基于残差平方和取舍常数项的地应力场反演分析

蒙 伟，何 川，严 健，吴枋胤，周子寒，寇 昊

(西南交通大学 交通隧道工程教育部重点实验室，四川 成都 610031)

近年来，随着复杂艰险山区高速公路、铁路等重大工程和西部大开发相继实施，在这些区域逐渐涌现出大量长大、深埋隧道等地下工程[1]。该区域大量地应力资料表明，竖直地应力小于水平地应力的现象非常普遍，即中国西部地区普遍存在挤压构造运动[2-3]。故对此区域隧道等地下工程岩体的初始地应力场进行反演分析具有重要的实际意义[4]。

在反演过程中，郭怀志等[5]给出了能够获得唯一解的反演分析方法，故此方法被广泛地沿用至今。文献[5-9]、[10-14]对岩体的初始地应力场进行反演分析采用的回归模型分别为

(1)

(2)

由式(1)、式(2)可知，依据初始地应力场回归模型是否含有常数项，回归模型可分为两类：①不含常数项的回归模型;②含常数项的回归模型。然而，截至目前，在岩体初始地应力场回归反演过程中，含与不含常数项回归模型被单独地用于反演岩体的初始地应力场。由于这两类回归模型的反演结果不相同，因此，仅采用其中一类回归模型进行反演会存在一个问题，即在这两类回归模型中，反演结果较差的回归模型很可能被用于反演。为解决此问题，提高岩体初始地应力场的反演结果，在反演过程中，有必要对如何选择最佳反演结果的回归模型进行分析研究。

本文为客观地选择最佳反演结果的回归模型，结合概率统计理论提出了衡量岩体初始地应力场反演结果的量化指标。然后在郭怀志教授提出的反演分析方法的基础上，在含与不含常数项回归模型中，基于量化指标提出了客观获得最佳反演结果的反演流程。最后，以某复杂艰险山区铁路隧道为依托对反演流程进行工程应用分析。

1 取舍回归模型常数项的反演机理

1.1 衡量反演结果的量化指标

由概率统计相关理论可知，残差平方和可描述随机因素影响引起的偏差，即残差平方和可反映实测值与计算值之间的偏差程度。

(3)

假设y由一个自变量组成，样本容量为5，展示含与不含常数项回归模型实测值与计算值的偏差程度见图1。

图1 含与不含常数项回归模型的偏差程度

由图1可知，由于实测值是固定不变的，当用含与不含常数项回归模型去匹配实测值时，实测值与计算值之间的偏差程度是不相同的。同理，在岩体初始地应力场回归反演过程中，采用含与不含常数项回归模型的反演结果也不相同，故有必要对如何选择最佳反演结果的回归模型进行分析研究。

通过对比式(3)与式(1)、式(2)可知，不含与含常数项回归模型的应力实测值与计算值的偏差程度分别为

(4)

(5)

应力是矢量，除了大小还具有方向，即合理的岩体初始地应力场反演结果应满足应力实测值与计算值在大小和方向上基本一致[9]。因此，需采用应力大小残差平方和与应力方向残差平方和两个指标表示应力实测值与计算值之间的偏差程度。应力方向残差平方和为

(6)

为综合评价应力实测值与计算值的大小和方向偏差程度，需把应力大小残差平方和与应力方向残差平方和两个指标放到一个计算式中，而应力大小残差平方和的单位(Pa2或MPa2)与应力方向残差平方和的单位[(°)2]不相同，如果采用应力大小残差平方和加应力方向残差平方和的形式则还需对其进行标准化处理。为方便工程人员以后应用此方法，故采取应力大小残差平方和乘以应力方向残差平方和的形式，计算式为

(7)

1.2 基于残差平方和取舍常数项的反演流程

在反演过程中，首先同时采用含与不含常数项回归模型进行反演，然后通过式(7)客观地选择其中最佳反演结果的回归模型(较目前仅采用其中一类回归模型的反演流程，可避免采用反演结果较差的回归模型进行反演)，最后采用该回归模型即可得到与实测原位地应力较为吻合的岩体初始地应力场。如果采用最小二乘法求解回归系数，则具体的实施方案见流程图2。

图2 基于残差平方和取舍常数项的地应力场反演流程

2 工程概况

某隧道全长约16 km，最大埋深1 350 m，为某新建复杂艰险山区铁路隧道的关键控制性工程，属于单洞单线隧道，处于中国地壳板块构造运动最强烈的区域之一(青藏高原)。隧址区主要分布的是闪长岩等硬岩，弹性模量为40 GPa，隧址区属构造剥蚀区，山势陡峻，沟谷深切，最大切割深度为2 100 m，典型高山峡谷地貌，加上地壳板块构造运动强烈，在隧道开挖过程中频繁发生中等及强烈岩爆。

在隧道DK186+183—192附近伴有强烈岩爆发生，为研究该区域附近复杂的初始地应力场，本文以DK183+000—188+000深埋峡谷段为例进行反演分析，隧道穿越的岩体物理参数见表1。该区域通过水压致裂法测得的原位地应力见表2[17]。此外，隧道还存在高地温，记录的最高地温为89.5 ℃，且地温高居我国之首[15-16]。

表1 岩体物理参数

表2 实测原位地应力

由表2可知，最大水平主应力大于垂直应力，故在该水压致裂钻孔附近存在强烈的水平构造运动，与隧道在拱顶附近发生岩爆现象基本一致，表明由于板块构造运动，该钻孔附近区域岩体的初始地应力场较自重应力场更加复杂，且DK-SZLSD-2钻孔5个测点的原位地应力较难反映该局部区域岩体初始地应力场的宏观分布规律，故有必要对该局部区域岩体的初始地应力场进行反演分析。

3 工程应用

3.1 建立三维数值模型

数值模型的网格及地形见图3，共计2 853 904个单元，3 863 303个节点，模型底面为隧道轴线标高以下1 000 m，顶面为岩体实际地形。为方便获取垂直于隧道轴线方向的应力，数值模型以隧道轴线方向为x轴，垂直隧道轴线方向为y轴，沿隧道埋深铅直方向为z轴。

图3 划分网格后的数值模型(单位：m)

3.2 定义边界条件

因原位地应力采用水压致裂法进行测量，见表2。受其原理限制，所测得的原位地应力无铅垂面内的剪切应力，故剪切构造荷载仅考虑水平xoy平面内的构造荷载[5, 11, 14]，施加的自重及构造荷载大小见图4。施加方法：除数值模型的顶部边界，其他非加载的边界均施加法向位移零的约束；对图4(a)数值模型施加竖向9.8 m/s2重力加速度；对图4(b) 数值模型施加沿x轴的1 MPa水平均布荷载；对图4(c) 数值模型施加沿y轴的1 MPa水平均布荷载；对图4(d) 数值模型施加x方向的1 kPa/m水平三角形荷载；对图4(e)数值模型施加y方向的1 kPa/m水平三角形荷载；对图4(f)施加水平剪切荷载，依据文献[18]提出的公式，可对该数值模型垂直于x轴的边界施加5 000 m/100 000=0.05 m的剪切位移，对该数值模型垂直于y轴的边界施加3 000 m/100 000=0.03 m的剪切位移，即可模拟得到该模型的水平剪切构造应力场。

图4 荷载分布

3.3 转换原位地应力坐标系

如图3所示，为方便获取垂直于隧道轴线方向的应力，数值模型以隧道轴线方向为x轴，因此数值模型采用的xoy坐标系与实测原位地应力采用的NoW坐标系不一致，故需要对实测原位地应力进行坐标变换。应力由旧坐标系转换到新坐标系时，弹性力学应力转轴式为

σu′v′=αu′uαv′vσuv

(8)

式中：σuv、σu′v′分别为转换前后的应力张量；αu′u、αv′v为转换系数。

含与不含常数项回归模型的反演结果见表3。

表3 反演结果

3.4 回归反演岩体初始地应力场

依据式(1)和式(2)的回归模型，结合图4的自重和构造荷载，以及温度应力[16]，隧道采用的不含与含常数项回归模型分别为

(9)

(10)

使应力实测值与计算值的残差平方和最小即可求解得到回归系数。在式(2)的含常数项回归模型中，即需使式(5)最小即可，为使残差平方和最小，则需求残差平方和分别关于回归系数的偏导数，并使偏导数为零，即有

s=1,2,…,n

(11)

化简式(11)有

(12)

式(12)为通过最小二乘法求解回归系数的正规方程组，为方便求解回归模型的回归系数，可将式(12)改写成直观的矩阵形式为

(13)

同理，对于式(1)的不含常数项回归模型，通过最小二乘法求解回归系数的矩阵形式正规方程组为

(14)

由式(13)、式(14)可知，在通过最小二乘法求解回归系数时，所采用的正规方程组是不相同的。因此，在求解回归系数的过程中需要注意，若采用含常数项回归模型，则应采用式(13)的正规方程组求解回归系数；若采用不含常数项回归模型，则应采用式(14)的正规方程组求解回归系数。而在目前的文献中，大都选择式(14)求解回归系数，却得到了含常数项的回归模型，这在今后岩体初始地应力场的反演中应避免。

通过回归模型及回归系数的显著性检验后得到的结果分别为

(15)

(16)

式中：σT≈0.016ZMPa。其中，Z=H-h，H为埋深，h为岩体恒温带离陆地表面的距离[16]。

回归得到式(15)的相关系数为0.997，式(16)的相关系数为0.992，由此可见，两类回归模型子构造应力场与岩体初始地应力场之间线性相关的密切程度都较高。

3.5 基于残差平方和取舍回归模型常数项

依据式(4)、式(5)可得不含、含常数项的回归模型应力实测值与计算值的大小残差平方和分别为

(17)

(18)

(19)

tanα=(σ1-σx)/τxy

(20)

(21)

由式(7)可得含与不含常数项回归模型应力的大小和方向偏差程度，分别为1 185.75、2 442.96。

由此可见，隧道采用含常数项回归模型的应力大小和方向偏差程度较不含常数项回归模型更小，故采用含常数项回归模型的反演结果更好。因此，从反演结果出发，本文建议采用含常数项回归模型对此局部区域的岩体初始地应力场进行反演分析。如果不采用图2的反演流程，则有可能采用不含常数项回归模型进行反演，而采用此类模型反演得到的反演结果较含常数项回归模型更差。含与不含常数项回归模型的反演结果对比见图5～图10。图中，εi=实测值-计算值。

图5 σx残差平方和

图6 σy残差平方和

图7 σz残差平方和

图8 τxy残差平方和

图9 含常数项回归模型方位角的残差平方和

图10 不含常数项回归模型方位角残差平方和

由图5～图10可知，两类回归模型应力实测值与计算值在量值上存在一定的差异。含常数项回归模型应力大小最大残差为1.5 MPa，应力方向最大残差为-7.8°；不含常数项回归模型应力大小最大残差为1.7 MPa，应力方向最大残差为10.3°。即采用含常数项回归模型反演得到应力实测值与计算值的大小和方向偏差程度更小，反演结果较不含常数项回归模型更好，与式(22)、式(23)的结果一致。

综上所述，在岩体初始地应力场反演过程中，采用含与不含常数项回归模型得到的反演结果是不相同的。因此，仅采用其中一类回归模型进行反演，很可能会导致反演结果较差的回归模型被用于反演，从而导致反演得到的岩体初始地应力场不是最优的。故为提高岩体初始地应力场的反演结果，在反演过程中，以应力的大小和方向偏差程度为衡量标准，采用图2的反演流程可客观地获得与实测原位地应力较为吻合的岩体初始地应力场。

4 结论

(3)在岩体初始地应力场反演中，若采用含常数项回归模型，则应采用式(13)的正规方程组求解回归系数；若采用不含常数项回归模型，则应采用式(14)的正规方程组求解回归系数。