基于DEM-SRM和盲数理论的边坡地震稳定可靠性分析

刘焕玉，魏汝明，胡文奎，程爱华，陈训龙

(1.济南市勘察测绘研究院，山东 济南 250000； 2.华中科技大学岩土与地下工程研究所，湖北 武汉 430074)

1 引 言

我国是一个山脉众多和地震频发的国家，因地震而产生的滑坡灾害给国家和人民的财产和人身安全造成了巨大的损失[1]。因此，对边坡工程的抗震分析和研究具有重要的意义。

现有的边坡抗震分析主要是基于拟静力法和动力时程分析法[2～4]。拟静力法的基本思想[5]是在静力计算的基础上，将地震作用简化为一个惯性力系附加在研究对象上，其核心是设计地震加速度的确定问题。由于拟静力法的计算简单，工程物理意义明确，且计算结果偏于安全，已成为边坡工程中最常用的计算分析方法[6]。动力时程分析法是考虑地震的动态过程性，将地震力以地震加速度时程的方式施加于研究对象上，以此来分析研究对象的动力响应及稳定性，其核心问题是边坡模型确定及地震输入的准确性[7，8]。随着计算机技术及数值分析方法的发展，动力时程分析法可以实现复杂边坡模型的建立和地震动力响应的快速计算和结果输出。因此，对于复杂边坡模型的地震稳定性分析，动力时程分析法是一个优选方法[9]。然而，现有基于以上两种分析方法的边坡地震稳定可靠性研究往往仅考虑了边坡岩土体参数的一种或两种不确定性，即随机性和模糊性[10，11]。事实上，边坡工程存有诸多不确定性因素，使得边坡岩土体参数也具有多种不确定性，如随机性、模糊性、灰色性和未确知性等。单独仅考虑其中的一种或两种不确定性易造成对边坡稳定可靠程度的高估，从而影响边坡的抗震分析[12]。盲数理论可以同时考虑参数的多种不确定性，但目前采用盲数理论分析边坡稳定性问题的研究还很少[13，14]，而将盲数理论应用于边坡抗震稳定可靠性分析的研究更是未见报道。

因此，本文基于盲数理论进行边坡的地震稳定可靠性分析，并提出了一种分析方法。首先，基于离散元强度折减法和拟静力法计算分析地震作用下边坡的安全稳定状态，并结合响应面法和蒙特卡洛法建立地震边坡的极限状态方程，其次，运用盲数理论对边坡岩土体参数进行盲数化处理，并进一步建立地震边坡的安全系数盲数模型，最后，通过MATLAB编写相应的计算程序实现安全系数盲数解的快速输出和分析。

2 方法与基本原理

2.1 离散元强度折减法基本原理

离散元强度折减法的基本思想是在采用离散元法分析边坡稳定性时，采用强度折减法计算出边坡的安全系数。

(1)离散元法基本原理

离散元法[15](Distinct Element Method，简称DEM)是20世纪70年代由Cundall提出的一种处理非连续介质问题的数值模拟方法，其离散单元能更真实地表达岩土体的几何特征，便于处理所有非线性变形等岩土体破坏问题，因此，被广泛应用于模拟边坡滑坡等大变形的力学过程。该方法的理论基础是结合不同本构关系的牛顿第二定律，采用动态松弛法求解方程，其基本方程由块体运动方程和物理方程组成。

在数值计算分析中，将块体离散化为有限差分三角形，网点由三角形的顶点组成，各网点的运动方程如式(1)表示：

(1)

物理方程表示块体间接触点的力和位移的关系。各块体之间的接触可归纳为角-边接触和边-边接触。力与位移的关系采用增量表示，具体如式(2)和式(3)所示：

△Fn=Kn△un

(2)

△Fs=Ks△us

(3)

其中：△Fn和△Fs分别为法向力增量和切向力增量；△un和△us分别为其对应的法向位移增量和切向位移增量；Kn和Ks分别为接触面上的法向刚度和切向刚度。

根据牛顿第二定律及中心差分法，可进一步得到：

(4)

由此，可计算出每一时步的变形、振动及其对应的变形速度、转动速度等参数。

(2)强度折减法基本原理

强度折减法[16]又称为抗剪强度折减系数法，其基本原理为：对边坡岩土体的抗剪强度指标，即黏聚力C和内摩擦角φ，同时除以一个折减系数K，如式(5)和(6)所示的形式进行折减，然后将这组折减后的值C′和φ′作为新的参数代入进行计算。当计算刚好不收敛时，此时坡体被认为达到极限稳定状态，对应的值称为边坡的安全系数。

C′=C/K

(5)

φ′=arctan(tanφ/K)

(6)

2.2 盲数理论

(1)盲数[12]的概念

(2)盲数的运算

设*为运算符号+，-，×，÷中的一种，A与B为2个盲数，其表达式分别为：

(7)

(8)

则A和B间的运算法则A*B可通过可能值和可信度值两个矩阵来实现。需要注意的是，在进行盲数运算时，A*B中的相同元素须进行合并，当成一个元素重新进行排列，同时，在可信度矩阵中须将相同元素对应的可信度值相加，所得之和即为这个元素的最终可信度值。

(3)盲数可信度

设A与B均为盲数，A*B在r处的可信度可表示为P(A*B≥r)，其表达式为：

(9)

其中：r为具体问题中要求值或标准值。

3 地震边坡稳定可靠性分析的盲数模型

3.1 地震边坡极限状态方程的建立

以一单层均质土坡为例，考虑地震力作用，选用拟静力法施加地震作用力，水平地震作用力Fe的表达式如式(10)所示：

Fe=Wαhξαi/g

(10)

其中：W为岩土体重力；αh为水平向设计地震加速度系数；ξ为地震作用效应的折减系数；αi为岩土体质心处水平地震加速度分布系数，可根据图1内插算得[17]，其中αm为分布系数最大值，当地震设计烈度为Ⅶ、Ⅷ、Ⅸ时，对应αm分别取3.0、2.5和2.0；g为重力加速度，一般取 9.8 m/s2。

图1 水平地震加速度分布系数参考取值

选取重度γ、黏聚力C和内摩擦角φ作为基本变量，则根据响应面法原理，建立地震作用下边坡极限状态的响应面函数为：

g(γ，C，φ)=Fs-1=α+ϖ1γ+ϖ2C+ϖ3φ+χ1γ2+χ2C2+χ3φ2

(11)

其中：Fs为地震边坡安全系数；α，ϖ1，ϖ2，ϖ3，χ1，χ2，χ3为响应面的待定系数，需要7个方程确定求解，其求解步骤如下[12]：

(1)假定初始迭代点Y(1)={γ(1)，C(1)，φ(1)}，一般选取均值点。

(2)采用拟静力法和离散元强度折减法，分别计算出地震边坡的极限状态函数值g{γ(1)，C(1)，φ(1)}、g{γ(1)±fσγ，C(1)，φ(1)}、g{γ(1)，C(1)±fσC，φ(1)}及g{γ(1)，C(1)，φ(1)±fσφ}，其中f为响应面参数，一般选取固定值为2，σγ，σC，σφ分别为对应参数γ，C，φ的标准差。

(3)将步骤(2)得到的迭代点和函数值代入式(11)中，计算得到7个待定系数，从而确定响应面函数。

(4)采用蒙特卡洛法，计算步骤(3)中得到的极限状态响应面函数的可靠度指标βk、可靠度Pk和验算点Y*(k)。

(5)计算验算点Y*(k)的极限状态函数值g{γ*(k)，C*(k)，φ*(k)}，并利用下式(12)计算下一步迭代点Y(k+1)。

(12)

(6)根据步骤(5)得到的迭代点Y(k+1)，重复步骤(2)～(5)进行下一步迭代计算，如满足下式(13)的条件，可认为计算收敛，则结束迭代计算，输出此时的可靠度指标βk以及可靠度Pk，此时得到的极限状态响应面函数即为所求。

|βk-βk-1|<0.001

(13)

3.2 地震边坡安全系数盲数模型的建立

(14)

其中：γmin为γ的下限，γmax为γ的上限，f1(γ)为γ的盲数；Cmin为C的下限，Cmax为C的上限，f2(C)为C的盲数；φmin为φ的下限，φmax为φ的上限，f3(φ)为φ的盲数。

由此，根据式(11)，易建立并得到地震边坡安全系数的盲数模型为：

Fs=α+1+ϖ1{[γmin，γmax]，f1(γ)}+ϖ2{[Cmin，Cmax]，f2(C)}+ϖ3{[φmin，φmax]，f3(φ)}

(15)

3.3 地震边坡安全系数盲数模型的快速计算

根据已知相关参数值及建立的盲数模型，通过MATLAB编写相应的计算程序，可以实现地震边坡安全系数盲数解的快速输出和分析，程序的具体计算步骤如下：

(1)根据式(15)，计算地震边坡安全系数盲数解，输出地震边坡的总体可信度概率分布图，同时，将计算结果以盲数标准形式输出，如式(16)所示的三段式盲数解，实际可根据具体需要划分区间段数：

(16)

其中：Fmin，F1，F2，Fmax为安全系数；P1，P2，P3为可信度。

(2)根据步骤(1)中盲数解以及前述响应面法算得的地震边坡可靠度Pk，进一步计算出边坡盲数安全系数Fps[12]。

(3)根据《建筑边坡工程技术规范》(GB50330-2013)[18]中规定的边坡安全最低安全系数值及式(16)，快速计算出边坡最小安全要求可靠度Ps[12]。

4 算例分析

4.1 算例概况及模型建立

选用文献[12]中的算例进行验证和分析，已知该边坡为一级边坡，边坡岩土体的物理力学参数如表1所示，其中γ，C，φ服从正态分布。该边坡的场区抗震设防烈度为Ⅶ度，地震的相关参数如下：αh=0.10g；ξ=0.25；αm=3.0。

边坡岩土体物理力学参数 表1

根据该边坡的剖面尺寸，采用离散元法建立边坡的三维模型如图2所示，其中Z轴方向长度取为 1 m。边界条件采用左右两侧水平约束，底部固定，上部自由。为简化计算，固定Z轴约束和变形，将边坡模型分析转化为平面应变问题，采用莫尔-库伦屈服准则的弹塑性条件进行计算和分析。

图2 边坡的离散元三维数值模型

4.2 盲数计算与分析

根据已知条件，采用本文方法首先计算出地震作用下边坡的极限状态方程，同时，得到响应面法迭代计算最终的可靠度Pk为58.32%。然后，进一步求解出地震边坡安全系数盲数解，具体结果如表2所示。

地震边坡安全系数盲数解 表2

表2给出了地震边坡安全系数取值区间及对应的可信度大小，其中可信度可直接理解为可靠度。由此，可以计算出任一安全系数值所对应的可靠度值，从而可以更全面地分析地震边坡的安全稳定状态，并根据实际工程及抗震要求来进行边坡抗震设计与加固处理。根据表2的结果，可进一步采用MATLAB软件快速计算并输出地震边坡的总体可信度概率分布，结果如图3所示。

图3 地震边坡总体可信度概率分布图

从图3中可以很直观地看出各安全系数值所对应的地震边坡可信度的大小。当安全系数取值为0.432时，对应地震边坡的可信度最大，此时可信度值为1，可认为边坡在地震作用下取此安全系数值最保险；而当安全系数取值为1.714时，对应地震边坡的可信度最小，可信度值为0，可认为边坡在地震作用下取此安全系数风险最大，风险概率为100%。

4.3 地震边坡稳定可靠性分析

根据表2和图3的结果，可进一步计算地震作用下边坡的盲数安全系数Fps和最小安全要求可靠度值Ps。根据响应面法迭代计算最终得到的可靠度Pk=58.32%，可计算出地震边坡盲数安全系数Fps=1.073。已知该边坡为一级边坡，其对应的边坡稳定安全系数要求值为1.30，由此，可计算出地震边坡最小安全要求可靠度值Ps=33.08%。由Fps=1.073<1.30，可知该边坡未达到一级边坡稳定安全系数的最低要求值，是偏不安全的，同时，由Ps=33.08%，可知该边坡在地震作用下具有较高的失效风险。

未验证本文方法的可行性和有效性，与现有的一般方法进行了对比，各方法的计算结果如表3所示，其中方法(1)为采用拟静力法考虑地震作用，M-P法计算均值安全系数以及蒙特卡洛法计算可靠度；方法(2)为采用拟静力法考虑地震作用，有限元强度折减法计算均值安全系数以及响应面法计算可靠度；方法(3)为不考虑地震作用，基于盲数理论计算的盲数安全系数和最小安全要求可靠度[12]。

边坡稳定可靠性计算结果比较 表3

从表3中可以看出：①在考虑地震作用的情况下，本文方法是有效可行的，且计算结果较现有方法(1)和(2)偏低，由于本文方法考虑了边坡岩土体参数的多种不确定性，分析因素更为全面，从而使得计算结果更为可靠且更具参考价值，对工程而言也更安全；②本文方法与方法(3)对比可知，地震作用对边坡的整体稳定可靠性具有明显影响，这与现有已知结论一致，因此，在边坡工程分析、加固处理及抗震设计时，应考虑地震作用。

5 结 论

本文采用拟静力法考虑地震作用，并采用离散元强度折减法、响应面法和盲数理论提出了一种边坡地震稳定可靠性分析方法。结合具体算例分析，得出如下结论：

(1)本文方法考虑了地震作用和边坡岩土体参数的多种不确定性，计算结果更为可靠且更具参考价值，对工程而言也更安全。

(2)根据本文方法计算出的地震边坡安全系数盲数解和边坡总体可信度概率分布，可以计算出任一安全系数值所对应的可靠度值，可以更全面地分析地震边坡的安全稳定状态，并根据实际工程等级情况及抗震要求来进行边坡抗震设计与加固方案设计。

(3)地震作用对边坡的整体稳定性具有明显影响，在边坡工程分析中应考虑地震的因素，同时，本文方法也为边坡地震可靠度分析提供了一个新的研究思路和参考。