基于GA-Sarma算法的边坡最不利滑裂面搜索方法

李生清

(山东省地质矿产勘查开发局八○一水文地质工程地质大队(山东省地矿工程勘察院)， 山东济南 250014)

0 引言

自然环境变迁以及人类社会过度开发，给土地、山体以及坡地带来较大的负面影响，不稳定的土壤条件、松散的地层以及活跃的大陆板块，使坡地出现了严重程度不同的滑裂面，影响自然生态安全。在过去的若干年间，国内外大量学者对于边坡稳定问题的研究一直很关注(祝志恒，2009；佘小年等，2011；樊成等，2020；邓宇轩等，2021；宋二祥，2021)。目前应用较多的方法为刚体极限平衡法，即通过简单计算将边坡稳定问题转化为二维问题。尽管一些方法能够快速获得滑裂面，但对于滑裂面的危险程度却难以区分。从20世纪中后期开始，越来越多的学者将边坡稳定性问题转化为三维问题，通过三维分析模型以及智能化、自动化系统搜索最不利滑裂面，但三维滑裂面本身较为复杂，因此该项工作仍存在较大困难(温树杰等，2018；韩林和蔡永昌，2020；石峰等，2021)。

针对几种常规搜索方法的不足，吴志轩等(2019)中提出利用Bishop方法和有限元方法搜索边坡最不利滑裂面，从无限长边坡基-填界面开挖台阶情况出发，利用Bishop 方法推导不同破坏模式下的安全系数公示，提出了理论解法并解释了稳定性强化机理。方宏伟等(2019)基于SLFT(滑移线场理论)设计边坡双折减系数强度失稳判据，根据极限状态下的边坡数据分析结果，讨论边坡局部位置是否开始出现裂纹。该方法利用整体安全系数定位发生问题的位置，实现对边坡稳定性、安全性的基本判断，并且相较于已有方法，提出了失稳判据的客观标准，排除了人为主观因素的影响。余鹏等(2022)基于XFEM构造模拟土质边坡失稳滑裂过程的方法，将边坡失稳过程概化为滑裂面的发生、扩展和贯通的过程。孔德志等(2022)基于极限平衡理论，建立目标函数与变量的关系表达式,通过求驻点的方法求解各条块的坐标值确定其相应的稳定性系数和滑裂面。卢坤林等(2009)基于滑面正应力修正模式的假定，根据边坡的平衡条件，导出含安全系数的平衡方程组,通过反复迭代求解安全系数。研究结果表明：该方法能够得到更小的安全系数，最危险滑面所对应的幂次随着滑体长度的增大先略微降低后逐渐增大。然而，上述传统方法在选取潜在滑裂面关键点时，未综合考量剪切滑裂点、受拉滑裂点对搜索过程的影响程度，导致搜索过程易陷入局部最优，从而丢失其他同等级搜索数据。

针对传统方法存在的不足，本研究针对GA-Sarma算法的遗传性和概率性特点，进一步优化边坡滑裂面搜索方法，通过GA-Sarma算法搜索边坡最不利滑裂面，为边坡支护、加固及改造等地质工作提供更加可靠的治理技术。

1 基于GA-Sarma算法的边坡最不利滑裂面搜索方法

1.1 计算稳定性系数提取潜在滑裂面

根据一般滑裂面的基本结构建立三维匀质边坡模型，以其中坡脚线中点为原点，设边坡水平面为x轴，与边坡水平面垂直的坡脚线设为y轴，过原点垂直于两轴的坐标轴为z轴。设置滑面方程，即存在：

z=f(x,y)

(1)

设置边坡土体粘聚力为a，内摩擦角为θ，土体重度为b，则可以做出三组假设：

(1)假设边坡发生变形时滑体为刚体(杨山奇等，2018；刘锋涛等，2019)；

(2)假设求解滑裂面存储势能时，变形后的微面变形大小可利用刚度系数dij模拟，则刚度系数与微面面积之间，存在正比例关系；

(3)假设外力作用产生虚位移能够减少滑裂面存储势能，则外力与虚位移分别为：

(2)

(3)

式中：λ表示随机微面单元外法线的单位向量，该值可通过公式(4)获得。由于每一个微面单元不一定是完全平整的曲面，为了便于计算与研究，默认微面单元为光滑平面，则根据微面坐标轴显示的坐标可得到该参数的取值：

(4)

(5)

确定土体单元受到的法向力计算剪应力S1(李珂等，2018)。由于该值的计算过程与上述弹性势能高度近似，因此不再重复说明，默认剪切势能为S2，则求解获得的总势能为：

(6)

(7)

式中：u、v、w表示函数求解参数。设置抗滑力参数为G1，设置边坡滑力参数为G2，则稳定性系数的一般计算表达式为：

G=G1/G2

(8)

G1和G2的值两组参数的计算表达式为：

(9)

式中：li表示第i处的活动距离；c表示滑动系数；α表示剪切力与虚位移之间夹角(孙加平等，2017)。公式(9)是在滑动系数的控制下，根据不同方向上的作用效果和边坡活动距离，对抗滑力和滑力参数展开计算。将公式(9)代入公式(8)中，求出边坡稳定性系数，然后可根据该系数提取潜在滑裂面。

1.2 选取潜在滑裂面关键点

随着大位移发生形成潜在滑裂面，无论在静力还是动力条件下，失稳时的坡体存在大面积破坏的情况(黄诗渊等，2018；冯晴枫等，2020)。已知岩土内部任意平面的剪应力超过极限剪应力，致使剪切应变增量发生突变现象，因此假设无限小的时间dt中，剪切应变增量为：

(10)

图1 搜索示意图

图1中的纵坐标表示坡高。将层数相同的单元合并为一组，利用单元ID标记每个单元。完成剪切滑点的选取工作后，选取受拉滑裂点。环境动态变化下，边坡受到水平惯性力，坡顶位置由于受拉而出现破裂(丁鑫品等，2020)。采用简化方法假设受拉裂缝为直向的细长裂缝，当坡顶出现连续受拉塑性单元时，叠加所有单元获得高度值，并将该值记为裂缝最大深度。受拉塑性单元与剪切塑性单元存在交互关系，针对这一研究结果定义受拉塑性单元高度为受拉裂缝深度，图2为受拉裂缝深度示意图。

图2 受拉裂缝深度示意图

图2中纵坐标表示裂缝深度，黑色区域表示拉塑性单元区，h表示受拉塑性单元与剪切塑性单元之间的深度(m)，也就是受拉裂缝深度，竖直提取该值以上单元坐标，作为受拉滑裂关键单元坐标。综合上述两组关键点，搜索边坡最不利滑裂面。

1.3 GA-Sarma算法搜索边坡滑裂面

GA-Sarma算法能够解决以折线形为滑面形态、滑裂路径可追踪、坡向不连续的岩质边坡最危险滑面的全局最优问题(张剑波等，2013；王紫杰等，2020)，其将滑裂面看作折线形，利用Sarma极限平衡法确定稳定性系数，通过遗传算法优化搜索最不利滑裂面的位置。很明显，GA-Sarma算法是具有搜索全局最优问题能力的遗传算法和传统力学模型的融合，为快速、合理得搜索边坡滑裂面提供了技术支持(毛谦等，2008；林永生和陈胜宏，2013；麻官亮等，2016)。根据预先模拟的边坡模型，该算法首先建立目标函数。当滑裂面由随机一点向坡顶扩展时，可能存在无数条变化路径，因此假设滑裂路径的方向为β。当坡体内存在顺坡向不连续结构面时，则路径沿不连续结构面变化(王章琼等，2018)。因此根据Sarma算法得到：

K=g(φn,γn)

(11)

式中：φn、γn表示方向为β时，选取的剪切滑裂点和受拉滑裂点。决策变量就是上述公式中参数的数量，因此将函数优化问题转换成算法搜索问题。根据其他方法中适应度函数的使用方法可知，该方法可以度量个体最优解。将求解稳定性系数作为滑裂面全局搜索的一个最小值问题，则得到适应度函数的计算公式：

h(φn,γn)=max(g(φn,γn))-g(φn,γn)

(12)

公式(12)的物理意义，描述了稳定性系数最小的路径适应度最大，当搜索过程中出现变异问题时，有更大机率被保存下来。算法利用选择算子评价个体，根据得到的适应度评价结果，强化算法的收敛性能(窦勇，2019)。已知随机采样方式下，适应度实际值与个体选中概率呈正相关，因此假设N表示整个搜索群体集合，用i表示个体排序情况，通过上述公式(12)获得，则概率结果为：

(13)

式(13)中：i∈N表示单个个体和个体总数量。交叉算子作为有效计算工具，能够产生新个体，提高整个群体中所有个体的多样性和动态变化特性，防止部分个体结果，提前成熟而影响最终的位置搜索(张子映等，2018)。为了防止破坏大量现有的较好模式，要求变异概率Pi的值要满足最小值要求，则变异操作可按照图3所示的方法执行。

图3 变异操作示意图

为摆脱个体多样性丧失、搜索进程陷入僵局的问题，GA-Sarma算法通过多次增大变异概率Pi的值改善这些现象，同时为了进一步加强搜索结果，还引入了灾变策略(CS)，大规模消减处理群体，并进行新后代生成操作，强化搜索最优这一目的。至此基于GA-Sarma算法，实现边坡最不利滑裂面搜索方法。

2 实验与分析

为验证上述设计的基于GA-Sarma算法的边坡最不利滑裂面搜索方法的应用性能，设计如下实验。

2.1 实验准备

通过地形测量技术、地质测量技术、物探技术、钻探技术以及岩矿分析等技术手段，对西北某地区中存在潜在滑裂面的峡谷、山坡进行勘察。此次勘察工作依据8项不同的城市治理规范、土地使用规范、滑坡防治规范以及工程勘察方案而进行。

实验的工程测量工作中，采用2000国家大地坐标系，建立仿真环境测试模型，高模型的基准为1985国家高程基准，设置比例尺为1:500，范围为滑坡体及周围外的50 m范围内，测量面积0.2 km2。利用无人机对整个边坡进行扫描，单次扫描面积设置为0.30 km2。采用高密度电法进行区域物探。在选择的勘测区域内设置260个物理点、15个质检点。计算各点均方相对误差，计算公式为：

(14)

图4 Ⅱ区、Ⅲ区山体滑动裂缝图(部分)

根据图4显示的实际调查结果可知，在实验选择的勘测区域中，存在大量威胁性较大的滑裂面以及潜在滑裂面。为此，使用不同的勘察手段分析测试区域内的岩土类型。由于篇幅的限制，本研究主要对Ⅲ区展开研究。表1为获得的Ⅲ区各个岩土层力学指标。

表1 测试区域岩土层主要力学指标勘测结果

为了检验本文方法的实际应用性能，设计对比应用实验。将本文方法作为实验组，将两组传统搜索方法-Bishop方法及基于SLFT(滑移线场理论)设计边坡双折减系数强度的失稳判据，然后根据极限状态下的边坡数据分析结果，讨论边坡局部位置是否开始出现裂纹方法(分别是方宏伟等(2019)和呈志轩等(2019)所得的方法)分别作为对照A组、对照B组，结合图4与表1显示的数据，对整个勘察地区进行边坡最不利滑面搜索工作，通过比较搜索结果是否最优，来确定不同方法之间的差异性。

2.2 最优解测试结果分析

利用仿真软件、CAD等软件模拟整个搜索区域的三维影像，已知仿真场景中存在3处最不利滑裂面、5处刚刚形成且危险程度偏低的滑裂面以及14处潜在滑裂面。因此，分别以最不利滑裂面、刚刚形成且危险程度偏低的滑裂面(一般滑裂面)和潜在滑裂面为搜索对象，利用不同方法对其展开搜索，根据搜索结果验证不同方法的实际性能，其比较准则为不同方法对仿真场景中的不同滑裂面搜索的吻合度和数量。图5为三组不同方法的最优滑裂面搜索结果。

图5 最优搜索结果对比

根据图5显示的测试结果可知，面对同样的搜索场景，实验组获得了所有最不利滑裂面的所在位置。对照组A遗漏了1处最不利滑裂面所在位置，同时将3处刚刚形成且危险程度偏低的滑裂面，当成了最不利滑裂面。对照B组只获得了一组最不利滑裂面所在位置，剩余的2处最不利滑裂面直至实验测试结束，也没有被发现，同时将4处刚刚形成且危险程度偏低的滑裂面，当成了最不利滑裂面。综合上述实验测试结果可以看出，两组常规边坡最不利滑裂面搜索方法中，使用的算法性能较差，过早收敛导致获得的结果不是最优。

在此基础上，为了进一步验证不同方法的应用性能、保证实验测试结果的可靠性，保证基本测试条件不变，利用3种方法对选择的测试场地中的最不利滑裂面展开20轮搜索，以此来更直观地观察不同方法的搜索效果。表2为20轮测试过程中得到的最优解测试结果。

表2 最不利滑裂面搜索最优解统计表

表2汇总结束后，发现对照B组从第11轮开始，一直到第17轮的测试结果均为0。为了保证实验测试数据完全可信，利用事先安装的监控软件，查询对照B组仪器设备以及软件的运行状态，发现对照B组正常完成了20轮搜索工作，并不存在任务系统或软硬件故障问题，因此第11轮至第17轮的测试结果可信。根据表2中的数据可以直观看到，实验组方法在20轮测试中，全部得到了3组最不利滑裂面；对照A组每次搜索最不利滑裂面时，总是存在遗漏位置的情况；对照B组则有35%的概率，难以获得最不利滑裂面的搜索结果。

为了方便比较，计算三组方法20轮测试的最优解平均值，分别为3处、1.4处以及0.8处。根据该计算结果可以得出：在同样的测试背景环境下，本文方法完全能够获得最优解，而两组传统方法获得最优解的平均值不足2处和1处，即面对3处最不利滑裂面所在位置，两组方法每次搜索时只能获得1.4处以及0.8处最不利滑裂面。由此可见，本文设计的基于GA-Sarma算法的边坡最不利滑裂面搜索方法更加满足设计预期。

3 结论

针对实际实验结果，本研究在选取潜在滑裂面关键点的基础上，引入了GA-Sarma算法，通过加强算法自身的计算能力、降低其他因子对算法的影响，并建立目标函数和适应度函数，在考虑异变的前提下搜索边坡最不利滑裂面。具体结论如下：

(1)工程测量工作中，采用2000国家大地坐标系，建立仿真环境测试模型，利用无人机对整个边坡进行扫描且采用高密度电法进行区域物探。勘测区域质量检查总均方相对误差在±3.42%之间，满足上述8项规范中不大于±5%的要求。

(2)本文引入GA-Sarma算法中，利用遗传算法对边坡最不利滑裂面进行搜索，显示出强大的全局最优搜索能力。遗传算法和Sarma算法结合是具有搜索全局最优问题能力的遗传算法和传统力学模型的融合，为快速、可靠得搜索边坡滑裂面提供了技术支持。

(3)在20轮测试中，对照A组每次搜索最不利滑裂面时，总是存在遗漏位置的情况；对照B组则有35%的概率难以搜索到最不利滑裂面；实验组方法全部得到了3个最不利滑裂面，取得了最好的搜索效果。

此外，在今后研究工作中，基于GA-Sarma算法进一步扩大搜索范围，对多种不同地质进行系统研究，优化特征提取或智能化处理步骤是值得研究的课题，为土地保护与治理提供更加快速、可靠的技术支持。