引入和声策略的果蝇算法求解边坡最小安全系数

李 冉，曾 宇，于广明，李 亮

(1. 青岛理工大学 土木工程学院，山东 青岛 266033；2. 云南省交通规划设计研究院有限公司，云南 昆明 650041)

边坡稳定性分析是岩土工程中常见的课题之一。用于评价边坡稳定性的比较完善的方法共有3种：第1种方法是极限平衡方法。该方法视滑动土体为刚体，不考虑土体的应力-应变关系，假定边坡可能的滑动面，然后利用该滑动土体失稳时的静力平衡条件求解边坡沿该滑动面的安全系数，利用各种优化算法变换不同的滑动面，并从中选择安全系数最小的滑动面作为设计滑动面[1-21]。第2种方法是数值分析方法，主要包括有限单元法和有限差分法。数值分析方法考虑复杂的边坡土层情况，考虑土体的应力-应变关系，通过与强度折减法结合，自动判断边坡最危险的滑动区域和相应的滑动安全系数[22-23]。第3种方法是极限分析方法。该方法基于塑性力学中极限分析的上、下限定理发展而来。文献[24]中结合上限定理发展了极限分析上限法，并基于虚功原理，推导了斜条分模式下边坡安全系数的计算方法。第1种方法虽然引入各种假定，无严格理论基础，但是计算简单，工程经验丰富，在土木工程各相关规范中得到了广泛应用。

目前针对极限平衡方法的研究热点集中在最小安全系数的搜索方面。自20世纪90年代以来，国内外学者尝试应用遗传算法、粒子群算法、蜂群算法、模拟退火算法、禁忌算法、鱼群算法、和声算法、萤火虫算法等启发式优化算法进行边坡稳定最小安全系数的搜索。针对复杂程度不同的边坡稳定分析问题，各种算法表现不同，每种算法各有优势和劣势，因此在已有算法基础上，融合2种或多种算法提出混合搜索算法具有较强的研究意义和应用价值。本文中结合和声搜索算法和果蝇算法，提出引入和声策略的果蝇算法，通过复杂边坡最小安全系数优化求解问题，对比分析基本果蝇算法和改进的果蝇算法的计算效率。

1 边坡稳定性优化模型

边坡稳定性优化模型示意图如图1所示。边坡稳定性优化模型为

minfs(g,p,f,s)，

g—边坡的几何形状；ρ—边坡土层的密度；c—边坡土层的黏聚力；φ—边坡土层的内摩擦角；f—边坡承受的荷载，如自重、堆载等；s—边坡滑动面。图1 边坡稳定性优化模型示意图

式中：fs为安全系数；g为边坡的几何形状；p为边坡土层的材料参数，包括密度ρ、内摩擦角φ与黏聚力c；f为边坡承受的荷载，如自重、堆载等；s为边坡滑动面。

边坡最小安全系数的确定原理如下：在给定边坡几何形状g、土层材料参数p、荷载f的条件下，不断变换可能的失稳模式s，从而找到最小的安全系数及其对应的滑动面。在极限平衡方法中，滑动面可以假定为圆弧或非圆弧，不失一般性地，本文中假定滑动面为非圆弧，由满足一定规则的点通过样条函数插值得到，限于篇幅，非圆弧滑动面的模拟步骤不再赘述，具体可参阅相关文献[15-16]。对于给定的非圆弧滑动面s，采用Spencer法[17]计算安全系数fs。

2 基本果蝇算法

受果蝇寻食行为启发，潘文超[25]于2011年提出基本果蝇算法。基本果蝇算法主要通过不同果蝇发现的食物浓度调整果蝇的飞行方向和距离，具有参数少、简单且易实现等优点。此外，相对于其他物种，果蝇具有敏锐的嗅觉与视觉，能够很好地搜寻到空气里漂浮的气味，甚至能够嗅到40 km以外的食物源。靠近食物位置后，能凭借视觉的敏锐度找到食物和同伴聚集的位置，然后飞向目标源。图2所示为基本果蝇算法的计算流程。基本果蝇算法的参数包括果蝇个数N、最大飞行距离L、最大迭代次数t。在算法迭代开始前，首先在搜索范围内随机产生N个果蝇个体，位置记为Xi=(xi1,xi2,…,xin)，i=1,2,…,N，n为正整数。x1,x2,…,xn为优化变量个数。在边坡稳定分析中，为了求解最小安全系数对应的滑动面，优化变量通常包括模拟滑动面所需的一系列变量。

图2 基本果蝇算法的计算流程

为了便于实现，令最大飞行距离为向量，即L=(l1,l2,…,ln)，其中li为在第i维度上的最大飞行尺度。果蝇位置的味道浓度T(Xi)类似于边坡稳定分析中的安全系数fs(Xi)。安全系数越小，该果蝇发现的食物浓度越大，对其他果蝇的吸引力越大。在基本果蝇算法中，当前果蝇群体中味道浓度最大的果蝇位置具有特别强的吸引力，其他果蝇个体均飞向该个体，然后在此基础上继续依靠嗅觉和视觉找寻味道浓度更大的位置，直至迭代次数达到指定的次数。基本果蝇算法的优势在于具有较强的局部寻优能力，在当前味道浓度最大的位置附近进行寻优；缺点在于忽视了果蝇个体在朝最大浓度位置飞行过程中的自主选择性。为此，本文中提出引入和声策略的果蝇算法，通过算例分析，研究果蝇个数N、最大飞行距离L对计算结果的敏感性，从而为基本果蝇算法的应用提供借鉴。

3 引入和声策略的果蝇算法

在基本果蝇算法中，果蝇发现浓度最大的个体后，会直接飞向该最优个体。实际上，果蝇在向最优个体飞行的过程中，仍然会保留原来个体的部分信息，并根据自己的嗅觉和视觉综合判断朝最优个体飞行的程度。鉴于此，本文中引入和声策略[4]模拟该飞行过程。设当前最优果蝇个体位置为Xbest=(xb1,xb2,…,xbn)，拟飞向Xbest的个体为Xi=(xi1,xi2,…,xin)，利用和声策略，确定第i只果蝇下一次飞行到达的位置。和声策略流程如图3所示。

图3 和声策略流程

在引入和声策略后，改进的果蝇算法计算步骤如下：

1)设定果蝇个数N，和声策略参数hr、pr，最大迭代次数t；

2)随机生成N个果蝇位置X1，X2，…，XN；

3)利用Spence法[17]计算N个果蝇位置X1，X2，…，XN相应的安全系数fs(X1)，fs(X2)，…，fs(XN)；

4)挑选安全系数最小的果蝇位置，记为Xbest；

5)利用图3所示的和声策略产生第i个果蝇位置，i=1,2，…，N；

6)迭代次数加1，转步骤3)继续进行搜索，直至迭代次数达到t。

4 算例分析

图4所示考虑空间变异特性的边坡。坡高为10 m，考虑不排水工况下边坡的短期稳定性，强度参数取为不排水强度cu，即考虑黏聚力与内摩擦角的综合强度。假设cu服从对数正态分布，均值为40 kPa，变异系数为30%，即标准差为12 kPa，cu的空间变异特性用一维(垂直)随机场模拟，采用厚度为1 m的水平土层离散整个边坡。该均质边坡被离散为12个水平的小土层，即不排水强度为cu(1)～cu(12)。在每个小水平土层内，土体参数完全相关。在不同的小水平土层之间，土体材料参数不同，但是具有一定的相关性。根据随机场理论，相关性随着小水平土层之间距离的增大而减小，例如，cu(1)与cu(2)之间的相关性比cu(1)与cu(12)之间的相关性强。相关性的计算公式为

(1)

cu(i)—第i水平土层的不排水强度，i=1,2,…，12。图4 考虑空间变异特性的边坡

式中：ρcu(i)cu(j)为cu(i)与cu(j)之间的相关系数；yi为cu(i)小水平土层中心的y坐标；yj为小水平土层cu(j)中心的y坐标；λ为波动范围，是描述土层空间变异特性的参数。λ越大，则不同水平土层之间的相关系数越大。当λ→+∞时，不同水平土层之间的相关系数均为1。

在得到相关系数矩阵后，可以利用Cholesky分解生成符合指定概率密度函数的一系列样本。每个样本实际上是指cu(1)～cu(12)共12个cu值。本文中的算例即针对随机生成的10个样本展开。表1所示为波动范围为1 m时，随机生成的10个样本值。本文中采用连续7个点的直线连接来模拟非圆弧滑动面，因此优化变量个数为7，即n=7。

表1 波动范围为1 m时生成的10个样本值 kPa

4.1 基本果蝇算法参数的敏感性分析

选取果蝇群体数N=10，20，30，40，50，最大飞行距离L=(10,10,1.0,1.0,1.0,1.0,1.0)，最大迭代次数t=500，针对每个样本值进行10次计算，得到10个最小安全系数值，利用10个安全系数值的平均值作为该样本条件下基本果蝇算法的结果。图5(a)所示为不同群体个数时，10个样本条件下10次计算的平均值。由图可知，对于某个样本，随着群体个数N的增加，基本果蝇算法的寻优能力逐渐增强，即安全系数平均值逐渐减小。图5(b)所示为不同群体个数时，10个样本条件下10次计算结果的标准差。由图可知，不同群体个数时，10次计算结果的标准差呈现不同的规律。当群体个数较小时(N=10)，对于第3样本，标准差最小，但是安全系数平均值最大。综合来看，当群体个数较小时，算法寻优效率低。由此可知，当果蝇群体个数取大于50时，可以获得较好的结果。

(a)安全系数平均值

算法参数取N=50，t=500，最大飞行距离分别取L=(10,10,1.0,1.0,1.0,1.0,1.0)、L=(5,5,1.0,1.0,1.0,1.0,1.0)和L=(1,1,1.0,1.0,1.0,1.0,1.0)，研究最大飞行距离对算法参数的影响。为了便于叙述，以上3种飞行距离依次称为飞行距离1、2、3。图6所示为不同飞行距离时安全系数平均值的变化。由图可知，当样本个数为10时，飞行距离3对其中6个样本给出了最小的安全系数平均值，飞行距离1对其中4个样本给出了较好的结果。由此可见，针对不同样本情况，需要采取不同的飞行距离进行计算。

图6 基本果蝇算法最大飞行距离L的敏感性分析

4.2 引入和声策略的果蝇算法敏感性分析

在引入和声策略的果蝇算法中，首先设定果蝇个数N=10，和声策略参数pr=0.1，分别选取和声策略参数hr=1.0，0.95,0.85，研究hr对计算结果的影响。图7(a)所示为10个样本时，不同hr参数条件下，安全系数平均值的变化。为了便于比较，图7(a)中还给出了基本果蝇算法在N=10条件下，飞行距离1时的结果。由图可知，当样本个数为9时，改进的果蝇算法的结果均优于基本果蝇算法的。不同的hr会导致不同的结果。hr取值越大，安全系数平均值越小，算法的寻优效率越高。建议hr取值大于0.95。

(a)声策略参数hr

设定hr=0.95，pr=0.1，分别研究N=10,20,30,40,50时，10个样本条件下的安全系数平均值变化情况。图7(b)所示为不同果蝇群体个数N对计算结果的影响。由图可知，与基本果蝇算法类似，随着N的增加，安全系数的平均值逐渐减小，说明改进的果蝇算法中群体个数不应取得过小，建议取为大于50。

4.3 2种果蝇算法对比分析

为了验证本文中改进的果蝇算法的优势，对2种果蝇算法进行对比分析。在基本果蝇算法中，取飞行距离为1。在改进的果蝇算法中，取hr=0.95，pr=0.1，分别对比N=10,20,30,40,50时，2种果蝇算法的寻优能力，结果如图8所示。由图可以看出，在N相同时，改进的果蝇算法的寻优能力均优于基本果蝇算法的，表明引入和声策略对基本果蝇算法的飞行起到了良好的促进作用，也证明了果蝇个体在飞向浓度最大的位置时，有选择保留原来位置的趋势，这种趋势会使群体获得更好的食物。

5 结论

在边坡稳定分析中，快速求解最小安全系数具有较强的工程意义和研究价值，尤其对于考虑空间变异特性的边坡而言更是如此。本文中综合和声策略和基本果蝇算法，提出改进的果蝇算法并进行边坡稳定分析，研究了基本果蝇算法中群体个数N和最大飞行距离L对计算结果的影响，其次还分析了和声策略参数hr、pr对计算结果的影响，得到以下主要结论：

1)基本果蝇算法和改进的果蝇算法中果蝇个数N对计算结果影响较大，建议N取为大于50。

2)最大飞行距离L对计算结果影响较复杂，针对不同问题建议采用不同飞行距离进行试算，最终找到最优结果。

3)在改进的果蝇算法中，hr应尽量取较大值，一般建议取为大于0.95，pr建议取为0.1。