基于Rhino精准曲面建模的常村煤矿原岩地应力场反演方法

祝凌甫，伊丙鼎，曲秋扬，胡 滨

(1.中煤科工开采研究院有限公司，北京 100013； 2.煤炭科学研究总院开采研究分院，北京 100013；3.煤炭资源高效开采与洁净利用国家重点实验室，北京 100013)

煤矿井下各种地质构造复杂，原岩地应力场分布极不均匀，并且受限于煤矿井下有限的操作环境和高昂的测试费用，很难通过实测的方式获取整个煤矿区的地应力场。而采用数值模拟的方法可以对煤矿区地应力场进行地应力反演，进而获取整个煤矿区的地应力场分布，该方法可以为煤矿井下开采和支护等现场实践提供更加全面和可靠的依据。

国内外学者在地应力反演和分析方面做了大量研究[1-11]。康红普等[12]针对华亭矿区的大型向斜构造附近的地应力场进行了反演；苏生瑞[13]通过离散元数值模拟反演了断层附近的地应力场；张勇慧等[14]采用数值模拟的方法对水电站地下厂房地质模型进行了计算，反演出地下厂房及周边应力场的分布情况；孙猛[15]采用有限元软件ANSYS和有限差分软件FLAC3D反演了平顶山矿区原岩地应力场;张宁博[16]通过数值模拟反演了东荣矿的地应力场，同时分析了断层对应力场和冲击地压发生的影响。

地应力反演结果的精确性取决于很多因素，包括地质模型与现场实际的相似程度、数值模拟过程中各层岩石力学参数赋值的准确性，以及模型计算过程中载荷施加情况等。其中，在三维地质建模方面，通过煤矿区大量的钻孔柱状图对不同岩层进行标识和拟合为目前公认准确度较高的方法，常用的三维建模软件有Rhino、Surfer和ANSYS等。一般煤岩层均为起伏不平的曲面，而部分建模软件在建立曲面方面效果一般，无法准确反映实际地层延伸情况，因此，选择一个能够精准建立曲面模型的建模软件对于地应力反演的准确度至关重要。

Rhino建模软件作为专业的建模软件具有以下几个特点：①Rhino建模方式多样化，效果好，其还带有分析功能和渲染功能，基本满足了3D建模的所有需求；②自动快速形成曲面，曲面精度高，与现场实际贴合性强；③兼容性佳，其默认3D模型保存格式是3dm，也可以多种格式保存，Rhino软件经过转换，可以把三维文件转换成二维图形和线条文件；④可操作性强，界面绘图指令丰富，既可以采用编程的方式绘图，还可以采用直接手绘的方式作图；并且Rhino建模软件建立的三维地质模型与FLAC3D、UDEC及3DEC等数值模拟软件具有非常好的兼容性，可以直接将三维模型导入这些软件进行数值计算。

选取潞安矿区常村煤矿+470 m水平煤岩层为研究对象，采用Rhino软件构建其三维地质模型，采用有限差分数值模拟软件FLAC3D进行地应力场的计算和反演,最后将地应力反演结果与现场实测结果进行对比，评估地应力反演的准确性,以及Rhino建模的适用性。

1 Rhino三维地质模型及数值模型的建立

1.1 工程概况

常村煤矿主采煤层平均厚度为5.98 m。矿井分为2个水平，笔者只针对+470 m水平开展建模和计算。该水平埋深在470～600 m内。煤层直接顶为粉砂岩，基本顶为中粒砂岩，老底为泥岩。井田内地质构造以褶曲为主，小断层相对比较多，大断层主要是aF9断层和FB-1断层，矿井三维地质模型的建立只考虑了这2个断层的存在。

1.2 Rhino三维地质模型的建立思路

由于煤矿井下地层结构的不可见性，无法了解从地表到煤层这部分地层的准确分布情况，因此采用一种简化的办法模拟地表到煤层的地层分布结构。总共分为2个步骤：首先对煤矿钻孔柱状图进行简化，选取不同标志层作为关键地层分界面；然后，在地质钻孔资料中将各标志层的标高进行标识，按照标高将这些标识点连接并拟合成曲面，该曲面与地层的延伸方向相似度比较高，可以表示数值模拟模型中地层的分布。

分析常村煤矿的钻孔柱状图等地质资料，从钻孔资料的精确性和分布的均匀性两个方面考虑，选取了15个地质钻孔，各钻孔标识点分布情况如图1所示。按照上述的简化思路和曲面拟合的原则，将选取的地质钻孔的标志层进行标识，连接拟合出精度最好的曲面来表示岩层的分布。

图1 各钻孔标识点分布图

将拟合出的曲面地层合并，通过Rhinoceros软件建立常村煤矿+470 m水平煤岩层的三维地质模型，如图2所示，其长×宽×高=7 300 m×3 800 m×750 m。常村煤矿煤层赋存最大深度为750 m，即褶曲构造的核心部位，为了更加准确模拟整个常村煤矿的煤层内部地应力分布情况，选择750 m作为所建立的地质模型的高度。

图2 Rhinoceros软件构建的常村煤矿三维地质模型

考虑反演的精度和计算机的能力，将地质模型进行单元划分，划分精度选取为40 m左右，总共生成了281 733个单元，180 026个节点。将划分好网格的模型导入FLAC3D数值模拟软件，生成FLAC3D数值计算模型，如图3所示。本构关系采用Mohr-Coulomb屈服准则。由于断层带附近岩石承载能力基本为零，断层带的模拟采取岩石力学参数弱化的方式实现。

图3 常村煤矿FLAC3D数值计算模型

2 数值模型的计算

2.1 计算原理

煤岩体水平应力一般由两部分组成：一是由自重应力引起的水平应力分量；二是由构造运动引起的水平应力分量[1]。最大水平主应力计算公式如下：

σH=σHg+σHt

(1)

(2)

(3)

式中：σH为最大水平主应力；σHg为由自重应力产生的水平应力分量；σHt为由构造应力场产生的水平应力分量；σV为垂直应力；μ为泊松比；εHt为构造应力引起的水平应变；E为弹性模量。

选取某个地应力测点为初始点进行应力反演(即S3胶带下山测点，其最大、垂直和最小主应力分别为12.87、9.33、6.71 MPa，最大主应力为水平构造应力，为走滑型应力状态)。由式(2)可知，由重力侧压效应所产生的水平应力分量为4.00 MPa；由式(1)计算可得，由构造应力场产生的水平应力分量为8.87 MPa，构造应力引起的水平应变εHt为1.66×10-3。因此，基于以上计算过程可知，为了还原地应力场原始赋存状况，需要施加给模型侧边的水平应变为1.66×10-3。

2.2 计算过程

常村煤矿+470 m水平三维模型总共分为5个边界，在应力场反演过程中需要分别对其施加不同的边界载荷条件或约束条件，这5个边界在计算模型中的数学坐标，以及表示方法如图4 所示。

图4 模型边界数学坐标表示方法

模拟过程分2个阶段：①自重应力场演化。固定模型的侧面和底部位移，施加重力加速度，模型在自重作用下达到平衡。②构造应力场演化。解除模型侧面的法向水平位移固定边界条件，施加速度边界条件，计算直至平衡。应力场演化模拟的边界条件如图5所示。

图5 应力场演化模拟的边界条件

3 地应力反演结果及分析

基于FLAC3D数值模拟计算结果，为了对比同一埋深下的地应力变化情况，采用固定埋深这一变量，输出不同坐标的地应力值，绘制三向主应力折线图，分析煤层内部和基本顶内应力场的分布规律。依据模型数学坐标，设置了4条勘探线来监测矿区内部应力场的变化。具体为：在沿x方向的煤层和基本顶内部，每隔300 m的距离设置程序输出三向主应力量值，绘制煤层和基本顶内部y=600 m、y=2 500 m、y=3 500 m三向主应力沿x轴的变化折线图，分析主应力的变化趋势；在沿y方向的煤层和基本顶内部，每隔100 m的距离设置程序输出三向主应力量值，绘制煤层和基本顶内部x=7 000 m三向主应力沿y轴的变化折线图，分析主应力的变化趋势。各主应力变化情况如图6～9所示。通过该方法，定量分析了常村煤矿+470 m水平井下煤岩应力场的分布及演化。同时，由于本文构建的三维地质模型为曲面模型，更接近地层真实赋存情况，地层起伏等情况导致在远离大型褶曲和断层的地理区域内，地应力反演的数据也会随地层起伏在一定范围内波动。

图6 y=600 m、z=470 m勘探线基本顶主应力随x坐标变化

图7 x=7 000 m、z=470 m勘探线主应力随y坐标变化

图8 y=2 500 m、z=470 m勘探线主应力随x坐标变化

图9 y=3 500 m、z=470 m勘探线主应力随x坐标变化

通过分析可以得出，常村煤矿+470 m水平基本顶最大主应力量值(构造应力)稳定在10～18 MPa内。远离断层区域内，在煤岩层褶曲构造附近，最大主应力达到10 MPa和18 MPa两个极值，其中，向斜构造轴部的应力值最高，达到18 MPa；其他范围内的最大主应力值基本稳定在12 MPa，最小主应力值稳定在6 MPa。因此，在向斜构造区域内进行开采活动时，一定要密切监测顶板应力值变化情况，当应力值产生突变或者远高于其他区域时，应采用加强支护或者卸压的方法维护采场和巷道的围岩稳定。在断层附近区域，地应力产生突变，断层端部会产生应力集中，这也符合断裂力学中裂纹扩展力学机制。但在断层核部则进入完全松弛区，地应力值降低为 6 MPa 左右或者进入完全卸压区，岩石承载力基本为0，地应力值接近于0。因此，在断层区域内进行开采活动时，一定要密切监测应力变化和岩石强度，一般情况下，应对断层附近区域采用加强支护或者预注浆的方法保证工作面围岩的强度和稳定性。

煤炭科学研究总院开采研究分院经过10余年的现场地应力测量工作，在潞安矿区常村煤矿进行了17个不同地点的地应力测量，获取了17组地应力测量数据，实测数据见表1。

表1 常村煤矿地应力实测数据

备注：①由于地应力现场实测过程中未列出地应力测点的地理坐标，无法在地应力反演结果中准确定位，未给出每个地应力测点处的地应力反演结果;② 9号测点的异常性有可能是由不可避免的测量误差、处于地质异常区等因素导致。

经过对比分析实测数据与通过数值模拟获取的地应力反演结果，可以发现，潞安矿区常村煤矿现场实测地应力最大主应力分布在10.17～13.57 MPa内，平均值为12.10 MPa；最小主应力分布在5.30～7.30 MPa内，平均值为6.40 MPa；垂直主应力分布在8.03～11.54 MPa内，平均值为9.20 MPa。数值模拟计算结果显示，在远离地质构造(断层和褶曲等)区域内，最大主应力基本在12 MPa左右波动，最小主应力基本在6 MPa左右波动，中间主应力变化范围相对比较大，但也基本在9 MPa左右上下浮动，与地应力实测的平均值和变化范围一致性非常高。综上分析，两者均显示，常村煤矿+470 m水平最大水平主应力稳定在12 MPa左右，最小水平主应力稳定在6 MPa左右，垂直主应力稳定在9 MPa左右，属于走滑型应力状态。因此，在地应力量值方面，地应力实测数据和地应力反演数据具有比较好的一致性。在三维地质模型建立与现场实际贴合度比较高的条件下，采用数值模拟的方法进行地应力反演具有可行性，能够真实反映煤矿井下原岩地应力场的分布，并且可以弥补在地质构造附近很难获取准确的地应力实测数据的不足。

4 结论

1)利用Rhino建模软件可以比较精确地对煤矿区的地质模型进行构建，以曲面代替岩层更加真实地反映了煤矿井下地层的实际分布情况，褶曲等地质构造可以被真实体现，原岩地应力场的数值模拟反演结果更加准确可靠。

2)大型地质构造区域原岩地应力场与其他区域有较大的差异，向斜轴部原岩地应力场集中程度最大，从断层中部到断层端部再到远离断层区域地应力场呈现“小—大—原”的分布规律。

3)通过对比常村煤矿现场实测的17个地应力测点量值与地应力反演量值，常村煤矿+470 m水平最大主应力稳定在12 MPa左右，最小主应力稳定在6 MPa左右，垂直主应力稳定在9 MPa左右，属于走滑型应力状态，地应力实测结果和地应力反演结果具有高度的一致性。