基于三维构造恢复的地应力预测技术

The Stress Prediction Technique Based on 3D Structural Restoration

摘要: 组合弹簧理论模型是预测地应力的一种常用方法，但模型忽略了对岩层各向异性的考虑，无法反映由于局部断层、裂缝发育等各向异性造成的应力局部扰动变化；对构造应力项的求取基于薄板理论，简单利用地层形变开展地应力预测，忽略了地质力学参数的空间变化。本文通过叠前弹性参数反演和精细构造解释对常规弹簧组合模型的构造应力项求取进行了优化，提出了基于三维构造恢复的地应力预测技术，从构造建模开展地质力学分析与构造恢复，求取现今地应力的古构造应力残余分量，提高了模型对地应力的预测精度，结果与实钻井具有较好的吻合性，并通过实际工区应用验证了该方法的有效性，对页岩气勘探选区评价、井位部署、钻后评价、水平井设计具有借鉴作用。

Abstract: The composite spring model is a commonly used method for predicting the in situ stress. However, the model ignores the consideration of the anisotropy of the rock mass, which cannot reflect the stress disturbance caused by the anisotropy due to local faults, fractures, and other factors. Moreover, the calculation of the tectonic stress term is based on the thin plate theory, which simplifies the geological deformation and neglects the spatial variation of the geological mechanical parameters. In this paper, the construction stress term of the conventional spring combination model was optimized by using pre-stack elastic parameter inversion and refined tectonic interpretation. A 3D tectonic reconstruction-based in situ stress prediction technology was proposed, which analyzed the geological mechanics and tectonic reconstruction from the perspective of structural modeling. The residual tectonic stress component of the present-day stress was obtained, improving the prediction accuracy of the model for in-situ stress. The results have good consistency with actual drilling, and the effectiveness of the proposed method was verified by application in a practical working area. This method has reference significance for shale gas exploration zone evaluation, well deployment, post-drilling evaluation, and horizontal well design.

1. 引言

页岩储层较常规油气储层具有“低孔低渗”的特点，压裂改造是实现页岩储层商业开发的必须条件。随着勘探的推进与深入，页岩气勘探从浅层走向到深层 [1]。随着埋深增加，地层温度、压力增大，页岩塑性增强，水平井压裂改造难度增大。研究表明地应力总体表现为随着埋深的增大而增大；同时受到构造挤压作用，盆缘复杂构造带靠近控盆断裂地应力梯度较高，盆内地应力梯度较低；宽缓构造区具有较低的地应力梯度，低水平应力差，是深层页岩气的有利甜点区。地应力变化特征及演化规律是深层页岩气可压性评价的关键参数，相对低的两向水平应力差有利于压裂形成复杂缝网 [2] [3] [4] [5]。由于地下介质受到挤压、扩张等构造运动及岩层的结构、温度、压力、孔隙度等各不相同，导致不同岩层的自重应力和构造应力显著不同。上覆岩石的重量和岩石孔隙、骨架等决定了地应力的垂向分量大小，而构造运动导致岩层发生位移，决定了地应力水平应力的大小，并共同构成了复杂的地应力场 [6]。

目前国内外针对地应力的研究和预测主要从工程、测井、数值模拟、地震预测等几方面展开，形成了基于工程的水力压裂法、应力恢复法、应力解除法等，基于声波测井、成像测井的应力计算方法，基于数值模拟的有限元法、边界位移调整法、应力函数法、位移反分析法，基于地震预测的反射系数法、曲率法、岩石物理模型法、弹簧组合理论法等 [7] [8] [9]。其中，基于工程和测井资料的地应力预测方法成本高、代价大，能较准确地反映出地应力的局部大小和方向；基于数值模拟的方法未能考虑地层的地质力学参数和各向异性特征，对地应力的预测精度不高；基于地震预测的方法缺乏考虑构造应力的影响，适用性较差。针对以上几类地应力预测方法存在的不足，本文基于弹簧组合理论，通过叠前弹性参数反演和精细构造解释对常规弹簧组合模型的构造应力项求取进行了优化，提出了基于构造恢复 [10] [11] [12] 的地应力预测技术，以提高地应力的预测精度，促进对页岩气选区评价、井位部署和水平井方位的优选。

2. 方法原理

基于弹簧组合理论，现今地应力的组成由上覆岩石压力、地层孔隙压力分量、各向异性应力扰动量、构造恢复古应力残余量构成。其中上覆岩石压力主要引起垂向应力的变化，地层孔隙压力分量、各向异性和构造应力的变化主要引起水平应力的变化。因此，在实际应力场描述中通常采用垂向应力、最大水平地应力、最小水平地应力进行地应力的表征，同时用水平应力差作为可压性评价的关键参数。相应应力项的求取公式主要如下 [13] ：

① 垂直应力

垂直方向上的应力为上覆岩层的重力，随埋深增加而增加：

$S_v={\displaystyle \int \rho \left(h\right)g\text{d}h}$ (1)

式中，S_v为垂直应力；h为地层埋深；ρ为地层的密度，随埋深变化而变化；g为重力加速度。

② 水平应力

水平方向上的应力较为复杂，组成因素较多。由于地下储层受到上覆岩石重力、压实作用及不同岩石骨架产生的孔隙压力、构造运动导致的构造应力等诸多因素的共同影响，造成地下储层各点的应力状态各不相同。其中，垂直应力对岩层进行挤压时，不仅会产生垂直应力，还会使岩石横向上产生形变，产生水平应力；岩层各向异性的存在导致不同的孔隙应力等其他应力的叠加；构造运动造成岩层发生位移，产生的构造应力是水平方向的主要部分。

假设水平方向的应变受限，可得到由垂向应力诱导的最大水平应力S_H1与最小水平主应力S_h1相等，可以表示为：

$S_{H1}=S_{h1}=\frac{v}{1-v}\left(S_v-\alpha P_p\right)+\alpha P_p$ (2)

式中，v表示地层的泊松比；P_p为孔隙压力，α为Biot系数；S_v为垂直应力。

构造运动产生的水平应力是由于地质板块和地壳构造运动等活动而产生的作用力，由于构造作用的部位不同且自然界影响因素较多，环境复杂多变，因此不同地区不同深度的构造应力大小和方向各异，根据胡克定律其表达式为：

$\{\begin{array}{l}{S}_{H2}=\frac{E{\epsilon }_H}{1-v^2}+v\frac{E{\epsilon }_h}{1-v^2}\\ S_{h2}=\frac{E{\epsilon }_h}{1-v^2}+\nu \frac{E{\epsilon }_H}{1-v^2}\end{array}$ (3)

式中，E地层的杨氏模量；ε_H为地层的最大主应变；ε_h为最小主应变。

水平主应力除了由上述垂直应力引起的水平应力、构造应力、孔隙应力构成之外，还受到其他应力的影响，如塑性泥岩和石膏等因素的影响造成地应力的软化现象、岩石中矿物的变化造成的局部应力的变化，由于其计算难度较大，在研究过程中尚未考虑。因此，水平方向上的应力由垂直应力引起的水平应力以及构造应力两者叠加而构成：

$\{\begin{array}{l}{S}_H=S_{H1}+S_{H2}=\frac{v}{1-v}\left(S_v-\alpha P_p\right)+\alpha P_p+\frac{E}{\left(1-v^2\right)}\left({\epsilon }_H+v{\epsilon }_h\right)\\ S_h=S_{h1}+S_{h2}=\frac{v}{1-v}\left(S_v-\alpha P_p\right)+\alpha P_p+\frac{E}{\left(1-v^2\right)}\left({\epsilon }_h+v{\epsilon }_H\right)\end{array}$ (4)

式中，S_H为最大水平主应力；S_h为最小水平主应力；v表示地层的泊松比；P_p为孔隙压力，α为Biot系数；E为地层的杨氏模量；ε_H为地层的最大主应变；ε_h为最小主应变。

常规弹簧组合模型中基于曲率的原理将求取的最小、最大曲率等效为地层受构造运动所产生的最小、最大主应变量，应变与曲率的关系可由下式表达：

$\{\begin{array}{l}{\epsilon }_H=\frac{h}{2}{K}_{\max }\\ S_h=\frac{h}{2}{K}_{\min }\end{array}$ (5)

式中，ε_H为最大主应变、ε_h为最小主应变。将式(1)和式(4)相加，带入式(7)即可得到常规弹簧组合模型的求取公式。

弹簧组合理论求取地应力的方法基于薄板理论，简单利用地层形变开展地应力预测，忽略了地质力学参数的空间变化。改进后的地应力预测方法优化了构造应力项最大、最小应变的求取方法。新方法结合地质认识和区域构造运动背景，在构造精细解释的基础上，根据研究区特点确定模型边界，建立区内边界、断层、层位之间的约束连接，结合层序地层学沉积初始状态为水平面的观点开展三维地质建模及恢复。应变的建模和计算采用动态松弛和有限元方法，将断层和层位进行网格化处理，并遵循最小应变原则，得到最小应变下模型拉平的最优解。模型内部的岩性和力学特征参数由叠前弹性参数反演产生，使得建立的模型体更加接近真实的目标地质体。通过提取拉平前后的形变量，可得到格林拉格朗日应变结果，通过对分量计算得到最大水平应变量和最小水平应变量。具体的求取公式如下：

在岩石力学中，应变被定义为岩石受力产生的形变，具体的公式可由下式表达：

$\overline{\stackrel{¯}{F}}=\frac{\partial \stackrel{\to }{x}\left(\stackrel{\to }{X},t\right)}{\partial \stackrel{\to }{X}}$ (8)

式中，X是受力前的初始位置，x是受力后的位置。二者的梯度比值表示应变量，但是用变形梯度进行应变量的表征对简单的平移或者旋转运动得到的应变量不为零，具有一定的局限性。

本文的做法是将应变视作空间中任意点上形变张量的对称部分，即格林拉格朗日应变，其矩阵形式可由下式表达：

$\bar{E}=\frac{1}{2}\left({\bar{F}}^T\cdot \bar{F}-\bar{I}\right)$ (9)

格林拉格朗日应变由9个分量组成，各分量之间可按照不同的计算公式得到最大水平主应变和最小水平主应变，具体计算公式如下：

$\{\begin{array}{l}{\sigma }_1=i{s}_0+j{s}_1+k{s}_2\\ {\sigma }_2=i{s}_3+j{s}_4+k{s}_5\\ {\sigma }_3=i\left(s_1{s}_5-s_2{s}_4\right)+j\left(s_2{s}_3-s_0{s}_5\right)+k\left(s_0{s}_4-s_1{s}_3\right)\end{array}$ (10)

式中，σ₁代表最大水平主应变，s代表格林拉格朗日分量，σ₃代表最小水平主应变。

根据组合弹簧理论地应力模型将构造应力项中最小、最大主应变求取方法用三维构造恢复的方法进行替代(式10)，结合岩石各向异性特征，得到本文的地应力预测新方程如下：

$\{\begin{array}{l}{S}_H=\frac{v\left(1+v\right)}{1+E{Z}_N-v^2}\left(S_v-\alpha P_p\right)+\alpha P_p+\frac{E}{1-v^2}\left(v{\sigma }_3+{\sigma }_1\right)\hfill \\ S_h=\frac{v\left(1+v\right)}{1+E{Z}_N-v^2}\left(S_v-\alpha P_p\right)+\alpha P_p+\frac{E}{1-v^2}\left({\sigma }_3+v{\sigma }_1\right)\hfill \end{array}$ (11)

式中，S_h、S_H为最小、最大水平主应力；ν为泊松比、E为杨氏模量，Z_N为法向柔度，S_v为垂向应力，α为Biot系数，P_p为为流体压力，其中垂向地应力S_v可通过密度、重力加速度在深度域积分实现，P_p为流体压力可通过系列成熟的压力预测方法(Eaton法、Fillippone法、Stone法等)实现，ν、E可通过叠前弹性参数反演获得；σ₃、σ₁为最小、最大主应变，可通过构造恢复进行求取。

3. 实际应用

根据上述新推导改进后的弹簧组合模型地应力预测新方程，以四川盆地某工区五峰–龙马溪组优质页岩层段为例，开展三维构造恢复和地应力预测，研究工区内不同实钻井在相同地质条件下压裂测试产量差异大的具体原因，分析地应力分布规律，验证方法精度(具体方法流程如图1所示)。该研究区内实钻井(Well_1井、Well_2井、Well_3井、Well_4井)在优质页岩厚度、TOC、孔隙度、含气量等地质条件相当的情况下压裂测试产量具有较大的差异，其中Well_1井和Well_2井具有高产的结果，而Well_3井和Well_4井的产量相对较差(实钻井统计表见表1)。

Figure 1. Technical flow of 3D structural restoration method

图1. 三维构造恢复方法技术流程图

Table 1. Comparison table of geological parameters of wells in the study area

表1. 研究区实钻井地质条件参数对比表

结合研究区构造特征，现今的构造面貌是在两期构造应力的叠合作用下形成的。因此，基于研究区OVT域叠前时间偏移地震数据，考虑研究区内受构造运动主要影响造成的大断层，针对五峰–龙马溪组目的层解释层位TP2l、TP1l、TO3w共三个层系，主要断层12条，参考已钻井4口(well_1井、well_2井、well_3井、well_4井)，完成深度域精细构造解释，为建立三维模型提供可靠数据(如图2所示)。

Figure 2. Fine structure interpretation results in the study area

图2. 研究区精细构造解释成果图

三维建模基于paradigm-Gocard软件平台，结合研究区的断层发育情况和交接关系，对原始断层数据进行编辑，确定断层间的主要断层和分支断层，使得断层之间的交接关系符合实际地质认识。然后，基于有限元网格原理，设置横向和纵向网格大小，将断面进行网格化处理，得到三维网格化插值后的断层模型(如图3所示)。

Figure 3. 3D structural modeling model of the study area

图3. 研究区三维构造建模模型图

基于断层网格模型，结合构造精细解释，建立地质网格模型(如图4所示)。通过建立的三维地质网格模型，结合研究区弹性参数反演得到的实钻井处的密度、泊松比、杨氏模量等参数(如图5所示)进行均值求取，得到目的层处的参数平均值。然后，将参数平均值带入到地质网格模型中，设置三维构造恢复时目的层处的地球物理参数。

反演结果表明研究区杨氏模量、泊松比随着埋深增加而增加，其中Well_1、Well_2井杨氏模量、泊松比较低，Well_3、Well_4井杨氏模量、泊松比较高。结合建立的地质网格模型及对古构造所受地应力方向的地质认识，设置三维构造恢复时的钉面。钉面在进行构造恢复时固定不动，其他面可自由滑动，得到模型三维构造恢复拉平前后的对比图(如图6所示)。

Figure 4. Geological grid model of the study area

图4. 研究区地质网格模型图

(a) (b)

Figure 5. Pre-stack elastic parameter inversion planar graph. (a) Inversion plan of Young’s modulus in the study area; (b) inversion plan of Poisson’s ratio in the study area

图5. 叠前弹性参数反演平面图。(a) 研究区杨氏模量反演平面图；(b) 研究区泊松比反演平面图

Figure 6. Comparison of 3D structure restoration and leveling of the target layer in the study area

图6. 研究区目的层三维构造恢复拉平前后对比图

通过提取构造恢复中地层的构造形变程度(即地层膨胀度)，可得到研究区受构造运动影响形成的地层形变大小(如图7所示)。其中，膨胀度大于1表示压缩，等于1表示古今构造无变化，小于1表示拉伸。从图中可以看出，well_4井受到的压缩最大，well_1井、well_2井和well_3井受到的压缩程度相当。

Figure 7. Expansion degree plan of the Wufeng-Longmaxi Formation of the target layer in the study area

图7. 研究区目的层五峰–龙马溪组膨胀度平面图

结合前述中(9)的公式，将构造恢复拉平前的三维模型体积与拉平后的三维模型体积求梯度，得到反映构造恢复形变大小的格林拉格朗日应变量矢量(如图8所示)。图中颜色代表应变的大小，从图中可以看出在well_4井处的应变程度最大，而在well_1井、well_2井、well_3井相当。

Figure 8. Strain vector plan of Wufeng-Longmaxi Formation in studyarea

图8. 研究区五峰–龙马溪组应变矢量平面图

进一步的将求取得到的格林拉格朗日应变的9个分量，按照(10)中的公式进行计算，可得到研究区受构造运动影响形成的地层应变大小，即构造运动产生的得到最大水平主应变和最小水平主应变(如图9所示)。

(a) (b)

Figure 9. Maximum and minimum horizontal principal strain plan of the Wufeng-Longmaxi Formation in the study area. (a) Maximum horizontal principal strain; (b) Minimum horizontal principal strain

图9. 研究区五峰–龙马溪组最大、最小水平主应变平面图。(a) 最大水平主应变；(b) 最小水平主应变

将构造恢复求取的最大最小主应变带入到新的地应力计算方程中，得到新方程的地应力计算结果(如图10所示)。

(a) (b)

Figure 10. Maximum and minimum horizontal principal strain plan of the Wufeng-Longmaxi Formation in the study area. (a) Minimum horizontal principal stress prediction plane; (b) Maximum horizontal principal stress prediction plane

图10. 研究区五峰–龙马溪组最大、最小水平主应变平面图。(a) 研究区最小水平主应力预测平面图图；(b) 研究区最大水平主应力预测平面图

通过分析研究区最大最小水平主应力预测结果可知，研究区地应力规律整体表现为随着埋深增加而增加的趋势，在构造紧闭部位(图中well_3井、well_4井附近)地应力高，在背斜翼部地层平缓，地应力较小(well_1井、well_2井附近)。为了进一步说明该方法的有效性，将最大、最小水平地应力作差与基于常规弹簧组合模型应用曲率原理求取得到的水平应力差做对比(如图11所示)。

(a) (b)

Figure 11. Comparison result of horizontal stress difference between old and new methods in the study area. (a) Horizontal stress difference prediction of traditional method; (b) Horizontal stress difference prediction of new method)

图11. 研究区新老方法水平应力差对比结果图。(a) 传统方法水平应力差预测平面图；(b) 改进后弹簧组合模型水平应力差预测平面图

通过对比图11中常规和改进后的弹簧组合模型求取地应力方法得到的水平应力差预测平面图可知常规弹簧组合理论求取的水平应力差在well_3井、well_4井处较小，而新方法水平应力差较高，解释了well_3井和well_4井压裂测试产量低的原因。

将新方法计算得到的水平应力差在各井点进行统计，预测well_1井和well_2井水平应力差15 MPa左右，well_3井和well_4井预测应力差大于17 MPa。结合实钻井压裂数据，well_1井和well_2井压裂裂缝以复杂缝为主，返排率低，测试产量高；well_3井和well_4井压裂以主缝为主，返排率高于50%，测试产量低。新方法对地应力的预测结果较老方法更加符合实际压裂情况，说明新方法具有更高的精度和适用性(如表2所示)。

Table 2. Statistical comparison table of horizontal stress difference to fracturing in the study area

表2. 研究区水平应力差于压裂情况统计对比表

4. 结论

本文基于弹簧组合模型，推导了基于三维构造恢复的地应力预测技术，将该方法应用到四川盆地某研究区得到了较好的应用效果，得到了以下几点认识和结论：

1) 现今地应力的组成复杂，主要由上覆岩石压力、地层孔隙压力分量、各向异性应力扰动量、构造恢复古应力残余量构成；

2) 常规组合弹簧理论模型求取地应力忽略了对岩层各向异性的考虑，无法反映由于局部断层、裂缝发育等各向异性造成应力释放的分布规律；

3) 新方法从构造建模、岩石力学参数反演开展地质力学分析与构造恢复，求取现今地应力的古构造应力残余分量，对弹簧组合模型进行了优化，解决了地层应变求取难题，提高了地应力预测精度，结果与实钻井具有较好的吻合性，为页岩气水平井部署提供了借鉴。