均质边坡稳定性极限曲线法

方宏伟，李长洪，李 波

（1.辽宁交通高等专科学校 道桥系，沈阳 110122；2.北京科技大学 土木与环境工程学院，北京 100083）

1 引 言

边坡稳定性分析主要是条分法和有限元法，包括计算安全系数和搜索临界滑裂面两方面[1]。条分法局限性是需假定滑裂面的位置[2]，有限元法分为滑裂面应力分析法和强度折减法两类，前者为有限元与极限平衡理论的结合，仍需假定和搜索临界滑裂面，后者克服了前者的不足，但对边坡临界破坏的失效判据学术界尚无统一的意见[3]，而且不同程序计算的结果也不一致，因此有限元法还不能够替代极限平衡法[4]。为解决上述问题，学者们一方面努力寻找准确可靠的搜索方法[5]，另一方面不断提出新观点，例如，朱大勇等[6]认为直接求解最小安全系数是解决问题的突破口，并提出了边坡临界滑动场法。以上都是从力的角度进行研究的，陈震[7]认为可按边坡变形发展过程定量估定稳定性，提出可用破坏坡度与实际坡度之比作为安全系数，并指出联系变形发展情况还可能有其他表达形式。

滑移线法在边坡中的应用主要是计算无重边坡极限荷载[8]，近年来已有学者认为滑移线具有更重要的意义，并且提出与有限元法结合确定滑裂面[9]，张天宝等[10]认为该理论对土石坝合理边坡形状设计具有启发意义，并与高广岩等[11]通过有限元论证堆石坝合理边坡形状是凹形曲面的结论相一致。基于以上研究，以边坡变形量评价其稳定性为出发点，本文提出应用滑移线场理论计算得到的极限坡面曲线分析边坡稳定性的新方法，称极限曲线法。计算极限坡面曲线有2 种方法即Sokolovskii[7]研究得出的特征线法差分方程组和Cehkob 根据试验得到的极限稳定边坡曲线方程近似公式[10]，限于篇幅，本文不列出相关公式。对应的极限曲线法也分为称S曲线法（S curve method，SCM）和称C 曲线法（C curve method，CCM），采用Matlab 编程计算，笔者可提供源程序。通过经典考题与典型算例及34个边坡样本[2]的验算考察方法合理性和结果可靠性，进行因素敏感性分析，并应用于露天矿边坡稳定性和最终边坡角的分析与确定。

2 极限曲线法简介

2.1 基本概念与公式

滑移线法已计算得到了无重度边坡极限荷载，对于有重度边坡则要求坡面为凹形曲面才能求得解析解[8]，本文采用该计算逆过程。对于有重度边坡，在极限荷载作用下坡面形状应为凹形曲面，二维坐标下为凹形曲线，需要说明的是，即使边坡所受荷载不是极限荷载也可以求得极限坡面曲线，这条曲线是将现有荷载定义为极限荷载后边坡极限平衡状态下的坡面形状。将地基线以上边坡体放入第一象限，以坡脚为坐标原点，设坡高为H,F1(x)为极限坡面曲线拟合二次函数，当F1(x)与正x 轴的交点x11＞0时，将坡面线与极限坡面曲线之间的面积S1和坡面线与正x 轴所成面积S2之比定义为安全度（degree of safety，DOS），见图1。设坡脚到坡顶横坐标 x22，S2=x22H/2，极限坡面曲线与正x 轴所成面积，则S1=S2-S3，DOS=S1/S2，DOS 越大稳定性越好，值域为(0,1)，以极限坡面曲线与坡面线相交为变形破坏准则，其交点横坐标x1与x22之比的负值定义为破坏度DOF，见图2。此时，x11＜0，DOF=-x1/x22，DOF 越小，稳定性越差，值域为(-1,0)。由以上分析可知，x11为变形破坏准则判断值。

强度折减法构筑一个与真实边坡相同轮廓的“虚拟”边坡，强度指标缩减，缩减的系数即安全系数。极限曲线法正好是这个方法的对偶过程，即强度指标不变，但坡面缩减变形，按其变形量评价稳定性，因此不必假设和搜索临界滑动面。

图1 DOS 计算示意图Fig.1 Calculation sketch of DOS

图2 DOF 计算示意图Fig.2 Calculation sketch of DOF

2.2 算法函数与流程

极限曲线法函数：

式中：γ为重度(kN/m3)；c为黏聚力(kPa)；φ为内摩擦角(°)；α为坡度(°)；N1为主动区边界步长数；N2为过渡区边界步长数；Δx为SCM（S 曲线法）主动区边界步长(m)，CCM（C 曲线法）时将Δx改为 Δy=H/N1，其中N1为H 剖分数，无N2；f1为SCM/CCM，F0(x)为边坡坡面函数（本文为过原点的一次函数）；f2为章节2.1 求DOS、DOF 的方法。

采用准确率最高的3 次样条差值[12]计算x11，该法可以回避插值问题的Runge 现象，又是连续光滑的。为保证F1(x)与x 轴有交点，要求F1(x)纵坐标最小值ymin＜-1。对SCM 法，从理论上讲，应x11＞0时，x1＜0；x11＜0时，x1＞0，但由于F1(x)随着取点变化而发生多项式摆动[13]，会产生x11＜0时x1＜0的情况，此时DOS/DOF=0，当坡面无荷载时，满足计算条件的坡顶最小荷载[7]Pmin=ccot φ·(1+sinφ)/(1-sinφ)，此时无过渡区，则N2=0，而CCM法不存在上述2个问题，计算流程见图3。

图3 算法流程图Fig.3 Algorithm flow chart

3 考题和算例验算与样本分析

王国体[2]建立了边坡稳定性应力状态法，引用了其他文献的大量例题，为便于比较，本文例题皆取之于该文献，不同之处是用陈祖煜[1]的经典考题b 代替c。

3.1 经典考题与典型算例的验算

考题a 的参数为γ=20 kN/m3，c=3 kPa，φ=19.6°，α=27°，H=10 m。考题b 的参数为γ=20 kN/m3，c=32 kPa，φ=10°，α=27°，H=10 m。检验有限元强度折减法的典型算例参数为γ=20 kN/m3，c=42 kPa，φ=17°，α为变化范围为[30°，50°]，均分5个区间计算，H=20 m。

由于SCM 法采用特征线差分法，而拟线性双曲型微分方程组没有惟一解[14]，相应的数值方法也具有类似的特点，即SCM 计算的DOS、DOF 会随着边界条件不同而发生变化。为考察该法稳定性，采用不同的Δx 和N1及节点数k 计算分析，结果见表1、2。随着Δx 的减小和N1与k 的增加，计算精度随之增加，计算结果表明SCM 法稳定。计算范围影响分析见表3、4。当Δx 不变，N1增加时k 变多，对函数拟合有影响，DOS、DOF 会发生变化，但x11不变，说明解的变化完全是由于拟合振荡产生的，对边坡稳定性的判断不会产生影响。由于计算机性能的限制，N1不可能取得很大，本文计算N1≤999。CCM 法没有稳定性问题。

表1 经典考题SCM 法稳定性分析Table 1 Stability analysis of SCM

表2 典型算例SCM 法稳定性分析(不同坡度)Table 2 Stability analysis of SCM(different slope)

表3 经典考题SCM 法计算范围的影响分析Table 3 Calculation influence range of SCM

表4 典型算例SCM 法计算范围的影响(Δx=0.02 m)Table 4 Calculation influence range of SCM(Δx=0.02 m)

不同坡度SCM 法(Δx=0.02 m)计算图如图4 所示，其中的三角形为原边坡体，阴影区域为计算得到的极限稳定状态边坡体，CCM 法(Δy=0.001 m)计算图例如图5 所示，其中上侧边线包含的三角形为原边坡体，下侧边线包含的是极限稳定状态边坡体，随着边坡角增大，安全系数逐渐变小，边坡稳定性变差，DOS、DOF 也在变小，x11由大于0 变为小于0，即极限坡面曲线与x 轴的交点逐渐向左偏移；当x11＜0时，极限坡面曲线与坡面相交，且安全系数越小，x11越小，极限坡面曲线与坡面相交点x1越大，由此证明了变形破坏准则的正确性。

DOS、DOF 与安全系数对比见表5。安全系数较大时，SCM 法、CCM 法与其具有可比性；当安全系数变小时，SCM 法、CCM 法偏于保守，且CCM法保守程度大于SCM 法，原因应该是相对于坡顶无荷载的初始条件，SCM 法计算时在坡顶附加了最小荷载值，而CCM 法将坡顶以下部分边坡体自重作为外荷载。

3.2 样本计算与正确率分析

图4 不同坡度的SCM 计算图(典型算例)Fig.4 SCM diagram in different angles(typical example)

图5 不同坡度的CCM 计算图(典型算例)Fig.5 CCM diagram in different angles(typical example)

取文献[2]中（去掉c=0 和φ=0 的实例）共34个边坡样本并按原序列重新排序，计算结果见表6。分析可知，安全系数法计算正确率为67.7%，应力状态法为73.5%，CCM 法为79.4%，SCM 法为70.6%，SCM 法判断错误的样本安全系数法也都判断错误，CCM 法判断错误最少。

4 因素敏感性分析

按王国体[2]的作法，分析SCM 法、CCM 法在不同坡度下各因素变化对DOS、DOF 的敏感性，增加了坡高H 一项，计算结果分别见表7、8。安全系数法中各因素对边坡稳定性影响趋势是[15]：单一参数变化而其他参数固定，重度γ、边坡角α、坡高H 增加时，边坡稳定性降低，而黏聚力c 和摩擦角φ增加时，边坡稳定性增大。由各个表格数据可知SCM 法、CCM 法因素敏感性结论与安全系数法完全一致。

5 露天矿边坡工程应用

鞍钢集团鞍山矿业公司大连石灰石矿（原名大连甘井子石灰石矿）是鞍钢主要的辅料基地，在进行-50 m 深部开拓和南帮扩采2 项工程中，南帮扩采后最终边坡角等参数的确定是保证矿山安全生产的前提，因此边坡稳定性及参数研究成为矿山公司立项研究的重要内容[16]。通过工程类比分析，临界安全系数值确定为1.20～1.25，考虑爆破和地下水时安全系数为1.40～1.45。石灰岩岩体参数：重度γ=26.5 kN/m3，黏聚力c=225 kPa，内摩擦角为φ=36°，坡高H 的敏感度分析选用135 m，设最终边坡角为α0，其敏感度分析选用50°、55°、60°、65°，计算结果见表9，其中SCM 法边界参数为Δx=0.5 和N1=100，CCM 法边界参数 Δy=0.013 5 和N1=10 000，SCM 法和CCM 法α=65°时为DOF 情形，则最终边坡角应α0＜65°，这与原报告初步确定最终边坡角α0=55°～60°的结论一致。选择不同高度4个剖面进行稳定性分析，计算结果见表10，可知SCM 法、CCM 法与安全系数法结论一致，即在上述参数条件下边坡处于稳定状态。

当岩体强度参数取最小值 cmin=200 kPa 和φmin=33°时，原报告推荐最终边坡角 α0=55°～58°，且局部会有滑坡风险。4个剖面在边坡角为55°～60°时，安全系数法和SCM 法计算结果见表11，可知当α=60°时，H=128 m 剖面对应SCM法为DOF 情形，说明最终边坡角应为α0＜60°，SCM 法结论与原报告一致，即α0=55°～58°。由表11 分析可知，按临界安全系数值为1.40，则由安全系数确定的最终边坡角，除H=106m 剖面为α0=58°，其余均应α0＜55°，CCM 计算结果见表12。由表可知不同高度剖面具有不同α0，如H=128 m时，α0=55°，H=116 m时，α0=57°，H=106 m时，α0=59°，H=120 m时，α0=56°，可见CCM 法更有利于工程实践，且各剖面α0普遍提高了1°～2°。

表5 DOS 法、DOF 法与安全系数的对比Table 5 Comparison of DOS/DOF and safety factor

表6 样本分析与结果对比Table 6 Sample analysis and result comparison

表7 SCM 参数敏感性计算(不同坡度)Table 7 SCM parameter calculation(different slopes)

表8 CCM 参数敏感性计算(不同坡度)Table 8 CCM parameter calculation(different slopes)

表9 坡高(135 m)与最终边坡角敏感性分析Table 9 Sensitivity calculation of slope high(135 m)and ultimate pit slope

表10 不同剖面稳定性计算Table 10 Stability calculation of different profiles

表11 最小参数不同剖面SCM 稳定性计算Table 11 Different profile stability SCM calculation of the minimum parameters

表12 最小参数不同剖面稳定性CCM 法计算Table 12 Different profile stability CCM calculation of the minimum parameters

6 讨 论

对滑移线场理论，有学者认为其计算复杂，通用性差，在边坡工程中无实用价值。实际上，数值方法普遍具有复杂性，这就需要计算机软件和编程语言的支持[17]。本文采用Matlab 编程，尽管节点数量巨大，但时间成本很低，对于SCM 法，采用有限差分法，在N1=999时第一次运算40 min 左右，由于计算机储存功能二次计算N1不超过第一次则1 min 以内便可完成。CCM 法为试验方程近似公式，当N1=104时运算时间也只要1 min 左右。近40 多年来，学者们一直没有停止过对滑移线场理论的研究，在不同工程领域不断提出新观点，如在土压力中的应用[18]，在地基承载力中的应用[19-20]等。本文通过大量例题和露天矿边坡工程的计算数据证明了基于滑移线场理论的均质边坡稳定性极限曲线法的合理性和实用性。

7 结 论

（1）以边坡坡面变形量评价其稳定性为出发点，提出了基于滑移线场理论的均质边坡极限曲线法，由理论分析可知，该法是求有重边坡极限荷载逆过程，也是强度折减法对偶过程。根据计算极限坡面曲线方法的不同，将其分为基于特征线法差分方程组的SCM 算法和基于试验极限稳定边坡曲线方程近似公式的CCM 算法。以极限坡面曲线与坡面相交为变形破坏准则，定义了安全度（DOS）和破坏度（DOF）2个稳定性评价指标。

（2）由经典考题和典型算例的验算说明SCM法稳定，典型算例的计算数据与图例证明了变形破坏准则的正确性，以上2个例题的计算结果表明，安全系数较大时，SCM 法和CCM 法与其具有可比性，由于计算时增加了外荷载，故安全系数变小时，SCM 法和CCM 法偏于保守。34个样本正确率：安全系数法为67.7%，应力状态法73.5%，CCM 法为79.4%，SCM 法为70.6%，表明该法正确率较高，计算结果可靠。

（3）SCM 法和CCM 法因素敏感性分析结论与安全系数法一致。应用于露天矿边坡稳定性和最终边坡角的分析与确定，两种计算结果与原报告结论一致，当参数变小时，CCM 法更有利于工程实践，表明该法具有一定的工程实用价值。

（4）SCM 法和CCM 法采用Matlab 编程，时间成本低廉，且不必假设和搜索临界滑动面。对c=0和φ=0 的情况尚需研究并编入程序，对非均质情形，笔者已完成了成层土质边坡稳定性极限曲线法的研究工作，对更加一般的非均质情形尚需进一步研究。本文的极限曲线法具有探索性，需要更多的工程实践检验。

[1]陈祖煜.土质边坡稳定性分析——原理方法程序[M].北京:中国水利水电出版社,2003.

[2]王国体.边坡稳定和滑坡分析应力状态法[M].北京:科学出版社,2012.

[3]陈力华,靳晓光.有限元强度折减法中边坡三种失效判据的适用性研究[J].土木工程学报,2012,45(9):136-146.CHEN Li-hua,JIN Xiao-guang.Study of the applicability of three criteria for slope instability using finite element strength reduction method[J].China Civil Engineering Journal,2012,45(9):136-146.

[4]CHENG Y M,LANSIVAARA T,WEI W B.Reply to“Comments on ‘Two-dimensional slope stability analysis by limit equilibrium and strength reduction methods’ by Y.M.Cheng,T.Lansivaara and W.B.Wei,”by J.Bojorque,G..De Roeck and J.Maertens[J].Computers and Geotechnics,2008,35(2):309-311.

[5]CHENG Y M.Location of critical failure surface and some further studies on slope stability analysis[J].Computers and Geotechnics,2003,30:255-267.

[6]朱大勇,张四维.边坡稳定性计算新方法[C]//第五届全国工程地质大会文集.河南:中国地质学会工程地质专业委员会,1996.

[7]陈震.散体极限平衡理论基础[M].北京:水利电力出版社,1987.

[8]陈祖煜,汪小刚,杨健.岩质边坡稳定性分析-原理方法程序[M].北京:中国水利水电出版社,2003.

[9]朱以文,吴春秋,蔡元奇.基于滑移线场理论的边坡滑裂面确定方法[J].岩石力学与工程学报,2005,24(15):2609-2616.ZHU Yi-wen,WU Chun-qiu,CAI Yuan-qi.Determination of slip surface in slope based on theory of slip line field[J].Chinese Journal of Rock Mechanics and Engineering,2005,24(15):2609-2616.

[10]张天宝,张立勇,敖天其.土坝和土石坝合理边坡设计研究[J].水力发电学报,1985,3(3):28-38.ZHANG Tian-bao,ZHANG Li-yong,AO Tian-qi.Study on the design of reasonable slope of earth and earthrockfill dams[J].Journal of Hydroelectric Engineering,1985,3(3):28-38.

[11]高广岩,张天宝.高堆石坝合理边坡形状的静力有限元分析[J].四川大学学报(工程科学版),2002,34(1):28-31.GAO Guang-yan,ZHANG Tian-bao.Static finite element analysis of reasonable slope shape of high rockfill dams[J].Journal of Sichuan University(Engineering Science Edition),2002,34(1):28-31.

[12]马莉.MATLAB 数学实验与建模[M].北京:清华大学出版社,2010.

[13]MATHEWS J H,FINK K D.数值方法(MATLAB版)(第4 版)[M].周璐,陈渝,钱方等译.北京:电子工业出版社,2005.

[14]胡祖炽,雷功炎.偏微分方程初值问题差分方法[M].北京:北京大学出版社,1988.

[15]廖珊珊,张玉成,胡海英.边坡稳定性影响因素的探讨[J].广东水利水电,2011,40(7):31-40.LIAO Shan-shan,ZHANG Yu-cheng,HU Hai-ying.Factors affecting the slope stability study[J].Guangdong Water Conservancy and Hydropower,2011,40(7):31-40.

[16]鞍钢集团鞍山矿业公司研究所.鞍矿大连石灰石矿南帮扩建边坡稳定性及参数研究[R].鞍山:鞍钢集团鞍山矿业公司研究所,2006.

[17]张鲁渝,欧阳小秀,郑颖人.国内岩土边坡稳定分析软件面临的问题及几点思考[J].岩石力学与工程学报,2003,22(1):166-169.ZHANG Lu-yu,OUYANG Xiao-xiou,ZHENG Ying-ren.Problems and thoughts of development of slop stability analysis software in China[J].Chinese Journal of Rock Mechanics and Engineering,2003,22(1):166-169.

[18]彭明祥.挡土墙被动土压力的滑移线解[J].岩土工程学报,2011,33(3):460-469.PENG Ming-xiang.Slip-line solution to passive earth pressure on retaining walls[J].Chinese Journal of Geotechnical Engineering,2011,33(3):460-469.

[19]SMITH C C.Complete limiting stress solutions for the bearing capacity of strip footings on a Mohr-Coulomb soil[J].Geotechnique,2005,55(8):607-612.

[20]MARTIN C M.New software for rigorous bearing capacity calculations[C]//Proceedings of the British Geotechnical Association International Conference on Foundations.Dundee:Thomas Telford,2003:581-592.