基于几种最优化方法的边坡稳定分析及在MATLAB中的实现

基于几种最优化方法的边坡稳定分析及在MATLAB中的实现姜健;余湘娟;毛尚礼
【摘要】以国内应用广泛的毕肖普条分法结合MATLAB优化工具箱中具有优良特性的拟牛顿法和单纯形法编程,同时实现用这几种最优化方法求解最小安全系数,并就有关方面予以比较并查看它们的求解差异,以期通过多种方法的比较尽可能地找到真正的危险滑裂面.考虑到初值的选取对优化计算的重要影响,适当地采用了一些合理的方法.
【期刊名称】《水运工程》
【年(卷),期】2011(000)004
【总页数】5页(P9-13)
【关键词】最优化方法;边坡稳定性;MATLAB;毕肖普法
【作者】姜健;余湘娟;毛尚礼
【作者单位】河海大学土木与交通学院,江苏南京210098;河海大学岩土工程研究所,江苏南京210098;河海大学岩土工程研究所,江苏南京210098
【正文语种】中文
【中图分类】TU476
采用条分法分析边坡稳定性已行之已久，目前仍然是该课题主要的分析手段。

利用条分法的首要任务就是确定最危险滑裂面，之后可以用任何一种条分法算得安全系数。

传统的方法是假定许多滑动面，用某种条分法分别计算各滑动面的安全系数，最终通过比较得到最小安全系数，从而确定其对应的滑裂面为最危险滑裂面。

这种
方法得到的所谓的“最小安全系数”往往不是真正的最小安全系数。

于是最优化方法被广泛应用于求解最小安全系数。

如果滑裂面的曲线为y（x），那么，求解最小安全系数这个问题具体化为寻找下列泛函的极值：F=F[y(x)][1]。

岩土工程中的边坡的几何形状各异，以及岩土材料的非均质性等复杂条件决定了纯解析的变分原理很难进行极值计算。

用最优化方法进行数值方法求解，是一个现实可行的途径。

1.1 概述
这里介绍拟牛顿法，因为它是MATLAB软件关于最优化求解的默认方法，该最优化方法的突出优点是收敛速度快，是一种集中了许多种算法优点的一种方法。

经理论证明和实践检验，拟牛顿法已经成为一种公认的比较有效的方法。

它不仅是求解无约束问题非常有效的算法，而且已经被推广用来求解约束极值问题。

由于它避免了牛顿法中计算二阶导数矩阵及其求逆过程，又比梯度法的收敛速度快，而且对高维问题具有显著的优越性，因而拟牛顿法获得了很高的声誉。

牛顿法最突出的优点就是收敛速度快，但是其缺点就是要用到导数，对于复杂的函数这是一个严重的缺陷，即难求其Hesse矩阵及其逆矩阵。

拟牛顿法恰好解决了这个问题，它的基本思想就是利用不包含二阶导数的矩阵Hk取代了牛顿法中的Hesse矩阵的逆矩阵[Δ2f(x(k))]-1，再沿d(k)=-Hkgk作一维搜索。

由于构造近似矩阵Hk的方法不同，因而有不同的拟牛顿法，其中最主要的两种为DFP法和BFGS法，其中后者具有更好的数值稳定性和实用性。

1.2 拟牛顿法简介
在最优化方法的牛顿法中，基本迭代公式是：
式中：x为所求函数的自变量向量；λk为最优步长；d(k)为搜索方向。

代替搜索方向表达式d(k)中的Hesse矩阵的逆矩阵；对于BFGS法，用
代替搜索方向表达式d(k)中的Hesse矩阵的逆矩阵。

具体请参考文献[2]。

另外用
直接搜索法中的单纯形法对拟牛顿法进行校核。

1.3 关于初值
令人遗憾的是，几乎所有的优化方法都存在陷入局部极值的可能，无法获得全局极值，初始估计值对优化问题的影响很大。

经验表明，在下面两种情况下，数值计算的收敛问题变得极为严峻：1）当问题包括较多的自由度时，可以想到，此时目标函数的自变量空间存在许多局部极值点，使得搜索局部极值的任务变得很困难。

2）当边坡断面较复杂时，此时也会出现许多安全系数的局部极值[1]。

为此，笔者认为可以适当采用以下方法：采用不同的初始值来进行优化。

一种方法是采用随机选取的方法，随机产生一组自变量然后将他们作为初值进行迭代，产生的随机数越多，得到合理结果的可能性越大，当然明显不合理的应该去掉。

如果随机产生的数组数量够大，那么可以作为一种独立的计算最小安全系数的方法，称为随机搜索法[1]，这完全得益于计算机计算水平的日益提高。

另一种方法是对自变
量采用黄金分割点选取[3]。

第3种方法可以使用相对简单的瑞典法进行初步计算，并把结果作为毕肖普法的初值。

其次，尽可能使用不同的优化方法检验所得到的结果。

2.1 条分法概述
条分法是以极限平衡理论为基础，假定土体是理想的塑性材料，把土体划分为若干竖向土条，并把土条看做刚体。

极限平衡法是建立在大家熟悉的摩尔-库伦强度准
则基础上的。

极限平衡法的基本特点是只考虑静力平衡条件和土的摩尔-库伦破坏
准则。

在大多数条件下问题是不静定的，其处理对策就是引入一些简化假定，使得问题变成静定的，这种处理使得问题的严密性受到了破坏，但是对计算结果的精度影响并不大。

毕肖普法是条分法的一种，它建立了竖向力的平衡方程和整体力矩平衡方程，另外假定土条竖向条间力为零，并回避了横向条间力，所以它并非严格的条分法。

但是，
其计算方程相对简单，收敛性好而且诸多工程实例证实其精度也能满足工程需要，所以目前在国内较多使用的就是毕肖普法。

毕肖普法对安全系数F采用如下定义:土的实际抗剪强度与保持平衡（指总体平衡）所需要的强度之比，即
对该安全系数的定义应该这样理解：保持平衡并不需要发挥全部的强度，只需要发挥一部分τ，发挥了这样大的强度整个土坡达到极限平衡。

所谓保持平衡所需要的强度就是实际发挥的强度，或者说，若滑面上各点的抗剪强度为τ=τf/Fs，则土坡达到极限平衡，当然此时各土条也处于极限平衡状态。

安全系数是对整个土坡而言的，对各土条均取这一相同值，意味着假定滑动体各部分强度发挥程度一致。

边坡形式及土条受力如图1所示。

安全系数公式[3]为
式中：b为土条宽度；c为粘聚力；φ为内摩擦角；α为滑弧线切线与水平线夹角；Y为竖向条间力；E为横向条间力；W为土条重力；N为滑面法向力；T为滑面剪切力；其中
（6）式就是广泛应用的毕肖普法。

2.2 关于圆弧形滑裂面
圆弧形滑动面形式最为简单，且建立了一套完整的工程指标和经验指标，所以应用也最为广泛。

对于均质的各向同性的黏性土质，它是一个比较好的近似。

毕肖普法并非通用条分法，所以其只适用于圆弧滑裂面这种假定。

假定有圆弧滑动面y=y(a,b,R),其中圆心坐标为（a,b）,圆的半径为R。

研究表明，安全系数是圆心坐标和半径的函数[4]。

于是，按照毕肖普法得到最终的安全系数
的表达式将是一个以[a,b,R,F]为自变量的函数，然后可用最优化方法求解最小安全系数。

当考虑其他实际因素时，例如孔隙水压力、多层土质、地震以及非线性强度关系等，
则要对安全系数的表达式做相应的变化，但总体的求解思路是不变的。

MATLAB语言是当今国际上科学计算领域最具影响力也最具活力的软件，它提供了强大的优化方法工具箱，其中的优化程序使用的是一些目前公认的有效算法。

3.1 3种最优化方法的同时实现
假设已经按照毕肖普的方法得到了安全系数的表达式：F=G(a,b,R,F),下面的代码会帮助同时实现BFGS法、DFP法、单纯性法求最小安全系数，另外还可以从输出中得知迭代次数，收敛情况等信息。

3.2 实例计算及方法、结果分析
1986年澳大利亚计算机应用协会（ACADS）对澳大利亚所使用的边坡稳定分析程序进行了一次调查。

笔者就以此次调查设计的考核题中的一个为例进行计算。

考核题为一个均质边坡，如图2所示。

材料特性如表1所示。

将安全系数表达成圆心坐标和半径的函数可参见文献[4,6]，以其中一个随机产生的初值x0= [35,40,18]为例演示,其中F=1作为第一个迭代值，每次以新求得的安全系数进行迭代，3～4次即可得到满意结果。

设想的随机数产生代码为
“a=unifrnd(m,n,1,k),b=unifrnd(m,n,1,k),R=unifrnd(m, n,1,k)”，k为产生的随机数的个数，[m,n]为3个自变量的可能取值范围，据实际情况自行设定。

此代码是按照均匀分布产生随机数，按照最优化理论，用位于同一个“凹”区内的随机数求得的极值是一样的，否则就可能有差异，这样做的目的就是确保全局最小值。

对于例题这种简单边坡，潘家铮认为其圆心的大体位置可这样确定：过边坡中点，分别以L/2和3L/4为半径画弧，与边坡中法线和中垂线形成一个范围，圆心即在其中[7]。

当然，确定了此范围也可使用有约束最优化进行求解，MATLAB中使用的是基于内部映射牛顿法的子空间置信域法、序列二次规划法以及BFGS法，结果是一样的。

其计算结果如表2所示。

相应的危险滑动面图2所示。

通过计算得知，3种方法都成功收敛且结果相差很小，但是DFP法得到的滑裂面是偏于危险的，另两种方法得到的滑裂面是偏于安全的。

另外，单纯搜索的迭代次数要远大于拟牛顿法，前者超过100次，后者差不多是其半数。

尽管存在这些差别，但还是尽量同时使用它们以便取长补短，确定最危险滑裂面，一法定安全是不可取的。

对于随机产生的不同的迭代初值中，除个别明显不合理的外，在3种方法中所得的最小安全系数值是相同的。

如果条件允许的话，通过瑞典法再计算一遍以便比较，这里用单纯形法的计算结果如下：
然后以上面初值进行毕肖普法的上述3种最优化计算，其结果非常接近表2中单纯搜索的结果，由此可知单纯搜索的结果是全局最小值的可能性最大。

对于这种简单边坡，瑞典法的结果往往是偏小的[8]。

综上所述，要选DFP法所计算的0.998 5，据此对相应的危险滑裂面进行加固。

对于此调查题，各个科研机构提供的标准答案是：1.000和0.990。

另外此题是比较简单的情况，旨在说明最优化方法和MATLAB在这方面的应用。

结合例题一方面产生均匀分布的随机数作为初值，另一方面用相对更简单的瑞典条分法计算结果作为初值，对毕肖普法优化计算得到了最小安全系数。

针对岩土边坡问题的复杂性，文章主导的思想是尽量在各个环节采取多途径进行电算，以确保全局最小值。

笔者编程水平有限，相信就此问题而言还可以有更方便、可靠、智能的程序设计。

滑坡灾害的研究有近百年的历史，有关滑坡方面的研究成果和文献资料数不胜数，然而滑坡灾害及其造成的损失却让人失望的与日俱增[9]。

殷切希望各位有志之士从各个角度致力于边坡研究，保护人类生命及财产。

【相关文献】
[1] 陈祖煜.土质边坡稳定分析——原理,方法,程序[M].北京:中国水利水电出版社,2003:87-97.
[2] 郭科,陈聆,魏友华.最优化方法及其应用[M].北京:高等教育出版社,2007:86-115.
[3] 殷宗泽.土工原理[M].北京:中国水利水电出版社,2007: 410-423.
[4] 张天宝.土坡稳定分析圆弧法的数值解研究[J].成都工学院学报,1978(1):99-101.
[5] 苏金明,张莲花,刘波,等.MATLAB工具箱应用[M].北京:电子工业出版社,2004:190-202.
[6] 杨庚宇.土坡稳定分析中条分法的解析计算[J].力学与实践,1995,17(2):59-60.
[7] 钱家欢,殷宗泽.土工原理与计算[M].北京:中国水利水电出版社,1996:308.
[8] 方玉树.边坡稳定性分析条分法最小解研究[J].岩土工程学报,2008,30(3):333-334.
[9] 张鲁渝,欧阳小秀,郑颖人.国内岩土边坡稳定分析软件面临的问题及几点思考[J].岩石力学及工程学报,2003, 22(1):166-169.。