含有倾斜薄裂缝孔隙地层中的井孔声场

阎守国，谢馥励，龚丹，章成广，张碧星

1中国科学院声学研究所声场声信息国家重点实验室，北京 100190

2长江大学油气资源与勘探技术教育部重点实验室，荆州 434102

1 引言

随着世界油气勘探的发展，裂缝性储集层已成为油气勘探的主要对象之一.在声波测井中，由于地层中的裂缝宽度一般很小，裂缝与井壁之间的位置关系也不确定，且往往是多条裂缝一起构成裂缝带，因此裂缝地层对阵列声波测井波形的影响十分复杂，这使得对裂缝地层的识别与评价成为当今测井所面临的难题之一.通过数值模拟计算，对含裂缝地层的井孔声场进行研究，掌握裂缝对井内声场的影响规律，对于声波测井数据处理和解释具有重要意义.

针对裂缝性地层的井孔声场问题，早期的研究工作采用平板状裂缝模型（Hornby et al.，1989），利用低频解析公式分析计算了斯通利波通过单一裂缝时的传播及反射特性，之后类似的研究（Tang，1990；Tang and Cheng，1993）指出，声波在经过渗透带处也会产生反射的斯通利波，后来也有学者通过不同的方法（Derov et al.，2009；Bakku et al.，2013）研究了裂缝存在时的井孔声场，但这些方法大都假定声场具有轴对称性；积分方程法（Spring and Dudley，1992）、有限差分法（Kostek et al.，1998；陈德华等，2004；丛健生，2004）和有限元方法（Matuszyk，2013）陆续被应用于裂缝性地层井孔声波场的研究中，但他们的计算都针对二维情况，并且水平层或裂缝宽度较大.由于有限元方法在地震波模拟问题上比有限差分法更耗时，因此对于井孔声场的模拟大都采用有限差分方法，三维交错网格有限差分法近年来被应用于大斜度井、倾斜分层地层及各向异性介质地层等非轴对称问题的数值模拟研究（Leslie and Randall，1992；Cheng et al.，1995；Sinha，2006；林伟军等，2006；阎守国等，2011），得到了较好的效果，但采用该方法对非轴对称裂缝地层的井孔声场研究还未见报道.这主要是由于当裂缝宽度较小时，需要划分细小的网格来满足薄裂缝的计算需求，而这将大大增加三维有限差分的计算量，导致模拟计算无法顺利进行.本文针对这一问题，对含有倾斜薄裂缝孔隙地层中的井孔声波场采用三维不规则交错网格有限差分法进行数值模拟计算.考察不同裂缝及地层参数情况下的井孔声场特性，为进一步研究利用井孔声场信息反演裂缝参数奠定基础.

2 物理模型及计算方法

2.1 物理模型

本文采用三维应力－速度交错网格有限差分法数值模拟计算含有倾斜裂缝地层的井孔声场问题.图1为含有倾斜裂缝地层的井孔模型示意图.假定声源位于笛卡尔坐标系的原点，井轴与z轴重合，井外裂缝与xz轴构成的平面垂直，与xy轴构成的平面倾角为α，裂缝下界面与z轴的交点与原点的距离为dz，裂缝的垂直厚度为H，裂缝间隔为h.

图1 倾斜裂缝地层中的井孔模型Fig.1 The configuration of borehole in porous formation with tilted thin fracture

2.2 变网格有限差分处理

由于需要计算的裂缝很薄，并受井孔模型大小的影响，如果采用很小的网格进行数值计算会导致计算速度过慢，甚至出现无法计算的情况.所以我们采用不均匀网格的办法在裂缝处采用小网格，在非裂缝处采用大网格，这样可以大大的减小计算量，在保证计算精度的同时提高计算速度.在使用变网格时，网格步长的变化可能会导致数值反射的出现.其原因是因为波场离散后，相速度是网格步长的函数，当相速度梯度较大时，即使速度和密度都没有变化，入射波的能量也会部分反射回来，导致数值反射现象.在网格变化的区域对波场进行插值计算虽可以在某种程度上压制数值反射，但是插值算法计算量大，且效果并不令人满意，所以本文没有采用插值计算.本文采用的算法是根据各点所对应的步长的变化计算其差分系数，此方法不会产生数值反射，并且减小了计算量（孙卫涛和杨慧珠，2004）.下面给出了具有四阶空间精度的不规则网格差分算子（图2）.

图2 变网格差分示意图Fig.2 The configuration of nonuniform finite－difference grid

节点i处的场量ψ对空间的一阶导数∂ψ／∂x可以表示为相邻节点的线性组合形式：

其中，

（1a）即为网格发生变化时的四阶空间精度差分算子.而当s1＝r1＝r2＝Δr时，即在均匀网格情况下，（1a）式可以退化为均匀网格四阶空间精度差分算子，此时

2.3 流体－孔隙介质界面处的差分方程统一

由于实际地层中往往是裂缝与孔隙同时存在，为了能够更加接近实际情况，同时也为了考察孔隙度渗透率等地层参数对裂缝产生的影响，本文在数值模拟时采用流体饱和孔隙介质作为井外地层介质，地层中的弹性波传播采用Biot孔隙弹性波方程组描述（Biot，1955；1956；1962a；1962b；王秀明等，2003）.孔隙介质中场量和参数更多，波动方程形式复杂，与流体波动方程形式不统一，不能直接应用参数平均法获得流体－孔隙介质交界面上场量的离散表达式，因此无法直接进行差分计算.关威等（Guan et al.，2009）曾利用井壁边界上的应力和位移连续性条件推导出了井壁网格点的差分离散表达式，实现了对孔隙地层中声波测井的有限差分模拟算法，但这种处理方法在算法的实现上相对麻烦，尤其是地层中存在裂缝时，需要利用界面处的边界条件重新推导差分公式.本文采用将裂缝中的流体利用孔隙介质方程参数取流体极限的办法，统一了弹性波动方程，简化了介质交界面处的处理过程.

依据Biot理论，孔隙介质的本构关系为

其中，Tij为总应力张量，S为流体应力，Pf为流体压强，分别为颗粒、骨架以及流体的体积模量，u为固体位移矢量，w为固体相对流体位移矢量，φ为孔隙度，μb是骨架剪切模量，λc是拉梅系数，λc＝H－2μb，H＝λb＋2μb＋α2M，M＝

将方程（2）两边对时间t求导，并进行分量展开，可写为

其中，Vx，Vy，Vz分别表示固体介质在x，y，z方向上的质点速度分量，Wx，Wy，Wz分别表示固体相对流体介质在x，y，z方向上的质点速度分量.

孔隙介质的运动方程为

式中，σij＝Tij＋φPfδij，ρ11＝ （1－φ）ρs－ （1－E）φρf，ρ12＝ （1－E）φρf，ρ22＝Eφρf，ρs，ρf分别为固体颗粒与孔隙流体密度，E为孔隙迂曲度，可近似表示为，K为渗透率，η为黏滞系数.

如果令：Q1＝－ρ11ρ22，R1＝ρ12／Q1，R2＝ρ22／Q1，R3＝ （ρ12＋ρ22）／Q1，P1＝ （ρ11＋ρ12）／Q1，P2＝ （ρ12＋ρ22）／Q1，P3＝ （ρ11＋2ρ12＋ρ22）／Q1，则孔隙介质中的速度分量可以写为

当场量和介质参数取极限的状态下，描述孔隙介质中弹性波传播的Biot方程组可以退化为单相的固体或流体介质中的弹性波方程组.本文仅考虑退化为流体的情况，当孔隙介质参数选取如下极限状态时，退化为流体介质：φ＝1，Kb＝0，Ks＝Kf，H＝M＝Kf，μ＝0，a＝1，b＝0，由于此时ρ11＝0，ρ12＝0，ρ22＝ρf，使得Q＝0，P，R的计算不再有意义，可根据流体中的波动方程直接定义R1＝1／ρf，R2＝0，R3＝0，流体中不存在相对渗流位移，可直接令W＝0.在流体与孔隙介质交界面处进行参数平均处理，密度取算数平均值，弹性系数取调和平均值.经过以上的处理，在计算中遇到流体和孔隙介质交界面时可以直接进行差分计算，非常适合处理含有裂缝的地层模型.

2.4 差分算法的实现

本文采用三维应力－速度交错网格有限差分，各物理量在网格上的离散规则如下：

其中，下标i，j，k表示为整网格点，下标i＋1／2，j＋1／2，k＋1／2为半网格点，声源类型为点源，时域脉冲选择变形的Ricker子波，表达式为f（t）＝，时域波形及频谱如图3所示.

图3 声源时域波形及频谱Fig.3 The waveform in time domain and spectrum in frequency domain of source pulse

3 计算结果及分析

为符合裂缝性储层的实际特点，本文选择低孔隙度、低渗透率地层进行数值模拟计算.所采用的介质参数如表1所示，其中孔隙介质参数Ks，ρs为孔隙固体颗粒的体积模量及密度，VrP，VrS为骨架纵、横波速度，η，φ，K分别为黏滞系数、孔隙度及渗透率，Vf，ρf为孔隙内流体的速度及密度.井孔半径a＝0.1m，计算过程中当声源频率不大于5kHz时，空间网格最大可选为2cm，变网格差分中可任意更改小网格大小，根据裂缝宽度不同，本文计算过程中最小网格选取为20μm.

表1 模型的介质参数Table 1 The parameters of the media

3.1 算法正确性验证

当井外为无裂缝的孔隙介质地层时，将有限差分结果与实轴积分法结果进行对比，对比结果如图4所示.声源频率为2.5kHz，接收器与声源间距为1m.图中黑线为有限差分法计算结果，红线为实轴积分法计算结果.两种方法计算结果基本吻合，证明了有限差分方法的正确性.图中幅度上的差异是由归一化造成的，不影响本文结果.

图4 有限差分与实轴积分法的结果对比图Fig.4 The comparison results of real axis integrate method and finite－difference method

采用完全匹配层（PML）来处理各方向吸收边界（张鲁新等，2010）.图5为井外为均匀孔隙介质情况下，xz平面上z方向应力（Tzz）的声场快照，图中所示时刻为1.35ms，x，y，z方向上的计算距离均为±2.0m，从图5可以看出各方向均无反射波存在，边界吸收效果良好.

图5 井外无裂缝时xz平面上Tzz的声场快照图Fig.5 The Tzzacoustical fields snapshot of xz plane when there is no fracture outside borehole

3.2 裂缝对井内声场的影响

裂缝引起的介质突变会使声波在界面处产生反射，并使穿过裂缝的声波能量减弱.图6是井外存在裂缝情况下xz平面上场量Tzz的快照图，计算范围为±2.0m，水平裂缝位于声源上方0.9m处，裂缝宽度为2mm.图中P，S，ST分别代表纵波、横波和斯通利波，下标d代表直达波，r代表由裂缝产生的反射波.由图6a可见，此时阵列波形图中可看到直达的纵波、横波及斯通利波，斯通利波只在井孔内沿井轴方向传播.纵波通过裂缝后直达场减弱，此时从图中还观察不到反射波；声波经过一段时间的传播，在图6b中可以看到反射纵波，并且由于孔隙介质的耗散作用，井内的斯通利波在传播过程中强度逐渐减弱；图6c中可以看到反射横波的出现，此时由于各波列之间的相互干扰，以及头波的影响，井外声场显得比较复杂.虽然从阵列波形图中可以明显地看到反射波的存在，裂缝的位置也可以很直观的反映出来，但是由于实际测井中只能记录井内的声场，因此需要针对阵列波形图进行研究.

图6 不同时刻，井外存在裂缝时xz平面上的声场快照图（Tzz）（a）0.51ms；（b）0.68ms；（c）0.85ms.Fig.6 The Tzzacoustical fields snapshot of xz plane when there are fractures outside borehole

图7为井外存在裂缝时井轴上的阵列波形图，横轴是接收时间，纵轴为接收点距离声源的距离，波列的幅度经过归一化处理.从图7中可以看到，井轴上的阵列波形图存在纵波（P）、横波（S）以及斯通利波（ST），由于声源激发频率及地层参数的影响，此时纵、横波激发幅度较小，斯通利波占主要成分，当正常传播的各种波型遇到裂缝时，就会产生相应的反射波，图7中可以观察到的反射波列有两条，按照其传播速度判断分别为横波及斯通利波反射.

图7 井轴上的阵列波形图Fig.7 The array waveform neceived on the borehole axis

如果能够在阵列波形图中观察到反射波列，即可发现裂缝的存在.但是由于全波列声波测井仪器的声源频率一般在15kHz以内，比较常用的频率在2.5～10kHz之间，而地层中裂缝的宽度一般都比较小，介于20～2000μm之间，裂缝产生的反射波强度往往很弱，因此在声波测井的频段范围内能否观察到细小裂缝的存在是我们首先要考察的问题.

图8是数值模拟计算得到的声源频率为2.5kHz时不同裂缝宽度情况下，井轴上接收到的阵列全波图.此时全波图的主要成分是斯通利波，纵、横波的激发幅度很小，从波列图中可以得到斯通利波的速度约为1459m·s－1.图中a、b分别是裂缝宽度为2cm和2mm的情况，从图8a中可以观察到裂缝反射的斯通利波，但随着裂缝宽度减小，在图8b中已经无法观察到反射波的存在.可以看出低频声源对单条裂缝的检测能力较弱.

图8 数值模拟计算得到的声源频率为2.5kHz时不同裂缝宽度情况下，井轴上接收到的阵列全波图（a）H＝2cm；（b）H＝2mm.Fig.8 The simulation results of full waveform array received on the borehole axis with different fracture width，the source frequency is 2.5kHz

图9 数值模拟计算得到的声源频率为5.0kHz时不同裂缝宽度情况下，井轴上接收到的阵列全波图（a）H＝2cm；（b）H＝2mm；（c）H＝0.2mm；（d）H＝0.02mm.Fig.9 The simulation results of full waveform array received on the borehole axis with different fracture width，the source frequency is 5.0kHz

图9是数值模拟计算得到的声源频率为5.0kHz时不同裂缝宽度情况下，井轴上接收到的阵列全波图.与2.5kHz时的情况相同，此时波列中占主要成分的依然是斯通利波，其速度约为1490m·s－1.从图中可以看到，当裂缝宽度较大时（图9a），存在明显的两条反射波列，按照其传播速度判断分别为地层中沿井壁滑行的横波及井内流体中斯通利波，当裂缝宽度减小后，反射波列中的斯通利波明显减弱，在裂缝宽度为2mm时（图9b）就已经很难观察到反射斯通利波波列，但反射波列中的横波却并没有减弱，并且当裂缝宽度减小到0.02mm时依然可以观察到明显的反射横波.这一现象有助于我们探测水平薄裂缝的存在.

实际储层中的裂缝往往不是单一出现的，而是由多条薄裂缝一起构成一条裂缝带，我们定义裂缝带密度为d＝H／（H＋h），它相当于裂缝在整个裂缝带范围内的占空比.图10给出了由单条宽度为2mm的裂缝构成的裂缝带在不同的裂缝带密度情况下的阵列波形图，图中裂缝带宽度均为6cm，图10（a，b）是频率为2.5kHz，裂缝带密度分别为0.1和0.25时的波形图，图10（c，d）是频率为5.0kHz时的对应情况.图8b及图9b是单条2mm的裂缝在这两个频率下的阵列波形图，通过对比可以看出声波通过裂缝带时会产生相比于单条裂缝更强的反射斯通利波，且随着裂缝带密度的增加反射波也随之增强，最终会形成相当于同等宽度的宽裂缝的反射效果.同时，对比2.5kHz与5.0kHz时的结果，发现高频声源在裂缝带密度相对较低时的反射波列更加明显.因此，在声波测井常用的频率范围内，高频声源对低密度裂缝带的检测能力更强，更适合裂缝带的检测.

在实际地层中裂缝与井孔的关系往往很复杂，可能出现与井孔相交或者不相交等不同情况，经过模拟计算发现，在本文所使用的声源类型及频率范围内，未与井孔相交的裂缝或裂缝带对井内声场产生的影响很难通过阵列波形图加以辨别，因此本文只考虑裂缝或裂缝带与井孔相交存在一定倾斜角度的情况.图11是声源频率为5.0kHz时，裂缝带相对井孔的倾斜角度分别为15°、30°和45°时的阵列波形图，裂缝带由单条宽度为2mm的裂缝构成，裂缝带密度为0.25，总宽度为6cm.通过与相同状态下水平裂缝的对比（图10d），可以发现此时裂缝带的反射横波消失了，即使倾斜角度较小时（图11a）也无法观察到；但裂缝倾斜角度对反射斯通利波的影响并不大，图中用符号ST标注出了裂缝产生的反射斯通利波波列，对比图11（a—c）可发现，即使在裂缝倾斜角度较大时依然可以观察到明显的反射斯通利波波列.

图10 单条宽度为2mm的裂缝构成的裂缝带在不同的裂缝带密度情况下的阵列波形图，裂缝带宽度均为6cm（a）f＝2.5kHz，d＝0.1；（b）f＝2.5kHz，d＝0.25；（c）f＝5.0kHz，d＝0.1；（d）f＝5.0kHz，d＝0.25.Fig.10 The simulation results of full waveform array with different fracture zone density，the fracture zone constituted by the fractures 2mm width，and the fracture zone width is always 6cm.

图11 不同裂缝倾斜角度时的阵列波形图（5.0kHz）（a）α＝15°；（b）α＝30°；（c）α＝45°.Fig.11 The simulation results of full waveform array with different tilted angle of the fracture zone and borehole axis，the source frequency is 5.0kHz.

图12 不同渗透率条件下的阵列波形图（a）K＝0.002μm2；（b）K＝0.005μm2.Fig.12 The simulation results of full waveform array with different permeability of formation

图12显示了孔隙介质渗透率改变时阵列波形图的变化情况，声源激发频率为5kHz，水平裂缝，图中标出了各波列的位置，P，S，ST分别对应纵波、横波及斯通利波.从图12可以看到，从声源发出的声波在未遇到裂缝时有纵波、横波及斯通利波三种成分，由于此时纵、横波激发强度较弱，很难从阵列波形图中观察到它们的反射波，而斯通利波幅度相对较大，当渗透率相对较小时，斯通利波衰减较弱，当其受到界面的影响时，产生较强的反射横波及斯通利波；随着渗透率的增大，斯通利波的衰减明显增强，使得斯通利波在到达裂缝时能量减弱，此时反射波还可以观察到反射的横波但反射的斯通利波则明显减弱，很难从波形图中观察到.

4 结语

本文采用有限差分法数值模拟了有裂缝及裂缝带存在时孔隙介质地层中的井孔声场问题，采用非均匀网格有限差分，所计算的最低裂缝宽度达到了20μm，利用统一形式的声波波动方程，通过将孔隙介质退化为流体的办法，简化了流体与孔隙介质交界面处的处理过程.通过对裂缝及裂缝带宽度、裂缝倾斜角度及孔隙介质渗透率等因素的考察总结出以下结论：

（1）当存在水平裂缝时，使用低频声源（2.5kHz）只能观察到反射斯通利波，而使用高频声源（5kHz以上）还会观察到反射横波波列，随着裂缝宽度的减小，反射斯通利波逐渐减弱，但反射横波强度变化不大，可以用来检测单条水平薄裂缝的存在.

（2）当多条裂缝构成裂缝带时，裂缝带的密度和宽度会对反射波的强度产生较大影响，当密度较小时，裂缝带无法形成整体作用，反射波较弱，但随着裂缝带密度的增大，反射波逐渐增强，最终会形成相当于同等宽度的宽裂缝的反射效果.

（3）与井孔相交的裂缝随着倾斜角度的增加，其反射横波会消失，但对斯通利波的影响不大，倾斜角度很大时依然存在反射斯通利波，这对倾斜裂缝的检测是有利的.

（4）斯通利波在传播过程中会向孔隙介质辐射能量，产生衰减，随着渗透率的增大，斯通利波的衰减也随之增大，会导致界面处斯通利波幅度减小，可能观察不到界面的反射斯通利波.

本文通过数值模拟的办法，定性研究了孔隙地层井孔中声波通过裂缝及裂缝带时的反射规律，更加深入的工作是找到声波的衰减及反射透射规律与裂缝参数之间更深层次的定量关系，最终实现利用阵列声波信息对地层中裂缝相关参数的反演，这是我们下一步将要进行的主要工作.

Bakku S K，Fehler M，Burns D.2013.Fracture compliance estimation using borehole tube waves.Geophysics，78（4）：D249－D260.

Biot M A.1955.Theory of elasticity and consolidation for a porous anisotropic solid.J.Apply.Physica，26（2），doi：10.1063／1.1721956.

Biot M A.1956.Theory of deformation of a porous viscoelastic anisotropic solid.JournalofAppliedPhysics，27（5），doi：10.1063／1.1722402.

Biot M A.1962a.Mechanics of deformation and acoustic propagation in porous media.J.Apply.Physica，33（4），doi：10.1063／1.1728759.

Biot M A.1962b.Generalized theory of acoustic propagation in porous dissipative media.TheJournaloftheAcoustical SocietyofAmerica，34（9）：1254－1264.

Chen D H，Cong J S，Xu D L，et al.2004.The simulation of borehole field in fractural formation.JournalofDaqing PetroleumInstitute，28（3）：4－6，13，116.

Cheng N Y，Cheng C H，Toksöz M N.1995.Borehole wave propagation in three dimensions.J.Acoust.Soc.Am，97（6）：3483－3493.

Cong J S.2004.Simulated for the acoustic field of the layer and fractures in and outside of a borehole using finite difference method［Master thesis］.Daqing：Daqing Petroleum Institute.Derov A，Maximov G，Lazarkov M，et al.2009.Characterizing hydraulic fractures using slow waves in the fracture and tube waves in the borehole.2009SEG Annual Meeting，4115－4119

Guan W，Hu H S，He X.2009.Finite－difference modeling of the monopole acoustic logs in a horizontally stratified porous formation.J.Acoust.Soc.Am，125（4）：1942－1950.

Hornby B E，Johnson D L，Winkler K W，et al.1989.Fracture evaluation using reflected Stoneley－wave arrivals.Geophysics，54（10）：1274－1288.

Kostek S，Johnson D L，Randall C J.1998.The interaction of tube waves with borehole fractures.PartⅠ：Numerical models.Geophysics，63（3）：800－808.

Leslie H D，Randall C J.1992.Multipole sources in boreholes penetrating anisotropic formations：Numerical and experimental results.J.Acoust.Soc.Am.，91（1）：12－17.

Lin W J，Wang X M，Zhang H L.2006.Acoustic wave propagation in a borehole penetrating an inclined layered formation.Chinese J.Geophys.（in Chinese），49（1）：284－294，doi：10.3321／j.issn：0001－5733.2006.01.036.

Matuszyk P J，Torres－Verdín C，Pardo D.2013.Frequency－domain finite－element simulations of 2Dsonic wireline borehole measurements acquired in fractured and thinly bedded formations.Geophysics，78（4）：D193－D207.

Sinha B K，Şimşek E，Liu Q H.2006.Elastic－wave propagation in deviated wells in anisotropic formations.Geophysics，71（6）：D191－D202.

Spring C T，Dudley D G.1992.Acoustic－wave propagation in a cylindrical borehole with fractures.J.Acoust.Soc.Am.，91（2）：658－669.

Sun W T，Yang H Z.2004.A3－D finite difference method using irregular grids for elastic wave propagation in anisotropic media.ChineseJ.Geophys.（in Chinese），47（2）：332－337，doi：10.3321／j.issn：0001－5733.2004.02.022.

Tang X M.1990.Acoustic logging in fractured and porous formations［Ph.D.thesis］.Cambridge：Massachusetts Institute of Technology.

Tang X M，Cheng C H.1993.Borehole stoneley wave propagation across permeable structures 1.GeophysicalProspecting，41（2）：165－187.

Wang X M，Zhang H L，Wang D.2003.Modelling of seismic wave propagation in heterogeneous poroelastic media using a highorder staggered finite－difference method.ChineseJ.Geophys.（in Chinese），46（6）：842－849.

Yan S G，Song R L，LüW G，et al.2011.Numerical simulation of acoustic fields excited by cross－dipole source in deviated wells in transversely isotropic formation.ChineseJ.Geophys.（in Chinese），54（9）：2412－2418，doi：10.3969／j.issn.0001－5733.2011.09.025.

Zhang L X，Fu L Y，Pei Z L.2010.Finite difference modeling of Biot′s Poroelastic equations with unsplit convolutional PML and rotated staggered grid.ChineseJ.Geophys.（in Chinese），53（10）：2470－2483，doi：10.3969／j.issn.0001－5733.2010.10.021.

附中文参考文献

陈德华，丛健生，徐德龙等.2004.裂缝性地层中的井孔声场模拟.大庆石油学院学报，28（3）：4－6，13，116.

丛健生.2004.利用有限差分法模拟计算具有分层和裂缝地层井内外声场［硕士论文］.大庆：大庆石油学院.

林伟军，王秀明，张海澜.2006.倾斜地层中的井孔声场研究.地球物理学报，49（1）：284－294，doi：10.3321／j.issn：0001－5733.2006.01.036.

孙卫涛，杨慧珠.2004.各向异性介质弹性波传播的三维不规则网格有限差分方法.地球物理学报，47（2）：332－337，doi：10.3321／j.issn：0001－5733.2004.02.022.

王秀明，张海澜，王东.2003.利用高阶交错网格有限差分法模拟地震波在非均匀孔隙介质中的传播.地球物理学报，46（6）：842－849.

阎守国，宋若龙，吕伟国等.2011.横向各向同性地层斜井中正交偶极子激发声场的数值模拟.地球物理学报，54（9）：2412－2418，doi：10.3969／j.issn.0001－5733.2011.09.025.

张鲁新，符力耘，裴正林.2010.不分裂卷积完全匹配层与旋转交错网格有限差分在孔隙弹性介质模拟中的应用.地球物理学报，53（10）：2470－2483，doi：10.3969／j.issn.0001－5733.2010.10.021.