韩城区块西北部三维地应力场反演分析

赵海峰 ，杨昌增 ，冯 堃 ，甄怀宾

（1.中国石油大学（北京） 石油工程学院，北京 102249；2.中国石油塔里木油田分公司，新疆 库尔勒 841901；3.中石油煤层气有限责任公司，北京 100028；4.中联煤层气国家工程研究中心有限责任公司，北京 100095）

煤层气是一种清洁高效的能源，加快对煤层气的开发与利用能够有效缓解国内的能源供给压力[1-2]。现有煤层气资源勘探调查结果显示，煤层气全球的可采储量约为256.3×1012m3，我国的煤层气资源含量十分丰富[3]。研究表明，韩城区块埋深在1 500 m 以内的煤层气储量达到2 007×108m3，埋 深1 500～2 000 m 之 间 的 煤 层 气 储 量 达 到 了934×108m3。这些数据表明，韩城区块具有较大的煤层气开采价值，因此，它是我国深层煤层气开发的重点对象[4]。

研究显示，现今地应力的分布情况和煤层渗透率的大小决定着煤层气勘探开发的成功与否[5-7]，对煤层气的高效经济开采具有重要的指导意义。在实施煤层顶板压裂的过程中，地应力控制着压裂裂缝从起裂到扩展再到延伸的整个过程，从而决定了压裂的成败[8]。同时，地应力对煤层的渗透率也有着较大的影响[9]。地应力的确定是进行煤层气开采的基础，目前通过钻孔进行地下应力的测量是一种公认的最为可靠的应力测量方式，但由于受地质环境和条件的限制，无法实现大面积的地下应力测量，而在实际的工程应用中需要的是整个区块的地应力分布[10-11]。

因此，根据研究区块的地质资料和实测数据，建立准确地计算模型并进行地应力场反演，对后续的煤层气开发具有重要的指导意义[12-13]。近年来，不同学者提出了多种用于地应力反演分析的方法[14-19]，其中，多元线性回归法原理简单且解唯一，被众多学者所沿用[20-23]。鉴于此，依据韩城北部某区域的地应力实测数据，采用三维有限元回归分析法对研究区进行应力场反演，进一步分析，获得了该研究区目前的地应力场分布规律，为该区域的煤层气的开采和增产提供了重要的数据基础。

1 研究区域概况

韩城区块地质构造为北东走向、北西倾向的单斜构造。位于鄂尔多斯盆地渭北隆起东北部，渭北隆起自三叠纪后经历了印支运动期NNW 向挤压、燕山运动中期NWW 向挤压和喜马拉雅运动期NE 向挤压与NW-SE 向伸展作用三期构造作用。印支、燕山、喜马拉雅三期构造运动叠加形成该区域构造形态，其中燕山运动奠定了本区的构造基础。

因此，韩城区块挤压构造带主要发育在东南边浅部地层，虽然伸展构造遍及全区，但主要集中在文家岭隆起以南地区，而且发育规模及强度从南向北、从边浅部向中深部慢慢变弱。南北区域构造发育的差异性表现在南区以断层为主，褶皱发育较少，北区则相反。在区域应力挤压的作用下，断层较为发育，受3 条断裂带的影响，韩城区块煤层气的富集差异大。研究区块构造纲要如图1。

图1 研究区块构造纲要图Fig.1 Outline diagram of research block structure

2 地应力场实测数据分析

实测地应力数据是现今地应力场分析极其重要的数据基础，在地应力反演过程中实测数据亦是必不可少的重要参数。因此，在进行地应力反演前需进行地应力实测数据分析。通过钻孔进行地下应力的测量是目前一种公认的最为可靠的应力测量方式，但是现场应力实测需要耗费巨额的资金，并且受地质环境和技术条件的限制，无法实现大面积地下应力测量。因此，利用水压致裂法获得的6 个测点的实测应力数据对研究区块的应力场进行反演分析。

2.1 水力压裂法原理

向地层封隔段中注入高压流体，记录地层破裂的整个过程，通过压裂曲线得出垂直于钻孔平面的最小水平主应力的大小，然后通过经验公式算出最大水平主应力。弹性力学是水压致裂法测量原理的基础，须满足以下3 个要求：①岩石是线性、均匀、各向同性的弹性体；②岩石是完整的且不具有渗透性，注入的流体符合达西定律；③钻孔轴方向与垂向主应力方向相同。

可以用3 个应力分量来表示地壳中某一点的应力状态，而钻孔中某处的应力状态可以通过压裂曲线计算得出。在压裂的过程中，封隔段内的液柱压力会不断增大，壁面上的切向应力将随之变小，最后慢慢演变为张应力。当压力超过岩石的抗拉强度 σt时，岩石会发生破裂。此时的液体压力被称之为破裂压力，其表达式为[24]：

式中：pb为 破裂压力；Sh为钻孔破裂面上最小主应力，MPa；SH为钻孔破裂面上最大主应力，MPa； σt为 岩石的抗拉强度；p0为孔隙压力。

计算最大、最小主应力的公式如下[24]：

式中：ps为 瞬时关闭压力；pr为裂缝重压力。

垂向应力Sv可以通过地层的密度计算得出。

2.2 测点布置

研究区域的地应力测量共选择6 个勘探钻孔作为应力测试孔。所有测点均布置在煤层顶板砂岩层中。除6 号钻孔外，其他钻孔深度均未超过1 000 m，最大测点深度为1 003 m。在水压致裂测量试验结束之后进行印模定向测试，所测得的破裂面的走向就是最大水平主应力方向。测点布置如图2，钻孔水压致裂法地应力测量结果见表1。

表1 钻孔水压致裂法地应力测量结果Table 1 Measurement results of in-situ stress by drilling hydraulic fracturing method

图2 测点布置图Fig.2 Measuring points layout

以正东向为x轴正向，正北向为y轴的正向，竖直向上为z轴的正向，建立三维直角坐标系。考虑到进行地应力场反演分析的过程中都是以应力分量的形式来表征地应力进行回归计算，而水压致裂法测得的地应力值是以主应力的平面方位和倾角给出，所以为了便于计算，将根据方向余弦将实测地应力值进行坐标变换。

首先，计算实测主应力与建立的坐标轴之间的方向余弦[25]经进一步推导得出：

式中：Li、Mi、Ni分别为 σi对x、y、z轴的方向余弦； σi为 实测主应力； βi为 σi与水平面之间的夹角； αix为 σi与x轴正方向间的夹角。

根据所测得的主应力的大小和方向余弦，再计算主应力转换坐标系下的应力分量[26]：

式中：σ1、σ2、σ3分 别为原始坐标系下x、y、z方向的主应力； σx、σy、σz分别为x、y、z方向的主应力； τxy、τxz分别为与y向力相关的剪应力和与z向力相关的剪应力； τyz为σy与z向力相关的剪应力。

坐标系下实测地应力分量值见表2。

表2 坐标系下实测地应力分量值Table 2 Measured geostress component values in the coordinate system MPa

综合上述6 个测点的实测地应力资料可以得出，地应力状态主要表现为Sv>SH>Sh，最小水平主应力小于16 MPa，三向主应力的这种分布关系表明现今地应力处于拉伸应力状态，工程类比判断韩城区块的最大主应力方向为N45°E，2 个钻孔水压致裂法实测均为N45°E 上下浮动3°～5°，最大主应力方向基本一致。根据现场测井解释结果显示韩城区块的平均地应力梯度为2.56 MPa/hm；6 个测点埋深范围为754～1 002 m 之间，测得的垂向应力基本符合测井解释的结果。

3 地应力场三维有限元回归分析

3.1 地应力场回归分析原理

为了全面掌握研究区域的现今地应力场分布规律，引用有限元应力回归分析法进行计算。该方法基于数值分析理论，通过计算机技术手段，把现有的地质资料和现场实测地应力资料及数理统计理论结合起来，充分考虑和分析影响地应力的因素，保证反演出的应力场最大限度地接近实际地层应力。其基本原理是：

1）根据研究区域的地质资料，建立起精确的有限元数值模型。

2）根据现今地应力成因理论，选定垂直方向的自重因素和水平方向的2 个构造因素作为待定因素，以单位载荷分别作用于先前建好的数值模型，得到各因素在单位载荷下的有限元计算值。

3）根据应力线性叠加原理，将地应力回归值σ0n作为因变量，把有限元计算求得的单因素的地应力计算值σin作为自变量，则回归方程的形式为[27]：

最小二乘法残差的平方和Q为[27]：

应力分量m取1～6，分别对应于6 个分量，根据最小二乘法的原理，求解使残差的平方和取得最小值的方程式[27]，即可得到k个回归系数D1、D2、···、Dk，将其代入式（13）即可求得模型区域任意位置的地应力回归计算值。

3.2 三维有限元计算模型的建立

为了提高有限元力学模型的精度，充分利用现有测井资料，在Petrel 中建立地质模型，并通过测井曲线将泊松比和弹性模量粗化到整个网格中，最终在Fracman 中完成地应力反演过程。建模思路如图3。

图3 建模思路图Fig.3 Modeling idea diagram

首先，从研究区域的井轨迹和分层数据中得出各个层段（顶板顶部、煤层顶部、煤层底部、底板底部）的各层面点数据；然后通过Petrel 的make surface 模块进行曲面拟合，形成三维空间的曲面，通过三维建模模块生成三维网格模型；最后，将测井曲线中的泊松比和弹性模量粗化至每个网格中。将得到的模型导入Fracman 中，通过裂缝建模生成煤层的割理裂隙，并将割理裂隙通过Oda地应力计算粗化至每个网格中以弹性模量和泊松比的变化体现出来[23]。该模型计算域为：东西向分别以19 433 040 m 和19 437 060 m 坐标网格为边界，南北向分别以3 933 605 m 和3 936 485 m坐标网格为边界，东西向宽4 020 m，南北向长2 880 m，顶板下表面深度取到海拔-245 m，上表面深度取-28 m，计算域内共划分178 640 个六面体单元，三维有限元计算模型如图4。

图4 三维有限元计算模型Fig.4 3D finite element calculation model

如果直接在模型加载应力边界条件，会产生边界应力集中，与实际应力情况有很大差别，且煤层和顶底板砂岩层起伏较大，直接加载应力也会产生较大误差。

由于自重作用和构造运动是引起地应力的主要因素，因此，模拟的边界条件只考虑这2 个因素。构造应力场力学模型为：2 个侧面分别施加x方向的I1的单位均布载荷和y方向的J1单位应力载荷，来模拟水平方向构造运动作用；在垂直方向加载K1的单位应力载荷来模拟上覆地层压力；在相应的未加载应力的对侧分别施加x、y、z方向的位移约束。

3.3 初始地应力场反演回归结果

基于6 个测点的测试数据和3 个地应力场计算结果，最终得到对应于自重应力场 σz、水平方向构造应力场 σx和 σy的3 个回归系数，分别为D1=20.2，D2=18.5，D3=12.5。因此，得到本次研究区的地应力场σ0回归方程为：

利用式（15）即可得到测点位置处的应力回归计算值，通过数学关系转换得到对应的最大和最小水平主应力，并与应力实测数据进行对比分析。实测数据与回归计算值对比见表3。

表3 实测数据与回归计算值对比Table 3 Comparison of measured data and regression calculated values

由表3 可知，最大、最小水平主应力和垂向应力都表现出随着埋深的增加而逐渐增大的趋势，且三者的大小关系基本保持为Sv>SH>Sh。现今地应力处于拉伸状态，本次反演结果与现场实测结果的应力特征基本保持一致。实测值与回归计算值之间的误差范围在0～5 MPa 之间。其中，垂向应力的误差较小基本在0～1 MPa 之间浮动。这是因为垂向应力与地层的深度基本呈线性变化，拟合结果较为准确；水平应力的误差偏大，在2～5 MPa 之间波动。由于水平方向的应力成因复杂，单纯用一种线性载荷来模拟是不够严谨的，因此误差偏大，但也在工程应用可接受最大不超过25%的范围之内。因此，本次反演结果能够满足工程精度的要求，反演得到的地应力场可表达研究区地应力场的主要特征。

煤层顶面水平最大主应力大小云图如图5，煤层顶面水平应力差云图如图6，煤层顶面水平最大主应力方向如图7。

图5 煤层顶面水平最大主应力大小云图Fig.5 Cloud diagram of horizontal maximum principal stress on top surface of coal seam

图6 煤层顶面水平应力差云图Fig.6 Cloud diagram of horizontal stress difference on top surface of coal seam

图7 煤层顶面水平最大主应力方向Fig.7 Direction of maximum horizontal principal stress on top surface of coal seam

从图5 可以看出，煤层的水平最大主应力整体趋势还是随深度的增加而增大，变化范围在10.7～25 MPa 之间，基本中部、北西向高，东部低的趋势；在研究区中部有高应力区的存在，这是可能是由于此处处于地层的褶皱区，在顶部受到了较大的拉伸应力造成的。

从图6 可以看出，煤层的整体应力差处于1～5 MPa 之间，在北西、西南和中南部存在较大的应力差，结合前面的应力分析，其原因可能是此处的地层起伏较大，构造应力较为复杂，导致水平最大主应力的大小超过了垂向应力，出现了高应力区。所以，建议在前期布井时应该尽量避免这3 个高应力差的地区，尽量选择在东北部布井。

从图7 可以看出，除了中部的高应力区外，整体的应力方向趋势呈N45°E。在这2 个高应力区的应力方向由水平向突变为竖直向，初步分析认为是此处的最大水平应力超过了垂向应力，而绘制这个应力方向矢量图用的是水平方向和倾角，此时的中间应力由原本的最大水平应力方向变为了垂向应力，所以会发生应力突变。

4 结 语

1）通过分析研究区域的构造特点、实测数据和已有的研究成果，经多次的反演计算之后得出，在垂向施加20 MPa、东西向施加18.5 MPa、南北向施加12.5 MPa 应力，所取得的反演效果最好。

2）测点位置处的三向应力大小均表现出随着埋深的增加而逐渐增大的趋势，而现场实测应力由部分与回归计算值有较大差异。这由于是测点的位置在煤层与顶板的交界处，该处的地质参数不稳定从而影响地应力大小产生不规律的变化。3）煤层的水平最大主应力处于10～25 MPa 之间，呈现中部、北西向高，东部低的趋势；煤层的整体应力差处于1～5 MPa 之间，在西、西南和中南部有较大应力差存在；究区整体应力方向趋势呈N45°E。建议在前期布井阶段应尽量避免这3 个高应力差的地区，尽量选择在东北部布井。