-
花岗岩残积土广泛分布于我国东南沿海地区。不同地区气候条件不同，其风化程度及内部结构也具有较大差异，在福建、广东等省份，花岗岩残积土层厚度最大可达40 m[1]。花岗岩残积土是一种结构性很强的特殊土，具有不均匀性和各向异性、软化和崩解性、扰动性等工程特性。花岗岩残积土常存中间颗粒缺失、级配不良、空隙比大、压缩性中等或者偏高、物理力学性质变异性强等特点[2]，同时具有显著的各向异性和空间变异性。
许多学者对花岗岩残积土的性质进行了分析和研究。梁仕华等[3]分析了广州东部地区花岗岩残积土不同物理力学指标之间的相关性关系。陈洪江等[4]根据我国东南沿海地区花岗岩残积土的土工试验资料，研究了花岗岩残积土内摩擦角等指标的概率分布特征。Hu等[2]基于广东省花岗岩残积土土工试验数据，在概率统计和回归分析的基础上对其工程性质进行了评价。这些研究为增进各地区花岗岩残积土物理力学性质的理解起到了重要的作用，为涉及花岗岩残积土的岩土工程设计和分析提供了有用的参考。
在雨季期间，雨水入渗容易引发花岗岩残积土中负孔隙水压力消散，导致土体抗剪强度的降低，诱发边坡失稳灾害。许多研究者对降雨条件下的边坡稳定性进行了研究，如张灿亨等[5]，蒋中明等[6]。由于花岗岩残积土力学性质变异性强，一些学者开始采用可靠度理论对其进行稳定性分析。例如，郑敏洲等[7]采用可靠度理论研究天然状态下土体变异性对花岗岩残积土边坡稳定性的影响。Zhang等[8]采用拉丁抽样法，率先研究了降雨条件下花岗岩残积土边坡的可靠性。Wang等[9]使用简化Bishop法，从系统可靠度角度对降雨条件下的花岗岩残积土边坡的稳定性进行了分析。上述研究对揭示土体性质不确定性对花岗岩残积土边坡稳定性分析的影响起到了重要的作用。不足的是，上述研究中均假设花岗岩残积土的不确定性参数为随机变量，没有考虑土体参数空间变异性对边坡稳定性的影响。
土体的空间变异性表明，即使在同一层土体中，任何两点土体的性质也未必完全相同，但两者的相似程度会随着距离的减小而增加[10]。土体的空间变异性对岩土系统的失效模式、失效概率都有较大影响。从数学角度来讲，土体的空间变异性可采用随机场理论来模拟。为考虑土体空间变异性对边坡稳定性的影响，基于随机场原理的边坡可靠度分析正在成为研究热点[11,12,13,14,15,16]。例如，Griffiths和Fenton[11]使用随机有限元和蒙特卡罗方法分析了空间变异性对边坡失效概率的影响。吴振君等[13]引入已知取样点参数的约束随机场，在此基础上进行了边坡可靠度随机有限元分析。蒋水华等[16]提出考虑土体参数空间变异性的边坡可靠度及其敏感性分析多重响应面法，比较了自相关函数对边坡可靠度的影响。上述文献为研究空间变异性对花岗岩残积土边坡稳定性的影响提供了良好的借鉴。不足的是，这些研究中没有考虑降雨对边坡可靠度的影响。
本文的目的是采用随机场理论研究花岗岩残积土强度参数的空间变异性对强降雨条件下花岗岩残积土边坡稳定性评价的影响。论文将首先介绍本文采用的边坡可靠度模型，包括强降雨条件下花岗岩残积土边坡中的孔隙水分布、强度参数变化模型以及边坡安全性分析方法。在此基础上，提出了基于有限差分软件FLAC的可靠度模型求解方法。最后，通过具体算例研究了强度参数空间变异性对强降雨条件下花岗岩残积土边坡可靠度的影响。
HTML
-
降雨条件下边坡中的孔隙水压力分布可通过非饱和渗流分析进行计算，需要土体的水土特征曲线、渗流函数等参数[17]。Lumb[18]在对香港地区花岗岩残积土边坡失效的研究中提出湿润锋的概念，指出强降雨(即降雨强度大于土体饱和渗透系数的降雨)条件下花岗岩残积土边坡中将产生明显的湿润锋，如图1和图2所示。
Figure 1. Pore pressure distribution in a slope under intense rainfall infiltration
Figure 2. Illustration of the wetting front
对于残积土边坡，降雨诱发的滑坡失稳通常为浅层滑坡，深度大多在1.0～1.5 m左右[19]。在湿润锋以上为饱和或接近饱和的土体，土体基质吸力对土体强度影响不大；湿润锋以下为保持初始含水量状态的非饱和土体，基质吸力对土体强度有较大贡献。本文将采用Lumb[18]提出的湿润锋模型模拟强降雨条件下花岗岩残积土边坡中的孔隙水分布。
-
非饱和土体的强度不但与有效应力参数有关，还与基质吸力有关。在强降雨条件下土体中既存在饱和区，也存在非饱和区。因此，在对降雨条件下花岗岩残积土进行分析时，需采用非饱和土强度模型。由于Vanapalli模型[20]可用来模拟非饱和土在吸力达1 MPa范围内的强度特征，适用范围广，本文采用Vanapali模型模拟花岗岩残积土的强度。其公式为：
((1)) ((2)) 式中：τ为抗剪强度；c′为有效黏聚力；ϕ′为有效内摩擦角；(σn-ua)f为净法向应力；cT为总黏聚力；θ为体积含水率；θs和θr分别为饱和体积含水率和残余体积含水率；(ua-uw)为基质吸力。
在使用Vanapalli模型时需获得土体基质吸力与含水量之间的关系，即水土特征曲线。Fendlund和Xing水土特征曲线模型[21]适用范围广。本文采用Fendlund和Xing模型模拟花岗岩残积土的水土特征曲线，如下所示：
((3)) 式中：a、m、n为模型参数；e为自然对数的底；ψ为吸力。由于本文的重点是研究强度参数空间变异性对边坡稳定性的影响花岗岩残积土的水土特征曲线参数模拟成确定性参数。参考文献[8]中的花岗岩残积土水土特征曲线的数据库，本文中水土特征模型参数取值为θs＝0.453，θr＝0.148，a＝7.21 kPa，m＝0.74，n＝0.85。
1.1 强降雨条件下花岗岩残积土中孔隙水分布
1.2 非饱和土强度模型
-
随机场为定义在空间上的连续随机过程。在本文中将土体的强度参数考虑为平稳随机场，即该随机过程为平稳过程，每点处的随机变量具有相同的分布。已有花岗岩残积土力学性质的概率统计分析表明，花岗岩残积土的强度参数服从对数正态分布[4]。此外，现有边坡可靠度分析中，土体强度参数也多假定服从对数正态分布[11,12,14]。使用对数分布模型，可避免内摩擦角、粘聚力等出现负值这种物理上不可接受的情况。因此，本文假定内摩擦角、粘聚力服从对数正态分布随机平稳场。多元对数正态分布。对于服从对数正态分布的随机变量x～LN(μln,
((4)) ((5)) 式中：μ和σ为x的均值和标准差；μln和σln为ln x的均值和标准差。
自相关函数定义为随机过程中的任意两点处随机变量的相关系数，在本文的对数正态分布假定下，取为各随机变量对数间的相关系数。本文采用了高斯自相关函数[15]：
((6)) 式中：δx和δy分别为水平方向和数值方向上的相关距离，量度了随机变量的空间变异性；Δx和Δy分别是两点间水平方向和竖直方向上的距离。
从式(6)可以得出，当距离减小，相关系数逐渐增大。当距离为0时，相关系数为1，这可以保证土体参数的空间连续性。距离越大，则相关系数越小。当Δx＝δx，Δy＝0时(或Δx＝0，Δy＝δy时)，相关系数ρ＝0.043。在相关距离外，可以认为相关性很小。根据式(6)可对土体中某一点对周围土体作出等相关系数线，可知它们一般为椭圆；当δx＝δy时，等相关系数线为圆，随机场变为各向同性随机场。
对于一般沉积土层，水平相关距离δx通常比δy要大很多，这是由沉积土层沉积形成的特性决定的[22]。对于一般的沉积土层，垂直相关距离约在0.2～2.0 m之间，水平相关距离一般介于20～80 m，有的甚至大于100 m[15]。由于花岗岩残积土由花岗岩原地风化[1]，可能会表现出不同于一般沉积的不均匀性和各向异性。目前，花岗岩残积土空间相关距离的研究还少见报道。本文将采用参数敏感性方法研究花岗岩残积土空间距离对边坡稳定性的影响。
-
随机场的离散方法主要有局部平均法、中心点法、形函数法、K-L正交展开法[14]等。本文采用K-L正交展开法。K-L正交展开法采用一组两两不相关的正交单位随机向量来离散随机场[14]。该方法具有可根据需要任意控制精度、离散效率高、且不受单元形状的限制等优点[22]。作为示例，图3给出了某边坡的粘聚力参数在δx＝10 m、δy＝2 m时根据K-L展开获得的一组离散样本，其中网格颜色灰度的大小代表了土体参数值的大小。由于上述相关距离参数土体在水平方向具有更多的相关性，图3中土体表现出一定的分层特征。
Figure 3. A realization of random field of cohesion
2.1 强度参数的随机场模型
2.2 随机场离散
-
本文基于FLAC有限差分软件中的FISH语言编程，使用强度折减法[23,24]计算边坡的安全系数。强度折减法中目前应用较广泛的失稳判据主要有3种[23]：(1)边坡位移发生突变；(2)计算结果不收敛；(3)潜在滑面塑性区贯通。FLAC软件中采用最大不平衡力比率来判断计算结果的收敛性。
设原强度指标为c′和ϕ′。计算时，使用强度折减系数F对强度指标进行折减，则假想的强度指标[24]分别为：
((8)) ((9)) 强度折减法求解边坡安全系数的过程即为搜索F值，使得边坡刚好位于稳定与不稳定的临界状态。该安全系数可采用二分法进行求解。
首先设定初始求解区间[F10，F20]，初始区间过大会使得求解过程过多，增大了求解时间，而初始区间过小则可能遗漏可能值。一般可以设为[0，2]，对于安全系数可能较高的边坡，则可以进行适当扩大。然后计算初始中值
((10)) 再使用式(8)和(9)计算假想强度指标，进行第一次循环求解。将假想强度指标代入边坡稳定性分析。如果在假想强度指标条件下边坡失效，则说明此时强度折减系数过大，下一次循环取下半区间的中值；如果在假想强度指标条件下边坡未失效，则说明此时强度折减系数过小，下一次循环取上半区间的中值，即：
((11)) 直到求解区间长度满足精度要求，取该区间中值即可作为安全系数值。
-
随机场离散后会产生较多的随机变量。由于蒙特卡罗法的精度不受随机变量数目的影响，本文采用蒙特卡罗法计算边坡的失效概率。令x代表边坡可靠度模型中的不确定性参数。定义函数：
((12)) 假设蒙特卡罗法中生成的样本数为N，则边坡的失效概率估计值为：
((13)) 本文使用Matlab进行随机场离散，采用FLAC软件求解边坡安全系数；通过建立两者之间的接口程序，自动实现边坡失效概率的自动求解。具体过程如下：
1)使用FLAC根据模型参数建立边坡模型，输出网格数据文件。
2)采用Matlab读取网格数据文件，结合随机场参数，对随机场进行离散，生成N个随机场样本(N为蒙特卡罗模拟次数)。
3)利用FLAC自带的FISH程序读入离散后的随机场样本，并对各单元参数进行赋值。
4)进行稳定性分析计算，保存分析结果。
5)重复步骤3)和步骤4)，直到读取、分析所有的随机场样本。
6)使用Matlab读取计算结果，对结果进行后处理。
3.1 强度折减法
3.2 失效概率计算方法
-
参考文献[9]，选取了如图4所示的花岗岩残积土边坡算例。图4中z为湿润锋深度，即饱水层的厚度。坡高为15 m，坡角为45°，边坡土体的重度为常量γ＝19.5 kN/m3。
Figure 4. Slope test problem in this paper
参考有关花岗岩残积土的统计分析结果[1,3,4]，粘聚力具有较高的变异性，而内摩擦角的变异性较小。各参数选取的均值、标准差、中位数如表1。其中cT的均值为基质吸力为12 kPa时，基于式(1)～式(3)计算。表中同时注明了各强度参数的中位数。
参数 c/kPa ϕ1＝ϕ0/(°) cT/kPa 均值 6.8 32.3 11.7 标准差 9.0 12 9.0 中位数 4.1 30.3 9.3 Table 1.
Soil strength parameters -
为研究湿润锋深度对边坡失效模式的影响，本节先不考虑土体的空间变异性，土体取强度参数中位数进行边坡稳定性分析。图5(a)、(b)分别给出了湿润锋深度为0.5 m和2.0 m条件下根据强度折减法计算得到的边坡失效时的剪应变云图。当湿润锋深度较小时(图5〈a〉)，此时的滑动面同时穿越表层饱和土层和湿润锋以下非饱和土层，且滑动面主要位于非饱和土层中。当湿润锋深度较大时(图5〈b〉)，此时滑动面位于饱水层内，受饱水层的土体控制，这与文献[9]的结果相似。
Figure 5. Shear strain contour under different failure types: (a) z＝0.5 m; (b) z＝2 m (Dash lines denote the wetting fronts)
-
图6(a)、图6(b)分别给出了在湿润锋深度z＝1 m、土体相关距离为δx＝10 m、δy＝2 m条件下两次蒙特卡罗模拟中边坡失稳时刻的剪应变图。图6(a)中，滑动面主要位于湿润锋以上的饱和土层中。在图6(b)中，滑动面主要位于湿润锋以下，且形成了两个潜在滑动面。这两者的差别主要是由于两次蒙特卡罗模拟中生成土体性质的差异性造成的。
Figure 6. Shear strain contour under different failure types: (a)δx＝10 m, δy＝2 m，z＝1 m；(b)δx＝10 m, δy＝2 m，z＝1 m (Dash lines denote the wetting fronts)
图7(a)、图7(b)分别给出了在湿润锋深度为1 m、土体相关距离分别为δx＝10 m、δy＝2 m和δx＝100 m、δy＝20 m条件下一次蒙特卡罗模拟中土体内聚力的分布图，其中网格颜色灰度的大小分别代表了粘聚力值高低。由图7可知，随着相关距离的减小，网格颜色的变异程度增加，表明土体强度参数变异程度的增加。图7中同时给出了失稳时变形网格图以及对应的滑动面位置。在边坡稳定性分析中，边坡总是沿着安全系数最小的滑动面滑动。当土体空间变异性较大时，更容易形成路径较短的软弱滑动面，导致滑坡规模较小。
Figure 7. Deformation of the slope at the moment of failure under different correlation distances：(a)δx＝10 m, δy＝2 m，z＝1 m; (b)δx＝100 m, δy＝20 m，z＝1 m (Dash lines denote the slip surfaces)
-
首先不考虑土体的空间变异性，采用随机变量方法对该边坡的可靠度进行分析。此时，土体对应的相关距离为δx＝∞、δy＝∞。图8给出了此种情况下边坡失效概率与湿润锋深度的关系。作为比较，图中还给出了不同湿润锋深度下边坡的的中值安全系数。由图可知，随着湿润锋的增加，边坡安全系数不断下降，而边坡的失效概率不断增加。当湿润锋深度为2.0 m时，边坡的失效概率达到45%以上。当考虑土体参数不确定性后，边坡的安全系数会存在不确定性，导致边坡是否稳定也存在不确定性。此时，采用失效概率衡量边坡的稳定性更能合理地反映不确定性对边坡稳定性评价的影响。
Figure 8. Results from RFEM under different correlation distances
为研究土体空间变异性的影响，图8中还给出了土体参数空间相关距离分别取为δx＝10 m，δy＝2 m条件下边坡失效概率与湿润锋深度之间的关系。在考虑空间变异性条件下，失效概率随着湿润锋深度的增加而增加，趋势与不考虑空间变异性的情况相同。不过，考虑土体空间变异性条件下失效概率对湿润锋深度变化更为敏感。其原因在于，考虑土体空间变异性后，土体不同点处的变异性可以相互抵消，导致滑动面土体的平均性质的变异性降低、安全系数的不确定性减小，在文献[14]中也有类似的结论。这也导致在安全系数较低时，不考虑空间变异性会低估失效概率。
图中还给出了δx＝100 m，δy＝20 m以及δx＝1 000 m，δy＝200 m条件下边坡的失效概率。随着相关距离的增加，土体不确定性的空间折减效应降低，导致边坡安全系数不确定性增加。此时，边坡失效概率对湿润锋深度的敏感性介于δx＝10 m、δy＝2 m与δx＝∞、δy＝∞之间。
-
图9为使用不同粘聚力和内摩擦角的相关系数时，随机有限差分蒙特卡罗方法的计算结果。对应湿润锋深度也标于图中。由图9可知，在湿润锋深度相同的情况下，边坡失效概率随着粘聚力和内摩擦角相关系数的增加而增加。
Figure 9. Probability of failure versus correlation coefficient of cohesion and friction angle
粘聚力和内摩擦角相关系数为正数时，二者可能同时取得较大值或较低值。相关系数为负数时，当其中之一取得较低的值时，另一个参数比较容易取到相对大值，降低了粘聚力和内摩擦角同时出现较低值的情况，使得边坡的失效概率降低。已有资料表明[25]，粘聚力和内摩擦角具有一定的负相关性。此时，不考虑粘聚力和内摩擦角的负相关性会高估边坡的失效概率，造成设计的浪费。文[14]中也观察到了类似的现象。
4.1 湿润锋深度对边坡失效模式的影响
4.2 空间变异性对滑坡失效模式的影响
4.3 空间相关距离对边坡稳定性的影响
4.4 粘聚力和内摩擦角相关性的影响
-
本文中的研究工作及结论可总结如下：
1)以强降雨条件下的分层花岗岩残积土边坡为研究对象，建立了考虑花岗岩残积土空间变异性的随机场有限差分模型，提出了基于现有软件的实现方法。
2)当湿润锋深度较浅时，边坡最危险滑动面可越过湿润锋，进入到非饱和土层中；随着湿润锋深度的增加，边坡最危险滑动面位于湿润锋以上的饱和土层中。
3)在湿润锋深度相同的条件下，土体空间变异性会导致边坡出现不同的失效模式。土体空间变异性越强，边坡失稳的规模越小。
4)随着土体空间变异性的增加，土体不确定性的空间折减效应增大，安全系数的不确定性降低，边坡失效概率对湿润锋变化的更加敏感。
5)边坡失效概率随粘聚力和内摩擦角之间的相关系数的增加而增加。粘聚力和内摩擦角间通常呈一定的负相关性。此时，忽视两者之间的负相关性会高估边坡的失效概率。
DownLoad:
