穿越多断层公路隧道初始地应力场反演分析

Inverse Analysis of the Initial Ground Stress Field of a Highway Tunnel through Multiple Faults

摘要: 本文以大丰华高速公路鸿图嶂隧道为依托，建立ANSYS三维数值计算模型，通过现场水压致裂法测得的数据，进行多断层下隧道周边区域初始地应力场反演计算分析；研究结果表明：隧道工程区地应力场的分布基本以构造应力为主，最大水平主应力和垂直应力有随隧道埋深的增大而增大的趋势，且在断层处出现应力释放现象。研究结果为后续设计和施工阶段提供基础数据。

Abstract: Based on HongTuZhang Tunnel of DaFengHua Expressway, this paper establishes the ANSYS 3D numerical calculation model and performs the inversion analysis of the initial ground stress field in the tunnel surrounding area under multiple faults through the data measured by hydro-fracturing ground stress measurement; the research results show that: the distribution of the ground stress field in the tunnel engineering area is basically dominated by tectonic stress, and the maximum horizontal principal stress and vertical stress have a tendency to increase with the increase of tunnel burial depth. The maximum horizontal principal stress and vertical stress have a tendency to increase with the increase of tunnel depth, and the stress release phenomenon appears at the fault. The results of the study provide basic data for the subsequent design and construction stages.

1. 引言

近年来，随着公路交通基础设施的不断发展，“多、长、大、深”的公路隧道数量逐年增加，随着越岭深埋特长公路隧道的不断出现，深部工程问题变得越发重要。岩体的初始应力场大小与隧道开挖过程中的围岩稳定性、地质灾害等的发生密切相关，因此加强隧道工程区的初始应力场的基础性研究，不仅对设计阶段起指导性意义，同时可为后续的施工期间工程的安全性研究提供基础性资料。目前国内外对于岩体初始地应力场分析已有大量研究 [1] [2] [3] [4]，但其内容主要针对于隧道轴线应力本身或高地温因素的影响，对隧道穿越多断层的模拟及对断层附近地应力变化规律的研究并无深入讨论，本文以广东省大丰华高速公路鸿图嶂隧道为工程依托，针对断层对隧道地应力分布的影响进行探讨分析。

2. 工程概况

大(埔)丰(顺)(五)华高速公路(简称“大丰华高速”)是广东省拟规划高速公路网中的加密联络线，是梅州市南部区县相互连接的交通主干线，也是梅州市高速路网的横向连接线。隧道为双洞独立式，双洞洞长分别为6332 m和6335 m，属特长隧道，洞身最大埋深达751 m，穿过地区山体植被较发育(隧道穿越洞身图见图1)。隧址区岩石大多为燕山期花岗岩及侏罗系安山玢岩。

隧址区位于深圳–五华断裂带、大埔–海丰断裂带之间，处于榕江断裂带南西侧。隧道穿越地段未见较大型断裂，但穿越小型断层数量较多(构造遥感解译示意图见图2) [5]。

Figure 1. Design cave diagram

图1. 设计洞身图

Figure 2. Tectonic remote sensing interpretation

图2. 构造遥感解译示意图

3. 现场地应力实测分析研究

现场测量地应力的方式较多，主要分为直接测量法和间接测量法两种。本项研究采用直接测量法中的水压致裂法进行地应力实测 [6]。

水压致裂法是目前较为成熟的一种测试深埋地应力的手段，具有操作简便、测量周期相对较短、测量深度大等优点。由于在计算中不需要考虑岩石材料参数，从而避免了弹性常数导致的误差。因此，选取水压致裂法在鸿图嶂隧道现场进行了地应力测试，为后文的地应力回归分析提供依据。

分别在K93+280位置钻孔CSZK14花岗岩段及K91+100位置，钻孔XSZK21安山玢岩段采用水压致裂法进行地应力测试，测试基本情况见下表1、表2 [7]。

Table 1. CSZK14 granite section ground stress test situation

表1. CSZK14花岗岩段地应力测试基本情况一览表

Table 2. XSZK21 porphyritic section ground stress test situation

表2. XSZK21玢岩段地应力测试基本情况一览表

由于现场实测的地应力资料按主应力平面方位给出的，而在地应力场回归分析中，以计算坐标系XYZ中坐标应力分量为基本对象，即正北方位与Y轴存在一定的夹角，与Y轴正向夹角为18˚，因此需要将实测地应力根据实际的夹角进行变换以获得真实数据。

将实测地应力场时的坐标系记为X、Y、Z数值仿真时的坐标系记为X’、Y’、Z’，实测坐标系与数值模拟系之间的余弦值记为 $l,m,n$，则当由现场的坐标系向模型坐标系转换时，各应力分量根据表3和式(1)进行计算。

Table 3. Directional cosine table

表3. 方向余弦表

$\left(\begin{array}{c}{\sigma }_{x^{\prime }}\\ {\sigma }_{y^{\prime }}\\ {\sigma }_{z^{\prime }}\\ {\tau }_{x^{\prime }{y}^{\prime }}\\ {\tau }_{y^{\prime }{z}^{\prime }}\\ {\tau }_{x^{\prime }{z}^{\prime }}\end{array}\right)=\left[\begin{array}{cccccc}{l}_1^2& m_1^2& n_1^2& 2l_1{m}_1& 2m_1{n}_1& 2l_1{n}_1\\ l_2^2& m_2^2& n_2^2& 2l_2{m}_2& 2m_2{n}_2& 2l_2{n}_2\\ l_3^2& m_3^2& n_3^2& 2l_3{m}_3& 2m_3{n}_3& 2l_3{n}_3\\ l_1{l}_2& m_1{m}_2& n_1{n}_2& l_1{m}_2+l_2{m}_1& m_1{n}_2+m_2{n}_1& l_2{n}_1+l_1{n}_2\\ l_2{l}_3& m_2{m}_3& n_2{n}_3& l_2{m}_3+l_3{m}_2& m_2{n}_3+m_3{n}_2& l_3{n}_2+l_2{n}_3\\ l_3{l}_1& m_3{m}_1& n_3{n}_1& l_3{m}_1+l_1{m}_3& m_3{n}_1+m_1{n}_3& l_3{n}_1+l_1{n}_3\end{array}\right]\left(\begin{array}{c}{\sigma }_x\\ {\sigma }_y\\ {\sigma }_z\\ {\tau }_{xy}\\ {\tau }_{yz}\\ {\tau }_{xz}\end{array}\right)$ (1)

转换所得模型计算坐标系下的地应力值如表4所示。

Table 4. In-situ ground stress conversion value

表4. 原位地应力转换值

注：S_x为X方向的水平应力；S_y为Y方向的水平应力；S_z为垂直主应力；S_xy为水平方向剪应力。

4. 鸿图嶂隧道初始应力场回归分析

对鸿图嶂隧道初始应力场反演回归分析的步骤如下：

1) 依据工程勘察建立ANSYS计算模型；

2) 分别施加自重应力场及五种构造应力场边界条件，提取六次模拟计算数据，将钻孔CSZK14测点数据将其作为已知数据基于多元回归方法建立多元回归方程。将钻孔XSZK21数据代入验证多元回归效果；

3) 基于最小二乘法原理，对残差平方和进行最小化，即可获得多元回归方程中各自变量(待定因素)系数的最优解，进一步可求得鸿图嶂隧道工程区的初始地应力场 [8]。

4.1. 数值模型建立

在平面图上选择了隧道线路附近约7200 × 1200 m的一块长方形区域进行模型建立，该区域包含鸿图嶂隧道全线并且范围包含了鸿图嶂隧道7条控制性断层，由于断层空间结构复杂，给数值建模及网格划分带来困难，故进行近似处理：在离隧道轴线较远处做近似直线延伸处理；断层厚度依据地勘资料各断层实际厚度范围取值，同一断层取相同平均厚度；老断层被新断层截断的外延部分在模型中不做模拟(例如图中F2-7断层)。如图3为断层近似处理示意图，图中红线为断层实际延伸线，黑色直线为近似处理延伸线。

Figure 3. Fault approximation processing map

图3. 断层近似处理示意图

Figure 4. Schematic diagram of grid division

图4. 网格划分示意图

本次三维模型离散化全部采用10节点四面体等参单元(Solid92) [9]。在划分有限元网格的时候，不同区域采用的网格精度不同，其划分规则为：在可能对结果产生较大影响的区域(如断层处)采用了更小的网格尺寸，在其余区域(如连续岩体)内网格采用了相对较大的网格尺寸。网格划分如图4所示，不同材料以颜色区分(断层较窄，图中颜色部分难以区分)。岩体材料计算参数的选取如表5所示。

Table5. Calculation parameters selection table

表5. 计算参数选取表

4.2. 边界条件及分析原理

鸿图嶂隧道水压致裂法所测得的地应力结果中可发现：隧道铅锤方向无对应的剪切应力，因此，计算中可忽略产生铅锤方向剪切应力的应力场，仅选取自重应力场和5个边界水平构造应力场作模拟分析(见图5)。

(a) 自重 (b) 东西向均匀挤压构造 (c) 南北向均匀挤压构造 (d) 东西向三角形荷载构造 (e) 南北向三角形荷载构造 (f) 水平面均匀剪切构造运动

Figure 5. Calculation of boundary conditions

图5. 计算边界条件

计算中具体的实现过程如下：

1) 自重应力场：在模型的四周及底面均加上约束，将通过表5岩体参数求得的不同岩石自重应力加于模型之上。

2) 构造应力场：均匀挤压构造及三角形荷载构造通过在底部施加约束后，于两侧施加约束，另两侧施加均布(三角形)荷载实现，剪切构造运动的实现采用模拟边界位移的方式，其中于长边界(与隧道轴线方向一致的边界)上施加0.012 m位移，在短边界(垂直于轴线方向的边界)上施加0.072 m位移。

以多元线性回归的方式设通过量测和分析获得的试验数据为：

$\left(y_i,x_{i1},x_{i2},\cdots ,x_{ip}\right),i=1,2,\cdots ,N$ (2)

式中，N为进行试验或实测的次数。

根据数理统计理论，可将上述表达式汇总为以下总式(3)：

$Y=X\beta +\epsilon$ (3)

式中，Y为观测值向量， $\beta$ 为参数列向量，X为可通过观测和分析而得到的系数矩阵，ε为N维随机向量，代表其余可能对地应力场造成影响的因素常量。

回归效果的检验采用复相关系数R即回归平方和 $S_R^2$ 在总离差平方和 $S_T^2$ 中所占比例大小衡量线性相关的密切程度。

$R=\sqrt{\frac{S_R^2}{S_T^2}}=\sqrt{\frac{{\displaystyle \underset{i=1}{\overset{N}{\sum }}{\left({\hat{y}}_i-\bar{y}\right)}^2}}{{\displaystyle \underset{i=1}{\overset{N}{\sum }}{\left(y_i-\bar{y}\right)}^2}}}$ (4)

采用偏相关系数V_i检验回归系数的显著性，计算公式为：

$V_i=\sqrt{1-\frac{S_c}{{S^{\prime }}_c}}$ (5)

4.3. 回归结果分析

采用2~5号点共4个测点实测资料(表1)和数值模拟所得地应力场结果，用最小二乘法进行多元线性回归分析 [10]，通过式(6)：

$\left[\begin{array}{cccc}{\displaystyle \underset{k=1}{\overset{m}{\sum }}{\displaystyle \underset{j=1}{\overset{6}{\sum }}{\left({\sigma }_{jk}^1\right)}^2}}& {\displaystyle \underset{k=1}{\overset{m}{\sum }}{\displaystyle \underset{j=1}{\overset{6}{\sum }}{\sigma }_{jk}^1{\sigma }_{jk}^2}}& \cdots & {\displaystyle \underset{k=1}{\overset{m}{\sum }}{\displaystyle \underset{j=1}{\overset{6}{\sum }}{\sigma }_{jk}^1{\sigma }_{jk}^n}}\\ 对& {\displaystyle \underset{k=1}{\overset{m}{\sum }}{\displaystyle \underset{j=1}{\overset{6}{\sum }}{\left({\sigma }_{jk}^2\right)}^2}}& \cdots & {\displaystyle \underset{k=1}{\overset{m}{\sum }}{\displaystyle \underset{j=1}{\overset{6}{\sum }}{\sigma }_{jk}^2{\sigma }_{jk}^n}}\\ & ⋮& ⋮& ⋮\\ & 称& & {\displaystyle \underset{k=1}{\overset{m}{\sum }}{\displaystyle \underset{j=1}{\overset{6}{\sum }}{\left({\sigma }_{jk}^n\right)}^2}}\end{array}\right]\times \left[\begin{array}{c}{b}_1\\ b_2\\ ⋮\\ b_n\end{array}\right]=\left[\begin{array}{c}{\displaystyle \underset{k=1}{\overset{m}{\sum }}{\displaystyle \underset{j=1}{\overset{6}{\sum }}{\sigma }_{jk}^1{\sigma }_{jk}}}\\ {\displaystyle \underset{k=1}{\overset{m}{\sum }}{\displaystyle \underset{j=1}{\overset{6}{\sum }}{\sigma }_{jk}^2{\sigma }_{jk}}}\\ ⋮\\ {\displaystyle \underset{k=1}{\overset{m}{\sum }}{\displaystyle \underset{j=1}{\overset{6}{\sum }}{\sigma }_{jk}^n{\sigma }_{jk}}}\end{array}\right]$ (6)

即可得到分别对应自重应力场、沿X轴方向的均布荷载引起的应力场、沿Y轴方向的均布荷载引起的应力场、沿X轴方向的三角形荷载引起的应力场、沿Y轴方向的三角形荷载引起的应力场和水平面内剪切应力的6个回归系数b₁ = 1.2466；b₂ = 49.0772；b₃ = 45.9447；b₄ = −6.6716；b₅ = −10.9647；b₆ = 175.5855，由此获得式(7)的鸿图嶂隧道的初始应力场回归公式：

$\sigma =1.2466{\sigma }_{自}+49.0772{\sigma }_{构1}+45.9447{\sigma }_{构2}-6.6716{\sigma }_{构3}-10.9647{\sigma }_{构4}+175.5855{\sigma }_{构5}$ (7)

式中：σ为初始应力场回归值；σ_自为岩体自重引起的应力场；σ_构1为由垂直X轴的面施加沿X方向的单位法向均布荷载引起的应力场；σ_构2为由垂直Y轴的面施加沿Y方向的单位法向均布荷载引起的应力场；σ_构3为由垂直X轴的面施加沿X方向的三角形荷载引起的应力场；σ_构4为由垂直Y轴的面施加沿Y方向的三角形荷载引起的应力场；σ_构5为由垂直X轴的面施加沿Y方向的剪切位移及垂直Y轴的面施加沿X方向的剪切位移引起的应力场。

通过公式(4)可以得到复相关系数r = 0.996232，其接近于1，表明反演数据与钻孔CSZK14处测值的变化规律较为吻合；回归系数显著性检验结果V₁ = 9.1605E−05；V₂ = 0.0069；V₃ = 0.0315；V₄ = 0.5814；V₅ = 0.0694；V₆ = 0.0008，其中V₄，V₅在显著性水平 $\alpha =0.05$ 上不满足显著性要求，考虑V₄，V₅代表仅有的三角形荷载回归系数显著性，故仅剔除显著性较差的回归系数σ_构4，保留显著性较好的σ_构5，重新进行回归分析。结果如下：

$\sigma =1.3519{\sigma }_{自}+40.9670{\sigma }_{构1}+44.2400{\sigma }_{构2}-10.5817{\sigma }_{构3}+173.6971{\sigma }_{构5}$ (8)

通过上述公式可以得到鸿图嶂隧道测点处的回归值，可算得测点值主应力值，提取各测点的最大值主应力方向换算为实际坐标系后如下表6所示。从表中可以看出回归值与实测值误差15%以内，反演得到的主应力变化规律与实测值相同。最大主应力角方向与实测结果较为接近，表明反演结果具有较高的可信度。

Table 6. Comparison table of main stress inversion of drill hole CSZK14

表6. 钻孔CSZK14主应力反演对比表

考虑到多元回归方程式是利用钻孔CSZK14处的测值求解得出，其回归值必然与实测值接近，从而导致其误差具有一定的局限性。因此，利用钻孔XSZK21处数据进行验证，将模型提取的应力值代入回归方程求解，得到的回归值与实测值(验证数据)进行误差比较如表7。其误差范围0.8%~24%，误差均在25%以内，符合程度较好，证明反演所得鸿图嶂隧道工程区应力场合理可信。

Table 7. Validation table of main stress inversion of drill hole XSZK21

表7. 钻孔XSZK21主应力反演验证表

5. 鸿图嶂隧道轴线地应力分布规律研究

通过提取模型中隧道沿线各分应力场数据，代入上节求得的回归方程，叠加5个子应力场，求得整个隧道轴线的地应力场，考虑到圣维南原理即存在有限元模型边界效应的影响，去除两端部分受荷载影响范围，所得隧道主力力图如下图6。

Figure 6. Main stress map from K90+200 to K94+900 at the axis of Hongtuzhang Roach Tunnel

图6. 鸿图嶂隧道轴线处K90+200~K94+900主应力图

通过前述应力分析表明，隧道区内以构造与自重应力场联合作用为主，结合上图隧道轴线主应力变化曲线图可得出隧道区内的地应力值有如下特征 [11]：

1) 鸿图嶂隧道轴线方向最大水平主应力在21 MPa以上，最大达40 MPa，根据地勘资料CSZK14钻孔岩芯饱和抗压强度90 Mpa，XSZK21钻孔岩芯饱和抗压强度80 Mpa，因此岩石强度比为2~4.29，可以判定鸿图嶂隧道为高地应力隧道。

2) 从现场实测的地应力点可以看出，鸿图嶂隧道区测点附近地应力场的分布基本以构造应力为主，在测量深度范围内三个主应力之间的大小关系依次为 $S_H>S_v>S_h$。

3) 从图中隧道低埋深区域可以发现其垂直应力远低于水平主应力的范围，说明此区域影响地应力的主要因素为构造应力。

4) 从隧道中部区域的曲线走势可以发现随着隧道埋深的增大，水平主应力和垂直应力均有不同程度的增加，其中垂直应力较水平主应力增幅更大，量值上接近最大主应力，并于K91+800附近超过最大主应力，说明在中部高埋深区域，构造应力影响不强，与垂直应力在同一量值水平；

5) 隧道全线的最小水平主应力较为稳定。

6) 从图中曲线走向和断层位置对比中可以明显发现，由于断层的影响，应力场在单个断层附近出现突变，产生了一定的应力释放现象。

7) 由于部分区域断层密集存在使得该区域的最大水平主应力要明显小于周边区域。相较于独立断层处的应力突变现象也更为明显。

6. 结语

通过上述对鸿图嶂隧道区应力场的研究，可以得出以下结论：

1) 通过水压致裂法原位地应力测试工作，从现场实测的地应力点可以看出，鸿图嶂隧道区测点附近地应力场的分布基本以构造应力为主，在测量深度范围内三个主应力之间的大小关系依次为 $S_H>S_v>S_h$。

2) 根据已有钻孔求得的鸿图嶂隧道工程区岩体的初始应力场回归公式经验证误差较小，其结果合理可信。

3) 随着隧道埋深的增大，水平主应力和垂直应力均有不同程度的增加，其中垂直应力较水平主应力增幅更大，量值上接近最大主应力，且在断层处出现应力场突变现象。并且由于部分区域断层密集存在使得该区域的最大水平主应力要明显小于周边区域。

岩体的初始应力场大小与隧道开挖过程中的围岩稳定性、地质灾害等的发生密切相关，因此加强隧道工程区的初始应力场的基础性研究，不仅对设计阶段起指导性意义，同时为后续的施工期间工程的安全性研究提供了基础性资料。

基金项目

国家自然科学基金项目(U2034205, 51878571)；广东省交通科技项目(2016-02-014)。