粒子群算法改进及内变量本构模型参数反演_袁克阔.pdf

第 45 卷 第 2 期 煤田地质与勘探 Vol. 45 No.2 2017 年 4 月 COAL GEOLOGY constitutive model; particle swarm optimization algorithm; variable inertia weight; variable learning factor 由于岩土介质的非均质、非线性等原因，使得 合理本构模型的提出和对应参数的正确给定成为数 值方法解决岩土工程问题最为重要和最为迫切的两 大难题，特别是对于各内变量模型，参数的确定更 是难以从测取数据得出。随着现代量测技术[1]和计 算能力的发展，反分析[2]为解决岩土本构模型及力 学参数的确定这两大影响岩土学科发展的瓶颈问题 打开了思路。 当前岩土工程的反分析主要基于数值方法展 开，即给定待求参数 Xx1,x2,,xn，通过数值计算 得到监测指标计算值 Uc  T ccc 12 ,,, n uuu ， 将其与量 测信息 Ut T ttt 12 ,,, n u uu 建立目标函数 X。其本 质是建立计算指标与量测指标的极值函数，通过调 节输入参数而求解计算指标，直至目标函数满足条 件，即将难以直接量测参数的反演问题转化为一个 易量测信息目标函数在约束条件下的寻优问题。量 测信息最为常见的是位移、应力、渗透系数等参数； 根据采用的输入信息的不同，可以命名相应的反演 方法，如以位移为输入参数的位移反演法，以渗透 ChaoXing 第 2 期 袁克阔 粒子群算法改进及内变量本构模型参数反演 113 系数、压力水头、流量等渗透性参数为输入参数的 渗流反演法等。关于反问题解的适定性、反问题的 求解方法等内容，已有相当的研究[3]。 反演理论已经由纯数值研究阶段发展到了智 能、系统的新型反演算法阶段；最初的纯数值求解 主要应用传统确定性优化技术，以单纯形法、 Nelder-Mead 法等为代表； 智能反演法采用随机性优 化算法，以模拟退火算法、遗传算法、人工神经网 络等为代表。由于岩土工程的复杂性及传统方法自 身的不足，使得应用传统优化算法进行参数的反演 计算难以得到满意的解答[4-5]。在众多的智能反演算 法中，粒子群算法以其概念简单、易于实现、算法 参数少、计算量小、收敛快并有深刻的仿生背景， 而且不需先验学习和训练、不需二进制运算或微分 导数操作、无需交叉变异、不需极度苛刻的限制条 件等优点成为反分析的首要选择[6]， 贾善坡等[7]基于 常规粒子群算法与混合罚函数法进行了反演优化算 法研究，编制了有限元优化反演分析程序，并完成 了盾构隧道盾尾空隙注浆充填等代层参数的反演计 算。然而该算法仍存在对于有多个局部极值点的函 数，容易陷入到局部极值点中等不足，需要进行改 进。苏国韶等[8]提出一种基于粒子群变惯性权重优 化与高斯过程的智能协同优化的反演新方法。李金 凤等[9]提出混沌直接搜索粒子群CHPSO-DS算法 进行了面板坝堆石体力学参数的反演。 本文在不引入新概念的条件下，通过基于自然 选择的自适应变惯性权重和异步变化学习因子的方 法改进粒子群算法，并通过 Matlab 软件予以程序实 现，反分析岩土材料内变量模型材料参数，显现了 该方法的先进。 1 粒子群算法基本原理 粒子群算法particle swarm optimization，PSO 是由 J Kennedy 和 R Eberhart 于 1995 年提出的一种 源于对鸟群捕食行为研究的群体智能迭代随机优化 算法[10-11]；该算法采用速度–位置搜索模型，即每个 粒子 Xi代表 n 维解空间内的一个候选解，其由速度 Vi 12 ,, iiim vvv（，和位置 Xi 12 ,,, iiim xxx用于本 身状态的调整，n 称为粒子的长度或者是群体规模， 表示待优化/待反演的参数个数，m 称为种群大小， 表示自变量的个数，取决于输入的标准值/量测值个 数量测点的个数。第i个粒子当前时刻搜索到的最 优位置称为个体极值，记为 12 ,, iii ppp, im p表 示目标函数在第i个粒子当前位置处的函数值；整 个粒子群当前时刻搜索到的最优位置为全局极值 表示位置，记为 12 ,,, ggggm pppp。开始时， 随机初始化粒子的位置和速度构成初始种群，初始 种群在解空间中为均匀分布；后续迭代计算中在找 到这两个最优值时，粒子根据如下公式来更新自己 的速度和位置[12]  1 1 12 2 kkkkkk ijijijijgjij vvc rpxc rpx   1 11kkk ijijij xxv   2 式中 1,2,,in、1,2,,jm；k为当前迭代次数； 1 c 和 2 c 为学习因子，也称加速常数acceleration constant，为非负常数；根据经验，通常 12 2cc。 1 r 和 2 r 为[0,1]范围内的均匀随机数[13]。 式1右边由 3 部分组成，第一部分为“惯性 inertia”或“动量momentum”部分，反映了粒子的 运动“习惯”， 代表粒子有维持自己先前速度的趋势； 第二部分为“认知cognition”部分， 反映了粒子对自 身历史经验的记忆或回忆，代表粒子有向自身历史 最佳位置逼近的趋势；第三部分为“社会social”部 分，反映了粒子间协同合作与信息共享的群体历史 经验，代表粒子有向群体或邻域历史最佳位置逼近 的趋势。粒子的速度 maxmax [,] ij vvv ， max v是常数， 由用户设定用来限制粒子的速度。有研究[14]给式1 第一项 k ij v添加了一个惯性权重参数，变为  1 1 12 2 kkkkkk ijijijijgjij vvc rpxc rpx   3 用以确定局部搜索能力和全局搜索能力的比 例， k ij v表示在多大程度上保留继承原来的速度。 较大， 全局收敛能力强， 局部收敛能力弱；较 小，局部收敛能力强，全局收敛能力弱。 2 粒子群算法的改进 优化算法的改进工作，无非基于寻优能力和寻 优速度两个方面。首先，从算法的进化方程看，标 准 PSO 算法依靠的是群体之间的合作和竞争，粒子 本身没有变异机制，因而单个粒子一旦受某个局部 极值约束后本身很难跳出局部极值的约束，此时需 要借助其他粒子的成功发现来跳出局部极值。变权 重和异步变化学习因子的 PSO较好地解决了全局与 局部、自身经验与社会经验的协调。其次，通过分 析算法的进化方程式1和式2发现， 粒子们在搜索 时，总是追逐当前全局最优点和自己迄今搜索到的 最优点，因此粒子们的速度很快降到接近于 0，导 致粒子们陷入局部极小而无法摆脱，这种现象被称 为粒子群的“趋同性”。这种“趋同性”限制了粒子的 搜索范围。要想扩大搜索范围，就要增加粒子群的 粒子数或者减弱粒子对整个粒子群当前搜索到的全 局最优点的追逐。增加粒子数将导致算法计算复杂 ChaoXing 114 煤田地质与勘探 第 45 卷 度增高，而且研究表明通过增加粒子数所获得的效 果并不明显，而减弱粒子对全局最优点的追逐又存 在算法不易收敛的缺点。采用自然选择策略来保持 粒子群体的多样性，解决算法的早熟问题，提高算 法的探索能力是一项较好的选择。 故为了平衡PSO算法的全局搜索能力和局部改 良能力，采用自适应惯性权重法，以实现惯性权重 系数是非线性动态变化，公式如下 maxmin minminavg avgmin maxavg , , ffff ff ff             ≤ 4 式中 maxmin 、分别为惯性权重 的最大值和最 小值； 对于大多数问题， 其取值范围为0.9, 0.2；f 表示粒子当前的目标函数值， avg f表示当前所有粒 子的目标函数平均值， min f为目标函数最小值。该 式确定的惯性权重随粒子的目标函数值而自动变 化，因此称之为自适应权重；该权重变化策略简单 易于计算且不与迭代步数直接关联，具有一定灵活 性。由式4可见，当各粒子的目标值趋于一致或者 趋于局部最优时，将使得惯性权重增加，而各粒子 的目标函数值比较分散时惯性权重减小；另外，目 标函数值优于平均值时惯性权重因子保持较小值， 从而保证了此类粒子进行更好的局部搜索； 相反地， 目标函数值差于平均值时则保持较大的惯性权重， 以进行全局搜索。 Asanga 等[15]研究表明， 随时间变化的学习因子 同样对算法的性能也有很大的影响；在优化的初级 阶段，粒子具有较大的自我学习能力和较小的社会 学习能力，进行着较强的全局搜索，而在优化的后 期，粒子则具有较大社会学习能力和较小的自我学 习能力，有利于收敛到全局最优解。从而可见，两 个学习因子在优化过程随时间进行不同而变化的异 步变化学习因子能够满足此类优化要求，采用如下 学习因子异步更新公式 1,fin1,ini 11,iniiter max 2,fin2, 22,iniiter max ini cc ccT T cc ccT T           5 式中 1,fin c和 1,ini c分别是 1 c 的最终值和初始值，2,finc 和 2,ini c分别是 2 c 的最终值和初始值，通常情况下， 1,fin 0.5c、 1,ini 2.5c， 2,fin 2.5c、 2,ini 0.5c； Titer 为当前迭代步数， max T为设定的最大迭代步数。 通过上述自适应变权重系数和异步变学习因子 的方法提高算法的寻优能力仍不够，将遗传算法中 的自然选择机理应用到粒子群算法中，得到基于自 然选择、自适应变惯性权重、异步变化学习因子的 改进粒子群算法，可非常好地保持粒子群体的多样 性，解决算法的早熟问题，提高算法的探索能力。 其基本思想为每次迭代过程中将整个粒子群按适应 度排序，用群体中最好的一半粒子的速度和位置替 换最差的一半的位置和速度，同时保留原来每个个 体所记忆的历史最优值。 3 改进算法的程序实现 基于自然选择、自适应变惯性权重、异步变化 学习因子的改进粒子群算法的程序实现流程如下 ①初始化一个规模为 N、维数为 M 的粒子群， 设定初始位置和速度。 ②计算每个粒子的适应值。 ③对每个粒子将其适应值和其经历过的最好位 置 pbest的适应值进行比较，若较好，则将其作为当 前的全局最好位置。 ④对每个粒子将其适应值和全局经历过的最好 位置 gbest的适应值进行比较，若较好，则将其作为 当前的全局最优位置。 ⑤更新每个粒子的速度和位置 按式2和式3进行粒子速度和位置的更新， 其 中权重按照式4来计算，学习因子 1 c 、 2 c 按式5 计算。 ⑥将整个粒子群按适应度排序，用群体中最好 的一半粒子的速度和位置替换最差的一半的位置和 速度，而保持 pbest和 gbest不变。 ⑦如果满足终止条件，则输出解，否则返回第 ②步。 4 改进算法测试 为了测试改进粒子群优化算法的性能，采用多个 经典基准测试函数[16]进行了计算和试验。这里给出 2 个具有代表性的基准函数的测试结果，Sphere 函数为 较为简单的单峰函数，在0 i x 处达到极小值 0； Rastrigrin 函数为复杂的多峰函数， 当0 i x 时达到全局 极小值 0，在其周围存在多个局部极小点。两测试函数 基本特性见表 1，两函数 2 维情形下的图像如图 1，两 函数在算法执行过程中各参数的设置如表 2。 两函数进化过程中适应度变化曲线如图 2 所示。 由图 2、表 3 可见，不论对于单峰函数还是多 峰函数，改进算法在寻优速度和寻优能力上均有巨 大的改进，效率实现了百倍的提高。因此，可以将 改进算法应用于岩土材料参数的反演分析中。 ChaoXing 第 2 期 袁克阔 粒子群算法改进及内变量本构模型参数反演 115 表 1 标准测试函数 Table 1 Standard test functions 函数 Sphere Rastrigrin 表达式 2 1 1 n i i f xx   2 2 1 10cos2π 10 n ii i fxxx    峰值特点 单峰 多峰 优化维数 30 30 变量范围 [ 100,100]n [ 5.12,5.12]n 最优值 0,,00f 0,,00f 图 1 测试函数图像n2 Fig.1 Functional images of the test 表 2 算法参数设置 Table 2 Parameter setting of PSO algorithm 参数 PSO 改进 PSO 惯性权重系数 0.9 max 0.9， min 0.2 学习因子 12 2.0cc 1,fin1,ini 0.52.5cc、 ， 2,fin2,ini 2.50.5cc、 最大迭代步数 1 000 1 000 图 2 适应度随进化代数变化情况 Fig.2 Variation of fitness with evolving algebra 函数最优值如表 3。 表 3 函数测试结果 Table 3 The test results 函数 PSO 改进 PSO Sphere 0.072 964 173 889 82 0.001 302 752 805 81 Rastrigrin 2.884 650 063 705 86310-5 0 5 基于改进算法的反演程序 基于量测数据进行岩土本构参数确定过程中， 选择目标函数 X在约束条件下求解 min X是核 心问题。本文采用目标函数表达式为  x1, x2,, xn 2 c t 1 1 n i i i u u       6 式中，直接自变量为 Uc，间接自变量 X 为内变量本 构模型中待反演材料参数，Uc UX。在大型有限 元程序 ABAQUS 的基础上， 以 Matlab 语言为平台， 结 合 改 进 粒 子 群 算 法 ，编 制 优 化 反 分 析程 序 GeoPSOInverse.m 以实现算法模块与有限元软件的 相耦合。程序流程框图如图 3。 反演程序实现步骤如下 ①依据具体问题建立有限元模型； ②通过 system 命令调用 ABAQUS 进行计算； ③读取 ABAQUS 结果文件， 获取计算所得 uc值； ④应用改进粒子群算法，对目标函数进行优化 计算，更新参数 X； ChaoXing 116 煤田地质与勘探 第 45 卷 ⑤读取有限元计算文件，赋予待反演参数更新 值，进行新一轮的计算。 图 3 岩土体参数反演程序框图 Fig.3 Flow chart of inverse analysis program 6 泥岩内变量蠕变模型参数反演 深埋煤矿巷道的失稳和破坏往往不是在开挖完 成后就立即发生，而是要经历应力与变形不断调整 的较长时期，伴随产生的底臌、两帮收敛可达米级， 并引起混凝土衬砌大面积开裂， U 型钢扭曲等现象， 严重影响构筑物的安全。因而在深埋高压软岩地下 工程的设计时， 不仅要考虑岩体的瞬时弹塑性变形， 而且还要高度重视其随时间发展的蠕变变形。 作者基于上下加载面修正剑桥模型建立了以不 可恢复应变为内变量的不显含时间项的富水泥岩蠕 变模型[17]，并完成了数值实现；其需反演辨识的参 数为 c A 、 c m 、 c n 3 个未知参数。 为考察所改进算法及反演程序的可行性，分别 对富水泥岩 4 800 kPa 和 6 400 kPa 应力水平下的单 轴蠕变和 1.52.5 MPa 应力水平下的三轴蠕变试验 进行反演，最后应用反演所得参数拟合两类试验最 后应力水平下的试验结果。对于单轴蠕变，有限元 模型底部与径向施加法向约束，三轴蠕变给定定常 围压 2233 ，试样底端轴向约束，径向施加接地 弹簧单元 SPRING；然后施加轴向载荷，据蠕变本 构计算轴向应变， 反分析计算蠕变过程中应变函数  y ccc ,,, ijij xA m n 7 反分析过程中采用应力加载方式，对于每一个 给定的蠕变时刻，数值计算结果和试验结果中的应 变存在一定的差值，以应变差值的绝对值趋向于最 小值为目标函数，即  x1, x2,, xn 2 c t 1 1min n i i i u u       8 式中 y x 为试验测得的已知材料参数， t i u为荷载作 用下不同蠕变时间下的轴向变形试验值， c i u是有限 元计算值。 反演过程与有限元计算过程相互作用，有限元 软件提供计算值 c i u，反演程序根据优化算法提供更 新后的有限元所用材料参数。反演过程是有约束条 件的，据经验和先期试算给定参数范围如表 4，最 终反演结果如表 5。 表 4 待反分析参数范围 Table 4 The parameter range for inversing 参数 Ac/10-21 nc mc 范围 1.010.0 2.010.0 –6.0–1.0 表 5 反分析计算结果 Table 5 The inverse results 参数 Ac/10-21 nc mc 固结蠕变 3.13 8.82 –1.76 反演值 三轴蠕变 1.65 8.46 –2.32 由反演结果可见， 不同蠕变试验单轴与三轴 下反演所得材料参数十分相近；蠕变曲线试验与 数值计算结果拟合程度很高图 4， 特别是应用两 图 4 数值结果与试验结果对比 Fig.4 Comparison of creep curves between simulated and experimental results ChaoXing 第 2 期 袁克阔 粒子群算法改进及内变量本构模型参数反演 117 类试验的部分试验数据反演获取的参数，分别 数值模拟 8 000 kPa 和 3.0 MPa 应力水平下的蠕 变特征，结果显示模拟结果与试验结果吻合度 非常高，证明了蠕变模型的正确性和反演程序 的有效性。 7 结 语 首先通过分析岩土工程反演现状分析，指出 了应用粒子群算法进行反分析的优越性和算法的 改进策略。其次基于自然选择、自适应变惯性权 重、异步变化学习因子的策略改进了粒子群算法 并进行了程序实现，通过 Sphere 和 Rastrigrin 两 函数的测试证实，改进算法无论寻优能力还是寻 优速度都大为增强。最后以 Matlab 语言为平台， 联合大型有限元程序 ABAQUS，编制优化反分析 程序 GeoPSOInverse.m；应用所编程序反演了以 不可恢复应变为变量的、不显含时间的泥岩新型 蠕变模型参数。结果表明，本文改进的粒子群算 法在岩土工程参数反演计算中体现出了可靠的反 演能力和很快的收敛速度。 参考文献 [1] GIODA G，SAKURAI S．Back analysis procedures for the interpretation of field measurements in geomechanics[J]. Int J Numerical Analytical s Geomechanics，1987，11 555–583． [2] 杨林德．岩土工程问题的反演理论与工程实践[M]．北京科 学出版社，1999. [3] 李守巨，刘迎曦，孙伟. 智能计算与参数反演[M]. 科学出版 社，200848–65. [4] 冯夏庭，张治强，杨成祥，等. 位移反分析的进化神经网络方 法研究[J]. 岩石力学与工程学报，1999，185529–533. FENG Xiating，ZHANG Zhiqiang，YANG Chengxiang，et al. Study on genetic-neural network of displacement back analysis[J]. Chinese Journal of Rock Mechanics and Engineering，1999，185529–533. [5] 高玮， 郑颖人. 一种新的岩土工程进化反分析算法[J]. 岩石力 学与工程学报，2003，222192–196. GAO Wei，ZHENG Yingren. New evolutionary back analysis algorithm in geotechnical engineering[J]. Chinese Journal of Rock Mechanics and Engineering，2003，222192–196. [6] 李宁. 粒子群优化算法的理论分析与应用研究[D]. 武汉华 中科技大学，2006. [7] 贾善坡， 伍国军， 陈卫忠. 基于粒子群算法与混合罚函数法的 有限元优化反演模型及应用[J]. 岩土力学， 2011， 32增刊 2 598–603. JIA Shanpo， WU Guojun， CHEN Weizhong. Application of finite element inverse model based on improved particle swarm optimization and mixed penalty function[J]. Rock and soil mechanics，2011，32 S2598–603. [8] 苏国韶， 张克实， 吕海波. 位移反分析的粒子群优化高斯过 程协同优化方法[J]. 岩土力学，2011，32增刊 2510–524. SU Guoshao， ZHANG Keshi, LYU Haibo. A cooperative optimization based on particle swarm optimization and gaussian process for displacement back analysis[J]. Rock and Soil Mechanics，2011，32 S2510–524. [9] 李金凤，杨启贵，徐卫亚. 基于改进粒子群算法 CHPSO-DS 的面板坝堆石体力学参数反演[J]. 岩石力学与工程学报， 2008，2761229–1235. LI Jinfeng， YANG Qigui， XU Weiya. Back analyzing mechanical parameters of rockfill based on modified particle swarm optimization CHPSO-DS[J]. Chinese Journal of Rock Mechanics and Engineering，2008，2761229–1235. [10] KENNEDY J，EBERHART R. Particle swarm optimization[C]//Proceedings of the IEEE International Conference on Neural Networks，1995，41942–1948. [11] EBERHART R C，KENNEDY J. A new optimizer using particle swarm theory[C]//Proceedings of the Sixth International Symposium on Micro Machine and Human Science，MHS95， Nagoya，Japan，199539–43. [12] LANGDON W B，POLI R. Evolving problems to learn about particle swarm and other optimizers[C]//2005 IEEE Congress on Evolutionary Computation，2005，181–88. [13] 谢晓锋，张文俊，杨之廉. 微粒群算法综述[J]. 控制与决策， 2003，182129–134. XIE Xiaofeng，ZHANG Wenjun, YANG Zhilian. Overview of particle swarm optimization[J]. Control and Decision，2003， 182129–134. [14] SHI Y ， EBERHART R C. A modified particle swarm optimizer[C]//Proceedings of the IEEE International Conference on Evolutionary Computation，199869–73. [15] ASANGA R，SAMNA K H. Self-organizing hierarchical particle swarm optimizer with time-varying acceleration coefficients[J]. IEEE Transactions on Evolutionary Computation，2004，38 240–255. [16] 张顶学， 关治洪， 刘新芝. 一种动态改变惯性权重的自适应粒 子群算法[J]. 控制与决策，2008，2311 1253–1257. ZHANG Dingxue，GUAN Zhihong，LIU Xinzhi. Adaptive particle swarm optimization algorithm with dynamically changing inertia weight[J]. Control and Decision， 2008， 2311 1253–1257. [17] 陈卫忠， 袁克阔， 于洪丹， 等. Boom Clay 蠕变特性研究[J]. 岩 石力学与工程学报，2013，32101981–1990. CHEN Weizhong，YUAN Kekuo，YU Hongdan，et al. Creep behavior of boom clay[J]. Chinese Journal of Rock Mechanics and Engineering，2013，32101981–1990. 责任编辑 张宏 ChaoXing