粒子群算法在时效变形参数反演中的应用

杨文东^{1, 2}，张强勇²，李术才²，王刚³，李勇²

(1. 中国石油大学(华东) 储运与建筑工程学院，山东 青岛，266580；

2. 山东大学 岩土与结构工程研究中心，山东 济南，250061；

3. 山东科技大学 山东省土木工程防灾减灾重点实验室，山东 青岛，266590)

摘要：以往的现场刚性承压板压缩蠕变试验中都假设岩体变形为黏弹性变形，然后推求岩体的蠕变变形公式，并结合最小二乘法反演获得蠕变参数，这种方法由于进行了黏弹性变形的假设从而无法考虑岩体的黏塑性变形，为此，采用粒子群智能算法的数值反演方法进行研究，并且提出将流变模型中控制瞬时变形和时效变形的2种参数分开反演的二次粒子群算法。研究结果表明：采用这种方法可以有效地减小反演的难度，提高了反演精度；将此方法应用于现场压缩蠕变试验的参数反演，拟合时效变形参数曲线与试验曲线较吻合，从而可获得相应的流变参数。

关键词：岩石力学；流变模型；参数反演；粒子群算法

中图分类号：O 319.56           文献标志码：A         文章编号：1672-7207(2013)01-0282-07

Application of particle swarm optimization in time-dependent parameters inversion

YANG Wendong^{1, 2}, ZHANG Qiangyong², LI Shucai², WANG Gang³, LI Yong²

(1. College of Pipeline and Civil Engineering, China University of Petroleum, Qingdao 266580, China;

2. Research Center of Geotechnical and Structural Engineering, Shandong University, Jinan 250061, China;

3. Shandong Provincial Key Laboratory of Civil Engineering Disaster Prevention and Mitigation,

Shandong University of Science and Technology, Qingdao 266590, China)

Abstract: In the previous research of in-situ compressive creep test using rigid bearing plate, the deformation of rock mass was always assumed to be visco-elastic deformation, and in this work the creep deformation equation is deduced, and the creep parameters are inversed with least-squares method. This method which assumed the deformation of rock mass in the compressive creep test is visco-elastic deformation does not consider the visco-plastic deformation. Based on the particle swarm intelligence algorithm, the numerical inversion method was used which was more in line with the actual situation.The secondary particle swarm optimization for rheological parameters inversion in which the instantaneous parameters and time-controlled parameters were separately inversed was proposed, and this method was applied in the in-situ compressive creep test. The results show that the methods can reduce the difficulty of the inversion and improve the accuracy. Example results verify the credibility of program compiling. The fitting curves are in good agreement with the experimental curves, which indicates that the model is correct and reasonable, and the rheological parameters are inversed.

Key words: rock mechanics; rheological model; parameter inversion; particle swarm optimization

岩石本构模型及力学参数的正确性是确保岩石工程设计安全可靠的重要保障，对流变问题更是如此。与室内的三轴流变试验相比，现场的大型承压板试验更能反映真实岩体的力学性质，但是，由于承压板下岩体的黏弹塑性蠕变变形解析公式难于表达，往往假设仅发生黏弹性变形^[1-3]，与真实情况相差较大。而采用嵌入智能反演算法的数值反演方法可以解决此问题。早期的岩石力学反演问题大多采用直接法、最小二乘法和优化反演法^[4-5]，但是，计算时间长，且当未知参数较多时，收敛速度慢，解的稳定性差，常常陷入局部最优解而不能得到全局最优解。随着数学算法的发展和计算机水平的提高，一些智能算法如神经网络法、遗传算法、模拟退火算法、支持向量机、蚁群算法和粒子群算法等^[6-17]逐渐产生并应用于岩石力学反问题中。目前，在压缩蠕变试验中，岩体蠕变参数的反演方法主要有2种：一是基于黏弹性变形解析公式的理论反演法，即通过推导出承压板下岩体的黏弹性变形公式，然后利用最小二承法反演获得蠕变参数；二是根据现场承压板下岩体的分布状况建立其网格模型，利用黏弹塑性流变模型进行数值计算，结合相应的反演算法获得其蠕变参数。在理论反演法中，承压板压缩蠕变试验中复杂黏弹塑流变模型的变形公式难于推导，只能简化为黏弹性变形；而采用数值反演法则可解决该问题，并可以反演获得复杂的黏弹塑性模型参数，与实际结果更相符，但是，数值反演法效率较低。以往的研究中大多采用理论反演法，无法处理复杂的黏弹塑性流变模型^[1-3]；而采用粒子群算法时也将所有流变参数统一对待，反演精度和效率有待提高^[12-13]。为此，本文作者针对现场压缩蠕变试验的流变参数反演过程中参数较多、反演难度较大、耗时较长等问题，提出将流变模型中控制瞬时变形和时效变形的2种参数分开反演的二次粒子群算法，既可以用于反演复杂流变模型的蠕变参数，又提高了反演效率。

1  时效参数反演的粒子群算法

受鸟和昆虫自然群聚和觅食行为的启发，Kennedy等^[18-19]提出了最早的粒子群算法(particle swarm optimization, PSO)。该算法通过群体的信息共享和个体自身经验的总结来修正个体行动策略，以适应度函数指导随机化搜索方向，最终找到复杂搜索空间中的最优区域。

1.1  基本粒子群算法

设粒子群规模为N，其中每个粒子在M维空间中的坐标位置表示为：

；i=1, 2, …, N     (1)

粒子i的速度定义为每次迭代中粒子移动的距离：

；i=1, 2, …, N     (2)

粒子i(i=1, 2, …, N)在第d(i=1, 2, …, M)维子空间中的飞行速度v_id根据下式进行调整：

          (3)

粒子通过下式调整自身的位置：

               (4)

式中：为当前粒子的历史最优位置记录；为整个粒子群的历史最优位置记录；为惯性权重；c₁和c₂为加速常数，c₁为认知参数，c₂为社会参数；和为(0, 1)之间的随机数。

粒子的运动由上述方程共同作用，其运动速度与其历史飞行经验和群体飞行经验相关^[16-17]。粒子群优化示意图见图1。

图1  粒子群优化示意图

Fig.1  Schematic diagram of particle swarm optimization

1.2  时效参数反演的粒子群算法

流变模型的参数往往较多，而参数反演过程中每增加1个参数将使工作量成倍增加、难度增大并且精确性下降。待求解参数的增加也使问题的多解性增加，因此，极易陷于局部最优解，而无法求得全局最优解。为解决这一问题，应尽可能地减少待反演参量，以减小后序工作的难度。

流变参数很多，但可分为2类：一类决定岩体的瞬时变形行为，另一类决定岩体的时效变形。这2类参数共同决定了岩体在外荷作用下的蠕变曲线。决定岩体的瞬时变形行为的参数称为瞬时参数，岩体的瞬时变形行为实际上也就是蠕变过程中时间取0时的变形行为。基于此，可以把流变参数反演分为2步：(1) 先按弹塑性计算确定岩体的瞬时参数；(2) 将求出的瞬时参数作为已知量代入蠕变计算式确定其他蠕变参数。在2步反演过程中，仍采用粒子群算法(PSO)进行求解，这可大大减少流变参数反演的难度，并提高反演精度。

1.3  时效参数反演中的粒子群算法流程

时效参数反演的粒子群算法流程是将粒子种群分为2个：一个为瞬时参数的种群1，另一个为蠕变参数的种群2，分步骤进行反演，先进行种群1的搜索，再进行种群2的搜索。图2所示为具体的时效参数反演的粒子群算法流程图，其中：N为粒子群规模；i为粒子数编号；v_i为粒子i的速度；X_i为粒子i的位置；(P_best)_i为第i个粒子飞行所经历的最好位置；G_best为所有粒子飞行经过的最好位置；Q为最大迭代次数；q为当前迭代次数编号；fitnessX_i为粒子的适应值；fitness(P_best)_i为P_i所对应的适应值；fitnessG_best为G所对应的适应值；ε为合适的适应值。

图2  时效参数反演的粒子群算法流程图

Fig.2  Flow chart of improved PSO calculation program for time-dependent parameters

1.4  粒子群边界处理方法

Croes^[20]针对出现的此类情况提出3种不同边界处理方法：

(1) 吸收墙(absorbing walls)。当粒子在某一维上越出解空间的边界时，在该维方向上的速度被置为0，粒子最终被牵引回可行解的区域，这好像边界吸收了粒子向外逃逸的能量一样，如图3(a)所示。

(2) 反射墙(reflecting walls)。当粒子在某一维上越出解空间的边界时，改变粒子在该维方向上的速度使粒子被反射回解空间，如图3(b)所示。

(3) 隐匿墙(invisible walls)。粒子允许越过物理限定区域，但是，对越界的粒子并不估计其适应度，如图3(c)所示。

图3  3种边界处理方法

Fig.3  Treatment method of three boundaries

本文采用反射墙的处理方法，即当某一粒子在飞行过程中超出解空间的边界时，把该粒子的位置置为边界值，并改变粒子速度方向，粒子最终被牵引回可行

解的区域。适应度函数为。

2  算例验证

通过小试块单轴压缩的流变计算算例对参数反演程序进行验证。

试件体积为10 cm(X向)×20 cm(Y向)×1 cm(Z向)，共划分200个单元和462个节点。模型在底部Y方向约束，顶部施加1个100 MPa的分布压力，如图4所示。

采用广义Kelvin模型，参数取值为：瞬时弹性模量E₀=50 GPa，泊松比μ=0.2，黏弹性剪切模量G₁=40 GPa，黏滞系数η₁=100 GPa·h；瞬时剪切模量；瞬时体积模量。取由这些参数计算所得的曲线为目标曲线。

假认泊松比μ为常量，则有E₀，G₁和η₁ 3个待反演变量，以此3个变量为待求参量进行反演。粒子群搜索空间为：40＜E₀＜60 GPa，30＜G₁＜50 GPa，50＜η₀＜150 GPa·h；粒子群规模N为10；最大搜索次数Q为8。经反演求得：E₀=49.394 GPa；G₁=41.094 GPa；η₁=98.116 GPa·h。3个参量广义Kelvin模型见图5。

图4  单轴压缩算例示意图

Fig.4  Schematic diagram of uniaxial compression example

图5  3个参量广义Kelvin模型

Fig.5  Generalized Kelvin model with three parameters

图6所示为适应度随迭代次数的变化结果，图7所示为反演曲线与目标曲线比较结果。

图6  适应度随迭代次数变化

Fig.6  Relationship between fitness and iteration times

图7  变形值反演曲线与目标曲线比较

Fig.7  Comparisons of time-dependent parameters between inversion curve and target curve

从图6和图7可见：反演曲线与目标曲线较吻合，反演参数值与原参数值相差不大；在粒子群迭代过程中，适应度随迭代次数增大而不断减小，说明参数随着迭代步不断优化逐渐接近最优解，由此验证了程序的正确性。

3  软弱岩带蠕变参数的反演

为获得某水电站坝基软弱岩带的流变力学参数，进行了现场的压缩蠕变试验^[21]。下面选取典型试验进行说明，E₀218-3承压板以下岩体地质描述见表1。

表1  E₀218-3承压板以下岩体地质描述

Table 1  Geological description of rock mass under bearing plate of E₀218-3

压缩蠕变参数数值反演计算范围(长×宽×厚)为16 m×16 m×16 m(该计算范围已远大于压缩荷载的影响范围)。

图8所示为PD218-3压缩蠕变参数数值反演计算网格。由于网格对称，为了计算方便，取1/4体进行计算。图中圆形承压板半径取50 cm，高度10 cm，辉绿岩脉软弱岩带取1.8 m，下部岩体为Ⅱ~Ⅲ₁类花岗岩，整个模型共剖分为12 308个单元和14 177个节点。图9所示为蠕变试验曲线。

图8  计算网格划分

Fig.8  Mesh sketch of numerical model

图9  蠕变试验分级曲线

Fig.9  Multi-stage creep curves

以p=0.5 MPa时的蠕变曲线为目标曲线进行反演，采用损伤流变模型(HKKP模型)^[21]，如图10所示。

假认泊松比μ为常量，则只有1个瞬时变量E₀，以此变量为待求参量进行反演。粒子群搜索空间为：0＜E₀＜1 000 MPa；粒子群规模N为5；最大搜索次数Q为10。由PSO1反演求得：E₀=56.21 MPa。

图10  损伤流变模型(HKKP模型)

Fig.10  Damage rheological model (HKKP model)

图11所示为瞬时变量反演过程中适应度随迭代次数的变化过程。从图11可见：随着迭代的进行，适应度逐渐减小，说明待求参量逐步向最优解逼近。

图11  瞬时变量反演过程中适应度随迭代次数变化

Fig.11  Relationship between fitness and iteration times in inverse process of instant parameters

有G₁，η₁，G₂，η₂和共5个蠕变变量，以此5个变量为待求参量进行反演(其中：G₁和G₂分别为第1个和第2个Kelvin体的黏弹剪切模量；η₁和η₂分别为第1个和第2个Kelvin体的黏滞系数；为损伤参数)。粒子群搜索空间为：0＜G₁＜100 MPa；0＜η₁＜100 MPa·h；0＜G₂＜100 MPa；0＜η₂＜100 MPa·h；0＜＜1；粒子群规模N为3；最大搜索次数Q为10。损伤部分的由如下公式计算^[22]：

           (5)

其中：E₀和分别为瞬时弹性模量和长期弹性模量；和分别为压缩蠕变的瞬时变形值和达到稳定后的最终变形值；q₀为施加的法向均布荷载；R₀为承压板半径。由PSO2反演求得：G₁=40.20 MPa；η₁=2.0 MPa·h；G₂=60.60 MPa；η₂=80.0 MPa·h；=0.55。

图12所示为蠕变变量反演过程中适应度随迭代次数的变化过程；图13所示为变形反演结果与试验结果的比较。从图13可见：变形反演曲线与试验曲线较吻合，说明反演效果较好。表2所示为反演获得各级曲线的流变参数(其中：p为法向均布荷载；ε_threshold为损伤阈值)。

图12  蠕变变量反演过程中适应度随迭代次数变化

Fig.12  Relationship between fitness and iteration times in inverse process of creep parameters

图13  蠕变反演曲线与试验曲线比较

Fig.13  Comparisons of creep between inversion curve and test curve

表2  反演获得的流变参数

Table 2  Rheological parameters obtained by inversion method

4  结论

(1) 以往的压缩蠕变试验的参数反演大多采用理论反演法，都假设岩体变形为黏弹性变形，本文采用粒子群智能算法的数值反演法，可以实现复杂流变模型的参数反演，能更准确地反映实际情况。

(2) 提出将流变模型中控制瞬时变形和时效变形的2种参数分开反演的二次粒子群算法，有效地减小了反演的难度，提高了效率。

(3) 通过算例验证了程序的正确性；利用现场压缩蠕变试验分级加载条件下的试验结果对其参数进行反演，获得了流变参数。

参考文献：

[1] 李云鹏, 王芝银, 丁秀丽. 流变荷载试验曲线的模型识别及其应用[J]. 石油大学学报, 2005, 29(2): 73-77.

LI Yunpeng, WANG Zhiying, DING Xiuli. Model identification for rheological load test curve and its application[J]. Journal of University of Petroleum, 2005, 29(2): 73-77.

[2] 徐平, 丁秀丽, 全海, 等. 溪洛渡水电站坝址区岩体蠕变特性试验研究[J]. 岩土力学, 2003, 24(增刊): 220-226.

XU Ping, DING Xiuli, QUAN Hai, et al. Testing study on creep behavior of rock mass at Xiluodu Dam site[J]. Rock and Soil Mechanics, 2003, 24(Suppl): 220-226.

[3] 杨文东, 张强勇, 张建国, 等. 刚性承压板下深部岩体压缩蠕变参数反演[J]. 岩土力学, 2009, 30(3): 762-768.

YANG Wendong, ZHANG Qiangyong, ZHANG Jianguo, et al. Inversion of compression creep parameters of deep rock under rigid bearing plate[J]. Rock and Soil Mechanics, 2009, 30(3): 762-768.

[4] 吕爱钟, 蒋斌松. 岩石力学反问题[M]. 北京: 煤炭工业出版社, 1998: 12-112.

L Aizhong, JIANG Binsong. Inverse problems of rock mechanics[M]. Beijing: China Coal Industry Press, 1998: 12-112.

[5] 杨林德. 岩土工程问题的反演理论与工程实践[M]. 北京: 科学出版社, 1996: 15-67.

YANG Linde. Inverse theory and engineering practice in geotechnical engineering[M]. Beijing: Science Press, 1996: 15-67.

[6] 冯夏庭. 智能岩石力学导论[M]. 北京: 科学出版社, 2000: 1-48.

FENG Xiating. Introduction to intelligent rock mechanics[M]. Beijing: Science Press, 2000: 1-48.

[7] Feng X, Chen B, Yang C, et al. Identification of visco-elastic models for rocks using genetic programming coupled with the modified particle swarm optimization algorithm[J]. International Journal of Rock Mechanics and Mining Sciences, 2006, 43(5): 789-801.

[8] 刘保国, 孙钧. 岩体粘弹性本构模型辨识的一种方法[J]. 工程力学, 1999, 16(1): 18-25.

LIU Baoguo, SUN Jun. An identification method of visco-elastic constitutive model of rock mass[J]. Engineering Mechanics, 1999, 16(1): 18-25.

[9] 高玮, 郑颖人. 采用快速遗传算法进行岩土工程反分析[J]. 岩土工程学报, 2001, 23(1): 120-122.

GAO Wei, ZHENG Yingren. Back analysis in geotechnical engineering based on fast-convergent genetic algorithm[J]. Chinese Journal of Geotechnical Engineering, 2001, 23(1): 120-122.

[10] 苏国韶, 冯夏庭. 基于粒子群优化算法的高地应力条件下硬岩本构模型的参数辨识[J]. 岩石力学与工程学报, 2005, 24(17): 3029-3034.

SU Guoshao, FENG Xiating. Parameter identification of constitutive model for hard rock under high in-situ stress condition using particle swarm optimization algorithm[J]. Chinese Journal of Rock Mechanics and Engineering, 2005, 24(17): 3092-3034.

[11] 陈炳瑞, 冯夏庭, 黄书岭, 等. 基于快速拉格朗日分析-并行粒子群算法的黏弹塑性参数反演及其应用[J]. 岩石力学与工程学报, 2007, 26(12): 2517-2525.

CHEN Bingrui, FENG Xiating, HUANG Shuling, et al. Inversion of viscoelasto-plastic parameters based on fast Lagrangian analysis of continuum-parallel particle swarm algorithm and its application[J]. Chinese Journal of Rock Mechanics and Engineering, 2007, 26(12): 2517-2525.

[12] 李志敬, 朱珍德, 周伟华. 基于CPSO算法的岩石蠕变模型非定常参数反演分析[J]. 河海大学学报: 自然科学版, 2008, 36(3): 346-349.

LI Zhijing, ZHU Zhende, ZHOU Weihua. Back analysis of non-stationary parameters of rock creep model based on chaos particle swarm optimization algorithm[J]. Journal of Hohai University: Natural Sciences, 2008, 36(3): 346-349.

[13] 罗润林, 阮怀宁, 朱昌星. 基于粒子群-最小二乘法的岩石流变模型参数反演[J]. 辽宁工程技术大学学报: 自然科学版, 2009, 28(5): 750-753.

RUO Ruilin, RUAN Huaining, ZHU Changxing. Parameter inversion of rock creep model based on PSO-least square method[J]. Journal of Liaoning Technical University: Natural Science, 2009, 28(5): 750-753.

[14] 刘文彬, 刘保国, 刘中战, 等. 基于改进PSO算法的岩石蠕变模型参数辨识[J]. 北京交通大学学报, 2009, 33(4): 140-143.

LIU Wenbin, LIU Baoguo, LIU Zhongzhan, et al. Parameter identification of creep constitutive model of rock based on modified PSO algorithm[J]. Journal of Beijing Jiaotong University, 2009, 33(4): 140-143.

[15] 王伟. 改进粒子群优化算法在边坡工程力学参数反演中的应用[D]. 南京: 河海大学岩土工程研究所, 2007: 17-21.

WANG Wei. The application of an improved particle swarm optimization in inversion of mechanical parameter of slope engineering[D]. Nanjing: Hohai University. Geotechnical Research Institute, 2007: 17-21.

[16] 杜好. 基于微粒群算法的堆石坝坝料参数反演分析[D]. 大连: 大连理工大学土木工程学院, 2006: 21-26.

DU Hao. Back analysis of rock-fill material parameters using particle swarm optimization[D]. Dalian: Dalian University of Technology. School of Civil Engineering, 2006: 21-26.

[17] 高玮. 基于粒子群优化的岩土工程反分析研究[J]. 岩土力学, 2006, 27(5): 795-798.

GAO Wei. Back analysis algorithm in geotechnical engineering based on particle swarm optimization[J]. Rock and Soil Mechanics, 2006, 27(5): 795-798.

[18] Kennedy J, Eberhart R. Particle swarm optimization[C]//Proc IEEE International Conference on Neural Networks. Piscataway, NJ: IEEE Service Center, 1995: 1942-1948.

[19] Eberhart R, Kennedy J. A new optimizer using particle swarm theory[C]//Proc on 6th International Symposium on Micromachine and Human Science. Piscataway, NJ: IEEE Service Center, 1995: 39-43.

[20] Croes G A. A method for solving traveling-salesman problems[J]. Operations Research, 1958, 6(6): 791-812.

[21] 杨文东. 复杂高坝坝区边坡岩体的非线性损伤流变力学模型及其工程应用[D]. 济南: 山东大学土建与水利学院, 2011: 179-185.

YANG Wendong. Nonlinear damage rheological mechanical model for complex rock of high dam slope and its engineering application[D]. Jinan: Shandong University. School of Civil Engineering, 2011: 179-185.

[22] 熊诗湖, 周火明, 钟作武. 岩体载荷蠕变试验方法研究[J]. 岩石力学与工程学报, 2009, 28(10): 2121-2127.

XIONG Shihu, ZHOU Huoming, ZHONG Zuowu. Study of methodology of plate-loading creep test of rock mass[J]. Chinese Journal of Rock Mechanics and Engineering, 2009, 28(10): 2121-2127.

(编辑  陈灿华)

收稿日期：2011-11-01；修回日期：2012-01-10

基金项目：国家自然科学基金资助项目(51279097，51109123)；中央高校基本科研业务费专项资金资助项目(11CX04051A)

通信作者：杨文东(1982-)，男，山东潍坊人，博士，讲师，从事岩石力学与工程的教学和研究工作；电话：15376713557；E-mail: wendongy@gmail.com