基于变形监测的软岩隧道有限元反演分析研究

周晓靖，罗红星，钟明文，秦雨樵

(1.云南楚大高速公路扩建工程有限公司， 云南 大理 671000； 2.云南大永高速公路建设指挥部， 云南 大理 671000；3.中国科学院 武汉岩土力学研究所， 岩土力学与工程国家重点实验室， 湖北 武汉 430071； 4.中国科学院大学， 北京 100000)

在软岩地层开挖隧道的过程中由于软岩自身的流变性质，围岩往往会出现较大的变形，可能会导致严重的工程事故。所以很有必要对开挖后的围岩进行及时的变形位移监测，并通过监测数据反演出岩体蠕变的各项参数，从而更好地服务于支护设计以及选择合适的加固时机。在这些方面，国内外学者做了大量的工作。Sakurai S等[1]提出优化方法的思想，采用逆解法反算隧道围岩的弹性模量以及周围的地应力。冯紫良等[2]结合位移数据反演得到岩体的初始应力场。王芝银等[3]使用有限元以及边界元方法分析了常见流变模型，并提出逆解优化以及回归的方法。叶飞等[4]分析了不同位移量测项目得到的时间位移曲线，并给出数据处理方法。朱合华等[5]采用参数反分析以及位移计算反复迭代的方法较为精确地预测了后续施工的位移状态。李德宏[6]则结合参数反演以及施工检测方法对结构设计进行指导。宋桂锋等[7]则利用监测数据分析围岩以及支护的变形受力特征。近些年来，信息理论、系统理论等多种新技术新思路也运用到实际工程中来。比如刘维宁等[8]在概率论的基础上研究了反分析结果的稳定性问题。袁勇等[9]则进一步完善了概率反分析方法。同时大量基于人工智能或者BP神经网络的研究也引入到参数反分析当中[10-14]。总而言之，基于位移监测的岩土参数反分析已经普遍运用到工程实际当中，具有良好的效果。

本文以云南大丽高速公路天井山隧道工程为依托，首先采用非接触测量对软岩隧道围岩变形进行监测，初步判断稳定性情况。然后结合有限元参数反演方法以及监测位移数据，对软岩蠕变参数进行反分析计算，并将得到的数值计算结果与实测数据进行比较。最后用反演得到的参数分析评价隧道围岩的稳定性，为工程施工提供建设性意见。

1 工程概况

天井山隧道是云南大丽高速公路毕节至威宁高速公路的一段，为一座分离式隧道。隧道左右线总长约为2 000 m，两侧中线距离平均为25 m。该隧道穿越洱海断陷盆地周边的中山区，该区地形切割较深，地势较陡峻，同时地形复杂，沟谷纵横。根据地质调查和钻孔资料，隧道左、右幅处于相同地貌单元。该区域需要关注的地质层位是砂岩、泥岩厚互层。这些层位的岩体破碎，岩质较软，特别是开挖之后，风化崩解，遇水软化。而强风化岩体的节理裂隙也十分发育。总体而言，工程区岩体较软，表现出较强的软化以及蠕变性质，对于隧洞的开挖以及长期稳定性均不利，需要进一步通过现场监测以及数值计算进行分析。

2 隧道围岩变形监测

2.1 监测手段

围岩的变形测量可以采用直接与围岩接触的接触变形测量方法，也可以采用全站仪设站的非接触变形测量方法。接触测量方法由于受到施工条件的限制以及人为测量因素的干扰，不能很好地满足长距离大断面隧道施工的变形监测要求[15]。而非接触变形测量方法是一种在不接触测量点的前提下测量该点位移的方法。这种方法(特别是运用全站仪的监测预报体系)能高效且精确地完成测量任务，并受人为因素影响较小。本文即采用非接触变形量测方法对软岩隧道围岩进行监测，所用观测仪器为瑞士徕卡TS30。

在非接触变形测量方法中，除了需要满足工程施工的要求以及关键点位以外，还需要尽可能地满足测网图形条件以及要求2个以上的能够由前视点转化为后视点的测点。所以这些测点在前后2个测点都能通视，并能被同一点位找准。对于莱卡全站仪而言，其测距范围在25 m～80 m之间。但是由于隧道施工中会产生较大的尘土，导致能见度很低，不适用于太长的测量范围。所以在实际测量中，必须要准备强光射灯用以照明[16]。

根据上述的各项要求，在天井山隧道左洞布置3个监测断面，分别为：ZK0+655，ZK0+665和ZK0+675；在天井山隧道右洞布置5个监测断面，分别为：YK0+670、YK0+660、YK0+650、YK0+640和YK0+630。同时在每个断面处布置3个沉降观测点(A1、A2、A3)，3条收敛观测线(AB、AC、BC)，如图1、图2所示。由于左右两洞基本对称，故在本文分析中仅选用左洞的三个监测断面进行分析。

图1 天井山隧道沉降观测点布置图

图2天井山隧道收敛观测点布置图

2.2 监测结果以及分析

隧道交叉段左洞的监测累计值见表1，同时在图3中画出了ZK0+675断面不同监测点沉降累计值以及监测线收敛累计值随时间变化的曲线。ZK0+665断面和ZK0+655断面监测点沉降累计值与监测线收敛累计值随时间变化曲线与ZK0+675断面相似，不再画图。从图表中可以看出，对于断面ZK0+675,监测沉降点A1、A2、A3点均在第15 d沉降稳定，与该断面的监测线收敛积累量稳定的时间一致，但是对于沉降量来说，左侧累计沉降量明显偏大，在侧线收敛累计值中也出现了同样的规律。对于断面ZK0+665，A1点在第11 d已经稳定，而A3点则在第16 d才稳定。同时与该侧点相关的两条监测线收敛累计量也晚于另一条监测线稳定，累计值也同样明显较大。对于断面ZK0+655，AB、BC、AC三条监测线的收敛位移值完全相同，但左侧A2点的沉降值明显小于其他两点。三个断面的围岩变形速率随时间的增加逐渐减少，最后达到收敛，并没有发生大的位移突变，说明此时围岩在开挖后均达到稳定状态。但是在这三个断面的围岩变形并不是对称的，而是或多或少地向某一方向侧集中。这种现象可能会导致隧洞另一侧出现拉应力，从而威胁到围岩的稳定性，需要特别注意。

表1 隧道交叉段左洞观测累计值(数据来自施工监测)

图3各断面沉降累计值以及监测线收敛累计值随时间变化曲线

3 基于变形监测的有限元参数反演方法

本文所关注的隧道岩体因为本身强度较低，表现出较强的非线性蠕变特性。这意味着需要用较为复杂地蠕变本构方程来进行描述。为了与实际工程情况相对应，有必要运用有限元参数反演的方法来标定模型的各种参数。而非接触测量除了可以监测围岩的变形收敛情况之外，也可以为参数反演过程提供重要的目标依据。由于监测了三个断面的变形随时间的变化情况，可以考虑通过其中一个断面的数据进行参数反演，并用实测结果进行验证。这样就能科学客观地得到围岩岩体的参数，方便进一步地分析计算。

3.1 有限元优化反演法的基本理论

有限元参数反演的目的是寻找一组参数X=(x1，x2，…，xm)T使目标函数趋近无穷小值。这个问题由于其复杂性，通常需要增加一些约束条件将其转化为约束优化问题来处理，才能保证得到稳定且唯一的解。

一般而言，非线性优化算法需要求得目标函数的导数。但是反演模型的目标函数很复杂，很难用解析表达式直接表达。这时，要直接求出目标函数的导数就有些不切实际。为了克服该问题，可以采用Nelder-Mead法。这种方法不必计算函数导数，算法比较直接而简单，能很好解决变量不是很多的方程极值问题[16]。其基本思想是：先给出几个待反演参数的初始值，并带入方程求解出目标函数值。然后根据这些函数值的大小关系找出目标函数下降的方向，在下降的方向上可以再找出一组新的反演参数值，再带入计算求出目标函数值。接着比较新值和上一次计算出的函数值，找出新的函数下降方向。如此不断搜索，直到最终满足收敛要求为止。虽然该方法只适于求解无约束的最优化问题，但可通过罚函数法将非线性约束最优化问题转化为无约束优化问题，从而满足了其要求。

3.2 蠕变损伤模型的数值实现

ABAQUS提供了基于扩展Drucker-Prager屈服准则的线性蠕变模型，其蠕变势函数采用与屈服势函数相同的形式，在p-q面上为双曲线函数：

(1)

岩石的蠕变过程总是伴随着其内部的损伤和弱化，因此，整个蠕变过程是非线性的。为精确模拟隧道围岩的蠕变力学行为，需考虑蠕变损伤的影响。

非线性蠕变损伤模型能较好地反映岩石蠕变的第一阶段和第二阶段，即瞬态蠕变和稳态蠕变，但不能反映加速蠕变阶段，而加速蠕变对工程来说是更有意义的。因此，基于验模型式，提出隧道交叉段软岩的蠕变损伤力学模型。

(2)

当考虑到蠕变损伤时，采用下面的损伤方程式进行描述：

(3)

式中：εc为蠕变应变；Dc为蠕变损伤；q为等效应力；t为时间；A，m，n，α为材料常数，是本次有限元反演优化计算的目标值。

3.3 有限元优化反演法的实现

基于上述提出的理论与方法，在大型有限元程序ABAQUS的基础上，以MATLAB语言为平台，结合精确罚函数法以及Nelder-Mead算法，实现了基于有限元的反演算法，其具体流程见图4：

图4岩土体参数反演程序框图

(1) 首先输入目标实测值、待反演参数初始值以及优化算法的控制精度。

(2) 然后读取上述提到的有限元程序命令流文件，调用ABAQUS软件进行有限元计算。

(3) 通过读取计算结果文件，获取测点的有限元计算值，计算目标函数F(X,μ)的值，并进行优化处理。

(4) 如果得到的函数值F小于收敛精度则求出了反演参数X，输出结果；若大于收敛精度则修正了反演参数X，并重新计算(2)-(3)步直至满足收敛精度。

4 软岩隧道稳定性分析

为了分析该软岩隧道的稳定性，首先按照本文提出的反演方法和思路，对天井山隧道开展平面有限元反演计算，得到其岩体蠕变本构模型的参数。然后分析其塑性区以及位移的分布情况，以此来分析该隧道的稳定性情况。

4.1 天井山隧道有限元反演结果

由于天井山隧道ZK0+655断面沉降观测数据以及BC线收敛观测数据具有较好的稳定性，选取这几个数据随时间的变化规律作为反演目标。由式(3)可知，需要反演参数为A、m、n、α共计4个蠕变参数。

根据隧道交叉段工程地质特征，建立平面有限元计算模型,尺寸为120 m×120 m，其中Y轴为重力方向。有限元计算模型见图5。在本模型中不仅模拟了隧道围岩，同时也模拟了隧道初衬。有限元计算模型加载模式为底面和侧边采用法向约束，顶部按围岩自重施工荷载；侧面按侧压力系数施加侧向应力，根据地实力地质构造特征，隧道交叉段侧压力系数取1.5。有限元加载模式见图6。根据工程地质条件和岩石力学试验结果，隧道交叉段围岩有限元计算力学参数如表2所示。

表2 地下洞室区岩体参数取值

遵循由3.3中提出的反演参数过程，采用ABAQUS以及MATLAB程序进行反演分析，在反演未知蠕变参数的同时，也对已知的弹性模量、泊松比等已知参量进行调整，以便得到与实际变形情况更为符合的结果。最后反演得到岩体参数结果见表3。

图5 平面有限元计算模型图

图6 有限元模型加载模式

为了验证反演参数的合理性，将反演结果与实测结果进行对比，见图7。从图7中可以看出，四个不同位置计算值的变化趋势与实测数据吻合较好，其变形速率均随时间的增长而降低，并最终趋于稳定。特别是在监测初期，两者得到的数据基本一致。从这些证据中可以看出，通过本文提出的参数反分析计算求得的参数能很好地反映隧道开挖后变形过程中的状态。

4.2 隧道围岩稳定性分析

通过这些参数以及反演计算的结果，可以进一步分析隧道围岩的稳定性情况。观测结束时围岩塑性区分布如图8所示。围岩塑性区主要分布在隧道的顶板以及底板附近，其中顶部塑性区范围达7.15 m，底部围岩塑性区范围达7.71 m。特别是在底部围岩处，塑性区域非常集中。而在隧道的两侧，塑性区区域相较而言并不太发育。同时从图9中可以看出，当监测结束，初次衬砌已全部进入塑性，主要集中区域也发生在底部，顶部的塑性也较大。而隧道位移分布图10也揭示出隧道开挖后即使加上初衬，底部隆起的位移也很大，同时在两侧，位移影响区域非常大，基本快达到模型边界。其中，隧洞顶部沉降值为4.4 cm，最大水平收敛值为4.5 cm。说明隧道上部也存在向内滑塌的可能。总体而言，隧道围岩在没有恰当支护措施的情况下，很难保持稳定，应提出合理支护调整方案和二衬支护时机。

5 结 论

本文介绍了利用全站仪设备对软岩隧道进行变形监测的方法，采用有限元计算分析方法对软岩蠕变参数进行了反演计算，并以此为依据分析隧道围岩的稳定性，得到以下结论：

(1) 根据现场实际监测数据,隧道围岩的变形趋势总体上是随时间变化逐渐收敛的，但变形呈现出较强的不对称性，需要着重注意。

图7 观测值与实测值对比图

图8 围岩塑性区分布云图

图9 初衬塑性区分布云图

图10隧洞位移分布云图

(2) 根据现场围岩沉降与收敛观测数据，采用有限元反演分析方法，得到隧道交叉段围岩相关流变参数，并以实测参数变化规律作为对比，验证了参数的有效性。

(3) 隧洞变形有限元反演结果表明，该洞段围岩软弱、流变特征显著，现有支护措施下，由于围岩流变塑性区范围将持续增长，最终会导致初砌破坏达不到支撑效果，因此，应提出合理支护调整方案和二衬支护时机。