基于改进狼群算法的边坡极限平衡稳定性分析方法研究

马金城 石 崇 张一平 张文浩 王乐荣

(1.河海大学 岩土力学与堤坝工程教育部重点实验室,南京 210098;2.河海大学 岩土工程科学研究所,南京 210098)

边坡稳定性分析是岩土工程十分重要的研究内容,涉及水利水电、公路、铁路、矿山工程等诸多领域,是判断边坡稳定性,是否需要加固及采取何种防护措施的主要依据.在进行边坡稳定分析时,关键是如何得出滑面并利用滑动力、抗滑力计算安全系数、优化得出最危险滑面.因此,对于危险滑面搜索方法进行理论和应用研究非常有价值.

在常规的条分法边坡稳定计算中,如瑞典条分法和毕肖普法假定滑坡体为不变形的刚体,利用自重计算滑面上滑动力与抗滑力,再利用抗滑力矩与滑动力矩之比来定义稳定性系数[1].而自然岩土体处于弹塑性变形状态,达到屈服时应力可向周围传递,因此条分法得到的各条分底部抗滑力、滑动力往往与实际偏差很大.而数值模拟方法综合考虑物理方程、平衡方程、几何边界条件,因此可以逼真地得到岩土体应力场,但不能直接给出稳定性系数.如果将这两者有效结合起来,利用数值模拟得到真实的应力场,采用矢量和法可定义稳定性系数,对边坡稳定性计算有很好参考价值.因此近些年很多学者致力于探索危险滑面搜索方法的工作,如遗传算法[2],蚁群算法[3],粒子群算法[4],模拟退火法[5]等.这些算法在边坡滑面搜索中得以运用,与早期使用的穷举法、二分法相比,上述方法的精度和收敛性提高了很多,对边坡稳定研究和实际工程应用起到了积极的作用,但蚁群算法容易出现停滞等现象[6],粒子群算法会过早陷入局部最优等缺点[7],模拟退火算法和遗传算法对参数选择要求较高,如若参数选择不当会使搜索陷入局部极值或搜索结果不令人满意[8].

狼群算法[9]是2011年提出的一种新型的智能优化算法,该算法在全局搜索和局部开发能力上具有良好的适用性,被广泛地运用到路径规划等实际问题中,具有鲁棒性好,全局寻优能力强,收敛速度快等优点,但目前尚未运用到边坡稳定性分析之中.本文基于传统的狼群算法,尝试结合边坡稳定性计算特点进行优化改进,建立狼群滑面位置搜索及稳定计算模型,并用已有成果验证其可行性和准确性.

1 极限平衡稳定性分析基本原理

1.1 抗滑力和下滑力的计算

对一个二维边坡中潜在的滑面,计算稳定时一般如图1所示将滑面划分为若干段,不同段相交点称为一个状态点,从起始点到终止点通过状态点首尾相连可形成一条路径,每条路径都是一条潜在的滑动面.

图1 滑面受力示意图及单元体图

设滑动方向向右,逆时针为下滑力正方向,则其中某一段上的抗剪强度τf和切向应力τ的计算公式如下：

第i段的平均滑动力：

第i段的平均抗滑力：

式中：σx、σy、τxy分别为单元形心沿x,y方向的正应力和切应力；φ和c分别为材料的内摩擦角和内聚力；σα为斜面的法向应力；α为斜面与水平面的夹角.

1.2 安全系数的定义

利用数值模拟的方法获取边坡真实应力场后,计算出每一段上的滑动力和抗滑力,然后采用矢量和法定义安全系数,安全系数定义为：

式中：A和B分别为边坡起始点和终止点；n为滑动面所经过的段的数量；τ为沿滑动面滑动方向的剪应力；dL为沿滑动面每一段的长度增量；τi为从第i-1个节点到第i个节点路径上的剪应力；τfi为从第i-1段到第i段路径上的抗剪强度；Δd i为从第i-1阶段到第i阶段路径的长度[10].

2 滑面搜索的改进狼群算法

对于已知滑面边坡,采用上述极限平衡稳定性分析方法可直接求出安全系数,但对于滑面未知的边坡,如何确定滑动面上点的位置,找到最危险滑面是非常困难的[11].从数学问题看,确定滑面实际上是图形路径的优化问题.近年来,狼群算法作为一种生物仿生算法被提出,可用于图形最优化问题的求解.因此,本文针对二维边坡的网格,将滑出滑入点当做头狼与母狼的位置,某点的应力当作当前位置的猎物气味浓度,网格节点当作猛狼,网格内部的某些点为探狼,通过不同种类狼之间的信息交互、“胜者为王”的头狼产生机制和“强者生存”的狼群更新机制来寻求全局路径的最优解,可以用于解决边坡最危险滑面的搜索.

2.1 经典狼群算法基本原理

经典狼群算法分为如下几个步骤：

1)头狼生成：在某一待寻优区域内,将位置最优的狼视为头狼.头狼是狼群的核心,负责狼群的指挥工作[12].

2)游走行为：选取除头狼外最佳的数匹狼作为探狼,执行游走行为,根据某种计算规则得到探狼处气味浓度,然后和当前头狼气味浓度进行比较,选择二者中较优者为头狼.

3)召唤行为：当头狼发出召唤行为后,猛狼以较大步长向头狼位置逼近.如果在奔袭过程中,猛狼位置优于头狼,则猛狼代替头狼发出召唤行为,否则猛狼继续奔袭,直至距离头狼位置小于围攻距离,此时转化为围攻行为[13].

4)围攻行为：当猛狼奔袭到距头狼较近位置时,猛狼联合探狼对猎物进行围攻.

5)“强者生存”的狼群更新机制：剔除掉较差的狼,只有适应性较好的狼才能生存下来.

2.2 狼群算法的改进

1)模型的建立：根据边坡形状及地层分布状况,结合边界条件,对边坡进行合理的有限单元网格划分,单元格数量由上述因素及所需结果精度共同确定.得到有限单元网格后,计算边坡的应力场,得到每个节点处的应力和坐标[14].将所得应力场视为草原平面,也即狼群活动的平面,将节点处的应力值视为当前位置的猎物气味浓度.

2)头狼生成：在狼群中,头狼和母狼处于支配地位,头狼负责狼群的指挥,母狼负责配合头狼工作.头狼和母狼对整个狼群的优劣起决定作用.如图2 所示,头狼和母狼在可能滑入点与滑出点的某一区域内,初始头狼与母狼分别为该区域内横坐标值最小的点.算法中引入母狼的原因是为配合头狼工作,两者共同把握整个狼群的位置,避免陷入局部最优.

图2 可能圆心位置与滑出滑入点切线斜率

3)狼群区域划分：如图3 所示,以滑入点(X0,Y0)为坐标原点建立搜索区域坐标系,用相同间距为x d=(X n-X0)/n的纵向线将待寻优区域划分成n部分,其中X n为滑出点横坐标.将每条纵向线邻域视为一个搜索区域,分别记为LQ1,LQ2,…,LQn.将各纵向线横坐标依次记为X1=(X n-X0)/n,X2=2(X n-X0)/n,…,X n=X n.在给定滑出滑入点位置后,需确定下一个点(X ML1,Y ML1)的搜索位置,潘家铮等人通过研究表明,对于简单均质土坡,作边坡中点的垂线和法线,以坡面中点为圆心,以1/4坡长和5/4坡长为半径画同心圆,最危险滑弧的圆心在如图2所示阴影部分ABDC范围内.据经验可知,圆心纵坐标越小,滑出点和滑入点与圆弧相切切线斜率越大(滑出点和滑入点切线的斜率仅用于确定滑出滑入点的下一个点搜索范围,并不对滑动面是否为圆弧面做限制).设滑出点处的切线最大斜率绝对值为k1,滑入点处最大斜率绝对值为k2.

图3 狼群搜索区域示意图

如图4(a)所示,LQ1的搜索范围为：

如图4(b)所示,在LQ i中,横坐标在i(X n-X0)/n的左右邻域内波动.

图4 狼群搜索区域及探索步长

由于滑动面点为一下凹曲线,因此点(X0,Y0)与点(X ML1,Y ML1)所在直线为LQ2搜索区域的下限,X=X1-δ,X=X1+δ是X=X1左右邻域.点(X0,Y0)与点(X n,Y n)所在直线为LQ2搜索区域上限.也即LQ2的搜索范围为：

根据上述规则,由前一搜索区域所得猛狼位置来推算下一狼群搜索区域,直至最后一匹猛狼产生.这样做可以有效缩小搜索区域的范围,提高搜寻效率,减少搜索时间.

4)游走行为：探狼从每个待寻优区域最上侧的点开始以固定步长y d向下游走探索.处于头狼和母狼附近的探狼分别以固定步长x d水平向右游走探索.

5)召唤行为：当头狼和母狼及探狼已处于各自狩猎位置后,头狼发出召唤行为,召唤探狼周围某一范围内猛狼.视探狼所处单元格内部为召唤区域.处于单元格节点处的猛狼在接收到召唤指令后,向附近区域探狼奔袭,通过召唤区域上猛狼处气味浓度综合确定探狼处气味浓度,也即由已知多边形节点的坐标和应力通过机器学习的方法确定待求点的应力状态(具体确定过程见(8)).

6)围攻行为：可能最优路径上的头狼和猛狼经上述方式产生后,对猎物展开围攻,形成一狼群.

7)“强者生存”“优胜劣汰”的狼群更新机制：将适应度D作为目标函数,则狼群算法可将求最小安全系数的问题转换为求D最大值问题,目标函数定义为：

适应度D为狼群对草原适应性的好坏,D越大,则狼群越优.在完成每次可能最优路径适应度D i计算后,将其适应度D i与当前适应度最好的Dlead进行比较,直至最后一条路径计算结束,将适应度最好的狼群视为优胜狼群,也即：

其中：Dmax为全局最优解,其对应的各种路径为最适宜路径.

8)机器学习方法确定猛狼处气味浓度：在Python 中使用sklearn库中的Gradient Boosting Re-gressor模块,将给定数据(结点坐标和应力值)分为训练集和测试集[15](将部分已知节点应力和坐标作为训练集,其余节点坐标和应力作为测试集,在训练集中找寻规律,在测试集中验证规律[16].可采用L1范数损失函数(loss function)来估量模型的预测值f(x)与真实值Y的不一致程度,其值的大小可以用来评价预测的准确性,该函数把目标值(Y)与估计值(f(x))的绝对差值的总和最小化[17]：

相比于传统的插值法求中间部分应力,机器学习方法可以很直观地通过S值的大小来判断对某点预测的准确性.

3 程序实现与算例验证

根据以上原理,形成改进狼群算法边坡滑面搜索方法,计算流程图如图5所示.

图5 改进狼群算法流程图

本文采用python语言编制程序,引入一经典算例来验证狼群算法在边坡滑面搜索以及计算安全系数的准确性.该算例为澳大利亚计算机协会于1987年委托BDonald和PGiam 设计ACADS边坡稳定性分析程序的考核题,该题的参考滑动面如图6(a)所示,推荐答案的安全系数为1.00[18],具体参数为：容重γ=20 k N/m3,泊松比μ=0.25,黏聚力c=3.0 kPa,内摩擦角φ=19.6°.由于该算例已经过大量学者的引用和对比分析验证,因此其结果视为标准值.

基于改进狼群算法理论的边坡临界滑动面搜索方法搜索出最危险滑动面如图6(b)所示,对应的安全系数为1.03675,计算所得安全系数与推荐值1.00误差约3.6%,采用本文方法得到的临界滑动面与推荐滑动面基本一致,这表明基于改进狼群算法的边坡稳定性分析是可行的,其值偏高的原因是数值应力结果比刚体极限平衡协调,因此有更高的承载能力.

图6 边坡滑面(单位：m)

4 工程案例研究

溪洛渡水电站左岸谷肩堆积体的典型地质剖面如图7(a)所示.本文采用flac-3D 软件建立如图7(b)所示数值模型,选用摩尔-库仑的本构模型,以及表1中岩土物理力学参数进行数值计算,得到其应力场后,采用Python语言编写程序,运用改进狼群算法搜索其最危险滑面,并计算其安全系数.

表1 岩土物理力学参数

如图7(a)所示,IN13钻孔倾斜仪与剖面较为接近,图7(c)表明,溪洛渡水电站左岸谷肩堆积体的最危险滑面位于滑坡堆积体与宣威组上部交界处.这与历史监测资料所揭示的滑动面位置一致.文献[19]在自然工况下采用改进动态规划算法,圆弧形滑面计算所得安全系数为1.100,而在相同工况下采用狼群算法得到的安全系数为1.044,如图7(d)所示.经对比滑面可知,狼群算法搜索得到的滑面可保证更多比例滑面位于宣威组,而现场勘查表明该边坡处于蠕滑状态,自然工况安全系数为1.00～1.05之间,计算结果接近蠕滑状态自然工况下安全系数的上限1.05,故采用本方法得到的结果更加接近于实际.

图7 工程案例计算结果

5 结论

基于数值计算得到的边坡应力场,提出一种基于改进狼群算法的边坡临界滑动面搜索方法,为边坡稳定分析提供一条新的有效思路.得到主要结论如下：

1)改进狼群算法是利用数值模拟手段获得真实的应力场后,采用矢量和法来定义安全系数.在边坡稳定性分析过程中可得到较为准确的边坡最危险滑动面和最小安全系数.

2)在运用改进狼群算法的过程中,用机器学习方法代替传统的插值法来计算指定坐标点处的应力,可以通过损失函数直观地看出计算结果的准确性.

3)通过对经典边坡算例的分析,本文所提方法的计算结果与典型算例参考答案具有较好的一致性,证明了该算法的合理性.