基于有限元分析和多元线性回归的初始地应力场反演分析及应用

Inversion Analysis and Application of Initial Geostress Field Based on Finite Element Analysis and Multivariate Regression

摘要: 本文根据玲珑金矿的工程地质条件及实测地应力资料，运用有限元和多元线性回归数值分析方法，对该矿指定区域内的初始应力场进行回归反演分析。首先，利用航拍资料建立了研究区域的有限元模型，通过三维有限元分析方法获得了多种设定工况的应力分布。然后，采用具有迭代求精的多元回归分析方法，得出最优回归系数，从而确定了初始地应力场三维回归分析模型。最后，根据反演分析结果，讨论了对于初始地应力分布的影响因素。通过对比实测点的应力计算值与现场实测值，可确定两者在量值上相当，且方向上接近，表明经过回归得到的地应力场是合理的，可以比较精确地反映该金矿指定区域内的初始地应力场，为金矿开挖过程中的岩爆灾害预测与防护提供了可用的工具。

Abstract: Based on the engineering geological conditions and measured geostress data of the Linglong Gold Mine, this paper uses the finite element and multivariate regression analysis method to carry out the regression inversion analysis of the initial geostress field in the specific region of the mine. Firstly, the finite element model of the study area was established using the aerial photography data, and the stress distribution of various setting conditions was obtained by the three-dimensional finite element analysis method. Then, using the multivariate regression analysis method with iterative refinement, the optimal regression coefficient is obtained, and the initial three-dimensional regression analysis model of geostress field is determined. Finally, according to the results of inversion analysis, the influencing factors on the initial ground stress distribution are discussed. Through the comparison between computed and measured geostress values of measuring points, it can be determined that the two values are identical in magnitude and close in direction, indicating that the measured geostress field obtained by regression is reasonable, and can reflect the initial measured geostress field in the designated area of the gold mine more accurately, which provides a useful tool for the prediction and protection of rock burst disaster during the gold mine excavation.

1. 引言

在进行岩体工程设计和稳定性分析时，初始地应力的分布特征是重要的影响因素之一 [1]。实测地应力是提供区域地应力场最直接的途径，但由于场地和经费等原因，不可能进行大量的测量 [2]。因此，需要在现场实测地应力资料的基础上，通过有效的计算分析方法，反演得到准确的适用范围较大的地应力场。

目前，国内外反演初始地应力场的方法主要可分为两类：一类是位移反分析法，这一方法是1971年由Kavanagh和Clough [3] 首先提出的，利用有限元法根据位移测值计算岩体力学性质参数。其基本思想为：根据大地构造变形观测或工程施工开挖中的位移监测，反演研究区域的应力场。这是一种间接方法，适合研究区域内缺乏地应力实测资料或实测地应力是扰动地应力的情况，该方法常用于地下工程局部范围岩体应力场的反演。另一类是应力回归分析方法，应力回归分析方法于1983年由天津大学郭怀志、马启超、薛玺成和王大年教授提出 [4]，其基本思想为：结合对区域地应力场产生条件的规律性认识，建立该区域地应力场三维地质概化模型。根据工程所在区域少量地应力实测资料进行回归计算，使得计算应力场与实测应力场达到最优拟合，以求得工程区域初始地应力场。当计算区域内已有初始地应力实测资料时，使用该方法较为高效、可靠，故该方法被广泛运用 [5] - [11]。

本文采用第二种反演分析方法，根据某金矿实测地应力资料，采用三维有限单元法，结合多元线性回归分析，进行初始地应力场的反演计算，以获得较为准确的初始地应力场，为开挖及长期稳定性分析提供了重要依据。

2. 初始地应力场反演回归分析方法

初始地应力场反演回归分析的基本思想为：1) 根据已知地质地形勘测试验资料，建立三维计算模型。2) 把可能形成初始地应力场的因素(自重、构造运动和温度等)作为待定因素，对每一种待定因素用数值计算获得已知点位置的应力值，然后在每一种待定因素计算的应力值与已知点实测地应力值之间，建立多元回归方程。3) 用统计分析方法(最小二乘法)，根据残差平方和最小的原则，求得回归方程中各自变量(待定因素)系数的最优解，同时，在求解过程中，可对各待定因素进行筛选，贡献显著的引进，不显著的剔除，从而获得区域初始地应力场分布规律。

2.1. 有限元模型的建立和基本应力场的模拟

首先，界定考察域的范围，范围选取得是否合理直接影响到计算结果是否准确。对于一定的工程考察域而言，初始应力场是客观存在的，是唯一的、确定的，如果取域不当，将会导致应力场的失真。因此，考察域的有效范围要能真实反应地应力场的地形地质条件。一般来说，考察域的边界宜选在有地形特征的部位。

初始地应力场主要由自重应力场和地质构造运动应力场两个部分组成，并通过三维有限元模型计算来模拟。首先，根据地应力实测资料、金矿地形地质条件、岩体的力学性能等因素，建立考察域的三维有限元计算模型。对自重应力的模拟，采用岩体实测密度进行计算。构造应力的模拟较为复杂，一般视为是以下5个子构造应力场的线性叠加：1) 东西向水平均匀挤压构造运动(本文计算坐标系中的x轴向)。2) 南北向水平均匀挤压构造运动(本文计算坐标系中的y轴向)。3) 水平面内的均匀剪切变形构造运动(本文计算坐标系中的xy平面)。4) 东西向垂直平面内的竖向均匀剪切变形构造运动(本文计算坐标系中的xz平面)。5) 南北向垂直平面内的竖向均匀剪切变形构造运动(本文计算坐标系中的yz平面)。实际所在面沿不同方向分布荷载的具体值是未知的，本文采用直接对计算域边界加相应单位分布位移，利用ANSYS程序对各种因素单独作用下的三维模型进行计算，从而得到6个基本应力场即自重应力场和5个子构造应力。

最后，通过形函数插值，提取各个基本应力场在地应力测点处的6个空间应力分量。若有N个实测点，那么将会产生6 × N组应力分量，作为数学回归模型的“观测值”。

2.2. 多元线性回归模型

根据多元回归法原理，将地应力回归计算值 ${\hat{\sigma }}_k$ 作为有6个应力分量的实测值视为因变量，把有限元计算求得的自重应力场和构造应力场相应于实测点的应力计算值 ${\sigma }_k^i$ 作为自变量(相应有6个应力分量)，则回归方程的形式为：

${\hat{\sigma }}_k={\displaystyle \underset{i=1}{\overset{n}{\sum }}{L}_i\cdot {\sigma }_k^i}$ (1)

式中：k为测点的序号， ${\hat{\sigma }}_k$ 为第k观测点的回归计算值，L_i为相应于自变量的多元回归系数，n为工况数。

对于每一个应力状态 ${\sigma }_k^i$ 都可确定一个回归计算值 ${\hat{\sigma }}_k$。用量测值 ${\sigma }_k^{*}$ 与回归计算值 ${\hat{\sigma }}_k$ 之差来表示量测值与回归方程的偏离程度，即该量测值的残差 $S_k={\sigma }_k^{*}-{\hat{\sigma }}_k$。假设有m个测点，则全部量测值与回归方程的偏离程度，用全部量测值 ${\displaystyle \underset{k=1}{\overset{m}{\sum }}{\displaystyle \underset{j=1}{\overset{6}{\sum }}{\sigma }_{jk}^{*}}}$ 与回归计算值 ${\displaystyle \underset{k=1}{\overset{m}{\sum }}{\displaystyle \underset{j=1}{\overset{6}{\sum }}{\sigma }_{jk}^i}}$ 的残差平方和表示，即：

$S_k={\displaystyle \underset{k=1}{\overset{m}{\sum }}{\displaystyle \underset{j=1}{\overset{6}{\sum }}{S}_{jk}^2}}={\displaystyle \underset{k=1}{\overset{m}{\sum }}{{\displaystyle \underset{j=1}{\overset{6}{\sum }}\left({\sigma }_{jk}^{*}-{\displaystyle \underset{i=1}{\overset{n}{\sum }}{L}_i\cdot {\sigma }_{jk}^i}\right)}}^2}$ (2)

式中： ${\sigma }_{jk}^{*}$ 为k测点j应力分量的量测值， ${\sigma }_{jk}^i$ 为i工况下k观测点j应力分量的有限元计算值。

根据最小二乘法原理，使得S_k为最小值的方程式为 $\frac{\partial S_k}{\partial L_i}=0\left(i=1,\cdots ,n\right)$，即：

$\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\cdot {\sigma }_{jk}^2}}& \cdots & {\displaystyle \underset{k=1}{\overset{m}{\sum }}{\displaystyle \underset{j=1}{\overset{6}{\sum }}{\sigma }_{jk}^1\cdot {\sigma }_{jk}^n}}\\ {\displaystyle \underset{k=1}{\overset{m}{\sum }}{\displaystyle \underset{j=1}{\overset{6}{\sum }}{\sigma }_{jk}^2\cdot {\sigma }_{jk}^1}}& {\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\cdot {\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]\cdot \left\{\begin{array}{c}{L}_1\\ L_2\\ ⋮\\ L_n\end{array}\right\}=\left\{\begin{array}{c}{\displaystyle \underset{k=1}{\overset{m}{\sum }}{\displaystyle \underset{j=1}{\overset{6}{\sum }}{\sigma }_{jk}^{*}\cdot {\sigma }_{jk}^1}}\\ {\displaystyle \underset{k=1}{\overset{m}{\sum }}{\displaystyle \underset{j=1}{\overset{6}{\sum }}{\sigma }_{jk}^{*}\cdot {\sigma }_{jk}^2}}\\ ⋮\\ {\displaystyle \underset{k=1}{\overset{m}{\sum }}{\displaystyle \underset{j=1}{\overset{6}{\sum }}{\sigma }_{jk}^{*}\cdot {\sigma }_{jk}^n}}\end{array}\right\}$ (3)

解此方程式(3)，可得n个待定回归系数 $L={\left[L_1,L_2,\cdots ,L_n\right]}^{\text{T}}$，则计算域内任一点P的回归初始应力，

由该点各工况有限元计算值迭加可得：

${\sigma }_{jk}={\displaystyle \underset{i=1}{\overset{n}{\sum }}{L}_i\cdot {\sigma }_{jp}^i}$ (4)

3. 工程实例

某金矿是我国最大的黄金矿山之一，近些年，该金矿开采速度迅速增加，目前开采深度已到1000 m，随着开采深度不断增加，矿体赋存状况、地质条件和围岩稳定性状态显著恶化，地压显现加剧。为了提出符合深部矿体赋存条件的高效、安全、可靠的开采设计方法，必须进行矿床深部岩石力学的综合调查、研究，地应力场的反演计算是其中最重要的研究内容之一。

3.1. 计算范围与三维地质模型

计算范围的确定应遵循以下2个原则：1) 几何范围必须包含全部工程影响区域，且适当增大，减少边界影响。2) 边界处的几何约束条件必须易于确定。通过分析金矿工程区的范围和该处工程地质以及地应力实测点的分布资料，确定计算区为：x，y轴的计算范围为1500 m × 1000 m；z轴竖直向上，从最低点向上1250 m。采用8节点等参单元，共划分单元1563095个、节点276334个。计算网格见图1。矿区岩体物理力学计算参数见表1。

Figure 1. Meshes for numerical calculation

图1. 数值计算网格图

Table 1. Mechanical paremeters of rock masses

表1. 岩石力学参数

3.2. 现场地应力测试

该金矿地应力测量共布置6个测点，分别位于第1测点：19中段，(−570 m)，第2测点：19中段，(−570 m)，第3测点：22中段，(−690 m)，第4测点：22中段，(−690 m)，第5测点：26中段，(−850 m)，第6测点：26中段，(−850 m)，分别位于−570、−690和−850 m三个中段水平上，每个水平2个测点。共计6个测点，见表2。

Table 2. The location of geostress measurement points in calculating reference frame

表2. 计算坐标下地应力测量点位置

根据主应力与全应力分量之间的转换公式和大地坐标系与任意空间坐标系之间的全应力转换公式 [12]，将表3中的主应力值转换成计算坐标系中的应力分量，其计算结果见表4。

Table 3. Table of measurement results of 3D geostress

表3. 三维地应力测量结果表

Table 4. The measured geostresses component values of each measuring point in calculating reference frame

表4. 计算坐标下各测点实测地应力分量值

3.3. 基本应力场的模拟

利用通用有限元软件ANSYS对各个基本应力场分别进行模拟计算。1) 自重应力场。采用岩体实测密度，计算在自重作用下形成的自重应力场，考察域模型的4个侧面和底面均加法向位移约束。2) x轴向挤压构造运动见图2，计算时以边界位移为控制对象，控制单位取为1 cm；3) y轴向挤压构造运动见图3，计算时仍以边界位移为控制对象，控制单位取为1 cm；4) xz平面剪切构造运动见图4，边坡面自由，边界剪切位移选为10 cm；5) yz平面剪切构造运动见图5，边坡面自由，边界剪切位移也选为10 cm。

Figure 2. Extrusion geologic techonics movement of x-direction

图2. x轴向挤压构造运动

Figure 3. Extrusion geologic techonics movement of y-direction

图3. y轴向挤压构造运动

Figure 4. Shearing geologic techonics movement of xz-plane

图4. xz平面剪切构造运动

Figure 5. Shearing geologic techonics movement of yz-plane

图5. yz平面剪切构造运动

本次利用有限元程序分别对上述5种因素单独作用下的三维数值模型进行计算，根据地应力各测点的回归计算值，建立初始地应力场与各种影响因素之间的回归方程为：

${\sigma }_{地}=L_{自}{\sigma }_{自}+L_1{\sigma }_{构1}+L_2{\sigma }_{构2}+L_3{\sigma }_{构3}+L_4{\sigma }_{构4}+e_k$ (5)

式中：σ_地为初始地应力；σ_自为由自重引起的应力；σ_构1为yz平面施加沿x轴向的1 cm单位均匀挤压位移引起的应力；σ_构2为xz平面施加沿y轴向的1 cm单位均匀挤压位移引起的应力；σ_构3为yz平面上施加沿y轴向的1 cm切向分布单位位移引起的应力；σ_构4为xz平面上施加沿x轴向的1 cm切向分布单位位移引起的应力；σ_构5为yz平面上施加沿z轴向的10 cm切向分布位移引起的应力；σ_构6为xz平面上施加沿z轴向的10 cm切向分布位移引起的应力；e_k为随机变量；L_自、L₁、L₂、L3和L₄均为回归系数。

3.4. 地应力场反演计算结果分析

先对式(5)确定的5种计算工况利用ANSYS分别进行计算，然后以实测地应力资料(见表4)中的6个测点的应力为回归目标。利用多元线性回归方法得到回归系数：L_自 = 1.1989，L₁ = 186.37，L₂ = 43.102，L₃ = −5.135，L₄ = −2.1504，e_k = −0.5345，则该金矿研究区域初始地应力场回归方程为：

${\sigma }_{地}=1.1989{\sigma }_{自}+186.37{\sigma }_{构1}+5.135{\sigma }_{构2}-5.135{\sigma }_{构3}-2.1504{\sigma }_{构4}-0.5345$ (6)

回归分析中的复相关系数为0.90167，表明回归公式的相关性好。

经应力场平衡计算，得到计算坐标系下测点的地应力实测值与回归值对比见图6和表5；测点的实测与回归主应力值对比见表6。

由表5和图6可知，大部分回归值与实测值之间的绝对误差小于6.73 MPa，有些测点有些应力值甚至接近一致。判断应力场回归优劣的又一重要指标是计算平衡后模型内的最大主应力和最小主应力的数值，从表6和图7可以看出，回归主应力与实测主应力值绝大部分在数值上很接近。因此，回归得到的该金矿研究矿区的应力场是合理的。

取研究区域的最大主应力及最小主应力等值线如图8所示。

由图8可知，说明该金矿的初始地应力以自重应力为主，断裂带对初始地应力场有一点的影响。

Table 5. Comparison of measured and regressive geostress components of each measurement point

表5. 各测点实测与回归应力分量对比

Table 6. Comparison of measured and regressive principal stresses of measuring points

表6. 测点的实测与回归主应力对比

Figure 6. Comparison histogram of measured and calculated geostress values of measuring points in calculating reference frame

图6. 计算坐标系下测点的地应力实测值与回归值对比直方图

Figure 7. Comparison histogram of measured and calculated principal stress values of measuring points

图7. 测点的主应力实测值与回归值对比直方图

4. 结语

1) 根据金矿矿区地应力实测资料，利用该模型对矿区初始地应力场进行反演分析。通过实测点的计算应力值与现场实测值的比较，两者在量值和方向上接近，表明该方法较好地反映了实际地应力场的分布规律，能够满足矿区设计和研究的需要。表5、表6和图6、图7的比较分析结果表明，采用多元线性回归与三维地质模型的有限元计算相结合的初始地应力反演回归方法能够得到合理的应力场分布，反演回归应力值与实测应力值拟合较好，为金矿长期稳定性分析提供了合理的三维初始地应力场。同时，也表明该反演回归方法是一种高效、准确、可靠的方法。

2) 回归结果的好坏与测点的实测应力值密切相关，因此，实测点数据要尽可能准确，否则会影响回归效果。另外，还与回归系数的求解方法密切相关。不同的求解方法如最小二乘法和偏最小二乘法所求得的系数不同，因而对回归结果产生一定的影响，但总的来说，不会影响地应力场分布规律的判断。

3) 通过本文的尝试，探索了一条通过ANSYS软件实现边界构造运动组合条件的施加来回归地应力场的方法。

Figure 8. Contour of principal stresses (unit: Pa)

图8. 主应力等值线(单位：Pa)

基金项目

北京建筑大学校设科学研究基金资助项目，项目编号：00331614048。

参考文献