基于模拟退火的果蝇优化算法的边坡临界滑面搜索

余秀玲， 申翃， 郑万栋

(武汉理工大学 土木工程与建筑学院， 湖北 武汉 430070)

边坡稳定性分析是岩土工程中的重要研究课题之一，其研究方法包括：工程地质类比法、极限平衡法、极限分析法、数值分析法等。其中极限平衡法的应用最为广泛，但是其必须依赖研究者的经验，试算一系列滑面寻找最小安全系数及其对应的临界滑面位置。

近年来，很多学者致力于临界滑面自动搜索技术的研究，提出了各种适用于圆弧滑面或非圆弧滑面的搜索方法，如蒙特卡罗法与有限元结合、临界化方法、遗传进化算法、模式搜索法、相遇蚁群算法、神经网络与多元回归等。与早期使用的穷举法、二分法相比，上述方法的精度和收敛性提高了很多，对边坡稳定研究和实际工程应用起到了积极的作用，但有时也会发生计算成果并非整体安全系数极值的情况。该文将果蝇优化算法和模拟退火算法相融合，在Matlab中编程实现基于模拟退火机制的果蝇优化算法，通过两个土坡算例验证，并与其他方法进行对比分析。以期为边坡临界滑面搜索提供新途径。

1 边坡临界滑面搜索理论模型

简化Bishop法是目前常用的一种滑坡稳定性分析方法，被认为是圆弧滑面计算分析的最佳极限平衡法。其力学模型如图1所示。假定条间力的方向水平，切向剪切力为零，根据滑面上土体的静力平衡和整体的力矩平衡可以得到滑面的安全系数：

Fs=

(1)

式中：Fs为对应滑面的安全系数；n为条分数；ci、φi分别为滑动体第i条块土层的黏聚力和内摩擦角；Wi为第i条块自重；bi为第i条块的宽度；θi为第i条块中心线与圆弧滑面交点处圆弧切线的水平角。

图1 简化Bishop法力学模型

若忽略条间力的影响，即为瑞典条分法，安全系数的表达式简化为：

(2)

式中：L为圆弧滑面的弧长。其他参数同式(1)。

式(1)为隐式函数，求解安全系数时需要迭代，式(2)为显式函数，可以直接求解。

当边坡的几何参数和力学参数已知时，安全系数是滑弧圆心坐标和半径的函数，即：

Fs=f(x0,y0,R)

(3)

式中：x0、y0、R分别为圆弧滑面的圆心坐标和半径；f为安全系数Fs与可行搜索域参数的函数。

则式(3)转化为一个优化问题：

(4)

式中:x(j)={x(1)=x0,x(2)=y0,x(3)=R}为优化变量；[a(j),b(j) ]为第j个变量的取值区间。

2 基于模拟退火的果蝇优化算法

2.1 基本果蝇优化算法(FOA)

基本果蝇优化算法(FOA)由台湾学者潘文超基于果蝇的觅食行为而开发，是一种新型的群体智能全局优化算法，近年来已被研究人员广泛运用到诸多领域解决优化问题。果蝇优化算法机理是基于果蝇种群觅食过程的模拟，利用果蝇良好的嗅觉和视觉，设计相应的感官搜索环节，通过不断迭代从而获得最优解。果蝇群体寻找食物的迭代过程如图2所示。

果蝇优化算法步骤如下：

(1) 果蝇种群相关参数的初始设定：位置X_axis、Y_axis，规模Sizepop，最大迭代次数Maxgen。

图2 果蝇群体觅食迭代过程

(2) 更新果蝇个体觅食过程的随机距离(RandomValue)与方位。

(5)

(3) 估量果蝇种群中个体与原点之间的距离。

(6)

(4) 将味道浓度判定值Si代入味道浓度判定函数，即适应度函数(Fitness function)，得到果蝇个体的味道浓度Smell(i)。

Smell(i)=Function(Si)

(7)

(5) 在果蝇种群中找出味道浓度最佳的个体。

[bestSmellbestindex]=min[Smell(i)]

(8)

(6) 记录并保存最佳味道浓度值bestSmell及其位置X和Y，此时其他果蝇飞向该位置。

(7) 重复步骤(2)～(5)，此次与前一次的最佳味道浓度做比较，若优于前一次的味道浓度，并且迭代次数未达到Maxgen，则执行步骤(6)，否则，结束算法。

基本果蝇优化算法能够快速高效地搜索出最优解。但是，由于在利用适应度函数的寻优过程中极易陷入局部最优解，导致果蝇种群停止迭代更新位置，因此得到的最优解具有局限性。

2.2 模拟退火算法(SA)

模拟退火算法是根据物理学中固体退火原理而建立的一种全局最优算法。固体加热升温时内部粒子排序紊乱，内能逐渐增大，降温时，粒子排列趋于有序，内能减小。假设在状态A时，系统内能为E(A)，受到外界扰动后，变为状态C，内能变为E(C)。基于Metropolis准则，粒子在温度T趋于平衡的概率为:

(9)

由式(9)可以看出：退火时，新状态使内能函数减小，系统一定接受该状态；新状态使内能函数增加时，系统以一定概率接受该状态，即在优化进程中，可以获得最优解和一定概率的差解。将内能函数E(X)演变为目标函数，退火温度T演变为控制参数，随着T值的变化，算法执行“新解-判断-舍弃或接受”的进化，直至求得问题的最优解。

模拟退火算法的步骤如下：

(1) 随机赋予初始状态，设置初始解x0、初始温度T0以及最大迭代次数Maxgen，计算目标函数初始值E(x0)。

(2) 判断所赋初始值是否满足终止条件，若是，该解即为最优解，结束算法；否则，在x0邻域内随机选取x*重新赋值，计算内能增量ΔE=E(x*)-E(x0)。

(3) 若ΔE<0，则x=x*；否则，以概率exp(-ΔE/T)作为x*的解。

(4) 进入迭代寻优过程，重复步骤(2)和(3)，直至寻得最优解或达到最大迭代次数。

2.3 基于模拟退火的果蝇优化算法(SA-FOA)

如前所述，基本果蝇优化算法(FOA)参数少，易调节，效率高，却易陷入局部最优，而模拟退火算法(SA)可以有效地跳出局部最优，因此该文将两种算法融合，形成基于模拟退火的改进果蝇优化算法。其步骤为：

(1) 执行基本果蝇算法的步骤(1)～(6)。

(2) 模拟退火步骤，确定初始温度T0

T0=Smellbest/log(5)

(10)

(3) 计算当下温度状态果蝇味道浓度判定值Si的适应值：

TF(Si)=

(11)

(4) 计算随机概率:

p=rand()

(12)

(5) 采用轮盘赌法策略，从Si中计算得到全局最佳的某个值(bestX，bestY)作为下一代果蝇的搜索中心位置。

(13)

(7) 计算果蝇个体的味道浓度Smelli*，若Smelli*

(14)

(8) 进行退温操作：

T=0.5T0

(15)

(9) 进入迭代寻优，重复执行(3)～(9)，直至寻得最优解或达到最大迭代次数。

将基于模拟退火机制的果蝇优化算法(SA-FOA)应用于边坡的临界滑面搜索，其流程如图3所示。

待优化的圆弧滑面圆心坐标和半径的取值区间通过地质勘查和工程类比相关资料确定。根据式(16)、(17)进行滑弧圆心坐标和半径的初始化：

(16)

(17)

图3 基于SA-FOA的边坡临界滑面搜索流程图

3 算例分析

基于该文提出的SA-FOA算法，结合瑞典条分法/简化Bishop法，在Matlab中编写了相应的边坡临界滑面搜索程序，并对两个算例进行了验证分析。

3.1 算例1：双层土体边坡

算例1是文献[8]中的非均质土坡，安全系数计算采用瑞典条分法。该文分别采用FOA和SA-FOA算法搜索该土坡的临界滑面，并与文献[8]的SA算法进行比较，参见表1及图4。3种算法中，FOA算法安全系数较大，SA-FOA算法安全系数最小，即基于模拟退火的果蝇优化算法寻优能力最强。

表1 不同算法的临界滑面搜索结果对比

图4 不同算法临界滑面形状和位置比较

在Matlab中，FOA、SA-FOA搜索最小安全系数的收敛过程如图5所示，对比可见，SA-FOA迭代曲线的斜率比FOA更大，即SA-FOA收敛速度更快；SA-FOA迭代200次后的安全系数稳定在0.857处，而FOA迭代200次后的安全系数还未达到0.93，因此SA-FOA的收敛稳定性更好，搜索能力较FOA显著提高。

图5 最小安全系数收敛进程

算例1证明将SA-FOA算法用于边坡临界滑面的搜索不仅是可行的，而且更高效、准确。

3.2 算例2：复杂的3层土体边坡

算例2取自澳大利亚计算机应用协会(ACADS)发布的考题，土层剖面如图6所示，土层参数见表2。

图6 ACADS考核题剖面图

该边坡土层厚度不均、材料参数相差较大，是一个复杂边坡，在最危险滑面搜索过程中会出现多个极值。SA-FOA算法结合简化Bishop 法，假设一个初始滑面，则局部危险滑面及最危险滑面的搜索结果列于表3。SA-FOA算法在初始滑面与最危险滑面(临界滑面)相差较大的情况下，能够跨越众多局部危险滑面，搜索到整体最危险滑面，证明了SA-FOA具有较强的全局搜索能力。

表2 各土层材料参数

表3 不同圆弧滑面的安全系数

针对该考题各国学者以及裁判推荐答案同时列于表4，各个答案均基于简化Bishop法。由表4可见，该文SA-FOA计算结果与评判员答案吻合度很高，且临界滑面的安全系数较其他答案略低，进一步验证了该文SA-FOA算法用于边坡临界滑面搜索中是行之有效的。

表4 该文计算结果与评判员答案比较

4 结论

(1) 该文所提出的SA-FOA算法与瑞典条分法结合，在简单土坡中能够快速、高效地搜索出临界滑面，与单纯采用SA算法或FOA算法相比，收敛速度更快，寻优结果更准确。

(2) SA-FOA算法与简化Bishop法结合，也适用于复杂边坡的临界滑面搜索，显示出较强的全局搜索能力，得到的最小安全系数和临界滑面更可靠。

(3) 将模拟退火算法与基本果蝇优化算法相结合，克服了单一算法的局限性，经验证具有较强的全局搜索能力和较快的收敛速度，为边坡临界滑面的搜索提供了一条高效的新途径。