联合运用改进的极限平衡法和动态规划法分析边坡稳定性
吴海真
河海大学水利水电工程学院
江西省水利科学研究院
顾冲时
河海大学水利水电工程学院
为了解决边坡关键滑动面全局优化处理问题，将动态规划理论引入基于有限元的改进极限平衡法。与工程上被广泛认可的一般极限平衡法比较表明，本文提出的方法可以获得更好的近似解甚至是整体最优解，关键滑动面更符合工程实际，有效地克服了传统边坡稳定分析方法的不足；与目前国内外同类方法相比具有明显的优越性，丰富了边坡稳定分析和确定关键滑动面的方法和理论。
1 研究背景
目前，工程界对于边坡稳定分析评价常采用刚体极限平衡法，而土坡稳定分析计算中又以垂直条分法应用最为广泛，岩坡稳定分析计算中以Sarma法应用最为普遍。这些经典极限分析方法概念清晰，使用时间长，积累了丰富经验，但在理论上存在诸多不足，主要有：①没有考虑岩土体内部的应力应变关系，边坡失稳与变形无确定性联系，滑动面底部应力可能失真；②不能考虑岩土体与支挡结构的共同作用及其变形协调；③假定条件较多，如土条条分假设，条间力及滑动面假定等。
随着计算机技术的迅速发展，特别是有限单元法提出后，有力推动了边坡稳定分析方法的发展。有限元法已广泛应用于边坡稳定性分析，与极限平衡法相比，有限元法有如下优点：①可同时进行稳定分析与变形分析；②可考虑土的非线性应力应变关系；③可保证整个边坡同时满足力和力矩平衡条件；④可考虑先期应力状态和模拟开挖或填筑等施工过程。有限元法分析边坡稳定性一般可分成两类：有限元强度折减法和有限元与优化搜索法相结合的方法。早在1975年Zienkiewicz等［1］就通过等比折减土的黏聚力和内摩擦角正切值的方法分析边坡的稳定性，把使边坡达到临界稳定状态的材料强度折减系数值作为边坡的最小安全系数，并把该方法命名为“强度折减技术”。郑颖人等［2］采用有限元强度折减法进行土坡和岩坡的稳定分析。张鲁渝等［3］对有限元强度折减系数法计算土坡稳定安全系数的精度进行了研究。但是，有限元强度折减法在计算分析前对土的强度参数进行了折减，取得的应力与变形结果无法反映边坡实际的受力状态，因此该方法只能限于边坡最小安全系数的评价，而在边坡失稳判据等方面仍值得深入研究。另一方面，确定临界滑动面位置及形状是边坡稳定分析的主要任务。一般的边坡稳定分析工具要求分析者在计算前输入滑动面的位置，这对工程技术人员的理论水平和工程经验提出了较高的要求。因此，从20世纪70年代后期开始，很多学者致力于临界滑动面自动搜索和优化的技术研究，提出了多种确定边坡圆弧或非圆弧临界滑动面的优化搜索方法。Zhu［4］提出了边坡临界滑动场的概念，陈昌富等［5］把自适应蚁群算法引入边坡稳定分析领域。但由于蚁群算法与临界滑动场本身无法直接应用于连续性问题，因此在一定意义上这两种方法难以找到边坡的全局最小安全系数。Zou等［6］提出了采用动态规划法或改进的动态规划法确定非圆弧临界滑动面及其对应最小安全系数的方法，但该动态规划法仍局限于常规的条分法计算。Goh［7］结合滑动楔体理论，运用遗传算法对多折线滑动面进行了搜索。Cheng［8］采用模拟退火法成功分析了多个较复杂边坡，但其对于如何利用模拟退火法搜索临界滑动面的方法介绍比较简单。
有限元法与临界滑动面搜索优化相结合的思想源于极限平衡法的临界滑动面搜索，以充分利用有限元法的优势对边坡的稳定与变形进行分析评价。本文采用动态规划理论与有限元相结合的方法确定边坡的临界滑动面，其中利用有限元法分析某状态下边坡的应力和孔隙水压力分布情况，按有效应力法计算给定滑动面的安全系数，同时引入动态规划理论，成功地解决了对边坡关键滑动面的优化处理问题，有效地克服了传统边坡稳定分析方法的不足。
本文发表于2007年。
2 改进的极限平衡法基本理论
传统的极限平衡法需要对边坡条块之间的作用力进行一定的假定，使之成为静定问题，由此计算的滑动面上应力并不精确。众所周知，有限元法能够同时满足力和力矩平衡，从而可以进行准确的应力分析。因此可利用有限元法分析的应力结果计算边坡的安全系数，也即基于有限元法得出边坡内部的应力场，然后基于传统极限平衡法安全系数的定义，采用优化方法搜索关键滑动面的位置并得到相应的最小安全系数。Farias等［9］提出了有限元法计算边坡安全系数的详细方法，本文基于Terzaghi有效应力原理对该方法进行了改进。
2.1 安全系数的定义
据毕肖普20世纪50年代重新定义的边坡稳定安全系数概念，每个土条的安全系数Fsi等于该土条底部的抗剪强度与土体的剪应力之比［10］：
式（1）即为边坡沿滑动面的局部安全系数。对于整个滑动面，边坡稳定安全系数定义为滑动面底部的抗滑力之和与滑动力之和的比值：
式（2）即为边坡沿滑动面的整体安全系数。对于圆弧滑动面，式（2）即成为沿滑动面的抗滑力矩与滑动力矩的比值。
2.2 边坡内任意点应力状态
首先采用弹塑性有限元方法计算边坡内的应力状态，得到单元高斯点应力（本文采用四节点等参单元）。为了得到土条底部中点的应力值，把单元高斯点应力投影成单元节点应力：
式中：{σ}e为单元节点应力；［N］为位移插值函数矩阵；{σ}g为单元高斯点应力。
对所有节点进行绕节点平均，从而得到边坡内部各单元的节点应力值，进而可根据计算点的局部坐标关系，计算边坡内任意点的应力状态和孔隙水应力。
2.3 边坡滑动面应力状态
将滑动面所包含的土体进行条分，设任意土条底部长为ΔLi，则滑动面上一点的抗滑力τf和滑动力τn可由该点有限元分析的应力值σx、σy和τxy和孔隙水应力p计算，公式如下：
式中：σx、σy和τxy为土中一点正应力和剪应力，以拉为正；φ′和c′分别为土体有效内摩擦角和有效黏聚力；σ′n和σn分别为沿滑动面上的有效正应力和总应力，以压为正；p为孔隙水应力；α为土条底部中点之切线与x轴正向夹角，逆时针方向转动为正。
由此可根据式（1）和式（2）计算边坡的局部安全系数和整体安全系数。
3 动态规划法确定边坡关键滑动面方法
3.1 动态规划法
动态规划法是20世纪50年代由Bellman和Dantzig发展起来的多阶段决策优化理论，其中的费用最小值问题表述为［11］：
其中x（0）有固定值C，系统方程为：
应用动态规划方法研究上述问题就是希望把实际问题嵌入到一类似的问题中去。一般情况下，直接求解上述问题是困难的，但可以通过连接这一类问题中各组成部分的关系式即可求极值的嵌入递推方程进行求解。用函数I（x，k）规定为最小费用，则嵌入递推方程为：
基于嵌入递推方程，根据Bellman的阐述［11］，一个最优策略有这样的特性：不论初始状态和初始决策如何，相对于第一个决策所形成的状态来说，剩下的决策必定构成一个最优策略，此即动态规划法的基本原理。
3.2 数学模型的建立
本文基于条分法和有限单元法相结合的基本原理，引进动态规划理论确定边坡的关键滑动面，其数学建模过程如下：
引进辅助函数G：
式中：R为抗滑力；S为滑动力。
可以设想，要使式（2）的安全系数Fs最小，等价于使式（11）中的G达到最小：
为引进动态规划，用适当多的垂直线（共n+1条，也即单元边界线）条分边坡，如图1所示。每一垂直线称为阶段（stage），初始阶段和最终阶段表示分析区域，滑动面与条分线的交点（单元节点）为状态点（state point），可以设想边坡临界滑动面必然是由所有各阶段的某些状态点连成的折线。取出任一两相邻阶段（即两条分线）i和i+1考虑，如图2所示，假定连接点j和k之连线jk是可能的临界滑动面，计算滑动面所穿过单元的Ri和Si，则
式中：
；
。
图1 动态规划法的网格划分
图2 滑动面穿越相邻阶段示意
引入一个“优化函数”Hi（j），表示从初始阶段到i阶段j状态即（i，j）之间G的最小值。根据Bellman的动态规划原理［11］，其嵌入递推方程为：
式中：Hi+1（k）为在i+1阶段状态点k获得的优化函数；Hi（j）为在i阶段状态点j获得的优化函数；
为滑动面穿越i阶段的状态点和i+1阶段的状态点k所获得的函数值。
在初始阶段：H1（j）=0，j=1，2，…，Np1，Npi为第i阶段的状态点数。
到最终阶段：
以上即为应用动态规划方法确定边坡关键滑动面的数学模型。
3.3 相关问题的处理
（1）关键滑动面与坡面交点的确定。据上述由动态规划法原理确定边坡关键滑动面的方法进行搜索，可以得到边坡的临界滑动面。同时因应力值已经由有限元方法计算得到，故安全系数Fs的计算是线性问题，不需要迭代。但在实际计算时，必须扩大计算范围，即在关键滑动面与坡面的可能交点部位的坡体以外也要划分足够多的状态点，此时只要令抗滑力R和滑动力S为零即可，而且初始及终了阶段的状态点都必须在坡外，这主要是为了确保得到完整的滑动面，见图3。
图3 确定关键滑动面与坡面的交点示意
（2）搜索关键滑动面的起点设置。关键滑动面搜索起点的确定主要决定于初始阶段状态，如果初始阶段状态点总数为Np1，则逐个以初始阶段状态点为起点进行分析计算，得到相应的安全系数和最小值，其相应初始阶段之状态点即为关键滑动面搜索起点。
（3）单元网格密度。单元网格密度决定垂直条分线之间的宽度和状态点的间距，从而直接影响计算的精度和搜索时间，这主要取决于边坡的地质条件、规模和工程本身的精度要求。通过大量算例证实，对于4节点矩形单元，当单元密度达到每10m2不少于4个节点时，计算精度较为理想，如果再增加节点，计算精度应还能提高，但此时耗费的机时也将成倍增长。
4 算例分析
根据以上理论和方法，编制了边坡稳定改进极限平衡法的计算程序，采用动态规划法确定边坡的关键滑动面位置和形状，可以考虑地下水和外荷载的影响，尤其适用于复杂地质条件下的边坡稳定性分析研究，程序计算速度快，收敛性好。
4.1 算例1
图4所示为算例1的边坡剖面图。首先采用Morgenstem-Price法进行计算，对如图4所示的物理力学参数和剖面，取土条数30，穷举法搜索关键滑动面，相应安全系数为1.170；基于改进的极限平衡法基本原理，取弹20MPa，泊松比0.33，得到边坡的整体安全系数为1.130；若根据上述动态规划法的原理和步骤在改进的极限平衡法基础上再次进行关键滑动面的搜索计算，得到的边坡安全系数则为1.05。3种计算方法得到的关键滑动面位置如图5所示。
图4 算例1的计算剖面
图5 算例1的计算成果
4.2 算例2
计算一坡内具自由水面、坡外有水的情况，图6所示为算例2的边坡剖面及土层分布图，各土层的物理力学参数取值见表1。
图6 算例2的计算剖面
表1 算例2的物理力学参数取值
采用Morgenstern-Price法进行计算，边坡相应的安全系数为1.14；采用基于改进的极限平衡法，动态规划法优化滑动面位置及形状，得到边坡的整体安全系数为1.06。两种计算方法得到的相应关键滑动面位置如图7所示。可见采用改进的极限平衡动态规划法可获得更好的优化滑裂面及相应的近似解甚至整体最优解，滑裂面位置和形状也更符合实际。
需要说明的是，上述两个算例中Morgenstem-Price法与本文方法计算结果的差异，并非Morgenstern-Price法本身所造成，而主要是由于优化处理后滑裂面形状和位置的不同所致。
图7 算例2的计算成果
5 结语
边坡稳定性分析一直是岩土工程中的重要研究课题，而极限平衡法又是边坡稳定分析理论中的重要内容，但在理论上存在诸多不足，其中边坡内部应力尤其是滑动面底部应力“失真”是突出问题之一；另一方面，如何合理的确定关键滑动面是提高边坡稳定分析准确性的关键。与传统的条分法相比较，本文提出的基于改进极限平衡的动态规划法的特点在于：
（1）把有限元与极限平衡结合起来，从而使该方法本身具有有限元法的一切优点，并在原有方法的基础上引入了有效应力原理，能够考虑坡内具自由水面、坡外有水的情况，可较为真实地反映边坡体内部的应力和变形。
（2）基于改进的极限平衡法，引入动态规划理论与之相结合，成功地解决了对边坡关键滑动面的优化处理问题。算例分析表明，与工程上被广泛认可的Morgenstern-Price法比较，本文提出的方法可以获得更好的优化滑裂面及相应的近似解甚至整体最优解。
（3）本文提出的方法和理论为复杂边坡稳定分析提供了新的思路，丰富了边坡稳定分析和确定关键滑动面的方法和理论。
参考文献
［2］郑颖人，赵尚毅.有限元强度折减法在土坡与岩坡中的应用［J］.岩石力学与工程学报，2004，23（19）：3381-3388.
［3］张鲁渝，郑颖人，赵尚毅，等.有限元强度折减系数计算土坡稳定安全系数的精度研究［J］.水利学报，2003，（1）：21-27.
［5］陈昌富，龚晓南，王贻荪.自适应蚁群算法及其在边坡工程中的应用［J］.浙江大学学报（工学版），2003，37（5）：566-569.
［10］钱家欢，殷宗泽.土工原理与计算［M］.2版.北京：中国水利水电出版社，1995.
［11］罗伯特·E拉森.动态规划法原理［M］.陈伟基，等，译.北京：清华大学出版社，1984.
本文发表于2007年。