地球介质广泛存在波动各向异性，地震各向异性主要表现在地震波传播速度是传播方向的函数、体波间的相互耦合、横波发生分裂等(吴国忱，2006).裂隙介质不仅是油气运移的重要通道，更是重要的油气存储空间，因此在油气勘探领域越来越受到重视(金抒辛和何樵登，2005；吴国忱和秦海旭，2014).裂隙中弹性波传播表现为方位各向异性特征，具有方位各向异性特征的介质包括HTI介质、正交各向异性介质和单斜各向异性介质等，方位各向异性特征主要表现在振幅随炮检距和方位角的变化特征及NMO速度随炮检距和方位角的变化特征等方面(王月英等2010).地震波在裂隙介质中传播时，其振幅、速度和走时等都会随测线方位角和入射角的变化而变化.运用地震波的动力学和运动学特征，通过数值模拟方法模拟弹性波在裂隙介质中的传播过程，可以帮助找出裂隙发育方向及其他物性，进而为油气田的评价和开发服务.

1980年Schoenberg给出了线性滑动界面的概念，描述了弹性波在线性滑动界面的传播现象(Schoenberg，1980)，并于1988年给出了单组平行排列直立裂隙弹性参数的计算方法(Schoenberg and Douma J，1988).Bakulin在Schoenberg的基础上于2000年给出了单斜对称性各向异性介质弹性参数的计算方法(Bakulin et al.，2000a，b，c).单斜对称系统由多组非正交直立裂隙组成，图 1展示了单斜裂隙介质岩芯图片和单斜裂隙模式示意图.

图 1

Fig. 1

图 1 (a)单斜裂隙介质岩芯; (b) 单斜介质示意图 Fig. 1 (a) Core of monoclinic fractured media; (b) Schematic of monoclinic medium

地震波正演模拟是模拟地震波在地球介质中的传播过程，并研究地震波传播特性与地下介质参数的关系(孙成禹，2007).1987年Jean Virieux给出波动方程一阶速度-应力方程并用交错网格模拟地震波在非均匀介质中的传播过程(Virieux，1986).交错网格差分法相对规则网格差分法网格频散显著减小，精度明显提高，而且可以选取较大的空间步长，提高计算效率(董良国等，2000).正演数值模拟无法在无限大空间进行，人为设定计算区域会导致边界处产生明显边界反射现象，严重干涉正常波场信息的分析.本文采用最佳匹配层法(the perfectly matched layer，简称PML)消除边界反射，其由Berenger于1994年在研究电磁波传播理论时提出(Berenger，1994).

本文利用线性滑动理论结合Hudson理论，建立不同裂隙物性条件下的单斜各向异性岩石物理模型，利用交错网格有限差分法模拟弹性波在单斜介质中的传播，分析不同模型的波场和炮集，以及不同偏移距下的方位角道集.

线性滑动理论指出含有裂隙岩石的柔量可以表达成裂隙柔量与围岩柔量之和(Schoenberg and Douma J，1988)，公式为

2 000年Bakulin考虑两组不同任意方位的裂隙(Bakulin et al.，2000c)，公式为

这里我们默认第一组裂隙沿x₁方向排列，其柔量矩阵表达式为

而另一组裂隙法向方向与x₁方向夹角为Φ_i，可以利用Winterstein于1990年提出的邦德变换(Winterstein，1990)：

1993年Hsu和Schoenberg给出切向柔量和法向柔量的表达式为(Hsu and Schoenberg，1993)

当裂隙满足Hudson理论的假设：

(1)裂缝是定向排列并且稀疏分布于背景介质中，裂缝尺寸远小于地震波波长.

(2)裂缝是相互独立的薄扁球体，流体不能在裂缝之间流动，裂隙纵横比较小.

(3)裂缝内的气体、液体以及其他物质的体积模量和剪切模量都比背景介质小.

则可以从Hudson理论中给出Δ_T，Δ_N的表达式如下(Hudson，1981)：

(1)当裂缝中填充较小体积模量和剪切模量的固体时

(2)当裂缝为干裂缝时

(3)当裂缝中填充无黏滞流体时

根据应力应变关系：

与常规网格不同，交错网格差分算法将不同变量定义在整网格点和半网格点上，单斜介质变量定义网格如下图所示：

其不同网格位置上定义的变量如表 1所示：

图 2

Fig. 2

图 2 单斜介质变量定义网格 Fig. 2 Variables defined grids of monoclinic medium

表 1(Table 1)

表 1 单斜介质弹性参数变量的空间位置表 Table 1 Table of the position of elastic elastic parameters variable of monoclinic medium 网格点 1 2 3 4 5 6 7 弹性波场分量和弹性参数 σxx,σyy,σzz,c11,c12,c13,c16,c22,c23,c26,c33,c36 σyz,c14,c14 σxz,c35,c55 σxy,c66 vx,ρ-1 vy,ρ-1 vz,ρ-1,c45

表 1 单斜介质弹性参数变量的空间位置表 Table 1 Table of the position of elastic elastic parameters variable of monoclinic medium

利用泰勒公式将上述一阶速度-应力方程在时间和空间方向分别展开，即可得其差分形式，完整差分形式见附录A.

如上文所述，采用有限差分进行地震波数值模拟过程中，模拟区域不可能做成无限大，人为限定模型计算区域的方法会引起边界问题，其常常比物理边界引起的反射波更强，在地震记录中会干扰有效信息，因此针对边界问题必须采用一定的边界条件消除这种边界反射.Berenger针对电磁波传播情况，给出了一种高效的完全匹配层(PML)吸收边界条件，并在理论上证明该方法可以完全吸收来自各个方向、各种频率的电磁波，而不产生任何反射(Berenger，1994).其主要原理是在计算区域外加一定厚度的吸收层，对传入吸收层的波动造成指数形式衰减从而消除边界反射.

将波动方程分裂推导其PML形式，这里只以一个分量为例，公式为

利用前文线性滑动理论结合Hudson理论求取弹性系数矩阵建立单斜介质的均匀模型，模型大小均为2000 m×2000 m×2000 m，背景介质砂岩和泥岩各占一半，干燥裂隙，裂隙纵横比0.01，裂隙密度0.04，裂隙夹角45°.在模型中心处激发震源，时间网格选取dt=1 ms，空间网格选取dx=dy=dz=10 m.过震源的三个相互正交平面截取三分量0.2 s时的波场快照如图 4所示.

图 3

Fig. 3

图 3 三维PML吸收层示意图 Fig. 3 Schematic of three-dimensional PML absorbing layer

图 4

Fig. 4

图 4 裂隙密度0.06三维单斜介质波场 Fig. 4 Three-dimensional monoclinic medium wave field when fracture density is 0.06

结合图 1所示，单斜介质由两组非正交的垂直裂隙组成，这里给定裂隙之 间夹角为45°，xoy面内波场 呈现椭圆形，这是由于垂直裂隙方向波场传播速度小于平行裂隙方向波场传播速度，椭圆对称轴向右倾斜约22.5°，其余两个面内波场为正椭圆形.

图 5为xoy面波场随裂隙密度变化情况，裂隙密度分别给定为0.02、0.04、0.06、0.08和0.1.随着裂隙密度增大，纵波波场椭圆长短半轴之比增大，说明波场垂直裂隙方向和平行裂隙方向速度差增大，横波畸变现象加重.图 6为U_xoy、V_xoy、W_xoy三个面波场随裂隙夹角变化情况，裂隙密度分别给定为0°、30°、45°、60°和90°.随着裂隙夹角增大，纵波椭圆对称轴向右倾斜明显，当夹角为0°时波场形状类似于HTI介质，当夹角为90°时波场形状类似于正交介质.

图 5

Fig. 5

图 5 单斜介质波场随裂隙密度变化 Fig. 5 Monoclinic medium wave field changing with fracture density

图 6

Fig. 6

图 6 单斜介质波场随裂隙夹角变化 Fig. 6 Monoclinic medium wave field changing with fracture included angle

图 7为一个块状裂隙介质模型，棕色代表各向同性介质，其右下角黄色代表的单斜裂隙介质.图 8为斜层裂隙介质模型，右下方为黄色代表的单斜裂隙介质，左上方为棕色代表的各向同性模型.图 9和图 10为这两种模型的U分量三个面的波场，震源均位于顶层中心位置.我们可以清楚看出各个面波场中的透射qPqP波(P波入射后产生入射波为P波，以此类推)、 透射qPqS波、 反射qPqP波、反射qSqP波和反射qSqS波等，以及由模型中绕射点产生的绕射现象.

图 7

Fig. 7

图 7 单斜介质块状模型图 Fig. 7 Block model of monoclinic medium

图 8

Fig. 8

图 8 单斜介质斜层模型 Fig. 8 Oblique layer model of monoclinic medium

图 9

Fig. 9

图 9 单斜介质块状模型U分量波场 Fig. 9 U component of the wave field of the monoclinic medium block model

图 10

Fig. 10

图 10 单斜介质斜层模型U分量 Fig. 10 U component of the wave field of the monoclinic medium oblique layer model

如图 11所示，这里设计一个两层裂隙介质模型，大小为1000 m×1000 m×1000 m，上层黄色为单斜裂隙介质，裂隙密度0.06，裂隙夹角45°，下层棕色为各向同性介质.按照图 12所示方式抽取方位角道集，测线方向与x轴方向的夹角分别选取为0°、30°、45°、60°、90°、120°、135°、150°和180°.

图 11

Fig. 11

图 11 单斜介质两层模型图 Fig. 11 Double layers model of monoclinic medium

图 12

Fig. 12

图 12 方位角道集抽取示意图 Fig. 12 Schematic of drawing azimuth gathers

图 13为两层模型U分量不同方位测线记录的共炮点道集，从中我们可以清晰分辨直达qP波、直达qS波、反射qPqP波、反射qSqS波和反射qPqS波等的同相轴.这里设定两组裂隙的发育方向与x轴夹角为45°和90°，红圈内(测线方向小于90°)反射波旅行时大于黄圈内(测线方向大于90°)，且相比黄圈内反射同相轴振幅较强，这都是由两组裂隙发育方向决定的，体现了单斜裂隙介质的方位特性.

图 13

Fig. 13

图 13 不同方位CSP道集 Fig. 13 CSP gathers of different azimuths

若每隔5°抽取固定偏移距的单道道集并从左至右平行排列可以得到三个分量的方位角道集，如图 14所示，这里给定偏移距为800 m.从图中我们可以分辨直达qP波、直达qS波、反射qPqP波、反射qSqS波和反射qPqS波等，不同方位旅行时成周期变化，椭圆指出旅行时较小的方位.

图 14

Fig. 14

图 14 三分量方位角道集 Fig. 14 Azimuthsgathers of three components

单斜介质由两组非正交的垂直裂隙等效而成，其相对于单组平行排列的垂直裂隙(等效为HTI介质)和两组正交垂直裂隙(等效为正交介质)更具一般性.本文利用线性滑动理论结合Hudson理论根据裂隙物性求取裂隙介质弹性矩阵的方法建立单斜裂隙介质多种模型，根据裂隙介质弹性动力学理论建立三维单斜介质一阶速度-应力方程，利用交错网格高阶有限差分法模拟地震波在裂隙单斜介质中的传播过程.比较不同裂隙物性均匀模型中的波场，分析裂隙物性对单斜介质地震波传播的影响.通过简单模型说明地震波在单斜介质中的传播过程及透射反射现象.最后通过层状介质分析不同方位测线记录道集的特点和差别，从而说明单斜裂隙介质的方位特性.

A.三维单斜介质一阶速度-应力方程差分形式

B.三维单斜介质一阶速度-应力方程波场分裂形式和PML边界条件波动方程