基于萤火虫-支持向量机方法的地下厂房围岩力学参数反分析研究

2020-09-21王秀菊张金龙

河南科学 2020年8期

王秀菊，张金龙，戴薇，石崇

（1.南京交通职业技术学院，南京 211188； 2.河海大学岩土工程科学研究所，南京 210098；3.中国电建集团中南勘测设计研究院有限公司，长沙 410014）

破碎围岩地下洞室群工程中往往存在复杂的破碎带及不连续面，对工程的建设有重要影响，但岩体地下工程围岩具有非均匀性、非线性、非连续性特点，如何对围岩力学参数取值是进行变形预测与稳定分析的关键之一，也是困扰岩土工程领域的一大难题. 近年来，随着位移反分析方法的进步，采用智能反演算法结合数值模拟方法，应用计算机技术与岩土工程理论和工程实际相结合的方法来获取围岩变形参数，成为力学参数获取的重要途径［1-5］.

在围岩变形反分析方面，许传华等［6］结合索风营工程现场实际，建立了支持向量机和模拟退火位移反分析模型，并反分析其力学参数. 聂卫平等［7］综合适当的改进粒子群算法和支持向量机，进行糯扎渡水电站的围岩力学参数和初始地应力场反演，获得了良好的效果. 文辉辉等［8］利用BP神经网络和正交试验方法，反演了隧道围岩物理力学参数. 袁振华等［9］提出岩土力学参数GA算法与有限元耦合应力反分析模型，该模型采用应力反分析方法，克服了常规算法反演过程中易陷入局部极优的难题. 周喻等［10］基于PFC3D和BP神经网络方法建立计算模型并反演出岩土体细观力学参数. 陈志超等［11］结合地下洞室裂隙化花岗岩实例，建立了BP神经网络和DEC离散元方法相结合的计算模型，反演出了区域应力场值. 李澄清等［12］采用开源的颗粒离散元程序LMGC开展了土体双轴压缩数值试验，建立了BP人工神经网络反演系统，也获得了较为精确的变形参数.

本文依托某电站地下洞室群工程，通过萤火虫智能分析方法结合数值模拟模型进行岩体参数反分析，并利用反演参数对典型地下洞室围岩的稳定性进行分析与探讨.

1 萤火虫算法-最小二乘支持向量机反分析方法研究

1.1 萤火虫算法

萤火虫算法（Firefly Algorithm，FA）［13-15］是一种新兴的群智能优化技术，它通过模拟自然界中萤火虫个体根据发光亮度寻找其他个体的特性，以发光亮度为依据，模拟萤火虫算法在模拟萤火虫个体之间相互吸引、相互靠近的过程，以达到寻优的目的. 基于萤火虫算法的最小搜索模型有三个假设：

①萤火虫个体之间没有任何区别，只考虑发光的亮度，亮度决定萤火虫个体的吸引力，亮度低的被亮度高的萤火虫吸引；②假如一定区域内的某个萤火虫个体为局部最优，则进行随机移动；③萤火虫个体的适应值函数决定其亮度.

上面三条件决定的萤火虫算法，萤火虫个体根据其他萤火虫亮度而发生追逐的行为即为最优解的求解过程. 在FA算法中，萤火虫局部最优搜索的关键是萤火虫的发光亮度和相互吸引度［16-17］. 假设所有个体的初始值相同，个体的亮度决定该个体对其他个体的吸引度，个体的目标函数值决定个体的发光亮度，在固定的种群范围内，当前最优解即为种群中亮度最高的个体，为防止陷入局部最优，处于最优解的个体再在一定范围内随机移动，其他个体会根据个体间的吸引度向吸引度大的个体靠近，从而找出全局最优值.

1.2 支持向量机算法原理

支持向量机算法［18］以统计学习理论为基础，根据有限学习样本构建复杂模型，从而获得更好的推广能力. 最小二乘支持向量机（Least squares support vector machine，简称LSSVM）是标准支持向量机的一种扩展［19］，它将二次规划问题转化为线性方程组求解，降低了计算的复杂性，提高了求解速度. LSSVM在训练的过程中，采取等式约束、经验风险偏差的二次方和线性方程组求解的方式进行，同时也需要遵循结构风险最小化原则. 最小二乘支持向量机的标准输出量可定义为：

式中：K(x,xi)为满足Mercer条件的核函数；αi为回归系数；b为回归偏差；x为样本变量；xi为第i个预测样本值；l为样本数目.

对于实际工程的设计和施工来说，通过待反演参数获得的计算值和实测值之间的差距越小，则说明位移反分析获得最有效果. 因此，目标函数通常取以下形式：

其中：x={X1,X2,…,Xk,…,Xm}，为m个参数构成的一个样本，Xk为待反演的岩体参数，如弹性模量E、内聚力c、内摩擦角φ等；m为待反演的岩体参数个数；fi(X)为第i个位移计算值，ui为位移实测值；n为测点数.

与传统SVM 相比，LSSVM 计算复杂性更低，求解完成的解答中包含多数甚至全部的训练样本，更加便于在线计算，能够通过扩展其他形式来处理动态问题，在求解过程中结合支持向量机、前向神经网络等进行讨论，与机器学习的本质更加贴合.

1.3 参数反分析模型

根据以上原理，建立萤火虫-支持向量机参数反分析计算步骤如图1所示，利用该流程可将数值计算结果作为样本进行学习，从而得到与监测位移相互吻合的岩体力学参数.

图1 参数反分析模型Fig.1 Parameter back analysis model

2 工程概况及数值模型建立

江苏溧阳抽水蓄能电站属破碎岩层大型地下厂房洞室群工程，该工程地质条件十分复杂，总体上岩体以完整性较差为主，岩体结构主要为镶嵌碎裂结构，部分属于层状结构，围岩稳定问题突出，施工安全压力大. 围岩类别：Ⅲ类占45%，Ⅳ类占52%，Ⅴ类约占3%.

根据溧阳抽水蓄能电站地下洞室群工程地质条件，建立地下洞室群块体离散元数值模型［18］，如图2 所示，地下洞室群工程数值计算模型共计约3 万个块体，采取对边界进行法向位移约束的形式，其中主厂房区域整体计算模型如图3所示.

图2 溧阳抽水蓄能电站地下洞室群布置图Fig.2 Underground cave group layout of Liyang pumped storage power station

图3 主厂房整体计算模型三维示意图Fig.3 Three-dimensional diagram of the overall calculation model of the main workshop

3 主厂房围岩力学参数反分析研究

地下洞室工程中，洞室的变形往往与侧向约束条件有很大关系，因此必须考虑初始应力场. 在考虑主厂房区域地应力影响时需要考虑两个水平向侧压力系数，并转化到模型坐标系下，因此，在反演研究区域初始地应力时，实际反演两个侧压力系数. 围岩变形的主要影响因素为岩体弹性模量E和泊松比μ. 地下洞室变形侧向变形主要取决于应力场的变化，而泊松比影响较小，但是弹性模量对于围岩变形有着较大影响. 因此本文在反演参数的选取中选择弹性模量E，而采用确定的泊松比μ. 岩体力学建议取值如表1、表2所示.

表2 结构面物理力学参数Tab.2 Physical and mechanical parameters of structural surfaces

根据正交设计原理设计试验方案，利用块体离散元数值模型计算得到各个方案下的计算位移值，可得LSSVM算法的训练样本如表3所示. 采用萤火虫算法优化模型参数，得到反映岩体参数和围岩位移之间映射关系的最优LSSVM模型.

表3 正交试验获得的正演训练样本Tab.3 Forward training samples obtained by orthogonal experiment

将方案1～16作为训练样本，基于萤火虫算法的最小二乘支持向量机算法，以实际监测位移作为样本输入到已经训练完成的LSSVM模型中，利用图1流程进行学习预测. 得出LSSVM最佳参数组合如表4所示.

表4 围岩参数反演结果Tab.4 Inversion results of surrounding rock parameters

将表4 中反演得到的参数代入数值计算模型，进行正分析，选取典型（C1）断面的6 个监测点进行现场实测变形和数值模拟监测点变形进行对比分析，如图4、图5 所示. 可见现场实测变形和数值模拟监测点变形具有较好的对应，各点变形与现场监测数据基本一致. 开挖过程中围岩变形最大位于拱顶，现场实测值为15.1 mm，数值模拟值为12.2 mm；围岩变形最小位于上游边墙，现场实测值为-0.5 mm，数值模拟值为2.2 mm；由于受开挖过程中母线洞开挖的影响，明显存在下游围岩变形大于上游的现象，这些现象都与工程实际相符合，因此该模型获得的反演参数较为准确，可为相关研究提供参考（图4 和图5 中MC1-1，MC1-2，MC1-3，MC1-4，MC1-5，MC1-6都为监测点编号）.

图4 现场实测变形曲线Fig.4 Deformation curves measured on site

图5 监测点变形曲线的数值模拟Fig.5 Numerical simulations of deformation curves of monitoring points

4 多洞联合作用下的变形协同分析

溧阳抽水蓄能电站地下洞室群，共有3大洞室，分别是主厂房、主变洞和尾闸室，其中，主厂房共分为7层开挖，主变洞共分为3层开挖，整个地下洞室群共分8步进行开挖，其中，第4层开挖为母线洞，第8层开挖为交通洞.

4.1 分层开挖下的围岩位移场分布特征

分层开挖步骤下的围岩位移场分布如图6所示. 由图6（a）可知，第一层开挖洞室顶拱部分后，洞室围岩主要产生竖直方向的变形，其中，拱顶以下沉变形为主，拱底位置处由于产生向上的回弹变形而出现“底臌”现象. 随着第二层进行开挖，洞室高边墙逐渐形成，此时洞室围岩变形逐渐转变为以水平变形为主. 第二层开挖后，开挖面逐渐远离顶拱部分，顶拱围岩变形受开挖扰动影响越来越小，围岩变形速度减慢. 此时，洞室上下游边墙的变形基本对称分布. 随着开挖的进行，主厂房下游侧围岩变形区和主变洞上游侧变形区发生交接，两洞室的围岩变形开始变为受2次开挖扰动影响，在双重开挖扰动作用下，主厂房边墙开始发生非对称变形，下游侧壁围岩变形明显增大. 洞室开挖完成后，从图6（h）中可以看出，主厂房变形最为严重，最大变形点位于下游侧岩锚梁处，主厂房顶拱围岩最大变形发生在上游侧拱顶处，边墙最大变形发生在下游侧边墙处；主变洞顶拱围岩最大变形同样发生在上游侧拱顶位置处，边墙最大变形发生在上游侧边墙.

图6 不同开挖步骤下的位移图Fig.6 Displacement diagrams under different excavation steps

总体上，洞室顶拱和底板以竖向变形为主，高边墙以水平变形为主. 顶拱变形受开挖面扰动影响越来越弱，围岩变形速度越来越慢，边墙围岩水平变形受开挖扰动影响较大，尤其是受母线洞的开挖影响，开挖母线洞时，主厂房、主变洞边墙变形发生明显的增大. 在实际开挖过程中，应加强支护工作，尤其是主厂房下游侧岩锚梁处的支护.

4.2 分层开挖下的围岩应力场分布特征

图7为地下洞室群在不同开挖步骤下的最大主应力云图. 由图可以明显地看出，由于地下厂房的开挖扰动作用，在洞室开挖面附近产生了明显的应力重分布现象，形成一定范围的拉应力区，在远离洞室开挖面一定距离，围岩处于原岩应力状态，处于压应力区.

由图7（a）可知，地下洞室群第一层顶拱部分开挖后，在拱顶、拱底产生较大的拉应力区，引起厂房顶拱发生拱顶下沉、拱底臌起现象，通过4.1节的洞室开挖位移场可以得到印证. 同时，在顶拱与边墙交接处形成了明显的应力集中现场. 由图7（b）可知，从第二层开挖开始，地下厂房的高边墙逐渐形成，洞室边墙临空部分越来越大，在应力场作用下，边墙产生以水平方向为主的变形，从图中可以看出在洞室边墙位置形成较大范围的拉应力区，在拉应力的作用下洞室围压产生向洞室方向的变形，产生“片帮”现场，其中，由图7（d）可知，在洞室第4步开挖母线洞后，在主厂房、主变洞与母线洞交接部位形成明显的应力集中现象. 随着开挖的进行，厂房顶拱部分距离开挖面越来越远，开挖扰动作用对顶拱的影响越来越弱，集中的应力逐渐得到释放，应力水平逐渐降低.

由图7（h）可知，整个地下洞室群开挖完成后，在3大洞室的边墙部位形成了较大范围的拉应力区域，应力集中主要发生在顶拱与边墙、边墙与母线洞交界处.

4.3 分层开挖下的围岩塑性区分布特征

图8为地下洞室群在不同开挖步骤下的各剖面塑性区分布图. 可以看出，地下厂房开挖后，在洞室表层围岩处产生一定范围的剪拉应力带.

由图8（a）可知，3大洞室顶拱部分开挖后，在拱顶部位形成了剪切应力带，第2步开挖开始后，厂房顶拱的应力集中现象更加明显，拱顶产生少量剪切破坏，随着开挖面的不断深入，顶拱部位围压逐渐趋于稳定，后续洞室开挖对顶拱塑性区分布影响不明显；洞室高边墙逐渐开挖形成后，在洞室围岩尤其是临空面附近围岩形成了较大区域的剪切和张拉塑性破坏区，其中，母线洞岩体开挖后，在母线洞与主厂房、母线洞与主变洞交接部位形成了明显的应力集中现象.

图8 不同开挖步骤下的典型剖面塑性区图Fig.8 Plastic zone diagrams of typical sections under different excavation steps

随着厂房分层开挖的不断进行，洞室围岩中产生的塑性区范围也在不断增大，所产生的破坏形式主要为剪切破坏和剪切张拉破坏，引起的破坏区主要集中在洞室上下游边墙部位、边墙与母线洞交接部位. 对比分析主厂房上下游侧的围岩塑性区范围可知，主厂房下游侧产生的塑性区范围要比上游侧塑性区大，主变洞上游侧产生的应力集中要比下游侧应力集中明显，主要原因是母线洞的开挖扰动对主厂房的下游侧、主变洞的上游侧影响比较大.