CN113431056B - 基于粒子群算法和离散元分析的锚固边坡安全性评价方法 - Google Patents

Abstract

本发明提供了一种基于粒子群算法和离散元分析的锚固边坡安全性评价方法，基于现场勘察和仪器监测收集岩质锚固边坡的基础数据；基于离散元分析软件3DEC构建岩质锚固边坡分析计算模型；根据岩体蠕变试验确定蠕变本构模型及蠕变参数范围；基于粒子群算法和岩质锚固边坡分析计算模型的响应确定岩质锚固边坡分析计算模型中蠕变本构模型的最优参数；基于设定的分析时间，提取岩质锚固边坡的特征变量。根据特征变量的变化情况进行岩质锚固边坡长期安全性评价。本发明的有益效果是：有效提高了岩质锚固边坡安全性评价的准确性和评价效率，为岩质锚固边坡滑坡的演化机理研究与动态安全性评价提供有利条件。

Description

基于粒子群算法和离散元分析的锚固边坡安全性评价方法

技术领域

本发明涉及地质灾害防治领域，具体涉及水电站岩质锚固边坡的动态长期安全性的评价方法，尤其涉及一种基于粒子群算法和离散元分析的锚固边坡安全性评价方法。

背景技术

自一系列国家水电发展规划实施以来，我国西南地区兴建了一大批水利水电工程，随着这些工程的建设和运行，开展了大量的岩质高边坡治理工程，锚固结构作为边坡治理工程中最常用的手段之一被大量应用。因此形成了大量的岩质锚固边坡，这些边坡的长期安全性不仅关系着整个水电工程的运行，还严重关系着库区人民的生命和财产安全。在目前评价滑坡稳定性最常用的规范中仅仅是将锚固结构等效为均匀分布的一个力进行平衡计算，难以考虑到锚固结构和边坡的动态相互作用，同时目前规范对于边坡稳定性评价是一种静态的评价方法，对岩质锚固边坡的长期安全状态难以作出准确判断。

离散单元法(Discrete element method,DEM)是一种对存在结构面和节理的岩体进行应力分析的数值方法，在1971年Cundall率先提出。可模拟非连续岩体的大位移，接触面滑移、分离等现象,同时可以更加真实的反应出结构面和节理岩体内部变形以及应力分布情况。该方法以牛顿第二定律作为理论基础。可运用于刚体(RigidBlock)和可变形体(Deformableblock)的分析研究。3DEC是三维离散元程序(3Dimension distinct elementCode)的简称，是由美国ITASCA公司开发研究的成熟离散元数值模拟软件。不仅可以模拟存在不连续情况的工程案例，还可以模拟岩体在动态或静态荷载下的受力和变形情况，主要适用于大变形和转动的分析计算,不但可以自动判别块体间各接触点还能分析动态或静态荷载下介质受力以及变形情况。目前应用在了边坡评价、地下工程和层理、节理或断层等结构对岩石地基的影响等领域。然而，目前大多数分析均采用Mohr模型，没能考虑到岩体材料蠕变特性导致的其临界应变应变状态改变和强度劣化；同时对于蠕变参数的选择也是通过试算和试验导致确定参数效率低下，岩质锚固边坡分析计算模型的响应与实际监测数据还存在一定的差异。这两个不足导致了目前对岩质锚固边坡的演化机理和长期安全性评价等研究存在偏差。

发明内容

为了解决上述问题，本发明提供了一种基于粒子群算法和离散元分析的锚固边坡安全性评价方法，借助于3DEC软件实现，主要包括以下步骤：

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

S2：根据所述基础数据，利用离散元分析软件3DEC构建岩质锚固边坡分析计算模型；

S3：根据岩体蠕变试验确定蠕变本构模型和蠕变参数范围；

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

S5：将步骤S3中得到的蠕变本构模型和步骤S4中得到的岩质锚固边坡分析计算模型中蠕变本构模型的最优参数带入步骤S2中构建的岩质锚固边坡分析计算模型，基于设定的岩质锚固边坡分析时间，提取岩质锚固边坡的特征变量；

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

进一步地，所述基础数据包括由现场勘察得到的形态数据以及仪器监测得到的锚固边坡的变形数据和锚固结构的受力数据，所述形态数据包括锚固边坡的高程、滑坡边界、结构面的形态和分布特征和边坡的变形特征。

进一步地，根据岩体蠕变试验确定蠕变本构模型和蠕变参数范围的过程如下：

首先，整理岩体的三轴蠕变试验数据，得到不同围压下轴向应变-时间关系曲线，然后根据该轴向应变-时间关系曲线特点，确定元件模型所包含的元件类型，最后选择包含所述元件类型的元件模型，将该元件模型的典型应变-时间曲线与三轴蠕变试验的试验结果进行对照，确定蠕变本构模型；

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

进一步地，确定岩质锚固边坡分析计算模型中蠕变本构模型的元件模型的最优参数的过程如下：

S4-1：通过python编写粒子群算法将3DEC软件中的岩质锚固边坡分析计算模型响应和粒子群算法的目标函数相结合；

S4-2：将步骤S1收集的仪器监测数据，导入步骤S4-1所述的python中，将3DEC软件中的多个坡表监测点在X方向的位移及其对应的实际监测数据的差值求和作为目标函数；

S4-3：借助3DEC软件内置的ipython编写环境，导入粒子群算法，设置对应参数，该对应参数包括：目标函数func、目标维数ndim、群体规模size_pop、最大迭代次数max_liter、参数下限lb、参数上限ub、惯性权重、个体记忆和群体记忆c1、c2，然后进行粒子群算法和岩质锚固边坡分析计算模型的分析计算的联合迭代计算，结合现场监测数据，不断获取适应度值，最终获得最符合实际监测数据的蠕变本构模型的最优参数。

进一步地，提取岩质锚固边坡的特征变量的过程如下：

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

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

S5-3：向3DEC软件内输入设定的岩质锚固边坡分析时间，将特征位置的监测变量和锚索的轴力值作为特征变量进行提取。

进一步地，进行岩质锚固边坡长期安全性评价的过程如下：

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

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

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

附图说明

下面将结合附图及实施例对本发明作进一步说明，附图中：

图1是本发明实施例中一种基于粒子群算法和离散元分析的锚固边坡安全性评价方法的流程图。

图2是本发明实施例中基于粒子群算法和离散元分析法的锚固边坡长期安全性评价方法的实现流程图。

图3是本发明实施例中常用的元件模型示意图。

图4是本发明实施例中粒子位置更新方式的示意图。

图5是本发明实施例中粒子群算法实现流程图。

图6是本发明实施例中3DEC软件计算分析流程图。

图7是本发明实施例中3DEC岩质锚固边坡分析计算模型。

图8是本发明实施例中围压为200kpa时岩体的应变-时间曲线。

图9是本发明实施例中Burgers模型的典型应变-时间曲线。

图10是本发明实施例中Burgers-mohr模型。

图11是本发明实施例中粒子群算法目标函数值随迭代次数的变化图。

图12是本发明实施例中坡表监测位置及锚索位置图。

图13是本发明实施例中坡表特征点X方向位移随时间变化图。

图14是本发明实施例中锚索轴力随时间变化图。

图15是本发明实施例中岩质锚固边坡位移场及锚索应力分布图。

具体实施方式

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

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

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

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

根据锚固岩质边坡现场勘察收集滑坡的形态数据例如：锚固边坡的高程、滑坡边界、结构面的形态和分布特征、边坡的变形特征等。根据现场的监测仪器收集锚固边坡的变形数据和锚固结构的受力数据等。例如：基于GPS监测墩收集坡表位移数据和基于锚索测力计收集锚索锚固力的变化情况。根据设计和施工资料确认锚索的安装角度、长度、注浆体的强度等。根据试验资料确认岩质锚固边坡的力学参数和蠕变试验资料等。

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

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

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

S3-1：确定蠕变本构模型。目前岩土进行蠕变分析最常用的蠕变本构模型分为两类，分别是基于试验和数据拟合形成的经验模型和通过理论分析得到的元件模型。但是经验模型缺乏理论依据，只是试验数据的拟合。而元件模型可以更加真实地反映岩土体本身力学性质，同时简单易用，因此在实际工程中使用的更加广泛。Kelvin模型、Maxwell模型、Poynting-Thomson模型、Burgers模型等是最常被使用的原件模型如图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软件中的岩质锚固边坡分析计算模型的响应和粒子群算法的目标函数相结合。同时借助锚固边坡的现场监测数据确定蠕变本构模型的最优参数。

粒子群算法是将所有的粒子在一个D维空间进行搜索，当前粒子位置的优劣是通过适应度函数进行判断的，而粒子群的适应度值是由目标函数决定的。粒子根据适应值以及记忆功能，可以记住粒子在迭代更新的过程中的最优解，以及所有粒子的群体最优解，在粒子进行位置更新时会受到当前的位置和速度的影响，同时也会根据个体最优和全局最优的情况进行调整种群粒子的更新位置和速度,从而保证了种群的整体最佳，找到全局的最优解。

如图4所示，粒子群优化算法求最优解的原理如下：

D维空间中，有N个粒子，粒子i的位置向量如下：

x_i＝(x_i1,x_i2,...,x_iD),i＝1,2,...n

粒子i的速度表示如下：

v_i＝(v_i1,v_i2,...,v_iD)

粒子i个体经历过的最好的位置即个体的最优值pbest表示如下：

p_i＝(p_i1,p_i2,...,p_iD)

全局最优值为gbest也就是种群历史最好的位置如下：

p_g＝(p_g1,p_g2,...,p_gD)

粒子位置的变化限定在[X_min,d,X_max,d]。

粒子i的第d维速度更新公式为：

v_id(t+1)＝v_id(t)+c₁r₁(t)(p_id(t)-x_id(t))+c₂r₂(t)(p_gi(t)-x_id(t))

粒子i的第d维位置更新公式为：

x_id(t+1)＝x_id(t)+v_id(t+1)

式中：d—粒子维数

i—第i个粒子

c₁—调节粒子飞向个体最优值的加速常量

c₂—调节粒子飞向全局最优值的加速常量

r₁—调节粒子与个体最优值之间关系的随机函数，取r₁∈[0,1]

r₂—调节粒子与全局最优值之间关系的随机函数，取r₂∈[0,1]

粒子位置的更新方式如下图4所示，粒子群算法实现过程：

(1)对粒子群进行初始化包括随机位置和速度。本方法选择岩体蠕变本构模型的参数作为目标参数；

(2)确定适应度函数，计算适应度函数值；

(3)确定个体历史最佳位置,在计算完粒子的适应度函数值后,进行比较更新后确定最新的个体最优值；

(4)确定全局最佳位置，得到粒子的适应度函数值之后，比较当前函数值与全局最佳位置从而更新群体的最优值；

(5)依据公式更新种群中粒子的速度与位置；

(6)进行循环迭代，如果没有达到最优结果，回到第二步再次进行计算,多次连续循环上述过程直到满足结束条件。

粒子群算法的实现过程如图5所示，为了便于计算使用本方法依据粒子群算法原理，使用python语言编写粒子群算法计算程序。

Python编写粒子群算法参数说明：

目标函数func：目标函数

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

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

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

参数范围lb和ub：是指参数的上下限，一般将优化问题的参数范围设置为粒子的取值范围，一般在每一维中，粒子的变化不可以超过设定的边界值，并且不同维数的粒子的范围也可以进行不同的设置的，粒子范围根据具体的问题设置不同的值。

惯性权重：通过调节可以限制粒子的搜寻范围，使得粒子有保持运动的惯性和空间搜索能力。促进全局和局部的平衡，在新区域搜索。设置较大值便于搜索领域，较小值便于局部搜索。合适的惯性权重ω将影响粒子群算法的搜索速度和精度。

个体记忆和群体记忆c1、c2：反映了两个方面的学习：一是自我的历史经验，二是群体中的最优粒子。该参数反映了粒子群的信息交流，在调整参数方面，不能设置较大或较小，较大的c1会发散,较大的c2过早收敛一般来说c1和c2取固定常数，并且取c1和c2相等的；在众多的资料中c1和c2取值范围为[0,4]。

S4-2：将步骤S1收集的工程实例监测数据，导入python中，将3DEC模型的多个坡表监测点X方向的位移与其对应的实际监测数据的差值求和作为目标函数。

S4-3：借助3DEC软件内置的ipython编写环境，导入粒子群算法，设置对应的参数。进行粒子群算法和岩质锚固边坡分析计算模型的分析计算的联合迭代计算，不断获取适应度值，最终获得最符合实际监测数据的最优的本构模型参数。如图2所示。

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软件计算分析过程如图6所示。

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

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

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

步骤S6-1：基于步骤5提取的特征变量，构建位移随时间的变化曲线，以及锚索轴力随时间的变化曲线。

步骤S6-2：观察位移变化曲线是否呈现收敛趋势，锚索的轴力变化曲线是否收敛。计算锚索轴力的最大变化幅度。进行岩质锚固边坡的长期安全性评价。

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

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

本实施例中，所述步骤1基于现场勘察、施工资料、设计资料、试验资料和仪器监测收集岩质锚固边坡的基础数据如下表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_为锚索剪切应变极限。

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

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

建立的岩质锚固边坡分析计算模型如下图7所示，本实施例中，所述步骤3根据岩体蠕变试验确定蠕变本构模型、确定蠕变参数范围的具体实现步骤如下：

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

步骤S2-2：Burgers-mohr模型共包含6个参数，c、

E_M、η_M、E_K、η_K，分别是粘聚力、内摩擦角、弹性模量、马克斯韦尔粘性系数、粘弹性模量和开尔文粘性系数。C和

为常规强度参数已经通过试验确定。因此还需确定四个参数。则Burgers-mohr模型的应力-时间关系可简化为：

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

依据试验曲线(图8)，可知模型的初始应变ε₀和稳定蠕变时的斜率m。当t＝0，基于式(5)可得ε₀＝σ/E_M；当t取较大数值时，

的值趋近于0，所以m＝1/η_M；基于最小二乘法就可以确定参数E_M、η_M。最终岩体本构模型的拟合参数如表3所示。

表3岩体蠕变参数

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

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

本实施例中，所述步骤S4基于粒子群算法和岩质锚固边坡分析计算模型的响应以及现场的实际监测数据确定岩质锚固边坡最优参数的具体实现为：

从表4中可知岩体蠕变参数的取值范围，借助3DEC软件自带的ipyhon编程工具，将其带入使用python编写的粒子群算法中，该计算参数在生成初始粒子群时使用。将3DEC模型的多个坡表监测点X方向的位移与其对应的实际监测数据的差值求和作为目标函数。因此粒子群算法的参数情况如下表5所示，带入参数后结合步骤2构建的3DEC分析模型进行粒子群算法和3DEC离散元算法的联合计算。目标函数值随迭代次数的变化情况如下图11所示。

表5工程实例粒子群算法参数表

如图11可知随着粒子群算法的迭代目标函数值越来越小，说明随着粒子群算法的迭代3DEC岩质锚固边坡分析计算模型的响应与真实监测数据越来越接近，即蠕变参数在逐步逼近最优值。而且在40代左右目标函数值趋于稳定，收敛很快，在40代左右就达到最优值。

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

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

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

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

由图13可知，特征变量随着时间增加从一开始的急速增加。到经历几个波动后最终趋于稳定，锚索轴向荷载的变化规律是从一开始由于锚索刚安装初步发挥锚固作用的急速增加，到增速放缓，经历波动后最终趋于稳定。稳定后锚索的荷载没有发生突变，波动的范围也小于10％。基于《岩土锚固工程的长期性能与安全评价》所提出的判断阈值，锚索荷载波动范围小于10％可判断为安全，结合特征变量的变化趋势分析判断该锚固边坡的长期安全性为：安全。

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

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

(1)粒子群算法依据粒子的速度完成搜索，具有记忆能力，收敛速度快参数设置少。

(2)粒子群算法原理比较简单，构造也比较简单，无需交叉和变异，具有记忆性，易操作实现。

(3)粒子群算法是有良好的特质，粒子是根据历史上最优的粒子和全局最优的粒子位置和速度进行自我更新，搜索效率高，展现了非常好的效果。

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

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

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

Claims (4)

1.一种基于粒子群算法和离散元分析的锚固边坡安全性评价方法，借助于3DEC软件实现，其特征在于：包括以下步骤：

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

S2：根据所述基础数据，利用离散元分析软件3DEC构建岩质锚固边坡分析计算模型；

S3：根据岩体蠕变试验确定蠕变本构模型和蠕变参数范围；

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

S5：将步骤S3中得到的蠕变本构模型和步骤S4中得到的岩质锚固边坡分析计算模型中蠕变本构模型的最优参数带入步骤S2中构建的岩质锚固边坡分析计算模型，基于设定的岩质锚固边坡分析时间，提取岩质锚固边坡的特征变量；

提取岩质锚固边坡的特征变量的过程如下：

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

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

S5-3：向3DEC软件内输入设定的岩质锚固边坡分析时间，将特征位置的监测变量和锚索的轴力值作为特征变量进行提取；

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

进行岩质锚固边坡长期安全性评价的过程如下：

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

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

2.如权利要求1所述的一种基于粒子群算法和离散元分析的锚固边坡安全性评价方法，其特征在于：步骤S1中，所述基础数据包括由现场勘察得到的形态数据以及仪器监测得到的锚固边坡的变形数据和锚固结构的受力数据，所述形态数据包括锚固边坡的高程、滑坡边界、结构面的形态和分布特征和边坡的变形特征。

3.如权利要求1所述的一种基于粒子群算法和离散元分析的锚固边坡安全性评价方法，其特征在于：步骤S3中，根据岩体蠕变试验确定蠕变本构模型和蠕变参数范围的过程如下：

首先，整理岩体的三轴蠕变试验数据，得到不同围压下轴向应变-时间关系曲线，然后根据该轴向应变-时间关系曲线特点，确定元件模型所包含的元件类型，最后选择包含所述元件类型的元件模型，将该元件模型的典型应变-时间曲线与三轴蠕变试验的试验结果进行对照，确定蠕变本构模型；

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

4.如权利要求1所述的一种基于粒子群算法和离散元分析的锚固边坡安全性评价方法，其特征在于：步骤S4中，确定岩质锚固边坡分析计算模型中蠕变本构模型的最优参数的过程如下：

S4-1：通过python编写粒子群算法将3DEC软件中的岩质锚固边坡分析计算模型响应和粒子群算法的目标函数相结合；

S4-2：将步骤S1收集的仪器监测数据，导入步骤S4-1所述的python中，将3DEC软件中的多个坡表监测点在X方向的位移及其对应的实际监测数据的差值求和作为目标函数；

S4-3：借助3DEC软件内置的ipython编写环境，导入粒子群算法，设置对应参数，该对应参数包括：目标函数func、目标维数ndim、群体规模size_pop、最大迭代次数max_liter、参数下限lb、参数上限ub、惯性权重、个体记忆和群体记忆c1、c2，然后进行粒子群算法和岩质锚固边坡分析计算模型的分析计算的联合迭代计算，结合现场监测数据，不断获取适应度值，最终获得最符合实际监测数据的蠕变本构模型的最优参数。

Images