刘 亚 坤

(西南交通大学土木工程学院，四川 成都 610000)

0 引言

岩土工程中分析边坡稳定性方法主要可分为传统的极限平衡法和有限元法。传统的极限平衡法以刚体极限平衡理论为基础，分析结果物理意义明确，但需事先假设滑动面的形状、位置，然后不断求解对应的安全系数，直到最小安全系数，此时对应的滑动面为临界滑动面。该方法不仅计算量大，而且难以保证临界滑动面的准确性，尤其是对复杂地质边坡的实用价值大打折扣；而有限元法不仅能够模拟边坡的变形过程及其滑动面的形状，搜索临界滑动面时无需事先假定滑动面形状和位置[1-3]。

Davidon[4]和CHEN Z Y & SHAO C M[5]基于最优控制理论，提出了能最快速搜索目标的变度量法，这为边坡临界滑动面快速搜索提供了新思路，但其计算精度必须建立在准确确定初始滑动面的基础之上。已有研究表明，有限元强度折减法能较为精确的确定初始滑动面。

为此，本文借鉴有限元强度折减法获取初始滑动面、局部优化法的单变量精确计算临界滑动面和最快速获得搜索目标的边度量法的优点，提出了一种能编程简单、搜索效率较高且能搜索任意形状边坡的临界滑动面交替变度量法搜索方法。

1 交替变度量法基本原理

1959年，D’esopo[4]把多变量求最优值问题转化成单变量轮流求最优值问题，从而提出了坐标轮换法——以节点坐标为变量，每次搜索只改变一个变量，其余变量保持不变，即沿坐标方向轮流进行搜索的寻优方法；1959年Davidon[4]提出变度量法。CHEN Zuyu, SHAO Changming[5]于1988年将最优控制理论中变度量法改进，成功的用于边坡滑动面的搜索中。本文借鉴坐标轮换法及变度量法的原理，将多变量问题转化成单变量问题，利用变度量法逐一搜索每个节点的最优值，从而提出了交替变度量法。

为确定临界滑动面，利用有限元强度折减法确定等效塑性应变区域，利用样条函数沿最大塑性应变区域拟合一条光滑曲线(最大塑性应变值连线为波浪形[6])，作为初始滑动面。在初始滑动面上确定若干节点，先令其余节点固定不变，使第一个节点(坡脚节点)沿边界方向进行搜索，获取最小安全系数，如图1a)所示；然后沿着最优化方向移动另一节点，获取最小安全系数，其余节点仍然保持固定，见图1b)；滑动面上所有节点均按照这种方法搜索最小函数值，直至得到一次循环的最小安全系数。

利用交替变度量法确定边坡最小安全系数和临界滑动面的具体步骤如下：

1)通过样条函数沿最大塑性应变值拟合的滑动面作为初始滑动面。

2)i为循环变量，第一次循环记为i=1。

3)在初始滑动面上确定若干节点，节点位置需分布均匀且满足运动学规律，即节点之间距离尽可能相等。节点间的距离过小则无法满足搜索灵活性，距离过大，则可能无法满足运动学规律。节点之间用直线段连接，初始滑动面上节点的坐标为Zi=(x1,y1,x2,y2,…,xm,ym)。

继续对第二个节点进行沿最优化方向搜索，其余节点不变，直至获得最小安全系数，最优化方向记为Sk，按式(1)确定。

(1)

其中，Hk为对称矩阵。本文节点移动的步长λ采用“黄金分割”法沿Sk方向进行一维搜索。第二次搜索后临界滑动面上的节点坐标变为：

5)按步骤4)对除最后一个节点外的剩余节点依次进行搜索，直至第(m-1)个节点搜索完毕，此时临界滑动面的坐标变为：

6)对最后一个节点(第m节点)沿边界移动搜索，直到获得对应的最小安全系数值。

7)重复4)～7)步，当最后两次循环的结果差值满足式(2)时，循环搜索结束，最后一次循环所得的坐标连线即为临界滑动面，所得最小函数值即为边坡安全系数。

2 算例分析

对澳大利亚计算机应用协会(ACADS)提供的均质边坡算例[7]进行计算验证，并与几种常用的极限平衡法计算结果进行了比较，验证了上述方法的正确性及优越性。

该例来自ACADS考核题EX1(a)，土体重度γ=20.0 kN/m3，黏聚力c=3.0 kPa，内摩擦角φ=19.6°，弹性模量E=1.0×104kPa，泊松比v=0.25，具体材料参数、几何参数见I.B.Donald和P.Giam[7]的研究。分别采用属于圆弧法的瑞典法、简化Bishop法以及属于非圆弧法的Morgenstern-Price法与本文方法进行最危险滑动面搜索和安全系数计算。

图2为初始滑动面位置及计算节点的布置。由文献[6]知，临界滑面是由沿深部方向上的等效塑性应变的极大点所组成，但最大塑性应变值的连线成波浪形，故本文采用样条函数沿最大塑性应变区域拟合一条光滑曲线作为初始滑动面。然后等间距选取8个计算节点，通过交替变度量法搜索每个节点最优值(最小安全系数值对应的节点坐标)。

图3为交替变度量法及几种常用传统极限平衡法搜索的临界滑面位置，可以看出，本文方法确定的临界滑动面总体上要比圆弧法确定的临界滑动面更深，但比Morgenstern-Price确定的临界滑动面更浅。本文所得结果为折线型滑动面，计算节点仅为8个，但已然很接近于圆弧状滑面，若增加滑动面计算节点数，则最终搜索结果基本等同于圆弧状滑面。

表1给出了几种传统极限平衡法及本文方法计算得到的安全系数，虽然本文方法计算结果小于算例给出的答案，但只相差1.2%，比瑞典法、简化Bishop法都接近算例答案，几乎与Morgenstern-Price法结果相同。

表1 安全系数计算成果汇总表

表1还总结了不同方法尝试搜索滑动面的次数。可以看出，传统极限平衡法中瑞典法搜索次数较少，为140次，Bishop法搜索次数为400次，Morgenstern-Price法搜索次数明显大于其他传统极限平衡法的搜索次数，达到了667次。从表1搜索时间还可以看出，本文方法仅用时3 min，相较瑞典法、Bishop法和Morgenstern-Price法分别降低了9.33，14.00和25.33倍。结合表1的安全系数可知，瑞典法的精度较差，安全系数为0.942，Bishop法精度适中，安全系数约为0.986，Morgenstern-Price法精度最高，达到了0.989，而采用本文方法，计算程序只需循环搜索15次，最后两次计算结果的差值就可满足误差要求|Fs(Zi+1*-Fs(Zi*)|<1×10-5，此时安全系数达到0.988。可见，与传统极限平衡法相比，本文方法在保证精度的同时能大大提高搜索效率。

传统极限平衡法假设在整个临界滑动面上安全系数相等，这并不符合实际情况，边坡的破坏并不是整体同时破坏而引起的，而是局部破坏的逐步扩展从而引起整体失稳，本文方法恰好能反映这一点。图4给出了临界滑动面上安全系数的变化，每点位置与图3坐标表示一致。由图4可以看出，边坡坡脚处较为稳定，沿着坡脚到坡顶安全系数逐渐降低，接近坡顶时稳定性又逐渐增加，与边坡实际破坏机理基本一致。

3 结语

1)通过有限元强度折减法确定初始滑动面，再利用交替变度量法可以快速、有效的搜索临界滑动面。通过ACADS的边坡算例验证，其在边坡滑动面的搜索中的应用是可行的，且搜索效率明显高于传统极限平衡方法。

2)利用交替变度量法搜索结果虽然为折线型滑动面，但只要增加计算节点，就能很好的拟合为圆弧状滑面。因此利用交替变度量法可以搜索任意形状的滑动面。为了提高计算效率，在保证准确性的情况下，计算节点数量以8个～15个为宜。

3)本文方法以有限元强度折减法作为计算基础，可以真实反映滑动面不同位置安全系数的变化，符合边坡破坏的实际机理。

参考文献：

[1] 郑 宏,李春光,葛修润.求解安全系数的有限元法[J].岩土工程学报,2002,24(5):626-628.

[2] 李春光,朱宇飞,刘 丰,等.基于下限原理有限元的强度折减法[J].岩土力学,2012,33(6):1816-1821.

[3] 杨光华,钟志辉,张玉成,等.用局部强度折减法进行边坡稳定性分析[J].岩土力学,2010,31(S2):53-58.

[4] 马昌凤.最优化方法及其Matlab程序设计[M].北京:科学出版社,2010.

[5] CHEN Zuyu,SHAO Changming.Evaluation of factor of safety in slope stability analysis[J].Canadian Geotechnical Journal,1988,25(4):735-748.

[6] 孙冠华,郑 宏,李春光.基于等效塑性应变的边坡滑动面搜索[J].岩土力学,2008,29(5):1159-1163.

[7] I.B.Donald, P.Giam. The ACADS slope stability programs review. Proc[A].The 6th International Symposium on Landslides[C].1992:1665-1670.