基于混合变异策略差分进化算法的边坡滑裂面搜索研究

张子映， 柴军瑞， 张书滨， 钱武文

(1.西安理工大学 西北旱区生态水利工程国家重点实验室培育基地， 陕西 西安 710048；2.江西省河道湖泊管理局， 江西 南昌 330009)

1 研究背景

边坡滑动面搜索是边坡稳定计算中的一个关键问题，其实质为安全系数[1-4]与最危险滑动面[5-7]的搜索问题。在这一领域出现了较多的优化搜索方法，模式搜索法、动态规划法、数学归纳法、变分法、随机搜索法等。随着优化方法研究的深入，人工智能搜索算法如遗传算法[8-9]、模拟退火算法、粒子群优化算法、蚁群算法[10]、神经网络以及蒙特卡洛法[11]等被用于边坡最危险滑动面搜索，且目前已成为滑动面搜索研究中的主流方法。这些智能寻优法各有其优点但也都存在不足之处，例如：遗传算法存在早熟收敛、算法性能对参数选取强依赖性而参数不易确定等缺陷[9]；模拟退火算法在实际应用中由于退火不充分并不能保证一定能搜索到全局最优值；粒子群优化算法[12-13]易陷入局部最优；蚁群算法[14-15]最初被设计用于组合优化问题，其缺少用于连续数值优化的操作算子，用于边坡滑动面搜索问题时，需要对边坡进行一定精度的离散,操作过程复杂，同时也存在易于陷入局部最优的问题，有待进一步改进。

差分进化(differential evolution，简称 DE)算法是由 Storn等[16]于1995年提出的一种智能优化算法，该方法具有操作简单、收敛性好、全局寻优能力强的优点[17-19]。钱武文等[20]受Nelder-Mead方法启发，提出了一种名为“反射变异”的新型变异策略，具有更明确的变异方向，可即使基向量拥有一定的随机性，但仍为一种贪婪搜索策略，为进一步增强其全局收敛性能，本文在此基础上引入一个基本变异策略形成混合变异策略，并将该优化方法用于边坡临界滑动面的搜索中，通过两个经典考核算例，验证了该方法应用到边坡稳定性研究是合理可行的。

2 改进的差分进化算法

算法的全局搜索能力和局部开采能力是相互对立的，变异策略高度影响着差分进化算法的性能。参考文献[20]提出了一种名为“反射变异”的新型变异策略，并与参数自适应调整方法结合，形成了一种具有良好全局探测能力和局部开采能力的差分进化算法(RMADE)。考虑到“反射变异”从某种程度上仍属于贪婪变异策略，仍易陷入局部最优，因此为进一步改进其全局收敛能力，本文在RMADE的基础上引入一种基本变异策略组成混合变异策略共同控制子代个体的生成。

2.1 混合变异策略

受Nelder-Mead方法启发，文献[20]提出了一种名为反射变异策略的新型变异策略，该策略不论求解问题的维数，始终选择种群中的4个任意个体组成单纯形，每完成一次变异操作须引入一个单纯形。该策略操作如下:

Vi,g=Xo,g+F(Xa,g-Xd,g)

(1)

式中：Vi,g为i个体在g代的试向量；Xo,g为单纯形中心；Xa,g和Xd,g分别为单纯形中的最好点和最坏点，F为变异缩放因子。

不同于Nelder-Mead方法，本文使用的单纯形中心Xo,g是根据单纯形顶点的函数值计算，并且除最坏点外的各个顶点均按其函数值权重分配到中心中。单纯形中心Xo,g的计算公式如下：

(2)

式中：w1，w2和w3分别为除最坏点外的其余单纯形顶点的权重；Xa,g，Xb,g和Xc,g为除最坏点外的剩余单纯形顶点。

这种变异策略类似于Nelder-Mead方法中的反射操作，具有明确的搜索方向。即使基向量拥有一定的随机性，但仍为一种贪婪搜索策略，为进一步增强其全局收敛性能，本文还引入一种基本子代生成策略(DE/current-to-rand/1)，其变异向量生成如下。

Vi,g=Xi,g+F(Xc,g-Xi,g+Xa,g-Xb,g)

(3)

式中：Xi,g为当前操作个体，其他同上。

引入DE/current-to-rand/1后，一个种群中将同时存在两种变异策略，本文根据各策略的成功率(成功更新子代个体的比率)为各个个体轮盘赌选择其对应的变异策略。成功率的计算公式如下：

(4)

式中：Nk为种群中使用策略k的个体总数；srk为策略k的成功率，在初始种群时将其设置为0.5，以保证算法开始时各变异策略拥有总体相等的个体数。

2.2 控制参数自适应调整机制

传统的差分进化算法采用固定的变异因子和交叉因子，而针对不同的优化问题，使用固定的控制参数将严重妨碍算法的性能。

参数自适应调整机制允许每一个体使用不同于其他个体的缩放因子F和交叉因子CR去生成子代个体，且能根据算法进化过程中成功个体的经验反馈调节控制参数。本文采用SHADE[19]算法中的参数自适应调整方法。

2.3 算法的基本思想和步骤

本文在文献[20]中算法的基础上，引入基本的变异策略DE/current-to-rand/1，从而进一步改进算法的全局搜索能力，并在降低早熟收敛可能性的同时保持原算法的收敛速度。使用自适应参数调节方法进一步提升了算法的整体性能。算法的具体步骤如下：

Step 1：初始化参数和初始种群；

Step 2：判定是否达到终止条件,若是转Step 7，若否则转Step 3；

Step 3：为每个个体分配变异因子和交叉因子；

Step 4：根据各变异策略的成功率使用轮盘赌方法为各个体分配变异策略；

Step 5：对每一靶个体，按公式(1)或(3)完成变异操作，然后完成交叉操作；

Step 6：计算试向量的函数值并与靶个体比较，选择较优个体替换靶个体；

Step 7：输出最优个体。

3 基于DE的边坡稳定分析

3.1 滑动面上任意点应力张量计算

为克服极限平衡法分析结果不考虑土体应力应变状态的缺陷，本文采用有限元软件进行边坡极限平衡分析，即先计算获得边坡高斯应力场，然后根据双线性插值计算得到各单元节点的应力值。

(5)

式中：Ni为双线性插值函数；ξ，η为自然坐标，本文采用笛卡尔坐标系。

生成滑动面y=y(x)，对滑动面离散成多个线单元。依次对滑动面上每个线单元的每个高斯点查找其归属的有限元单元，由这个单元的各个节点的方向应力形函数插值得到高斯点的方向应力，即σx，σy，τxy。

(6)

3.2 安全系数计算

(8)

式中：c、φ分别为材料的凝聚力和内摩擦角。σn、τ分别为曲线上任意一点土体方向应力和这点位置沿曲线的剪应力。

若将滑动面视为由一系列连续线段组成，将滑裂面离散化处理，应用高斯积分得边坡安全系数为:

(9)

3.3 任意滑裂面构造方法

滑裂面由n个控制点{(x1,y1) , (x2,y2),…, (xn,yn) }组成的n-1个小片段近似构造。由于滑入点和滑出点都位于边坡轮廓线上，其纵坐标y1和yn都可据已知的轮廓线求出，因此通过事先固定中间控制点的横坐标，进一步减少搜索变量的个数，即由n个控制点近似的滑裂面的搜索变量为S={x1,y2,…,yn-1,xn}。

合理的滑裂面并非这些点在平面内任意分布的组合，通常将这些片段组成的滑裂面向上凹时视为合理，如图1所示，这意味着滑裂面通常需要满足以下特征[21]，即：

α1≥α2≥…≥αn-1

(10)

式中：αi为线段i的倾角；n为滑裂面控制点的个数。

式(10)亦可用线段斜率ki表示，即：

k1≥k2≥…≥kn-1

(11)

使用点坐标表示后，有：

(12)

则优化问题变为带有n-2个约束的优化问题，优化变量为S={x1,y2,…,yn-1,xn}，即：

minF(S),S={x1,y2,…，yn-1,xn}

(13)

j=1,2,…,n-2;li≤Si≤ui；

i=1,2,…,n

在智能优化领域，约束的存在易造成可行域的减少进而使得解的搜索变得更加复杂。边坡临界滑裂面的搜索属于有约束优化问题，考虑到一般情况不合理的滑裂面对应的安全系数比合理滑裂面更大的特性，可将有约束优化直接按无约束优化处理。

图1 合理的边坡滑动面

3.4 求解最危险滑裂面

求解最危险滑裂面的程序流程如下(程序框图见图2)：

(1)根据边坡地层分布情况及边坡剖面界限，运用有限元软件进行合理的有限元网格划分，并计算边坡应力场。

(2)差分进化算法应用于滑裂面搜索，每个个体代表一个滑裂面，每个滑裂面由n个控制点组成的曲线，故每个个体就是一个n维的向量Xi=[x1,y2,…,yn-1,xn]。

(3)在算法搜索前，根据剖面的几何模型对个体的每个维度设定合理的搜索区域。其控制方式为：第一个滑入点和最后一个滑出点的搜索范围位于坡表，中间控制点纵坐标为搜索区间随机初始化值y2,y3,…yn-1。

(4)构建初始种群{X1,X2,…,XNP}后，使用有限元极限平衡法计算个体Xi的安全系数Ki，用差分进化算法搜索出最佳滑裂面Xbest。

图2 基于差分进化算法的有限元极限平衡法的流程图

4 算例验证

本算例使用8个滑裂面控制点(滑入点、滑出点以及6个中间控制点)，因此只需搜索滑入点和滑出点的X值以及3个中间控制点的Y坐标值。为了使中间控制点更好地反映滑面的形状，尽可能地将中间控制点的横坐标均匀分布在搜索域内。种群规模设置为40，最大迭代次数设为150。

4.1 均质边坡

引用澳大利亚计算机应用协会(ACADS)通用考题1，坡高为10 m，坡角30°，土容重γ=20 kN/m3，黏聚力c=3 kPa，土的弹性模量E=104kN/m2，内摩擦角φ=19.6°，泊松比ν=0.25。推荐的安全系数为1.0。应力计算模型如图3所示，假定区域左右边界为法向约束，底部边界为固定端约束。

建立有限元模型，使用理想弹塑性本构模型计算基于Mohr-Coulomb强度破坏准则的边坡土体的应力场分布，进而以本文方法确定最危险滑动面的位置和其对应的安全系数。为了进行比较，将结果与使用传统极限平衡法spencer、Bishop和Janbu法稳定性分析的结果比较，将各方法所得到的安全系数和最危险滑裂面汇总绘制于图4中。由图4可知，圆弧滑动搜索结果为KJanbu=0.966，KBishop=0.996，KSpencer=0.995，推荐安全系数为1.0。本文方法搜索的边坡最危险滑面安全系数为1.055，仅与推荐的安全系数相差5.5%，在合理范围，本文方法正确。

图3 边坡几何模型

图4 各方法搜索的最危险滑裂面对比

4.2 非均质边坡

引用澳大利亚计算机应用协会(ACADS)一通用考题。坡体由3种土料构成，坡高为 10 m，坡比为1∶2，边坡几何特征如图5所示，其材料特性如表1所示，推荐的安全系数为1.0。

图5 边坡几何模型

表1 非均质边坡材料性质

建立有限元模型，使用理想弹塑性本构模型计算基于Mohr-Coulomb强度破坏准则的边坡土体的应力场分布，进而以本文方法确定最危险滑动面的位置和其对应的安全系数。同上，将各方法所得到的安全系数和最危险滑裂面汇总绘制于图6中。

由图6可知，圆弧滑动搜索结果为KJanbu=1.341，KBishop=1.442，KSpencer=1.419，推荐安全系数为1.39。本文方法搜索的边坡最危险滑面安全系数为1.375，比推荐的安全系数小1.08%，在合理范围。本文计算结果处于Janbu法与Bishop法和Spencer法之间，且略小于Bishop法和Spencer法，说明本文方法的准确性。基于有限元弹塑性变形的极限平衡法滑面起始部分较圆弧法略有“上翘”，尾部略有“下沉”，考虑到土体的弹塑性变形这种现象是可以解释的。说明是否考虑塑性变形对边坡的滑面位置和安全系数影响较大。

图6 各方法搜索出的最危险滑裂面对比

4.3 收敛性能比较

使用上节非均质边坡作为计算模型，将本文改进的差分进化算法与单纯形法、DE、RMADE进行比较，研究不同算法对边坡临界滑动面搜索的影响。鉴于单纯形法为无约束优化算法，而滑面受到边坡坡体区域的限制，因此本文对单纯形法中的初始单纯形生成方法进行了略微调整，即引入差分进化算法中的初始种群技术生成初始单纯形。

文中使用的单纯形法的反射系数、扩张系数以及收缩系数分别为1.0、1.3和0.7，DE的缩放因子F=0.8，交叉因子CR=0.8。种群规模NP=40，最大迭代次数为200。为了研究算法搜索临界滑面的稳定性，各算法均独立运行10次，记录最大值、平均值、最小值以及标准差于表2中，各算法的典型收敛过程线绘于图7中。

表2 独立10次计算下各算法性能比较

图7 各算法的收敛过程

由表2可知，单纯形法结果最为发散，说明其已陷入局部最优，本文算法的平均值和标准差最小，说明其性能较其他方法要好且更稳定。由图7可知，本文算法拥有较快的收敛速度。综合可知，本文算法性能优于其他几种方法，说明本文对RMADE的改进具有实际意义。

5 结 论

(1)本文在RMADE的基础上引入一个基本的变异策略DE/current-to-rand/1而形成了一种改进的差分进化算法，该方法根据各变异策略的子代成功率使用轮盘赌选择自动选择个体变异策略，并与其他优化方法进行比较验证了本文提出的算法具有更好的寻优性能。

(2)将本文改进的差分进化算法与有限元极限平衡法结合，应用于边坡稳定性分析中，通过两个算例验证了该法解决边坡稳定性问题的有效性。