均质边坡稳定状况快速辨识及参数反演的一种实用方法
张科，曹平
(中南大学 资源与安全工程学院，湖南 长沙，410083)
摘要：基于塑性极限分析的运动单元法(kinematical element method)提出一种边坡稳定状况快速辨识方法。研究岩土体参数组合(坡高H、黏聚力c、内摩擦角
和重度γ)与安全系数F和滑动面位置间的内在联系。在此基础上，绘制边坡稳定分析图，求解不需要任何迭代，可快速求解安全系数和滑动面位置。根据边坡稳定分析图，提出一种新的参数反演方法。结果表明，无量纲参数
控制着滑动面位置和
。λ不变，则滑动面位置和保持不变；λ越大，边坡失稳模式由浅层破坏变为深层破坏，滑动面坡顶滑出点离坡肩越来越远。当λ＞0.2时，λ与符合线性关系，可用拟合公式快速计算安全系数。根据原始边坡形态和滑动面坡顶滑出点离坡肩的距离，即可快速反算滑带岩土体的黏聚力和内摩擦角。
关键词：边坡稳定；运动单元法；临界滑动面；稳定分析图；参数反演
中图分类号：TU 443 文献标志码：A 文章编号：1672-7207(2013)07-2855-06
A practical method for fast identifying stability state and back calculating parameters in homogeneous slopes
ZHANG Ke, CAO Ping
(School of Resources and Safety Engineering, Central South University, Changsha 410083, China)
Abstract: A new fast identification of slope stability state based on KEM (kinematical element method) was presented. The internal relationship between combination of rock and soil parameters (slope height H, cohesion c, friction angle and unit weight γ) and stability state (factor of safety F and location of critical slip surface) was studied. The slope stability charts were developed, and a rapid and reliable way was provided to calculate the factor of safety and the location of critical slip surface without iteration. Based on the above slope stability charts, a new parameters inversion method was presented. The result shows that the location of critical slip surface and are determined by dimensionless parameter . If λ remains constant, the location of critical slip surface and remain the same. As λ increases, the failure mode changes from shallow slip to deep slip, and the horizontal distance from the toe of the slope to the failure surface at crest is much farther. When λ＞0.2, the relationship between λ and is linear, and the factor of safety can be calculated fast by the fitting formulas. With the slope geometry and the horizontal distance from the toe of the slope to the failure surface at crest in a specific cross-section, strength parameters of landslide slip can be obtained through back-calculation.
Key words: slope stability; KEM; critical slip surface; slope stability chart; parameter inversion
边坡稳定性一直是岩土工程的重要研究内容之一，其关键问题是确定边坡的最小安全系数及相应的最危险滑动面。边坡稳定计算属于超静定问题，无法直接由静力平衡条件得出边坡的安全系数。为了使问题静定可解，极限平衡法对条间力作出适当的假定，可分为瑞典法、简化Bishop法、Janbu法、Spencer法、Morgenstern-Price法和Sarma法等[1]。这些假定带来一定的计算误差难以估计，直接影响结果的精度。Gussmann[2]创立了塑性极限分析的运动单元法(kinematical element method)，该方法具有严格的数学力学基础。曹平等[3]进一步发展运动单元法，并开发了相应的计算软件。罗晓辉等[4]应用运动单元法，结合遗传算法提出最危险滑动面与地基极限承载力的求解方法。Gussmann[5]应用运动单元法求解有限荷载作用下的边坡稳定性。李峰等[6]采用运动单元法进行基坑坑底稳定性分析。滑带岩土体强度参数的取值正确与否，直接影响边坡稳定评价和工程设计的结果[7]。目前，抗剪强度参数的确定方法主要有工程类比法、试验测试和反分析[8-10]。工程类比法是一种经验估算方法，由于滑坡成因、几何边界、岩土体性质和研究者的工程经验存在着差异性，工程类比法难以准确得出抗剪强度参数。室内试验的试样尺寸小，取样时易对岩土体造成扰动，试验结果往往具有较高的离散性；而原位试验所需经费多，试验周期长。已发生的滑坡是非常难得的边坡现场原位试验，采用反分析方法求解抗剪强度参数更为经济、简便和可靠。反分析有2种方法[7-9]：(1) 对某个剖面进行反算。前人研究认为，采用该方法求解c和
2个未知数，其解具有不确定性，所以需要假定其中一个容易确定其取值范围的变量，结合敏感性分析，反算另一个变量。(2) 对多个剖面联立方程进行反算，选取的剖面必须相似。本文作者基于塑性极限分析的运动单元法提出一种边坡稳定状况快速辨识方法，分析抗剪强度参数与安全系数和滑动面位置间的内在规律，绘制了边坡稳定分析图。在此基础上，提出了一种新的参数反演方法。
1 运动单元法基本原理
运动单元法的理论体系包括单元运动分析、单元静力分析和多变量目标函数极值的优化求解。滑移边界上正应力σ和剪应力τ满足Mohr-Coulomb准则。运动单元法将近似假设的塑性滑动区离散为有限个单元，称为运动单元。运动单元法边坡稳定分析模型如图1所示。其中：L为坡趾到坡顶滑出点的水平距离；D为坡趾到滑动面最低点的深度。
根据单元运动协调条件，建立单元运动方程组：
(1)
式中：[K]为运动方程系数矩阵；{V}为各单元未知位移的列向量；
为已知位移的列向量。
求解上述方程组得各单元的位移，就可以确定各单元间的相对运动方向，求得作用在单元边界上切向力的作用方向。
将作用在单元边界上的法向力作为未知量，建立单元静力平衡方程组：
(2)
式中：[Ks]为静力方程组系数矩阵；{N}为各边界上未知法向力的列向量；{F}为已知力向量，与边界的体力和黏聚力有关。
求解静力方程组后，迭代求出中间安全系数。将安全系数作为目标函数，将滑动面节点坐标作为变量，利用多变量优化求极值方法搜索临界滑动面，并确定相应的最小安全系数。
图1 运动单元法边坡稳定分析模型
Fig.1 Slope stability analyses model on KEM
2 边坡稳定状况快速辨识方法
2.1 基本原理
在极限平衡法中，某一滑动面的安全系数被定义为此滑动面上岩土体抗剪强度指标与剪应力的比值[11-13]。一般采用摩尔库仑准则来表征岩土体的剪切强度，所以安全系数的定义为：
(3)
式中：
，
；c和分别为边坡岩土体的黏结力和内摩擦角；
和
分别为边坡处于极限平衡状态下的黏结力和内摩擦角。
为研究岩土体参数组合与安全系数和滑动面位置间的内在联系，设某边坡黏聚力为c，内摩擦角为，最小安全系数为F；将式(3)进行变换，得
。另有一边坡黏聚力为
，内摩擦角为
，其他参数都相等，设最小安全系数为F'。若令
，得
，即
，说明这2个边坡在极限平衡状态下计算参数相等。代入式(3)可得，F'=wF，从而得
，
。由此可见，对于相同几何形状和重度的边坡，
控制着临界滑动面位置、
和
。
众多研究者整理了坡高H、坡角β、与岩土体参数c，和γ之间的关系，以此为基础绘制了大量的边坡稳定分析图，一般都需要迭代步骤[14]。而Bell[15]提出的方法是最便捷的，不需要任何迭代。Bell认为为λ的函数，定义无量纲参数λ如式(4)所示。
(4)
为进一步研究岩土体参数对临界滑动面位置和安全系数的影响规律，选取文献[15]中的无量纲参数λ，将文献[16]的匀质边坡作为分析对象，坡高H=10 m，坡角β=45°，岩土体参数为：重度γ=20 kN/m3，黏聚力c=12.38 kPa，内摩擦角=20°。
改变黏聚力和内摩擦角，令折减系数
，研究λ为定值情况下，对边坡安全系数和滑动面位置的影响，计算结果如表1和图2所示。由表1和图2可知，λ不变，则滑动面位置和保持不变。
分别改变λ，研究λ对边坡安全系数和滑动面位置的影响，计算结果如表2和图3所示。由表2和图3可知，无量纲参数λ控制着滑动面位置和。λ越大，越大，滑动面越深，失稳模式由浅层破坏变为深层破坏，滑动面坡顶滑出点离坡肩的距离也越来越远。
表1 无量纲参数λ和安全系数F的关系
Table 1 Relationship between dimensionless parameter λ and factor of safety F
2.2 边坡稳定分析图
对于特定几何形状的边坡，临界滑动面的位置只与λ有关，因此，在本文绘制的边坡稳定分析图中定义无量纲参数D/H和L/H，用于表征临界滑动面的位置。
图2 无量纲参数λ和临界滑动面位置的关系
Fig.2 Relationship between dimensionless parameter λ and location of critical slip surface
表2 无量纲参数λ和安全系数F的关系
Table 2 Relationship between dimensionless parameter λ and factor of safety F
图3 无量纲参数λ和滑动面位置的关系
Fig.3 Relationship between dimensionless parameter λ and location of critical slip surface
坡角分别选取15°，30°，45°，60°，75°，90°，稳定分析图如图4所示。对于非特定的坡角，可以通过线性插值的方法来求解安全系数和滑动面位置。
假设某匀质边坡，边坡高度为8 m，坡角为30°；重度、黏聚力和内摩擦角分别为18 kN/m3，20 kPa和15°。计算
。由图4可得：=6.5，D/H=0.15，L/H=2.4。求得F=1.74，D=1.20 m，L=19.2 m。
图4 边坡稳定分析图
Fig.4 Stability charts for slope
从图4可以看出，当λ＜0.2时，λ和间呈非线性关系。当λ≥0.2时，λ与呈线性关系，采用最小二乘法拟合，拟合公式为
，也可以写成
。
拟合结果见表3，说明两者符合线性关系。对于λ＞0.2的边坡，可以用表3的拟合结果快速计算安全系数。
选取本节算例，
＞0.2，采用拟合公式计算安全系数。将表3中的系数代入拟合公式，得
。
表3 λ-F拟合结果
Table 3 Fitting result of λ and F
3 抗剪强度参数反演新方法
3.1 基本原理
边坡滑带岩土体抗剪强度参数反分析就是根据滑坡体的几何尺寸，结合现场勘察得出的边坡安全系数，建立数学模型，然后利用该模型反算岩土体的c和。安全系数建议根据边坡变形状态参考表4[17]选取。
研究表明，对于特定几何形状的边坡，λ控制着临界滑动面位置和。通过现场勘查得到滑体形状参数(滑动面坡顶滑出点离坡肩或坡趾的距离比较容易确定，故本文选取该参数)，根据图4(c)中
关系曲线就可以找到对应的无量纲参数λ，其临界滑动面与真实滑动面吻合。然后根据图4(a)中
关系曲线以及待分析边坡的安全系数，即可求出内摩擦角，再应用式(4)求出黏聚力。对于非特定的坡角，可以通过线性插值求解c和。
滑带岩土体抗剪强度参数反演方法流程见图5。
表4 滑坡不同阶段的安全系数[17]
Table 4 Stability coefficient for different development stages of landslide[17]
图5 强度参数反算流程
Fig.5 Flow chart of strength parameters inversion method
3.2 算例验证
选取2.2节中的算例，坡高H=8 m，坡角β=30°，重度γ=18 kN/m3，将安全系数和滑动面位置视为已知，即F=1.74，L/H=2.4；而假定抗剪强度参数(c和)未知。
根据L/H=2.4，由图4(c)可得λ=0.52，由图4(a)可得=6.5，已知F=1.74，可得
；根据式(4)，求得c=20.05 kPa。反演所得的结果与精确解几乎一致。
3.3 实例分析
某露天矿山边坡7'与9'勘探线之间发生了多台阶中型滑坡，塌滑体垂直高度H约112 m，原坡角约42°，滑动面后缘距坡肩距离
约10 m，即L'/H=0.09，塌滑体重度为18 kN/m3。
当坡角β=30°，根据L/H=
+0.09=1.82，读取图4(c)可得λ=0.03，读取图4(a)可得=2.3，取F=1.0，可得
；根据式(4)，求得c=27.03 kPa。
当坡角β=45°，根据L/H=
+0.09=1.09，由图4(c)可得λ=0.02，由图4(a)可得=1.4，取F=1.0，可得
；根据式(4)，求得c=29.60 kPa。
通过线性插值可得坡角β=42°边坡的综合抗剪强度参数为：
，c=29.09 kPa。
4 结论
(1) 提出了一种基于运动单元法的边坡稳定状况快速辨识方法，绘制了边坡稳定分析图。
(2) 无量纲参数λ控制着滑动面位置和。λ不变，则滑动面位置和保持不变。λ越大，边坡失稳模式由浅层破坏变为深层破坏，滑动面上缘离坡肩越来越远。当λ＞0.2时，λ与符合线性关系，可用拟合公式快速计算安全系数。
(3) 基于边坡稳定分析图，提出一种新的抗剪强度参数反演方法。根据某个剖面上的原始边坡形态和滑动面上缘离坡肩的距离，可快速反算滑带岩土体的黏聚力和内摩擦角。
参考文献：
[1] 陈祖煜. 土质边坡稳定分析: 原理·方法·程序[M]. 北京: 中国水利水电出版社, 2003: 18-19.
CHEN Zuyu. Stability analysis of soil slopes: Theory, methods and programs[M]. Beijing: China Water Power Press, 2003: 18-19.
[2] Gussmann P. Kinematical elements for soils and rocks[C]// Proceeding of 4th International Conference on Numerical Methods in Geomechanics. Edmonton, 1982: 47-52.
[3] 曹平, Gussmann P. 运动单元法与边坡稳定性分析[J]. 岩石力学与工程学报, 1999, 18(6): 663-666.
CAO Ping, Gussmann P. Kinematical element method and slope stability analysis[J]. Chinese Journal of Rock Mechanics and Engineering, 1999, 18(6): 663-666.
[4] 罗晓辉. 用遗传算法确定地基承载力[J]. 岩石力学与工程学报, 2001, 20(3): 394-398.
LUO Xiaohui. Determination of subsoil bearing capacity using genetic algorithm[J]. Chinese Journal of Rock Mechanics and Engineering, 2001, 20(3): 394-398.
[5] Gussmann P. Effective KEM solutions for the limit load and the slope stability problem[J]. International Journal for Numerical and Analytical Methods in Geomechanics, 2000, 24(14): 1061-1077.
[6] 李峰, 王晓睿, 罗晓辉, 等. 基坑坑底稳定性的机会约束评估方法[J]. 岩土力学, 2010, 31(12): 3867-3874.
LI Feng, WANG Xiaorui, LUO Xiaohui, et al. Assessment methods of chance constrained on bottom stability of foundation pit[J]. Rock and Soil Mechanics, 2010, 31(12): 3867-3874.
[7] 郑明新. 论滑带土强度特征及强度参数的反算法[J]. 岩土力学, 2003, 24(4): 528-532.
ZHENG Mingxin. Research of strength characteristic of landslide slip and revised counter calculation method[J]. Rock and Soil Mechanics, 2003, 24(4): 528-532.
[8] 高德军, 徐卫亚, 郭其达. 长江三峡大石板滑坡计算参数反分析[J]. 河海大学学报: 自然科学版, 2006, 34(1): 74-78.
GAO Dejun, XU Weiya, GUO Qida. Back analysis of calculation parameters for Dashiban landslide at Three Gorges Reservoir on Yangtze River[J]. Journal of Hohai University: Natural Sciences, 2006, 34(1): 74-78.
[9] 高丽娜, 严明. 反算法在某边坡稳定性评价中的应用[J]. 水土保持研究, 2006, 13(2): 75-77.
GAO Lina, YAN Ming. The application of the back calculation method in the evaluation of the slope stability[J]. Research of Soil and Water Conservation, 2006, 13(2): 75-77.
[10] 陈骏峰. 降雨型堆积层滑坡抗剪强度参数反演分析[J]. 华中科技大学学报: 城市科学版, 2008, 25(4): 249-252.
CHEN Junfeng. Back analysis of the shear strength parameters of rainfall-induced colluvial landslides[J]. Journal of Huazhong University of Science and Technology: Urban Science Edition, 2008, 25(4): 249-252.
[11] Michalowski R L. Stability charts for uniform slopes[J]. Journal of Geotechnical and Geoenvironmental Engineering, 2002, 128(4): 351-355.
[12] Jiang J C, Yamagami T. Charts for estimating strength parameters from slips in homogeneous slopes[J]. Computers and Geotechnics 2006; 33: 294-304.
[13] 郑颖人, 陈祖煜, 王恭先, 等. 边坡与滑坡工程治理[M]. 北京: 人民交通出版社, 2007: 95-96.
ZHEN Yingren, CHEN Zuyu, WANG gongxian, et al. Engineering treatment of slope and landslide[M]. Beijing: China communication press, 2007: 95-96.
[14] Duncan J M, Wright S G. Soil strength and slope stability[M]. Hoboken: John Wiley & Sons, 2005: 105-107.
[15] Bell J M. Dimensionless parameters for homogenous earth slopes[J]. Journal of Soil Mechanics and Foundation Engineering Division, ASCE, 1966, 92(5): 51-65.
[16] Dawson E M, Roth W H, Drcscher A. Slope stability analysis by strength reduction[J]. Geotechnique, 1999, 49(6): 835-840.
[17] 徐汉斌, 王军. 反算法中滑坡稳定系数的取值问题[J]. 四川地质学报, 1999, 19(1): 86–89.
XU Hanbin, WANG Jun. Selection of landslide stability factor in inverse calculation[J]. Acta Geologica Sichuan, 1999, 19(1): 86-89.
(编辑 赵俊)
收稿日期：2012-07-28；修回日期：2012-10-10
基金项目：国家自然科学基金资助项目(10972238，51174228)；高等学校博士学科点专项科研基金资助项目(20120162110009)；湖南省研究生科研创新项目(CX2012B069)
通信作者：张科(1986-)，男，浙江舟山人，博士研究生，从事边坡工程理论、试验和数值计算研究；电话：15697485024；E-mail: zhangke_csu@163.com
