基于实测结果及线性规划法的地应力场DDM反演方法

摘要

本发明提供一种基于实测结果及线性规划法的地应力场DDM反演方法，该方法依据工程现场工程地质勘查报告，并通过实验获得所需参数信息；依据得到的信息建立位移不连续法DDM数值模型；分别在竖向施加重力荷载、施加水平x方向单位荷载、施加水平z方向单位荷载，经过DDM处理，得到区域应力场，并提取地应力实测点的应力数据处理值；依据得到的各测点的应力数据处理值，通过线性规划法处理各工况对应的回归系数，然后依据回归系数对各工况处理得到的应力场进行组合，求和得到反演地应力场。

说明书

技术领域

本发明涉及的是隧道和地下空间工程领域中的地应力场反演方法，具体是一种采用位移不连续法（DDM）作为计算工具的基于实测结果及线性规划法的地应力场反演方法。

背景技术

地应力是天然地质体（岩土体）中的未受工程扰动的应力，也称岩体初始应力或原岩应力。地应力的组成主要有五个方面：岩体自重、地质构造运动、地形、剥蚀作用和封闭应力。

在隧道、水电站地下硐室和地下空间工程中，工程区域的地应力是设计与施工必须掌握的工程条件。地下工程的岩体开挖是在地应力作用下进行的，开挖卸载将引起岩体变形与应力重分布，严重时会发生失稳、垮塌和破坏。因此，工程区域地应力的获取关系到地下工程安全与稳定。

地应力的实测方法主要有水压致裂法、应力解除法、应力恢复法和声发射法等，但这些实测方法获得的地应力是具体某一个点的地应力，而整个工程区域的地应力场需要依靠若干个点的实测地应力值，通过一定的数值方法进行反演获得。

既有的地应力场反演方法有数理统计、灰色建模、人工智能等方法，所依靠的数值计算工具一般为有限元和有限差分法。这些方法的不足主要有两点：一是反演方法复杂，计算过程繁琐，在数学上进行了优化反演，但力学概念上并不清晰；二是所采用的计算工具，如有限元和有线差分法，需要建立完整的三维实体模型，而工程区域的地形和地质条件可能十分复杂，建模可能十分困难，且三维实体模型需要划分大量的单元，对计算机的运算速度与能力要求较高。

如中国专利公开号为1601304A、申请号为200410043992.2的发明专利申请，该专利公开偶极横波测井地应力多频反演方法，其特征在于：将交叉偶极阵列声波测井数据进行角度旋转后得到两主时间序列，在此基础上提取两主方向的弯曲波频散曲线，利用两个主偏振方向弯曲波频散曲线上的n个频率点的数据根据公式(5)和公式(6)建立n个等式，其中公式(5)和公式(6)均由井孔模式波声弹性方程反演而来，其中n大于或等于6，通过求解这n个等式就可以得到地层的最大主应力SH和最小主应力Sh。但该专利方法也存在着明显的不足：首先，该方法通过声发射法测量岩石传播波速反演地应力，仅能获得测点处的地应力值，无法得到整个工程区域的地应力分布情况；其次，该技术得到的地应力值仅有两个水平最大和最小地应力值，而真实的地应力是三维条件下的，竖向地应力也是工程必需的。而本发明在这两点上具有明显的优势，反演得到的地应力场涵盖整个工程区域，对于区域中任意点的值均可获得，而且得到的地应力场为三维空间分布，涵盖水平和竖向地应力，为工程提供了完整的地应力资料。

发明内容

针对现有技术中的缺陷，本发明的目的是提供一种采用DDM作为计算工具的基于实测结果与线性规划法的地应力场反演方法，克服上述技术背景所缺乏的如何方便建模、有效计算的地应力场反演问题。

本发明是通过以下技术方案实现的，包括如下步骤：

第一步，依据工程现场工程地质勘查报告，并通过实验获得所需参数信息，包括地形地貌、各断层的位置、地应力测点的位置、各测点的地应力实测值、岩石密度ρ、岩石弹性模量E、岩石泊松比v、各断层的粘聚力c和内摩擦角

地应力实测采用孔径变形法，通过三孔交汇获得三向应力状态。在岩体中钻取3个不同方向（相邻孔夹角不小于30°）的交汇孔，对3个孔分别进行孔径变形测量，最后获得岩体三维应力状态。

其余参数可直接采用地质勘查报告中的数据，但岩石弹性模量E须进行现场实验测定：在清除爆破松动影响的岩体表面后加工试点，对粘贴承压板（直径为50.5mm的圆形钢板）处的岩石表面进行人工凿平。岩面起伏差控制在5mm内。承压板以外影响带的岩体表面加工至大致平整，并清除松动的岩块和碎石。一般液压千斤顶加载。加载方式为逐级一次循环法，载荷分5级施加，按预期岩体强度进行均分。根据现场试验实测变形资料，绘制压力～变形关系曲线，并分析岩体变形性质、特征，依据弹性理论，按半无限空间弹性体公式求得岩体变形（弹性）模量。

$>E=\frac{\pi }{4}·\frac{(1-v^2)\mathrm{PD}}{W}$>

式中：

E—岩石弹性模量（GPa）；

W-岩石变形（cm）；

P—按承压板面积计算的压力（MPa）；

D—承压板直径（cm）；

v—泊松比。

第二步，依据第一步的信息，建立位移不连续法（DDM）数值模型。

所述的位移不连续法（DDM）数值模型：在工程区域的地表曲面划分若干单元，在各断层面划分若干单元，并输入岩石密度、岩石弹性模量、岩石泊松比、各断层的粘聚力和内摩擦角。

具体的操作为：首先采用几何建模模块画出地表轮廓线，并沿延长线拖拽成曲面，断层面一般为平面，直接画出几何边界即可；然后采用网格划分模块，必须首先对地表面进行网格划分，在曲面中布点并给点编号，每相邻三点组成一个三角形单元，给每个单元编号，每个单元三个点的排序必须为逆时针，地表单元网格划分完成后，依次对各断层进行网格划分，顺序不可改变；然后在参数输入模块输入各参数，同样，必须先输入地表参数，一般地表处粘聚力和内摩擦角为0，各断层粘聚力和内摩擦角须依次输入，顺序不可改变，最后输入岩石密度、弹性模量和泊松比；此时可启动模型信息生成模块，输出可用于数据处理的计算文件，其中包含了模型几何信息、网格信息和参数信息。

所述计算文件是一个文本文件，用于后续的数据处理，之后的反演等处理过程直接读取该文件内的信息，该文件中数据排列与顺序必须保持不变，计算文件中依次包含了岩石的力学参数、断层的数目及其力学参数、每个网格节点的坐标、每个网格单元的点号以及待处理区域的范围。

第三步，在竖向施加重力荷载，经过DDM处理，得到区域应力场，并提取地应力实测点的应力数据处理值。

具体的操作为：利用第二步生成的计算文件，在DDM数据处理模块中输入荷载，重力荷载的施加须指明方向，在地表处取为0，沿模型竖向逐渐增大；启动DDM数据处理模块，待处理完成后得到结果文件；利用结果查看功能，输入实测点坐标，得到该点处的地应力数据处理值。

所述的DDM数据处理模块，采用一种人机交互的数据处理过程与方法，要求提供一个可供读取的输入文本文件，文件中的信息必须按照严格的格式与顺序排列，依次为岩石的力学参数、断层的数目及其力学参数、每个网格节点的坐标、每个网格单元的点号以及待处理区域的范围，将该文本文件提供给DDM数据处理模块，并在该模块中输入荷载，即已满足模块运行条件；模块的数据处理是依据位移不连续法，以断层单元的位移不连续量为待求解的基本未知量，依据输入的断层信息数据，将其处理为矩阵数值信息，通过应力影响系数求得边界应力值，该边界应力值符合给定的荷载边界应力条件，对输入的荷载值进行坐标转换处理，由此进行内核分析处理，得到基本未知量，继而完成输入数据的处理，得到一个结果文件；结果文件是DDM数据处理模块完成后生成的一个文本文件，由两部分内容组成：待处理的工程区域的各个点的地应力数据处理值（分为6个应力分量表达）、待处理区域的各个点的位移值（包含大小和方向）。

第四步，施加水平x方向单位荷载，经过DDM处理，得到区域应力场，并提取地应力实测点的应力数据处理值。

具体的操作为：利用第二步生成的计算文件，在DDM数据处理模块中输入x方向单位荷载，此处的荷载为应力荷载，不需指出方向，在地表处可按坐标转换进行应力分量计算处理并自动施加；启动DDM数据处理模块，待处理完成后得到结果文件；利用结果查看功能，输入实测点坐标，得到该点处的地应力数据处理值。

第五步，施加水平z方向单位荷载，经过DDM处理，得到区域应力场，并提取地应力实测点的应力数据处理值。

具体的操作为：利用第二步生成的计算文件，在DDM数据处理模块中输入y方向单位荷载，此处的荷载为应力荷载，不需指出方向，在地表处可按坐标转换进行应力分量计算处理并自动施加；启动DDM处理模块，待处理完成后得到结果文件；利用结果查看功能，输入实测点坐标，得到该点处的地应力数据处理值。

第六步，依据建立的DDM数值模型，进行数据处理，并反演求和。

依据第三步至第五步得到的各测点的应力数据处理值，通过线性规划法得到各工况对应的回归系数，然后依据回归系数对各工况处理得到的应力场进行组合，求和得到反演地应力场。

将第三步至第五步得到的三组实测点的地应力数据处理值依次输入到回归数据处理模块，然后输入实测地应力值和测点数目，各组输入数据的格式须统一为：测点编号、x方向正应力、y方向正应力、z方向正应力、xy剪应力、yz剪应力和zx剪应力，中间用空格间隔。

待求解回归系数的方程：

其中，m为测点数目，为k观测点j应力分量的观测值；为i工况下k观测点j应力分量的数值模型（DDM）处理值；应力分量j=1,2,…,6，对应6个初始应力分量。

启动回归数据处理模块求解以上方程组，得n个待定回归系数L=(L₁,L₂,…,L_n)^T。在反演求和模块输入第二至第五步的地应力场数据处理结果，并输入n个回归系数，回归系数与输入地应力场须一一对应，输入地应力场的格式须为：x方向正应力、y方向正应力、z方向正应力、xy剪应力、yz剪应力和zx剪应力，中间用空格间隔。启动反演求和模块，这样计算域内任一点p的回归初始应力可由该点各工况数值模型（DDM）处理值叠加而得

$>{\sigma }_{\mathrm{jp}}=\underset{i=1}{\overset{n}{\Sigma }}{L}_i{\sigma }_{\mathrm{jp}}^i$>

其中，j对应初始地应力的各个分量，j=1,2,…,6。

所述的回归数据处理模块，依据输入的若干组实测地应力值，进行回归处理，得到各工况对应的回归系数。需要输入实测地应力值和测点数目以及各工况的经DDM数据处理模块处理得到的对应于实测点的应力结果，各组输入数据的格式须统一为：测点编号、x方向正应力、y方向正应力、z方向正应力、xy剪应力、yz剪应力和zx剪应力，中间用空格间隔；回归处理过程依据DDM数据处理模块得到的各工况在实测点的应力结果与实测结果的相关关系，给予不同的荷载工况各自不同的权重，使各种工况组合形成的应力值与实测应力值最接近；经过处理得到的结果即为各工况对应的回归系数，可用于后续的反演处理。

所述的反演求和模块，依据输入的各工况在整个区域的经DDM数据处理模块得到的地应力结果，与经回归数据处理模块得到的各工况对应的回归系数，进行反演求和处理，将各工况分别得到的地应力场按照回归结果进行叠加，组合得到最终的反演地应力场。

本发明采用DDM作为数据处理工具的基于实测结果与线性规划法的地应力场反演方法，线性规划法将待反演目标按组成因素不同进行分解，再依据目标函数进行回归，求得各组分的回归系数，再叠加得到反演目标，具有清晰的力学概念，且数学处理简单。位移不连续法（DDM）属于间接边界元法，具有降维的特点，工程区域的三维模型只需在地表和断层等结构面上划分单元，具有建模简单和计算效率高的特点。

与现有技术相比，本发明具有如下的有益效果：

本发明首次利用DDM数据处理模块、基于线性规划法的回归数据处理模块与反演求和模块，是一种简易高效的地应力反演方法。本发明所采用的地应力实测值采用三孔交汇孔径变形法测得，可得到三维状态分布的地应力，岩石的弹性模量通过现场承压板试验实际测得，可反应区域岩体的整体力学性质，采用DDM作为数据处理工具，模型简单，几何建模模块和网格划分模块可以方便的将工程现场的物理模型建立为DDM所需的数值模型，DDM数据处理模块的计算效率优于其他计算手段，可快速得到计算结果，回归数据处理模块只需求解一个方程组，得到一组表达十分简单的回归系数，最后采用线性规划法的反演求和模块一次性得到最终反演地应力场，整个操作过程具有标准的流程化特点，易于实际应用。

附图说明

通过阅读参照以下附图对非限制性实施例所作的详细描述，本发明的其它特征、目的和优点将会变得更明显：

图1为DDM数值模型图；

图2为待开挖区域竖向地应力场图；

图3为待开挖区域水平向地应力（左右）图；

图4为待开挖区域水平向地应力（垂直纸面）图。

具体实施方式

下面对本发明的实施例作详细说明，本实施例在以本发明技术方案为前提下进行实施，给出了详细的实施方式和具体的操作过程，但本发明的适用范围不限于下述的实施例。

实施例

某水电站工程地下厂房硐室群处于河谷一侧山体中，山体表面坡度为55°，硐室区域水平埋深在230m~430m，竖向埋深在285m~480m，岩石为具有较高强度的中粗粒花岗岩。工程区域主要受4条断层影响，编号为fc-1、fc-11、fc-13和fc-23。有6个地应力实测点。拟采用基于实测结果及线性规划法的DDM反演方法获得地应力场，本实施例中没有详细说明的部分，按照发明内容描述的技术方案进行操作，发明内容中没有说明的部分，采用本领域常规操作进行。

本实施例具体施工步骤简述如下：

第一步，依据工程现场工程地质勘查报告，可确定地形地貌、断层位置等工程资料，该工程地表为55°倾角的坡面，有4条断层（fc-1、fc-11、fc-13和fc-23）,以及岩石密度ρ、岩石泊松比v、各断层的粘聚力c和内摩擦角详见表1。

岩石弹性模量的现场实验测试：在清除爆破松动影响的岩体表面后加工试点，对粘贴承压板（直径为50.5mm的圆形钢板）处的岩石表面进行人工凿平。岩面起伏差控制在5mm内。承压板以外影响带的岩体表面加工至大致平整，并清除松动的岩块和碎石。一般液压千斤顶加载。加载方式为逐级一次循环法，载荷分5级施加，按预期岩体强度30MPa进行均分。根据现场试验实测变形资料，绘制压力～变形关系曲线，并分析岩体变形性质、特征，依据弹性理论，按半无限空间弹性体公式求解岩体变形（弹性）模量。

$>E=\frac{\pi }{4}·\frac{(1-v^2)\mathrm{PD}}{W}$>

式中：E—岩石弹性模量（GPa）；W—岩石变形（cm）；P-按承压板面积计算的压力（MPa）；D—承压板直径（cm）；v-泊松比。

计算得到的弹性模量见表1.

表1  岩石与断层参数

为获得现场地应力，采用孔径变形法进行监测。在地层中钻取3个不同方向的交汇孔，并分别进行孔径变形测量，即可获得该位置处的三维应力状态。该工程的地应力监测分布及测点的实测地应力值见表2。

表2  地应力实测结果

注:主应力量值单位为MPa。α为主应力在水平面上投影方位角，β为主应力倾角，仰角为正。

第二步，根据工程的空间关系建立位移不连续法（DDM）数值模型。DDM数值模型要求将工程区域的地表曲面划分为若干单元，断层的处理同理，并根据地质勘察报告输入相应的岩石的密度、弹性模量、泊松比，以及断层的粘聚力和内摩擦角。对于本工程案例，根据地质勘察结果，采用几何建模模块勾画出地表轮廓线，进一步得建立三维模型，而断层面一般为平面，直接画出几何边界即可；然后进行模型的网格化，利用网格划分模块的网格划分功能，首先对地表表面进行网格化分，并对每个单元进行编号记录，接着对各断层进行网格划分。最后是，根据参数输入模块的要求，依顺序进行模型参数的输入，先输入地表参数，其中地表处的粘聚力和内摩擦角为0，接着是输入断层的粘聚力和内摩擦角，最后才输入岩土的物理参数。若此时启动模型信息生成模块，可输出用于模拟分析的计算文件，包括模型的几何参数、物理参数等。本工程的位移不连续法（DDM）数值模型如附录中的图1。

第三步，根据第二步生成的计算文件，在DDM数据处理模块输入重力荷载。地表处荷载应力为0，随着竖向深度的增加荷载逐渐增大，启动DDM数据处理模块，经DDM处理，可得到区域应力场，利用结果查看功能提取地应力实测点的应力数据处理值，见表3所示。

表3  重力荷载作用下实测点的应力数据处理值

第四步，根据第二步生成的计算文件，在DDM数据处理模块中输入x方向的单位荷载，该荷载为应力荷载。启动DDM数据处理模块，可得到区域应力场。利用结果查看功能，输入测点坐标，即可得到该处的地应力数据处理值。各测点的应力数据处理值见表4。

表4  x方向单位荷载作用下实测点的应力数据处理值

第五步，施加水平z方向的单位荷载，操作过程同第四步。在DDM数据处理模块中输入z方向的单位荷载，该荷载为应力荷载。启动DDM数据处理模块，可得到区域应力场。利用结果查看功能，输入测点坐标，可得到该处的地应力数据处理值。各测点的应力数据处理值见表5。

表5  z方向单位荷载作用下实测点的应力数据处理值

第六步，依据前面步骤所建立的DDM数值模型进行处理，并反演求和。

依据第三步至第五步得到的各测点的应力数据处理值，通过线性规划法处理各工况对应的回归系数，然后依据回归系数对各工况处理得到的应力场进行组合，求和得到反演地应力场。

将第三步至第五步得到的三组实测点的地应力数据处理值依次输入到回归数据处理模块，然后输入实测地应力值和测点数目，根据回归系数对各工况处理的应力场进行组合，进而反演求和得到地应力场。各组输入数据的格式须统一为：测点编号、x方向正应力、y方向正应力、z方向正应力、xy剪应力、yz剪应力和zx剪应力。

待求解回归系数的方程：

其中，m为测点数目6，为k观测点j应力分量的观测值；为i工况下k观测点j应力分量的数值模型（DDM）处理值；应力分量j＝1,2,…,6，对应6个初始应力分量。

启动回归数据处理模块求解以上方程组，得n个待定回归系数L=(13107,7.7899,19.1199)^T。在反演求和模块输入第二至第五步的地应力场数据处理结果，并输入n个回归系数，回归系数与输入地应力场须一一对应，输入地应力场的格式须为：x方向正应力、y方向正应力、z方向正应力、xy剪应力、yz剪应力和zx剪应力，中间用空格间隔。启动反演求和模块，这样计算域内任一点p的回归初始应力可由该点各工况数值模型（DDM）处理值叠加而得

$>{\sigma }_{\mathrm{jp}}=\underset{i=1}{\overset{n}{\Sigma }}{L}_i{\sigma }_{\mathrm{jp}}^i$>

其中，j对应初始地应力的各个分量，j=1,2,…,6。

各测点的反演处理结果见表6。

表6  各测点的反演应力值

开挖区域的反演地应力场见图2~图4。

本实施例介绍了某地下工程区域地应力场反演。该区域有6个地应力实测点，4条断层，一个坡面，基于6个测点的地应力实测结果，建立DDM数值模型并进行3种荷载条件下的处理，通过线性规划法处理3种情况的回归系数，最后得到了反演应力场，整个反演过程调理清晰，计算简单，整个操作过程具有标准的流程化特点，易于实际应用。

以上对本发明的具体实施例进行了描述。需要理解的是，本发明并不局限于上述特定实施方式，本领域技术人员可以在权利要求的范围内做出各种变形或修改，这并不影响本发明的实质内容。