0 引言

地下介质中应力不仅是控制油气运移^[1]、压实成岩作用^[2]、构造演化^[3]等的关键因素，同时也对储层压裂^[4]、注水开发^[5]、钻井井壁稳定性^[6]、钻井安全^[7]等工程作业影响较大。因此，地应力分析受到业界的关注。Fillippone^[8]基于地层速度与地层有效应力的关系，综合不同实际资料给出了地层速度与地层压力之间的经验关系式，提出了一种利用地震资料预测地层压力的有效方法，该方法在油气勘探中得到进一步发展和应用^[9-10]；Dillen^[11]建立了反射系数与应力的关系，直接利用地震资料预测地应力，并逐渐形成技术体系^[12-13]。

地层异常压力预测一般基于以下认识：当地层中孔隙流体在特定地质条件下承担了部分上覆地层重力时，流体压力增大，同时地层有效应力相对降低，导致地层压实作用减弱，速度减小。另外，根据岩石力学相关认识可知，在地层中存在构造变化(起伏或断裂)时，上覆介质重力会不均匀作用于地层^[14]，此时地层有效应力发生横向变化，必然影响地层速度。文献[15]初步探讨了上述问题，指出断层下降盘的滚动背斜的速度急剧变化可能部分源于应力变化引起的差异压实作用，实例分析证明这种速度变化严重影响同类构造的勘探、开发。截至目前，未见地应力横向变化对地层速度及构造解释影响的定量分析文献。为此，本文基于测井数据，参考了土木工程、固体矿产、钻井工程等领域的研究成果^[16-19]，采用有限差分法模拟应力场，分析局部地应力对地层速度的影响机理，以提高构造解释精度。

1 地应力与地层速度的关系 1.1 测井数据统计

在沉积环境相对简单、构造平缓、无断裂影响的情况下，同一位置、不同深度的垂向地应力差ΔP_z主要来源于重力，其与地层密度ρ的关系为

$ \Delta P_z=\int_{z_1}^{z_2} \rho(z) g \mathrm{~d} z $

(1)

式中：z为深度，z₁、z₂分别为数据统计层段的起、止深度；g为重力加速度，取9.81m/s²。实际应用中，可基于测井密度数据估算ΔP_z，其离散表达式为

$ \Delta P_z \cong \sum\limits_{i=1}^m \rho_{z_i} g \Delta z $

(2)

式中：Δz为样点的深度间隔；m为样点数；ρ_zi为深度z_i处的地层密度。

若以特定深度为基准点，可利用式(2)及测井数据统计速度随垂向地应力的变化。针对地应力与地层速度之间的关系，选取A区典型钻井数据进行分析。以HZ-A1井为例，该井在深度大于1200m后井径变化较小，对应数据较可靠，因此将1200m深度作为该井数据的参考点。为便于进一步统计、分析，对测井数据进行整理，包括：①通过重采样将不同数据深度对齐(采样率为2m)，并且剔除异常点；②利用密度数据求取任意深度点相对参考点的垂向地应力差，将时差数据转换为速度数据；③利用泥质含量区分砂(泥质含量小于25%的碎屑岩)、泥(泥质含量大于75%的碎屑岩)岩；④绘制地层纵波速度随相对地应力变化曲线(图 1a)，可见，地层速度近似与垂向地应力呈线性正相关关系，而且砂、泥岩差别较小，斜率(应力对地层速度的影响系数，以下简称影响系数)约为40m·s^-1·MPa^-1。

A区其他井的砂、泥岩纵波速度随相对地应力的变化与HZ-A1井类似，如HZ-B1井(图 1b)、PY-A1井(图 1c)、XJ-A1井(图 1d)。表 1为A区影响系数统计表，可见：①同一井点，地应力变化对砂、泥岩速度的影响接近；②不同井点，地应力变化对地层速度的影响略有差异；③影响系数分布于30~40m·s^-1·MPa^-1。

对于埋深较小、范围较小的靶区，可忽略沉积环境、胶结作用等其他因素的影响，上述认识也适用于横向地层。速度与垂向地应力之间的线性正相关关系可简述为：速度增量Δv与应力差ΔP_z之间呈正比例关系，即

$ \Delta v=\alpha \Delta P_z $

(3)

式中α为影响系数，本文取分析井点的平均值(约为37m·s^-1·MPa^-1)。对于地质情况复杂或大尺度区域，可以认为α值随空间位置变化，但式(3)仍然成立。

1.2 地应力与速度换算

基于式(3)，假设未压实碎屑岩对应的平均纵波速度为v₀，可估算地下任意点介质的速度期望值

$ v=v_0+\alpha P_z $

(4)

式中：v为层速度；P_z为垂向地应力。进一步通过

$ \frac{Z}{v_{\mathrm{a}}}=\frac{T_Z}{2}=\int_0^z \frac{1}{v_0+\alpha P_z} \mathrm{~d} z $

(5)

获得平均速度。式中：T_Z为深度Z处的地震波双程旅行时；v_a为0~Z深度段的平均速度。式(5)的离散形式为

$ \frac{N_Z}{v_{\mathrm{a}}}=\sum\limits_{i=1}^{N_Z} \frac{1}{v_0+\alpha P_{z, i}} $

(6)

式中：N_Z为地面至深度Z处的样点数；P_{z, i}为第i个样点处的垂向地应力。

基于式(6)将地应力场转化为平均速度场，可以指导精细构造解释。

2 地应力模拟

采用连续介质条件下的有限差分法阐述二维介质地应力场模拟的基本原理及实现过程。

2.1 基本原理

一般情况下地下介质同时具有弹性和塑性，在平衡状态下，介质的塑性无法恢复，其应力源于弹性形变，应力和应变之间服从胡克定律。假设介质为各向同性介质，在二维情况下存在以下关系

$ \left\{\begin{array}{l} \sigma_x=(\lambda+2 \mu) \varepsilon_x+\lambda \varepsilon_z \\ \sigma_z=(\lambda+2 \mu) \varepsilon_z+\lambda \varepsilon_x \\ \tau_{x z}=\mu \varepsilon_{x z} \\ \lambda=K-\frac{4}{3} G \\ \mu=G \end{array}\right. $

(7)

式中：λ、μ为拉梅常数；K、G分别为介质的体积模量、剪切模量；σ_x、σ_z分别为x、z方向的正应力；τ_xz为切应力；ε_x、ε_z分别为x、z方向的正应变；ε_xz为切应变。弹性介质在外力作用下会发生位移，其在x、z方向的位移分量分别为u(x, z)和w(x, z)(下文简写为u和w)，应变与位移之间有如下关系

$ \left\{\begin{array}{l} \varepsilon_x=\frac{\partial u}{\partial x} \\ \varepsilon_z=\frac{\partial w}{\partial z} \\ \varepsilon_{x z}=\frac{\partial u}{\partial z}+\frac{\partial w}{\partial x} \end{array}\right. $

(8)

另外，在平衡状态下，地下介质受应力、重力共同作用，应力分量满足运动平衡方程

$ \left\{\begin{array}{l} \frac{\partial \sigma_x}{\partial x}+\frac{\partial \tau_{x z}}{\partial z}=0 \\ \frac{\partial \tau_{x z}}{\partial x}+\frac{\partial \sigma_z}{\partial z}=-\rho g \mathrm{~d} z \end{array}\right. $

(9)

根据式(7)~式(9)，得到以下二阶偏微分方程组

$ \left\{\begin{array}{l} (\lambda+2 \mu) \frac{\partial^2 u}{\partial x^2}+\mu \frac{\partial^2 u}{\partial z^2}+(\lambda+\mu) \frac{\partial^2 w}{\partial x \partial z}+ \\ \frac{\partial(\lambda+2 \mu)}{\partial x} \frac{\partial u}{\partial x}+\frac{\partial \lambda}{\partial x} \frac{\partial w}{\partial z}+\frac{\partial \mu}{\partial z}\left(\frac{\partial w}{\partial x}+\frac{\partial u}{\partial z}\right)=0 \\ (\lambda+2 \mu) \frac{\partial^2 w}{\partial z^2}+\mu \frac{\partial^2 w}{\partial x^2}+(\lambda+\mu) \frac{\partial^2 u}{\partial x \partial z}+ \\ \frac{\partial(\lambda+2 \mu)}{\partial z} \frac{\partial w}{\partial z}+\frac{\partial \lambda}{\partial z} \frac{\partial u}{\partial x}+\frac{\partial \mu}{\partial x}\left(\frac{\partial w}{\partial x}+\frac{\partial u}{\partial z}\right)=-\rho g \mathrm{~d} z \end{array}\right. $

(10)

求解式(10)得到u和w，依次代入式(8)和式(7)得到地应力场。在地应力场尺度下地下岩石一般为非均质介质，弹性常数的偏导数不为0，因此需要保留式(10)的一阶微分项。

2.2 边界条件

基于应力仅受重力和构造影响且处于平衡状态的假设条件，给定模型边界条件如下：

(1) 上边界设定为自由边界，即横向位移、垂向位移的法向导数均赋0，为第二类边界条件；

(2) 左、右边界设定为横向固定边界，即横向位移赋0，垂向位移的法向导数赋0，为混合边界条件；

(3) 下边界设定为垂向固定边界，即垂向位移赋0，横向位移的法向导数赋0，为混合边界条件。

2.3 模型参数

建立二维地质模型需要ρ、K、G共3个参数模型，各参数取值参考A区HZ-B1井测井数据(井筒长度接近4km)，参数模型建立过程如下。

(1) 绘制ρ、K、G随深度变化曲线，并拟合趋势线(图 2)。数据中存在个别异常点，部分可能与特殊岩性有关，但对趋势线影响较小，这里不详细讨论。

(2) 利用拟合趋势线对应关系，给定模型中每个深度点的参数值作为基础参数值。

(3) 基于基础参数设计构造格架，对每套地层在一定范围内给定增量(变化量为幅值的10%)。

(4) 基底为致密岩石，ρ赋常值，K、G分别为50、30GPa。

为消除重力变化的影响，所有模型的ρ均采用基础参数值。以水平层状介质为例，建立水平层状地质模型(图 3)。可见，该模型中每一质点上覆地层重力随深度逐渐增加而增大，当埋深为4000m时，重力场最大(约87MPa，图 3d)。

3 典型构造影响分析 3.1 披覆背斜构造 3.1.1 实钻速度特征

披覆背斜是在古隆起上沉积的地层经压实作用形成的背斜构造，一般由深至浅构造幅度逐渐减小，且倾角与基底相差较大，但走向、倾向基本一致。实钻结果表明，该类构造翼部的地层速度通常小于高点。如HZ-A油田已进入开发后期，井控程度较高，基于井点等T₀值(图 4b)与顶面实钻深度(图 4b的蓝色实心点的红色数值)获得各井点位置v_a(图 4c的紫色数值)，可见翼部的井点v_a显然小于高部位井点，这种速度变化由不同时深关系^[20]体现。利用单井时深关系(图 5a，方法1)及H1油藏局部沿层时深关系(图 5b，方法2)预测的构造含油面积分别为1.22(图 6a)、3.95km²(图 6b)，对应的预测储量分别为X万m³、2.72X万m³，而根据该油藏实际生产动态反推H1油藏储量约为2.85X万m³，因此方法2的结果与实际更吻合。显然，对于少井油田，由方法1预测的构造含油范围更小、深度预测误差更大，会明显低估油田规模(图 6a、图 6b)。将平面上距油田5km范围内的探井数据纳入方法2统计结果(图 5c，区域时深关系——方法3)，预测储量为1.04X万m³(图 6c)，与方法1的结果接近，也低估了油田规模。

方法1的结果反映了正常压实情况下地层速度与深度的关系，即速度随埋深增加而增大。方法2时深关系斜率更小，反映了局部范围内速度变化较快，构造翼部速度减小。方法3与方法1的时深关系斜率接近，表明较大区域范围内沿层地层速度变化受正常压实规律控制。本文认为局部范围内的速度变化主要受局部应力场变化的影响，且对勘探初期的油气田(特别是海上油气田，常用方法1或方法3计算储量)储量评估影响较大。

3.1.2 应力模拟分析

参考图 4a建立地质模型模拟应力，弹性参数按前文方法给定。图 7为HZ-A油田地应力模拟及速度预测结果。可见：HZ-A1井至HZ-A2井之间基底介质的起伏幅度约为80m(图 7a、图 7b)，根据σ_z局部剖面(图 7c)，采用式(4)~式(6)获得v_a局部剖面(图 7d)及z=2460m的v_a切片(图 7e)，其中HZ-A2井的v_a较HZ-A1井约小8m/s，与HZ-A1井(图 7f)和HZ-A2井(图 7g)层位标定结果估算的速度差(约7m/s)接近，误差较小。因此，基于应力模拟预测的速度变化与实钻数据基本吻合。

3.1.3 机理模型分析

为了解披覆背斜的应力变化特征、控制因素，设计了模型2(图 8a、图 8b)、模型3(图 8c、图 8d)进行地应力模拟。

模型2地应力模拟结果(图 9)表明：①σ_z在构造高点集中且数值较大，在翼部减小且变化较快(图 9a)；σ_z在浅层变化较缓、深层变化较快，在古构造附近较集中。②σ_x在基底附近集中，在古潜山的高点较大，向翼部快速减小，在外围逐渐增大(图 9b)。③τ_xz在古潜山翼部集中，最大约为5MPa(图 9c)。

模型3地应力模拟结果与模型2类似，文中不再赘述。图 10为z=1700m、z=2500m的σ_z切片。由图可见，模型2、模型3的σ_z横向变化基本一致，说明该类构造局部应力变化主要受古构造控制。综合模型2、模型3的地应力模拟结果可知：古构造高点上覆地层在τ_xz作用下分担了较多的重力，造成σ_z较高；古构造翼部上覆地层重力被其他位置地层分担，导致σ_z降低。

分别利用式(4)和式(6)将σ_z(图 9a)转化为v(图 11a)和v_a(图 11b)(v₀=2000m/s)，并提取z=1700m(图 11c)、z=2500m(图 11d)的v_a切片。可见，该类构造在地应力影响下导致构造高部位v及v_a较翼部大，且埋深越大变化越快。

3.2 滚动背斜构造 3.2.1 实钻速度特征

除披覆背斜外，发育在正断层上盘的小型背斜构造(滚动背斜或逆牵引背斜)也可能受压实作用影响，文献[15]分析了南海东部该类构造局部范围的速度变化规律，得到如下经验关系式

$ V_{x, y} \approx V_0-k H D_{\mathrm{F}} \ln \frac{L_{x, y}}{L_0} $

(11)

式中：V_{x, y}为水平方向平均速度；H为深度；V₀为参考井平均速度；L₀为参考井距断层垂直距离；L_{x, y}为(x, y)点距断层垂直距离；D_F为断层最大断距；k为经验常数(63.38m·s^-1·km^-2)。参考A区滚动背斜油田HZ-D实钻数据(图 12a)，设定H、V₀、L₀、D_F分别为2510m、2497.12m/s、890m、205m以及2990m、2632.62m/s、572m、205m预测V_{x, y}变化，统计井点位置V_x，y变化并叠合显示(图 12b、图 12c)，预测与实钻结果基本一致。可见，在A区断层上盘的V_{x, y}具有向断层方向增大、埋深越大且越靠近断层变化越快的特点。

3.2.2 机理模拟分析

断层影响致使应力场特征复杂，难以准确获取应力场分布。受技术水平限制，本文基于理想假设条件，采用数值模拟方法分析应力场机理：①剩余应力、区域应力在研究范围对地层的影响是均匀的，故模拟中不予考虑；②正断层引起断层两盘地层之间的剪切应力释放，相互作用减小；③将断层等效为弹性模量小于围岩且具有一定宽度的弹性介质(传递介质)；④将整个下盘作传递介质，以突出(滚动背斜发育的)上盘的应力分布特征。

基于上述假设条件建立滚动背斜地质模型(图 13a、图 13b)，加载重力后模拟地应力得到σ_z(图 13c)。结合式(10)、式(12)得到v(图 14a)和v_a(图 14b)，并提取z=2500m(图 14c)和z=3000m(图 14d)的v_a切片，可读取同样1000m距离内的速度变化量分别约为70m/s(图 14c)和101m/s(图 14d)。可见，正断层上盘断面附近σ_z快速增大，导致v_a具有向断层方向增大、埋深越大且越靠近断层变化越快的特点。模拟的速度场特征与实钻数据吻合，表明地应力场可能是滚动背斜速度变化的主要因素之一。

3.3 讨论与分析 3.3.1 地应力模拟假设条件

在沉积盆地中应力场与构造运动、压力场、温度场相互作用^[21]，无论在区域尺度还是局部尺度上准确获取应力场分布均是难点问题^[16-18]，本文结合地震数据研究地应力，基于以下较理想的假设条件：①不考虑地层异常压力等复杂情况；②针对研究区的拉张应力背景^[22]，将区域应力场或古应力场视为均匀分布，不影响速度横向变化，因此在构造圈闭尺度上重力可能是局部构造应力变化的主要因素。实例分析表明：在没有异常地层压力情况下，区域沿层时深关系(图 5c)与单井时深关系一致(图 5a)；局部沿层时深关系斜率更小(图 5b)，即局部存在压实变化，间接说明研究区应力分布满足上述假设。

3.3.2 地层速度的应力影响因素

实例分析表明，局部地层速度变化主要取决于应力变化，且是由地层压实作用间接形成的。因此，当地层应力减小时，压实作用减弱，在地层速度减小的同时，储层物性变好(孔隙度增大)，反之亦然。HZ-A油田HZ-A2井的地层速度较小，储层有效孔隙度整体较HZ-A1井高(图 15)，佐证了上述观点。

目前人们认为，影响纵波速度横向变化的应力因素主要为地层异常流体压力。本文综合地应力模拟、速度预测及其实例分析认为，在圈闭或局部构造尺度上，应力变化引起的速度变化可能会严重干扰地震资料解释。因此，有必要将与构造有关的应力变化作为一种影响地层速度的因素予以关注。

4 结束语

文中利用多口钻井数据，统计、分析了研究区碎屑岩地层在压实作用下速度与垂向地应力之间的关系，获得统计学意义上速度随应力变化的趋势，进一步结合地应力模拟讨论了两种典型局部构造的速度分布特征及成因机理，形成以下认识：

(1) 地应力通过压实作用直接影响地层速度，地应力受局部构造影响而变化较快，进一步影响局部速度分布规律及地震资料的构造解释，应予以关注；

(2) 在古构造影响下，披覆背斜构造翼部地应力相对较小，导致翼部速度较低，实际构造幅度可能小于常见构造预测方法的预测结果；

(3) 在断层影响下，正断层下降盘的构造应力向断层方向增大，导致远离断层速度减小，真实构造高形态与常规方法预测结果出现差异。

尚需指出，本文基于二维模型研究局部地应力对地层速度的影响机理，仅可定性或半定量分析圈闭或局部构造尺度的速度变化，若要获得准确的应力—速度关系，还需考虑三维区域应力、应力变化史等因素。