考虑网格自适应的边坡可靠性随机有限元极限分析研究

许晓亮 张家富 曾林风 徐健文 史为政

(三峡库区地质灾害教育部重点实验室(三峡大学)， 湖北 宜昌 443002)

边坡作为常见的岩体之一，具有显著的不确定性，其中一个主要的表现便是岩土参数的空间变异性[1]．为了更好地考虑岩土体参数空间变异特性，部分学者利用随机场[2]来进行刻画，也有学者通过机器学习[3]与贝叶斯压缩感知[4]等方法利用有限现场实测数据对参数进下概率反演，从而更好地表征岩土参数变异性．岩土参数空间变异性与可靠性分析结合受到越来越学者的关注和探究[5-6]．

边坡确定性分析方法主要包括极限平衡(LE)、极限分析法(LA)和有限元法(FEM)等．LE法需要提前假设一个破坏机制，需人为地假定一系列平衡条件，同时没有考虑岩土体的应力应变关系，故理论严密性不足，而LA法将真实解限制在下限解[7]和上限解[8]之间，可信度要高于LE法，但是对于复杂的岩土工程问题，传统LA法很难人为地构建静止许可应力场和机动许可速度场，这导致了其应用受到较大限制．直到有限元和极限分析法相结合形成了有限元极限分析法[9](finite element limit analysis，FELA)，该方法既克服了LA法在复杂工程下无法人为构建应力场与速度场问题，也避免了LE法理论严密性不足等问题．

如前所述，将随机场理论与可靠性分析相结合，是考虑岩土体参数空间变异性影响的一个重要途径．在边坡工程的研究中，Cho[10]将整个边坡区域生成随机场，并进行极限平衡分析，但当参数空间变异性较大时，得到的滑移面不尽合理；Griffiths等[11]较早地将FEM与随机场有效结合，研究了空间变异性对边坡稳定性的影响．同时，基于FELA的优势，将FELA与随机场相结合的研究也自然受到了重视．Krabbenhoft等[12]将强度折减法应用于FELA，并结合随机场理论得到了基于FELA的随机场方法，实现了考虑参数变异性下的边坡FELA强度折减稳定性计算；Huang等[13]将有限元极限分析应用于滑坡问题，通过Monte-Carlo模拟建立了一种新的滑坡概率分析框架；陈朝晖等[14]基于FELA分析了参数空间变异性对边坡稳定性的影响，并与经典极限平衡法进行比较，发现通过传统极限平衡法计算安全系数再进行可靠性分析，会高估边坡失效概率．上述研究表明，除在边坡稳定性的确定性分析外，FELA在边坡可靠性分析中也逐渐得到了应用．然而，在FELA计算中，为了得到更为接近真解的上下解，需要增大单元数量，压缩上下限区间至极小状态，进而要开展大量的网格剖分，导致计算规模大幅增加，因此部分学者[15]通过网格自适应加密，通过对塑性区等应变率较大区域继续加密处理，从而提高计算效率得到更小的上下限区间．目前，考虑网格自适应的FELA法在边坡确定性分析中已经取得了一定的进展[16-17]，但在考虑空间变异性的随机场条件下，采用网格自适应的FELA对边坡稳定性及可靠性分析结果的影响并不清楚．

为此，本文结合随机场理论和FELA方法，利用Matlab编制了基于Karhunen-Loeve随机场离散的二维随机场程序以及Matlab与极限平衡、有限元软件交互程序，提出了基于网格自适应有限元极限分析边坡可靠性分析方法，针对一边坡算例，将网格自适应前后的FELA法以及传统极限平衡法进行对比分析，并基于网格自适FELA边坡可靠性分析结果，探讨了岩土参数变异系数和波动范围对安全系数均值、可靠度指标以及边坡潜在滑移面等因素的影响．

1 随机场理论

1.1 土体空间变异性

为了更真实地表征土体参数的空间特征，需在对边坡进行稳定性分析时考虑参数空间变异性．根据随机场理论，随机场可通过参数均值、变异系数、波动范围来表征空间变异性，为生成满足要求的平稳随机场，采用被广泛使用的指数型自相关函数来表示土体的空间变异性[18]，表达式如下：

(1)

式中：(x，y)分别为两点空间坐标(x1，y1)与(x2，y2)；δh和δv分别为水平波动范围与竖向波动范围，波动范围越大，表明土性参数分布越平稳．

1.2 随机场离散

随机场模型一般不能在有限元分析中直接使用，需将随机场离散从而与有限元相结合．随机场离散的方法有协方差法、局部平均细分、Karhunen-Loeve(K-L)级数展开法等，本文采用K-L展开法[19]生成满足要求的平稳正态随机场如下：

(2)

(3)

由于岩土体物理参数往往不能取负值[20]，需要将标准正态随机场转化成对数正态、Beta和截尾指数等分布的非高斯随机场来表征参数的空间特征．考虑到岩土体参数大多服从正态和对数正态分布[20]，本文采用对数正态随机场，将任一点的参数X作为随机变量并服从对数正态分布，其K-L离散结果表达式如下：

(4)

式中：ulnXi和σlnXi分别为相应高斯参数随机场中lnXi均值和方差，由对数正态变换式(5)～(6)计算得：

(5)

(6)

上述二维随机场进行离散时，对于可分离型自相关函数，二维的特征值、特征函数可简化为一维方向上的特征值乘积、特征函数乘积[21]；

λi=λxi·λyi

fi(x，y)=fi(x)·fi(y)

(7)

一维方向上特征值与特征函数求解步骤如下，以x方向为例(y方向同理)：

1)设x区域随机场范围为Ω=[xmin，xmax]，用[-1，1]区间的标准Legendre多项式函数Lj-1(x)表示[xmin，xmax]区间的正交基函数：

(8)

其中：a和b分别为缩放参数和平移参数，a=(xmax-xmin)/2，b=(xmax+xmin)/2．

2)分别由式(7)生成x方向矩阵Ax(Ay同理)．

(9)

式中：ρ(x，y)为前述指数型自相关函数；通过积分换元以及参数的缩放平移可进一步改进为：

Lj-1(u)Lk-1(v)dudv

(10)

3)利用矩阵方程求解两个方向上特征值λx(λy对应Dy)

AxDx=DxΛ

(11)

式中：Dx为矩阵Ax特征矩阵；Λ为对角矩阵，将其主对角元素从大到小排列即为λxi．

4)利用正交基函数hi(x)[22]求解特征函数fi(x)(fi(y)对应hk(y))

(12)

上述特征值以及特征函数的大小均取决于展开项数M，M越大离散越精确，为了控制离散精度和计算量，引入期望能比率因子ε[22]≥0.95，来控制离散误差从而确定合适的展开项数M．

(13)

式中：S为二维离散区域面积；ε为期望能比率因子．

2 网格自适应有限元极限分析

2.1 有限元极限分析原理

有限元极限分析方法[12-13]基于塑性极限理论，结合广义原理和混合有限元理论．有限元极限分析通过有限单元将岩土体内应力场和机动场离散化，从而在离散后的场内构建满足条件的约束方程，并将岩土体内总的能量耗散(上限分析)或外力荷载(下限分析)作为上、下限目标函数，最终通过数学规划的手段求解目标函数，从而得到上、下限解，基本思想与方法如下：

对于弹塑性材料，其广义变分原理如下式：

s.t.F(α，κ)≤0

(14)

式中：α为荷载因子；σ为任一点应力状态；b为重力矢量；u为位移矢量；t代表边界应力．F(α，κ)为屈服函数，当κ=0即为理想弹塑性材料．

FELA下限解(LB)假定应力场内满足平衡条件，应力边界条件和屈服条件且理想弹塑性材料不破坏，此时的应力场即为静力许可应力场，对应的外部荷载即为破坏荷载，下限解即满足静力许可条件的最大荷载，对应的目标函数优化问题可如下表述：

maxα

PTσ=αtonsσ

f(σ)≤0

(15)

FELA上限解(UB)假定机动场满足应变-位移关系、关联流动法则和位移边界条件，在弹塑性材料破坏时，内力功大于等于外力功．则对应的目标函数优化问题可表述为：

(16)

2.2 网格自适应加密

在上述优化问题求解过程中，为了缩小上下解区间，需增加网格单元数，若对网格均匀划分的话，要达到满足的精度需要极大的计算规模，但对塑性区外的网格进行划分并不会提高精度，因此采用网格自适应加密，只在塑性区等应变率较大的区域进行加密，计算规模较小并能提高计算精度．

文中上限自适应采用一种基于能量耗散的自适应加密方法，即计算出每个单元的能量耗散，然后通过反复加密能量耗散较大的三角单元来降低破坏区域内单元间能量耗散的差异，从而在破坏区形成能量耗散密集的单元，而在其他区域单元呈疏松状态，具体加密过程见文献[16]．

下限自适采用一种基于节点应力的自适应加密策略[15]，将Mohr-Coulomb准则改写成一个等式约束和一个二阶锥约束，将前面优化问题转化为锥优化的形式求解，最终通过衡量屈服准则残余以及等效变形来指导网格进行剖分，从而使破坏区域网格更密集．

3 考虑网格自适应的边坡可靠性随机有限元极限分析方法

3.1 参数随机场模拟

通过前述K-L离散方法生成平稳随机场，为验证所生成随机场的合理性，以黏聚力为例，生成不同均值度下不同波动范围随机场，如图1所示．

图1 不同波动范围下的随机场

由图1中(a)、(b)可知，当δh/δv=1时，所模拟的黏聚力呈块状的分布，且随着波动范围的增大，块状分布更大，低黏聚力以及高黏聚力分布越集中；由图中(c)、(d)可知，当δh/δv=10时，由于水平波动范围较大，竖向波动范围较小，不同波动范围下的随机场中，黏聚力分布呈水平层状，且随着波动范围的增大，水平层状带更长，这一规律与已有研究一致[23]．

3.2 考虑网格自适应的边坡可靠性FELA流程

本文通过Matlab生成N组对数正态参数随机场，结合Monte-Carlo法，将每组参数随机场作为随机变量样本，结合OPTUM-G2对边坡进行N次确定性分析并结合强度折减[12]得到N组边坡上、下限安全系数，代入式(17)所表征的功能函数g(X)，并由式(18)可计算得出边坡失效概率Pf与可靠度指标β．

g(X)=Fs(X)-1

(17)

β=Φ-1(1-Pf)

(18)

实现流程如图2所示，该过程由4个步骤组成：

Step1：建立边坡计算模型，读取单元中心坐标．首先在OPTUM-G2建立边坡模型并进行网格划分，导出对应单元信息．

Step2：生成参数随机场．结合岩土参数统计信息，选取自相关函数与波动范围，利用Matlab编制对数正态随机场程序，基于Step1中单元节点信息生成N组岩土参数随机场．

Step3：批量确定性分析．将上述随机场参数批量导入原始边坡模型中，生成N组边坡计算文件，利用Matlab调用OPTUM-G2进行批量网格自适应FELA确定性分析并结合强度折减计算上、下限安全系数．

Step4：计算失效概率及可靠度指标．统计Step3中N组确定性分析计算结果，结合式(18)计算出失效概率与可靠度指标．

图2 基于网格自适应FELA边坡可靠性分析流程图

4 算例分析

以一简单土坡为算例，如图3所示，坡高为8 m，边坡顶部宽度为4 m，临空坡面的坡比为1∶2．边坡模型均匀划分为1 000个三角单元，单元尺寸为0.46 m，单元尺寸与竖直波动范围之比[24]为0.153<0.25，满足要求，边坡土体采用使用摩尔-库仑准则和非关联流动法则的理想弹塑性本构模型，计算模型底部采取水平与竖直约束，左侧采取水平约束，右侧无约束；模型加自重荷载．

图3 边坡模型图

确定性分析时，不考虑参数空间变异特征，黏聚力c、内摩擦角φ、弹性模量E、重度γ、泊松比ν等参数采用均值见表1．

表1 土性参数取值表

在可靠性分析时，仅考虑c和φ的空间变异特征，其余参数视为定值，弹性模量、泊松比、剪胀角等参数对边坡稳定性状态分析结果的影响可忽略不计[25]．算例中，c和φ均服从对数正态分布，变异系数、波动范围等参数参考文献[26-27]选取，初步设立随机场参数取值如下，δh=30 m，δv=3 m．

4.1 确定性分析结果

当不考虑空间变异特征进行确定性分析时，计算参数采用均值，见表1．借助GEO-SLOPE采用Morgenstern-Price(M-P)法得出的安全系数为1.412．由FELA计算的塑性乘子云图如图4(a)、(c)所示，上限解(UB)安全系数为1.453，下限解(LB)为1.393；而采用网格自适应3次迭代后FELA计算云图如图4(b)、(d)所示，上限安全系数为1.419，下限为1.395．需要指出的是，随着迭代次数的增大，计算所需时间逐渐增加，但经试算发现，当迭代次数大于3以后，迭代次数的增加带来的上下限区间缩小幅度逐渐减小，为提高计算效率，文中采取3次迭代的计算结果．

图4 网格自适应前后FELA法确定性分析图

通过塑性乘子云图对比分析，自适应后的等效塑性区潜在滑移面网格分布较未自适应更加密集，网格加密区所反映的破坏形态更为清晰，滑动带明显更窄；且FELA法在网格自适应后计算的上限安全系数较小，下限也有所增大，上下限的区间范围明显更小；若要增加单元数量要达到相同效果，由表2可知，以上限解分析为例，需要增加单元至8 000个，计算时间也将由网格自适应所需的8.91 s增加至19.65 s．

表2 时间对比

可见，若不考虑网格自适应，在一次计算中需通过增加单元数来缩小上下限区间，其用时更长，而对于后续可靠性分析，成千上万次的计算将大幅增加时长，而考虑网格自适应则能在保障上、下限解逼近真解的同时显著提高计算效率．

4.2 随机性分析结果

考虑土体参数的空间变异特征进行可靠性分析时，参数取值见表1，采用网格自适应的FELA方法，利用Monte Carlo法模拟2 000次．结果表明，土体参数变异性会导致边坡滑移面发生显著变化，同一安全系数往往会对应多个滑移面形式，对应于4.1节中自适应下限解为1.395的情况，在可靠性分析中共出现11个滑面(如图5所示)，滑移面与图4(b)中确定性分析潜在滑移面有较大区别，部分滑移面不再经过坡脚而经过坡面，滑移面积差异性也较大，其中最小滑移面积28.43 m2，最大滑移面面积64.87 m2，这表明参数空间变异性下，坡体潜在滑移面位置和规模也呈现出较为显著的不确定性．

图5 同一安全系数滑移面分布

为了进一步研究参数随机场分布与潜在滑面的关系，以黏聚力c的随机场为例，图6(a)、图6(b)分别给出了安全系数最小与最大时对应的潜在滑移面以及随机场，图6(c)、图6(d)为滑移面积最小与最大时对应的潜在滑移面和随机场．可看出潜在滑移面更靠近且穿过土体参数较小的区域，且基于M-P法得出的边坡安全系数及潜在滑移面积均大于网格自适应FELA的结果．

图6 典型样本滑移面及随机场

对比图7中M-P法、FELA法和自适应后的FELA所得的安全系数概率密度曲线可知，采用FELA法，不论其上限或下限解，在考虑网格自适应和未考虑时其安全系数变异性均比较接近，且都明显小于M-P法的计算结果．同时，FELA法所得安全系数概率密度曲线均比M-P法“高瘦”，无论是否采用网格自适应，FELA法所得失效概率均小于M-P法的失效概率，故由传统极限平衡法进行可靠性分析往往会高估边坡的失效概率．

图7 概率密度曲线规律对比

分别对FELA法和自适应后的FELA生成的2 000次安全系数进行统计分析，图8给出了相应的累计分布曲线和频率分布直方图．由图8可知，自适应前后得到的边坡安全系数均服从对数正态分布，自适应后，安全系数上限均值从1.408减小至1.378，下限均值从1.331增大至1.338，上下限均值区间明显减小；同时，上限可靠度指标从1.903减小至1.805，下限可靠度指标从1.594增加至1.630，上下限更为趋近，由此计算出来的可靠度指标明显更接近与真实状态．且由图8中的概率密度曲线可知，无论是否采用网格自适应，2 000次FELA法计算出来的安全系数中有约65%以上都低于确定性分析安全系数，这表明在该边坡算例中，忽略岩土体参数空间变异性将会得出较高的安全系数，造成对边坡稳定性的不保守估计．

图8 网格自适应前后FELA法安全系数统计规律

4.3 随机场参数敏感性分析

变异系数与波动范围是表征随机场的重要参数，为研究其对计算结果的敏感性，仍采用网格自适应FELA方法，基于已有研究[26-27]设计计算参数，见表3．计算中，改变黏聚力变异系数ccov，内摩擦角变异系数φcov，水平波动范围δh和竖向波动范围δv，研究随机场参数对边坡的安全系数均值、可靠度指标、安全系数概率密度曲线以及滑移面分布范围的影响规律．

表3 敏感性分析参数取值表

4.3.1 变异系数影响

为研究土体参数变异系数对边坡的影响，取表3中组合1、2，分别研究c、φ变异系数的影响规律．以组合1为例，φcov=0.2，ccov取值范围为[0.2,0.8]，波动范围取值为δh=30 m和δv=3 m．

图9和图10分别为不同变异系数ccov、φcov下的可靠性指标β和安全系数均值uF值，结果表明β和uF有着相同的变化趋势，均随着变异系数的增大而减小．当ccov从0.2增加到0.8时，相应的uF上、下限解降幅分别达到8.1%和10%；β上、下限解降幅分别为68.6%和72.4%；而当φcov从0.05增加到0.25时，相应的uF上、下限解降幅为3.1%和3.9%；β上、下限解降幅则为45.9%和45.4%．对比数据，β比uF变化幅度明显更大，说明相对于安全系数，可靠度指标更为敏感，这与类似研究的结论一致[27]．

图9 β和uF随ccov变化曲线

图10 β和uF随φcov变化曲线

以下限法为例，图11和图12分别给出了安全系数概率密度曲线随ccov和φcov的变化趋势．

图11 不同ccov下安全系数概率密度曲线

图12 不同φcov下安全系数概率密度曲线

不难看出，随着参数变异系数增大，坡体物理力学性质的离散性变大，薄弱区更集中，更容易出现局部强弱化区，在概率密度曲线上整体呈现出向左偏移的同时逐渐扁平趋于“矮胖”状的特点；相应的安全系数均值下降但变异性增大，坡体失效概率逐渐增大，可靠度指标减小．对比ccov和φcov对概率密度曲线的影响规律可见，两者虽然有相同趋势，但后者相较于前者的影响更为显著，即参数φ的变异性更为敏感．

仍以下限法为例，图13和14为不同ccov和φcov下的坡体潜在滑移面分布范围，图中同一颜色的两条滑面分别为某一变异系数下计算得到的滑面上下界限．当ccov和φcov增加时，浅层滑面位置上移，最小滑移面面积减小，深层滑面下移，最大滑移面面积增大，即坡体潜在滑面的分布范围更大，滑移面分布更为离散．

图13 不同ccov下滑移面分布图

图14 不同φcov下滑移面分布图

4.3.2 波动范围影响

为研究随机场参数波动范围对边坡可靠性的影响，取表3中参数组合3～5进行分析．其中组合3，控制竖向波动范围为3 m，水平波动范围取值范围[10 m，50 m]，主要研究水平向波动范围的影响；组合4中，控制水平波动范围为30 m，竖向波动范围取值范围[1 m，5 m]，主要研究竖向波动范围影响；由于一般水平波动范围比竖向波动范围高一个数量级，为研究两者共同影响，控制δh/δv=10，即表3中组合5．图15和图16反映了随机场波动范围δh、δv对β和uF的影响，与前述变异系数的影响有所不同，随着波动范围的增大，β呈现减小而uF呈现增大趋势．当δh从10 m增加50 m时，相应的uF上、下限解增幅分别为0.9%和0.6%，β上下限降幅分别为11.9%和13.4%；当δv从0.05增加到0.25时，uF上、下限解增幅分别为3.0%和5.7%；而β上、下限解降幅分别为27.9%和28.5%．

图15 β和uF随δh变化曲线

图16 β和uF随δv变化曲线

从图15～16可以看出，相对安全系数，可靠度指标对波动范围更敏感，且竖向波动距离δv比水平向波动距离δh对安全系数和可靠性指标影响更大；此外，值得注意的是，随着波动范围增大，坡体安全系数均值变大的同时，其可靠性指标减小，这与一般地“越安全，越可靠”的认知相斥．

为了探究上述安全系数与可靠性指标变化趋势互异的原因，图17和图18给出了不同δh和δv下的安全系数概率密度曲线．

图17 不同δv下安全系数概率密度曲线

图18 不同δh下安全系数概率密度曲线

可以看出，随着波动范围的增大，安全系数概率密度曲线整体呈现出向右偏移且趋于“矮胖”状的特点，且随竖向波动范围δv的变化中更为明显．概率密度曲线整体右移将使得安全系数均值有所增大，趋于“矮胖”状意味着安全系数变异性增大，进而坡体失效概率增大，相应的可靠度指标减小．由于波动范围表示随机场中任意两点力学参数不在相关的临界距离，不难理解，波动范围越大，表明土性参数分布越平稳，其变异程度减弱，可进一步解释安全系数均值有所增大的原因．

此外，结合图15～16可以发现，无论δv还是δh，当增加到一定值时，安全系数均值与可靠度指标变化趋于平缓，受波动范围影响渐渐减弱．这是由于波动范围较小时，边坡土体参数可能更容易出现局部弱强度区，当波动范围增大，弱强度区范围也逐渐增大，但是增加到一定值时，土体参数整体分布更趋于均值，对坡体可靠性指标等影响不再显著．

图19给出了δh和δv同步增大时2 000次蒙特卡罗模拟中下限解对应的潜在滑移面分布范围．与坡体抗剪强度参数变异系数增大引起滑移面分布更为离散相反，当波动范围整体增大时，土性参数分布越平稳，其变异程度呈减弱趋势，最小滑移面与最大滑移面之间区间范围减小，滑移面分布范围越来越集中．

图19 不同波动范围下滑移面分布图

5 结 论

提出了考虑网格自适应的边坡可靠性随机有限元极限分析方法，通过均值边坡算例，探讨了网格自适应前后的FELA法的确定性和可靠性分析结果，同时在网格自适应FELA可靠性分析的基础上，从安全系数均值、可靠度指标、安全系数分布规律、潜在滑移面分布范围等角度探讨了参数变异系数和波动范围的影响，主要得出以下结论：

1)同等计算精度下，若不考虑网格自适应，需通过大幅增加单元数来缩小FELA的上下限区间，其用时更长，而考虑网格自适应的FELA上下限更为接近，潜在滑面位置更加明确，在更易于逼近真解的同时可显著提高计算效率，尤其对于需大量抽样计算的可靠性分析，其优势更突出．

2)对比基于M-P和网格自适应的FELA的随机性分析结果，前者更易得出更为深层的潜在滑面和偏大的失效概率；如忽略岩土体参数空间变异性，采用确定性分析则会得出更大的安全系数，进而引起对边坡稳定性的过高估计．

3)参数变异系数对边坡可靠性有显著影响，随着ccov和φcov的增大，安全系数概率密度曲线在整体向左偏移的同时逐渐略趋于扁平，安全系数均值下降但变异性增大，可靠度指标减小，同时，坡体潜在滑面的分布范围更大，更为离散；相对于黏聚力c，可靠性分析结果及潜在滑面位置对内摩擦角φ的变异性更为敏感．

4)波动范围越大，安全系数概率密度曲线在整体向右偏移的同时逐渐略趋于扁平，坡体安全系数均值与变异性均增大，进而可靠性指标减小，且竖向波动范围的影响更显著；但当波动范围增加到一定程度时，波动范围对可靠度指标的影响也渐渐减弱，主要是由于波动范围越大土性参数分布越平稳，其变异程度减弱，相应得出的潜在滑移面分布则更加集中．