地球物理学报  2015, Vol. 58 Issue (4): 1290-1304   PDF    
引用本文
杜启振, 郭成锋, 公绪飞. 2015. VTI介质纯P波混合法正演模拟及稳定性分析. 地球物理学报, 58(4): 1290-1304, doi: 10.6038/cjg20150417   
DU Qi-Zhen, GUO Cheng-Feng, GONG Xu-Fei. 2015. Hybrid PS/FD numerical simulation and stability analysis of pure P-wave propagation in VTI media. Chinese J. Geophysics, 58(4): 1290-1304, doi: 10.6038/cjg20150417   
VTI介质纯P波混合法正演模拟及稳定性分析
杜启振, 郭成锋, 公绪飞    
中国石油大学(华东)地球科学与技术学院, 青岛 266580
投稿时间:2014-09-02    最后修改时间: 2014-12-16 .
基金项目：国家自然科学基金(41174100)、国家科技重大专项(2011ZX05019-008-08)和中国石油天然气集团公司(2014A-3609)联合资助.
作者简介：杜启振, 教授,主要从事地震弹性波传播理论与成像. E-mail:multi-wave@163.com
摘要：各向异性介质纯P波方程完全不受横波的干扰,在一定程度上可以减缓由于介质各向异性引起的数值不稳定,本文推导了具有垂直对称轴的横向各向同性(VTI)介质纯P波一阶速度-应力方程.由于纯P波方程存在一个分数形式的伪微分算子,无法直接采用有限差分法求解.针对该问题,本文采用伪谱法和高阶有限差分法联合求解波动方程,重点分析了混合法求解纯P波一阶速度-应力方程的稳定性问题,并给出了混合法求解纯P波方程的稳定性条件.数值模拟结果表明纯P波方程伪谱法和高阶有限差分混合法能够进行复杂介质的正演模拟,在强变速度、变密度的地球介质中仍然具有较好的稳定性.
关键词： VTI介质     纯P波方程     稳定性条件     有限差分法     伪谱法    
Hybrid PS/FD numerical simulation and stability analysis of pure P-wave propagation in VTI media
DU Qi-Zhen, GUO Cheng-Feng, GONG Xu-Fei    
School of Geosciences, China University of Petroleum (East China), Qingdao 266580, China
Abstract: Seismic anisotropy is undergoing rapid expansion as a result of recent advances in parameter inversion and seismic processing. The simplest and most commonly used anisotropic model is the transverse isotropy (TI) with a vertical symmetry axis (VTI) or with a tilted symmetry axis (TTI). TI models can provide good representation of intrinsic anisotropy of the earth media. Vertical and tilted transverse isotropy have become an integral part of velocity fields in prestack depth migration algorithms for imaging complicated anisotropic structures, especially those algorithms based on the wave equation such as reverse time migration (RTM). Reverse time migration in anisotropic media based on the full elastic anisotropic wave equation demands tremendous computation and requires an anisotropic model that simultaneously contains P- and S-wave velocity information. However, the seismic data recorded with conventional vertical geophones often mainly contain P-wave reflections. So an alternative way to reduce parameter requirement and computational cost for P-wave modeling and reverse time migration in TI media is to apply an acoustic approximation. Acoustic approximation yields a kinematically accurate P-wave propagation for TI media, however, this method can't completely get rid of shear-wave noise. Moreover, acoustic approximation equations become unstable for anisotropic models with parameters ε<δ and for heterogeneous models with highly varying dip and azimuth angles. For these above mentioned reasons, the pure P-wave equation for modeling and migration in TI media has attracted more and more attention. The pure P-wave equations are also called decoupled equations which are derived from the exact dispersion by using the first-order Taylor series expansions. Here, we propose the first-order form of the pure P-wave equation with stress and particle velocities as wavefield variables for VTI media. Our main work is restricted to pure P-wave mode for practicality. Just like most other anisotropic pure wave mode equations, the proposed first-order pure P-wave equation also involves pseudo-differential operators in space which are difficult to handle by using finite-difference (FD) method directly. The most natural choice to implement these equations is pseudospectral (PS) method. Although pseudospectral method requires only two nodes per wavelengths, it is still computationally expensive due to several forward-backward Fourier transforms in wavefield extrapolation at each time step. Here we adopt a hybrid pseudospectral and finite-difference scheme to solve the first-order pure P-wave equation. The hybrid PS/FD scheme requires only two Fourier transforms at each time step in 2D case and thus it is faster than the pseudospectral method alone. In addition, pure P-wave equation is completely free from shear-wave, so one can use large discretization step in space to reduce the computational cost further. The stability of the hybrid PS/FD scheme on regular and staggered grids is investigated by von Neumann's method. The analysis shows that the staggered-grid hybrid PS/FD scheme is more stable than regular-grid hybrid PS/FD scheme. The hybrid PS/FD scheme on regular grid is unstable when the anisotropy parameteris ε greater than δ and this instability problem is restricted to wavenumbers near Nyquist. To deal with these high frequency spatial instabilities, a low pass filter is applied to the wavenumber operators in our implementation.We test the hybrid PS/FD scheme with various kinds of VTI model. The numerical simulations provide insight into the kinematic and dynamic accuracy of the pure P-wave VTI equation in comparison with the elastic wave equation. The numerical examples testify that:(1) Pure P-wave equation is completely free from shear-wave compared to elastic wave equation and conventional acoustic VTI equation; (2) The proposed pure P-wave equation can provide good kinematic and dynamic approximations to the elastic wave equation for homogeneous media with weak-to-moderate anisotropy and is more efficient than full elastic wave equation; (3) The hybrid PS/FD scheme on regular grid suffers from some numerical instabilities when the anisotropy parameteris ε greater than δ and the low pass filtering approach do improve the stability; (4) The hybrid PS/FD scheme on staggered grid for pure P-wave equation has better stability than the regular-grid PS/FD scheme. Even in complex heterogeneous media, it still can provide stable and accurate modeling of P-wave propagation. The numerical examples demonstrate the hybrid PS/FD scheme is an effective and efficient method to solve the first-order velocity-stress pure P-wave equation. At the same time, it should be noted that the staggered-grid hybrid PS/FD scheme is more numerically stable than regular-grid hybrid PS/FD scheme, however, its stability condition still is stricter than that of acoustic equation.
Key words: VTI media     Pure P-wave equation     Stability criterions     Finite-difference method     Pseudo-spectral method    
1 引言

随着地震勘探技术的不断发展，地震各向异性被证实普遍存在于地下介质中，特别是在沉积岩地区.通常情况下各向异性是由定向排列的垂直裂缝或周期性薄互层等因素引起的，该类各向异性可以用横向各向同性介质(TI)模型来近似，包括具有垂直对称轴(VTI)和倾斜对称轴(TTI)的横向各向同性介质.近年来，随着计算机技术的快速发展和计算能力的大幅提高，各向异性介质弹性波正演模拟和逆时偏移愈来愈受到人们的重视.但是，由于介质各向异性的复杂性，基于TI介质全波场正演模拟的计算成本仍然很高.为降低计算成本，TI介质正演模拟以及逆时偏移主要基于各向异性声学近似.

一般情况，这些方法可以分为两大类.一类是具有人为横波干扰的耦合波动方程(Alkhalifah，2000；Zhou et al., 2006;Hestholm，2009;Fletcher et al., 2009;Fowler et al., 2010;Duveneck and Bakker， 2011; Zhang et al., 2011)，另一类是完全没有横波干扰的纯P波方程(Klié and Toro，2001; Etgen and Brandsberg-Dahl， 2009; Liu et al., 2009;Crawley et al., 2010;Pestana et al., 2011;Chu et al., 2011，2013; Zhan et al., 2011，2013;黄翼坚等，2011).第一类方法最早由Alkhalifah(2000)提出，由于横波速度对纵波的影响很小，可以人为设定沿对称轴方向上的横波速度为零，由此得到的声学近似方程可以很好地近似描述纵波.声学近似是对各向异性介质波场传播的一种很有效的简化，但是声学近似方程并不能完全去除横波，仍然会残留菱形SV波(Grechka et al., 2004)，并且在介质各向异性参数ε<δ的情况下会引起数值计算的稳定性问题；另外，当TTI介质中对称轴的倾斜角度发生剧烈变化时，声学近似方程同样会出现数值不稳定(李博等，2012).由于以上问题，各向异性纯P波方程逐渐受到众多学者的关注.但是，纯P波方程存在一个分数形式的伪微分算子，直接采用高阶有限差分法求解是十分困难的.伪谱法作为一种地震数值模拟的常用方法(Kosloff and Baysal， 1982;Reshef et al., 1988a，1988b)，可以很好地处理这类具有复杂微分算子的方程.很多学者在伪谱法的基础上提出了一些求解纯P波方程的谱类方法(Etgen and Brandsberg-Dahl， 2009;Crawley et al., 2010;Chu et al., 2013).伪谱法正演模拟过程中需在每个时间迭代过程中进行多次傅里叶变换，计算成本昂贵.为进一步提高效率，Zhan等(2013)联合伪谱法和有限差分法求解二阶纯P波波动方程，并提出了进一步提高效率的计算策略.

为了能够处理变速度、变密度等复杂介质，并结合现有的高阶有限差分方法，本文推导了适用性更广泛的纯P波方程的一阶速度-应力方程形式.借鉴Zhan等(2013)提高效率的计算策略，本文采用伪谱法和高阶有限差分法联合进行一阶速度-应力方程的求解.

对于地震波场数值模拟来说，数值计算的稳定 性是一个重要问题.Virieux(1986)给出了一阶速度-应力方程关于P波和SV波的稳定性条件，Lev and er(1988)、Moczo等(2000)、Carcione等(2002)学者陆续对弹性介质波动方程有限差分解法的稳定性进行了探讨分析.董良国等(2000)对一阶弹性波方程交错网格高阶差分解法的稳定性进行了系统研究；牟永光和裴正林(2005)给出了弹性介质波动方程高阶差分格式的稳定性条件.Furumura等(2002)用并行有限差分法与伪谱法混合方法模拟了我国台湾集集地震的地面强震运动，但是并未对混合法的稳定性问题做详细的讨论.由于复杂微分算子的存在，纯P波波动方程的稳定性与常规波动方程的稳定性并不是完全相同，常规各向异性方程的稳定性条件不再适用.目前关于采用混合法求解纯P波方程的稳定性问题的研究很少，为此本文针对纯P波一阶速度-应力方程的稳定性问题进行了重点讨论，分别分析了中心规则网格和交错网格两种情况下混合法求解此波动方程的稳定性条件.通过理论分析，我们发现在ε>δ的介质中若直接采用中心规则网格混合法求解纯P波方程会引起严重的稳定性问题.针对此问题，我们给出了一种简单有效的解决方法，并对相应的数值正演结果进行了分析.

2 VTI介质纯P波波动方程

首先从精确的VTI介质P-SV相速度出发，二维情况下可以表示为(Tsvankin，1996)

其中，和VS0分别是P波和SV波沿对称轴方向的相速度，θ是相角，ε和δ是各向异性参数(Thomsen，1986)，正号代表P波相速度，负 号代表SV波相速度.根据声学近似假设(Alkhalifah，2000)，令沿对称轴方向上的横波速度VS0=0，那么f=1，因此上式可以简化为

取上式正号项，在弱各向异性条件下对式中的根号项进行一阶泰勒近似，得到纯P波相速度公式：

图 1分别给出了VTI介质弹性P波和纯P波的相速度曲线，以及两者之间的相对误差剖面.其中，图 1a是各向异性参数为ε=0.2，δ=－0.1时的相速度对比图；图 1b是各向异性参数为ε=0.05，δ=0.3时的相速度对比图；图 1c是各向异性参数ε由0变化到0.4，δ由-0.2变化到0.4的弹性P波和纯P波的相速度最大相对误差图.对于纵横波速度比我们选取了常用的1.73，纵波速度取V=3500 m·s-1.图 1a和1b中，黑色点线指示的是弹性P波速度曲线，灰色实线指示的是纯P波速度曲线，从图中可以看到两者吻合得很好，误差非常小.图 1c中的各向异性强度范围基本上处于弱各向异性到中等各向异性之间，从中可以看到在弱各向异性条件下(各向异性参数小于0.2)各向异性纯P波方程的误差基本在1%以内，即使在中等各向异性条件下误差也小于3%.图 2更直观地展示了各向异性纯P方程近似精度与各向异性参数η(Alkhalifah and Tsvankin， 1995)的关系，概括来讲，各向异性纯P波方程在弱各向异性条件下具有非常好的精度，在中等各向异性情况下精度有所下降，但仍具有一定的适用性.

由关系式为角频率，kx和kz为x和z方向上的视波数，可以得到如下的频散公式：

方程(4)即为VTI介质纯P波频散公式，地震波传播过程中没有虚假横波的干扰.附录中，我们给出了声学近似方程和纯声波方程更加直观的联系.显然，上式并不能直接采用有限差分法求解，其中一个方法是通过引入辅助变量，将方程(4)转化为一个耦合的二阶偏微分方程组进行求解，但是在每个时间迭代过程中需要求解一个大型的线性方程组，计算效率比较低(Chu et al., 2011)；相对于求解线性方程组的方法，伪谱法可以很好地处理类似方程(4)的具有复杂微分算子的偏微分方程，并且具有很高的数值精度.

将方程(4)进行傅里叶反变换到时间-空间域，有如下形式：

式中，F[]和F－1[]分别为傅里叶变换和傅里叶逆变换，是虚数单位，本文通过引入一个辅助项ψ，将方程(5)改写为如下的纯P波一阶速度-应力方程:

其中p是应力场，u和w分别为x和z方向质点速度，ψ为辅助项，V(x)是地震波沿z方向传播速度，ρ(x)是介质密度.一阶速度-应力波动方程的主要优点是对于非均匀介质无须对密度和波速进行平滑(Virieux，1984).当各向异性参数为零时，上式退 化为声波方程.与VTI介质弹性波方程相比，纯P波方程形式更简洁，各向异性的物理意义更为直观，无横波干扰.为了提高计算效率，我们联合伪谱法和高阶有限差分法求解(6)式，即傅里叶变换部分采用伪谱法求解，空间导数部分采用高阶有限差分法求解.在数值模拟的过程中由于没有横波，可以取较大的空间步长.

图 1 VTI介质弹性P波和声学近似纯P波的相速度曲线，以及两者之间的相对误差剖面. (a)各向异性参数ε=0.2，δ=－0.1；(b)各向异性参数ε=0.05，δ=0.3；(c)最大相对误差剖面图 1a和1b中黑色点线指示的是弹性P波速度曲线，灰色实线指示的是纯P波速度曲线 Fig. 1 Comparison of phase velocity between elastic wave equation(black dotted line) and pure P-wave equation(gray solid line)：(a)anisotropy parameters ε=0.2，δ=－0.1;(b)anisotropy parameters ε=0.05，δ=0.3; (c)maximum relative error for a wide range of ε and δ
3 稳定性分析

对于地震波场数值模拟来说，数值计算的稳定性是一个重要问题.差分算法离散求解波动方程，其参数的设置对最终的模拟结果有直接的影响，若参数选择不当会导致计算溢出而终止.由于伪微分算子的存在，方程(6)的稳定性与常规波动方程的数值求解的稳定性是不完全相同的，本节主要针对伪谱法和高阶有限差分联合求解方程(6)的稳定性问题进行分析讨论.

图 2 各向异性纯P波方程精度与各向异性参数η的关系图，其中η=(ε－δ)/(1+2δ)Fig. 2 Maximum relative error of pure P-wave equation versus the variation of η，where η=(ε－δ)/(1+2δ)

方程(6)是具有如下形式的一阶双曲型波动方程：

式中，，那么微分算子A为

对(7)式两边进行时间求导，可得

式中

将方程(9)时间导数项采用二阶中心差分代替，并对时间和空间进行傅里叶变换，可得

式中表示单位矩阵，如下式所示:

其中，.

由于方程为双曲型方程组，因此必然存在一个非奇异矩阵S使得成立，其中Λ是由矩阵 的特征值λi构成的正对角矩阵(董良国等，2000)，于是有

令，显然矩阵G和矩阵T相似，那么有det(G)=det(T).要使(11)式有非零解，其系数矩阵行列式必须为零，那么有

由于，所以

通过计算得，那么有.

当且仅当式(13)中的两个不等式同时满足，数值解法才是稳定的.我们首先考察，由于Km的符号不确定，我们取Km的绝对值，那么有

对于中心规则网格有限差分格式，采用泰勒展开式，函数u(x)的一阶导数的2L阶差分近似可以表示为

其中，cl(l=1，2，…，L)为中心规则网格有限差分系数.将差分近似变换到波数域，有ikx=，那么有.分数形式的算子采用伪谱法求解，空间离散化后，取到最大的空间波数——Nyquist波数时，有和，进而可以得到，则

当Δx=Δz=Δh时，有

对于交错网格有限差分法有

其中al(l=1，2，3，…，L)为交错网格有限差分系数.交错网格中参数的定义方式有很多种，图 3为本文交错网格的定义方式，整网格用实线表示，半网格用虚线表示.应力场p，辅助变量ψ，各向异性参数ε和δ，地震波速度V定义在整网格点处，质点震动速度和密度定义在半网格点处.

图 3 交错网格划分示意图Fig. 3 Staggered-grid geometries

由于与中心规则网格情况讨论类似，这里直接给出结果：

接下来，对这一约束条件进行分析.由前面的讨论可知

若Km≥0，则该约束条件无条件成立；若Km＜0，该约束条件未必一定成立.因此下文只针对Km＜0 这种情况进行分析.记

对于中心规则网格有限差分格式，有

对于交错网格有限差分格式，有

直接对(23)式和(24)式进行讨论是很困难的，因此我们首先对两种差分格式的频散关系进行了分析.图 4所示为不同阶数的差分算子在波数域的频散曲线，图 4a为中心规则网格有限差分算子对精确微分算子的逼近，图 4b为交错网格有限差分格式，横坐标是归一化波数，最大值为Nyquist波数.图 5是φ(kx，kz)在差分阶数为2阶，各向异性参数ε=0.3、δ=－0.1的特定介质下的数值模拟图，图 5a是采用中心规则网格有限差分计算出的数值解，图 5b是交错网格有限差分格式.从图 4a可以发现对于中心规则网格有限差分格式，不论取多少阶近似，在靠近Nyquist波数的区域，差分算子响应都接近为0.同时，又由于Lxz＞0，而Km＜0，因此当x和z方向同时取到Nyquist频率时，有φ(kx，kz)＜0， 如图 5a所示，其中四个边角处的值均出现了小于零的情况.在计算过程中这部分高波数区域会引起计算溢出导致波场值呈指数增长，计算终止.因此在ε>δ的情况下不论取多少阶近似，直接联合伪谱法和中心规则网格有限差分法求解方程(6)都无法满足稳定性条件，在正演模拟过程中将出现严重的稳定性问题.

图 4 有限差分算子频散曲线Fig. 4 Dispersion curves for different finite-difference(FD)schemes to approximate first-order derivative

对于交错网格有限差分格式，从图 4b可以看出，当差分阶数提高时差分算子单调递增地逼近微分算子，即对2L阶交错网格有限差分，有

因此我们只需考察L=1的情况即可.当L=1时，a1=1，方程(24)可以写为

其数值解如图 5b所示，其中并没有出现小于零的区域.虽然该数值解是在特定的介质参数下求取的，但是仍能说明联合伪谱法和交错网格有限差分法求解方程(6)有可能是稳定的.当然，这种稳定并不是无条件的，下面我们给出严格的分析.
图 5 φ(kx，kz)的数值模拟(ε=0.3，δ=－0.1，L=1)Fig. 5 Numerical solution of φ(kx，kz) with ε=0.3，δ=－0.1，L=1

空间离散化后，取到最大的空间波数——Nyquist波数时，，方程(26)可以写成

方程(27)给出了一个与时间迭代步长和介质速度无 关，只与介质参数和空间采样间隔有关的稳定性条件：

当Δx=Δz时，上式与空间采样间隔无关，可进一步化简为

不等式(28)实际上限定了采用混合法求解方程(6)所能适用的各向异性强度范围.实际地震勘探中介质一般为弱各向异性，因此上式给出的稳定性条件是可以满足的.另外，提高差分阶数在一定程度上是 可以提高稳定性的，但是效果不是很明显.Alkhalifah(2000)给出的声学近似方程要求ε≥δ，这显然要比不等式(29)要求更严格.在沉积岩层中，各向异性参数ε<δ的情况也是普遍存在的(Thomsen，1986)，此时Alkhalifah(2000)给出的声学近似方程不再适用，而本文的纯P波方程却是稳定的.由此可见，纯P波方程降低了对介质各向异性的要求，扩大了各向异性声学近似的适用范围.

综上讨论，直接联合伪谱法和中心规则网格有限差分法求解一阶各向异性纯P波方程会引起较为严重的稳定性问题，但是联合伪谱法和交错网格有限差分却可以得到相对较好的稳定性，其稳定性条件由不等式(19)和不等式(28)共同给出.事实上，与VTI介质弹性波方程相比较，采用混合法求解纯P波方程的稳定性条件是更为严格的，其中主要的原因在于分数形式的伪微分算子的存在.但是，纯P方程要比VTI介质弹性波方程更简洁，各向异性的物理意义更为直观，无横波干扰.

4 提高稳定性的方法

针对直接联合伪谱法和中心规则网格有限差分法求解一阶各向异性纯P波方程所出现的问题，我们给出了一种解决方案.从图 5a可以看出不等式(21)不满足的区域主要是在高波数区域，那么我们可以通过对高波数区域进行加权衰减来改善稳定性.在实际的计算中，波数的取值范围是从负Nyquist波数到正Nyquist波数，Nyquist波数的大小由空间步长和网格点数通过采样定理来确定.地震信号为带限信号，地震波场中有效信号的绝大部分能量都集中在主频附近的一个窄的频带范围内，这个有效的频带范围仅是正、负Nyquist频率之间很小的一个区间，在这个区间内波场延拓算子的精度决定了波场延拓过程的计算精度.所以对高波数区域可以进行加权处理.值得注意的是，任何带通滤波都会不可避免地产生截断效应，带来一定的噪音污染.这里我们对谱算子Lxz进行窗函数约束.方程(6)可以改写为

其中，D[]表示加权约束.

本文利用有限差分法和Z变换构造波数域窗函数(Yan and Sava， 2009).我们以二阶中心差分格式为例，说明窗函数的构造方法，对于二阶中心差分格式有如下的Z变换形式:

将上式变换到波数域，有

根据上式我们可以定义窗函数

其中，k是归一化波数.其他阶数的窗函数可同样构造，这里列举了4阶、6阶和8阶的窗函数公式：

图 6所示为作用在谱算子Lxz的不同阶数的滤波函数.利用有限差分和Z变换构造窗函数的优点在于其灵活性，即可以根据不同的介质各项异性参数来选择不同阶数的窗函数，当|ε－δ|的值较大时选用较低阶的窗函数，反之选用高阶窗函数，尽量保存高波数区域的信息.

图 6 滤波函数Fig. 6 Filter function
5 数值模拟

为了验证理论分析的正确性以及对比分析各向异性纯P波方程的精度、计算量以及内存需求等问题，我们进行了三组数值模拟测试，前两组数值模拟测例为简单均匀介质，第三组为强变密度、变速度复杂介质.在进行数值模拟的过程中，可以通过下面的处理降低计算成本，方程中的辅助变量ψ可以重新被定义为，相应地，其他定义不变.如此便可将原来的4个傅里叶变换减少为2个，通过简单的推导可以发现改进后的波动方程与未改进的方程具有相同的稳定性条件，附录B给出了说明.

首先，在各向异性参数为ε=0.2，δ=-0.05和ε=0.05，δ=0.2的均匀介质中进行数值模拟，同时与各向异性介质弹性波方程和Duveneck等(2008)提出的各向异性拟声波方程进行对比.在数值模拟时，时间采用2阶有限差分格式，空间采用10阶有限差分格式，边界条件采用完全匹配层(PML)边界条件(Berenger，1994).模型计算区域大小为350×400的四边形区域，网格间距为Δx=Δz=10 m，沿垂直对称轴速度V=3000 m·s-1，时间步长dt=0.001 s.震源子波采用雷克子波，子波中心频率fm=25 Hz，震源置于模型中间.图 7是第一组VTI介质(ε=0.2，δ=-0.05)中弹性波场垂直分量波场快照(图 7a)、各向异性拟声波垂直分量波场快照(图 7b)和各向异性纯P波(图 7c、d和e)垂直分量波场快照，其中，图 7c采用了规则网格有限差分混合法，图 7d对规则网格混合法中的伪微分算子进行了8阶窗函数约束，图 7e采用了交错网格有限差分混合法.从图 7可以看出：(1)相对于弹性波方程和各向异性声学近似方程，各向异性纯P波方程完全无SV波干扰；(2)在ε>δ的VTI介质中，直接联合伪谱法和规则网格有限差分法求解一阶各向异性纯P波方程会引起较为严重的稳定性问题，波场呈指数增长，计算溢出，对伪微分算子进行加窗约束后稳定性得到明显改善；(3)交错网格有限差分混合法对一阶各向异性纯P波方程具有较好的稳定性.图 8是第二组VTI介质(ε=0.05，δ=0.2)中弹性波场(图 8a)和各向异性纯P波(图 8b和c)的垂直分量波场快照，其中，图 8b采用了规则网格有限差分混合法，图 8c采用了交错网格有限差分混合法.在ε<δ的VTI介质中，由于不满足VTI 介质拟声波方程的稳定性条件，是无法进行数值模拟的(Alkhalifah，2000)，而纯P波波动方程可以得 到稳定的数值模拟结果，如图 8b和图 8c所示，这与前面的理论分析是一致的，即在ε<δ的介质中规则网格有限差分混合法和交错网格有限差分混合法对一阶各向异性纯P波方程都具有较好的稳定性.概括来说，交错网格有限差分混合法对一阶各向异性纯P波方程具有较好的稳定性，而规则网格有限差分混合法会出现数值稳定性问题，需要对伪微分算子进行窗函数约束.

图 7 各向异性VTI介质(ε=0.2，δ=-0.05)波场快照：(a)弹性波场垂直分量；(b)各向异性声学近似方程模拟波场；(c)规则网格混合法求解纯P波方程波场；(e)伪谱算子滤波后规则网格混合法求解纯P波方程波场；(d)交错网格混合法求解纯P波方程波场 Fig. 7 Wavefield snapshots in a VTI medium with ε=0.2，δ=-0.05：(a)Vertical component synthesized by elastic wave equation;(b)Wavefield synthesized by anisotropic acoustic wave equation;(c)Wavefield synthesized by pure P-wave equation using hybrid PS/CFD scheme;(d)As for(c)，except filtering the pseudo-differential operator;(d)As for(c), except using hybrid PS/SGFD scheme

图 8 各向异性VTI介质(ε=0.05，δ=0.2)波场快照：(a)弹性波场垂直分量；(b)规则网格混合法求解纯P波方程波场；(c)交错网格混合法求解纯P波方程波场 Fig. 8 Wavefield modeling in a VTI medium with ε=0.05，δ=0.2：(a)Vertical component synthesized by elastic wave equation;(b)Wavefield synthesized by pure P-wave equation using hybrid PS/CFD scheme;(c)As for(b), except using hybrid PS/SGFD scheme

图 9是前两组VTI介质中弹性波和纯P波的波形对比图，其中黑色实线是弹性波，灰色虚线是纯P波.从图中可以看出，除了没有SV波之外，VTI介质纯P波方程与弹性波方程模拟波场的波形振幅与传播规律十分吻合，在小偏移距内两者几乎完全相同，这说明纯P波方程对均匀介质中P波的运动学和动力学特征的描述是比较准确的.然而，值得注意的是，纯P波方程是不存在能量转换的，因而纯P波方程是无法处理非均匀介质分界面处转换波的问题，如何保证其动力性特征还需要进一步地深入探索研究.

图 9 弹性波(黑色实线)和纯P波(灰色虚线)波形对比图Fig. 9 Comparison between elastic(solid black line) and pure P wave(dashed gray line)seismograms

为了验证本文纯P波一阶应力-速度方程对复杂介质的适用性及长时间数值正演中的稳定性，我们采用变密度Hess VTI模型进行正演模拟测试.图 10所示为Hess VTI模型(部分).我们对该模型进行了重采样，模型计算区域大小为750×1200，网格间距为Δx=Δz=12.5 m，时间步长dt=0.001 s.震源子波中心频率fm=25 Hz，震源位于(50，650)处.图 11a所示为各向异性弹性波方程单炮垂直分量波场快照；图 11b是交错网格混合法模拟的各向异性纯P波垂直分量波场快照，图 12是与之对应的总时间为6s的共炮点道集.从波场快照和地震记录可以看出，对于复杂介质各向异性纯P波方程仍旧可以得到比较稳定的模拟结果，强速度分界面处的传播也十分稳定，并且相对于弹性波，纯 P波波场更加清晰，频散误差更小.这对于以垂直分 量为主的常规单分量地震资料偏移成像是非常有利的.

图 10 VTI Hess模型Fig. 10 VTI Hess model

图 11 Hess模型不同方程正演模拟波场快照： (a)VTI介质弹性波方程正演结果；(b)VTI介质纯P波方程正演结果Fig. 11 Wavefield snapshots of vertical component for Hess Model synthesized by(a)elastic wave equation and (b)pure P-wave equation using hybrid PS/SGFD scheme

图 12 Hess模型不同方程正演模拟地震记录：(a)VTI介质弹性波方程正演结果；(b)VTI介质纯P波方程正演结果；(c)和(d)分别是(a)和(b)的矩形区域标出的部分Fig. 12 Shot gathers generated form Hess model：(a)Seismic profile with elastic wave equation；(b)Seismic profile with pure P-wave equation by using hybrid PS/SGFD scheme；(c)Zoom of rectangle area in(a)；(d)Zoom of rectangle area in(b)

为了对比各向异性纯P波的计算量和内存需求，针对Hess模型(部分)我们还进行了各向异性拟声波方程的数值测试.由于模型大小、边界条件的选取、算法的优化、有限差分的阶数以及计算机硬件设备均会影响到计算效率，因此我们除了采用的方程不同之外，其他参数设置完全相同，计算设备为4核Intel·CoreTM i5处理器，观测时间为6000个步长.表 1给出了三种不同波动方程单炮平均耗时和内存消耗.纯P方程同各向异性拟声波方程在计算区域需要4个变量，弹性波方程需要5个变量，边界区域由于采用分裂PML边界，弹性波方程比纯P波方程和拟声波方程会消耗更多的内存.从测试结果可以看出，各向异性纯P波方程的计算效率比各 向异性声学近似方程要低一些，但相比全弹性波方程，计算效率明显得到提高.实际上，相比弹性波方程和各向异性声学近似方程，各向异性纯P波方程的计算优势在于没有横波频散限制，可以取到更大的空间步长，以此减少计算量，提高计算效率，降低内存需求.

表 1 三种不同波动方程的单炮平均耗时和内存消耗Table 1 The average time costs and memory requirements of three different equations
6 结论

基于VTI介质，推导出了纯P波一阶速度-应力方程，分析了VTI介质纯P波方程在运动学和动力学上的近似精度.在数值模拟的过程中采用伪谱法和高阶有限差分法联合求解该波动方程，重点讨论了中心规则网格和交错网格两种差分格式下混合法的稳定性条件，给出了其适用的各向异性强度范围.纯P波方程降低了对介质各向异性的要求，扩大了各向异性声学近似的适用范围，在介质各向异性参数ε<δ的情况下仍然具有很好的稳定性.

针对中心规则网格稳定性问题，本文给出了一种解决方案，通过对谱算子进行窗函数约束来提高数值稳定性，数值模拟结果表明了该方法的正确性.与规则网格有限差分混合法相比较，交错网格有限差分混合法具有更好的稳定性和更高的精度.复杂介质数值模拟结果表明了纯P波一阶速度-应力方程可以有效处理强变速、变密度介质.

致谢 感谢Hess公司与SEG提供二维VTI模型，感谢评审专家反馈的意见和建议，本文曾得到中国石油大学(华东)王书亭教授以及国土资源部中国地质调查局青岛海洋地质研究所方刚博士的帮助，谨此致谢.

附录A 两类声学近似的关系

这个附录简要说明了两类声学近似的关系.Alkhalifah(2000)给出的频散关系如下：

式中，VP0是P波沿对称轴速度，ω是角频率，VNMO=√1+2δ VP0 是动校正速度，η=(ε－δ)/(1+2δ)是非椭圆率，当η=0时P波波前面为椭圆形状，kx，ky和kz是沿x，y和z轴的视波数.方程(A1)是关于ω2的二次方程，其解的形式与表达式(2)类似，分别对应着一个P波和菱形SV波解.经过简单代数运算，可以把方程(A1)写成下式：

纯P波方程的3D形式如下所示：

对比上面两个式子可以发现，对于分数形式算子的分母项存在下面的近似：

这说明纯P波方程是声学假设的一种椭圆再近似.实际上，纯P波方程中可以存在不为零的横波速度，只是横波速度对纵波的影响非常小，故本文沿用了Alkhalifah(2000)的做法将其舍去.通过对公式(1)做泰勒近似同样可以得到纯SV波的波动方程，只是一阶近似的精度较低，可以采用较为稳定的高阶Padé连分式逼近来提高近似精度，通常情况下二阶近似就可以得到非常好的效果.

附录B 改进的纯P波一阶速度-应力方程

在数值计算的过程中，我们对方程做了进一步改进，将辅助变量ψ重新定义为，其他定义不变，相应有.虽然改进后的方程并不是严格意义上的一阶方程了，但是如此处理可以将原来的4个傅里叶变换减少为2个，降低了计算成本.下面对改进的P波一阶应力速度方程的稳定性做简要说明.改进后的P波一阶应力速度方程仍然可以采用方程(7)的写法，只是其中的微分算子A稍有不同，记为

那么

通过计算可知(A′)2和A 2具有相同的特征值，即，因而改进后的方程和方程(6)具有相同的稳定性条件.

参考文献
[1] Alkhalifah T, Tsvankin I. 1995. Velocity analysis for transversely isotropic media. Geophysics, 60(5):1550-1566.
[2] Alkhalifah T. 2000. An acoustic wave equation for anisotropic media. Geophysics, 65(4):1239-1250.
[3] Berenger J P. 1994. A perfectly matched layer for absorption of electromagnetic waves. Journal of Computational Physics, 114(2):185-200.
[4] Carcione J M, Herman G C, Ten Kroode A P E. 2002. Seismic modeling. Geophysics, 67(4):1304-1325.
[5] Chu C L, Macy B K, Anno P D. 2011. An accurate and stable wave equation for pure acoustic TTI modeling. 81th Annual International Meeting. SEG, Expanded Abstracts, 179-184.
[6] Chu C L, Macy B K, Anno P D. 2013. Pure acoustic wave propagation in transversely isotropic media by the pseudospectral method. Geophysical Prospecting, 61(3):556-567.
[7] Crawley S, Brandsberg-Dahl S, McClean J, et al. 2010. TTI reverse time migration using the pseudo-analytic method. The Leading Edge, 29(11):1378-1384.
[8] Dong L G, Ma Z T, Cao J Z. 2000. A study on stability of the staggered-grid high-order difference method of first-order elastic wave equation. Chinese J. Geophys. (in Chinese), 43(6):856-864.
[9] Duveneck E, Milcik P, Bakker P M, et al. 2008. Acoustic VTI wave equations and their application for anisotropic reverse-time migration. 78th Annual International Meeting. SEG, Expanded Abstracts, 2186-2190.
[10] Duveneck E, Bakker P M. 2011. Stable P-wave modeling for reverse-time migration in tilted TI media. Geophysics, 76(2):S65-S75.
[11] Etgen J T, Brandsberg-Dahl S. 2009. The pseudo-analytical method:Application of pseudo-Laplacians to acoustic and acoustic anisotropic wave propagation. 79th Annual International Meeting. SEG, Expanded Abstracts, 2552-2556.
[12] Fletcher R P, Du X, Fowler P J. 2009. Reverse time migration in tilted transversely isotropic (TTI) media. Geophysics, 74(6):WCA179-WCA187.
[13] Fowler P J, Du X, Fletcher R P. 2010. Coupled equations for reverse time migration in transversely isotropic media. Geophysics, 75(1):S11-S22.
[14] Furumura T, Koketsu K, Wen K L. 2002. Parallel PSM/FDM hybrid simulation of ground motions from the 1999 Chi-Chi, Taiwan, earthquake. Pure Appl Geophys, 159(9):2133-2146.
[15] Grechka V, Zhang L B, Rector J W Ⅲ. 2004. Shear waves in acoustic anisotropic media. Geophysics, 69(2):576-582.
[16] Hestholm S. 2009. Acoustic VTI modeling using high-order finite differences. Geophysics, 74(5):T67-T73.
[17] Huang Y J, Zhu G M, Liu C Y. 2011. An approximate acoustic wave equation for VTI media. Chinese J. Geophys. (in Chinese), 54(8):2117-2123, doi:10.3969/j.issn.0001-5733.2011.08.019.
[18] Klié H, Toro W. 2001. A new acoustic wave equation for modeling in anisotropic media. 71st Annual International Meeting. SEG, Expanded Abstracts, 1171-1174.
[19] Kosloff D D, Baysal E. 1982. Forward modeling by a Fourier method. Geophysics, 47(10):1402-1412.
[20] Levander A R. 1988. Fourth-order finite-difference P-SV seismograms. Geophysics, 53(11):1425-1436.
[21] Li B, Li M, Liu H W, et al. 2012. Stability of reverse time migration in TTI media. Chinese J. Geophys. (in Chinese), 55(4):1366-1375, doi:10.6038/j.issn.0001-5733.2012.04.032.
[22] Liu F Q, Morton S A, Jiang S S, et al. 2009. Decoupled wave equations for P and SV waves in an acoustic VTI media. 79th Annual International Meeting. SEG, Expanded Abstracts, 2844-2848.
[23] Moczo P, Kristek J, Bystricky E. 2000. Stability and grid dispersion of the P-SV 4th-order staggered-grid finite-difference schemes. Studia Geophysica et Geodaetica, 44(3):381-402.
[24] Mu Y G, Pei Z L. 2005. Seismic Numerical Modeling for 3-D Complex Media (in Chinese). Beijing:Petroleum Industry Press, 33-34.
[25] Pestana R C, Ursin B, Stoffa P L. 2011. Separate P-and SV-wave equations for VTI media. 12th International Congress of the Brazilian Geophysical Society, 1227-1232.
[26] Reshef M, Kosloff D, Edwards M, et al. 1988a. Three-dimensional acoustic modeling by the Fourier method. Geophysics, 53(9):1175-1183.
[27] Reshef M, Kosloff D, Edwards M, et al. 1988b. Three-dimensional elastic modeling by the Fourier method. Geophysics, 53(9):1184-1193.
[28] Thomsen L. 1986. Weak elastic anisotropy. Geophysics, 51(10):1954-1966.
[29] Tsvankin I. 1996. P-wave signatures and notation for transversely isotropic media:An overview. Geophysics, 61(2):467-483.
[30] Virieux J. 1984. SH-wave propagation in heterogeneous media:velocity-stress finite-difference method. Geophysics, 49(11):1933-1942.
[31] Virieux J. 1986. P-SV wave propagation in heterogeneous media:Velocity-stress finite-difference method. Geophysics, 51(4):889-901.
[32] Yan J, Sava P. 2009. Elastic wave-mode separation for VTI media. Geophysics, 74(5):WB19-WB32.
[33] Zhan G, Pestana R C, Stoffa P L. 2011. An acoustic wave equation for pure P wave in 2D TTI media. 12th International Congress of the Brazilian Geophysical Society, 993-997.
[34] Zhan G, Pestana R C, Stoffa P L. 2013. An efficient hybrid pseudo-spectral/finite-difference scheme for solving the TTI pure P-wave equation. J. Geophys. Eng., 10(2):025004.
[35] Zhang Y, Zhang H Z, Zhang G Q. 2011. A stable TTI reverse time migration and its implementation. Geophysics, 76(3):WA3-WA11.
[36] Zhou H, Zhang G, Bloor R. 2006. An anisotropic acoustic wave equation for VTI media. 68th EAGE Conference and Exhibition incorporating SPE EUROPEC 2006.
[37] 董良国, 马在田, 曹景忠. 2000. 一阶弹性波方程交错网格高阶差分解法稳定性研究. 地球物理学报, 43(6): 856-864.
[38] 黄翼坚, 朱光明, 刘池洋. 2011. 一个近似的VTI介质声波方程. 地球物理学报, 54(8): 2117-2123, doi: 10.3969/j.issn.0001-5733.2011.08.019.
[39] 李博, 李敏, 刘红伟等. 2012. TTI 介质有限差分逆时偏移的稳定性探讨. 地球物理学报, 55(4): 1366-1375, doi: 10.6038/j.issn.0001-5733.2012.04.032.
[40] 牟永光, 裴正林. 2005. 三维复杂介质地震数值模拟. 北京: 石油工业出版社, 33-34.