0   引言

地应力是存在于地壳中未被人工扰动的天然应力，其在指导矿山地下开采工程、土木建设工程、边坡修建等工程中的作用越来越重要，是各种工程开挖设计必须要做好的前提工作.国内外对获取初始地应力进行了大量研究，并取得了许多研究成果^[1-2].目前初始地应力的计算方法大致有以下几种：边界荷载调整、应力函数法、有限元数学模型回归分析法、三维有限元反演分析法等^[3-5].在地应力反演方法上，主要采用的是以分析开挖工程区域小范围内的初始地应力场的数学计算方法，包括：位移反分析方法，由监测位移反演局部区域的应力分布^[6]；应力反演方法，由现场有限个点的地应力实测值，通过应力函数或数值计算方法回归分析得到研究区域的应力.

本文在前人的研究成果基础上，假定初始地应力场主要由自重应力场和地质构造运动应力场两部分组成，将各运动模式分解成应力基本运动模式，并结合多元线性回归原理，采用Fortran编写相关程序，以及FLAC3D强大的数值模拟分析功能^[7]，以此建立地应力场反演多元线性回归模型，发展和完善地应力场反演方法.

1   多元线性回归反演理论

多元线性回归法反演初始地应力场是当前常用且比较精确的方法之一^[8]，其基本思路是：通过在岩体的局部位置进行实际测量得到一定的地应力值，然后根据应力与荷载之间存在的相互关系，分析研究得到岩体的应力测量信息，再采用多元线性回归法展开计算分析，以此来反演地应力场.

1.1   基本影响因素模拟

由于初始地应力主要受岩体自重和地质构造运动两方面影响，进而把初始地应力场认为主要由自重应力场和地质构造运动应力场两个部分组成^[9]，通过用FLAC^3D对自重应力场和地质构造应力场进行数值模拟.通常情况下，地质构造运动应力场由以下5种子构造应力模式线性叠加而成：

（1）左右方向水平挤压构造运动产生的应力场；

（2）前后方向水平方向水平挤压构造运动产生的应力场；

（3）水平面内剪切构造运动产生的应力场；

（4）左右垂直平面内的竖向剪切构造运动产生的应力场；

（5）前后垂直平面内的竖向剪切构造运动产生的应力场.

通过以上分析，就可以从影响地应力的众多因素中提取出6个主要因素，便于建立模型.

1.2   多元线性回归模型

基于多元线性回归原理，把所要求的初始地应力场作为因变量，形成地应力的各种基本因素作用作为自变量，采用弹性工作状态下的线性叠加原理写出初始应力场的回归模型^{[6, 10]}为：

式(1)即为回归方程，C_i为回归系数，n代表自重和构造运动基本运动模式的种类(数值模拟时亦可称为工况).式（1）中，σ_jk⁰为k测点j应力分量的初始地应力计算值；C_i（i=0, 1, 2, …, 6）为7个待估参数，分别代表每类运动模式（自重运动模式及子构造运动模式）引起的应力分量大小；σ_jk⁰、σ_jk¹、σ_jk²、σ_jk³、σ_jk⁴、σ_jk⁵、σ_jk⁶分别为自重和5个子构造运动作用下的k测点j应力分量；ε为模型误差.

上述回归模型的基本假定是：

（1）ε是因变量σ_jk⁰的误差，是相互独立的随机变量，没有系统性，其数学期望全为零，即E（ε_n）=0；

（2）每次观测相互独立，并有相同的精度，即ε_n之间的协方差关系式为

（3）ε_n服从正态分布.

1.3   回归系数C_i的确定

在使用回归方法反演区域应力场时，记第k个测点的实测地应力j分量值为，(j=1, 2, …，L；k=1, 2, …，m)，m代表地应力总测点的个数；L为应力分量个数，平面问题L=3，空间问题L=6.回归的目的是使m组实测值和数学模型计算值σ_jk⁰相等，对于每一个应力状态σ_jkⁱ，可确定一个回归计算值σ_jk⁰，观测值与回归值的偏离程度用观测值与回归计算值σ_jk⁰之差，即残差ε=ε_jk= -σ_jk⁰表示.则由式(1)知第k个地应力测点的第j个应力分量的实测值与对应测点的计算值残差平方和为

运用最小二乘法原理，使得残差平方和为最小值，即式(2)对C_i取偏导数并令其为零，则

对式(3)整理后得回归系数C_i的线性方程组矩阵：

由此方程组可求解出回归系数，它的解是唯一的，可以得到n+1个待定回归系数C=（C₀，C₁, C₂，…，C_n）^T.根据各基本运动模式（计算工况）应力回归方程，反演研究区域内任意点的应力，从而得到反演后区域内的整体应力场空间分布.

2   算例分析

为了验证上述所提出的初始地应力场反演与构建方法的可行性，利用Fortran语言编写了多元线性回归初始地应力场反演程序，并基于Flac3D有限差分软件，对一个简单平面应变问题的初始地应力场进行反演与构建.

如图 1所示为一地质剖面计算模型，模型左侧施加P=3.0MPa的均布荷载，右侧水平约束底部垂直约束，材料均质，弹性模量E=28GPa，泊松比μ=0.23，容重r=2.8×10³kg / m³.

图  1  地应力反演算例模型示意图

下载: 全尺寸图片 幻灯片

算例在方案设计时，选取如图 1所示计算模型计算区域内的6个点（A, B，…，F）作为假想的初始地应力实测点，各点在模型中相对应的单元编号分别为1070、1080、670、732、320、332.首先在Flac3D中按算例模型边界条件及自重条件下计算得到如表 1所示6个观测点的应力值，这些应力值作为各测点的初始地应力实测应力值；然后将自重因素视为未知因素，构造运动基本运动模式采用应力边界条件基本模式，分别计算单位自重、模型左侧边界单位均匀分布载荷、三角形分布载荷和剪应力分布载荷等4种工况的基本运动模式的应力场，再由这些已知的地应力值，利用所编写的多元线性回归初始地应力场反演程序，得到各工况基本运动模式的回归系数.至此，可根据各基本运动模式（计算工况）应力回归方程，反演研究区域内任意点的应力，从而得到反演后区域内的整体应力场空间分布.

表  1  地应力实测值

下载: 导出CSV 

| 显示表格

表 2列出了各运动模式单位荷载作用下应力回归系数，实际只考虑了自重和均布左侧边界加载，对模型的外部实际加载与反演加载相对误差分别为0.103%和1.4433%，计算表明反演非常成功，尤其是自重因素对形成初始地应力影响极其稳定.

表  2  各运动模式单位荷载作用下应力回归系数

下载: 导出CSV 

| 显示表格

图 2列出了各测点实测应力值与反演应力值，通过比较可以看出，各测点σ_x和σ_z实测应力值与反演应力值相对误差极小，均在0.5%以下，而剪应力的实测应力值与反演应力值误差略大一些，但相对误差也均在2.2%以下，由于区域的应力以水平和竖直方向为主，从工程的角度而言是可以接受的.

图  2  测点实测应力值与计算应力值相对误差

下载: 全尺寸图片 幻灯片

由此说明：反演应力与“实测应力”吻合较好，回归效果理想，采用多元回归方法反演线弹性地应力场是可行的，并且有较高的精度.

3   结论

（1）基于多元线性回归原理，把所要求的初始地应力场作为因变量，形成地应力的自重基本因素和构造运动的子构造应力基本因素作为自变量，采用线性叠加原理形成因变量与自变量的表达关系式，建立了多元线性回归地应力场反演数学模型，完善了多元线性回归地应力场反演方法.该方法保证了解的唯一性，使引进和剔除形成初始地应力场的各种因素有了根据，有利于对形成初始地应力复杂因素的认识.

（2）多元线性回归地应力场反演方法数学推导严密、物理意义明确，直接针对形成应力场的物理成因展开模拟，考虑了自重、挤压构造运动及剪切变形构造运动等共同作用.标准算例验证表明，计算应力值与观测应力值极其吻合，相对误差均在2.2%以下，回归效果理想，说明本文的地应力场反演方法可行，计算程序可靠，并且有较高的计算精度.