孟凡顺，张 亮，李景岩，李洋森

（1.中国海洋大学 海洋地球科学学院，山东 青岛 266100；2.中国海洋大学 “海底科学与探测技术”教育部重点实验室，山东 青岛 266100；3.中海石油（中国）有限公司 湛江分公司，广东 湛江 524057）

准P波方程紧致交错网格井间地震波场模拟及边界条件

孟凡顺1，2，张 亮3，李景岩1，2，李洋森1，2

（1.中国海洋大学 海洋地球科学学院，山东 青岛 266100；2.中国海洋大学 “海底科学与探测技术”教育部重点实验室，山东 青岛 266100；3.中海石油（中国）有限公司 湛江分公司，广东 湛江 524057）

研究井间地震波场的形成过程以及波场的传播机理、规律，对于指导实际井间地震勘探有着重要的意义。基于具有垂直对称轴的横向各向同性（VTI）介质中的一阶准P波方程，应用具有无条件稳定性质的紧致交错网格隐式差分格式求解该方程。重点研究了紧致交错网格求解该方程的完全匹配层（PML）吸收边界条件，在此基础上实现了VTI介质中一阶准P波方程的井间地震波场的正演模拟。数值算例表明：紧致交错网格能精准模拟VTI介质中准P波的传播过程，得到高精度的正演结果。一阶准P波方程能以足够的精度描述VTI介质中准P波特征。完全匹配层吸收边界能有效地解决人工边界问题，是一种高效的边界吸收算法。

井间地震，紧致交错网格；VTI介质；准P波方程；完全匹配层

0 前言

井间地震是将震源与检波器都置于井中进行地震波观测的一种新的物探方法。该方法可以获得较高分辨率的地震信号，与地面地震互补可以大幅度地提高复杂陆相储层的描述精度［1、2］。但是当地下地质比较复杂时，尤其是地质呈现各向异性时，井间地震波场类型多，而且各种类型波之间相互干涉或叠加，形成了非常复杂的井间地震波场，这对于井间地震波场的识别和分离以及成像带来了很大的麻烦［3］。因此应用正演模拟技术研究复杂地层中井间地震波场特征和传播规律，具有非常重要的理论和实际意义。

目前，各向异性介质中地震波的正演模拟，主要通过求解弹性波方程实现［4～6］。该方法用三分量矢量场描述地震波场，每个分量都包含纵波和横波，因此能模拟比较丰富的波场信息。但也具有十分明显的缺陷［7］，主要表现在：各向异性介质弹性波方程包含多个分量，这就必然导致计算速度慢、效率低，并且这些弹性参数物理意义不明确，在实际生产中这些参数也很难得到。目前在实际生产中仍以纵波勘探为主，即使在各向异性介质中，也仍然只记录纵波信息。在目前技术的情况下，多分量理论数据的应用受到了限制。为克服以上缺陷，Alkhalifan［8］假设横波速度为零，推导出各向异性介质中的准P波方程。该方程对弹性波方程能有效的近似，更加有利于研究各向异性介质中准P波的传播规律。随后，Hestholm［9］推导了该方程的一阶速度～应力方程形式，并采用高阶差分实现了该方程的数值模拟。

在数值模拟领域中，由于有限差分法计算效率高，容易解决边界问题，编程简单等优点被广泛地应用。Madariage［10］首次提出了交错网格有限差分算法，极大地提高了有限差分数值模拟的精度。Lele［11］提出了高阶紧致网格有限差分算法，紧致网格是一种隐式差分格式，与显式差分格式不同的是其具有无条件稳定的性质，其主要优势在于应用较少的网格节点就可以得到高精度的模拟记录，而且较容易解决边界问题［12～14］。Nagarajan［15］应用紧致交错网格有限差分成功解决了大涡问题。Bendiks［16］应用紧致交错网格有限差求解了Navier－Stokes方程，结果表明紧致交错网格有限差分能够很好地压制数值频散。

作者在本文根据VTI介质中的一阶准P波方程，研究了该方程在紧致交错网格中求解的高阶有限差分算法。重点研究了紧致交错网格一阶准P波方程的PML边界方程，在此基础上实现了VTI介质中准P波方程的正演模拟，得到了高精度的合成记录。

1 基本原理

1.1 一阶准P波方程

二维VTI介质中准P波方程的一阶速度～应力方程为［9］：

其中 Vx、Vz分别为x、z方向的质点速度；ρ为密度；ψ、κ、ζ为中间变量。

1.2 紧致交错网格时间差分近似

1.3 紧致交错网格空间差分近似

设f（x）是连续函数，并且存在n阶偏导数，令

2 N点2 N＋2阶紧致交错网格有限差分格式为［18］式（4）。

其中 Δx是空间采样间隔，待定系数a、bn可以通过将处Taylor公式展开后解方程组求得。8点10阶紧致交错网格的待定系数为：a＝0.257 89、b1＝0.889 87、b2＝ 0.216 12、b3＝－0.004 7、b4＝0.000 15。例如，在x方向应力σxx的空间一阶差分8点10阶近似为：

其中 A和B是n×n阶系数矩阵，Fj和fj是n×1阶列向量。

向前差分格式：

z方向空间一阶差分8点10阶近似为：

其中 Fi＝ ［Fi，1Fi，2…Fi，n－1Fi，n］；fi＝ ［fi，1fi，2…fi，n－1fi，n］是1×n维行向量，AT是A 的转置，向前差分B＝B1，向后差分B＝B2。

1.4 紧致交错网格差分格式的构造

紧致交错网格单元的空间交错策略见下页图1，利用公式（5）、公式（6）及网格的空间分布即可得VTI介质中一阶准P波方程的紧致交错网格时间二阶，空间8点十阶有限差分格式为：

式中 Δx、Δz为空间离散间隔；Δt为时间离散间隔；i、j为空间离散点序号；n为时间离散点序号。

1.5 PML吸收边界条件

根据PML分裂原理［19］针对二维VTI介质中的一阶准P波方程，应用紧致交错网格有限差分导出了适应该方程的吸收边界方程，见式（8）～式（16）。

其中 式（12）和式（13）分别对应x方向和z方向，其余各式在x方向和z方向均可使用；Px、Pz分别为P在x方向和z方向上的分量；d（x）、d（z）为x与z方向的吸收因子。

作者以式（14）～式（16）为例给出其时间二阶精度，空间8点十阶精度的紧致交错网格有限差分格式，即：

其中 di、dj为x与z方向吸收因子的离散值。

2 数值模拟

2.1 均匀VTI介质模型

首先设计一个400×400的均匀VTI介质模型，来检验紧致交错网格有限差分数值模拟算法的可行性。合成计算的空间步长和时间步长分别为Δx＝ Δz＝10m，Δt＝1ms，介质参数v ＝2 401m／s、vv＝ 2 450m／s、η ＝ 0.518、ρ ＝2 150kg／m3。震源为主频30Hz的Ricker子波位于模型的中间，如下页图2所示。

采用时间二阶精度，空间八点十阶精度对上述模型进行模拟，得到的波场快照如图3所示（见下页）。波场快照中外层传播较快的即为准P波，内层呈菱角状传播较慢的波即为由计算误差产生的P2波［20］。

由分析波场快照可以得出：紧致交错网格隐式格式有限差分模拟的地震记录数值频散较小，模拟的效果比较有利于波场的识别、速度分析、偏移成像等。紧致交错网格隐式格式准P波方程的PML吸收边界，能有效地吸收人工边界产生的虚假反射波，而且与内层网格的高阶格式相匹配，避免了吸收边界区与内层网格交界面由于差分阶数的突变而产生的伪反射，是一种高效的边界吸收算法。

为了验证准P波方程对VTI介质中纵波的近似程度，作者在本文做了如下试验，震源点位于模型中间，接收点位于震源右下角水平和垂直距离均为1 000m处，其中弹性波方程正演模拟采用的参数和一阶准P波方程相一致，可通过计算得出［21］。并且弹性波方程模拟时采用纵波源激发，弹性参数同上。

后面图4为在相同弹性参数的情况下，分别应用一阶准P波方程和弹性波方程得到的理论地震图，分析地震图可以得到以下认识：

（1）紧致交错网格有限差分算法产生的数值频散很小，主要表现在接收的振动曲线平滑很少出现剧烈抖动。

（2）一阶准P波方程在动力学上对VTI介质中的准P波具有很高的近似精度。主要表现在一阶准P波方程和弹性波方程在准P波的起跳时间上一致，而且两种方程得到振动曲线准P波部份重合度很高。

图1 紧致交错网格空间单元交错示意图Fig.1 Elementary cells of compact staggered grids

（3）弹性波方程模拟时，即使采用纵波源激发VTI介质中也并非只产生qP波，还将产生与qP波耦合的qSV波；准P方程产生与弹性波方程同样的qP波，还将产生由计算误差引起的P2波。

2.2 层状介质模型

模型如图5所示（见下页），模型计算区域为400×400，空间采样间隔Δx＝Δz＝10m，时间采样间隔Δt＝1ms，激发井和接收井分别位于水平1 500m和2 500m处，井间距1 000m，道间距10m，共400道。记录时间长度t＝1 500ms，炮点深度为1 200m和2 000m，介质的分界面位于垂直深度1 500m、2 500m处，介质“1”为各向同性介质，介质“2、3”为VTI介质。介质的相关参数见后面表1。其它模拟参数同均匀VTI介质模型。

表1 层状介质模型参数Tab.1 The physical parameters of layer media

通过分析图6、图7井间地震层状介质模型的共炮点道集记录特征可以得出：

（1）紧致交错网格有限差分方法能高精度的模拟井间地震波场。各向同性介质中井间地震直达波同向轴为双曲线，而VTI介质中井间地震直达波的同向轴不再是双曲线形态。因为当地层表现为各向同性时，地震波传播的速度在不同方向上为恒定不变的，时距曲线方程为双曲线方程。当地层表现为VTI介质时，地震波传播的速度在不同方向上是变化的，此时方程就不是双曲线方程了，表现在共炮点道集记录上直达波就不是双曲线形态。

（2）震源位于各向同性介质中时，只产生纵波向四周传播，波传播到模型交界面处将产生透射纵波和反射纵波（见图6标注），界面处不产生转换横波。这是一阶准P波方程模拟和弹性波方程模拟井间波场时的主要区别，弹性波方程模拟时将产生各种转换横波［22］。对比结果表明，一阶准P波方程能够对VTI介质中准P波进行很好的近似。

（3）震源位于VTI介质中时，不仅产生纵波还产生实际介质中并不存在的P2波，此波向上传播到各向同性介质交界面处将产生反射波（见图7中标注I波）；当向下传播到VTI介质的交界面处将产生反射波（见图7中标注II波）和透射波（见图7中标注III波）。得知P2波只产生于VTI介质中，在各向同性介质中不会产生P2波，因此实际地质条件下，我们只需将震源放在沉积相对比较稳定可以近似认为各向同性的地层中进行激发，就可以有效地消除P2波的影响。

3 结论

应用紧致交错网格有限差分对井间地质模型进行数值模拟，结果表明：

（1）紧致交错网格能准确模拟VTI介质中准P波的传播过程，得到的井间地震波场齐全、清晰。克服了弹性波场正演模拟波场复杂、计算效率低、弹性参数不明确等缺点。

（2）完全匹配层吸收边界能够有效解决人工边界问题，并且和计算区域交界面处产生虚假反射很小，是一种理想的吸收边界条件。

（3）地震各向异性的研究，使地震学理论向实际地球介质波动理论研究迈出了一大步。随着对地震精度要求的不断提高，基于各向异性的处理解释技术也受到人们越来越多的重视。一阶准P波方程能高精度地描述VTI介质中准P波的特征，具有较高的近似精度，在纵波勘探领域有着广阔的应用前景。

（4）在井间二维VTI介质中准P波的数值模拟方面，将震源置于各向同性介质中激发可以消除计算误差产生的P2波影响，但是还需进一步研究更好的去除P2波的方法，使震源位于VTI介质时也能去除P2波的影响。

［1］ 周建宇，何惺华，李安夏，等.罗家地区井间地震方法与效果［J］.石油地球物理勘探，2001，36（6）：745.

［2］ 孔庆丰，王延光，左建军，等.井间地震有限角度叠加方法研究与应用［J］.石油地球物理勘探，2007，42（3）：256.

［3］ 李万万，裴正林.井间地震弹性波传播特征数值模拟［J］.物探与化探，2008，32（2）：207.

［4］ ZHANG J F.VERSCHUUR D J.Elastic wave propagation in heterogeneous anisotropic media using the lumped finite－element method［J］.Geophysics，2002，67（2）：625.

［5］ 祝贺君，张伟，陈晓非.二维各向异性介质中地震波场的高阶同位网格有限差分模拟［J］.地球物理学报，2009，52（6）：1536.

［6］ 汪利民，徐义贤，江贵.井中和井间地震波场正演模拟［J］.石油物探，2009，48（2）：146.

［7］ 韩令贺，何兵寿.VTI介质一阶准P波方程正演模拟及边界条件［J］.石油地球物理勘探，2010，45（6）：819.

［8］ ALKHALIFAH T.An acoustic wave equation for anisotropic media［J］.Geophysics，2000，65（4）：1239.

［9］ HESTHOLM S.Acoustic VTI modeling using high－order finite－differences［J］.Geophysics，2009，74（5）：T67.

［10］MADARIAGA R.Dynamics of an expanding circular fault［J］.BSSA，1976，66（3）：639.

［11］LELE S K.Compact finite difference schemes with spectral－like resolution［J］.J.Comput Phys，1992，103（1）：16.

［12］王书强，杨顶辉，杨宽德.弹性波方程的紧致差分方法［J］.清华大学学报：自然科学版，2002，42（8）：1128.

［13］马廷福，金涛，曹福军，等.三维双调和方程的高阶紧致差分格式及其多重网格方法［J］.兰州理工大学学报，2010，36（3）：142.

［14］葛永斌，刘国涛.三维Helmholtz方程的高阶隐式紧致差分方法［J］.工程数学学报，2010，27（5）：853.

［15］NAGARAJAN S.LELE S K，FERZIGER J H.A robust high－order compact method for large eddy simulation［J］.Comput.Phys，2003，191（2）：392.

［16］BENDIKS.A staggered compact finite difference formulation for the compressible Navier－Stokes equations［J］.Comput.Phys，2005，208（2）：675.

［17］董良国，马在田，曹景忠，等.一阶弹性波方程交错网格高阶差分解法［J］.地球物理学报，2000，43（3）：411.

［18］DU Q Z，LIN B，HOU B.Numerical modeling of seismic wavefields in transversely isotropic media with a compact staggered－grid finite difference scheme［J］.Applied Geophysical，2009，6（1）：42.

［19］COLLINO F，TSOGKA C.Application of the perfectly matched absorbing layer model to the linear elastodynamic problem in an isotropic heterogeneous media［J］.Geophysics，2001，66（1）：294.

［20］何兵寿，张会星.VTI介质中准P波方程的数值解法［J］.煤炭学报，2006，31（4）：446.

［21］TARIQ A.An acoustic wave equation for anisotropic media［J］.Geophysics，2000，65（4）：1239.

［22］吴律.层析基础及其在井间地震中的应用［M］.北京：石油工业出版社，1997.

P 631.4

A

10.3969／j.issn.1001－1749.2012.05.02

1001—1749（2012）05—0510—08

2011－12－06 改回日期：2012－05－17

孟凡顺（1960－），男，博士后，教授，主要从事地震波数值模拟与层析成像研究。