引言

全球发现的油气藏中有大量的裂缝性油气藏，包括碳酸盐裂缝型油气藏、 致密砂岩裂缝型油气藏及基岩裂缝型油气藏(毛宁波等，2008; 杨晓等，2010)． 裂缝通常是裂缝性油气藏的储存空间或运移通道，因此裂缝介质中地震波场的研究越来越受到关注(谢桂生，1994; 金抒辛，何樵登，2005; 潘保芝等，2006)． 研究发现，在含有定向分布的裂缝介质中地震波的传播呈现各向异性的特征(张中杰，何樵登，1989； 郭建，1993； Schoenberg，Sayers，1995)． 因此，利用地震波场来研究裂缝发育带位置、 确定储层物性并进行储层预测具有重要的意义．

目前研究裂缝的方法之一是将裂缝看成散射体，用散射理论研究裂缝介质的地震响应(刁顺等， 1994a，b)，但此方法的正演效率非常低. 本文用等效理论将裂缝介质等效为各向异性介质，运用各向异性介质理论研究裂缝的地震响应． 各向异性是指介质的弹性参数随方向而变化，该方法效率高，能得到较好的效果． 裂缝对地震波的影响已研究了多年． 例如： Hudson(1981)提出裂缝介质中P波和S波速度与衰减的公式； Thomsen(1986)提出利用弹性弱各向异性理论研究定向裂缝与等径孔隙度对地震速度的综合影响； Schoenberg和Sayers(1995)使用裂缝法向与切向柔度推导裂缝介质中的有效弹性常数; Bakulin等(2000a，b，c)回顾了几种不同的裂缝等效理论，并比较了它们对裂缝参数估计的可行性．

地震波场数值模拟是地球物理学领域的一个研究热点，既可以提高我们对各种复杂介质中地震波传播规律的认识，又可以作为检测工具检验各种方法技术的应用效果． 前人研究结果表明，地震波在定向裂缝介质中的传播具有横向各向同性的特征(Hudson，1981； Thomsen，1986； Bakulin等(2000a，b，c． Crampin(1985)提出构造应力产生的具有水平对称轴的横向各向同性(horizontal transverse isotropy，简称为HTI)介质，可用来模拟高陡倾角裂缝介质(孔丽云等，2012)． 目前常用的裂缝介质正演模拟方法主要是交错网格有限差分方法(董良国等，2000； 裴正林， 2004，2006； 裴正林，王尚旭，2005； 李雪静，刘定进，2008)． 普通交错网格能有效地模拟地震波在裂缝介质中的传播，但也存在一定的缺点. 首先，当模拟非均匀性比较强的介质或各向异性介质时，为了获得满意的计算结果，需要对介质参数进行平均或内插; 其次，在遇到自由界面时，必须对自由界面进行单独的显式处理(严洪勇，刘洋，2012)． 旋转交错网格将相同的物理量定义在相同的网格点上，计算其沿对角线的差分，然后利用对角线差分组合计算其沿坐标轴的差分. 该方法不需要对介质参数进行平均或内插，在遇到自由边界时不需要单独的显示处理(Saenger et al，2000)，从而解决了普通交错网格存在的缺点，是目前正演中比较准确的方法． 本文采用旋转交错网格正演模拟方法对裂缝性储集层进行正演模拟.

1.   裂缝介质的等效理论

运用各向异性理论进行裂缝介质正演模拟，需要将裂缝介质等效为各向异性介质． 目前广泛存在的裂缝等效理论主要有Thomsen等效理论、 Hudson等效理论及线性滑动理论等． Thomsen等效理论相当于低频模型，假设流体流动无限慢，该方法限制条件较多，不利于研究(Thomsen，1986)． Hudson等效理论与Thomsen等效理论正好相反，相当于高频模型，其假设流体流动无限快，也不利于实际研究(Hudson，1999)． 线性滑动理论(Bakulin等(2000a，b，c)是将裂缝看作地层中的一个柔性面，该柔性面可以用裂缝的切向柔度和法向柔度进行表示，柔度符合线性滑动边界条件便于求解. 对于多组裂缝可以对裂缝柔度进行矢量叠加便于实际求解． 因此本文采用线性滑动等效理论将裂缝介质等效为各向异性介质．

含有裂缝岩石的总柔度等于裂缝柔度与围岩柔度的和，即

式中，S为岩石总的柔度矩阵. S_f为裂缝的柔度矩阵. S_b为围岩的柔度矩阵. 其中柔度矩阵为弹性矩阵 C 的倒数，即 S = C^－1. 根据线性滑动理论，裂缝介质的弹性矩阵为

式中，λ和μ为背景岩石的拉梅常数，ΔN和ΔT分别为裂缝介质的法向柔度与切向柔度. 当介质中存在几组裂缝时，先求出每组裂缝介质的切向柔度与法向柔度，然后对切向柔度与法向柔度进行矢量叠加求出裂缝介质总的柔度矩阵． 裂缝介质等效示意图如图 1所示．

图  1  裂缝介质等效示意图图中ΔN和ΔT分别为裂缝介质的法向柔度与切向柔度.(a)单缝模型；(b)裂缝组合模型

Figure  1.  The sketch of fracture equivalenceΔN denotes the normal weakness， and ΔT denotes the tangential weakness(a)Sigle fracture model;(b)Combination model of fractures

下载: 全尺寸图片 幻灯片

由于裂缝介质的法向柔度与切向柔度需要满足线性滑动边界条件，不方便求解． 而采用Hudson等效理论与线性滑动理论相结合，既弥补了Hudson等效理论的局限性，又可以方便计算裂缝的法向柔度与切向柔度. 实际计算中采用线性滑动理论与 Hudson理论相结合的等效理论将裂缝介质等效为各向异性介质． 此时，裂缝的法向柔度和切向柔度分别为

式中，e为裂缝密度，k′和μ′分别为裂缝充填物的体积模量与切变模量，k为背景介质的体积模量，a为裂缝的长度，c为裂缝的宽度，g=v_S²/v_P²． 弹性介质模型的性质由刚度矩阵 C 确定，C 确定了应力与应变之间的关系，但其确定弹性波动方程系数的物理意义很不直观． Thomsen(1986)提出了一整套表征各向异性的参数，具体表达如下：

式中V_P0，V_S0分别为qP波和qS波沿各向异性介质对称轴方向的相速度； ε，δ和γ是表示HTI介质各向异性强度的3个无量纲因子． 根据式(1)-(4)可得裂缝参数与各向异性参数之间的关系为

将线性滑动理论与Hudson等效理论相结合，对4组典型的裂缝介质进行简单的等效，其结果如表 1所示． 裂缝层埋深3 000 m，岩层厚度为100 m，背景岩石的纵波速度为6 200 m/s，横波速度为3 500 m/s，密度为2 800 kg/m³．

表  1  裂缝介质参数

Table  1.  Fracture parameters

下载: 导出CSV 

| 显示表格

2.   HTI介质一阶速度应力方程

式(2)—(3)可将裂缝参数等效为表征裂缝介质的弹性矩阵参数，即将裂缝介质等效为各向异性介质． 下面推导HTI介质的一阶速度应力方程． 由广义胡克定律、 运动微分方程及几何方程可得HTI介质的波动方程．

广义胡克定律

运动微分方程

几何方程

This page contains the following errors:

Below is a rendering of the page up to the first error.

3.   裂缝介质旋转交错差分方程

3.1   时间2阶差分方程近似

3.2   空间高阶差分方程近似

式中，Δx为x方向的采样间隔，c_n为有限差分系数，具体如表 2所示．

表  2  几种常用有限差分系数

Table  2.  Several kinds of finite difference coefficients

下载: 导出CSV 

| 显示表格

3.3   一阶速度应力差分方程

交错网格主体网格定义的变量与交错网格半定义的变量不同，速度、 应力和弹性参数的定义规则如图 2所示. 其中变量分别定义在网格1，2，3，4点上.

图  2  标准交错网格示意图

Figure  2.  The sketch of st and ard staggered grid

下载: 全尺寸图片 幻灯片

在普通交错网格有限差分中，通常将介质参数定义在法向应力所在的网格点上． 本文中法向应力及所有的背景参数都定义在主体网格1点上，z方向速度定义在2点上，x方向速度定义在3点上，切应力定义在4点上. 应用普通交错网格需要定义1套主网格点和3套半网格点，在模拟裂缝介质传播时需要内插相关的弹性参数． 应力速度及背景参数的定义规则如表 3所示． 其中: 正应力τ_xx，τ_zz以及弹性参数c₁₁，c₁₃，c₅₅，c₃₃，ρ定义在主体网格1点上； 切应力τ_xz定义在4点所示的半网格点上，z方向速度v_z定义在2点所示的半网格点上； x方向速度v_x定义在3点所示的半网格点上； 半网格3点的密度需要用密度差值移位或者取平均得到．

表  3  标准交错网格中弹性波场分量和弹性参数的空间位置

Table  3.  The positions of elastic wavefield components and elastic parameters in st and ard staggered grid

下载: 导出CSV 

| 显示表格

当介质非均质非常强时，普通交错网格有限差分得不到满意的结果，此时需要用旋转交错网格模拟． 旋转交错网格首先将图 2沿对角线旋转，使得所有的应力分量位于同一点上，所有的速度位于同一点上． 本文的变量定义规则如图 3所示． 每个点代表的变量含义如表 4所示.

图  3  旋转交错网格示意图

Figure  3.  The sketch of rotated staggering grid

下载: 全尺寸图片 幻灯片

表  4  旋转交错网格中弹性波场分量和弹性参数的空间位置

Table  4.  The positions of elastic wavefield components and elastic parameters in rotated staggering grid

下载: 导出CSV 

| 显示表格

旋转交错网格中所有速度分量和应力分量均位于相同的网格点上，因此介质密度与弹性参数均位于速度所在的网格点或者应力所在的网格点上． 旋转交错网格中只需要定义一套主体网格和一套半网格，在模拟裂缝介质地震波场传播时不需要对弹性参数进行内插，从而提高了稳定性． 旋转交错网格中速度应力参数及背景参数如表 4所示．

表 4给出了应力速度及背景参数的定义规则. 其中正应力τ_xx和τ_zz、 切应力τ_xz以及弹性参数c₁₁，c₁₃，c₅₅，c₃₃，ρ均定义在主体网格1点上，x方向速度v_x与z方向速度v_z定义在2点所示的半网格点上． 旋转交错网格将相同的物理量定义在相同的网格点上，计算其沿对角线的差分，再利用对角线差分组合计算其沿坐标轴的差分. 此方法不需要对介质参数进行平均或者差值． 下面推导二维情况下旋转交错网格沿坐标轴方向的表达式.

根据矢量合成有

This page contains the following errors:

Below is a rendering of the page up to the first error.

将式(12)、(13)带入式(9)，得到旋转交错网格下裂缝介质的一阶速度应力方程与一阶应力速度方程为

式(14)中只需要用到对角线方向的差分格式，将对角线的差分格式进行线性组合即可得到最终的结果. 将式(11)带入式(14)得

This page contains the following errors:

Below is a rendering of the page up to the first error.

3.4   正演模拟中需要注意的问题

3.4.1   稳定性条件

对于裂缝介质一阶速度应力方程与一阶应力速度差分方程而言，得到其稳定性的解析式比较困难. 通过矩阵特征分析(Igel et al，1995)可得到其稳定性的一种表达形式. 通过矩阵特征分析法得到裂缝介质一阶速度应力弹性波方程时间(2阶)、 空间(10阶)的稳定性条件为

式中，Δt为时间采样间隔，C_max为弹性矩阵最大值，ρ为密度，Δx和Δz为空间采样间隔，c_i为差分系数(表 2)． 若选取参数不正确，程序不稳定，会使得正演模拟结果错误，消除的方法是根据模型参数通过式(16)选取合适的时间空间采样间隔. 须说明，时间采样点不是越小越好，时间采样点越小则需要的时间点数越多，累计误差越大，也会引起不稳定且计算时间太长，不利于实际研究．

3.4.2   边界条件

采用有限差分进行裂缝介质正演模拟的过程中，由于模拟区域不是无限大，这种人为限定模型计算区域的方法会引起边界问题． 这种边界引起的反射波常常比物理边界更强，在地震记录中会对有效信息产生干扰，因此针对边界问题必须采用一定的边界条件消除这种边界反射(何燕，2008)． 消除边界反射的方法有很多，其中具有代表性的有： Clayton和Engquist(1977)基于旁轴近似提出了吸收边界条件，其在特定的入射角和频率范围内具有较好的吸收效果； Reynolds(1978)提出了透明边界条件，其特点是零角度入射时反射系数为零； Marfurt(1984)提出了海绵边界条件； Berenger(1994)针对电磁波传播情况提出了最佳匹配层(perfectly matcher layer，简称PML)吸收边界条件，其在海绵吸收边界的基础上作了进一步改进，并从理论上证明了该方法可以完全吸收来自各个方向的各种频率波，且不产生任何反射，是目前应用最广泛的边界条件. 本文采用PML边界条件消除边界反射． 该方法的原理是在研究区域的边界上增加一个吸收层，使边界上传入吸收层的波随传播距离的增加成指数衰减，而不产生任何边界反射，可以得到裂缝介质一阶速度应力带有PML边界的差分格式． 旋转交错网格的PML边界条件与普通交错网格的PML边界条件一致. 下面以普通交错网格方程推导PML边界条件，然后再应用于旋转交错网格． 以公式(9)中其中一个方程(省略震源)为例，

根据复数坐标变换

对式(17)两边进行傅里叶变换得

根据波场分解原理，式(19)可分解为

将式(18)带入式(20)，并进行傅里叶反变换得

式(9)中其余方程的推导过程如上，不再一一重复．

3.4.3   网格频散

采用有限差分法存在固有缺陷，即在利用有限差分算子逼近微分算子时会产生误差项，它使得具有不同频率的地震波表现为不同的相速度，从而使地震波发生频散，影响数值模拟的精度和偏移成像效果． 数值频散实质上是一种因离散化求解波动方程而产生的伪波动，它是指组成地震波的不同频率分量以不同的速度传播，随着传播时间的推移地震波扩展为更长波列的现象(孙成禹等，2009)． 时间网格频散表现为相速度超前，空间网格频散表现为相速度滞后． 本文采用通量传输校正方法(Book et al，1975; Boris，Book，1973)消除频散，因为使用空间高阶有限差分方法会压制频散且计算量非常大． 通量传输校正方法主要分为以下几步：

1)输运计算，计算待修正值

2)求n时刻的扩散通量

式中η₁ 为扩散系数.

3)扩散计算

4)反扩散通量计算

式中η₂ 为反扩散系数.

5)限制反扩散通量条件

式中，Δ_i－1= Uⁿ⁺¹_i－ Uⁿ⁺¹_i－1，S =sign(Fⁿ⁺¹_i).

6)反扩散计算

4.   裂缝介质正演模拟分析

采用式(14)、(15)的差分方程，编写程序进行裂缝介质正演模拟，可以得到裂缝介质的正演模拟记录. 其中模拟过程中需要注意边界条件、 频散压制及稳定性条件. 裂缝引起的反射波能量相对于背景介质较小，因而正演模拟中裂缝介质引起的反射系数变化经常被淹没在背景介质中. 为了研究裂缝介质的波场，本文设计一个均匀模型裂缝位于均匀模型中间，如图 4所示． 模型尺寸为4000 m×4000 m，纵向、 横向网格间距均为10 m，纵波速度为3000 m/s，横波速度为1700 m/s，密度为2000 kg/m³. 黑色区域为裂缝所在的区域，其中: 裂缝长度为5 m，宽度为1 mm，密度为2条/m，裂缝充填物为油; 裂缝横向延伸为4000 m、 纵向延伸为100 m. 采用式(14)、(15)的差分方程进行正演模拟. 正演模拟采用雷克子波，子波的主频为25 Hz，时间采样率为1 ms，时间长度为5 s. 正演模拟结果如图 5所示．

图  4  单层模型示意图

Figure  4.  The sketch of single layer model

下载: 全尺寸图片 幻灯片

图  5  不同时刻地震波场示意图.(a)t=0.75 s;(b)t=1.0 s;(c)t=1.6 s

Figure  5.  The sketch of seismic wavefield in different moments.(a)t=0.75 s;(b)t=1.0 s;(c)t=1.6 s

下载: 全尺寸图片 幻灯片

图 5中由于裂缝引起的各向异性较小，qSV波和qSH波在地震波场中无法区分，统一用qS波表示. 从图 5a中可以看出，t=0.75 s时波还没有传播到裂缝处，地震波场中含有两种波，即qP波和qS波. 其中qP和qS分别为准P波和准S波． 从图 5b可以看出，t=1.0 s时qP波已经穿过裂缝层，地震波场中含有qS波、 qPRqP波(qP波入射qP反射波)、 qPRqS波(qP波入射qS反射波)、 qPTqS波(qP波入射qS透射波)及qPTqP波(qP波入射qP透射波). 地震波场比较复杂，可以看出裂缝介质顶底都会产生反射波、 透射波．从图 5c中可以看出，t=1.6 s时qS波已经穿过裂缝层，地震记录中含有qPRqP波(qP波入射qP反射波)、 qPRqS波(qP波入射qS反射波)、 qSRqP波(qS波入射qP反射波)、 qSRqS波(qS波入射qS反射波)、 qPTqS波(qP波入射qS透射波)、 qPTqP波(qP波入射qP透射波)、 qSTqS波(qS波入射qS透射波)及qSTqP波(qS波入射qP透射波). 由此可见，t=1.6 s时地震波场非常复杂，裂缝介质顶底都会产生反射波、 透射波. 地震波穿过裂缝介质时地震波场非常复杂，便于我们利用地震波场的信息反演裂缝的参数．

同样采用图 4模型示意图，模型尺寸改为400 m×400 m，模型采样间隔为10 m×10 m. 采用主频为75 Hz的雷克子波，时间采样率为0.1 ms，时间长度为0.5 s，炮点位于模型中间. 模型纵波速度为3000 m/s，横波速度为1700 m/s. 黑色区域为裂缝所在的区域，其中: 裂缝长度为5 m，宽度为1 mm，密度为2条/m，裂缝充填物为油； 裂缝横向延伸为400 m，纵向延伸为1—24 m． 从地震记录中提取第50道的PP波和SS波记录，如图 6所示.

图  6  不同裂缝延伸长度的单道地震记录.(a)PP 波;(b)SS波

Figure  6.  Seismic records of single trace in different crack growth.(a)PP wave;(b)SS wave

下载: 全尺寸图片 幻灯片

根据正演采用的参数，纵波的波长为40 m，横波的波长为23 m. 由图 6a可知，当裂缝延伸为纵波半个波长时，顶底反射完全分开； 当裂缝延伸为纵波1/4波长时，可以找出顶底反射界面，此时顶底反射波尚未完全分开． 由图 6b可知，当裂缝延伸为横波半个波长时，顶底反射完全分开； 当裂缝延伸为横波1/4波长时，可以找出顶底反射界面，此时顶底反射波尚未完全分开．

为了更清楚地了解裂缝介质对地震波场的影响，从不同的裂缝参数出发研究不同裂缝参数对地震波场的影响. 下面分别研究裂缝密度、 裂缝纵横比、 裂缝充填物等裂缝参数对地震波场的影响．

4.1   裂缝密度对地震波场的影响

采用图 4所示的模型示意图. 裂缝长度为2 m，宽度为1 mm，裂缝充填物为油； 裂缝密度分别为0.5，1，2，4和5条/m； 采用雷克子波进行正演模拟，子波的主频为25 Hz，时间采样率为1 ms，时间长度为5 s. 正演模拟结果如图 7所示． 可以看出，裂缝密度越大，引起的反射波能量越强．

图  7  不同裂缝密度(e)正演模拟的地震记录(a)e=0.5条/m;(b)e=1条/m;(c)e=2条/m;(d)e=4条/m;(e)e=5条/m

Figure  7.  Seismic records in the media with different fracture densities(a)e=0.5 m^-1;(b)e=1 m^-1;(c)e=2 m^-1;(d)e=4 m^-1;(e)e=5 m^-1

下载: 全尺寸图片 幻灯片

4.2   裂缝纵横比对地震波场的影响

采用图 4所示的模型示意图. 裂缝密度为1条/m，裂缝充填物为油，裂缝纵横比分别为1000，2000，5000，7000和10000； 采用雷克子波进行正演模拟，子波的主频为25 Hz，时间采样率为1 ms，时间长度为5 s. 正演模拟结果如图 8所示． 可以看出，裂缝纵横比越小，引起的反射波能量越强．

图  8  不同裂缝纵横比(a/c)正演模拟的地震记录(a)a/c=1 000;(b)a/c=2 000;(c)a/c=5 000;(d)a/c=7 000;(e)a/c=1 0000

Figure  8.  Seismic records in the media with different fracture aspect ratios(a)a/c=1 000;(b)a/c=2 000;(c)a/c=5 000;(d)a/c=7 000;(e)a/c=1 0000

下载: 全尺寸图片 幻灯片

4.3   裂缝充填物对地震波场的影响

采用图 4所示的模型示意图. 裂缝密度为1条/m，长度为2 m，宽度为1 mm，裂缝充填物分别为油、 水和泥岩； 采用雷克子波进行正演模拟，子波的主频为25 Hz，时间采样率为1 ms，时间长度为5 s. 正演模拟结果如图 9所示. 可以看出，裂缝充填物速度与背景速度差别越大，引起的反射波能量越强．

图  9  不同裂缝充填物正演模拟的地震记录(a)油；(b)水；(c)泥岩

Figure  9.  Seismic records in the media with different fracture fillings(a)Oil;(b)Water;(c)Mudstone

下载: 全尺寸图片 幻灯片

5.   讨论与结论

本文提出一种高陡倾角裂缝介质的正演模拟方法． 利用裂缝等效理论将高陡倾角裂缝介质等效为横向各向同性(HTI)介质，推导出HTI介质旋转交错网格一阶速度应力方程及所对应的差分方程，实现了裂缝介质的正演模拟，分析了裂缝介质中的地震波场响应特征，讨论了不同裂缝密度、 裂缝纵横比和裂缝充填物的地震响应特征． 结果表明： 裂缝的存在相当于人为增加反射界面； 裂缝密度越大，裂缝纵横比越小，裂缝充填物与背景介质弹性性质差别越大，引起的反射波能量变化越大. 本文模拟结果为利用地震数据进行裂缝介质参数反演与储层识别及预测提供了依据． 需要说明的是： 本文结果仅仅针对一组高陡倾角裂缝介质，而实际地层中裂缝的存在形式是多种多样的，需进一步研究实际裂缝介质中的地震波场特征； 本文仅仅对裂缝充填物、 裂缝密度、 裂缝纵横比等裂缝物性参数进行了正演模拟，需要进一步研究裂缝介质其它物性参数的地震响应特征； 本文使用弹性波进行正演模拟，而实际地震记录通常是纵波记录，因此需要进一步研究各向异性介质中纵横波分离的问题．