基于改进粒子群优化算法的地应力场反演.pdf
基于改进粒子群优化算法的地应力场反演 ① 李 欣1, 谢学斌1, 杨承祥2, 罗海霞1 (1.中南大学 资源与安全工程学院,湖南 长沙 410083;2.安徽铜冠(庐江)矿业有限公司,安徽 铜陵 231561) 摘 要 为提升粒子群优化算法反演地应力的效率和精度,采用改进粒子群算法,并对随机数进行针对性设置,联合使用 MATLAB 与 FLAC3D,以沙溪铜矿区的声发射数据为基础进行了算例分析。 结果表明,操作流程简单易行,改进后的算法精度高、效率高,对 类似工程有一定的参考性。 关键词 初始应力场; 地应力反演; 改进粒子群算法; 算例分析 中图分类号 TU443文献标识码 Adoi10.3969/ j.issn.0253-6099.2015.01.007 文章编号 0253-6099(2015)01-0023-04 Inversion Analysis of Initial Stress Field Based on Modified Particle Swarm Optimization LI Xin1, XIE Xue⁃bin1, YANG Cheng⁃xiang2, LUO Hai⁃xia1 (1. School of Resources and Safety Engineering, Central South University, Changsha 410083, Hunan, China; 2. Anhui Tongguan (Lujiang) Mining Co Ltd, Tongling 231561, Anhui, China) Abstract In order to enhance the efficiency and accuracy of the inversion analysis of initial stress field, the analysis of calculation examples was performed based on the acoustic emission data of Shaxi Copper Mine, by adopting the modified particle swarm optimization(MPSO) with a special set of random numbers, together with MATLAB and FLAC3D. Result shows that the MPSO is characterized by high accuracy and efficiency, as well as simple operation process. It is of certain reference value to the similar engineering. Key words initial stress field; inversion analysis of initial stress field; modified particle swarm optimization(PSO); analysis of calculation examples 初始地应力场是引起采矿、土木建筑、水利水电等 各种岩土体工程变形及破坏的重要因素,也是地下工 程设计施工的初始条件[1-5]。 由于成本和工作环境的 影响,仅能对岩体部分位置进行测量,要了解整个工区 的地应力场分布情况,必须运用合理的数值模型和优 化算法来进行地应力的反演。 力学参数反演问题为最 小化问题,本文认为适应度函数值越小的粒子对应的 适应度越好,所在的位置越好,所代表的一组反演参数 越优。 一般来说,多元回归分析方法、遗传算法和神经 网络方法在反演地应力时,都需要大量的正分析计算。 与其它反演方法相比,粒子群算法对较小的群体个数 也能充分搜索解空间,可有效解决工程问题,同时避免 过多的适应度估计,在地应力反演等其它岩土力学工 程中应用越来越多[6-14]。 基本粒子群算法在搜索后期 会出现粒子在全局最优解附近“震荡”的现象,影响了 反演的效率和精度,为了避免这个问题,可在迭代过程 中线性的减小惯性权重[15]。 另外,由于粒子群算法初 始样本较少,而粒子越界经常发生,为减少粒子的损 失,避免算法进入局部最优解,本文在算法中设置隐匿 墙[14],并对随机数进行针对性值域限定,保证算法收 敛的稳定性,提出了一种基于改进粒子群优化算法的 地应力反演方法。 将 MATLAB 与 FLAC3D相结合,运 用改进后的粒子群算法对地应力进行反演,工程实例 表明,这种方法简单易行且精度较高。 1 粒子群优化算法 粒子群优化是以邻域原理为基础进行操作的,群 体中的个体相互学习,而且基于获得的信息移动到更 相似于它们的较好的邻近区域[16]。 每次迭代过程不 是完全随机的,如果找到较好的解,将会以此为依据来 寻找下一个解。 其实质是所有粒子在确定的搜索空间 内,不停的追踪全局最优解和个体最优解,以获得满足 ①收稿日期 2014-09-02 基金项目 安徽省科技攻关计划项目资助(12010402148) 作者简介 李 欣(1987-),男,山东人,硕士研究生,主要从事岩土工程相关研究和工作。 第 35 卷第 1 期 2015 年 02 月 矿矿 冶冶 工工 程程 MINING AND METALLURGICAL ENGINEERING Vol.35 №1 February 2015 要求的最优解。 1.1 基本粒子群优化算法 系统初始化一组随机粒子,每个粒子都有其随机 速度在解空间中搜索,并且粒子具有记忆功能,在每一 次迭代过程中,粒子都向两个“最优值”靠近一个是 个体最优解,一个是全局极值。 一个由 m 个粒子组成 的群体在 D 维搜索空间中,在某时刻,第 i 个粒子的位 置表示为 Xi,第 i 个粒子的速度表示为 Vi,第 i 个粒子 经过的历史最好点表示为 Pi,群体内所有粒子所经过 的最好点表示为 Pg,如式(1)所示 Xi =(x i1,xi2,,xid,,xiD) Vi =(v i1,vi2,,vid,,viD) Pi =(p i1,pi2,,pid,,piD) Pg =(p g1,pg2,,pgd,,pgD) (1 ≤ i ≤ m, 1 ≤ d ≤ D) (1) 式中 xid为粒子某项特性的位置;vid为其变化趋势。 为 了改善算法收敛性能,Shi 和 Erhart 在 1998 年引入了 惯性权重的概念[17],将速度和位置更新,方程改为 vk+1 id = ωvk id + c 1ξ(p k id - x k id) + c2η(p k gd - x k id) xk+1 id = x k id + γv k+1 id { (2) 式中 c1和 c2称为学习因子,又称为加速系数,它们使粒 子具有自我总结和向群体中优秀个体学习的能力,从而 向自己的历史最优点以及群体内或邻域内的历史最优 点靠近。 c1与 c2之和最好接近 4.0[18],可取 c1 =c 2 = 2。 η 和 ξ 是[0,1]区间内的随机数。 粒子每一维的变化 量都会限制一个最大值 vmax,如果某维更新后的速度 超过了预先设定的 vmax,那么该维的速度就被限定为 vmax。 式中 ω 称为惯性权重,其大小决定了对粒子当 前速度继承的多少,合适的选择可使粒子具有均衡的 探索和开发能力。 1.2 改进粒子群优化算法 为提升粒子群算法的性能,目前对粒子群所做的 改进主要集中在粒子越界的处理方式和惯性权重的设 置上。 本文在考虑以上两方面的前提下,又对随机数 进行了针对性设置。 1.2.1 粒子越界的处理 粒子越界是 PSO 算法不可 避免的问题,在粒子群算法早期进化阶段,粒子越界的 数量和范围是非常可观的。 由于粒子群算法所采取的 样本数量较少,而粒子越界会造成相当大比例的粒子 流失,从而造成算法最终结果陷入局部最优解,而无法 达到精度要求。 目前有 4 种常规的方法来解决粒子撞 击边界的问题吸收墙、反射墙、衰减墙和隐匿墙,如 图 1 所示。 图 1 4 种规则算法 (a) 吸收墙;(b) 反射墙;(c) 衰减墙;(d) 隐匿墙 吸收墙将越界粒子速度置为 0,粒子能量被吸收; 反射墙将粒子速度置反,大小不变。 这两种方法都能 使粒子无法穿越边界。 衰减墙结合了吸收墙和反射墙 的特性,将越界粒子速度取反,随机的减少速度大小。 隐匿墙粒子飞出搜索空间,粒子的位置和速度不人为 的改变,保证其随机性,但该粒子不参与极值更新,在 迭代过程中重返解空间。 前 3 种方法均是将粒子置于边界上或某一范围内, 算法简易可行,但会造成多个越界粒子聚集在边界上, 如果边界或附近存在局部最优解,那么这些粒子容易陷 入局部最优点,从而导致算法早熟收敛。 隐匿墙操作既 节省了工作量,提高了反演效率,又不影响粒子群的自 然移动,在一定程度上削弱了粒子越界对算法进化的影 响,又不易造成局部最优解。 根据粒子速度更新公式可 知,在粒子群算法中,指导粒子的运动方向与步长的是 粒子个体最优位置和群体最优位置,这种运动引导机 制,保证了最优位置对各粒子的引力。 因此只要保证最 优位置在约束空间内,就可以把约束空间外的粒子通 过速度更新和迭代吸引到约束空间里面来[14]。 1.2.2 惯性权重的设置 为了避免基本粒子群算法 在解空间内后期搜索时粒子在全局最优解附近“震 荡”的现象,可随着迭代进程线性的减小惯性权重,如 式(3)所示 ω = ωmax- ωmax - ω min Nmax N(3) 其中 N 为当前迭代步数;Nmax为最大迭代步数;ωmax为 最大惯性权重;ωmin为最小惯性权重。 这种方法要使 用合理的迭代步数,如果迭代步数过少,算法即使搜索 次数不足也会过早进入收敛陷入早熟;如果迭代步数 过多,则要等惯性权重变小后才会搜索极值,造成等待 时间过长。 42矿 冶 工 程第 35 卷 1.2.3 随机数的值域 一般情况下,η 和 ξ 在[0,1]区 间内随机取值。 在迭代前期,惯性权重 ω 较大,各个 粒子速度更新值也较大,这样会频繁造成粒子超速和 越界问题,如图 2 所示。 如果此时随机数在较小值域 中取值,那么不仅会使各个粒子平稳过渡,还能保证粒 子飞行的随机性、自主性。 而算法后期,惯性权重 ω 变小,各个粒子当前位置与全局最优粒子和各自历史 最优粒子间的区别也越来越小,如果此时 η、ξ 取值较 小,算法可能会提前进入收敛。 因此,为了保证粒子群 的搜索空间和后期的速度有效更新,在迭代前期,η 和 ξ 取值范围在[0,0.1]区间,在迭代中后期,η 和 ξ 取值 范围在(0.1,1]区间。 图 2 错过全局最优解 2 岩体初始地应力反演应用 岩体应力反演应以现场实测资料为主要依据。 沙 溪铜(金)矿区由铜泉山矿段和凤台山矿段组成,在矿 山开采设计进行岩石力学研究阶段,利用声发射方法 对原岩应力进行了测量。 模型坐标轴 X 轴近东西方 向,模拟长度为 2 000 m;Y 轴与 X 轴垂直,近南北方 向,长度为 2 000 m;Z 轴正方向铅直向上,模拟范围为 0~1 300 m。 三维地质模型划分四面体和六面体等单 元共 18 579 个。 根据地质资料,不考虑山脉因素,利 用 FISH 函数生成伪随机表面。 模型采用 3 个地质分 层,范围从上至下分为三区第一区为石英闪长斑岩, 第二区为石英闪长玢岩,第三区为泥质粉砂岩。 三维 地质模型如图 3 所示。 图 3 初始应力场计算区域及网格划分 2.1 岩体反演参数选取 本文以实测原岩应力和地质勘探钻孔资料为基础 进行地应力反演。 利用 MATLAB 进行粒子群迭代,正 分析过程利用 FLAC3D。 依据 200 余块声发射试件得 出的沙溪铜矿原岩应力分布大致特征,选取应力边界 条件。 构造应力、自重应力是地应力反演所获得的主 要未知量,剪切应力与正应力相比小得多,因此本文不 考虑剪切应力的影响,将反演边界条件设置为坐标轴 X、Y、Z 三个方向正应力及其梯度。 具体反演参数为 南北方向应力 SYY,南北方向应力梯度 grad(SYY),东西 方向应力 SXX,东西方向应力梯度 grad(SXX),铅直方向 应力 SZZ,铅直方向应力梯度 grad(SYY)。 反演参考点 选取沙溪铜矿区原岩应力场声发射法测量研究报 告[19]中埋深分别为-320、-470、-555、-720、-740 m 的 5 个点。 地应力参照点选取坐标具体为S1(650, 650,470),S2(950,950,320),S3(1 100,1 100,740), S4(1 400,1 600,555),S5(1 600,1 900,720)。 本构模 型采用莫尔⁃库伦模型。 在确定地应力反演工区以及参演参数后,构造地 应力反演的目标函数为 minF = 1 NM{∑ N i = 1 ∑ M j = 1 [S′ij(P) - Sij] 2 Sij2} 1/ 2 (4) 其中 S′(p)为测点在模型参数为 P 时所得的地应力 值;Sij为第 j 个实测地应力值分量。 2.2 粒子群计算参数设定 式(4)中 F 作为粒子群算法中粒子的适应值函 数,来评价给定解位置的优劣。 采用设置隐匿墙和线 性递减惯性权重的改进粒子群算法,ωmax为 0.9,ωmin为 0.4。 为防止迭代计算无休止的进行下去,不以反演精 度而以算法迭代步数作为终止条件,本文反演迭代最 大次数 Nmax为 40。 反演终止条件以最小适应度小于 0.04 为准。 2.3 反演流程 粒子群算法反演地应力改进流程如图 4 所示。 在 这个计算流程中,粒子群计算过程是开放的,可以随时 查看算法状态。 如果发现程序进入早熟或者算法不收 敛,则返回第一步。 2.4 反演结果分析 经过反演计算,算法达到了设定次数 40 次,并且 最小适应度函数小于 0.04,取得了较好的反演结果。 绘制全局最优粒子适应函数值随迭代次数的进化曲线 图,如图 5 所示,其中 minF 为最优粒子适应度,N 为迭 代次数。 由图 5 可看出,粒子在早期迭代收敛速度较 快,后期趋于定值。 由于对随机数进行了设置,迭代前 期曲线较为平滑,说明算法较为稳定的进入收敛状态, 52第 1 期李 欣等 基于改进粒子群优化算法的地应力场反演 图 4 粒子群算法流程 图 5 最优粒子收敛过程 减少了更新随机粒子群重新进行计算的概率。 在第 25 步以后收敛速度减慢,最小适应函数的最优值保持 在 0.038 不变,且曲线基本为稳定状态,表明算法具有 良好的收敛性。 图 6 为粒子平均适应度函数值收敛过 程,其中 AVEF 为粒子平均适应度函数值,N 为迭代次 数。 据图 6 所示,当前粒子适应度函数平均值基本等 于最优粒子适应度函数值,这表明算法中粒子具有良 好的学习性。 图 6 粒子平均适应度函数值收敛过程 表 1 为实测、反演参数值及相对误差,其中 σ1、 σ2、σ3分别对应 3 个主应力。 由表可看出,各测点应 力分量值的模拟效果较为理想。 埋深-320 m 处第三 主应力相对误差较大,然而其绝对误差却只有 3.08 MPa,另外全局平均相对误差仅为 8.44%,所得的地应 力场模拟精度完全满足工程,所反演的应力边界条件 能较好的模拟构造运动。 表 1 测点地应力实测值、反演值和相对误差 埋深 / m σ1σ2σ3 实测 / MPa 反演 / MPa 相对误 差/ % 实测 / MPa 反演 / MPa 相对误 差/ % 实测 / MPa 反演 / MPa 相对误 差/ % -320 -20.92 -24.03 14.87 -17.08 -16.225.04-8.92 -12.0034.5 -470 -27.29 -23.27 14.73 -22.11 -21.542.58-13.97 -12.858.02 -555 -33.26 -33.120.42-23.45 -23.340.47 -16.85 -15.955.34 -720 -38.61 -36.714.92-27.71 -28.061.27 -20.38 -19.524.22 -740 -44.04 -37.23 15.46 -32.41 -28.91 10.80 -24.98 -24.003.92 3 结 语 1) 基本粒子群算法容易陷入局部最优解,改进后 的粒子群算法降低了粒子进入局部最优解的概率,提 高了收敛速度和反演精度。 2) 工程实例计算结果表明,改进粒子群优化算法 反演结果可靠,精度较高,收敛速度快,是一种很好的 初始地应力场反演方法。 3) 采用 MATLAB 进行粒子群矩阵迭代运算并在 FLAC3D中进行力学分析,简单易行,效率较高。 参考文献 [1] 蔡美峰,乔 兰,于 波,等. 金川二矿区深部地应力测量及其分 布规律研究[J]. 岩石力学与工程学报,1999,18(4)414-418. [2] 康红普,林 健,颜立新,等. 山西煤矿矿区井下地应力场分布特 征研究[J]. 地球物理学报,2009,52(7)1782-1792. [3] 李彦兴,董平川. 利用岩石的 Kaiser 效应测定储层地应力[J]. 岩 石力学与工程学报,2009,28(z1)2802-2807. [4] 罗超文,李海波,刘亚群,等. 深埋巷道地应力测量及围岩应力分 布特征研究[J]. 岩石力学与工程学报,2010,29(7)1418-1423. [5] 孟 文,陈群策,杜建军,等. 新加坡地应力测量[J]. 地球物理学 报,2012,55(8)2611-2619. [6] 张伟杰,兰思栋. 高应力下卸压巷道围岩破坏机理及卸压过程数 值分析[J]. 矿冶工程,2014(4)34-38. [7] 许梦国,王明旭,王 平,等. 不同巷道支护方式的数值模拟分析 [J]. 矿冶工程,2013(6)19-23. [8] 夏海燕,付建军,时 凯. 基于 BP 网络的深部软岩巷道围岩力学 参数反演研究[J]. 矿冶工程,2013(5)25-29. [9] 倪绍虎,肖 明,王继伟,等. 改进粒子群算法在地下工程反分析 中的运用[J]. 武汉大学学报(工学版),2009,42(3)326-330. [10] 陈炳瑞,冯夏庭,李庶林,等. 基于粒子群算法的岩体微震源分层 定位方法[J]. 岩石力学与工程学报,2009,28(4)740-749. [11] 杨文东,张强勇,李术才,等. 粒子群算法在时效变形参数反演中 的应用[J]. 中南大学学报(自然科学版),2013,44(1)282- 288. (下转第 30 页) 62矿 冶 工 程第 35 卷 表 1 不同采矿方法技术经济指标对比 采矿方法 生产能力 / (td -1 ) 工效 / [t(人班) -1 ] 贫化率 / % 损失率 / % 炸药单耗 / (gt -1 ) 采矿综合成本 / (元t -1 ) 全分段预裂挤压一次爆破法51277.85.18.90.18255.35 机械化盘区上向分层充填法24223.117.08.00.36349.00 对比分析表明,因采矿成本降低,企业由此获得的 经济效益为 1 490.15 万元,因矿石贫化率降低,企业由 此获得的经济效益为 10.55 万元,因采准切割工程减 少,企业由此获得的经济效益为 24.52 万元,三项合 计,试验采场新增经济效益为 1 525.22 万元。 4 采矿法适用条件与关键技术 4.1 适用条件 全分段预裂挤压一次爆破采矿法是一种高效采矿 方法,运用时要求① 矿体和围岩中等稳固以上;② 急 倾斜矿体中厚以上,缓倾斜和倾斜矿体为厚大矿体; ③ 无轨设备开采;④ 采空区尾砂胶结或尾砂充填。 4.2 采矿方法特征及关键技术 全分段预裂挤压一次爆破采矿法也是一种强化采 矿方法,具有 5 个特征① 安全的作业环境;② 单循 环的作业方式;③ 高效的采掘设备利用;④ 高强度的 回采强度;⑤ 较低的生产成本。 强化开采在两个方面 进行① 将分层开采变为分段开采;② 将多循环作业 变为单循环作业。 为此,必须解决好几个关键技术① 采场轮廓控 制技术。 采场轮廓涉及采场顶板、两帮、上盘和下盘, 主要控制技术为预裂爆破、光面爆破。 ② 高效支护技 术。 采场支护涉及采场顶板、两帮和上盘,主要支护技 术为预应力树脂锚杆+金属网+喷浆联合支护,支护工 作量大,需要使用锚杆台车和喷浆台车。 ③ 大规模爆 破技术。 全分段一次回采落矿,崩矿量大,炸药消耗量 大,为减少单耗药量,降低爆破振动危害,需要使用高 精度雷管和大规模爆破网路。 5 结 语 1) 新城金矿在不同时期、不同开采深度采用不同 的采矿方法,浅部开采采用分层充填采矿法,中深部开 采采用机械化盘区分层充填采矿法,深部开采采用全 分段预裂挤压一次爆破采矿法。 从分层充填采矿法转 变为机械化盘区分层充填采矿法,新城金矿的采矿方 法实现了机械化、无轨化,从机械化盘区分层充填采矿 法转变为全分段预裂挤压一次爆破采矿法,新城金矿 的采矿方法实现了无轨设备的配套化、高效化。 2) 全分段预裂挤压一次爆破采矿法通过全分段 一次凿岩、一次落矿爆破、集中出矿和集中充填成功实 现了深部矿体的高效高强度低成本开采,该方法值得 推广应用。 参考文献 [1] 刘益福,刘永会,蔚登峰,等. 新城金矿采矿技术现状与展望[J]. 黄金,2000,21(5)26-28. [2] 李启月,王树海,范作鹏,等. 盘区阶梯式无间柱连续充填采矿法 试验研究[J]. 矿冶工程,2010,30(3)16-19. [3] 王立君,葛富英,王剑波,等. 盘区上向高分层连续回采充填采矿 法试验研究[J]. 黄金,2002,23(1)17-21. [4] 赵国彦. 机械化盘区上向分层充填采矿法试验研究[J]. 金属矿 山,1999(6)1-4. [5] 王旭东,王中亮,彭永明,等. 新城金矿床Ⅰ号、Ⅴ号矿体特征对比 及找矿预测[J]. 中国科技信息,2012(17)42-44. [6] 王龙振,赵 海,刘国富,等. 新城金矿床Ⅰ号矿体深部探矿方案 探讨[J]. 有色金属(矿山部分),2002,54(1)17-19. (上接第 26 页) [12] 陈炳瑞,冯夏庭,黄书岭,等. 基于快速拉格朗日分析⁃并行粒子 群算法的黏弹塑性参数反演及其应用[J]. 岩石力学与工程学 报,2007,26(12)2517-2525. [13] 罗润林,阮怀宁,黄亚哲,等. 岩体初始地应力场的粒子群优化反 演及在 FLAC3D中的实现[J]. 长江科学院院报,2008,25(4)73 -76. [14] 田明俊,周 晶. 岩土工程参数反演的一种新方法[J]. 岩石力 学与工程学报,2005,24(9)1492-1496. [15] 张利彪,周春光,马 铭,等. 基于粒子群算法求解多目标优化问 题[J]. 计算机研究与发展,2004,41(7)1286-1291. [16] 蔡自兴,徐光祐. 人工智能及其应用[M]. 北京 清华大学出版 社,2010. [17] Shi Y H,Eberhart R. A modified particle swarm optimizer[C] ∥ IEEE Congress on Evolutionary Computation. Anchorage,AK,USA, 19989-73. [18] Kennedy. Thinking is social experriments with the adaptive culture model[J]. Journal of Conflict Resolution,1998,42(1)56-57. [19] 中南大学资源与安全工程学院,安徽铜冠(庐江)矿业有限公 司,铜陵有色金属集团控股有限公司技术中心采矿研究所. 沙溪 铜矿区原岩应力场声发射法测量研究[R]. 2011. 03矿 冶 工 程第 35 卷