具体实施方式

为了对本发明的技术特征、目的和效果有更加清楚的理解，现对照附图详细说明本发明的具体实施方式。

本发明的实施例提供了一种基于遗传算法和离散元分析法的锚固边坡安全性评价方法。如建筑边坡工程技术规范GB50330-2013和岩土锚杆与喷射混凝土支护工程技术规范GB50086-2015所述目前岩质锚固边坡的设计和常用的评价方法如简单条分法、毕肖普法和瑞典条分法等都是基于是静态分析只是把锚索的作用等效为一个均匀分布的力，计算整个边坡的静态稳定性。而且都没有考虑到岩体的蠕变效应。而本方法有两个创新：第一，可以充分利用离散元分析法考虑到锚索和岩体的相互作用，以及岩质边坡结构面(不连续面)的作用。第二，可以考虑岩石蠕变对于岩质锚固边坡长期安全性的影响。将安全性和时间挂钩进行动态分析，充分考虑了随着岩石蠕变边坡安全性的变化。第三，充分利用实验数据和现场动态的监测数据，将遗传算法和离散元分析软件相结合，实现了自动准确的确认岩质锚固边坡分析计算模型参数。使得整个岩质锚固边坡分析计算模型的计算响应更加符合实际工程情况，进而使得评价结果更加准确。该方法的实施不是简单各自计算后仅在结果进行结合寻优，而是基于python语言编程将遗传算法的每一代计算都与离散元分析计算过程相结合，动态的结合离散元分析法进行遗传算法的寻优计算，该结合方法首次提出。

请参考图1-2，图1是本发明实施例中一种基于遗传算法和离散元分析法的锚固边坡安全性评价方法的流程图，图2是本发明实施例中基于遗传算法和离散元分析法的边坡长期安全性评价方法实现流程图，具体包括如下步骤：

S1：基于现场勘察和仪器监测收集岩质锚固边坡的基础数据。

根据锚固岩质边坡现场勘察收集滑坡的形态数据例如：锚固边坡的高程、滑坡边界、结构面的形态和分布特征、边坡的变形特征等。根据现场的监测仪器收集锚固边坡的变形数据和锚固结构的受力数据等。例如：基于GPS监测墩收集坡表位移数据和基于锚索测力计收集锚索锚固力的变化情况。

S2：基于离散元分析软件3DEC构建岩质锚固边坡分析计算模型。

基于步骤S1收集到的数据，在3DEC离散元分析软件中构建岩质锚固边坡的概化岩质锚固边坡分析计算模型。忽略次要的影响因素，考虑主要的影响因素从而提高软件的计算效率。

S3：根据岩体蠕变试验确定蠕变本构模型和蠕变参数范围，步骤S3的具体实现步骤如下：

S3-1：确定蠕变本构模型。目前岩土进行蠕变分析最常用的蠕变本构模型分为两类，分别是基于试验和数据拟合形成的经验模型和通过理论分析得到的元件模型。但是经验模型缺乏理论依据，只是试验数据的拟合。而元件模型可以更加真实地反映岩土体本身力学性质，同时简单易用，因此在实际工程中使用的更加广泛。最常被使用的原件模型如图3所示，图3中，a为Maxwell元件模型，b为Kelvin元件模型，c为Burgers元件模型，d为Poynting-Thomson元件模型；

本方法采用三轴蠕变试验和元件模型曲线类比的方法，确定蠕变本构模型的基本种类。具体过程为首先整理岩体的三轴蠕变试验数据，得到不同围压下轴向应变-时间关系曲线。然后根据曲线特点，确定模型所应该包含的元件类型，如弹性元件、粘性元件等；最后从上述模型中选择包含这些元件的模型，将其典型应变-时间曲线与试验结果进行对照，确定蠕变本构的基本类型。

S3-2：基于三轴蠕变试验数据，估算本构模型的参数范围。

Maxwell本构模型如图3中a所示，其应力应变关系式：

式中k为弹性模量，η为粘滞系数，σ为应力，ε为应变。

Kelvin本构模型如图3中b所示，其应力应变关系式：

式中k为弹性模量，η为粘滞系数，σ为应力，ε为应变。

Burgers本构如图3中c所示，其应力应变关系式：

式中k为弹性模量，其中k₁为元件模型中Kelvin体的弹性模量，k₂为元件模型中Maxwell体的弹性模量，η为粘滞系数，其中η₁为元件模型中Kelvin体的粘滞系数，η₂为元件模型中Maxwell体的粘滞系数，σ为应力，ε为应变。

Poynting-Thomson本构模型如图3中d所示，其应力应变关系式：

式中k₁为元件模型中Maxwell体的弹性模量弹性模量，k₂为元件模型中弹簧的弹性模量，η为粘滞系数，σ为应力，ε为应变。

S4：在3DEC软件中将步骤S3确定的蠕变本构模型赋值给步骤S2中构建的岩质锚固边坡分析计算模型，将步骤S2确定的蠕变参数范围作为遗传算法的输入参数，将遗传算法和岩质锚固边坡分析计算模型的响应相结合，进行联合的迭代计算来确定岩质锚固边坡分析计算模型中蠕变本构模型的最优参数；具体实现步骤如下：

S4-1：通过python编写遗传算法将3DEC软件中的锚固边坡的模型模拟计算响应和遗传算法的目标函数相结合。同时借助锚固边坡的现场监测数据确定蠕变本构的最优参数。

遗传算法实现过程：

(1)确定被优化的目标参数。本方法选择岩体蠕变本构的参数作为目标参数；

(2)对目标参数进行二进制编码；

(3)适宜度函数选择

适宜度值决定了种群中个体的优劣情况，决定了当前的一代中的优秀的个体是否会被选择继承到下一代。适宜度函数值非负，而且值越大表示该个体越优秀。

(4)遗传算子

包括选择算子、交叉算子和变异算子

选择算子，根据不同个体对环境的适应度，给出被选择或者被淘汰的概率，在根据一定的规则对个体进行选择，常用的方法有轮盘赌选择、排序选择和联赛选择。本评价方法使用轮盘赌选择。

假设种群中包含n个个体分别为x1,x2···,xn,每个个体的适应度为f(x_i)，由于个体都有一定被选择的概率，因此适应度应满足式(1)：

1≥f(x_i)≥0 (1)

定义p_i为分布概率，表示个体的生存概率，根据式(2)进行求解

分布律应满足式(3)和(4)，

p_i≥0 (3)

当个体的适应度值高时，其分布概率也高，因此优良个体被选择的可能性就好加大。同时表现低劣的个体也有一定的概率被选择，可以保证后代种群个体的多样性。

交叉算子，是指指定基因交换规则使得父代基因互相交换从而生成全新的子代，该算子有一定的概率获得优于父代的个体。常见的交叉算子有单点交叉、双点交叉、多点交叉、部分映射交叉、顺序交叉等。

本方法使用的是单点交叉，首先选择一个交叉点，将两个个体自交叉点之后的二进制编码进行互换操作，如图4所示。变异算子，变异算子可以使得遗传算法中的种群保持一种良好的多样性，仿照了生物进化中的基因变异。实现过程如下：第一步按照设定的变异概率从种群中选择某些个体，对选中个体中的随机一位进行0-1的变换，具体过程如图5所示，整个遗传算法的实现过程如图6所示。

S4-2：将步骤S1收集的工程实例监测数据，导入python中作为目标函数的一部分；

S4-3：进行遗传算法的迭代计算和岩质锚固边坡分析计算模型的分析计算，结合现场的监测数据，不断获取适应度值，最终获得最符合实际监测数据的最优的本构模型参数。

Python编写遗传算法所需参数说明：

目标函数func：目标函数。

目标维数ndim：目标函数的维度，一般为参数的数量。

群体规模size_pop：种群规模表示进行计算时需要的染色体个数，染色体越多,染色体间交流的信息也就越多，因此优化性能就越好,但是当染色体个数达到一定数量时,在增加数目对于算法的提升就很小啦,并且还会影响到遗传算法的效率。

最大迭代次数max_liter：求解最优值过程的最大迭代次数。

变异概率prob_mut：遗传算法进行变异算子操作时发生变异的概率。

参数范围lb和ub：是指参数的上下限，一般将优化问题的参数范围设置为该范围。

精准度precision：进行参数寻优时求解所得参数的精度。

S5：基于设定的岩质锚固边坡分析时间，提取岩质锚固边坡的特征变量；

步骤S5的具体实现步骤如下：

步骤S5-1：将步骤S3和步骤S4中确定的蠕变本构模型以及最优参数赋值给步骤S2中构建的岩质锚固边坡分析计算模型。

3DEC软件计算分析基本原理：

离散单元法只要满足平衡方程即可吗，对于变形体还须符合本构方程即可。

(1)物理方程

在离散元分析中相互作用的两个块体需要满足如下的物理本构方程：

ΔF_n＝K_nΔU_n

ΔF_s＝K_sΔU_s

F_n＝F_n1+K_nΔU_n

F_s＝F_s1+K_sΔU_s

式中K_n、K_s为节理面的法向和切向接触刚度，ΔF_n、ΔF_s为块体间法向作用力增量和切向作用力增量，ΔU_n、ΔU_s为块体间法向位移增量和切向位移增量，F_n、F_s为块体间法向作用力和切向作用力，F_n1、F_s1为块体间法向力和切向力的初值或上次迭代的终值。

(2)运动方程

离散单元法的基本原理式牛顿第二定律，作用在离散块体上的力会使块体产生运动，块体形心上应该满足以下方程：

a＝F_合/m

式中，F_xi、F_yi为块体受到其他块体的作用力，F_x、F_y为块体所受的水平和竖向合力。M为块体所受的合力矩，x_i、y_i为作用点坐标，x、y为块体形心的坐标，a为加速度，m为岩块质量，为角加速度，I为转动惯量。

(3)采用向前差分格式进行数值积分，求得块体沿合力方向的速度和位移以及块体的转动量。

V(t₁)＝V(t₀)+aΔt

U(t₁)＝U(t₀)+V(t₁)Δt

式中，t₀为初始时间，Δt为计算时步，U(t₁)、V(t₁)分别为块体在t₁时间的位移和速度，θ(t₁)是t₁时间的位移和速度，θ(t₁)是t₁时间的转动惯量。

3DEC软件计算分析过程如图7所示。

步骤S5-2：依据步骤S1收集到的数据，基于现场监测仪器布置情况，设置岩质锚固边坡的岩质锚固边坡分析计算模型的特征位置进行变形监测，同时监测锚索的轴力变化情况。

步骤S5-3：基于设定的岩质锚固边坡分析时间，将特征位置的监测变量和锚索的轴力值作为特征变量进行提取。

S6：根据特征变量的变化情况，进行岩质锚固边坡长期安全性评价。

S6的具体实现步骤如下：

S6-1：基于步骤S5提取的特征变量，构建位移随时间的变化曲线和锚索轴力随时间的变化曲线以及锚索轴力的最大变化幅度；

S6-2：基于位移变化曲线和锚索的轴力变化曲线以及锚索轴力的最大变化幅度，进行岩质锚固边坡的长期安全性评价，当位移变化曲线和锚索轴力变化曲线均呈现收敛时且锚索轴力的最大变化幅度小于等于10％时，该边坡评价为安全，否则为不安全。

本方法基于遗传算法和离散元分析法，在3DEC软件中构建岩质锚固边坡蠕变岩质锚固边坡分析计算模型，采用遗传算法结合现场监测数据确定蠕变本构模型的最优参数，通过python编写程序将将两种分析方法相结合，提高了评价的准确率和效率的同时，揭示了蠕变条件下岩质锚固边坡长期变形特征与演化机理，实现了对岩质锚固边坡的长期安全性的动态评价，为同类型边坡的预测预警与防治提供了基础。

本方法选取的实例为西南地区某水电工程岩质锚固边坡。该边坡位于水电站坝址上游约11.5km的江右岸，滑坡纵长约700m，宽约320～400m，面积约0.25km²，总方量约1300万m³。

本实施例中，所述步骤S1基于现场勘察、施工资料、设计资料、试验资料和仪器监测收集岩质锚固边坡的基础数据如下表1所示：

表1滑坡岩体相关参数取值

其中ρ为岩体的密度，G是岩体的剪切模量，K是岩体的体积模量，C是岩体的粘聚力，dip是岩层的倾角，dip-direction是岩层的倾向，spaceing是岩层间隙，num是岩层的数量，stiffness-normal是岩层的法向刚度，stiffness-shear是岩层的切向刚度。

锚索的相关参数如下表2所示：

表2锚索相关参数取值

其中area为锚索的横截面积，e为锚索的弹性模量，grout_stiff_为注浆体刚度，cable_yield为锚索轴向屈服强度，cable_strain_limit_为锚索轴向应变极限，dowel_stiffness_为锚索剪切变形刚度，dowel_yield_为锚索剪切屈服强度，dowel_strain_limit_为锚索剪切应变极限。

本实施例中，所述步骤S2基于离散元分析软件3DEC构建岩质锚固边坡分析计算模型的具体实现步骤如下：

基于步骤S1中收集到的数据构建数值岩质锚固边坡分析计算模型长为700m，高为400m，倾角约55°。因为需要开展锚固岩质边坡剖面的二维计算分析，因此宽度即Y方向仅取5m。节理模型取倾角为70°倾向270°锚索长度取50m间隔20m和安装角度约35°。在3DEC中计算所需要的其余计算参数按表3和表4所示的值赋予岩质锚固边坡分析计算模型，建立的数值岩质锚固边坡分析计算模型如图8所示。

本实施例中，所述步骤S3根据岩体蠕变试验确定蠕变本构模型和蠕变参数范围的具体实现步骤如下：

S2-1：根据开展的岩石样品蠕变试验，得到的试验数据，如图9是围压为200kPa时得到的试验曲线。曲线呈现出较2个较为明显的特点：1.在荷载施加后的较短时间内，试样会产生明显的瞬时弹性变形，所以该模型应该包含弹性元件；2.随着试验进行，试样呈现了明显的减速蠕变现象，经过一定时间的减速蠕变，蠕变速率保持不变，呈现出粘弹性的特征，所以该模型应包含粘性元件。经过与常用的蠕变本构模型对比分析，发现Burgers模型的典型应变-时间曲线(图10)和该曲线最相似。因此初步确定该岩体的本构模型定为Burgers模型。由于在进行长期安全性评价中需要考虑到岩质边坡的破坏情况，因此串联一个M-C元件，十使得其复合摩尔库伦破坏准则。进而该模型变为Burgers-mohr模型(图11)。

S2-2：Burgers-mohr模型共包含6个参数，c、E_M、η_M、E_K、η_K，分别是粘聚力、内摩擦角、弹性模量、马克斯韦尔粘性系数、粘弹性模量和开尔文粘性系数。C和为常规强度参数已经通过试验确定。因此还需确定四个参数。则Burgers-mohr模型的应力-时间关系可简化为：

式中，σ为应力；ε为应变。

依据试验曲线(图10)，可知模型的初始应变ε0和稳定蠕变时的斜率m。当t＝0，基于式(5)可得ε0＝σ/EM；当t取较大数值时，的值趋近于0，所以m＝1/ηM；基于最小二乘法就可以确定参数EM、ηM。最终岩体本构模型的拟合参数如表3所示。

表3岩体蠕变参数

依据表3确定的取值范围如下表4所示：

表4岩体本构模型参数范围

由图12可知随着遗传算法的进行岩质锚固边坡分析计算模型的坡表位移响应与现场监测数据越来越接近，而且55代左右目标函数基本趋于稳定，收敛很快，参数值很快逼近最优值。

计算结果显示当E_M＝19958.7661MPa、η_M＝1.35E+07MPa·h、E_K＝41200.39718MPa、η_K＝254000.3466MPa·h时，岩质锚固边坡分析计算模型特征位置特征变量的响应与监测值接近，选取该组参数作为岩质边坡本构模型的参数。

本实例中，如所述步骤S5基于设定的岩质锚固边坡分析时间，提取岩质锚固边坡的特征变量的具体实现步骤如下：

将步骤S3中确定的蠕变本构模型和步骤S4中确定的最优参数带入步骤S2中构建的岩质锚固边坡分析计算模型。输入设定的蠕变分析时间，本案例设置120天为如蠕变分析时间，提取了坡表4个特征位置的水平方向的位移和锚索荷载随时间的变化数据。坡表特征位置和锚索位置如图13所示，提取的特征变量和锚索的数据如图14和图15所示。

本实施例中，所述步骤S6根据特征变量的变化情况，进行岩质锚固边坡长期安全性评价的具体实现步骤如下：

由图14和图15可知坡表特征点X方向位移随时间呈现出一开始波动随后趋于稳定的趋势。锚索轴力一开始由于刚安装锚索，随着坡表变形锚索轴力呈现出急速增加，随后由于锚索和岩体的相互作用，随着时间锚索轴力呈现出增速减慢随后下降的现象，最终锚索轴力呈现出稳定的趋势。而且趋于稳定后锚索轴力的变化率小于10％，基于《岩土锚固工程的长期性能与安全评价》一文所述的判断阈值即锚索轴力的变化幅度小于10％锚固结构判断为安全。结合特征变量以及锚索轴力的变化趋势分析，判断该锚固边坡的长期安全性为：安全。

从图12中可以看出，模拟计算模型响应的准确性随着遗传算法的迭代不断增加最终趋于稳定，这表明遗传算法可以依据实际监测情况快速准确的确定岩体蠕变本构的参数。而特征变量随时间的变化如图14所示，这表明岩质锚固边坡的长期安全性是一个动态的变化过程，而目前常规评价方法均是静态的，是不准确的。而本方法所提出的方法既考虑考虑到锚固结构与岩体的相互作用，也考虑到岩质边坡岩土体自身特性(蠕变、结构面等)的影响，将岩质锚固边坡的安全性与时间挂钩，实现了对岩质锚固边坡安全性的动态评价，同时还可以得到岩质锚固边坡的在指定分析时间后的位移分布情况和锚索的应力分布情况如图16所示。

遗传算法特点：

(1)自组织、自适应和智能性。在遗传算法设计中不需要提前描述出问题的全部特点并说明针对问题的各个不同特点算法应采取的措施。所以，遗传算法的鲁棒性很强，能够解决复杂的非结构化问题。

(2)在搜索的过程中使用的是基于目标函数值的评价信息，不需要满足优化函数的连续性及可导性约束的要求。

(3)易于并行化，可降低由于使用超强计算机硬件所带来的昂贵费用。

(4)算法的基本思想非常简单，并且有着极其规范的算法运行方式及实现步骤，能够直接用以具体的使用。

(5)遗传算法同时使用多个搜索点的搜索信息。遗传算法从由很多个体所组成的一个初始群体开始最优解的搜索过程，而不是从单一的个体开始搜索。是对这个群体进行选择、交叉、变异等运算，产生出的乃是新一代的群体。这些信息可以避免搜索一些不必要的搜索点，所以实际上相当于搜索了更多的点，这是遗传算法所持有的一种隐含并行性。

(6)遗传算法使用了概率搜索技术。遗传算法属于一种自适应搜索技术，其选择、交叉、变异等运算都是以一种概率的方式进行，从而增加了其搜索过程的灵活性。

目前离散元分析中本构模型的参数确认方式主要是通过室内试验，赋予一个确定的值。在结合一现场的监测数据进行反演分析，主要是通过人工试算看岩质锚固边坡分析计算模型响应是否符合现场监测和变形迹象等，效率低下且准确率低。本方法将遗传算法和离散元分析法相结合实现了岩质锚固边坡分析计算模型本构参数的自动寻优，提高工作效率的同时还增加了准确性。充分利用了遗传算法以下的优点：

(1)自组织、自适应和智能性。在遗传算法设计中不需要提前描述出问题的全部特点并说明针对问题的各个不同特点算法应采取的措施。所以，遗传算法的鲁棒性很强，能够解决复杂的非结构化问题。

(2)在搜索的过程中使用的是基于目标函数值的评价信息，不需要满足优化函数的连续性及可导性约束的要求。

(3)易于并行化，可降低由于使用超强计算机硬件所带来的昂贵费用。

(4)遗传算法同时使用多个搜索点的搜索信息。遗传算法从由很多个体所组成的一个初始群体开始最优解的搜索过程，而不是从单一的个体开始搜索。是对这个群体进行选择、交叉、变异等运算，产生出的乃是新一代的群体。这些信息可以避免搜索一些不必要的搜索点，所以实际上相当于搜索了更多的点，这是遗传算法所持有的一种隐含并行性。

(5)遗传算法使用了概率搜索技术。遗传算法属于一种自适应搜索技术，其选择、交叉、变异等运算都是以一种概率的方式进行，从而增加了其搜索过程的灵活性。

本方法使用python语言编程通过3DEC软件自带的ipython平台将遗传算法和离散元分析法相结合。本方法提出的评价方法不是简单的将两种评价方法各自独立计算后仅在结果上进行结合。而是将每一代计算的目标函数都动态的从离散元分析中提取分析结果，每一代新的染色体的生成都动态的把新的参数值赋给离散元进行动态计算如图2所示。最终在遗传算法和离散元分析法的迭代计算后确定最符合现场监测数据的本构模型参数。目前该方法还是首次应用。

本发明的有益效果是：可以有效提高岩质锚固边坡安全性评价的准确性和评价效率，该方法可以为岩质锚固边坡滑坡的演化机理研究与动态安全性评价提供有利条件。

以上所述仅为本发明的较佳实施例，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。