李同录，王刘华，张常亮，李 萍

（长安大学地质工程与测绘学院，陕西西安710054）

0 引言

1915—1916年，Peterson在对瑞典Goetebborg采石场边坡破坏形态研究的基础上，首先提出了将均质土坡滑动面近似看作是一个圆弧的边坡稳定性分析方法，此后瑞典岩土工程学会对大量边坡特别是铁路工程边坡的调查也进一步证实了该观点［1］。目前已经有多种基于圆弧滑动面的边坡稳定性分析方法被提出，并在工程中广泛应用［2］。

然而，对于一个特定边坡，计算模型是在给定滑动面的情况下建立的。稳定性评价还需要从无数个给定的圆弧滑动面中找出稳定系数最小的一个作为临界滑动面，这是一项繁琐而复杂的工作。早期由于运算能力低，边坡稳定性分析主要采用图解法。Fellenius首先提出了圆弧条分法，并利用图解法绘制成图表，以便实际应用［3］；Taylor提出的摩擦圆法也是基于图解法制成图表，这些图表在中国的土力学教材中都有介绍［4］。

随着计算机的应用，处理大量的滑弧滑动面已经不再是难题。当圆心和半径大概的范围给定，就可以对圆心区域网格化，对半径按步长增加，计算每个圆心网格点和给定半径的对应圆弧滑动面的稳定系数。根据所有圆心网格点的稳定系数便可绘制稳定系数等值线，由此可得出稳定系数的最小值及与之相对应的滑弧参数。如果斜坡是均质土坡，那么只有一个极值，但如果由不止一种土层组成，且各土层性质存在明显区别，那么可能存在几个极值，其中最小的一个就是所期望的值，目前普遍使用的商业软件Geoslope就是采用这种算法［5－6］。一个圆有3个自由度，一般是圆心的x、y坐标和半径R。在某些情况下，如粘聚力c趋于0时，半径R可能趋于无穷大，x和y坐标远离斜坡，因此不容易得到一个包括临界圆中心的限制区域（图1）。在这种情况下，以上方法是不严密的。针对该方法的不足，李同录等提出了一种更为严格的二维搜索方法［7］，笔者将该方法进一步推广到三维滑动面进行搜索。

1 临界滑动面搜索的优化算法

图1 滑弧中心无穷远与半径无穷大的特殊情况Fig.1 Special Case That the Center and Radius of Slip Circle Is Far from the Slope

三维临界滑动面的搜索以二维为基础，因此先介绍二维方法。建立如图2所示的直角坐标系，假定x轴的正方向指向坡内，z轴垂直，设点O（xO，zO）是给定滑弧的圆心（其中xO和zO分别代表点O的横坐标和纵坐标，其他点类推），点A（xA，zA）和B（xB，zB）分别是剪入点和剪出点，通过A点做圆弧的切线并与x轴相交于点T（xT，0），首先，将圆心坐标xO，zO和半径RO等3个参数用A、B和T等3个点的坐标表示。

图2 由A、B两点和切线AT确定的圆弧滑动面Fig.2 Illustration of Slip Surface Determined by Two Points A，Band Tangent Line AT

则圆弧滑面高程z的表达式为

滑坡坡面函数zs为

假定xA、xB给定，则zA、zB可以通过式（5）求得，且zT＝0。因此，滑弧实际是由xA、xB与xT3个参数控制的，且它们都可以控制在合理的有限域内。从数学和几何方面而言，xA和xB的定义域分别为xA∈（xE，＋∞）、xB∈（－∞，xF），然而此处不难给出一个合适的xA的右边界值和xB的左边界值。假定这两个边界值分别为xD、xG，则xA和xB的定义域分别为xA∈（xE，xD）、xB∈（xG，xF）。对于xT，由于滑弧不可能越过最右侧的AI线，所以它的右边界值应该为xI＝xA，而它的左边界值为xH。由于滑弧上的任何一点都不可能出现在边坡坡面的上方，基于该限制条件，很容易确定xH的横坐标，并得到xT的定义域xT∈（xH，xA）。对于任何圆弧条分法，稳定系数F都可以表示为通用函数式中：γ为土的重度；c为粘聚力；φ为内摩擦角。 这3个参数反映土的属性，可通过试验结果或者经验确定；坐标xA、xB与xT为变量。接下来找到xA、xB与xT的一个确定值，使稳定系数极小，该问题用黄金分割法来求解。

如果边坡土体是均匀的，F是凹函数，那么在整个循环过程中对应每一个xA、xB与xT分别应用黄金分割法，在循环运行结束后，就可以得到使F达到极小值所对应的xA、xB与xT的值，黄金分割法的整个过程如图3。最后，由式（1）～（3）计算临界滑弧的圆心和半径。

图3 用黄金分割法搜索土坡临界滑动面的流程Fig.3 Process for Searching the Critical Slip Surface by Golden Block Optimization Method

如果边坡由不同土层组成，那么F可能会存在局部极值。在这种情况下，黄金分割法可能会找到某一个极值点，但是并不一定是F的最小值。因此，首先按一定步长改变xA、xB与xT，找出所有局部极值域，然后在每一个极值域中运用黄金分割法求极小值，所有极小值中的最小值就是最终的结果。

基于上述讨论，对于给定的土质边坡，搜索其临界滑动面的具体步骤。

（1）确定xA、xB与xT的合理区间，即xA∈（xE，xD）、xB∈（xG，xF）、xT∈（xH，xA）。

（2）按xA、xB、xT的顺序，采用黄金分割法由外向内划分黄金分割点，当到最里圈时就可以确定xA、xB与xT的一个设定。

（3）将xA、xB与xT代入式（1）～（4）计算xO、zO以及滑弧的半径RO，并检验滑弧方程式z。

（4）将式（1）～（3）及相关土性参数代入式（6），计算F。

（5）重复以上（2）到（4），每一次循环都可以得到3个F值，通过比较这些值缩减相应变量的区间。

（6）最后，当xA、xB与xT的所有区间长度减小到给定误差以内时循环结束。此时得到F的最小值及相应的xA、xB与xT值。

（7）将xA、xB与xT代入式（1）～（3）计算xO、zO和RO，并画出相应的圆。

该算法与具体方法无关，对所有基于圆弧滑动面的方法都适用（如简化Bishop法）。

2 三维临界滑动面搜索

二维极限平衡法适合侧向延伸远且坡面形态简单的边坡，该类边坡属于平面应变问题，取其中一个单位宽度的窄条来分析。二维方法没有考虑边坡的侧向约束力，代表了最保守的情况。如果边坡地形条件复杂，如发育在沟谷或凸出的山梁上，边界效应十分突出，采用二维方法可能低估了其稳定性，此时应该采用三维方法。

对于均质土坡，三维临界滑动面一般假定为旋转椭球体。这一方面是大量破坏面观察近似反映出该种形态特征，另一方面该假定保持了与二维圆弧滑动面的假定一致性，因为旋转椭球的任一纵断面都是圆（图4），旋转椭球体的方程可用下式表示

令式（7）的y＝0，则简化为中轴面的圆弧方程。式（7）比圆弧方程多了一个椭球体的横向半径Rb，由此椭球滑动面的方程为

图4 三维边坡破坏面形态Fig.4 Fracture Sufrace Pattern of Three－dimensional Earth Slope

对于任一三维稳定性计算方法，其稳定系数可表示为

由式（9）可见，三维滑动面的搜索比二维多了一个自由度Rb，即旋转椭球的横向半径。式（9）适合具有旋转椭球滑动面的任何三维极限平衡法。

以一种具体方法为例，说明三维临界滑动面的搜索。目前的三维极限平衡法，基本上是基于二维极限平衡法的扩展。张常亮等建立了三维通用极限平衡法［8］，给定特定的假定条件，可以简化为具体方法。为了和二维简化Bishop法比较，将通用法简化为三维Bishop法［9］其计算公式为

式中：dW为微条柱的体力；dAz为微条柱底滑面的面积；dNz为微条柱底部法向力；u为孔隙水压力；FM为基于整体力矩平衡的稳定系数；αx为xOz平面上滑动面与x轴的夹角；αy为yOz平面上滑动面与y轴的夹角；αx、αy分别用滑动面方程求导来确定，即

Ry为椭球在垂直y轴、对应y取值的截面圆的半径，由式（7）求得，即

当y＝0时

将以上参数代入式（10）就简化为二维简化Bishop法的公式，即

式中：mα＝cosαx＋tanφsinαx／FM。

式（10）的分子和分母项都可以用数值积分，如Gauss积分求解，xL、xH、yL、yH分别代表x和y方向的积分边界。数值积分时，可不考虑潜在滑动区的实际边界，积分边界取x和y方向边界（即椭球和地面的交线）的最大、最小值，相当于在一个矩形域求解。计算时对每一积分点做判断，若地面高程小于椭球面，则不参与计算，其结果等效于按实际边界积分。

二维圆弧滑动面搜索时，FM随着xA、xB、xT的变化总能找到极小值。通过文献［10］的一个算例来分析FM随Rb的变化规律。该坡体为一均质土坡，土性和几何参数如图5，并假定横向无限延伸。首先利用二维简化Bishop法找出二维临界滑动面，其圆心为（36.6，27.4），半径RO为24.4m，稳定系数为2.08。

图5 用于三维计算的一个算例Fig.5 Worked Case for Three－dimentional Calculation

分别取Rb为RO的0.5、1.0、1.5、2.0、2.5、3.0、4.0、6.0倍，计算三维稳定系数。给定Rb后，式（9）的自由变量与式（6）相同，除计算公式采用式（10）外，与二维算法相同（图6）。

图6 算例和案例中椭球长短轴之比与稳定系数的关系Fig.6 Correlation Between the Ratios of Major and Minor Axises and the Factors of Safety for Worked and Real Cases

由图6可见，边坡的三维稳定系数随着滑动面横向宽度的增大而减小，当长轴与短轴之比Rb／RO小于3时，三维效应十分明显，三维稳定系数随着Rb／RO的增大而急剧减小；Rb／RO大于3时，稳定系数的变化趋于平稳；当Rb→∞时，三维FM趋近于二维FM。李同录等利用三维Sarma法也得出了类似的结果［11］。

以上结果表明，三维稳定系数随着滑动面宽度增大而单调减小，不存在极小值，但存在极限值，当宽度无限大时，极限值逼近二维稳定系数，二维稳定系数是三维的下限。但三维稳定系数随滑动面宽度变化具有明显的转折点，若用椭球滑动面的横纵半径之比大于3时，边坡稳定性分析可简化为二维问题；小于3时，应当按三维分析。前者如道路工程中开挖的长边坡、河流侵蚀形成的直线岸坡等；后者如发育在沟槽或梁峁处的自然边坡或者基坑的拐角处等。

横向无限延伸的边坡实际属于二维问题，所以三维结果逼近二维。具有三维属性的边坡，边坡破坏范围可结合边坡实际边界条件加以界定，再利用椭球滑动面拟合，以便计算稳定系数。给定约束条件后，也可能存在极值。

3 工程实例

现以西安咸阳国际机场污水排放口黄土高边坡为例说明空间滑动面搜索方法的具体应用。该边坡位于泾阳县高庄镇傅家村黄土塬的北缘，咸阳机场的污水经过净化处理后，用暗管引到此处塬顶，然后从塬顶引到坡脚，再经明渠流入泾河。该边坡高85m，坡度40°～50°，最陡处可达65°。斜坡北侧为黄土塬，南侧坡脚为泾河一级阶地；东侧为人工开挖的弧形高边坡，坡脚建有泵站，用于汲取泾河水，坡顶为废弃的砖瓦厂；西侧为天然冲沟，宽30余米，深22～50m，沟头周围为墓地（图7）。

图7 西安咸阳国际机场排污口边坡的卫星影像Fig.7 Satellite Image of the Slope near the Outlet of the Sewage of Xi’an Xianyang International Airport

咸阳国际机场污水渠于2003年开始运营，此前未做详细勘查。2006年3月春灌时，有流水加泥从排放口西约60m处坡脚的裂缝中流出，淹没塬下部分麦田，在斜坡上也形成冲沟。随后到塬顶查看，发现大量灌溉水流入一张开的裂缝，裂缝总长近100m，局部宽度达120cm，平行边坡近东西向展布，距坡顶边缘43m。调查当地村民得知，该裂缝在2001年3月春灌时就已出现，呈锯齿状延伸，出露长度约100m，宽10～60cm，停止浇地后裂缝未继续发展，并因农耕被掩埋。其后，该裂缝在每年灌溉时，仍时断时续出现。为了防止灌溉水的进一步灌入，2006年春灌后，当地村民沿裂缝开挖沟槽，用隔水土工布和灰土夯填进行处理。但2008年春灌后，发现该裂缝再次张开，说明该斜坡处于不断蠕变中。该斜坡滑动可能直接造成污水排放设施的损毁，因此引起了机场管理部门的高度重视，要求对其做稳定性评价。

构成该边坡的地层主要为第四系风积黄土夹古土壤层。黄土呈棕黄色，结构均匀，具大孔隙，垂直节理发育；古土壤呈褐红色，具有粒状结构。自坡顶到坡底可见L1—S9层黄土－古土壤序列，下伏河流冲积相的致密黏土。坡面一般覆盖10cm至数米的坡积黄土状土，边坡地层剖面如图8。

图8 西安咸阳国际机场排污口边坡地质断面Fig.8 Geological Profile of the Slope near the Outlet of the Sewage of Xi’an Xianyang International Airport

该边坡的变形主要由灌溉水下渗引起。黄土下部的致密黏土是相对隔水带，灌溉水沿坡顶裂缝集中下渗，到该层顶面被阻滞，使地下水位大幅度上升，图8斜坡上钻孔揭露的地下水位高程位于坡脚高程以上，灌溉时水位可能会更高。因此该边坡剪出口应位于潜水面附近的软弱带。

根据钻孔试样的物理力学性质测试结果，土的重度和有效强度参数取值为：γ＝18.7kN／m3；γsat＝19.8kN／m3；c′＝45kPa；φ′＝30°。

为了分析灌溉水对边坡稳定性的影响，分地下水位低和高两种情况计算。前者是指地下水位在致密黏土中，即剪出口以下；后者指地下水位在坡脚以上的黄土中，按实测水位计算，地下水的影响范围为灌溉区域，两侧的砖瓦厂和墓地没有地表水补给，可认为地下水位深。

首先根据地形和地表裂缝展布确定中轴面，再给定不同的椭球长轴和短轴的比例计算稳定系数（图5）。根据滑坡两侧地形条件的限制和后缘裂缝的展布范围，当RO为138.8m、Rb为267.5m时，滑动面椭球和地面的交线与后缘裂缝和侧缘边界最为吻合，因此应取Rb／RO为2.5对应的稳定系数为实际稳定系数值。由图6可见，不考虑灌溉的情况下，边坡稳定系数为1.12，处于基本稳定状态；考虑灌溉导致地下水位升高的情况下，边坡的稳定系数为1.04，由于钻孔所测水位不在灌溉期，灌溉期水位可能更高，边坡稳定系数在灌溉期更低，此时边坡接近临界状态。由此可见，该边坡的变形主要由农田灌溉引起。

为了分析潜在破裂面宽度对稳定性的影响，图6给出了不同的Rb／RO对应的稳定系数，可以看出，不考虑灌溉时，稳定系数的变化和图5的算例规律相同，随Rb／RO的增大，稳定系数减小，Rb／RO大于2.5时，趋于一个稳定值1.12。考虑灌溉时，总体规律也一样，但在Rb／RO为2.5时，达到极小值1.04，随后缓慢增大，此时的滑动面与实际边界条件吻合最好，这说明在考虑实际地下水和地形等边界条件的情况下，三维稳定性计算也存在极小值。与不考虑地下水的情况比较表明，该极小值和地下水位局部上升有关，灌溉引起的高地下水位分布在张裂缝带所在的有限范围，潜在滑动面底部与地下水位面达到最大吻合时，稳定系数达到极小，潜在滑动面再扩大，两侧水位降低，滑动面在水位线以上，则整体稳定性提高，稳定系数增大。由此可见，三维边坡稳定性评价时，现场确定边界条件尤为重要。

4 结语

（1）均质土坡三维潜在滑动面可按旋转椭球体来近似拟合，这和实际观测基本吻合，也和二维圆弧滑动面的假定一致，由此可将二维滑动面的搜索方法推广到三维滑动面。

（2）对三维旋转椭球滑动面边坡稳定系数的计算结果表明，侧向无限长边坡，三维稳定系数没有极小值，但有极限值。因此，当边坡侧向延伸远，使旋转椭球长短轴之比大于3时，边坡稳定性分析可简化为二维问题；当侧向受地形条件或岩土条件限制，边坡局限于一定范围时，应按三维分析。

（3）按三维分析时，应结合边坡的实际边界条件，包括地形条件、岩土性质和地下水动态等，以便使边坡潜在破裂面与边界条件吻合，此时可能存在稳定系数的极小值。

［1］ Taylor D W.Stability of Earth Dams［J］.Journal of Boston Society of Civil Engineering，1937，24（3）：197－246.

［2］ Bishop A W.The Use of Slip Circle in the Stability Analysis of Slopes［J］.Geotechnique，1954，5（1）：7－17.

［3］ Fellenius W.Calculation of the Stability of Earth Dams［C］∥United States Government Printing Office.Transaction of the 2ndCongress for Large Dams.Washington DC：United States Government Printing Office，1936：445－462.

［4］ 陈希哲.土力学地基基础［M］.第2版.北京：清华大学出版社，1989.

［5］ 郑 涛，张玉灯，毛新生.基于Geoslope软件的土质边坡稳定性分析［J］.水利与建筑工程学报，2008，6（1）：6－8.

［6］ 王 胜，李云华，苏 谦.基于Geoslope的高陡路基稳定性分析［J］.路基工程，2008，20（5）：138－140.

［7］ 李同录，邓宏科，李 萍，等.搜索简单土坡潜在滑动面的一种新方法［J］.长安大学学报：地球科学版，2003，25（3）：56－59.

［8］ 张常亮，李同录，李 萍.三维极限平衡法通用形式的建立及应用［J］.地球科学与环境学报，2010，32（1）：98－105.

［9］ Hungr O.An Extension of Bishop's Simplified Method of Slope Stability Analysis to Three Dimensions［J］.Geotechnique，1987，37（1）：113－117.

［10］ Zhang X.Three－dimensional Stability Analysis of Concave Slopes in Plane View［J］.Journal of Geotechnical Engineering，1988，114（6）：658－671.

［11］ 李同录，王艳霞，邓宏科.一种改进的三维边坡稳定性分析方法［J］.岩土工程学报，2003，25（5）：611－614.