基于改进差分进化算法的三维边坡锚固位置优化

尹志凯; 卢坤林; 石峰; 朱大勇

doi:10.11779/CJGE202007016

基于改进差分进化算法的三维边坡锚固位置优化 English Version

1.
合肥工业大学土木与水利工程学院,安徽　合肥　230009
2.
浙大宁波理工学院土木建筑工程学院,浙江　宁波　315100

基金项目:

国家自然科学基金项目 41402256

详细信息

作者简介:
尹志凯(1992—),男,硕士研究生,主要从事边坡稳定性和三维锚固优化方面的研究工作。E-mail：yin.zk@formail.com

通讯作者:
卢坤林, E-mail：lukunlin@hfut.edu.com

中图分类号: TU43
计量
- 文章访问数: 249
- HTML全文浏览量: 19
- PDF下载量: 138
出版历程
- 收稿日期: 2019-10-06
- 网络出版日期: 2022-12-05
- 刊出日期: 2020-06-30

Optimization of anchorage location of three dimensional slopes based on improved differential evolution algorithm

1.
School of Civil and Hydraulic Engineering, Hefei University of Technology, Hefei 230009, China
2.
School of Civil and Architectural Engineering, Ningbo Institute of Technology, Zhejiang University, Ningbo 315100, China

摘要

摘要: 基于改进的差分进化算法,在三维锚固力计算基础上考虑锚固力的三维空间分布状态,对三维边坡锚固位置进行合理优化,从而有效节约锚杆数量。利用集中力作用下的半无限弹性空间内一点应力解答,将锚固力等效叠加到边坡滑面正应力和剪应力上,由此推导了关于锚固位置坐标的锚固力函数,并以此作为锚固位置优化的目标函数。对于标准差分进化算法（DE）的改进,将自适应差分进化算法与一般约束乘子算法（PHR）结合,使算法能够解决锚固边坡约束优化问题。算例和工程实例表明,在边坡加固安全系数相同且单根锚杆受力保持不变的情况下,通过本文方法优化锚固位置能节约20%～25%的锚固点数量,该方法可为工程应用提供参考。
- 边坡 /
- 极限平衡 /
- 三维稳定性 /
- 安全系数
Abstract: Based on the improved differential evolution algorithm, considering the 3D spatial distribution of anchoring force through calculation of the 3D anchoring force, the anchorage position of 3D slope is optimized reasonably to save the number of bolts. By using the stress solution of one point in the semi-infinite elastic space under the action of concentrated force, the anchoring force is superimposed on the normal stress and shear stress of the slope surface. Then the function for anchoring force about the coordinates of the anchorage position is derived, and taken as the objective function of the optimization of the anchorage position. For the improvement of the standard differential evolution algorithm (DE), the adaptive differential evolution algorithm is combined with the general constraint multiplier algorithm (PHR) to solve the problem of constraint optimization of anchoring slopes. Numerical and engineering examples show that with the same safety factor of slope reinforcement and unchanged single anchoring force, through the optimization of anchorage position, the number of anchorage points can be saved by 20% ~ 25%, and the proposed method can provide reference for engineering application.
- slope /
- limit equilibrium /
- 3D stability /
- safety factor

HTML全文

0. 引言

锚杆（索）加固作为一种常用的高边坡稳定加固技术,广泛应用于各专业工程的边坡安全治理中。锚固施工技术和施工方案在锚固边坡理论研究的发展中不断改进,其中锚固力等效模型和锚固参数优化是边坡锚固理论研究的两个重要方向,也是国内外学者研究的热点问题。锚固边坡稳定性分析在极限平衡理论基础上展开了研究,二维稳定分性分析现已建立其统一的求解体系^[1]。基于极限平衡的锚固力计算方法是将锚固力作为集中力来处理并对锚固力计入方式进行探讨^[2-8]。雷军等^[9]基于弹性力学Boussinesq解和Cerruti解考虑了锚拉力扩散效应,实现了二维锚固边坡分析方法的另一种途径。

三维锚固边坡稳定性分析需要开展三维锚固力等效模型研究。近年来,卢坤林等基于经典土力学假设建立了锚固力等效三维滑面正应力模型,进行锚固边坡稳定性分析,并应用于工程实践取得了良好的效果^[10]。锚固边坡参数优化通常借助最优化方法进行研究,对锚固边坡的锚固角度、锚固长度、锚固间距、锚杆直径等锚固边坡设计要素进行优化,可以得到更合理的加固方案^[11-14],并且能够提高锚杆（索）加固边坡的锚固效果和经济性能。遗传算法、蚁群算法在边坡工程中已被应用于搜索边坡最危险滑面^[15]。差分进化算法适用性广,算法简单,易于编程改进,适用不同领域工程优化问题^[16-17]。

针对锚固边坡优化问题,本文将结合弹性力学理论^[18],建立一种考虑锚固力空间效应的线性方程组,求解带锚固位置坐标的锚固力函数。通过改进差分进化算法,对锚固边坡加固位置进行优化,保证在边坡加固后安全系数相同且单根锚杆受力基本不变的情况下,能够节约一定锚杆用量,提高锚固经济性。

1. 锚固力等效滑面内力

1.1 半无限空间内一点的应力状态

如图1（a）所示,预应力锚固作用效果可以等效于集中力P作用下边坡体内各点的应力状态。本文研究的一般边坡可近似看成半无限弹性空间体。如图1（b）,弹性体内某一点M(x,y,z)的6个应力分量 $σ_{x}$ , $σ_{y}$ , $σ_{z}$ , $τ_{x y}$ , $τ_{y z}$ , $τ_{x z}$ 在弹性体微分面上的分布状态,R为该点到作用点的距离。

图 1 弹性空间体内一点的应力状态

Figure 1. Stress state at a point in elastic space

下载: 全尺寸图片幻灯片

如图2,将空间内一点M的应力状态等效四面微分体的受力状态,微分面ABC的外法线为N,方向余弦为（l,m,n）。在研究边坡锚固问题时,需要对锚入角度 $α$ 进行坐标变换,以便将锚固力叠加到滑面应力上。若原坐标系绕y轴旋转 $α$ 角后得到的新坐标轴为Ox'y'z',转换后的坐标为（ $z \sin α + x \cos α$ ,y, $z \cos α - x \sin α$ ）。

图 2 斜微分面的应力状态

Figure 2. Stress state of oblique differential plane

下载: 全尺寸图片幻灯片

新坐标轴对于原直角坐标系各坐标轴的方向余弦如表1所示。

表 1 新坐标轴的方向余弦

Table 1. Direction cosines of new axis

坐标	x	y	z
x'	l₁	m₁	n₁
y'	l₂	m₂	n₂
z'	l₃	m₃	n₃

下载: 导出CSV

| 显示表格

设锚固力作用点坐标为(t₁,t₂,t₃),结合空间内一点的应力状态,得到带位置关系T_i（i=0～15）的滑面一点应力分布状态,T_i的系数表达式如式（2）,其中的 $σ_{x i}$ , $σ_{x i}$ , $σ_{x i}$ , $τ_{x y i}$ , $τ_{y z i}$ , $τ_{x z i}$ 为经典Boussinesq解答^[18]推导式。

(1)

锚固力在滑面上一点的应力状态为

$\begin{array}{l} σ_{n x} = \sum_{i = 0}^{15} σ_{n x i} T_{i}, σ_{n y} = \sum_{i = 0}^{15} σ_{n y i} T_{i}, σ_{n z} = \sum_{i = 0}^{15} σ_{n z i} T_{i}, \\ τ_{n x y} = \sum_{i = 0}^{15} τ_{n x y i} T_{i}, τ_{n y z} = \sum_{i = 0}^{15} τ_{n y z i} T_{i}, τ_{n x z} = \sum_{i = 0}^{15} τ_{n x z i} T_{i} 。 \end{array}}$

(2)

式中 T₀=1；T₁=t₁；T₂=t₂；T₃=t₃；T₄=t₁t₂；T₅=t₁t₃；T₆=t₂t₃；T₇=t₁²；T₈=t₂²；T₉=t₃²；T₁₀=t₁²t₃；T₁₁=t₂²t₃；T₁₂=t₁t₃²；T₁₃=t₂t₃²；T₁₄=t₁³；T₁₅=t₁t₂t₃。

锚固边坡滑面上一点全应力沿3个坐标轴方向的分量为

$\begin{matrix} x_{n} = l σ_{n x} + m τ_{n x y} + n τ_{n x z}, \\ y_{n} = l τ_{n x y} + m σ_{n y} + n τ_{n x z}, \\ z_{n} = l τ_{n z x} + m τ_{n y z} + n σ_{n z} 。 \end{matrix}}$

(3)

1.2 滑面正应力增量模型

将 $x_{n}$ , $y_{n}$ 和 $z_{n}$ 投影到法线N上,得到微分面ABC上的正应力增量 $σ_{n}$ ,并由式（1）～（3）可以建立滑动面上某一点正应力增量关于锚固位置坐标的滑面正应力增量函数：

$σ_{n} = l x_{n} + m y_{n} + n z_{n} = \sum_{i = 0}^{15} σ_{n i} T_{i}$ ,

(4)

式中, $σ_{n i} = l^{2} σ_{n x i} + m^{2} σ_{n y i} + n^{2} σ_{n z i} + 2 l m τ_{n x y i} + 2 n l τ_{n x z i} +$ $2 m n τ_{n y z i}$ 。

1.3 滑面剪应力增量模型

将锚固力作用产生的应力等效到垂直作用于水平微分面上的一点对该平面各处的剪切效应,同理,由式（1）～（3）得到滑面上某一点剪应力增量关于锚固位置坐标的滑面剪应力增量函数：

$τ_{n} = \sum_{i = 0}^{15} [l m (σ_{n y i} - σ_{n x i}) + (l^{2} - m^{2}) τ_{x y i}] T_{i}$ 。

(5)

2. 三维边坡锚固力的目标函数

2.1 概化模型

如图3（a）所示的锚固边坡与三维滑面,其中坡面方程为 $z = g (x, y)$ ,滑面方程为 $z = s (x, y)$ ,土体抗剪强度指标为 $N = σ d A$ 和 $S = τ d A$ ,重度为 $d A$ ,则滑体的单位重量为 $w = γ (g - s)$ 。图3（b）为锚固边坡主断面受力模型,锚固力为P,滑面正应力和剪应力分别为 $σ$ 和 $τ$ ,滑面上一点正应力与剪应力方向余弦分别为

图 3 概化模型

Figure 3. Generalized model

下载: 全尺寸图片幻灯片

$\begin{array}{l} (n_{σ}^{x}, n_{σ}^{y}, n_{σ}^{z}) = (l, m, n), \\ (n_{τ}^{x}, n_{τ}^{y}, n_{τ}^{z}) = (1 / Δ', 0, s_{x} / Δ^{'}) 。 \end{array}}$

(6)

式中, $Δ = \sqrt{1 + s_{x}^{2} + s_{y}^{2}}$ , $Δ' = \sqrt{1 + s_{x}^{2}}$ 。

2.2 基本假设

（1）锚固边坡滑面形式为固定滑面。

（2）滑动方向与oxz面平行,即滑面剪应力与滑动方向平行。

（3）沿用卢坤林等^[10]的假设

$σ = σ_{0} (λ_{1} f_{1} + λ_{2} f_{2} + λ_{3} f_{3})$ ,

(7)

式中 $λ_{1}$ , $λ_{2}$ , $λ_{3}$ 为待定参数；f₁,f₂,f₃为已知函数^[19]。则锚固边坡滑面正应力为

$σ = σ_{0} + σ_{n} λ_{p} (n - l s_{x} - m s_{y}) / Δ^{2}$ ,

(8)

式中, $σ_{0} = w [1 + (λ_{1} + λ_{2} s_{x} + λ_{3} s_{y}) (g' - s')] / Δ^{2}$ , $λ_{p}$ 为锚固力系数。

2.3 三维锚固平衡方程组

考虑3个力平衡和一个绕垂直滑动方向（y轴）的力矩平衡方程得

$\begin{array}{l} \iint σ m d A = - λ_{p} \iint σ_{n} m d x d y, \\ \iint (σ l + τ / Δ^{'}) d A = - λ_{p} \iint σ_{n} l d x d y, \\ \iint (σ n + τ s_{x} / Δ^{'}) d A - W = - λ_{p} \iint σ_{n} n d x d y, \\ \iint (σ n + τ s_{x} / Δ^{'}) x - (σ l + τ / Δ^{'}) s d A = M + λ_{P} M_{P} 。 \end{array}}$

(9)

式中,W= $\iint w d x d y$ ,M_P= $\iint σ_{n} (n x + l z) d x d y$ ；M= $\iint w x d x d y$ 。

根据莫尔-库仑准则,有

$τ = (c + σ \tan φ + τ_{n}) / F_{s}$ 。

(10)

对于锚固后安全系数已知的情况下,将式（1）～（8）,（10）代入式（9）简化为求解4个未知量 $λ_{p}$ , $λ_{1}$ , $λ_{2}$ , $λ_{3}$ 的线性方程组：

$\begin{array}{l} λ_{1} A_{11} + λ_{2} A_{12} + λ_{3} A_{13} + λ_{p} A_{14} = B_{1}, \\ λ_{1} A_{21} + λ_{2} A_{22} + λ_{3} A_{23} + λ_{p} A_{24} = B_{2}, \\ λ_{1} A_{31} + λ_{2} A_{32} + λ_{3} A_{33} + λ_{p} A_{34} = B_{3}, \\ λ_{1} A_{41} + λ_{2} A_{42} + λ_{3} A_{43} + λ_{p} A_{44} = B_{4} 。 \end{array}}$

(11)

式中系数A₁₁～A₁₃,A₂₁～A₂₃,A₃₁～A₃₃,B₁～B₄不受锚固力增量影响,与文献[10]一致；A₁₄叠加了相应的滑面正应力增量,A₂₄,A₃₄,A₄₄分别叠加了相应的滑面正应力增量和剪应力增量：

$\begin{array}{l} A_{14} = \iint A_{1} σ_{n} d x d y, \\ A_{24} = \iint A_{2} σ_{n} d x d y + \iint A_{5} τ_{n} d x d y, \\ A_{34} = \iint A_{3} σ_{n} d x d y + \iint A_{6} τ_{n} d x d y, \\ A_{44} = \iint A_{4} σ_{n} d x d y + \iint A_{7} τ_{n} d x d y 。 \end{array}}$

(12)

令 $s_{1} = - s_{x} + \tan φ \frac{Δ}{Δ^{'}} / F, s_{2} = 1 + s_{x} \tan φ \frac{Δ}{Δ^{'}} / F, s_{3} =$ $\frac{Δ}{F_{s} Δ^{'}} g^{'} - s^{'}$ ,则式（12）中, $A_{1} = \frac{s_{y}}{Δ} (1 + \frac{s_{x}^{2} + s_{y}^{2} - 1}{Δ^{2}})$ , $A_{2} = \frac{s_{x}}{Δ} + s_{1} \frac{s_{x}^{2} + s_{y}^{2} - 1}{Δ}$ , $A_{3} =$ $\frac{1}{Δ} +$ $s_{2} \frac{s_{x}^{2} + s_{y}^{2} - 1}{Δ^{3}}$ , $A_{4} = (\frac{1}{Δ} x^{'} + \frac{s_{x}}{Δ} z^{'}) + (x^{'} + s s_{x} -$ $s_{x} x^{'} + s)$ $s_{3} \frac{s_{x}^{2} + s_{y}^{2} - 1}{Δ^{3}}$ , $A_{5} = \frac{Δ}{F_{s} Δ^{'}}$ , $A_{6} = s_{x} \frac{Δ}{F_{s} Δ^{'}}$ , $A_{7} = (s_{x} x^{'} - s) \frac{Δ}{F_{s} Δ^{'}}$ 。

2.4 目标优化函数

三维边坡锚固位置优化的目标函数实际为包含位置参数的锚固力系数表达式,因此需求解有4个未知变量的线性方程组式（11）,由克拉默法则可得考虑锚固力空间分布的锚固系数 $λ_{p}$ 的表达式：

$λ_{p} = D_{n 4} / D_{0}$ ,

(13)

式中,

$D_{0} = | \begin{matrix} A_{11} & A_{12} & A_{13} & A_{14} \\ A_{21} & A_{22} & A_{23} & A_{24} \\ A_{31} & A_{32} & A_{33} & A_{34} \\ A_{41} & A_{42} & A_{43} & A_{44} \end{matrix} |, D_{n 4} = | \begin{matrix} A_{11} & A_{12} & A_{13} & B_{1} \\ A_{21} & A_{22} & A_{23} & B_{2} \\ A_{31} & A_{32} & A_{33} & B_{3} \\ A_{41} & A_{42} & A_{43} & B_{4} \end{matrix} | 。$

将系数行列式按第4列展开得

$D_{0} = - A_{14} D_{1} + A_{24} D_{2} - A_{34} D_{3} + A_{44} D_{4}$ ,

(14)

式中,

将式（4）,（5）,（12）代入到式（13）得到锚固位置优化目标函数表达式：

$f {(x)}_{λ_{p}} = \frac{D_{n 4}}{\sum_{i = 0}^{15} \iint [s_{4} σ_{n} + (A_{5} - A_{6}) τ_{n}] d x d y}$ ,

(15)

式中, $s_{4} = - A_{1} D_{1} + A_{2} D_{2} - A_{3} D_{3} + A_{4} D_{4}$ 。

编程求解各参数得到关于锚固位置坐标的锚固力函数,即为本文优化问题所需的目标函数。

3. 锚固优化自适应差分进化算法

3.1 自适应差分进化算法

对于DE算法而言,初始种群数量,迭代次数,变异算子,交叉算子,选择算子以及变异策略的选择会不同程度影响算法的寻优能力。在DE算法进化的不同时期选择合理的交叉、变异参数可以解决算法早熟和后期收敛速度慢的问题,同时针对锚固边坡约束条件和目标函数适应度优化变异策略,使算法精确地逼近最优解。

（1）群体初始化

在解空间内随机产生初始种群：

$x (i, :) = l b + (u b - l b) . * r a n d (s i z e (l b))$ ,

(16a)

$F i t n e s s (i) = F u n (x (i, :))$ 。

(16b)

式中,lb,ub分别是位置坐标的一维取值边界,rand表示在0到1之间的随机数,D为解空间维度。

（2）自适应差分变异模式

不同的差分进化策略,进化过程和全局搜索能力不尽相同。为了弥补单一进化模式的缺陷,进化前期以概率ω采用DE/rand/1/bin变异模式提高DE的全局搜索能力避免算法陷入局部最优,进化后期以概率1-ω按DE/best/1/bin模式变异,改善DE后期迭代精度和收敛速度。因此,ω随进化代数线性递增满足自适应变异模式需求。

$v_{i j} (g + 1) = {\begin{cases} x_{r 1 j} (g) + f (x_{r 2 j} (g) - x_{r 3 j} (g)) (rand ＜ ω) \\ x_{bestj} (g) + f (x_{r 1 j} (g) - x_{r 2 j} (g)) (rand \geq ω) \end{cases}$ 。

(17)

其中,r₁≠r₁≠r₁≠i且在种群中随机选择,变异算子F取值为0～2,用于控制偏差的放大作用。

（3）交叉操作

交叉操作保证了种群的多样性,利用交叉因子C_r控制D维空间内各子代个体的活跃度,C_r取值区间为（0,1）,一般C_r取0.3～0.6。

$u_{i j} (g + 1) = {\begin{cases} v_{i j} (g + 1) (rand (0, 1) \leq C_{r}^{}) \\ x_{i j} (g) (其它) \end{cases}$ 。

(18)

（4）选择操作

计算每次迭代个体的适应度函数的值,通过选择操作选择适应度更高的个体作为新的子代。选择操作：

$x_{i} (g + 1) = {\begin{cases} u_{i} (g + 1) (f (u_{i} (g + 1)) \leq f (x_{i} (g))) \\ x_{i} (g) (其它) \end{cases}$ 。

(19)

3.2 自适应调整参数

（1）参数C_r的自适应调整

C_r的取值大小决定了种群的多样性,C_r取值较小会导致种群多样性减小,子代参与度不高,算法容易收敛于局部最优解。C_r取值过大,收敛速度快,且易造成群体扰动反而导致收敛变慢。适应度值即为目标函数值,设置当代种群自适应交叉系数C_i,C_i的值随着子代中最优与最劣适应度值的差值改变而相应调整。 $C_{r_{1}}$ 、 $C_{r_{u}}$ 为C_r的最小值和最大值。因此,可通过式（20）策略调整交叉参数。

$C_{r_{i}} = {\begin{cases} C_{r_{1}} + C_{i} (f_{i} - f_{\min}) (f_{i} \geq \overline{f)} \\ C_{r} (其它) \end{cases}$ 。

(20)

（2）参数F的自适应调整

变异参数F的大小对算法收敛速度有直接的影响。自适应调整变异参数,在不断地迭代进化过程中,适当降低变异程度,能够保存优良的种群信息,提高算法后期的寻优能力。设置自适应算子 $λ$ ,G_m为最大迭代次数,G为当前迭代次数,F₀为变异系数的初始值,则有

$F_{i} = F_{0} 2^{λ}$ ,

(21)

式中, $λ = e^{u}$ ,u=(1-G)/(1+G_m-G)。

3.3 锚固约束处理

对于某一确定的滑坡体,锚固区域需要满足一定的边界条件（lb,ub）,并且锚杆加固应保持合适的间距避免发生群锚效应。根据锚杆（索）设计规范,对锚固边坡锚杆布置要求、锚固间距d以及各锚固位置横向间距v_d、纵向间距h_d等约束条件^[20],对应的约束方程为

$\begin{matrix} C_{1} (x) = - {(x_{i} - x_{j})}^{2} + 2 v_{d} \leq 0, \\ C_{2} (x) = - {(x_{m} - x_{n})}^{2} + 2 h_{d} \leq 0, \\ C_{3} (x) = x_{i} + l b \leq 0, \\ C_{4} (x) = - x_{i} - u b \leq 0, \\ C_{5} (x) = x_{m} + l b \leq 0, \\ C_{6} (x) = - x_{m} - u b \leq 0, \\ C_{7} (x) = - {(x_{i} - x_{j})}^{2} - {(x_{m} - x_{n})}^{2} + d^{2} \leq 0, \end{matrix}}$

(22)

式中, $(i, j = 1, 3, 5, \dots, i \neq j ； m, n = 2, 4, 6, \dots, m \neq n) 。$

标准差分进化算法往往不适用于解决带约束的优化问题,现通过改进DE,将算法程序中植入约束判断程序以此处理锚固边坡约束。一般约束问题乘子法（PHR算法）克服了罚函数法的病态性质,收敛速度快。通过PHR方法来处理锚固约束条件,由于锚固边坡参数约束条件中没有等式约束,故以下仅考虑不等式约束。构造的增广Lagrange函数为

$L (x, λ, σ) = f (x) + \frac{1}{2 σ} \sum_{i = 1}^{l} {{[\max (0, λ_{i} - σ c_{i} (x))]}^{2} - λ_{i}^{2}} 。$

(23)

Lagrange罚函数乘子向量 $λ$ 和罚参数向量 $σ$ 一般难以确定,现采用最优化方法 $λ$ 、 $σ$ 进行初始化和修正。首先设置初始点 $x_{0}$ 、初始乘子向量 $λ_{1}$ 、初始惩罚因子 $σ_{1}$ ,放大系数c（c>0）以及小于1的正数 $θ$ 。乘子修正公式为

${(λ_{k + 1})}_{i} = \max [0, {(λ_{k})}_{i} - σ c_{i} (x_{k})] (i = 1, 2, 3, \dots, l)$ 。

(24)

终止判断准则由 $φ_{k}$ 和控制误差 $ε$ 决定,具体判断标准为

$φ_{k} = {\sum_{i = 1}^{l} {[\max (c_{i} (x_{k}), \frac{{(λ_{k})}_{i}}{σ})]}^{2}}^{0.5} \leq ε$ 。

(25)

PHR算法步骤：①给定初始值[ $x_{0}$ , $λ_{1}$ , $σ_{1}$ ,c, $θ$ ]令k=1；②以x_k-1为初始点,不考虑约束条件下构建Lagrange增广函数,求解 $\min L (x, λ, σ)$ 得到最优解{x_k}；③计算 $φ_{k}$ ,若 $φ_{k} \leq ε$ ,则{x_k}为当前最优种群,进入下一步差分进化,否则转④；④当 $φ_{k} / φ_{k - 1} \leq θ$ 时,转⑤,否则令 $σ_{k + 1} = c σ_{k}$ ,转⑤；⑤修正乘子向量 $λ_{k}$ ,

${(λ_{k + 1})}_{i} = \max [0, {(λ_{k})}_{i} - σ c_{i} (x_{k})] i = 1, 2, 3, \dots, l$ 。令k=k+1,转②。

3.4 算法流程图

算法流程图如图4所示。

图 4 算法流程图

Figure 4. Flow chart of algorithm

下载: 全尺寸图片幻灯片

4. 算例分析

已知某均质干燥边坡,坡面与水平方向夹角为45°,坡高h为10 m,坡面宽度b为10 m。土层参数为：土体重度 $γ$ =18 kN/m³,边坡土体抗剪强度指标c=33 kPa, $φ$ =29°。计算模型如图5所示。坡面方程z=x+17.32,假设此均质边坡的三维滑动体为椭球,滑面方程为

图 5 计算模型简图

Figure 5. Schematic diagram of computational model

下载: 全尺寸图片幻灯片

$\frac{{(x + 17.32)}^{2}}{20^{2}} + \frac{y^{2}}{14^{2}} + \frac{{(z - 20)}^{2}}{20^{2}} = 1$ 。

(26)

经三维极限平衡方法^[19]求得该边坡稳定系数为1.23,现拟采用锚杆对该边坡进行加固,加固安全系数为1.30。根据岩土边坡锚固规范^[20]设定锚固参数的取值及范围：锚杆与水平方向夹角为15°,锚杆布置横纵间距分别为3.0,2.5 m,锚固点横纵坐标的取值范围为-15.32≤x≤-5.82,-3.5≤y≤3.5。

图6为本文算法优化过程中目标函数值与迭代次数曲线图,差分进化代数为300代,由于算例较为简单约进化到120代时收敛到最优状态。本文方法优化后需要12根锚杆,单根锚杆受力为391.3 kN,所需总锚固力为4695.2 kN,优化布置结果如图7所示。若采用常规加固方案使边坡安全系数达到1.3,在合理锚固范围内需要布设15根锚杆,单根锚杆受力为390.5 kN。

图 6 目标函数迭代曲线

Figure 6. Iterative curve of objective function

下载: 全尺寸图片幻灯片

图 7 锚固方案对比图

Figure 7. Comparison of anchorage schemes

下载: 全尺寸图片幻灯片

如图7所示,与本文方法相比,在安全系数和单根锚杆受力基本不变的情况下,本文优化方案节省3根锚杆。优化后坡体中部和坡体上部锚杆加固较为集中,可能因为均质滑体重心在坡面中心以上,重心部位锚杆布设相对集中有利于提高锚固效果。

5. 工程应用

5.1 工程概况

G216线民丰—克里雅山口公路建设项目K102+230—K104+294段为中山地貌,本项目推荐线K102+000—K127+000段以及比较方案BK101+000—BK110+200展布于中山带,基岩或直接裸露,或被粉砂所覆盖,形成山地草甸,降雨较多,植被繁茂,是该地区主要牧场分布区。K102—K104段土体材料主要为分为两层,上层为灰岩,灰岩层的土体材料重度 $γ$ =25.5 kN/m³,黏聚力c=40 kPa,内摩擦角ϕ=45°；下层为钙质砂岩,钙质砂岩的土体材料性质为重度 $γ$ =22 kN/m³,黏聚力c=30 kPa,内摩擦角ϕ=30°。原始边坡形态如图8所示,K102+230—294段左侧边坡防护平面图,以桩号K260所在位置为原点坐标,K260—K280方向为y轴正方向,过桩号K260垂直线路方向为x轴。K102+240—K102+280段左侧原始边坡可近似看成关于K102+260对称。

图 8 K102+230—K102+294段左侧边坡防护平面图

Figure 8. Plan of slope protection at left side of sections K102+240-K102+280

下载: 全尺寸图片幻灯片

5.2 确定三维滑体方程

采用本课题组分片光滑技术提取二维、三维最危险滑面^[21]。坡体各特征点高度来源于原始勘测数据,分段插值不影响坡体整体形态,与原始坡体保持良好的一致性,可直接应用于理论计算。滑面提取结果如表2,图9（a）,（b）分别为K102+230—K102+294段三维边坡坡面计算模型和三维滑面计算模型。

表 2 滑面及滑体方程提取

Table 2. Extraction of equations for sliding surface and sliding body

计算断面	滑面方程提取
K102+240	${(x + 1.71)}^{2} + {(y - 16.74)}^{2} = {16.96}^{2}$
K102+240	坐标点提取：(13.05,10.28),(17.10,12.40),(20.57,15.05),(23.21,19.05),(24.91,24.10）
K102+260	${(x + 18.86)}^{2} + {(y - 44.65)}^{2} = {45.55}^{2}$
K102+260	坐标点提取：(12.12,10.24),(19.68,14.12),(25.24,18.50),(31.02,25.34),(34.70,31.96)
K102+280	${(x + 8.41)}^{2} + {(y - 46.80)}^{2} = {47.17}^{2}$
K102+280	坐标点提取：(13.04,10.44),(20.16,14.04)、(26.24,18.56),(31.80,25.06),(35.54,32.04)
滑体方程	$\frac{x^{2}}{900} + \frac{y^{2}}{1089} + \frac{{(z - 21.72)}^{2}}{471.76} = 1$

下载: 导出CSV

| 显示表格

图 9 边坡三维计算模型

Figure 9. 3D model for slope

下载: 全尺寸图片幻灯片

5.3 自适应DE优化计算

经三维极限平衡方法计算该边坡安全系数为2.29,综合考虑该段边坡需达到的稳定性和安全性,加固后安全系数须达到2.5。K102+260断面边坡为滑体中性面,计算简图如图10所示,AC段坡面做相应削坡处理,AB为边坡平台,滑体最大高度为21.72 m。根据提取的滑体方程和计算简图得到算法优化必要条件见表3。

图 10 K102+260断面边坡计算简图

Figure 10. Calculation of K102+260 section slope

下载: 全尺寸图片幻灯片

表 3 相关方程及初始参数

Table 3. Equations and initial parameters

方程	初始参数
滑面方程	$z = 21.72 - {({21.72}^{2} - 0.773 x^{2} - 0.524 y^{2})}^{0.5}$
坡面方程	z=0.879x(0 m<x<8.986 m)
	z=7.92(8.986 m<x<11 m)
	z=1.31x-5.28(11 m<x<20 m)
	z=0.4x+11.72(x>20 m)
土层方程	土层1：(0<z<0.2143x+9.66)
土层方程	土层2：(0.2143x+9.66<z<0.8793x)
参数初始化	种群规模NP=30,初始变异系数0.5, $c r_{l}$ =0.1, $c r_{u}$ =0.6；迭代次数G=8000
锚固优化参数取值	锚固角度15°；v_d≥2.5 m,h_d≥2.5 m；2.5 m≤x≤16.5 m,-26 m≤y≤23 m

下载: 导出CSV

| 显示表格

利用本文第3节方法编程计算得到单锚作用下的锚固力函数 $f (x, y, z)$ ,设打入锚杆数量为n,根据叠加原理得到目标函数表达式：

$\min F (Q) = \sum_{i = 1}^{n} f (x_{i}, y_{i}, z_{i})$ 。

(27)

通过本文优化方法使该边坡安全系数达到2.5。差分进化算法优化过程如图11,总进化代数为8000,算法进化4000代时基本收敛于最优解。优化结果如图12,需要的锚固点数量n=67,单根锚杆施加锚固力为489.2 kN,总锚固力为32773.6 kN。

图 11 DE目标函数迭代曲线

Figure 11. Iterative curve of DE objective function

下载: 全尺寸图片幻灯片

图 12 三维边坡锚固位置方案

Figure 12. Optimization scheme of anchorage position of 3D slope

下载: 全尺寸图片幻灯片

5.4 优化方法与常规方法对比

常规方案采用4 m×3 m框格布置,加固后安全系数达到2.5,单根锚杆受力为485.5 kN,需要打入90根锚杆。采用本文方法优化同样满足安全系数为2.5,单根锚杆施加锚固力为489.2 kN与常规加固方案基本保持一致,但通过对锚固位置合理优化调整后仅需67根锚杆便能达到防护要求。

具体布置方案对比如图12所示,该边坡对应里程桩号为K102+240—K102+280区域内,锚杆加固较为密集,图12锚固加强区域,计划打入40根锚杆,约占总锚杆数量的59.7%,可能由于此区域对应三维滑体中心区域,为滑体质量集中的部位,与其它区域相比锚固能取得更好的效果。因此在满足安全系数为2.5和单根锚杆受力基本保持一致的情况下,采用本文方法布设锚杆可节约23根锚杆,大约比常规方案节省25.5%的锚杆。

6. 结论

（1）提出了一种考虑锚固力空间分布的三维锚固力分析方法。在三维锚固力计算方法基础上结合空间内一点的应力状态,推导了既定锚固安全系数下,固定滑面的三维边坡锚固力计算方法,对三维边坡锚固工程的锚固力计算提供了理论基础。

（2）实现了固定滑面下基于改进差分进化算法的锚固体位置优化方法。通过改进标准差分进化算法,加入锚固约束及对算法变异策略和参数自适应调整,使改进的差分进化算法适用于三维边坡锚固位置的优化,编程优化计算简便。

（3）算例和工程实例验证了本文方法对锚固边坡优化的适用性,在锚固安全系数相同且单根锚固体受力基本保持不变的情况下,通过本文方法可节约20%～25.5%的锚杆。

图 1 弹性空间体内一点的应力状态

Figure 1. Stress state at a point in elastic space

下载: 全尺寸图片幻灯片

图 2 斜微分面的应力状态

Figure 2. Stress state of oblique differential plane

下载: 全尺寸图片幻灯片

图 3 概化模型

Figure 3. Generalized model

下载: 全尺寸图片幻灯片

图 4 算法流程图

Figure 4. Flow chart of algorithm

下载: 全尺寸图片幻灯片

图 5 计算模型简图

Figure 5. Schematic diagram of computational model

下载: 全尺寸图片幻灯片

图 6 目标函数迭代曲线

Figure 6. Iterative curve of objective function

下载: 全尺寸图片幻灯片

图 7 锚固方案对比图

Figure 7. Comparison of anchorage schemes

下载: 全尺寸图片幻灯片

图 8 K102+230—K102+294段左侧边坡防护平面图

Figure 8. Plan of slope protection at left side of sections K102+240-K102+280

下载: 全尺寸图片幻灯片

图 9 边坡三维计算模型

Figure 9. 3D model for slope

下载: 全尺寸图片幻灯片

图 10 K102+260断面边坡计算简图

Figure 10. Calculation of K102+260 section slope

下载: 全尺寸图片幻灯片

图 11 DE目标函数迭代曲线

Figure 11. Iterative curve of DE objective function

下载: 全尺寸图片幻灯片

图 12 三维边坡锚固位置方案

Figure 12. Optimization scheme of anchorage position of 3D slope

下载: 全尺寸图片幻灯片

表 1 新坐标轴的方向余弦

Table 1 Direction cosines of new axis

坐标	x	y	z
x'	l₁	m₁	n₁
y'	l₂	m₂	n₂
z'	l₃	m₃	n₃

下载: 导出CSV

表 2 滑面及滑体方程提取

Table 2 Extraction of equations for sliding surface and sliding body

计算断面	滑面方程提取
K102+240	${(x + 1.71)}^{2} + {(y - 16.74)}^{2} = {16.96}^{2}$
K102+240	坐标点提取：(13.05,10.28),(17.10,12.40),(20.57,15.05),(23.21,19.05),(24.91,24.10）
K102+260	${(x + 18.86)}^{2} + {(y - 44.65)}^{2} = {45.55}^{2}$
K102+260	坐标点提取：(12.12,10.24),(19.68,14.12),(25.24,18.50),(31.02,25.34),(34.70,31.96)
K102+280	${(x + 8.41)}^{2} + {(y - 46.80)}^{2} = {47.17}^{2}$
K102+280	坐标点提取：(13.04,10.44),(20.16,14.04)、(26.24,18.56),(31.80,25.06),(35.54,32.04)
滑体方程	$\frac{x^{2}}{900} + \frac{y^{2}}{1089} + \frac{{(z - 21.72)}^{2}}{471.76} = 1$

下载: 导出CSV

表 3 相关方程及初始参数

Table 3 Equations and initial parameters

方程	初始参数
滑面方程	$z = 21.72 - {({21.72}^{2} - 0.773 x^{2} - 0.524 y^{2})}^{0.5}$
坡面方程	z=0.879x(0 m<x<8.986 m)
	z=7.92(8.986 m<x<11 m)
	z=1.31x-5.28(11 m<x<20 m)
	z=0.4x+11.72(x>20 m)
土层方程	土层1：(0<z<0.2143x+9.66)
土层方程	土层2：(0.2143x+9.66<z<0.8793x)
参数初始化	种群规模NP=30,初始变异系数0.5, $c r_{l}$ =0.1, $c r_{u}$ =0.6；迭代次数G=8000
锚固优化参数取值	锚固角度15°；v_d≥2.5 m,h_d≥2.5 m；2.5 m≤x≤16.5 m,-26 m≤y≤23 m

下载: 导出CSV

参考文献(21)

[1]	郑颖人, 杨明成. 边坡稳定安全系数求解格式的分类统一[J]. 岩石力学与工程学报, 2004, 23(16): 2836-2841. doi: 10.3321/j.issn:1000-6915.2004.16.031 ZHENG Ying-ren, YANG Ming-cheng. Classification unification of calculation formulation for safety factor of slope stability[J]. Chinese Journal of Rock Mechanics and Engineering, 2004, 23(16): 2836-2841. (in Chinese) doi: 10.3321/j.issn:1000-6915.2004.16.031
[2]	康景文, 朱大勇. 一种边坡锚固力新的计算方法[J]. 四川建筑科学研究, 2005, 31(5): 72-74. https://www.cnki.com.cn/Article/CJFDTOTAL-ACZJ200505021.htm KANG Jing-wen, ZHU Da-yong. A new method for calculating the anchoring force of slope[J]. Sichuan Building Science, 2005, 31(5): 72-74. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-ACZJ200505021.htm
[3]	赵洪领, 台佳佳, 朱大勇. 边坡稳定性计算中锚固力计入方式讨论[J]. 中国水利水电科学研究院学报, 2006, 4(4): 271-276. doi: 10.3969/j.issn.1672-3031.2006.04.006 ZHAO Hong-ling, TAI Jia-jia, ZHU Da-yong. Discussion on calculation method of anchorage force in slope stability calculation[J]. Journal of China Institute of Water Resources and Hydropower Research, 2006, 4(4): 271-276. (in Chinese) doi: 10.3969/j.issn.1672-3031.2006.04.006
[4]	WANG Wen-jie, SONG Qian-qiang. Mechanical behavior of fully grouted GFRP rock bolts under the joint action of pre-tension load and blast dynamic load[J]. Tunneling and Underground Space Technology, 2018, 38(3): 82-91.
[5]	ZHANG R, ZHAO J, WANG G X. Stability analysis of anchored soil slope based on finite element limit equilibrium method[J]. Mathematical Problems in Engineering, 2016(20): 1-8.
[6]	刘华丽, 李宏伟, 朱大勇. 锚固力在边坡加固中的工程应用[J]. 岩土力学, 2015, 36(2): 509-514. https://www.cnki.com.cn/Article/CJFDTOTAL-YTLX201502032.htm LIU Hua-li, Li Hong-wei, ZHU Da-yong. A method for calculation of anchor force in slope reinforcement[J]. Rock and Soil Mechanics, 2015, 36(2): 509-514. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTLX201502032.htm
[7]	HARIYADI S M, YOSHIMI S. Experimental analysis of anchor bolt in concrete under the pull-out loading[J]. Procedia Engineering, 2017, 171: 926-933. doi: 10.1016/j.proeng.2017.01.391
[8]	LAURA B M, TIJANI M, HADJ-HASSEN F. A new analytical solution to the mechanical behaviour of fully grouted rock bolts subjected to pull-out tests[J]. Construction & Building Materials, 2011, 25(2): 749-755.
[9]	雷军, 肖世国. 考虑滑体中锚拉力扩散效应的锚固边坡稳定性分析方法[J]. 岩土工程学报, 2019, 41(9): 1724-1730. https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC201909019.htm LEI Jun, XIAO Shi-guo. Stability analysis method for anchored slopes considering stress dispersion effect in sliding mass due to anchor tensile forces[J]. Chinese Journal of Geotechnical Engineering, 2019, 41(9): 1724-1730. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC201909019.htm
[10]	卢坤林, 王运敏, 朱大勇. 三维滑体锚固力计算方法及工程应用[J]. 岩土力学, 2017, 38(2): 501-506, 516. https://www.cnki.com.cn/Article/CJFDTOTAL-YTLX201702026.htm LU Kun-lin, WANG Yun-min, ZHU Da-yong. A method for calculating reinforcing forces of 3D sliding mass and its engineering applications[J]. Rock and Soil Mechanics, 2017, 38(2): 501-506, 516. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTLX201702026.htm
[11]	LI Z, WEI J, YANG J. Stability calculation method of slope reinforced by prestressed anchor in process of excavation[J]. The Scientific World Journal, 2014(10): 1-7.
[12]	曹国金, 姜弘道, 熊红梅. 一种确定拉力型锚杆支护长度的方法[J]. 岩石力学与工程学报, 2003, 22(7): 1141-1145. https://www.cnki.com.cn/Article/CJFDTOTAL-YSLX200307018.htm CAO Guo-jin, JIANG Hong-dao, XIONG Hong-mei. Calculation method of support length for stretched bolts[J]. Chinese Journal of Rock Mechanics and Engineering, 2003, 22(7): 1141-1145. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YSLX200307018.htm
[13]	YANG G, ZHONG Z, ZHANG Y, et al. Optimal design of anchor cables for slope reinforcement based on stress and displacement fields[J]. Journal of Rock Mechanics and Geotechnical Engineering, 2015, 7(4): 411-420.
[14]	LI C, WU J, WANG J, et al. Layout and length optimization of anchor cables for reinforcing rock wedges[J]. Bulletin of Engineering Geology and the Environment, 2016, 75(4): 1399-1412.
[15]	石露, 李小春, 任伟, 等. 蚁群算法与遗传算法融合及其在边坡临界滑动面搜索中的应用[J]. 岩土力学, 2009, 30(11): 3486-3492. https://www.cnki.com.cn/Article/CJFDTOTAL-YTLX200911057.htm SHI Lu, LI Xiao-chun, REN Wei, et al. Hybrid of ant colony algorithm and genetic algorithm and its application to searching critical slope slip surface[J]. Rock and Soil Mechanics, 2009, 30(11): 3486-3492. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTLX200911057.htm
[16]	HUANG Q J, ZHANG K, SONG J C, et al. Adaptive differential evolution with a Lagrange interpolation argument algorithm[J]. Information Sciences, 2018, 9(4): 180-202.
[17]	徐东方, 郭战伟, 侯春娟. 约束尺度和算子自适应变化的差分进化算法[J]. 数学的实践与认识, 2017(16): 151-158. https://www.cnki.com.cn/Article/CJFDTOTAL-SSJS201716018.htm XU Dong-fang, GUO Zhan-wei, HOU Chun-juan. Differential evolutionary algorithm of constraint scale and operator adaptive change[J]. Mathematics in Practice and Theory, 2017(16): 151-158. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-SSJS201716018.htm
[18]	王凯. 层状弹性体系的力学分析与计算[M]. 北京: 科学出版社, 2009. WANG Kai. Mechanical Analysis and Calculation of Layered Elastic System[M]. Beijing: Science Press, 2009. (in Chinese)
[19]	卢坤林, 朱大勇, 甘文宁, 等. 一种边坡稳定性分析的三维极限平衡法及应用[J]. 岩土工程学报, 2013, 35(12): 2276-2282. https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC201312023.htm LU Kun-lin, ZHU Da-yong, GAN Wen-ning, et al. 3D limit equilibrium method for slope stability analysis and its application[J]. Chinese Journal of Geotechnical Engineering, 2013, 35(12): 2276-2282. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC201312023.htm
[20]	岩土锚杆与喷射混凝土支护工程技术规范：GB 50086—2015[S]. 2015. Technical code for engineering of ground anchoranges and shotcrete support: GB 50086—2015[S]. 2015. (in Chinese)
[21]	基于分片光滑技术三维边坡安全系数计算软件V2.0(登记号：2019SR0071874)[CP]. 3D slope safety factor calculation software based on piecewise smoothing technique V2.0 registration. No. 2019SR0071874[CP]. (in Chinese)