CN114117593B - 岩体力学参数联合求解方法及装置 - Google Patents

Abstract

本发明提供基于自适应粒子群算法的岩体力学参数联合求解方法及装置，能够实现对多区域岩体参数的准确获取以及对工程运行状态的及时评估。方法包括：S1、获取监测数据；S2、自适应粒子群算法求解：基于有限元数据和获取的监测数据构建变形模型；初始化粒子群算法；将变形模型作为目标函数，将每个粒子的位置信息输入目标函数可获得适应度值，其值越小代表与岩体实际变形、渗流的接近程度越高，岩体参数越接近真实值；根据粒子在搜索空间中的位置，按照相对距离划分子种群，然后确定子种群个数；进行粒子群算法的迭代计算；S3、岩体参数的获取：取出整个种群中具有最小适应度值的粒子，该粒子的位置即为求解的所有岩体力学参数。

Description

岩体力学参数联合求解方法及装置

技术领域

本发明属于岩土工程技术领域，具体涉及一种基于改进的自适应粒子群算法的岩体力学参数联合求解方法及装置。

背景技术

岩体的渗流、变形和受力状态是岩土工程分析的重点内容，准确的岩体力学参数是进行物理试验和数值模拟分析的基础。由于岩体的非均匀性，室内试验或现场测试获得的参数具有较大的离散性和不确定性，同时工程在正常运行状态下无法进行物理试验确定内部岩体的参数，更无法捕获参数随时间变化的信息。根据岩体变形或渗流的监测数据进行反分析，挖掘其中的岩体参数信息，可以有效解决该问题。

当前在岩土工程中使用的反分析算法，一般需要反复手动调整参数，此外，由于本构模型的非线性程度高以及待求解参数的数量较多，算法易陷入局部最优解，一般仅能求解少量参数，不能考虑工程多区域的变形协调。这些反分析算法在过多人为介入的情况下搜索效果不佳，不仅费时费力，而且采用这些参数进行物理试验和数值模拟计算的结果与监测数据相差较大，无法准确获取工程应力、变形和渗流情况，更无法进行工程预测。

发明内容

本发明是为了解决上述问题而进行的，目的在于提供基于改进的自适应粒子群算法的岩体力学参数联合求解方法及装置，对于高非线性、高维度问题具有求解精度高、运行速度快、性能稳定的特点，能够实现对多区域岩体参数的准确获取以及对工程运行状态的及时评估。

本发明为了实现上述目的，采用了以下方案：

<方法>

本发明提供的基于改进的自适应粒子群算法的岩体力学参数联合求解方法，其特征在于，包括以下步骤：

S1、数据获取：

S1.1、选择岩土工程中一块或多块研究区域，获取该区域的变形或/和渗流监测数据；

S1.2、根据岩土本构模型或/和渗流模型选取待求解的弹塑性参数或/和渗透参数；

S2、自适应粒子群算法求解：

S2.1、基于有限元数据和步骤S1.1获取的监测数据构建变形模型：

式中，x₁，x₂，...，x_p为待反演的岩体力学参数；m为监测断面的个数；n为断面上监测点的个数；u_ij为断面上第i个监测断面上第j个监测点的有限元计算值；u_ij ^*为相应的实测值；w_i为第i个监测断面的权重系数；

S2.2、初始化粒子群算法：设置初始化参数，包括目标函数、维度D、种群大小N、最大迭代次数T_max、拓扑结构种类、子种群最大数量Z_max、子种群最小粒子数n、拓扑结构重组间隔T、搜索范围[X_min，X_max]、位置调整机制、粒子速度范围[-V_max，V_max]、速度调整机制、每个粒子的位置、每个粒子的速度、惯性权重初始值w_start、惯性权重终止值w_end、第一加速系数初始值c_1，start、第一加速系数终止值c_1，end、第二加速系数初始值c_2，start、第二加速系数终止值c_2，end、第三加速系数初始值c_3，start、第三加速系数终止值c_3，end；

将S2.1构建的变形模型作为目标函数，用于度量计算值与实际值的接近程度，将每个粒子的位置信息输入目标函数可获得适应度值，其值越小代表与岩体实际变形、渗流的接近程度越高，岩体参数越接近真实值；维度D即待求解的所有岩体力学参数的个数；种群大小N指粒子的个数，一个粒子在搜索空间中的位置代表一组岩体参数解；搜索范围[X_min，X_max]由具体的待求解参数确定；

S2.3、根据粒子在搜索空间中的位置，按照相对距离划分子种群，根据Calinski-Harabasz指标确定子种群个数；

S2.4、粒子群算法的迭代计算：根据粒子在搜索空间中的位置及目标函数计算每个粒子的适应度值，取一次迭代步中整个粒子群最小适应度值所对应的粒子位置作为全局最优位置，得到一组最优岩体参数；迭代次数加1，更新每个粒子的速度和位置，更新算法参数w、c₁、c₂和c₃，更新全局最优位置，判断迭代次数t是否为拓扑结构重组间隔T的整数倍，若是，则重新划分子种群，否则继续运算。随着迭代的进行，整个粒子群获得的最小适应度值呈不断下降的趋势，也即愈发接近工程实际，当整个种群的最小适应度值达到设定的精度要求，或迭代次数达到最大迭代次数T_max，算法停止；

在S2.4中的粒子速度按照下列公式(2)计算：

式中，

表示第i个粒子第d维度在第t次迭代的位置分量；

表示第i个粒子第d维度在第t次迭代的速度分量；w^t表示第t次迭代时的惯性权重；

表示第t次迭代中第i个粒子的第一加速系数，

表示第t次迭代中第i个粒子的第二加速系数，

表示第t次迭代中第i个粒子的第三加速系数，

表示第t次迭代中第i个粒子的个体最小适应度值；

表示第t次迭代中第z个子种群的子种群最小适应度值；

表示第t次迭代中的全局最小适应度值，r₁、r₂、r₃为[0，1]范围内的随机数；

S3、岩体参数的获取：

在S2.4算法停止后，取出整个种群中具有最小适应度值的粒子，该粒子在搜索空间中的位置即为求解的所有岩体力学参数，该组参数后续用于工程变形、渗流试验及模拟计算，控制工程运行状态。S3中获取的结果包含所有研究区域的所有待求解岩体参数(弹塑性参数和渗透参数)，单个粒子包含的待求解参数为所有参数，即最终求解结果将由1个粒子的位置表示。

优选地，本发明提供的基于改进的自适应粒子群算法的岩体力学参数联合求解方法，还可以具有以下特征：在S1.1中，监测数据包含变形数据和渗流数据中的至少一种；

在S1.2中，根据岩土本构模型或/和渗流模型的至少一种选取待求解的弹塑性参数或/和渗透参数；岩土本构模型包括弹性模型、弹塑性模型、粘弹塑性模型。

优选地，本发明提供的基于改进的自适应粒子群算法的岩体力学参数联合求解方法，还可以具有以下特征：在S2.1中，有限元数据为通过使用多组不同岩体力学参数的有限元计算获得的与监测点位置对应的变形、渗流数据。

优选地，本发明提供的基于改进的自适应粒子群算法的岩体力学参数联合求解方法，还可以具有以下特征：在S2.2中，粒子的速度和位置受约束，数值超过限定的范围将触发调整机制，具体包括：

粒子速度范围限制在[-V_max，V_max]，当速度过大将触发调整机制，包括2种方式：①将速度调整为上一迭代步的速度；②速度方向反转并使数值减小一半；

粒子位置范围限制在[X_min，X_max]，当位置超出搜索范围将触发调整机制，包括5种方式：①粒子位置调整至边界；②在搜索区域随机生成；③收缩粒子速度10倍，使其落入搜索范围内；④根据边界镜像反射，使位置落入搜索范围内；⑤将粒子位置定义到上一迭代步位置与边界的中点；

在初始化时应选择1种速度调整机制和1种位置调整机制。

优选地，本发明提供的基于改进的自适应粒子群算法的岩体力学参数联合求解方法，还可以具有以下特征：在S2.2中，拓扑结构种类包含金字塔形结构、环形结构、星形结构和冯诺依曼结构，子种群的拓扑结构与全局拓扑结构一致，在初始化时需选择其中1种，具体包括：金字塔形结构是一种多层结构，顶层粒子少，底层粒子多，表现为粒子的速度迭代仅受上、下层粒子影响；环形结构表现为在算法开始前对所有粒子进行编号，粒子的速度迭代仅受编号在其前、后的两个粒子影响；星形结构即每个粒子的速度迭代受全局其他所有粒子影响；冯诺依曼结构与环形结构类似，但影响速度迭代的粒子数目由维度确定。

优选地，本发明提供的基于改进的自适应粒子群算法的岩体力学参数联合求解方法，还可以具有以下特征：在S2.3中，子种群划分方法为基于相对距离的划分，具体包括：

在空间中将随机生成Z＝Z_max个中心位置，根据每个粒子至每个中心的距离划分子种群，判断每个子种群中的粒子数是否达到最小值n，若否，则中心位置的数目Z减1，继续划分子种群，否则重新计算中心位置并判断与上次生成的中心位置相比是否发生变化。当发生变化则重新划分子种群，否则算法结束，获得划分结果。该过程将持续20次以获得多个划分结果。

优选地，本发明提供的基于改进的自适应粒子群算法的岩体力学参数联合求解方法，还可以具有以下特征：在S2.3中，Calinski-Harabasz(CH)指标由公式(3)进行计算，在获得的多个子种群划分结果中，挑选出获得CH指标最高分数的划分结果进行使用；

式中，n_z为第z个子种群的粒子数，s_z为第z个子种群的中心位置，s为全局中心位置。

优选地，本发明提供的基于改进的自适应粒子群算法的岩体力学参数联合求解方法，在S2.4中，将产生个体最优位置、子种群最优位置、全局最优位置，以及个体最小适应度值、子种群最小适应度值、全局最小适应度值，具体包括：每个粒子在每次迭代中，均会根据自身所在位置获得适应度值，单个粒子在迭代过程中获得最小适应度值时的位置称为个体最优位置，共N个；子种群内的最小适应度值所对应的位置称为子种群最优位置，共Z个；整个粒子群内的最小适应度值对应位置称为全局最优位置，共1个；以上3个层次的最小适应度值及最优位置随每次迭代进行更新。

优选地，本发明提供的基于改进的自适应粒子群算法的岩体力学参数联合求解方法，在步骤S2.4中，算法参数w、c₁、c₂和c₃将随计算效果进行自适应调整，具体包括：

第i个粒子在第t次迭代过程中，定义该迭代步的适应度值f_i ^t与上一迭代步的适应度值f_i ^t-1的差值为寻优梯度，根据寻优梯度按照公式(4)～(7)更新惯性权重w、第一加速系数c₁、第二加速系数c₂和第三加速系数c₃；

寻优梯度减小代表粒子搜索能力较弱，将增大w和c₁，减小c₂和c₃；寻优梯度增大代表粒子搜索能力较强，将减小w和c₁，增大c₂和c₃。

优选地，本发明提供的基于改进的自适应粒子群算法的岩体力学参数联合求解方法，还可以具有以下特征：在S2.4中，粒子位置按照公式(8)计算

式中，

表示第i个粒子第d维度在第t+1次迭代的位置分量；

表示第i个粒子第d维度在第t+1次迭代的速度分量。

<装置>

进一步，本发明还提供了能够自动实现上述<方法>的岩体力学参数联合求解装置，其特征在于，包括：

数据获取部，获取岩土工程中研究区域的监测数据，并根据岩土本构模型或/和渗流模型选取待求解的弹塑性参数或渗透参数；

自适应粒子群算法求解部，采用S2.1～S2.4构建变形模型并进行计算求解；

S2.1、基于有限元数据和步骤S1.1获取的监测数据构建变形模型：

式中，x₁，x₂，...，x_p为待反演的岩体力学参数；m为监测断面的个数；n为断面上监测点的个数；u_ij为断面上第i个监测断面上第j个监测点的有限元计算值；u_ij ^*为相应的实测值；w_i为第i个监测断面的权重系数；

S2.2、初始化粒子群算法：设置初始化参数，包括目标函数、维度D、种群大小N、最大迭代次数T_max、拓扑结构种类、子种群最大数量Z_max、子种群最小粒子数n、拓扑结构重组间隔T、搜索范围[X_min，X_max]、位置调整机制、粒子速度范围[-V_max，V_max]、速度调整机制、每个粒子的位置、每个粒子的速度、惯性权重初始值w_start、惯性权重终止值w_end、第一加速系数初始值c_1，start、第一加速系数终止值c_1，end、第二加速系数初始值c_2，start、第二加速系数终止值c_2，end、第三加速系数初始值c_3，start、第三加速系数终止值c_3，end；

将S2.1构建的变形模型作为目标函数，用于度量计算值与实际值的接近程度，将每个粒子的位置信息输入目标函数可获得适应度值，其值越小代表与岩体实际变形、渗流的接近程度越高，岩体参数越接近真实值；维度D即待求解的所有岩体力学参数的个数；种群大小N指粒子的个数，一个粒子在搜索空间中的位置代表一组岩体参数解；搜索范围[X_min，X_max]由具体的待求解参数确定；

S2.3、根据粒子在搜索空间中的位置，按照相对距离划分子种群，根据Calinski-Harabasz指标确定子种群个数；

S2.4、粒子群算法的迭代计算：根据粒子在搜索空间中的位置及目标函数计算每个粒子的适应度值，取一次迭代步中整个粒子群最小适应度值所对应的粒子位置作为全局最优位置，得到一组最优岩体参数；迭代次数加1，更新每个粒子的速度和位置，更新算法参数w、c₁、c₂和c₃，更新全局最优位置，判断迭代次数t是否为拓扑结构重组间隔T的整数倍，若是，则重新划分子种群，否则继续运算。随着迭代的进行，整个粒子群获得的最小适应度值呈不断下降的趋势，也即愈发接近工程实际，当整个种群的最小适应度值达到设定的精度要求，或迭代次数达到最大迭代次数T_max，算法停止；

在S2.4中的粒子速度按照下列公式(2)计算：

式中，

表示第i个粒子第d维度在第t次迭代的位置分量；

表示第i个粒子第d维度在第t次迭代的速度分量；w^t表示第t次迭代时的惯性权重；

表示第t次迭代中第i个粒子的第一加速系数，

表示第t次迭代中第i个粒子的第二加速系数，

表示第t次迭代中第i个粒子的第三加速系数，

表示第t次迭代中第i个粒子的个体最小适应度值；

表示第t次迭代中第z个子种群的子种群最小适应度值；

表示第t次迭代中的全局最小适应度值，r₁、r₂、r₃为[0，1]范围内的随机数；

岩体参数获取部，在S2.4算法停止后，取出整个种群中具有最小适应度值的粒子，该粒子在搜索空间中的位置即为求解的所有岩体力学参数；以及

控制部，与数据获取部、自适应粒子群算法求解部、岩体参数获取部均通信相连，控制它们的运行。

优选地，本发明提供的岩体力学参数联合求解装置，还可以包括：输入显示部，与数据获取部、自适应粒子群算法求解部、岩体参数获取部、控制部均通信相连，根据用户输入的操作指令，显示相应的信息。

优选地，本发明提供的岩体力学参数联合求解装置，还可以具有以下特征：输入显示部能够对岩体参数获取部求解得到的所有岩体力学参数以列表形式进行显示，并能够将这些岩体力学参数在相应的岩体三维或者平面图中相应位置处进行显示，还能够以静态曲线图或者动态变化图的方式按照时间顺序显示出岩体力学参数的变化过程，从而直观地反应工程应力、变形和渗流情况。

发明的作用与效果

1.本发明以工程实际监测数据作为输入，该数据能够代表一块区域内岩体的真实材料性能，具有平均意义，克服了室内试验和现场测试测点离散、代表性不强的缺陷。

2.本发明采用的监测数据可随时间逐渐延长，从而可以分析出岩体参数随时间变化的趋势，能够判断岩体是否达到稳定，是室内试验和现场测试无法做到的。

3.本发明能够根据搜索状态高频、自动地调整参数和种群拓扑结构，对粒子的全局搜索和局部挖掘能力进行动态平衡，避免了人力的过度干预，保证了整体搜索性能的稳定。

4.本发明对速度公式进行了改进。首先，引入新的局部极值和学习因子，这有利于粒子之间的信息交流。其次，对各个参数均进行非线性参数调整，考虑到了搜索的不同状态，为粒子前期探索和后期开发分配不同的参数大小。最后，参数的非线性调整策略融入了前后迭代步间适应度值的变化梯度信息，避免了粒子群的过早局部收敛，保障了整个搜索过程的有效性。整个速度改进方案将有效提高算法的准确度和鲁棒性，保障算法面对不同的求解情形时具备稳定的求解能力和求解精度。

5.本发明在考虑粒子间信息交流的同时，为每个粒子赋予了独立的搜索参数，不会过早形成统一的运动趋势，在很大程度上避免了早熟收敛，在整个算法运算过程中保证广泛搜索，获得的岩体力学参数也会更接近真实岩体状态。

6.本发明将整个种群进行了动态划分，各个子群内的信息交流保障了整个种群不会向一个方向快速聚集，该方法能够有效保障搜索区域的多样性，对求解空间的探索更加充分，改善了粒子群算法自身的缺陷。对于单个粒子，其在自身最优位置、子群最优位置和全局最优位置三个信息的平衡引导下进行搜索，在不同的搜索阶段将处于不同子群内，改善了传统粒子群算法的信息获取途径。

7.本发明方法能够在同一时间求解多个区域的多个模型参数，对于高维度、高非线性问题的求解效率高，与常规优化算法相比求解能力整体提升，获得的岩体参数也更接近实际情况。

综上，本方法能够为岩土工程的变形、受力、渗流试验和数值模拟快速提供高准确度的岩体力学参数，保障工程全生命周期的安全运行。

附图说明

图1为本发明涉及的基于改进的自适应粒子群算法的岩体力学参数联合求解方法的流程图；

图2为本发明实施例中某岩土工程的材料分区及监测点分布图；

图3为本发明实施例中Calinski-Harabasz指标确定子种群数目过程图；

图4为本发明实施例中自适应粒子群算法的全局最小值迭代曲线；

图5为本发明实施例中通过本发明提出的方法获得的岩体参数与试验参数的有限元计算结果和实际监测数据的对比图；

图6为比较例中本发明方法在Ackley函数下的全局最优适应度值的演化曲线与其他粒子群算法的对比图；

图7为比较例中本发明方法在Rastrigin函数下的全局最优适应度值的演化曲线与其他粒子群算法的对比图；

图8为比较例中本发明方法在Discus函数下的全局最优适应度值的演化曲线与其他粒子群算法的对比图；

图9为比较例中本发明方法在Sphere函数下的全局最优适应度值的演化曲线与其他粒子群算法的对比图。

具体实施方式

以下结合附图对本发明涉及的基于改进的自适应粒子群算法的岩体力学参数联合求解方法及装置的具体实施方案进行详细地说明。

<实施例>

如图1所示，本实施所提供的基于自适应粒子群算法的岩体力学参数联合求解方法，包括如下步骤：

S1、数据获取：

S1.1、在本实施例中，参见图2，选择某堆石坝的堆石区B区和堆石区C区作为研究区域，根据布置在坝体上的监测仪器获取变形监测数据。

S1.2、根据弹塑性本构模型中的邓肯-张E-B模型，选择求解弹塑性参数K、n、K_b和m，由于研究区域有2部分，因此共需求解8个岩体参数，限制4种参数的搜索区域分别为[1000，3000]、[0.1，0.5]、[1000，3000]和[0.01，0.2]。

S2.1、基于监测数据和有限元数据构建的变形模型如公式(A1)所示，该变形模型也即自适应粒子群算法的目标函数，通过该函数获得的适应度值越小则越接近实际变形情况。

式中，x₁，x₂，...，x_p为待反演的岩体力学参数，该实施例中为K₁、n₁、K_b1、m₁、K₂、n₂、K_b2、m₂；m为监测断面的个数，该实施例中为1；n为断面上监测点的个数，该实施例中为15；u_ij为断面上第i个监测断面上第j个监测点的有限元计算值；u_ij ^*为相应的实测值；w_i为第i个监测断面的权重系数，该实施例中为1。

S2.2、进行粒子群算法的初始化操作。目标函数已在步骤S2.1中列出；维度D即待求解的岩体参数总个数，该实施例中为8；种群大小N设为50；最大迭代次数T_max设为2000；拓扑结构种类设为星形结构；子种群最大数量Z_max设为7；子种群最小粒子数设为6；拓扑结构重组间隔T设为200；搜索范围已在S1.2列出；位置调整机制设置为②；粒子速度范围设置为[-5，5]；速度调整机制设置为④；惯性权重初始值w_start设置为3；惯性权重终止值w_end设置为2；第一加速系数初始值c_1，start设置为2；第一加速系数终止值c_1，end设置为1；第二加速系数初始值c_2，start设置为1；第二加速系数终止值c_2，end设置为2；第三加速系数初始值c_3，start设置为1；第三加速系数终止值c_3，end设置为2。在搜索空间中随机生成粒子位置和速度。

其中，拓扑结构种类包含金字塔形结构、环形结构、星形结构和冯诺依曼结构，设置多种拓扑结构用于应对不同的非线性模型，初始化时需选择其中1种。金字塔形结构是一种多层结构，顶层粒子少，底层粒子多，表现为粒子的速度迭代仅受上、下层粒子影响；环形结构表现为在算法开始前对所有粒子进行编号，粒子的速度迭代仅受编号在其前、后的两个粒子影响；星形结构即每个粒子的速度迭代受全局其他所有粒子影响；冯诺依曼结构与环形结构类似，但影响速度迭代的粒子数目由维度确定。

速度调整方式包含2种，初始化时需选择其中1种：①将速度调整为上一迭代步的速度；②速度方向反转并使数值减小一半。

位置调整方式包含5种，初始化时需选择其中1种：①粒子位置调整至边界；②在搜索区域随机生成；③收缩粒子速度10倍，使其落入搜索范围内；④根据边界镜像反射，使位置落入搜索范围；⑤将粒子位置定义到上一迭代步位置与边界的中点。由于位置代表了待求解的岩体参数，因此设置多种位置调整方式以应对不同大小的搜索范围。

S2.3、进行子种群划分：在空间中将随机生成Z＝Z_max＝9个中心位置，根据每个粒子至每个中心的距离划分子种群，判断每个子种群中的粒子数是否达到最小值n＝5，若否，则中心位置的数目Z减1，继续划分子种群，否则重新计算中心位置并判断与上次生成的中心位置相比是否发生变化。当发生变化则重新划分子种群，否则算法结束，获得划分结果。运行多次获得多个划分结果。Calinski-Harabasz(CH)指标由公式(A2)进行计算，在获得的多个子种群划分结果中，挑选出获得CH指标最高分数的划分结果进行使用，子种群数目选择过程参见图3。

式中，n_z为第z个子种群的粒子数，s_z为第z个子种群的中心位置，s为全局中心位置。

子种群划分完毕后进入S2.4；

S2.4、对粒子群算法进行迭代：

步骤S2.4.1、根据粒子在搜索空间中的位置及目标函数计算每个粒子的适应度值，适应度值的数值越小代表粒子所在位置的岩体参数越接近实际参数。单个粒子在迭代历史中的最小适应度值对应的位置为个体最优位置，对于粒子数为N的种群，个体最优位置和个体最优位置有N个，该实施例中共50个；子种群内的粒子在迭代历史中的最小适应度值对应的位置为子种群最优位置，共Z个子种群最优位置和子种群最小适应度值，在该实施例中不超过9个；整个粒子群最小适应度值所对应的粒子位置为全局最优位置，也即一组最优岩体参数，该实施例中共1个。

步骤S2.4.2、迭代次数加1，采用公式(A3)～(A4)更新每个粒子的速度和位置。粒子速度范围限制在[-5，5]，当速度过大将触发调整机制②；粒子位置范围按照步骤S1.2中的搜索范围进行限制，当位置超出搜索范围将触发调整机制④。

式中，

表示第i个粒子第d维度在第t+1次迭代的位置分量；

表示第i个粒子第d维度在第t+1次迭代的速度分量；w^t表示第t次迭代时的惯性权重；

表示第t次迭代中第i个粒子的第一加速系数，

表示第t次迭代中第i个粒子的第二加速系数，

表示第t次迭代中第i个粒子的第三加速系数，

表示第t次迭代中第i个粒子的个体最小适应度值；

表示第t次迭代中第z个子种群的子种群最小适应度值；

表示第t次迭代中的全局最小适应度值，r₁、r₂、r₃为[0，1]范围内的随机数。

S2.4.3、采用公式(A5)～(A8)更新算法参数w、c₁、c₂和c₃。第i个粒子在第t次迭代过程中，定义该迭代步的适应度值f_i ^t与上一迭代步的适应度值f_i ^t-1的差值为寻优梯度，根据寻优梯度更新惯性权重w、第一加速系数c₁、第二加速系数c₂和第三加速系数c₃。寻优梯度减小代表粒子搜索能力较弱，将增大w和c₁，减小c₂和c₃；寻优梯度增大代表粒子搜索能力较强，将减小w和c₁，增大c₂和c₃。在该实施例中，w_start＝3，w_end＝2，c_1，start＝2，c_1，end＝1，c_2，start＝1，c_2，end＝2，c_3，start＝1，c_3，end＝2。

S2.4.3、根据新的位置计算每个粒子的适应度值，更新3个层次的最小适应度值和最优位置。

S2.4.4、判断迭代次数tt是否为拓扑结构重组间隔T＝200的整数倍，若是，则重新划分子种群，若否，则算法继续。在该实施例中，算法开始时进行1次划分，在开始迭代后运行至达到最大迭代次数T_max＝2000的过程中进行9次划分(t＝2000时算法结束，不进行划分)，共计10次。

S2.4.5、判断整个种群的最小适应度值是否达到设定的精度要求，或迭代次数是否达到最大迭代次数T_max，若否，则返回S2.4.1继续计算，若是，则算法停止。算法的全局最优适应度值收敛曲线参见图4。

S3、岩体参数的获取

根据步骤S2.4.5，在算法停止后，取出整个种群中具有最小适应度值的粒子，该粒子在搜索空间中的位置即为求解的岩体参数，即所有研究区域的所有力学参数，详见表1。

表1求解结果

该参数用于工程变形、渗流试验及模拟计算，控制工程运行状态。通过以上参数进行有限元模拟计算，与实际监测变形数据和试验参数的有限元计算结果进行对比，参见图5，可以看到，由该方法获取的参数进行模拟计算与监测数据吻合，证明该方法获取的岩体参数与实际情况基本一致。

进一步，本实施例还提供了基于以上方法自动进行岩体力学参数联合求解的装置，该装置包括数据获取部、自适应粒子群算法求解部、岩体参数获取部、输入显示部以及控制部。

数据获取部执行上文S1描述的内容，获取岩土工程中研究区域的监测数据，并根据岩土本构模型或渗流模型选取待求解的弹塑性参数或渗透参数。

自适应粒子群算法求解部执行上文S2描述的内容，采用S2.1～S2.4构建变形模型并进行计算求解。

岩体参数获取部执行上文S3描述的内容，在S2.4算法停止后，取出整个种群中具有最小适应度值的粒子，该粒子在搜索空间中的位置即为求解的所有岩体力学参数。

输入显示部与数据获取部、白适应粒子群算法求解部、岩体参数获取部均通信相连，根据用户输入的操作指令，显示相应的信息。例如，输入显示部能够对岩体参数获取部求解得到的所有岩体力学参数以列表形式进行显示，并能够将这些岩体力学参数在相应的岩体三维或者平面图中相应位置处进行显示，还能够以静态曲线图或者动态变化图的方式按照时间顺序显示出岩体力学参数的变化过程，从而直观地反应工程应力、变形和渗流情况。

控制部与数据获取部、自适应粒子群算法求解部、岩体参数获取部、输入显示部均通信相连，控制它们的运行。

<比较例>

本比较例通过实验对上述实施例中的自适应粒子群求解方法的效果进行验证，并与其他现有技术方法进行比较。

实验测试环境为Intel(R)Core(TM)i7-10700 CPU@2.90GHz，RAM 32GB，windows10操作系统。将本发明提出的自适应粒子群算法(Patent-PSO)与标准粒子群算法(PSO)、繁殖粒子群算法(BreedPSO)、线性变化加速系数的粒子群算法(HPSO-TVAC)、随机加速系数的粒子群算法(RandPSO)、线性变化惯性权重的粒子群算法(WPSO)进行对比分析。

在比较前对各个算法进行调试，获得每个算法性能最优时的参数。对于PSO，参数设置为w＝0.9，c₁＝2，c₂＝2；对于BreedPSO，参数设置为w＝0.9，c₁＝2，c₂＝2，Pc＝0.9，Sp＝0.2；对于HPSO-TVAC，参数设置为w＝0.9，c₁由2线性递减至0.5，c₂由0.5线性递增至2；对于RandPSO，参数设置为w＝0.9，c₁＝1.5*rand[0，1]，c₂＝1.5*rand[0，1]；对于WPSO，参数设置为w由0.9线性递减至0.4，c₁＝2，c₂＝2；对于Patent-PSO，参数设置为w_start＝3，w_end＝2，c_1，start＝2，c_1，end＝1，c_2，start＝1，c_2，end＝2，c_3，start＝1，c_3，end＝2。

在本实施例中，将最大迭代次数T_max设置为2000，维度设置为50维，以下表2中列出的几种目标函数为例，进行适应度值的计算，搜寻全局最小值，在列出的目标函数中全局最小值均为0，但由于函数较为复杂，算法只能逼近0却无法到达该数值，算法的接近程度代表了该算法的搜索能力。对于每个目标函数，每个算法运行30次取平均值。

表2目标函数

图6～图9列出了几种算法的运行结果，可以看到，Patent-PSO在高维度条件下的搜索能力显著高于其他对比算法，搜寻精度获得提高，很好地避免了其他算法易陷入局部最优解的缺陷。因此，本发明提供的Patent-PSO可以精确地获取真实岩体参数，服务工程实际。

以上实施例仅仅是对本发明技术方案所做的举例说明。本发明所涉及的岩体力学参数联合求解方法及装置并不仅仅限定于在以上实施例中所描述的内容，而是以权利要求所限定的范围为准。本发明所属领域技术人员在该实施例的基础上所做的任何修改或补充或等效替换，都在本发明的权利要求所要求保护的范围内。

Claims (10)

1.岩体力学参数联合求解方法，其特征在于，包括以下步骤：

S1、数据获取：

S1.1、选择岩土工程中一块或多块研究区域，获取该区域的监测数据；

S1.2、根据岩土本构模型或渗流模型选取待求解的弹塑性参数或渗透参数；

S2、自适应粒子群算法求解：

S2.1、基于有限元数据和步骤S1.1获取的监测数据构建变形模型：

式中，x₁，x₂，...，x_p为待反演的岩体力学参数；m为监测断面的个数；n为断面上监测点的个数；u_ij为断面上第i个监测断面上第j个监测点的有限元计算值；u_ij ^*为相应的实测值；w_i为第i个监测断面的权重系数；

S2.2、初始化粒子群算法：设置初始化参数，包括目标函数、维度D、种群大小N、最大迭代次数T_max、拓扑结构种类、子种群最大数量Z_max、子种群最小粒子数n、拓扑结构重组间隔T、搜索范围[X_min，X_max]、位置调整机制、粒子速度范围[-V_max，V_max]、速度调整机制、每个粒子的位置、每个粒子的速度、惯性权重初始值ω_start、惯性权重终止值ω_end、第一加速系数初始值c_1，start、第一加速系数终止值c_1，end、第二加速系数初始值c_2，start、第二加速系数终止值c_2，end、第三加速系数初始值c_3，start、第三加速系数终止值c_3，end；

将S2.1构建的变形模型作为目标函数，用于度量计算值与实际值的接近程度，将每个粒子的位置信息输入目标函数可获得适应度值，其值越小代表与岩体实际变形、渗流的接近程度越高，岩体参数越接近真实值；维度D即待求解的所有岩体力学参数的个数；种群大小N指粒子的个数，一个粒子在搜索空间中的位置代表一组岩体参数解；搜索范围[X_min，X_max]由具体的待求解参数确定；

S2.3、根据粒子在搜索空间中的位置，按照相对距离划分子种群，根据Calinski-Harabasz指标确定子种群个数；

S2.4、粒子群算法的迭代计算：根据粒子在搜索空间中的位置及目标函数计算每个粒子的适应度值，取一次迭代步中整个粒子群最小适应度值所对应的粒子位置作为全局最优位置，得到一组最优岩体参数；迭代次数加1，更新每个粒子的速度和位置，更新算法参数w、c₁、c₂和c₃，更新全局最优位置，判断迭代次数t是否为拓扑结构重组间隔T的整数倍，若是，则重新划分子种群，否则继续运算；随着迭代的进行，整个粒子群获得的最小适应度值呈不断下降的趋势，也即愈发接近工程实际，当整个种群的最小适应度值达到设定的精度要求，或迭代次数达到最大迭代次数T_max，算法停止；

在S2.4中的粒子速度按照下列公式(2)计算：

式中，

表示第i个粒子第d维度在第t次迭代的位置分量；

表示第i个粒子第d维度在第t次迭代的速度分量；w^t表示第t次迭代时的惯性权重；

表示第t次迭代中第i个粒子的第一加速系数，

表示第t次迭代中第i个粒子的第二加速系数，

表示第t次迭代中第i个粒子的第三加速系数，

表示第t次迭代中第i个粒子的个体最小适应度值；

表示第t次迭代中第z个子种群的子种群最小适应度值；

表示第t次迭代中的全局最小适应度值，r₁、r₂、r₃为[0，1]范围内的随机数；

S3、岩体参数的获取：

在S2.4算法停止后，取出整个种群中具有最小适应度值的粒子，该粒子在搜索空间中的位置即为求解的所有岩体力学参数。

2.根据权利要求1所述的岩体力学参数联合求解方法，其特征在于：

其中，在S1.1中，监测数据包含变形数据和渗流数据中的至少一种；

在S1.2中，根据岩土本构模型或/和渗流模型的至少一种选取待求解的弹塑性参数或/和渗透参数；所述岩土本构模型包括弹性模型、弹塑性模型、粘弹塑性模型。

3.根据权利要求1所述的岩体力学参数联合求解方法，其特征在于：

其中，在S2.1中，所述有限元数据为通过使用多组不同岩体力学参数的有限元计算获得的与监测点位置对应的变形、渗流数据。

4.根据权利要求1所述的岩体力学参数联合求解方法，其特征在于：

其中，在S2.2中，粒子的速度和位置受约束，数值超过限定的范围将触发调整机制，具体包括：

粒子速度范围限制在[-V_max，V_max]，当速度过大将触发调整机制，包括2种方式：①将速度调整为上一迭代步的速度；②速度方向反转并使数值减小一半；

粒子位置范围限制在[X_min，X_max]，当位置超出搜索范围将触发调整机制，包括5种方式：①粒子位置调整至边界；②在搜索区域随机生成；③收缩粒子速度10倍，使其落入搜索范围内；④根据边界镜像反射，使位置落入搜索范围内；⑤将粒子位置定义到上一迭代步位置与边界的中点；

在初始化时应选择1种速度调整机制和1种位置调整机制。

5.根据权利要求1所述的岩体力学参数联合求解方法，其特征在于：

其中，在S2.3中，子种群划分方法为基于相对距离的划分，具体包括：

在空间中将随机生成Z＝Z_max个中心位置，根据每个粒子至每个中心的距离划分子种群，判断每个子种群中的粒子数是否达到最小值n，若否，则中心位置的数目Z减1，继续划分子种群，否则重新计算中心位置并判断与上次生成的中心位置相比是否发生变化；当发生变化则重新划分子种群，否则算法结束，获得划分结果；该过程将进行多次，保存多个子种群划分结果。

6.根据权利要求1所述的岩体力学参数联合求解方法，其特征在于：

其中，在S2.3中，Calinski-Harabasz指标由公式(3)进行计算，在获得的多个子种群划分结果中，挑选出获得Calinski-Harabasz指标最高分数的划分结果进行使用；

式中，n_z为第z个子种群的粒子数，s_z为第z个子种群的中心位置，s为全局中心位置。

7.根据权利要求1所述的岩体力学参数联合求解方法，其特征在于：

其中，在S2.4中，将产生个体最优位置、子种群最优位置、全局最优位置，以及个体最小适应度值、子种群最小适应度值、全局最小适应度值，具体包括：

每个粒子在每次迭代中，均会根据自身所在位置获得适应度值，单个粒子在迭代过程中获得最小适应度值时的位置称为个体最优位置，共N个；子种群内的最小适应度值所对应的位置称为子种群最优位置，共Z个；整个粒子群内的最小适应度值对应位置称为全局最优位置，共1个；以上3个层次的最小适应度值及最优位置随每次迭代进行更新。

8.根据权利要求1所述的岩体力学参数联合求解方法，其特征在于：

其中，在步骤S2.4中，算法参数w、c₁、c₂和c₃将随计算效果进行白适应调整，具体包括：

第j个粒子在第t次迭代过程中，定义该迭代步的适应度值f_i ^t与上一迭代步的适应度值f_i ^t-1的差值为寻优梯度，根据寻优梯度按照公式(4)～(7)更新惯性权重w、第一加速系数c₁、第二加速系数c₂和第三加速系数c₃；

寻优梯度减小代表粒子搜索能力较弱，将增大w和c₁，减小c₂和c₃；寻优梯度增大代表粒子搜索能力较强，将减小ω和c₁，增大c₂和c₃。

9.岩体力学参数联合求解装置，其特征在于，包括：

数据获取部，获取岩土工程中研究区域的监测数据，并根据岩土本构模型或渗流模型选取待求解的弹塑性参数或渗透参数；

自适应粒子群算法求解部，采用S2.1～S2.4构建变形模型并进行计算求解；

S2.1、基于有限元数据和步骤S1.1获取的监测数据构建变形模型：

式中，x₁，x₂，...，x_p为待反演的岩体力学参数；m为监测断面的个数；n为断面上监测点的个数；u_ij为断面上第i个监测断面上第j个监测点的有限元计算值；u_ij ^*为相应的实测值；w_i为第i个监测断面的权重系数；

S2.2、初始化粒子群算法：设置初始化参数，包括目标函数、维度D、种群大小N、最大迭代次数T_max、拓扑结构种类、子种群最大数量Z_max、子种群最小粒子数n、拓扑结构重组间隔T、搜索范围[X_min，X_max]、位置调整机制、粒子速度范围[-V_max，V_max]、速度调整机制、每个粒子的位置、每个粒子的速度、惯性权重初始值ω_start、惯性权重终止值ω_end、第一加速系数初始值c_1，start、第一加速系数终止值c_1，end、第二加速系数初始值c_2，start、第二加速系数终止值c_2，end、第三加速系数初始值c_3，star、第三加速系数终止值c_3，end；

将S2.1构建的变形模型作为目标函数，用于度量计算值与实际值的接近程度，将每个粒子的位置信息输入目标函数可获得适应度值，其值越小代表与岩体实际变形、渗流的接近程度越高，岩体参数越接近真实值；维度D即待求解的所有岩体力学参数的个数；种群大小N指粒子的个数，一个粒子在搜索空间中的位置代表一组岩体参数解；搜索范围[X_min，X_max]由具体的待求解参数确定；

S2.3、根据粒子在搜索空间中的位置，按照相对距离划分子种群，根据Calinski-Harabasz指标确定子种群个数；

S2.4、粒子群算法的迭代计算：根据粒子在搜索空间中的位置及目标函数计算每个粒子的适应度值，取一次迭代步中整个粒子群最小适应度值所对应的粒子位置作为全局最优位置，得到一组最优岩体参数；迭代次数加1，更新每个粒子的速度和位置，更新算法参数w、c₁、c₂和c₃，更新全局最优位置，判断迭代次数t是否为拓扑结构重组间隔T的整数倍，若是，则重新划分子种群，否则继续运算；随着迭代的进行，整个粒子群获得的最小适应度值呈不断下降的趋势，也即愈发接近工程实际，当整个种群的最小适应度值达到设定的精度要求，或迭代次数达到最大迭代次数T_max，算法停止；

在S2.4中的粒子速度按照下列公式(2)计算：

式中，

表示第i个粒子第d维度在第t次迭代的位置分量；

表示第i个粒子第d维度在第t次迭代的速度分量；w^t表示第t次迭代时的惯性权重；

表示第t次迭代中第i个粒子的第一加速系数，

表示第t次迭代中第i个粒子的第二加速系数，

表示第t次迭代中第i个粒子的第三加速系数，

表示第t次迭代中第i个粒子的个体最小适应度值；

表示第t次迭代中第z个子种群的子种群最小适应度值；

表示第t次迭代中的全局最小适应度值，r₁、r₂、r₃为[0，1]范围内的随机数；

岩体参数获取部，在S2.4算法停止后，取出整个种群中具有最小适应度值的粒子，该粒子在搜索空间中的位置即为求解的所有岩体力学参数；以及

控制部，与所述数据获取部、所述自适应粒子群算法求解部、所述岩体参数获取部均通信相连，控制它们的运行。

10.根据权利要求9所述的岩体力学参数联合求解装置，其特征在于，还包括：

输入显示部，与所述数据获取部、所述自适应粒子群算法求解部、所述岩体参数获取部、所述控制部均通信相连，根据用户输入的操作指令，显示相应的信息。

Images