土质边坡的单元失效概率与失效模式研究

彭 普，李 泽，张小艳，申林方，许 芸

(1.昆明理工大学建筑工程学院，云南，昆明 650500；2.昆明理工大学电力工程学院，云南，昆明 650500；3.长江勘测规划设计研究有限责任公司，湖北，武汉 430010)

边坡的稳定性是与众多随机参数息息相关的，其中土体抗剪强度参数以及地下水的影响尤为重要[1]。天然土体的抗剪强度参数在空间上呈现变异性，研究表明土体参数空间变异性使得边坡的失效模式和安全性能均具有一定的不确定性，因此需要采用可靠度方法进行研究[2-3]。地下水从两个方面影响着边坡的稳定性：一是其使得边坡内部的渗流场发生变化，渗流场的变化导致边坡孔隙水压力发生改变，从而影响着边坡的安全性能；二是其降低边坡土体的材料抗剪强度参数，从而影响着边坡的安全性能[4]。通常情况下地下水位并不是一个确定值，其受众多因素影响。因此，边坡的稳定性问题是一个与土体参数空间变异性和地下水位随机性相关的不确定性问题，在进行边坡稳定性分析时应综合考虑这两个因素带来的影响。

近年来，边坡系统可靠性问题受到了广大学者的关注。基于极限平衡分析方法(LEM)可以直接得到边坡稳定性系数的数学分布以及边坡的失效概率，但LEM 在进行边坡稳定性计算时首先需要假定滑裂面，然后对滑裂面上的土体进行条分并假定条间力的大小及方向将超静定问题转化为静定问题进行求解[5-8]；基于有限元分析方法(FEM)可以得到边坡的应力-应变分布，还可以获得比LEM方法更合理的稳定性分析结果，但计算成本较高[9-11]；基于有限元塑性极限分析方法可以在不考虑边坡加载历史和材料本构关系的情况下，高效、准确地得到边坡处于极限状态下的稳定性系数和失效模式[12-18]。该分析方法主要通过构建静力许可应力场(下限法(LBM))和机动许可速度场(上限法(UBM))，通过塑性力学的上、下限定理可知，真实解必然位于上限解与下限解之间。张小艳等[19]、李亮等[20]将土体抗剪强度参数设为随机变量，提出了基于塑性极限分析的边坡可靠度上、下限法；陈朝晖等[21]提出了考虑参数空间变异性的边坡可靠度上、下限分析方法。以上研究促进了有限元塑性极限分析方法在边坡可靠度中的应用。然而，随机地下水位对边坡可靠度的影响并未得到系统的研究。

此外，传统的边坡失效概率分析采用整体失效概率法，其主要采用边坡稳定性系数的阈值判断边坡是否失效，并据此计算边坡的失效概率，该方法隐含假定边坡具有单一的失效模式，这与边坡存在多种失效模式不相符。为解决这一问题，李典庆等[22]、蒋水华等[23]提出了基于代表性滑动面的边坡失效概率计算方法，然而在进行边坡失效模式统计时，随着计算样本数的增加计算量会急剧增大，其计算效率和计算精度均较低[24]；张小艳等在文献[19]中提出了基于上限法的单元失效概率计算方法，利用边坡稳定性系数和速度场进行边坡单元失效概率的计算，为边坡系统失效概率分析提供了一种新思路。

鉴于此，本文将有限元离散技术、上限法理论、相关非高斯随机场模拟方法以及随机规划理论结合起来提出了随机地下水位作用下考虑参数空间变异性的边坡可靠度分析方法。

1 随机地下水位作用下考虑参数空间变异性的边坡可靠度分析

1.1 随机地下水位的边坡稳定渗流场

当前，地下水位通常被设定为某一确定值，在进行边坡可靠性分析时，边坡的极限状态函数不包含地下水位这一随机变量。然而由于土体颗粒和孔隙在边坡中随机分布使得边坡渗流场具有不确定性，此外降水、地表径流、地下汇流等对地下水的供给以及蒸发、抽水等对地下水的排出的不确定性使得边坡地下水位具有不确定性，不确定的地下水位将导致不确定的渗流场。近年来，对于渗流场的研究主要集中在土体材料变化导致的渗透系数变化，进而导致渗流场的变化上[25-26]。对于随机地下水位对渗流场的影响以及对边坡可靠性的研究成果依然较少。

本文将地下水位定义为随机变量，如图1 所示，假定其服从正态分布[27]，并假定地下水位的上界和下界，地下水位位于上界和下界之间随机变化。图1 中：为水位的下界；为水位的上界；为水位的均值。为了简化计算，本文假设随机地下水位作用下的渗流场为饱和稳定渗流场，并且不考虑水位突然升高或降低引起的超孔隙水压力以及对土体抗剪强度参数的影响，且位于浸润线以上的孔隙水压力p=0，位于浸润线以下的孔隙水压力p=γh，其中 γ为水的容重，h为浸润线上与计算点处于同一条等势线上的点到计算点的垂直距离。采用二维稳定渗流公式进行渗流分析[28]，具体公式如下：

图1 边坡随机地下水位作用示意图Fig.1 Schematic diagram of random groundwater level of soil slope

式中：kx为土体材料x方向的渗透系数；ky为土体材料y方向的渗透系数；Hr为土体内各点的随机水头函数。

1.2 描述参数空间变异性的随机场模型

当前普遍采用自相关函数表征土体的空间变异性，常用的自相关函数类型有指数型、高斯型、对数型、二阶自回归型、指数余弦型以及三角型[29]。考虑边坡的可靠性对自相关长度较为敏感，对自相关函数的形式不敏感[30-31]，本文选用数学表达式更为简单的指数型自相关函数进行研究，具体可采用式(2)表示：

式中：xa、xb为空间坐标；xaa、xba、xab、xbb为空间坐标的坐标分量；Lh为水平坐标方向的波动范围；Lv为竖直坐标方向的波动范围[31]。

由于土体参数之间存在一定的互相关性和土体参数本身存在一定的自相关性，因而涉及相关非高斯场的离散过程，本文采用类似文献[32]的方法，采用基于乔列斯基分解的中点法进行随机场的生成，对于指数型自相关函数，土体黏聚力、内摩擦角的相关对数随机场可表示为：

式中：m=c,φ，c为土体黏聚力，φ为土体内摩擦角；(xa,ya) ∈Ω；Hm(xa,ya)为相关非高斯随机场；为相关标准高斯随机场；σlnm=为对数正态分布m的均值， σm为对数正态分布m的标准差；ulnm为相应正态变量lnm的均值；σlnm为相应正态变量lnm的标准差。

1.3 边坡可靠度分析的随机规划模型

有限元塑性极限分析兼具FEM 法和极限分析的优点，成功克服了结构体在进行极限分析时静力许可速度场的困难。SLOAN 等[14-15]将土体采用三角形单元进行离散，通过构建三角形单元塑性流动约束条件、公共边速度不连续约束条件、三角形单元速度边界约束条件从而构建机动许可速度场。由上限定理易知，机动许可速度场和外荷载是一一对应的关系，其中最小外荷载无限接近于极限荷载，因此上限法可理解为求解极小值的数学规划问题。本文在SLOAN 等[14-15]、张小艳等[19]、王均星等[33]研究的基础上，建立随机地下水作用下考虑土体参数空间变异性的边坡可靠性分析随机规划模型如下：

式中：Z为边坡可靠度计算的极限状态函数；km为边坡稳定性系数；kγ为土体容重超载系数；We为有限单元的内功功率；Wd为有限单元公共边的内功功率；WG为自重在有限单元节点速度上所做的外功功率；为孔隙水压力在有限单元连续体内所作的外功功率；Wdp为孔隙水压力在有限单元公共边上所作的外功功率；λe为塑性乘子；λd为决策变量；为有限单元e的塑性流动约束条件矩阵；为相邻有限单元公共边d的塑性流动约束条件矩阵；ab为有限单元速度边界约束条件矩阵；ue为有限单元e的速度矢量；ud为相邻有限单元公共边d的速度矢量；ub为边界上有限单元的速度矢量。

采用三角形单元进行边坡离散，在随机地下水位作用下考虑土体参数空间变异性时，不同空间位置的单元将具有不同的抗剪强度参数，式(4)中有限单元的内功功率为：

式中：ce为有限单元e形心处的黏聚力；φe为有限单元e形心处的内摩擦角；Ae为有限单元e的面积，e=(1,2,···,Ne)，Ne为所有有限单元数量。

有限单元公共边的内功功率为：

式中，ld为公共边的长度，d=(1,2,···,Nd)，Nd为所有有限单元公共边的数量。

自重所作的外功功率为：

孔隙水压力所作的外功功率为：

2 随机规划模型的求解策略

式中：tw=(1,2,···,nw)，nw为边坡地下水位蒙特卡洛随机数的数量；为边坡地下水位的第tw个随机数；μw为边坡地下水位的均值； σw为边坡地下水位的标准差；Random为截尾正态分布随机数生成函数；Normal为地下水位随机数符合正态分布；Hlb为地下水位的下界，取最低水位；Hub为地下水位的上界，取最高水位。

2)生成边坡材料抗剪强度参数的随机场。假设土体材料的黏聚力cr和内摩擦角φr的自相关函数类型为指数型，采用乔列斯基分解的中点法进行随机场的生成：

本文采用容重超载的方式使边坡达到极限状态，随机地下水位作用下考虑参数空间变异性时，边坡的容重超载系数为：

式中：kγ(tm,tw)为第tw个地下水位作用下对应着第tm个随机场的边坡的容重超载系数；、为土体处于极限状态时的容重，其与以及相关； γa为土体实际容重；为强度折减以前第tm个非高斯随机场中有限单元e形心处的黏聚力和内摩擦角。

随机地下水位作用下考虑参数空间变异性时，边坡稳定性系数如下：

式中：tw=(1,2,···,nw)；tm=(1,2,···,nm)；、为强度折减以后第tm个非高斯随机场中有限单元e形心处的黏聚力和内摩擦角。

3)计算边坡稳定渗流场。将步骤1)生成的地下水位的随机数tw=(1,2,···,nw)，逐次代入边坡稳定渗流场的计算公式，进而得到离散单元e三个节点的孔隙水压力值：，e=(1,2,···,Ne)，tw=(1,2,···,nw)。

图2 边坡可靠度上限法流程图Fig.2 Flow chart of reliability analysis using upper bound method for slope

5)计算边坡稳定性系数并绘制相关曲线。

3 边坡系统失效概率和可靠度指标

边坡工程中普遍采用边坡稳定性系数来评价边坡的稳定性。当边坡稳定性系数km≥1，则认为边坡稳定；边坡稳定性系数km<1，则认为边坡失稳[34]。因此，边坡的失效功能函数按式(14)进行表示：

式中：tw=(1,2,···,nw)；tm=(1,2,···,nm)；I(tw,tm)为第tw个地下水位作用下对应着第tm个随机场的边坡失效功能函数。

根据式(14)边坡失效功能函数，可以得到边坡在第tw个地下水位作用下的失效概率：

式中：tw=(1,2,···,nw)；为边坡在第tw个地下水位作用下的整体失效概率。

式(15)被广泛应用于边坡失效概率的计算，然而隐含地假定了边坡仅具有单一失效模式，这与边坡存在多种失效模式是不相符的[35-36]，此外该式只考虑了边坡稳定性系数的影响并未考虑边坡的失效后果的影响。为克服上述不足，张小艳等[19]提出了单元失效概率法用以分析边坡的失效概率。该方法充分发挥有限元离散技术在构建速度场上的优越性，认为在每一个失效模式对应的速度场中，当单元速度大于0 时代表此单元已发生塑性流动(失效)、当单元速度等于0 时代表此单元未发生塑性流动(安全)。在此基础上，本文首次提出根据边坡中单元的失效信息来估算边坡的系统失效概率。边坡系统是由无数个潜在滑动面组成的串联系统，其系统失效概率近似等于各个代表性失效模式的失效概率之和。当前，采用LEM 法和FEM 法进行边坡系统可靠度计算时任意潜在失效模式(设为N)均会导致边坡发生失稳破坏，因此边坡系统失效的功能函数为：

式中：Zs为边坡系统失效的功能函数；Zi为边坡的第i种潜在的失效模式，i=(1,2,···,N)。

根据边坡系统失效功能函数的定义，可以得到边坡的系统失效概率估算公式为：

式(17)在进行边坡系统失效概率计算时需要在大量的潜在失效模式中识别最具有代表性的失效模式，因此，所识别的代表性失效模式对计算结果有较大的影响。塑性极限分析上限法可以直接得到边坡极限状态下的失效模式，每种失效模式均是由若干个失效单元组成的区域，鉴于此，本文定义边坡单元的失效功能函数如下：

式中：Ie(tw,tm)为第tw个地下水位作用下对应着第tm个随机场中单元e的失效功能函数；uec(tw,tm)为第tw个地下水位作用下对应着第tm个随机场中单元e的速度。

从分类角度看，采用边坡稳定性系数的阈值和有限单元的节点速度信息可以识别边坡所有失效模式。AP(Affinity Propagation)聚类算法非常适合于本文边坡失效模式的识别，进而进行边坡系统失效概率的估算，该算法的主要思想是将tw×tm个速度场的失效面积都当作潜在的聚类中心，各失效面积之间连线构成一个网络，再通过网络中各条边的消息传递计算出各样本的聚类中心。第m、j个速度场的失效面积Am、Aj，Aj作为Am聚类中心的能力，记为S(Am,Aj)，S(Am,Aj)越大，两个点距离越近；定义吸引度R(Am,An)为An作为Am聚类中心的适合程度，如图3(a)所示为Am选An的过程；定义归属度B(Am,An)为Am选择An作为其聚类中心的适合程度，如图3(b)所示为An选Am的过程。R(Am,An)与B(Am,An)的和越大，表示An作为聚类中心的可能性越大，Am隶属于An为聚类中心的可能性也越大。具体的，吸引度迭代公式为：

图3 AP 聚类算法选点过程Fig.3 AP clustering algorithm selection process

式中：

Rt+1(Am,An)为更新后的R(Am,An)；Rt(Am,An)为更新前的R(Am,An)；λ 为阻尼系数。

归属度迭代公式为：

式中：

Bt+1(Am,An)为更新后的B(Am,An)；Bt(Am,An)为更新前的B(Am,An)。

定义单元失效概率为：

根据边坡可靠度分析上限法数学规划模型的求解策略中所得nw×nm个边坡稳定性系数km(tw,tm)，通过统计可以得到边坡在第tw个地下水位作用下稳定性系数的均值、标准差以及在所有可能发生的地下水位作用下边坡稳定性系数的均值、标准差。具体为：

式中：tw=(1,2,···,nw)；μk(tw) 、σk(tw)分别为第tw个地下水位作用下边坡nm个稳定性系数的均值和标准差；μk、 σk分别在所有可能发生的地下水位作用下nw×nm个边坡稳定性系数的均值和标准差。

4 算例分析

根据本文提出的方法，编制了相应的上限法程序，对一个经典的边坡算例进行了计算分析并与极限平衡法、有限元法计算结果对比，验证了本文计算方法的正确性。

4.1 均质边坡基本信息

本文选取文献[37]的均质边坡算例，来验证本文计算方法的正确性。图4 为均质边坡计算模型，边坡高度为10 m，坡比为1∶1，坡顶宽度为10 m，诸多文献对此边坡的失效概率进行过分析[29,32,37]。本文采用三角形有限单元对边坡进行离散，共得到885 个有限单元、2655 个有限元节点、1279 条公共边。边坡左侧随机地下水位为，右侧坡脚的地下水位与地表齐平，设置P1、P2、P3 三个孔隙水压力关键点(坐标分别为：P1(5, 5)、P2(10, 5)、P3(15, 5))；设置有4 个有限单元失效概率监测单元，分别为坡顶的E1 单元、边坡中部的E2、E4 单元、坡脚的E3 单元，其中单元E1、E2、E3 均在边坡的临空面上，E1、E4 在同一竖直面上，E2、E4 在同一高程上。设定土体容重为γ=20 kN·m-3、渗透系数为K=5×10-7m/s，两者均为确定值，其他计算参数见表1。本算例的计算目的是：1)计算随机地下水位作用下考虑参数空间变异性的边坡稳定性系数分布类型；2)计算随机地下水位作用下考虑参数空间变异性的边坡稳定性系数概率密度曲线和累积概率密度曲线；3)采用AP 聚类算法进行边坡失效模式的识别并进行边坡系统失效概率的估算；4)计算边坡的整体失效概率以及单元失效概率与地下水位之间的关系。

表1 土体参数统计特性Table 1 Statistics of soil parameters

图4 均质边坡计算模型 /mFig.4 Homogeneous slope calculation model

4.2 随机地下水位和渗流场

图5 地下水位随机数分布直方图Fig.5 Histogram of random distribution of groundwater level

图6 边坡稳定渗流场Fig.6 Steady seepage field of slope

图7 关键点处孔隙水压力随机变化情况Fig.7 The pore water pressure of the monitoring points

4.3 土体抗剪参数空间变异性和随机场

本文假定土体黏聚力c和内摩擦角φ服从对数正态分布，采用指数型自相关函数模拟土体抗剪强度参数的空间变异性。土体参数随机场的数量nm取1000；根据式(10)、式(11)，生成边坡材料抗剪强度参数的1000 个随机场。图8 为tm=500时边坡抗剪强度参数随机场，在该随机场下土体黏聚力的最大值为17.7961 kPa，最小值为3.5555 kPa；内摩擦角的最大值为46.8418°，最小值为17.8413°。此外，土体黏聚力和内摩擦角在空间上呈现一定的负相关。

图8 边坡抗剪强度参数随机场Fig.8 Random Field of Slope Shear Strength Parameters

4.4 上限法与LEM 法、FEM 法的对比分析

文献[37]采用150 项K-L 级数展开的方法离散土体的随机场，通过5×104次MCS 模拟，得到边坡的可靠度；文献[32]采用基于乔列斯基分解的中点法(CD)模拟随机场，采用拉丁超立方抽样(LHS)产生土体的黏聚力和内摩擦角样本，进而进行边坡可靠度的计算。为验证本文计算方法在无地下水作用时的有效性，边坡计算模型以及土体参数均与文献[32, 37]一致，根据边坡可靠度分析随机规划模型式(4)和随机规划模型的求解策略。当无地下水作用时，离散单元三个节点的孔隙水压力值：均为0，根据第4.3 节土体参数随机场数量nm=1000，可以得到1000×1=1000 个计算样本，基于开发的并行计算程序，在一个小型工作站(处理器：AMD 锐龙3970X、32 核心，物理内存：128 GB)上计算1000 个样本花费了大约0.43 h(单线程计算约17.2 h)，平均每个样本的计算时间为1.56 s(单线程计算约61.92 s)。此外，使用有限元极限分析软件OptumG2 的FEM法对该边坡进行了稳定性分析，计算相同样本共花费了大约23.4 h，平均每个样本计算时间为84.24 s。本文并行计算程序相较于单线程计算效率提高约38.69 倍，相较于FEM 法计算效率提高约53 倍。计算结果如表2 所示，本文方法与K-L 法、LHS法相比，所得边坡稳定性系数的均值、标准差基本相同，所得失效概率略微偏小；本文方法与FEM法相比，所得边坡稳定性系数的均值、标准差基本相同，所得失效概率略微偏大。

表2 不同方法的边坡可靠度结果对比Table 2 Comparison of slope reliability results of different methods

为验证本文计算方法在地下水作用时的有效性，本文选取低、中、高三个地下水位进行对比分析。表3 为本文计算方法与LEM 法、FEM 法计算结果对比情况，其中选取的水位分别为tw=1、tw=25、tw=50。图9 为三种水位下边坡稳定性系数分布。

表3 三种水位作用下边坡的可靠度指标Table 3 Calculation results of reliability index

图9 边坡稳定性系数分布特征Fig.9 Distribution characteristics of slope safety factor

计算结果表明：

1)采用UBM 计算所得边坡稳定性系数均值的上限解大于LEM 法、小于FEM 法所得计算结果，但误差较小。随着水位的升高三种计算方法所得边坡稳定性系数均降低。在tw=1、tw=25、tw=50三种情况下，UMB 法所得边坡稳定性系数均值均大于相同条件下LEM 法所得的均值，符合上限解的特征。此外，随着地下水位的升高，边坡整体失效概率随之增加，这与实际是相符的。在三种地下水位情况下，采用UBM 法计算所得边坡整体失效概率均小于LEM 法所得计算结果，但误差较小，符合上限解的特征，根据上限定理：机动许可速度场对应的荷载是极限荷载的上界，即采用上限法所得边坡稳定性系数必然大于真实边坡稳定性系数，因此在统计边坡的整体失效概率时，会略微低估边坡的失效概率。然而，在tw=1、tw=25、tw=50三种情况下，UMB 法所得边坡稳定性系数均值均小于相同条件下FEM 法所得的均值，其原因主要是使用有限元极限分析软件OptumG2在进行边坡稳定性分析时，依据计算不收敛、塑性区贯穿、计算点位移发生突变作为边坡处于极限状态的判断标准，因此FEM 法所得边坡稳定性系数是严格的上限解还是下限解无法验证。

2)图9(a)为tw=1、tw=25、tw=50三种情况下，UBM 法与LEM 法、FEM 法所得边坡稳定性系数累积概率密度曲线，不难看出，UBM 法与LEM法所得边坡稳定性系数累积概率密度曲线都非常接近，误差较小；UBM 法所得边坡稳定性系数累积概率密度曲线都位于FEM 法左侧，其原因主要是UBM 法和FEM 法在进行边坡稳定性系数求解时所采用的方法不同，前者通过求解边坡可靠度上限法数学规划模型，后者依据边坡处于极限状态的判断标准。此外，随着地下水位的升高，累积概率密度曲线逐步向左移动。图9(b)、图9(c)、图9(d)分别为tw=1、tw=25、tw=50三种情况下UBM 法所得边坡稳定性系数频数分布直方图，均呈现中间高两侧低的规律，与地下水位频数分布直方图相似。在tw=1时边坡稳定性系数在1.2 附近内出现的频数最多，最高频次为172 次；在tw=25时边坡稳定性系数在1.15 附近内出现的频数最多，最高频次为164 次；在tw=50时边坡稳定性系数在0.9 附近内出现的频数最多，最高频次为165 次。

4.5 边坡稳定性系数分布规律

根据边坡可靠度分析随机规划模型式(4)以及随机规划模型的求解策略，共得到50 个地下水位在1000 个随机场下的边坡稳定性系数。图10、图11分别为边坡稳定性系数的50 条概率密度曲线以及累积概率密度曲线；图12 为边坡稳定性系数的均值与地下水位的关系。通过分析可得到如下规律：

图10 边坡稳定性系数概率密度曲线Fig.10 Probability density curve of safety factor of slope

图11 边坡稳定性系数累积概率密度曲线Fig.11 Cumulative probability density curve of safety factor of slope

1)边坡稳定性系数的分布与正态分布一致，且随着地下水位的升高，边坡稳定性系数的均值逐渐降低，边坡稳定性系数的概率密度曲线和累积概率密度曲线均向左移动。此外，随着地下水位的升高，边坡稳定性系数的标准差逐渐降低，边坡稳定性系数的概率密度曲线分布范围逐渐变窄，累积概率密度曲线逐渐变陡。

2)边坡稳定性系数的均值与地下水位直接相关，采用多项式拟合得到边坡稳定性系数的均值与地下水位的量化关系式为：

式(27)表明边坡稳定性系数的均值与地下水位成负相关。

3)对50 个地下水位在1000 个随机场时所得的50 000 个边坡稳定性系数进行统计分析，得到考虑土体参数空间变异性以及地下水位随机性的边坡稳定性系数的分布特征，图13 为随机地下水位作用下边坡稳定性系数的频率分布直方图。随机地下水位作用下边坡稳定性系数的均值为1.127，标准差为0.097。

图13 边坡稳定性系数的频数分布直方图Fig.13 Frequency distribution histogram of safety factors

4.6 边坡系统失效概率分析

通过50000 次的模拟，LEM 法默认边坡仅存在一种失效模式(如图14 中的LEM 滑动面)，FEM法得到边坡的4 种失效模式(如图14 中的FEM 滑动面)。本文UBM 法相较于LEM 法、FEN 法最大的优势在于可以根据边坡中单元的失效信息采用AP 聚类分析方法捕捉到边坡所有的失效模式(如图14 中的失效区域)，进而进行边坡系统失效概率的估算，图14 为边坡的6 种失效模式AP 聚类结果，表4 列举了边坡的6 种失效模式以及对应的失效概率，表5 为tw=1、tw=25、tw=50时边坡失效模式以及对应的失效概率。由结果可知，失效模式1、失效模式2 失效区域面积较小，属于浅层滑坡，失效面积为5.45 m2～15.52 m2(失效次数为1193)；失效模式5、失效模式6 失效区域面积较大，属于深层滑坡，失效面积为89.45 m2～112.65 m2(失效次数为192)；失效模式3、失效模式4 介于浅层滑坡和深层滑坡之间，失效面积为48.72 m2～72.42 m2(失效次数为4430)。

表4 边坡失效模式以及对应的失效概率Table 4 Slope failure mode and the corresponding failure probability

表5 不同地下水位作用下边坡失效模式以及对应的失效概率Table 5 Slope failure modes and corresponding failure probabilities under the action of different groundwater levels

图14 边坡失效模式AP 聚类结果Fig.14 Schematic diagram of AP clustering results of slope failure mode

地下水位tw=1时，边坡仅具有单一失效模式，具体为失效模式3(失效次数为8，失效概率为0.8%)；地下水位tw=25时，边坡的失效模式变为两种，具体为失效模式3(失效次数为13，失效概率为1.3%)、失效模式4(失效次数为12，失效概率为1.2%)；地下水位tw=50时，边坡的6 种失效模式均可能发生。其主要原因是：随着地下水位的升高，边坡浸润线相应上移，边坡土体中饱和区域的面积增加，从而增加了失效模式4 的发生概率。此外，由于边坡较陡，随着地下水位的进一步升高，边坡土体中饱和区域的面积进一步增加，地下水从坡面溢出可能性增加，从而进一步增加了失效模式1、失效模式2、失效模式5、失效模式6 的发生概率。

LEM 法所得的破坏模式仅与本文结果的失效模式3 相吻合，其主要原因在于：采用LEM 法进行边坡稳定性计算时，事先假定初始滑裂面，并根据抗剪参数的均值搜索得到临界滑裂面，继而再基于唯一的临界滑裂面计算各个抗剪参数样本对应的边坡稳定性系数；因此，从理论上讲LEM法所得到的临界滑裂面与各个抗剪参数样本对应的真实滑裂面存在差距。FEM 法所得的破坏模式与本文结果的失效模式2、失效模式3、失效模式4、失效模式5 相吻合，并未获得失效模式1 和失效模式6，其主要原因在于：采用FEM 法进行边坡稳定性计算时，通过构建同时满足土体静力平衡条件、应变相容条件以及本构关系的有限元模型，可获得相较于LEM 法更为合理的结果；由于FEM 法在进行边坡稳定性计算时需要考虑材料的本构关系，在已知土体本构关系的基础上可获得相较于上限法更为精确的结果，但在进行边坡稳定性系数以及失效模式判断时包含较多的人为因素。此外，在计算效率上相较与上限法大幅降低。上限法通过有限单元构建边坡可靠度上限法数学规划模型，然后采用随机数学规划的方法搜索得到边坡的失稳区域，其结果更加接近于真实情况。计算结果表明：在多种随机参数作用下，传统的LEM 法会忽略边坡的多种失效模式，可能会错估边坡的失效风险；FEM 法在已知材料本构关系的情况下可以获得边坡所有的失效模式，但土体本构关系复杂，当前并没有统一的本构模型进行表征。本文上限法可以忽略材料的本构关系，获得边坡稳定性系数的上限解和破坏机构，其适用性和计算效率均较高。

需说明的是，本文同时计算了地下水随机数的数量nw取30、50、70 三种情况。当nw=30时，会损失部分失效模式(失效模式1 和失效模式2)，计算精度较低；当nw=50时，可同时保证计算精度和计算效率；当nw=70时，计算精度较nw=50时虽略微提高，但计算量大幅增加、计算效率大幅降低。限于篇幅，本文仅给出了nw=50的计算结果。

4.7 边坡单元失效概率分析

根据式(15)计算了50 个随机地下水位作用下边坡的整体失效概率。随着地下水位的逐渐升高，边坡整体失效概率从0.80%增加到92.70%，边坡的安全性能大幅度降低。通过分段3 次埃尔米特插值得到边坡整体失效概率与地下水位的关系曲线(如图15 所示)。地下水位在5.5 m～11.5 m范围内，整体失效概率变化幅度较小；地下水位在11.5 m～14 m 范围内，整体失效概率变化幅度较大；地下水位在14 m～15 m 范围内，整体失效概率变化幅度又变小。边坡的整体失效概率随地下水位的升高整体上呈现“平稳增加急剧上升平稳缓慢上升”的三段式特征。

图15 边坡整体失效概率与地下水位的关系Fig.15 The relationship between slope failure probability and groundwater

根据式(22)计算了50 个地下水位作用下边坡的单元失效概率，图16 为tw=1、tw=25、tw=50时的单元失效概率等值线；图17 为特征点的单元失效概率与地下水位的关系曲线。对以上图形进行分析，可以得到：

图16 单元失效概率等值线Fig.16 Element failure probability contour

图17 单元失效概率与地下水位的关系Fig.17 The relationship between the element failure probability and groundwater

1) 单元失效概率等值线可以非常直观的看出土质边坡每个部位的失效概率，地下水位tw=1、tw=25、tw=50时边坡单元失效概率等值线与该水位下边坡的失效模式轮廓线基本一致，但在失效可能性大小上存在明显差距，其主要原因在于AP聚类分析根据边坡中单元的失效信息将具有同一特征的失效面积聚为同一类，而单元失效概率等值线通过对已知单元的失效情况拟合得到。从理论上来说，如果边坡只有一种失效模式，那么边坡各单元的失效概率是相同的；如果边坡具有多种失效模式，则边坡各单元的失效概率会存在明显差异，单元失效概率等值线图也将存在差异；所有失效模式中重叠区域的单元的失效概率最大(如图16 所示)。

2) 监测单元的单元失效概率与地下水位的关系如图17 所示。在相同的地下水位作用下，不同的单元失效概率不同；在不同的地下水位作用下，同一单元的失效概率不同，不同位置处的单元失效概率随着地下水位的变化情况也存在明显差异。E2 单元和E3 单元随着地下水位的升高失效概率呈现增加的趋势，然而E1 单元和E4 单元随着地下水位的升高失效概率呈现先增大后减少的趋势，其主要原因在于所有失效模式均包含E2 单元，5 种失效模式(失效模式2、失效模式3、失效模式4、失效模式5、失效模式6)包含E1 单元，4 种失效模式(失效模式3、失效模式4、失效模式5、失效模式6)包含E3 单元，3 种失效模式(失效模式4、失效模式5、失效模式6)包含E4 单元，因此单元失效概率：E2≥E1≥E3≥E4。

5 结论

本文将有限元离散技术、上限法理论、相关非高斯随机场模拟以及随机规划理论结合起来，提出了随机地下水位作用下考虑参数空间变异性的边坡可靠度分析方法，并编制了相应的计算程序，获得了随机地下水位作用下考虑参数空间变异性的边坡稳定性系数、失效概率分布规律，为边坡可靠度分析提供了一种新方法。其主要结论如下：

(1)传统的LEM 法、FEM 法仅根据边坡稳定性系数来进行边坡稳定性的判断，该方法仅反映了边坡失效概率的程度。本文采用单元失效概率法对边坡的失效概率进行研究，不仅可以反映边坡失效概率的程度，而且可以根据单元的位置信息准确定位单元的失效情况。

(2)采用基于蒙特卡洛的迭代方法求解边坡可靠度分析的随机规划模型时需要进行成千上万次的模拟，本文基于Matlab 编制了高效的上限法并行程序，大大提高了计算效率。

(3)基于AP 聚类算法得到50 种水位下考虑土体参数变异性的边坡6 种失效模式以及对应的失效概率，得到边坡的系统失效概率为11.6180%，然而在多种随机参数作用下，传统的LEM 法会忽略边坡的多种失效模式，可能会错估边坡的失效风险，这是因为LEM 法在计算边坡临界滑裂面时其计算原理为：首先通过蒙特卡洛方法生成计算所需随机场，然后对每个随机场的抗剪强度参数进行平均化取值，再通过计算得到边坡的临界滑裂面。因此，LEM 法并没有得到边坡所有的失效模式。而UBM 法通过寻求构建机动许可速度场的最小值得到边坡的失效概率，所得结果更加符合真实情况。

(4)地下水位对考虑参数空间变异性的边坡整体失效概率和单元失效概率影响均较大，其中整体失效概率随着水位的升高不断增大，边坡单元失效概率等值线可便于岩土工程师更加直观的了解边坡各个位置的失效情况。