一种实现工程岩体力学参数反演的改进粒子群优化方法与流程

本发明属于岩体变形和受力状态分析的岩土工程问题，具体涉及一种实现工程岩体力学参数反演的改进粒子群优化方法。

背景技术：

工程岩体变形和受力状态分析的难题之一是如何恰当地估计岩土体的力学模型参数。最优估计的岩土体力学模型参数是通过比较现场观测到的信息数据与理论模型得到的。通过定义目标函数，将参数识别反问题转化为优化问题处理。伴随着计算机技术的发展，正分析的理论和计算方法逐渐成熟，观测仪器的精度也逐步提高，根据现场观测数据进行岩土力学模型参数反演具有良好的应用前景。目前粒子群算法中的适应值在算法前期改善较快，而在中后期增加缓慢，迭代次数过多，占用较多运行时间，效率较低。

技术实现要素：

发明目的：本发明的针对现有技术中存在的问题，提供一种实现工程岩体力学参数反演的改进粒子群优化方法。

技术方案：本发明所述的一种实现工程岩体力学参数反演的改进粒子群优化方法，包括下述步骤：

(1)根据边坡工程岩体的已知性能参数设置改进粒子群优化算法参数；

(2)随机初始化每个粒子的速度与位置，定义反演适应度函数，计算各粒子的适应值，根据适应值评价各粒子位置的优劣程度；

(3)进入快速粒子群搜索阶段，不断更新当前粒子的个体极值及全局极值，获取粒子个体最优位置和粒子种群最优位置；

(4)检验是否满足直接搜索引入准则，如满足，进入直接搜索阶段，否则，返回到上一步骤继续快速粒子群搜索迭代；

(5)进入直接搜索阶段，将快速粒子群搜索阶段得到的最优解p0作为直接搜索阶段的初始解，粒子种群最优适应值pbest作为直接搜索阶段的初始代最佳适应值，进行精度搜索；

(6)引入改进算法终止判断准则，如满足，迭代优化终止，此时输出的最优解p0即为反演所得最优的边坡工程力学参数；如不满足，返回上一步骤继续搜索。

上述步骤(1)中，设置的改进粒子群优化算法参数包括粒子种群大小m、惯性权重ω、加速因子c1、c2、解的维数d、每一维最大速度vimax、每一维上的解上下限ximax、ximin、初级精度ε1、高级精度ε。

步骤(2)中，适应度函数为式中，fi(x)为第i个量测点上发生的位移量的计算值，ui为相应的位移量的实测值。f(x)的值越小，粒子适应度越高，位置越好，解的品质越高，也就是说该粒子对应的一组力学参数更优。

较优的，在计算各粒子的适应值前，先判断粒子是否违反约束，对于违反约束的粒子采用隐匿墙处理；否则通过适应度函数计算各粒子的适应值。

优选的，步骤(3)中，针对较差粒子施以“加时飞行策略”，根据下式更新粒子新位置：

xid^t+1＝xid^t+ti*νid^t+1；

νid^t+1＝ωνid^t+c1r1(pid^t-xid^t)+c2r2(p^tgd-xid^t)；

式中：

ti为第i个粒子的飞行时间，fi为较差粒子的适应值，farv为粒子群体的平均适应值；

vid^t+1为每个粒子的新速度，t为当前进化代数；d＝1，2，...，n，n为粒子空间搜索的维数；i＝1，2，...，m，m为种群规模；w为惯性权重；

c1、c2为加速因子，c1调节粒子向自身最优位置飞行的步长，c2调节粒子向全局最优位置飞行的步长；

r1、r2为分布在[0，1]之间的随机数；

pid为第i个粒子d维上的个体最优位置；xid为第i个候选解中第d个待优化力学参数变量。

同时，针对最优粒子采取“保持最优位置策略”，群体中除较差粒子和最优粒子以外的粒子采用正常的位置更新公式xid^t+1+1＝xid^t+νid^t+1计算其新位置。

进一步的，本步骤中，更新粒子个体极值时，若粒子的当前适应值较其历史最好个体极值小，则对其进行更新操作，否则个体极值不变；更新当前粒子的全局极值时，比较其当前适应值与全局极值，如前者小则用其更新全局极值，且全局极值对应的参数向量x为粒子种群的当前代最优解。

上述步骤(4)中，直接搜索引入准则包括如下两种情形：

(1)在指定的迭代次数内解的质量无明显改进，即个体最优位置几乎不变，出现了解的收敛停滞现象；

(2)粒子的适应值满足时，ε1为设定的初级精度；

当满足上述两个条件之一时，算法可进入直接搜索阶段。

具体的，步骤(5)中，在最优解p0对应的最优反演参数组合的附近空间内生成m-1个新的力学参数组，利用直接搜索法进行精度搜索，range线性递减，并计算每一个力学参数组合方案的适应值，若该值比种群最优适应值pbest小，则用该值更新pbest，最终将适应值最小的一组参数更新为最优解p0；其中，range为优化参数变化范围矢量，线性递减range＝(1-α)×range，其中α为优化参数的领域搜索范围线性递减参数。

上述步骤(6)中，算法终止判断准则包括下述两种情形：

(1)最大直接搜索代数k或在指定的代数内解的质量无明显改进；

(2)粒子适应度函数值ε1为设定的初级精度；

当满足上述两个条件之一时，算法即可终止。

有益效果：与现有技术相比，本发明的优点在于：本发明对粒子群优化算法进行改进，先通过快速粒子群搜索获取粒子个体最优位置和粒子种群最优位置，针对不同位置的粒子实行不同的位置更新方式，然后进一步通过精度搜索、range线性递减得到最优解，最终反演得到最优的边坡工程力学参数；将该方法应用于工程岩体力学参数反演，可以提高优化反演效率，减少迭代次数，减少正算计算量。

附图说明

图1为实施例中边坡开挖前后的二维离散模型，其中，(a)为开挖前的二维离散模型，(b)为开挖后的二维离散模型；

图2为实施例中应用本发明的改进粒子群优化方法反演边坡工程岩体力学参数过程中粒子的迭代收敛图。

具体实施方式

下面结合附图对本发明的技术方案作进一步说明。

本发明的一种实现工程岩体力学参数反演的改进粒子群优化方法，包括下述步骤：

(1)设置改进粒子群优化算法参数；根据工程岩体的已知性能参数设置改进粒子群优化算法的参数，设置的参数包括粒子种群大小m、惯性权重ω、加速因子c1、c2、解的维数d、每一维最大速度vimax、每一维上的解上下限ximax、ximin、初级精度ε1和高级精度ε；

(2)随机初始化每个粒子的速度与位置，定义反演适应度函数，计算各粒子的适应值，通过适应值评价各粒子位置的优劣程度；

反演适应度函数为：fi(x)为第i个量测点上发生的位移量的计算值，ui为相应的位移量的实测值。f(x)的值越小代表粒子适应度越高，位置越好，解的品质越高，也就是说该粒子对应的一组力学参数更优。

本步骤中，在计算各粒子的适应值前，要判断粒子是否违反约束，粒子越出解空间的边界即判定粒子违反约束；对于违反约束的粒子采用隐匿墙处理，隐匿墙处理是指粒子允许越过物理限定区域，但是越界的粒子并不估计其适应度；否则通过适应度函数计算各粒子的适应值。

(3)进入快速粒子群搜索阶段：不断更新当前粒子的个体极值及全局极值，获取粒子个体最优位置和粒子种群最优位置；

具体而言，根据粒子的适应值判断粒子的优差，适应度数值越小，认为其对应的适应度越好，所在位置就越好，所代表的一组反演参数就越优。

针对较差粒子施以“加时飞行策略”，根据下式更新粒子新位置xid^t+1＝xid^t+ti*νid^t+1，

其中，第i个粒子的飞行时间fi为较差粒子适应值，farv为粒子群体平均适应值；

每个粒子的新速度νid^t+1＝ωνid^t+c1r1(pid^t-xid^t)+c2r2(p^tgd-xid^t)，式中，t为当前进化代数；d＝1，2，...，n，n为粒子空间搜索的维数；i＝1，2，...，m，m为种群规模；w为惯性权重；c1、c2为加速因子，c1调节粒子向自身最优位置飞行的步长，c2调节粒子向全局最优位置飞行的步长；r1、r2为分布在[0，1]之间的随机数；pid为第i个粒子d维上的个体最优位置；xid为第i个候选解中第d个待优化力学参数变量。

针对最优粒子采取“保持最优位置策略”，群体中除此以外的粒子采用正常的位置更新公式xid^t+1+1＝xid^t+νid^t+1计算其新位置；

得到每个粒子的新速度和新位置后，再次计算其适应度值，根据适应度值更新当前粒子的个体极值及全局极值；其中，更新粒子个体极值时，若粒子的当前适应值较其历史最好个体极值小，则对其进行更新操作，否则个体极值不变；更新当前粒子的全局极值时，比较当前适应值与全局极值，如前者小则用其更新全局极值，且全局极值对应的参数向量x为粒子种群的当前代最优解。

(4)检验是否满足直接搜索引入准则，如满足，进入直接搜索阶段，否则，返回到上一步骤继续快速粒子群搜索迭代；

当满足下述两条件之一时，代表满足直接搜索引入准则，算法就可以进入直接搜索阶段：第一种，在指定的迭代次数内解的质量无明显改进，即个体最优位置几乎不变，出现了解的收敛停滞现象；第二种，粒子适应度值时，ε1为设定的初级精度，是一个较小的正数。

(5)进入直接搜索阶段：将快速粒子群搜索阶段得到的最优解p0作为直接搜索阶段的初始解，粒子种群最优适应值pbest作为直接搜索阶段的初始代最佳适应值，进行精度搜索；

具体而言，在最优解p0周围一定范围range内随机生成m-1个解(即在前期优化m迭代所得最优反演参数组合的附近空间内生成m-1个新的力学参数组开始利用直接搜索法进行精度搜索，range线性递减)，并计算每一个力学参数组合方案的适应值，与种群最优适应值pbest比较起来较小，则用其值更新pbest，将适应值最小的一组参数更新为最优解p0。

每一组反演参数组合包括下述参数：粒子种群大小m、惯性权重ω、加速因子c1、c2、解的维数d、每一维最大速度vimax、每一维上的解上下限ximax、ximin、初级精度ε1、高级精度ε。range为优化参数变化范围矢量，线性递减range＝(1-α)×range，其中α为优化参数的领域搜索范围线性递减参数。

(6)引入改进算法终止判断准则，如满足，迭代优化终止，此时输出的最优解p0即为反演所得最优的边坡工程力学参数；如不满足，返回上一步骤继续搜索。

当满足下述两条件之一时，代表满足算法终止判断准则，算法可以停止：第一种，最大直接搜索代数k或在指定的代数内解的质量无明显改进；第二种，粒子适应度值

实施例

如图1，边坡分三步开挖(边坡分上下两种介质)，高为50m，上两级坡比为1﹕1，开挖高度为25m，平台为5m；第三级开挖坡比为1﹕0.5，开挖高度为20m。上层介质30m，下层介质20m。

考虑分步开挖卸荷对边坡变形的影响，图1中(a)图和(b)图分别给出了其开挖前后的离散模型，反分析的基础数据为设在各坡角处观测点1，2，3的竖向位移。本构模型采用摩尔-库仑弹塑性模型，边坡材料参数值如表1所示。

表1材料参数表

先对该边坡用flac进行正演计算，以三个测点1、2、3的开挖完成时的竖向位移(向上为正)3.90mm，12.8mm，3.72mm作为真实观测值，进行参数反演计算，待反演的参数共2个，分别为上层介质变形模量e1和下层介质变形模量e2，其参数反演范围分别为0.2～0.8gpa，4.0～8.0gpa。

以上表1中内聚力c，内摩擦角φ、泊松比μ和介质容重γ作为初始参数，采用本发明的改进粒子群优化方法对边坡工程的力学参数进行反演，反演过程中粒子的迭代收敛图如图2，反演得到的一组最优边坡工程力学参数如下表2。

表2边坡岩体变形模量反演结果

同时采用标准粒子群优化算法和线性递减惯性权重的粒子群优化算法对实施例中的边坡工程力学参数进行反演，三种粒子群优化算法的运行情况如下表3，其中，fpso-ds代表本发明的改进粒子群优化算法，pso-ldw代表标准粒子群优化算法；spso代表线性递减惯性权重的粒子群优化算法。

表3：三种粒子群优化算法的运行情况

由表3可知，与现有的标准粒子群优化算法和线性递减惯性权重粒子群优化算法相比，本发明的实现边坡工程力学参数反演的改进粒子群优化算法可以有效提高优化反演效率，减少迭代次数，减少正算计算量。