基于最大最小蚁群优化算法的缓层土质三维边坡临界滑动面搜索

张江辉，张 勤，吕加贺*，邹俊鹏，焦玉勇

(1.中国地质大学(武汉)工程学院，湖北 武汉 430074；2.中国市政工程中南设计研究总院有限公司，湖北 武汉 430010)

滑坡是现实生活中最为常见的一种地质灾害，对水利、道路、矿山等工程建设以及人民的生命和财产安全都有着巨大的威胁。解决滑坡问题的关键在于通过对边坡稳定性设定合理的评判标准，从而准确地预测滑坡的产生。因此，对边坡稳定性的分析是解决滑坡问题的关键，而边坡稳定性分析的要点在于边坡临界滑动面的搜索以及边坡安全系数的计算。

现实中发生的滑坡现象，是一个典型的三维问题。然而，目前大多数边坡稳定性分析主要是采用二维分析方法，通过对简化的二维模型进行分析，搜索边坡的临界滑动面并计算边坡的安全系数，但这样计算出的边坡安全系数是偏于保守的。而三维边坡稳定性分析能够考虑边坡在空间内的发育特征，计算出的边坡临界滑动面以及边坡安全系数更加贴合实际，更具有指导意义，故本文考虑构建一种三维边坡临界滑动面智能搜索方法用于三维边坡的稳定性分析。

常用的边坡临界滑动面搜索方法可分为变分法、固定模式搜索法、数学规划法、动态规划法、随机搜索法和智能搜索法6类[1]。其中，智能搜索法以其显著的优势[2-3]，成为目前边坡临界滑动面搜索方法中的主流算法，而蚁群算法[4]以其特有的组合优化方式，在智能搜索法中具有重要的地位。蚁群算法最初是用于解决优化领域的问题，如旅行商问题、二次分配问题，其具有优化速度快、全局搜索能力强、易得出最优值等优点，且适用于非线性、离散型的优化问题，这些特点使其对于求解边坡临界滑动面搜索问题也极为适用。一些学者[5-8]已将蚁群算法应用于二维边坡临界滑动面搜索问题研究中，并取得了良好的效果。本文在研究了蚁群算法的机理后，发现其同样适用于三维边坡，故将蚁群算法扩展到三维边坡临界滑动面搜索问题研究中。然而，传统的蚁群算法也存在解的质量较低、易陷入局部极值等缺点，故本文采用近年来提出的一种求解效果较好且实现简便的改进算法——最大最小蚁群优化算法[9]来搜索三维边坡的临界滑动面。最大最小蚁群优化算法相比于传统蚁群算法的主要优势在于：①能够更加有效地利用搜索过程中找到的最优解，并且将蚂蚁搜索引导到高质量的解；②避免蚂蚁搜索的过早收敛。

最大最小蚁群优化算法搜索边坡临界滑动面的过程主要是基于边坡的应力场进行的，对边坡安全系数的求解采用有限元法，搜索过程中考虑了边坡岩土体的非线性本构关系、复杂的几何结构和材料性质。使用有限元法[10-12]求解边坡安全系数的主要方法为强度折减法，但是该方法需要多次验算，由于弹塑性模型在计算到接近极限状态时收敛缓慢，所以其在三维问题中计算效率较低，并且通过强度折减法反复试算所求得的边坡安全系数可能要比真实的边坡安全系数稍高[13]。故本文考虑采用有限元法中的滑动面应力法[14]，该方法只需要对边坡进行一次弹塑性计算，在计算效率上较为出色，且避免了通过强度折减法反复试算所带来的误差，得到的边坡安全系数计算结果更接近真实值。

鉴于采用传统有限元法在处理复杂边坡时，往往会由于网格畸变而造成应力失准，为了求得准确的边坡应力场，本文考虑采用一种改进的有限元法——光滑有限元法[15]。光滑有限元法相较于传统有限元法，有以下优势[16]：①采用线性点插值法来表示形函数，避免了坐标变换，从而降低了网格畸变对应力精度的影响；②软化了刚度矩阵，可有效避免许多问题中产生的体积锁定；③该方法中的相容应变场可改善传统有限元法中的应变场不连续问题，从而明显提高了位移解和应力解的精度；④该方法在低阶单元中可以得到较高的计算精度。

1 光滑有限元法

光滑有限元法的本质是将传统有限元法与无网格法中的光滑应变技术结合，改进有限元结构刚度的一种方法。Liu等[15]在建立了光滑有限元法的理论基础后，又提出了一系列光滑有限元数值算法，主要包括光滑子单元域有限元法、结点光滑有限元法、边光滑有限元法和面光滑有限元法等。其中，处理三维问题主要采用面光滑有限元法和边光滑有限元法。边光滑有限元法具有空间离散稳定性和时间响应稳定性，相比于传统有限元法模型中偏硬的刚度，边光滑有限元法模型中的刚度更接近结构的实际刚度，可以对动态问题给出精确、稳定的结果，对静态问题给出超精确的结果，并且其在静态和动态分析中得到的结果精度要优于面光滑有限元法[17]。显然，对于三维边坡问题，边光滑有限元法是适用的。

1. 1 光滑有限元法的原理

与传统有限元法相同，光滑有限元法也是通过构造形函数来建立问题域的位移场，其表达式如下：

(1)

对于光滑有限元法中常用的三角形单元和四面体单元，可以直接使用传统有限元法中的三角形单元和四面体单元的线性形函数，在此不作赘述。

利用应变/梯度光滑技术创建光滑应变场，当相容应变场易得到时，光滑应变场可通过修正相容应变场获得：

(2)

(3)

光滑应变场也可通过对位移场求导获得：

(4)

式中：Ld为微分算子矩阵。

(5)

(6)

(7)

式中：nx、ny分别为Ln(x)在x轴和y轴的外法向分量。

光滑域的划分形式以三维的边光滑有限元法为例，它以单元边为基准划分光滑域，其光滑域由基准边上的两个结点、相邻单元面的形心点和相邻单元的形心点连接而成，见图1。假设单元边总数为Neg，光滑域总数为Nx，单元边与光滑域为一一对应关系，即：Ns=Neg。

图1 边光滑有限元法中的光滑域划分

将公式(1)代入公式(6)，可得：

(8)

(9)

(10)

(11)

光滑有限元法的总体平衡方程可通过使用光滑伽辽金弱化形式得到：

(12)

将公式(1)和公式(8)代入公式(12)并化简，可得标准的离散方程组：

(13)

(14)

1. 2 光滑有限元法的精度分析

本文通过一个线弹性模型对边光滑有限元法相较于传统有限元法和面光滑有限元法的优势进行分析。如图2所示是一个上侧受均布压强、左侧固支的三维悬臂梁模型，将对该模型选取27 217个十结点四面体单元网格计算所得的应变能7.081×104N·m作为参考解。模型采用的材料参数分别为杨氏模量E7.2×104Pa、泊松比υ0.33、均布压强P10 N/m2。

图2 三维悬臂梁模型示意图

由于光滑有限元法主要以低阶单元为研究对象，所以本文采用线性四面体单元来划分模型。在本节中以不同的网格尺寸(以自由度来表示)对三维悬梁臂模型进行划分，模型网格划分形式见图3。图3中三维悬梁臂模型的自由度依次为552、1 188、1 395、2 901、3 384。

图3 三维悬梁臂模型网格划分形式

本文分别采用边光滑有限元法、面光滑有限元法和传统有限元法对三维悬梁壁模型进行分析，得到3组模型应变能，并结合参考解，对3种方法下模型应变能随自由度的变化趋势进行了对比分析，如图4所示。

图4 3种方法下三维悬梁臂模型应变能随自由度 变化趋势的对比

由图4可见：随着自由度的增加，3种方法所求得的模型应变能整体变化趋势都是向着参考解收敛，说明3种方法都具有良好的收敛性，且边光滑有限元法在较稀疏的网格状态下，就能取得较高的计算精度，说明在同等自由度情况下其解的精度显然要高于面光滑有限元法和传统有限元法。

2 最大最小蚁群优化算法

在自然界中，蚁群具备找到食物与巢穴之间最短路径的能力，这种能力不是单只蚂蚁拥有的，而是依靠蚁群之间的联系实现的。蚂蚁在其经过的路径上会留下一种挥发性分泌物(信息素)，在选择前进路径时，会根据先行的蚂蚁留下的信息素来选择要走的路径，且选择的概率与该路径上信息素的浓度有关，两者之间是一种正相关的关系，即某一条路径上经过的蚂蚁越多，信息素的浓度便越高，后面的蚂蚁选择该路径的可能性就越大，进而形成一种正反馈机制。正是依靠这种机制，在一定数目的蚂蚁经过之后，蚁群便可以找到一条最短的路径。Dorigo等[4]根据蚂蚁的这一行为提出了一种模拟自然界蚂蚁行为的智能优化算法——蚁群算法。

2. 1 滑动面应力法实现思路

实现蚁群算法之前，首先应明确边坡滑动面安全系数的求解方法，本文采用边坡稳定性有限元法中的滑动面应力法。滑动面应力法是将有限元法分析中产生的应力场与传统的极限平衡法相结合，来求取边坡滑动面安全系数的一种方法。下面以本文算法的实现流程为例，对滑动面应力法的实现思路做简要的介绍。

由于本文在搜索边坡临界滑动面时，是以二维线条来拟合三维的面，所以边坡滑动面安全系数的求解方法本质上与二维是一致的，故在此以xy平面中边坡滑动面安全系数的求解为例来介绍滑动面应力法。在求得边坡的应力场和滑动面之后，如图5所示将滑动面的一个剖面划分为n个小段Δli(小段取蚂蚁行走路径上相邻两个条带结点之间的路径)，每一小段上的应力取该段中点处的应力，表示为σxi、σyi、τxyi，则作用在小段Δli上的法向应力σni和剪应力τni分别为

图5 边坡剖面上的滑动面参数示意图

τxyisin2θi

(15)

(16)

在本文中由于对边坡基于光滑有限元法进行弹塑性分析时采用的是莫尔-库伦屈服准则，所以在小段Δli上岩土体的抗剪强度τji应为

τji=ci-σnitanφi

(17)

式中：ci、φi分别为小段Δli上岩土体的黏聚力和内摩擦角。

按照上述方法，将边坡滑动面所划分的所有小段上岩土体的剪应力和抗剪强度都求出之后，可得边坡滑动面安全系数FS为

(18)

对于公式(18)，在程序实现时，采用逐步累加的方式进行求解，在累加中若存在边坡滑动面安全系数的值为负值，则将其取为0，再进行后续的累加，直到得到最终结果。

2. 2 最大最小蚁群优化算法实现思路

在使用最大最小蚁群优化算法搜索边坡临界滑动面时，首先应对边坡进行离散化处理(见图6)，即假定将边坡划分为m×n条条带，每条条带上有k个离散点，从而形成了一组条带结点网络，离散点的位置便可由m、n、k三个参数值确定；然后根据边坡本身的性质，设置蚁群在边坡上行走的进出口，让蚁群在这些离散点上行走，寻找最优路径。

图6 边坡临界滑动面搜索方式示意图

蚂蚁在三维条带中的一次行走，即蚁群算法的一次迭代主要分两步：①在边坡滑入口区域随机确定一个点，作为出发点，此时保持n的坐标不变，让蚂蚁在mk平面中寻找出一条行走路径(见图6中紫色线条)；②在蚂蚁在mk平面中寻找路径时，每当到达一个新的点，认为蚂蚁有3种可能的行走方向，即首先在mk平面中行走，其次保持m坐标不变，在nk平面中向两侧依次生成蚂蚁的行走路径(见图6中红色线条)。在蚂蚁行走结束后，便可将所有搜索出的路径拟合成一个滑动面。

(19)

式中：Mr+1为第r+1条条分带上的离散点数；τ[(r,i),(r+1,j)](t)]为处于时刻t时，路径[(r,i),(r+1,j)]上的信息素浓度；η[(r,i),(r+1,j)](t)为处于时刻t时，路径[(r,i),(r+1,j)]上的信息素能见度；υ和β分别为信息素浓度和信息素能见度在蚂蚁选择路径时的相对重要程度(υ≥0,β≥0)。

图7 蚁群行走方式

在初始时刻，设定各种路径上的信息素浓度相等。与自然界中蚂蚁的信息素浓度变化趋势相同，随着迭代的进行，边坡中各路径上的原有信息素浓度会逐渐衰减，给定衰减系数ρ，则1-ρ代表信息素浓度的衰减度。如果有n只蚂蚁在一次迭代中找到了一个边坡的滑动面，则会在其经过的路径上释放新的信息素，故迭代后边坡中各种路径的信息素浓度可调整为

τ[(r,i),(r+1,j)](t+1)

(20)

(21)

若蚂蚁k没有经过路径[(r,i),(r+1,j)]，则：

(22)

在边坡稳定性问题研究中，对于信息素的能见度，可将其定义为一个与路径间距有关的量，对于路径[(r,i),(r+1,j)]，其信息素的能见度为

η[(r,i),(r+1,j)]=1/dij

(23)

式中：dij代表第r条条分带上的第i个离散点(r,i)与第r+1条条分带上的第j个离散点(r+1,j)之间的距离。

由公式(20)可知，在蚂蚁搜索出的滑动面所包含的路径上，信息素浓度将随着迭代的进行而增大，而其他路径上的信息素浓度会有一定的衰减。因此，随着迭代的进行，在边坡安全系数最小的滑动面所包含的路径上，信息素浓度将随着迭代次数的增加而不断增大，而其他路径上的信息素浓度将一直衰减。

对于一般的蚁群算法，由于信息素的更新可能会造成某些路径被过分强调，使得在几次迭代后，搜索便会过早地收敛，这样的收敛具有一定的偶然性。最大最小蚁群算法针对这一现象，给路径上的信息素浓度划定了一个范围[τmin,τmax]，超出这个范围的信息素浓度被设定为τmin或τmax。

此外，在蚁群算法的实现中，需要满足以下几个基本的条件：①滑动面的进出口位于边坡的表面；②滑动面的水平、垂直坐标范围应在边坡的水平、垂直坐标范围内；③滑动面应表现为下凸形，即满足θi+1≤θi；④滑动面与水平面的夹角变化应较为平缓，即θi-θi+1≤Δδ，其中Δδ为一个常量。由最大最小蚁群优化算法的搜索流程可以看出，蚁群算法中各参数的选择对其影响极大，主要表现在对收敛速度和求解质量的影响上，但目前针对这些参数的选取没有系统的理论依托，大多是依靠经验来选取。假定边坡滑入区与滑出区之间共有n条条带，经过程序的试算，本文算法中的一些基本参数设定如下：υ=2，β=1，τ0=0.5，τmin=0.1，τmax=5.0，θ0=π/2，Q=1.5，ρ=0.5，Δδ=π/8,蚂蚁总数NA为n。利用最大最小蚁群优化算法搜索边坡临界滑动面的具体流程，如图8所示。

图8 基于最大最小蚁群优化算法搜索边坡临界滑 动面流程图

3 边坡算例和工程实例分析

3. 1 均质边坡算例分析

为了证明本文构建的计算模型的可行性，本文首先将Zhang[18]采用的三维均质边坡作为算例1,对其进行了分析。

该算例是一个典型的均质边坡，其剖面的几何信息见图9，将该剖面沿其边坡走向拉伸为一个三维均质边坡模型，模型的长、宽、高为80 m×62.7 m×18.3 m，模型的网格划分形式见图10，模型中采用的材料参数见表1。

图9 算例1均质边坡剖面的几何信息

图10 算例1三维均质边坡模型的网格划分形式

表1 算例1三维均质边坡模型的材料参数

对该模型进行分析时，首先利用结合莫尔-库伦准则的光滑有限元法对该均质边坡进行弹塑性分析，得到边坡的应力场;然后利用最大最小蚁群优化算法及滑动面应力法求取边坡的临界滑动面及对应的边坡安全系数，得到该边坡的临界滑动面，见图11。

图11 算例1三维均质边坡的临界滑动面

该算例曾被许多研究人员用来验证自身算法的可行性，采用不同算法得到的该算例边坡安全系数见表2。此外，在边坡的二维剖面上将本文算法所得的边坡临界滑动面与Zhang[18]所得的结果进行了对比，见图12。

表2 不同算法得到的算例1三维均质边坡安全系数的对比

图12 不同算法得到的算例1二维均质边坡临界滑动 面的对比

由表2和图12可知：本文算法与其他几种算法所求得的该边坡临界滑动面形态和边坡安全系数均较为接近。

3. 2 非均质边坡算例分析

本文将一个双层材料边坡作为算例2，对其进行了分析，以验证本文算法对非均质边坡的适用性。

该算例边坡是一个非均质边坡，其剖面的几何信息见图13，将该剖面沿其边坡走向拉伸为一个简单的三维非均质边坡模型，模型的长、宽、高为50 m×45 m×30 m，模型的网格划分形式见图14，模型中采用的材料参数见表3。

图13 算例2非均质边坡剖面的几何信息

图14 算例2三维非均质边坡模型的网格划分形式

表3 算例2三维非均质边坡模型的材料参数

采用与算例1相同的分析方式，得到该边坡的临界滑动面，见图15，可知滑动面形态较为复杂，有较多起伏，表现为图中滑动面上的黑色阴影线条。

图15 算例2三维非均质边坡的临界滑动面

本文利用商业软件GEO-SLOPE使用Bishop法对该边坡的二维剖面进行计算[22]，得到二维极限平衡法即二维剖面上的临界滑动面，同时对三维的临界滑动面在相同剖面中进行切片，得到了一个二维的临界滑动面，见图16。为了验证本文算法的可行性，还利用商业软件ABAQUS对模型使用三维强度折减法得到了该边坡的安全系数。采用不同算法得到的算例2三维非均质边坡安全系数，见表4。

图16 不同算法得到的算例2二维非均质边坡临 界滑动面的对比

表4 不同算法得到的算例2三维非均质边坡安全系数的对比

由图16和表4可知：在二维剖面上，本文算法与二维极限平衡法所得的边坡临界滑动面整体形态比较接近，但其临界滑动面所截取的滑体体积略大一些；三维边坡中的临界滑动面对应的安全系数也与预想一致，较二维极限平衡法的计算结果要偏大一些，与三维强度折减法的计算结果较为接近。

3. 3 实际边坡工程算例分析

为了验证本文算法的工程实用性，将某公路路堑边坡工程实例作为算例3，利用本文算法对其进行了三维边坡临界滑动面搜索与分析。

该滑坡区属构造剥蚀低山河谷地貌，前缘为正在施工的某农村公路工程，区域内河流从东侧水库流入，向西南侧流出，河谷走向与滑坡坡向近乎直交。该区域整体地形受地质构造和河流切割双重作用控制，区域构造为北高南低的单斜构造，滑坡区整体坡角为10°～45°，其中北面山体坡角为35°～50°，山体下部缓坡地带坡角为8°～15°；区域东西向受地形河流深切后，地势受岩性控制呈东高西低、陡缓相间的阶梯状地形。该边坡由于受河谷切割和河谷岸边乡村公路修建开挖的影响，形成了近直立的坡体，坡体主要由崩坡积层和残积层组成，直接裸露于地表，见图17。

图17 某公路路堑边坡滑坡区全貌

图18 算例3实际边坡工程地质剖面图

本节中的计算模型，为该公路路堑边坡中截取的滑坡地带，其内部地层条件比较复杂。该模型的长、宽、高约为64 m×60 m×42 m，模型网格划分形式见图19，模型中采用的材料参数见表5。

图19 算例3三维边坡模型的网格划分形式

表5 算例3三维边坡模型的材料参数

在算例3的实际边坡上，采用本文算法进行分析，得到的三维边坡临界滑动面见图20，图中滑动面上的黑色线条同样表示滑动面中存在的起伏。

图20 算例3三维边坡的临界滑动面

以蚁群算法中第一步所寻找到的平面作为基准面，剖分出一个二维平面作为二维分析的边坡模型，在二维剖面上的边坡临界滑动面见图21，采用不同算法得到的算例3边坡安全系数，见表6。

图21 不同算法得到的算例3二维边坡临界滑动面的对比

表6 不同算法得到的算例3三维边坡安全系数的对比

由图21和表6可知：在该边坡工程实例中，本文算法与二维极限平衡法所得二维边坡临界滑动面的整体形态比较接近，但其临界滑动面所截取的滑体体积要偏大一些；本文算法求得的三维边坡安全系数相比二维极限平衡法偏大，比三维强度折减法略小。

由该边坡的实际滑坡区域(见图22)可以看出，箭头指示区域有变形产生，与本文求得的该边坡临界滑动面位置大致相同，说明该计算结果符合实际的滑坡发展趋势。

图22 某公路路堑边坡的实际滑坡区域

4 结 论

本文主要采用光滑有限元法，并结合最大最小蚁群优化算法和滑动面应力法，构建了搜索三维边坡临界滑动面的计算模型。该计算方法首先使用光滑有限元法计算出边坡的应力场；然后利用最大最小蚁群优化算法搜索出边坡的一系列滑动面，并采用滑动面应力法计算得到蚁群搜索出的滑动面对应的边坡安全系数；最后将边坡安全系数最小的滑动面确定为边坡的临界滑动面。

将本文算法应用于边坡算例和工程实例，结果表明：利用本文算法所搜索到的边坡临界滑动面和边坡安全系数与其他方法较为接近，且在边坡工程实例中，本文算法搜索到的三维边坡临界滑动面与实际的滑坡发展趋势较为吻合，故可以判定，本文将最大最小蚁群优化算法用于三维边坡临界滑动面搜索是可行和可靠的。此外，本文算法形式较为简单，便于程序实现，适用于大部分缓层土质三维边坡，在根据边坡的形态、构造特征给定滑入区和滑出区后，便能得到边坡中的任意形态的空间滑动面及对应的安全系数，可为缓层土质边坡的开挖和支护设计提供依据。