基于帕德逼近响应面法的边坡可靠度分析

梁 斌，高乐星，何 杰，吴 政，刘 杰

（湖南工业大学 土木工程学院，湖南 株洲 412007）

崩塌、泥石流、滑坡等边坡灾害会危害人们的生命及财产安全，因此，进行边坡稳定性分析和评估相当重要。由于边坡工程的破坏模式和岩土参数都存在不确定性和变异性，故常运用可靠度分析方法来评价边坡的稳定性，该方法已得到了广泛认可。常见的可靠度分析方法有直接积分法、矩分析法、响应面法、Monte-Carlo 法等[1-2]。直接积分法是可靠度分析中最基础的计算方法，但当计算结构失效概率时，其在实际工程应用中将遇到极大障碍。矩分析法在进行可靠指标求解时，通常需求解功能函数对基本变量的导数，这对高度非线性隐式函数来说非常困难。Monte-Carlo 法虽然能得到精确解，但计算量巨大，收敛速度较慢。响应面法的适用性较强，但是当功能函数的非线性程度较高时，使用传统的响应面法可能会出现精度及效率低的问题。因此，本文拟基于帕德逼近有理多项式提出一种新的代理模型响应面方法（response surface method based on Pade approximation，PRSM）。首先，采用拉丁超立方试验构造样本，并将样本数据代入简化Bishop 算法中，得到边坡功能函数值；然后，通过求解得到代理模型的待定系数，从而实现隐式功能函数显式化；最后，通过两个经典的边坡算例，验证该方法的可行性、计算的高效性及准确性，为边坡工程高度非线性隐式功能函数的可靠度求解提供一条切实可行的途径。

1 边坡稳定性分析功能函数

1.1 极限平衡法理论

极限平衡法、极限分析法以及滑移现场法[3]等都是常用的边坡稳定分析方法。目前国内外普遍使用的分析方法为极限平衡法。它以摩尔-库仑强度准则为基础，通过计算边坡已知滑移面的静力平衡，最终得到边坡稳定安全系数，以此评价边坡的稳定性。极限平衡法不仅逻辑严谨而且简单适用，至今已形成多种计算方法，如Bishop 法、简布法、摩根斯坦-普瑞斯法、萨尔玛法[3]等。极限平衡法的计算结果精度很高，能满足工程要求。所以，在诸多边坡稳定性分析方法中，极限平衡法的地位相当重要。

1.2 简化Bishop 法

目前，极限平衡法中的瑞典圆弧法（Fellennius法）和简化Bishop 法仍为边坡工程中两种常用的稳定性分析方法。相较于用瑞典圆弧法计算得到明显偏低的安全系数Fs，经过优化土条间作用力假设的简化Bishop 法更为合理，计算更简单。因此，本研究将采用极限平衡法中的简化Bishop 法建立边坡的稳定性功能函数。简化的Bishop 法大大改进了传统的瑞典圆弧方法[4]，除了给出安全系数Fs的定义之外，还忽略了条间剪力，并假定土条之间的力是水平的，则土条底面上的法向力由在竖直方向上的静力平衡条件求出。Fs借助边坡整体力矩平衡方程来确定。

简化Bishop 法中，滑体被划分为n个竖直条块，滑体及条块间作用力如图1所示。计算参数如下：R为滑面半径，Wi、bi、αi分别为第i个条块的重力、宽度、网弧底面倾角，ui、Ni、Ti分别为第i个条块的孔隙水压力、法向作用力、圆弧底面剪力，Ei、Ei+1、Xi、Xi+1分别为第i个与第i+1 个条块的横向作用力和纵向作用力，ci为黏聚力，φi为滑面内摩擦角。

图1 滑体及条块间作用力示意图Fig.1 Schematic diagram of the force between the sliding body and the bars

由太沙基有效应力原理和摩尔-库仑准则可得：

竖直方向上条块间力的平衡方程如下：

由式（1）和（2）可得：

式（3）（4）中：

仅考虑对圆心的力矩平衡，不考虑滑体条块水平方向力的平衡，则有：

将式（3）（4）代入式（6），得

忽略滑体条间剪力作用，则简化Bishop 法的安全系数Fs[5]的表达式如下：

建立边坡稳定性功能函数[6-7]，即

由式（8）（5）可知，Fs为隐式函数，同时也是关于ci、φi、Wi的函数，所以式（9）表示的边坡功能函数属于高度非线性隐式函数，此时采用传统的可靠度计算方法，如直接积分法，将会遇到极大困难。

2 基于帕德逼近的响应面法

2.1 抽样方法

边坡稳定性分析中，许多参数都存在随机性和变异性，故变量的样本数据选取需要采取试验设计抽样方法来实现，因为试验要想获得所有的参数样本非常困难。蒙特卡洛抽样、正交试验设计和拉丁超立方抽样都是常见的抽样方法。样本点的选择不仅要反映基本参数与功能函数之间的映射关系，还必须最大程度地代表整个变量空间的性质，并且试验次数还要尽量小。作为一种分层随机抽样方法，拉丁超立方抽样方法可用较少数量的样本反映总体的变化趋势，避免了重复抽样，有效减少了抽样次数。因此，本文参数样本的选取采用拉丁超立方抽样完成。

从理论上分析，样本点越多，计算结果精度越高，功能函数的模型代理表达式越精确，但相应的计算代价越大。根据目前已有研究成果[8]得知，模拟时一般选取30 组样本较为合适。

2.2 帕德逼近

帕德逼近是由法国数学家亨利·帕德发明的，它是一种非线性逼近方法，也是有理函数逼近中比较特殊的一种方法。与截断的泰勒级数相比，帕德逼近通常更加精确，并且帕德逼近在泰勒级数不收敛的情况下依旧可行，所以其在许多物理问题上甚至实际问题中都有较为普遍的应用。假设一个分式，它的分子为多项式，分母也为多项式，就等于一个无穷级数，如1/(1-x)=1+x+x2+···+xn+o(xn)，帕德逼近就是找到无穷项中收敛较慢的级数，然后将其变成有限多项式除法运算。以泰勒级数多项式为基础，帕德逼近构造了一类近似的有理多项式函数以替代原函数的方法。函数f(x)的有理逼近可通过其泰勒展开式获得[9]。

设f(x)∈CN+1(-a,a)，若有理函数

（式中Pn(x)与Qm(x)无公因式）满足：

则将Rnm(x)称为函数f(x)在x=0 处的(n,m)阶帕德逼近，表示为R(n,m)。

由此可得如下逼近公式：

大量研究结果[10-11]显示，当n=m，且n+m为常数时，帕德逼近具有最佳精度和最快收敛速度。

2.3 边坡功能函数模型代理

边坡的功能函数为一个高度非线性的隐式函数，首先考虑将其显式化后再进行可靠度计算，依据帕德逼近方法，可取如式（13）所示的功能函数近似显式表达式，其中n=m，边坡的各参数如ci、φi、γi（容重）可视为随机向量x=[x1,x2,…,xj]，可得近似的边坡功能函数Z的表达式为

式（14）即为边坡PRSM 代理模型。

然后，将边坡参数训练样本X=[x(1),x(2),…,x(m)]，容量为m，代入简化Bishop 功能函数式（9），经迭代后获得样本的真实响应值Y=Z=[Z(1),Z(2),…,Z(m)]；再将X、Y代入PRSM 代理模型式（14），求得Pk、Qk等待定系数，最后获得清晰的PRSM 边坡功能函数表达式。

2.4 可靠度计算步骤

利用PRSM 模型代理边坡的功能函数后，选取可靠度分析方法来计算可靠度与失效概率，本研究采用一次二阶矩法中的JC 法[1-2]编制Matlab 程序进行可靠度计算。

3 代理模型方法精度验证

3.1 实例1

一均质边坡[12]，其横截面如图2所示，高H=20 m，坡比为1∶1，其土层参数如下：泊松比v=0.3，弹性模量E=20 MPa，内摩擦角φ=20o，黏聚力c=40 kPa，剪胀角ψ=20o，容重γ=20 kN/m3。c、φ相互独立，均服从正态分布，变异系数都为0.1。

图2 实例1 边坡横截面示意图Fig.2 Schematic diagram of the slope cross section in example 1

边坡上部边界自由，底部边界为刚接，两侧为水平滑动支承，同时视边坡土体为摩尔-库仑理想弹塑性模型。上述土层参数中，容重γ变异较小，变异性相对较大的为内摩擦角φ和黏聚力c。因此容重γ的变异性不纳入考虑范围，视为定量参数进行分析，把c、φ视为随机参数。将其表示为随机向量x={c,φ}={x1,x2}的形式。通过拉丁超立方试验设计抽样得到30 组样本后，再由简化Bishop 法获取相应的功能函数值，其计算结果如表1所示。

表1 实例1 样本点及功能函数计算结果Table 1 Calculation results of sample points and function functions in example 1

在获得表1中的30 组数据和其相应的功能函数值之后，根据式（14）可求得其中的有理多项式系数Pk、Qk、Rk、Sk、Uk、Vk，从而构建了PRSM 边坡代理模型。然后，采用JC 法进行计算，选取均值点(40,20)作为验算点初值。取(1,1)阶帕德逼近响应面代理模型，经迭代后，求得可靠度指标β(1,1)=3.151 8，相应的失效概率Pf=8.114 1×10-4。取(2,2)阶帕德逼近响应面代理模型，求得可靠度指标β(2,2)=3.294 4，相应的失效概率Pf=4.931 5×10-4。

在文献[13]中，结合ANSYS 和响应面法所求解的可靠度为β=3.420，Monte-Carlo 法求解的可靠度为β=3.234，验算点法求解的可靠度为β=3.215，传统响应面法求解的可靠度为β=3.214。本文方法的计算结果与上述结果均相差很小，证明了本文方法是可行的。并且对比Monte-Carlo法求得的可靠度指标3.234，文献[13]中的方法误差为5.75%，本研究所用方法的误差仅为1.87%，可见本研究提出的方法拥有更高的精确度。在文献[13]中，ANSYS 的求解原理为强度折减，且程序一次只能求解一个结果，非常繁琐、缓慢。再者有理多项式的拟合精度比二次多项式的明显要高。在进行可靠度分析时，尤其是当随机变量的数量增加时，采用基于帕德逼近的响应面方法在效率与精度方面均优于上述其他方法。

3.2 实例2

图3为一双层边坡[14]，高15.24 m，坡比为1∶2。其土层参数如下：弹性模量E=100 MPa，泊松比v=0.3，黏聚力c1=38.31 kPa，c2=23.94 kPa，内摩擦角φ1=0o，φ2=12o，容重γ=20 kN/m3。c1、c2、φ2是互为独立的正态随机变量，变异系数分别为0.2,0.2,0.1，其余参数为定值。

图3 实例2 边坡横截面示意图Fig.3 Schematic diagram of the slope cross section in example 2

以上土层参数中，把c1、c2、φ2作为随机参数。同样30 组样本通过拉丁超立方抽样生成后，再由简化Bishop 法获取相应的功能函数值，其计算结果如表2所示。然后利用本文所提出的方法，分别采用基于(1,1)阶和(2,2)阶帕德逼近有理多项式的响应面方法进行边坡可靠度分析，所得计算结果见表3。

表2 实例2 样本点及功能函数计算结果Table 2 Example 2 sample points and function function calculation results

将表3中的数据与文献[14]和文献[15]中所计算的结果进行比较，可以发现可靠度指标结果比较接近，因而进一步证实了本文方法的可行性。但是本实例中，(1,1)阶PRSM 方法比(2,2)阶PRSM 方法的计算精确度更高。其原因可能是：当随机变量数量增加时，拉丁超立方抽样方法抽取的基本随机变量样本数量过少，在拟合PRSM 代理模型的多项式系数时没有取到最优值。

表3 不同计算方法的可靠度结果Table 3 Reliability results calculated by different methods

4 工程实例分析

某高速段路堑边坡坡高约为15 m，坡比为1:2。层状构造，岩体主要是炭质砂岩。边坡上部地层为红褐色粉质黏土，硬塑性，厚为0～1 m；中部地层为灰色弱风化碳质砂岩，厚度大于3 m；下部地层为青灰色全风化碳质砂岩，厚度大于10 m。进行土工试验后得知，该路段各地层容重γ变异较小，黏聚力c、内摩擦角φ的变异性较大，所以把c、φ作为随机参数，不考虑γ的影响。该工程实例的地层力学参数和统计特征如表4所示。

表4 工程实例地层力学参数及统计特征Table 4 Ground mechanics parameters with the statistical characteristics of engineering examples

把上部地层、中部地层、下部地层的物理力学参数c1、φ1、c2、φ2、c3、φ3写成随机向量的形式x={c1,φ1,c2,φ2,c3,φ3}={x1,x2,x3,x4,x5,x6}，所得计算结果如表5所示。

表5 工程实例样本点及功能函数计算结果Table 5 Project example sample points and function function calculation results

根据表5中的30 组样本和其功能函数值构建基于帕德逼近有理多项式的边坡功能函数的响应面模型；再利用JC 法，经计算，得到β=5.024 2，Pf= 2.528 1×10-7。计算结果表明，该路堑边坡出现失稳状况的可能性较低。

5 结论

1）本文采用极限平衡法中的简化Bishop 法，利用拉丁超立方试验设计构造样本，构建了基于帕德逼近有理多项式代理模型的边坡可靠度计算方法（PRSM），为边坡工程高度非线性隐式功能函数可靠度求解提供了一条新的途径。

2）PRSM 法针对传统响应面方法中的响应面函数进行改进，用基于帕德逼近有理多项式的形式替代传统的二次多项式形式，其计算结果精度更高，且当随机变量增加时，(2,2)阶PRSM 的精度比(1,1)阶PRSM 的更高。但也会出现(2,2)阶PRSM 的精度不及(1,1)阶PRSM 的情况，推测其原因为抽取的基本随机变量样本数量过少或没有在拟合PRSM 代理模型的多项式系数时取到最优值等，有待进一步研究。

3）通过实例验证，表明本文方法是可行的，且当计算精度接近时，其工作量小于蒙特卡洛法，具有一定的工程实用性。