引言

地应力是地壳内部应力的总和，一般认为，地应力主要由上覆岩层重力、地层压力、构造应力等方面构成^[1-2]。地应力影响油气运聚^[3-4]、断层封闭性^[5-6]、构造裂缝与压裂裂缝预测^[7-8]、井壁稳定性^[9-10]及井网部署^[11-13]，是石油天然气领域不可或缺的数据^{[10, 14-16]}。

地应力的概念从20世纪初被提出以来，人们对地应力的认识与研究经历了由浅到深的过程。地应力的测量方法、技术以及研究手段不断地改进与发展。目前，基于岩芯或井孔的地应力测量是明确地应力状态最直接、有效且成熟的手段，基于测井资料的地应力计算能够获取全井段连续的地应力数据，是地应力测量的重要补充^[17]。近年来，数值计算法在三维地应力的研究中应用愈加广泛，以井点地应力值为约束，采用一定的算法，包括遗传算法、神经网络法及模糊算法等^[18-19]，开展区域应力场预测。然而，受岩芯样品和测试成本的局限，地应力测量无法大规模开展，测井资料所计算的单井地应力状态也难以反映盆地或油藏整体应力场分布特征。虽然数值计算法在理论上具有便捷、经济且准确的三维预测结果，但往往受限于样本数量的不足，对地质条件复杂的地区不具备较强适用性^[17]。

目前，基于单井约束的应力场有限元数值模拟法是一种应用广泛、适用性强、精度较高的预测方法，在复杂断块区、非均质较强的深层浊积砂体以及物性特征多变的煤田区块都取得了成功应用^{[17, 20-21]}，经历数十年来的发展，学者们进行了多层面的改进与创新，实现了从简单二维均质力学模型到复杂三维非均质力学模型的转变。然而，在进行应力场预测的工作中，为了得到一种最优的加载方式，来达到预测的应力场井点处大小和方向与实测吻合，需要人为地反复尝试确定边界条件。这一过程往往花费较大的精力与时间，特别是对于大型复杂模型，影响研究效率。因此，针对这一问题，笔者基于ANSYS中的APDL语言编译自适应边界程序，将人为调整的过程自动化，实现应力场的智能预测，提高研究效率，并以渤南油田义176区块沙四上亚段的3、4砂组为例，预测其现今地应力场的三维空间分布。

1 研究区概况

渤南油田位于济阳拗陷沾化凹陷中部的渤南洼陷内。渤南洼陷是一个北陡南缓向北掀斜的箕状断陷盆地，受几条近东西向盆倾断层切割，从南到北依次形成南部断鼻带、南部断阶带、中部洼陷带和北部陡坡断裂带^[22]。义176块位于渤南油田义176—渤深4断阶带，面积49.8 km$^2$，受两条边界断层控制，区内断层走向基本为近东西向。

研究目的层为沙四上亚段的3、4砂组，各划分为9个小层，埋深3 200$\sim$4 500 m，发育扇三角洲前缘亚相和滨浅湖亚相，岩性以膏岩、碳酸盐岩、泥岩为主，夹杂灰褐色、褐色的含砾砂岩、细砂岩级粉砂岩，平均孔隙度约11.17%，平均渗透率为3.15 mD^[22]，属于低孔隙、低渗透的薄互层砂砾岩油藏，在平面和垂向上呈现出较强的非均质性(表 1)。

2 井点现今地应力特征 2.1 现今地应力方向

基于“横波分裂原理”，根据多极子阵列声波测井资料中的快横波方向能够确定现今地应力方向。由于横波信号在各向异性的地层中传播时会分裂成速度较快的快横波和速度较慢的慢横波，而快横波的传播方向主要沿岩石中天然裂缝的延伸方向或最大主应力方向，通常与地层各向异性相对应，故快横波的传播方向一般就是地层中现今水平最大主应力方向^[23-24]。因此，可以通过测量快横波的方位角来确定现今水平最大主应力方向。

通过对义176区块9口具有多极子阵列声波测井数据的井进行测井解释，确定了这9口井在目的层的现今地应力方向(图 1)。可以看出，现今地应力方向主要呈NWW—近EW向，受边界断层控制明显，最大水平主应力方向与断层走向基本一致。

2.2 地应力数值

声发射与差应变法是目前相对成熟且较为简单易行的地应力测试方法。

声发射法的原理是脆性材料对曾经的受载作用具有记忆性，通过观察岩样在加载过程中发出的声信号变化，判断Kaiser效应点，进而获得地应力值^[25]。为了克服岩石非均质性造成的测试结果误差，从岩芯样品4个方位钻取柱塞样(图 2)，分别计算各方位的应力分量，再用式(1)求取主应力的大小。共对义176区块6口井的6组(24块)样品开展了声发射测试，测试结果见表 2。

差应变法测试现今地应力的原理是，地下岩石取出后发生应力解除从而产生许多张性微裂缝，当对岩石进行应力重新加载，使张开的微裂缝重新闭合，计算微裂缝闭合过程中不同方向上对应的应变张量，即可计算地应力的大小^[26]。共对义176区块4口井的4块样品开展了差应变测试，测试结果见表 2。

从井点实测的地应力数据可见，义176区块现今地应力变化范围较大，最大水平主应力为55.21$\sim$71.08 MPa，最小水平主应力为46.55$\sim$ 58.99 MPa，垂向主应力为78.61$\sim$93.66 MPa。然而，离散的点数据难以全面反映义176区块整体的地应力场分布特征，因此，以这些井点数值为约束，开展应力场的三维预测。

$ \left\{ \begin{array}{l} {S_{\rm v}} = {\sigma _ \bot }\\ {S_{\rm H}} = \dfrac{{{\sigma _x} + {\sigma _y}}}{2} + \sqrt {{{\left( {\dfrac{{{\sigma _x} - {\sigma _y}}}{2}} \right)}^2} + {{({\sigma _{x45y}} - \dfrac{{{\sigma _x} + {\sigma _y}}}{2})}^2}} \\ {S_{\rm h}} = \dfrac{{{\sigma _x} + {\sigma _y}}}{2} - \sqrt {{{\left( {\dfrac{{{\sigma _x} - {\sigma _y}}}{2}} \right)}^2} + {{({\sigma _{x45y}} - \dfrac{{{\sigma _x} + {\sigma _y}}}{2})}^2}} \end{array} \right. $

(1)

式中：

$\sigma _ \bot$—$z$方向试样的应力，MPa；

$\sigma _x$—$x$方向试样的应力分量，MPa；

$\sigma _{x45y}$—45°方向试样的应力，MPa；

$\sigma _y$—$y$方向试样的应力，MPa；

$S_{\rm H}$—最大水平主应力，MPa；

$S_{\rm h}$—最小水平主应力，MPa；

$S_{\rm v}$—垂向主应力，MPa。

3 基于井点约束的三维应力场智能预测 3.1 三维应力场智能预测流程

第一步：构建三维岩石力学场。由于义176区块目的层在平面及垂向上具有较强的非均质性，传统平面均质、纵向分层的力学模型无法满足研究要求。本文通过地震层速度确定三维岩石力学参数，计算方法为

$ {E_{\rm d}} = {\rho _{\rm b}} v_{\rm s}^2 \dfrac{{3v_{\rm p}^2 - 4v_{\rm s}^2}}{{v_{\rm p}^2 - v_{\rm s}^2}} $

(2)

$ \mu_{\rm{d}} = \dfrac{{v_{\rm p}^2 - 2v_{\rm s}^2}}{{2(v_{\rm p}^2 - v_{\rm s}^2)}} $

(3)

式中：

${E_{\rm d}}$—杨氏弹性模量，MPa；

${\mu_{\rm d}}$—泊松比，无因次；

${\rho_{\rm b}}$—岩石密度，kg/m$^3$；

${v_{\rm p}}$—纵波波速，m/s；

${v_{\rm s}}$—横波波速，m/s。

第二步：建立有限元模型。为了提高建模精度，采用Petrel与ANSYS联合建模的方法，首先，利用Petrel进行地质建模，包括层面模型和断层模型；再将层面与断层数据导出，借助autoCAD重构曲面及曲线，接着利用一套自编的插件将模型转化为IGES格式，使ANSYS能够直接识别；然后，进行有限元网格的划分，综合考量模拟精度和运算效率，将断层与地层的网格步长分别设置为150和300，共划分出335 422个节点，1 944 398个单元；最后，将第一步中得到的三维岩石力学参数(包括弹性模量、泊松比及岩石密度)赋予到每一个单元网格中(图 3)。

第三步：最优边界条件确定。首先，根据义176区块所处的大地构造背景初步确定边界条件，设置应力边界条件的初始值；再以图 2和表 2中实测现今地应力方向与大小为约束，编写一种自适应边界程序，不断自动修正边界条件，使模拟结果与实测结果之间的误差达到最小(图 4)。最终以此边界条件进行应力场的数值模拟。

3.2 自适应程序的编译

井点实测值$\mathit{\boldsymbol{S}}$与预测值$\mathit{\boldsymbol{x}}$存在如下关系

$ {\mathit{\boldsymbol{S}}_i} = {\mathit{\boldsymbol{A}}_i}{\mathit{\boldsymbol{x}}_i} + {\mathit{\boldsymbol{e}}_i} $

(4)

式中：

$\mathit{\boldsymbol{A}}_i$—系数矩阵；

$\mathit{\boldsymbol{S}}_i$—井点实测应力矩阵；

$\mathit{\boldsymbol{x}}_i$—预测值矩阵；

$\mathit{\boldsymbol{e}}_i$—实测值与预测值之间的误差；

$i$—约束井点的数量。

在整个计算过程中，以井点预测值与实测值的差值函数作为目标函数

$ \phi ({x_j}) = \sum\limits_{j = 1}^M {e_i^2} $

(5)

$\phi$—目标函数；

$j$—误差累积的数量；

$M$—达到临界精度的最大约束井数量。

当目标函数的精度达到要求时，确定合理的边界条件；当目标函数未达到精度要求时，在构造应力宏观分析结果的容忍范围内，对边界条件以一定步长进行不断自动修正，直到目标函数精度达到要求为止。

程序主要包含边界条件加载、数值对比及判断语句循环3个模块。本次研究中，设置初始边界为：北边施加50 MPa压力，南边施加40 MPa压力，东边施加80 MPa压力，西边施加10 MPa的左旋剪切，模拟郯庐断裂的影响。顶面除上覆岩体的重力外，再加30 MPa压力。将表 2的井位坐标和实测应力值输入程序，最大误差设置为10%，应力边界条件为10$\sim$60 MPa，变化量设置为2 MPa。

通过了30次循环，耗时3 h 25 min，每个井点的相对误差控制在10%以内，达到了预设精度。此时的边界条件为：北边为46 MPa的压力，南边为23 MPa的压力，东边为90 MPa压力，西边的剪切力和顶面的压力保持不变，仍分别为10 MPa和30 MPa(图 5)。

4 预测结果及分析

以自适应边界程序最终修正结果为边界条件，通过模拟计算得到义176区块三维应力场分布，包括最大水平主应力($S_{\rm H}$)、最小水平主应力($S_{\rm h}$)、垂向主应力($S_{\rm v}$)及水平应力差($S_{\rm H}-S_{\rm h}$)的方向与数值，且与井点实测值的误差控制在了10%以内。

4.1 现今地应力数值

义176区块应力场整体呈带状分布，应力值有北高南低的趋势，主应力数值分布并不均匀，离散程度较高，多见“牛眼状”的零星展布。最大水平主应力主要在50$\sim$90 MPa，最小水平主应力为35$\sim$75 MPa，垂向主应力为55$\sim$100 MPa。靠近断层，3个主应力值减小，且减小的程度与断层规模有关。在南北两条大型边界断层附近，主应力值大约减小40%，而内部小断层附近，主应力值约减小20%。义176区块水平应力差总体不超过25 MPa，低值主要呈带状分布于区块中部，这也是井位较为密集的部位(图 6)。较低的水平应力差有利于压裂作业的开展^[27]。

剖面上，3个主应力与埋深呈线性关系(图 7，式(5))，且表现为$S_{\rm v}>S_{\rm H}>S_{\rm h}$，属于I$_{\rm a}$类地应力，其中，垂向主应力与埋深的关系最为理想，相关程度最高，这是由于垂向应力只与上覆岩石的重力有关，而水平应力除了与区域构造背景有关，还受岩石的非均质性、构造形态等因素影响，故水平应力的离散性较高，其相关系数较垂向应力的相关系数低。

$ \left\{ \begin{array}{l} {S_{\rm h}} = 0.0128H + 18.3514\\ {S_{\rm H}} = 0.0163H + 31.7129\\ {S_{\rm v}} = 0.0209H + 41.6177 \end{array} \right. $

(6)

除了3个主应力外，应力系数也是描述地应力状态的重要参数，包括侧压力系数($k$)、最大和最小水平主应力比值($k_{\rm{H/h}}$)、最大水平应力系数($k_{\rm H}$)、最小水平应力系数($k_{\rm h}$)等^{[10, 28-30]}。其中侧压力系数$k$和最大和最小水平主应力比值是最常用的两个参数，其表达式如下

$ \left\{ \begin{array}{l} {k_{\rm{H/h}}} = {S_{\rm H}}/{S_{\rm h}}\\ k = ({S_{\rm H}} + {S_{\rm h}})/(2{S_{\rm v}}) \end{array} \right. $

(7)

侧压力系数是岩体侧向平均应力与垂向应力的比值，是决定岩体稳定性的一个重要参数，其值一般小于1，且随埋深增大而降低。最大和最小水平主应力比值能够反映岩体所受水平应力的状态，可以体现井壁的水平荷载，与水平应力差类似，其数值越小(越近于1)，岩体越稳定。

义176区块的侧压力系数呈西低东高的分布趋势，除义180井处于侧压力系数的高值区，其余井点均位于侧压力系数的低值区，其值普遍低于0.8。最大和最小水平主应力比值呈条带状分布，具有南北高、中部低的规律，义176区块内的井点均位于$k_{\rm{H/h}}$的低值区，数值一般不高于1.4(图 8)。从地应力分布的角度来看，较低的$k$与$k_{\rm{H/h}}$是井位部署的优势区域，对井筒的稳定有积极影响，另外，若要在这种低渗透油藏开展压裂作业，具有低$k$与$k_{\rm{H/h}}$的部位易于形成复杂网状裂缝，也是压裂作业的有利区，该区可称为“应力甜点”。

4.2 地应力方向

义176区块最大水平主应力方向总体为近EW—SEE，西部主要呈SEE向，向东逐渐转为近EW。得益于自适应边界程序的应用，预测结果与实测结果的吻合程度较高。可以看出，靠近边界断层，最大水平主应力方向顺断层走向偏转，而研究区块内部，应力方向变化微弱，偏转一般不超过10°(图 9)。最小水平主应力与最大水平主应力垂直，方向总体为近NS—NNE。

地应力方向会受断层、岩相及构造形态的控制，一般断层的影响较为显著，其规模、走向、充填物以及形态都会引起地应力方向的偏转^[30]，并且往往会造成应力方向的杂乱，而岩相的非均质性及构造的起伏对应力方向影响比较微弱且均匀。

5 结论

(1) 义176区块现今最大水平主应力方向总体呈近EW—SEE向，现今地应力数值变化较大，在目的层范围属于I$_{\rm a}$类地应力，平面上呈带状展布。

(2) 义176区块中部存在一片条带状展布的“应力甜点”，其水平应力差和$k_{\rm{H/h}}$值较低，是井位部署和压裂作业的优势区域。

(3) 自适应边界程序大幅缩短了人工反复尝试边界条件的时间，提高了应力场模拟的效率，结合Petrel与ANSYS的联合建模法及三维非均质力学模型，进一步完善了三维应力场预测的方法，形成了一套应力场的智能预测技术。