黄土边坡可靠度的随机有限元分析

李 萍，白健忠，GRIFFITHS D V，李同录，郑亚楠

(1.长安大学 地质工程与测绘学院，陕西 西安 710054; 2. Department of Civil and Environmental Engineering, Colorado School of Mines, Colorado Golden 80401； 3.天水金泉矿业有限公司西安分公司，陕西 西安 710065)

0 引 言

强度参数变异性来源于土体固有的不确定性、取样的不确定性及试验的不确定性等。这些不确定性对岩土稳定性评价结果影响显著，需采用统计方法对参数进行处理和估计。如Lumb对海积土物理力学指标进行了不确定性估计，认为各指标线性回归后的残差服从正态分布[1]；冷伍明等给出了考虑各种主要不确定性因素和土性空间平均性质的土工参数空间变异系数的综合计算式[2]。近几十年来采用考虑土性变异性的可靠度理论解决岩土体稳定性问题，受到众多研究者的关注[3-6]。边坡稳定状态可以以边坡稳定系数(Fs)来表达，其表达式为

Z=Fs-1

(1)

式中：Z为状态函数，Z>0说明边坡处于稳定状态，Z=0说明边坡处于极限状态，Z<0说明边坡处于失稳状态。

岩土强度参数空间分布可分为确定场[7-9]和随机场[10]两种可靠度计算模型。确定场模型在一次稳定系数计算中，同一地层的随机变量(如黏聚力(c)、内摩擦角(φ))取值各处相同；随机场模型中，假定土性参数在一定距离内相关，引进了一个新的参数——自相关距离(δu)。自相关距离越小，反映了土体性质变化越快，随机性越显著；反之，自相关距离越大，反映了土体性质变化越缓慢，随机性越不显著。当自相关距离趋于无穷，土体强度参数的空间分布趋近于确定场。闫澍旺等分析认为天津塘沽地区软土的静力触探(CPT)锥尖阻力和侧摩阻力是一个平稳且具有各态历经性的随机场[11]。黄土为天然形成的非饱和土，强度不仅与固有的结构性有关，而且对含水状态敏感。李新生等通过对西安黄土的静力触探数据分析，显示黄土的自相关距离多在0.7 m以内，与其厚度相比，具有显著的空间随机性[12-13]。由此可见，利用随机场模型，能更准确地分析黄土工程的可靠度。

随机场模型可分为标准差折减法与随机场离散法。张健等分别对两类方法在边坡中进行了应用[14-15]。标准差折减法是将一层土点上获取的土性参数估计的标准差进行折减后，按确定场的方法估计边坡可靠度。其计算公式为

随机场离散法有中心点离散法、局部平均离散法及随机场插值法。近年来随机场离散网格与有限元网格相结合，应用较为广泛。Griffiths等采用局部平均离散法与有限元相结合，对高度10 m的简单土坡进行了可靠度分析[16]。但采用有限元进行随机场可靠度分析，存在计算效率低、商业有限元软件难以置入随机场等问题，李典庆等在处理该问题时进行了较多尝试[17-19]。

李萍等采用确定场分析黄土边坡可靠度，已经展开了一定的工作[20-23]，而采用随机场分析黄土边坡可靠度，目前还是空白。本文采用Fenton等提出的局部平均细分法[24]作为随机场离散法，建立黄土高边坡的随机场模型，对黄土高原大样本的强度参数变异性与概率分布进行统计，采用深静力触探孔(达第5层古土壤)的测试结果计算黄土垂直向自相关距离，采用弹-理想塑性有限元计算稳定系数，采用Monte-Carlo法计算失效概率，分析强度参数随机场对黄土边坡失效概率的影响，为随机场可靠度分析黄土边坡稳定性提供基础数据。

1 黄土强度参数变异性及其概率分布

本文收集了河南、山西、陕西、甘肃高速公路及铁路勘察报告中共4 018组黄土物理力学指标作为样本，统计结果见表1。数据东自山西吕梁山，西至甘肃永登，南自陕西西安，北至陕西靖边，基本覆盖了黄土高原的主要范围。由表1可见：密度及液、塑限的变异性较小，变异系数在0.10左右；而含水率、饱和度、内摩擦角和黏聚力的变异系数较大，都在0.25以上，尤其是黏聚力变异系数达0.73。将整个黄土高原的数据放在一起统计，会增大各指标的变异性。而按工点进行统计，更能代表场地特征。

表1 黄土物理力学指标统计结果Tab.1 Statistical Results of Loess Physical and Mechanical Indexes

图件引自文献[25]

侯晓坤等统计了黄土地区公路大桥、隧道、车站等79个工点的黏聚力和内摩擦角[25]。每个工点的统计样本数大于6个，自东向西按工点的地理位置排列，绘制变异系数随地域的变化(图1)。由图1可见，各工点黏聚力的变异系数(Vc)较大，内摩擦角的变异系数(Vφ)较小。黏聚力的变异系数多在0.30以上，内摩擦角的变异系数多在0.20以下。黄土强度参数变异性在黄土高原可明显地分成3段：甘肃马衔山以西地区多为砂黄土或粉土，变异性较小，黏聚力和内摩擦角的变异系数平均值分别为0.03、0.04；甘肃陇西、陇东，及陕西关中、富县、甘泉等地变异性较大，黏聚力和内摩擦角的变异系数平均值分别为0.38、0.14；而河南、山西等地变异性最大，黏聚力和内摩擦角的变异系数平均值分别为0.57、0.26。从西向东，黄土强度参数变异性增大，这主要是因为土的级配和矿物成分是决定黄土强度参数变异性的主要因素。黄土高原由西向东，细粒含量增高，尤其是黏粒含量增高，黄土强度对水的敏感性增强，随之变异性增大。图1中，Vc>0.80的样本多是山西西部的午城黄土，粒径小于0.005 mm的黏粒含量(质量分数)多超过30%，含水率低时黏聚力高，但随着含水率增高，黏聚力大幅减小，此时变异性最大。总之，黄土强度参数变异性具有区域性特点，这与其形成环境和黄土的成分、结构及对水的敏感性有关；而且不同时代的黄土由于其固结历史和固结应力不同，其性质也有差异。

由于工点的样本量少，难以完成概率分布的假设检验，所以根据黄土高原地貌和土性差异划分为10个区域[26]，根据垂直向按时代划分为Q1、Q2、Q3等3个地层，选取样本数大于50的区域地层黏聚力和内摩擦角进行假设检验，获得15组样本，检验结果如表2。采用K-S检验法进行正态分布和对数正态分布的假设检验。将n个样本测试值x的累积频率(Fn(x))与假设的理论概率分布(F(x))相比较，其最大差值(Dmax)满足条件可接受假设概型。其满足条件为

(4)

15组样本中，黏聚力有10组接受正态分布，8组接受对数正态分布，而内摩擦角有14组接受正态分布，仅有6组接受对数正态分布。由此可见，内摩擦角采用正态分布是合理的，黏聚力采用正态分布或对数正态分布接受的比例仅过半。

2 黄土自相关距离

表2 黏聚力和内摩擦角概率分布检验结果Tab.2 Test Results of Probability Distribution of the Cohesion and Internal Friction Angle

根据Vanmarcke建立的土体随机场模型[27]，计算自相关距离需要连续测试数据，一般采用静力触探结果。计算方法主要有递推空间法和相关函数法。闫澍旺等认为当样本容量足够大时，理论上这两种方法的计算结果是大致相同的[28-29]。假定黄土强度参数在空间上为齐次平稳的随机场，自相关距离计算公式为

δu=ΔzΓ2(Δz)h∞

(5)

δu=nΔz0Γ2(nΔz0)

(6)

因为地层厚度充分大时，自相关距离为一常数，Vanmarcke提出的递推空间法[27]建议采用绘制Г(n)-n曲线；当n取某一值后，折减系数趋向于平稳，找出平稳点n*；以该点为计算点，利用式(6)即可求出自相关距离。本文以陕西洛川地区所做的两个深静力触探孔数据[30]进行自相关距离计算。1#孔深度为37.1 m，2#孔深度为37.2 m，灌入达第5层古土壤。静力触探采用双桥探头，锥底直径35.7 mm，锥底面积10 cm2，摩擦筒表面积200 cm2，锥角60°，实测贯入速度0.5 m·min-1，每次贯入10 cm。图2为两个孔的静力触探曲线及其地层剖面。由于孔的样本数不可能无限大，在Г(n)-n曲线上一般难以找到平稳点n*。张梅等认为当地层厚度为有限值时，递推乘积的最大值应最接近自相关距离[29]，因此，ΔzΓ2(Δz)-Δz曲线的最大值接近自相关距离。本文采用这一方法，绘制ΔzΓ2(Δz)-Δz曲线(图3)。获得的自相关距离见表3。

图件引自文献[30]

李小勇等利用静力触探数据计算了太原和杭州粉土的垂直向自相关距离(δuv)，在1.0 m内；利用波速测试数据计算了水平向自相关距离(δuh)，在40 m左右[31]。高大钊利用静力触探数据对上海地区的典型地层垂直向自相关距离进行了统计，均值在1.60 m以内[32]。李新生等采用西安的静力触探数据计算了西安地区黄土地层的垂直向自相关距离，黄土和古土壤都在0.70 m以下[12]。杨勇等计算得到西安地区黄土的自相关距离平均值在0.30 m左右，最大值仅0.70 m[13]。表3列出的黄土垂直向自相关距离也多在0.70 m以内，仅有L1、L4和S0-S5的合并层大于0.70 m，最大值是采用2#孔S0-S5合并层的锥尖阻力计算值，达2.99 m。对地层厚度在10 m以内的地层，自相关距离在0.7 m以内进行考虑是可靠的；厚度在40 m以内的地层，自相关距离在3.0 m以内考虑是可靠的。

3 有限元随机场模型

以黏聚力和内摩擦角为随机变量，弹-理想塑性有限元与强度折减法结合计算稳定系数，弹性模量(E)取1×105kPa，泊松比(μ)取0.30，剪胀角(Ψ)取0，重度为确定值，取18.0 kN·m-3。有限元网格为正方形，8节点单元(图4)。坡高(H)的左边约束在距坡口线1.6H处，右侧约束在距坡脚0.4H处，底部约束在距坡脚下0.4H处。Monte-Carlo法计算失效概率(Pf)，随机抽样2 000次。采用Fenton等提出的局部平均细分法[24]进行随机变量离散，选用 Markovian 协方差函数。随机变量离散单元与有限元网格相匹配。图5、6分别为黏聚力和内摩擦角的离散结果及边坡变形，浅色网格代表黏聚力和内摩擦角低，深色网格代表黏聚力和内摩擦角高。水平向自相关距离采用40 m，垂直向自相关距离的分析采用了无量纲参数Θ。其表达式为

(7)

4 结果分析

4.1 强度参数概率分布类型

图7表示了黏聚力和内摩擦角随机分布的4种组合类型，其中横坐标是水平向与垂直向自相关距离的比值，取δuh/δuv>5的数据进行比较(黄土地层多符合这一条件)。由图7可见，黏聚力和内摩擦角都采用正态分布(normal(c)-normal(φ)组合)，将评估出较大的失效概率，比同时采用对数正态分布(log norm(c)-log norm(φ)组合)的失效概率大5倍以上。而黏聚力采用正态分布、内摩擦角采用对数正态分布的组合，以及黏聚力采用对数正态分布、内摩擦角采用正态分布的组合，居于二者之间，可见概率分布类型对可靠度分析有很大的影响。本文后续分析选用了黏聚力采用对数正态分布、内摩擦角采用正态分布的组合，这样一方面接近检验的结果，另一方面不会导致估计的失效概率过大或过小。

图8表示了稳定系数与失效概率的关系，其中取垂直向自相关距离为1.6 m，δuh/δuv=25。由图8可见：在Fs≤1.10条件下，3种参数概型组合下计算的失效概率差异较小；当Fs>1.20条件下，随着稳定系数的增长，失效概率显著分异，双对数正态分布型(log norm(c)-log norm(φ)组合)失效概率接近0，而双正态分布型(normal(c)-normal(φ)组合)失效概率高达41%。Griffiths等认为土的黏聚力和内摩擦角宜用对数正态分布，以避免随机抽样过程中负值的出现[33]。但本文通过对黄土强度参数的大样本假设检验，认为至少内摩擦角用正态分布为宜。《水电水利工程边坡设计规范》(SL 386—2007)对稳定系数的要求阈值在1.3左右[34]。在这个稳定系数条件下，参数概率分布类型对失效概率计算结果的影响不能忽略。本文选取黏聚力采用对数正态分布、内摩擦角采用正态分布的组合，这个组合有3个益处：①符合前述对黄土样本的假设检验结果；②估计的失效概率不会过大(与双正态分布组合比较)或过小(与双对数正态分布组合比较)；③避免了随机抽样过程中参数出现负值。由于黏聚力变异系数大，采用正态分布，将出现较多负值，而内摩擦角变异系数小，采用正态分布，出现负值的几率也小。

表3 递推空间法求解的垂直向自相关距离Tab.3 Autocorrelation Distances Calculated by the Recursive Space Method

图3 zГ2(z)-z曲线Fig.3 Curves of zГ2(z)-z

图4 有限元模型Fig.4 Finite Element Model

图5 随机变量黏聚力分布及破坏边坡变形Fig.5 Distribution of the Cohesion Which Is Random Variable and the Deformation of Unstable Slope

图6 随机变量内摩擦角分布及未破坏边坡变形Fig.6 Distribution of the Internal Friction Angle Which Is Random Variable and the Deformation of Stable Slope

图7 不同参数概型分布对失效概率的影响(Fs=1.50)Fig.7 Distribution Influence of Different Parameter Models on Failure Probability (Fs=1.50)

图8 不同参数概型下稳定系数与失效概率的关系 Fig. 8 Relationship Between Stability Factor and Failure Probability Under Different Parameter Models

4.2 自相关距离

图9 自相关距离对失效概率的影响(Vc=Vφ=0.50)Fig.9 Influence of Autocorrelation Distance on Failure Probability (Vc=Vφ=0.50)

计算取坡高为10 m，水平向自相关距离为40 m，黏聚力和内摩擦角变异系数都为0.50的边坡，其坡度为1∶1。自相关距离与失效概率的关系如图9所示。当Fs=1.05时，边坡接近极限状态，确定场(Θ=∞)的失效概率为49.4%，而随机场的失效概率随着垂直向自相关距离的增大而减小，在Θ接近0时，随机场的失效概率接近100%，当Θ>1.0时，即垂直向自相关距离大于坡高时，接近确定场的分析结果，随着垂直向自相关距离的继续增大，随机场的失效概率趋于稳定。在Fs=1.50时，代表边坡较为稳定，确定场的失效概率为17.3%，而随机场的失效概率随着Θ的增大而增大，当Θ接近0时，随机场的失效概率也接近于0%，当Θ>1.0时，随机场的失效概率接近确定场分析结果。黄土边坡高，而垂直向自相关距离小，Θ多小于0.2，可见垂直向自相关距离对黄土高边坡的可靠度影响显著。

图10绘制了在不同Θ下稳定系数与失效概率的关系，考虑了不同的内摩擦角情况，黏聚力和内摩擦角的变异系数取自表1。当失效概率大于41.0%时，相同的稳定系数条件下，Θ越小，失效概率越大；而当失效概率小于32.5%(这个界限因内摩擦角的不同略有差异)时，相同的稳定系数条件下，Θ越小，失效概率也越小；失效概率在32.5%～41.0%范围内，垂直向自相关距离对失效概率影响较小，与确定场的结果基本一致。Griffiths等采用随机有限元法分析边坡的失效概率与稳定系数的关系，认为其交点处的失效概率约为35.0%[35]。基于表1的平均变异系数水平，这个范围相对应的稳定系数多在1.10～1.20之间。由此可见，如果工程上选用的稳定系数在这个范围内，黄土的随机场特性对失效概率的影响并不大，而稳定系数大于1.20，随机场随着Θ越小，对失效概率影响越大。

图10 考虑随机场的黏聚力-内摩擦角边坡稳定系数与失效概率的关系Fig.10 Relationships Between Stability Factor and Failure Probability of c-φ Slope Against RFEM

4.3 变异系数

图11 不同自相关距离下的黏聚力变异系数与失效概率关系(Fs=1.30)Fig.11 Relationships Between Coefficient of Variation of the Cohesion and Failure Probability Against Varying Θ (Fs=1.30)

取坡高50 m，稳定系数为1.30的边坡，绘制图11的不同自相关距离下黏聚力变异系数与失效概率的关系，其中内摩擦角变异系数取了黏聚力变异系数的一半。由图11可见，黏聚力变异系数由0.10增加至1.50，失效概率由0%增加至100%，增长的趋势呈“S”型，即在变异系数小的条件下，失效概率增长速率较慢，随后快速增长，在接近100%失效概率条件下增长速率又变慢。自相关距离的Θ对“S”型失效概率的第一个转折点(起弯点)有较大的影响：Θ小，在黏聚力变异系数大的条件下起弯；Θ大，在黏聚力变异系数小的条件下起弯，起弯后失效概率随着黏聚力变异系数的增大而快速增大。当Θ=0.01时，黏聚力变异系数达0.75后起弯；当Θ=0.50时，黏聚力变异系数在0.30后就起弯。黄土的垂直向自相关距离小，表3中各层按地层厚度计算的Θ值最大为0.21， 37.0 m厚的S0—L5地层Θ值最大为0.08，可见黄土高边坡的Θ值非常小；而图1显示黄土的黏聚力变异系数多在0.30以上。上述两项指标相结合，失效概率多处于起弯段附近。黄土高边坡的这一土性特点意味着高边坡的稳定性具有极大的变异性，失效概率对土性的自相关距离和变异系数都极为敏感。图11中同时提供了确定场(即Θ=∞)的失效概率，发现在大部分黄土的变异系数与自相关距离水平下，确定场的失效概率大于随机场的失效概率，如Vc=0.50，Θ=0.10时，确定场失效概率为18.5%，而随机场的失效概率仅4.5%。由此可见，采用确定场进行可靠度分析，会大幅高估黄土高边坡的失效概率。但在极端变异系数水平下，如Vc=1.0，Θ=0.10时，确定场失效概率为35.4%，而随机场的失效概率达98.0%，采用确定场进行可靠度分析，又会大大低估黄土高边坡的失效概率。低估与高估之间有一个失效概率界限值，该值因随机变量的变异性及自相关距离不同而变化。变异系数越大，界限值越高；自相关距离越小，界限值越高。由此可见，考虑土性的变异性与自相关距离，对黄土高边坡可靠度的准确评估具有不容忽视的作用。

5 结 语

(1)黄土强度参数黏聚力和内摩擦角的变异性有系统差异。不同区域和不同时代黄土的黏聚力变异系数大多在0.30以上，内摩擦角变异系数小，多在0.20以下；且黄土高原从西向东，黏聚力和内摩擦角的变异性增大。

(2)强度参数(黏聚力和内摩擦角)的概率分布类型对边坡失效概率的影响极大。通过对不同地区和不同时代黄土的统计检验，认为黏聚力采用对数正态分布、内摩擦角采用正态分布为宜。该分布组合一方面接近检验的结果，另一方面不会导致估计的失效概率过大或过小。

(3) 采用两个孔深约37.0 m的黄土静力触探锥尖阻力和侧摩阻力计算了黄土地层的垂直向自相关距离，其值多在0.7 m以内，与计算厚度的比值多小于0.2。

(4) 在大部分黄土的变异系数与自相关距离水平下，采用确定场进行可靠度分析，将大大高估了黄土边坡的失效概率；在极端高变异系数水平下，采用确定场进行可靠度分析，又将大大低估黄土高边坡的失效概率。在黄土的平均变异系数水平下，失效概率在32.5%～41.0%范围内，对应的稳定系数在1.10～1.20之间，可不考虑随机场的影响，否则，考虑土性的变异性与自相关距离在评价黄土边坡可靠度时应受到足够重视。