随着地震勘探的目标日趋复杂，与对勘探精度要求的越高，各向异性介质已成为当今地球物理领域的研究热点.而正演模拟是研究各向异性介质中地震波场传播规律的重要工具，其结果也是检验地震属性研究、反演、地震资料处理和解释等结论的主要依据.

Alkhalifah提出了一种VTI介质下的近似频散关系式，即将沿对称轴方向的横波速度设置成零，忽略了横波的物理现象，却非常精确地描述了 qP 波在 VTI 介质中传播的波动特征(Alkhalifah，1998).相对于求解相对于求解VTI 介质全弹性波方程，利用声波近似方程的更为实用高效的方法，也更接近实际生产需求，引起了学者们的广泛关注(Klie and Toro， 2001；Du et al．，2008；Duveneck et al．，2008；Fowle et al．，2010；张岩和吴国忱，2013a).

然而该方法的主要问题是会产生假象波场，这种假象由人为设定横波速度为零时产生的伪SV波引起，因为简单地令横波速度等于零不会使横波相速度在声学近似VTI介质中的每一处都等于零(Grechka et al．，2004)，由此还产生了数值计算的稳定性问题，即需满足Thomsen参数满足ε>δ时，才能使此部分伪SV波分量稳定传播.有学者提出了多种解决方法：如当震源位置放置在各向异性介质之上的各向同性介质中，假象便会消失(吴国忱和梁锴，2005)；或者令ε=δ，简化为椭圆VTI介质，也可以起到压制SV波假象的作用，这类方法显然限制了声波近似方程的适用范围；另一种思路是采用算子滤波的方法对波场进行滤波压制伪SV波分量(Zhang，2009；张岩和吴国忱，2013b)，但上述这些方法仍不能彻底解决稳定性问题.

因此近年也有学者提出了将方程解耦提出纯P波方程的方法.最初是对VTI介质准P波方程的频散关系式直接进行近似展开分解后得到了纵横波解耦的近似方程(Liu et al., 2009)，但其中的横波分量仍需满足ε>δ才可适定.最近有学者则在VTI介质精确频散关系式基础上(Tsvankin，1996)，重新导出了解耦合的VTI介质纯P波与纯SV波近似方程(Pestana et al．，2011)，其推导过程中无需设置横波速度为零，能够较精确地描述准P波的运动学特征，同时解决伪横波干扰问题，更具有一般性.

利用有限差分方法，很多学者已经实现了计算高效的各向异性VTI介质方程的正演模拟算法(何兵寿和张会星，2006；Zhou et al．，2006；吴国忱和梁锴，2007；Hestholm，2007).但对于VTI解耦合声波方程中的分母偏导项及参数变量使得有限差分方法难以直接应用.伪谱法的在空间频散较低，很适合计算方程中空间混合偏导，实现过程直观简便，可直接利用频率域下的方程进行递推，但伪谱法在时间递推上为二阶差分精度不高，使得仍有频散和不稳定现象产生，且计算效率不高.伪解析法(pseudo-analytical method)波场延拓理论在伪谱法的基础上进行了改进(Etgen and >Br and sberg-Dahl， 2009).它利用伪谱法本身在波数域下运算的特点，通过时间差分的波动方程在波数域下推导出一种伪拉普拉斯算子，使得时间差分所导致的误差能够得到补偿，提高了伪谱法的时间递推精度.

本文首先研究推导了VTI介质纵横波解耦方程，并在前人研究的基础上提出采用伪解析法建立VTI介质解耦合声波方程的正演数值模拟方法，解决了数值计算的效率问题.该方法相当于改进的伪谱法，通过在频率域求解解决空间微分，同时利用伪拉普拉斯算子提高对时间导数项的逼近精度，递推产生的数值频散更小，稳定性条件要求低，因此可以采用更大的延拓步长，较传统伪谱法有明显的计算效率和精度优势.脉冲响应试验和Hess模型的正演模拟结果证实了方法的正确性和有效性.

VTI介质精确频散关系式为(Tsvankin，1996)

P波：

SV波：

P波：

SV波：

根据泰勒展开式的收敛条件，VTI介质下的纯P波、纯SV波方程，当P波和SV波频散关系满足以下条件时，对精确公式有比较好的近似.

图 1

Fig. 1

图 1 正演模拟波场快照对比(ε=0.25，δ=0.1) (a)准P波方程伪谱法正演波场快照；(b)准P波方程伪解析法正演波场快照；(c)纯P波方程伪谱法正演波场快照； (d)纯P波方程伪解析法正演波场快照. Fig. 1 Impulse responses modeling wavefield snapshots comparison(ε=0.25，δ=0.1) (a)Quasi P-wave equation wavefield snapshot with pseudo-spectral method；(b)Pure P-wave equation wavefield snapshot with pseudo-spectral method；(c)Quasi P-wave equation wavefield snapshot with pseudo-analytical method；(d)Pure P-wave equation wavefield snapshot with pseudo-analytical method.

图 1

Fig. 2

图 2 正演模拟波场快照对比(ε=0.1，δ=0.25) (a)qP-qS波方程伪谱法正演波场快照;(b)qP-qS波方程伪解析法正演波场快照; (c)纯P波方程伪谱法正演波场快照;(d)纯P波方程伪解析法正演波场快照. Fig. 2 Impulse responses modeling wavefield snapshots comparison(ε=0.1，δ=0.25) (a)Quasi P-wave equation wavefield snapshot with pseudo-spectral method;(b)Pure P-wave equation wavefield snapshot with pseudo-spectral method;(c)Quasi P-wave equation wavefield snapshot with pseudo-analytical method; (d)Pure P-wave equation wavefield snapshot with pseudo-analytical method.

伪解析法波场延拓理论主要的思路是，将二阶时间差分的波动方递推式变换到波数域，推导出伪拉普拉斯算子，即通过修改拉波动方程中的普拉斯算子项 Δ ²，补偿时间二阶差分的所带来的误差.

首先，各向同性声波方程可用如下时间二阶精度差分递推式求解

对于VTI介质纯P波方程而言，只需要(13)式中的ω替换成(6)式推导的频散关系式即可，但考虑速度及各向异性参数的变化，这里通过分别针对k²_z，k²_r，k²_rk²_z/(k²_z+Fk²_r)三个部分，求得不同的F(k)，再经过傅氏反变换求出三部分波场值加以组合得到VTI介质纯P波方程的延拓方法.

相对于传统的伪谱法法波场延拓方法，伪解析法主要有以下几个特点：

(1)通过求取包含了时间递推误差伪拉普拉斯算子，替代原来伪谱法中的拉普拉斯算子，从而抵消了时间差分递推误差，在速度均匀的情况下，可以达到近乎解析解的精度，由此称其为伪解析法.

(2)由于伪解析法是在速度均匀假设下导出的，为了适应各向异性变速介质，可以通过先求取参考相速度下的F(k)值(如最大、最小速度)，再利用线性插值方法获得各个位置不同相速度下的波场值.由于研究表明F(k)值与速度变化的关系不大(Chu 2010)，该方法能够保证一定的精度.

为了验证伪解析法VTI介质解耦合声波方程数值模拟的有效性，首先进行了点源响应试验，参数设置为：v_po=3000 m/s；各向异性参数ε=0.25，δ=0.1；空间采样间隔15 m，时间采样间隔为2 ms，波场快照时间0.4 s.

这里分别进行了VTI介质准P波方程伪谱法与伪解析法正演模拟，采用的是Du等提出的耦合形式方程(Du et al．，2008)：

从图 1a和图 1c中，可以看出VTI介质准P波方程与纯P波方程波场快照中的波形基本一致，但准P波方程波场快照图 1a、图 1b的结果中间存在菱形的伪SV波假象，而图 1c、图 1d所示的纯P波方程结果则彻底地消除了这些假象.另一方面，分别对比图 1a、c与图 1b、d，可以看出伪解析法数值模拟结果明显较伪谱法数值频散更小，伪解析法解决了伪谱法由于时间递推误差导致的频散问题.最后，上述四种波场快照结果中的P波波形和振幅都具有很好的一致性.

其次，为了对比VTI介质纯P波方程与准P波方程的稳定性，采用之前相同的参数和方式进行正演模拟，仅将其中的各向异性参数改为ε=0.1，δ=0.25.此时，由于各向异性参数ε<δ，未满足VTI介质准P波近似方程的稳定性条件，伪SV波频率呈指数增长，导致伪谱法与伪解析法数值模拟结果不稳定，已经无法得到其波场快照.为了验证正确性，这里采用了不设置SV速度为零的VTI介质qP-qS波方程(17)(Fowler et al．，2010)，与VTI介质纯P波方程的正演结果进行对比.正演设置垂向横波速度v_so=1000 m/s.而若v_so=0则(17)式方程与(16)式是完全等价的.正演同样获得了四种不同的波场快照，结果如图 2.

从波场快照结果中可以看到，VTI介质纯P波方程图 2c、图 2d与qP-qS波方程图 2a、图 2b模拟的qP波波形是完全一致的.一方面，VTI介质纯P波消除了伪SV波干扰，其近似方程仍能够保证稳定地传播.同时无论是伪谱法还是伪解析法，纯P波方程的模拟结果波场频散都比qP-qS波方程结果更小.

图 3

Fig. 3

图 3 Hess模型点源正演模拟波场快照对比 (a)伪谱法波场快照;(b)伪解析法波场快照 Fig. 3 Impulse responses wave field snapshots comparison in Hess model (a)Pseudo-spectral wavefield snapshot; (b)Pseudo-analytical wavefield snapshot.

图 4

Fig. 4

图 4 Hess模型不同方程正演记录 (a)VTI介质准P波耦合方程正演结果；(b)VTI介质纯P波方程正演结果. Fig. 4 Hess model modeling profiles with different equations (a)VTI medium quasi P-wave coupling equations modeling result;(b)VTI medium pure P-wave equation modeling result.

综合前述两组点源响应试验对比，可以看出本文推导的基于伪解析法的纯P波近似方程正演模拟方法能够较准确地描述VTI介质中qP波的运动学及动力学特征，且完全消除了伪SV波假象，保证了波场在不同各向异性介质参数情况下的稳定传播，数值频散也更小.

为了进一步检验在非均匀复杂介质情况下的伪解析法VTI介质解耦合声波方程数值模拟效果.首先这里采用了SEG提供的Hess VTI各向异性模型进行正演模拟试验，模型横向距离22050 m，纵向深度9144 m.首先进行了点震源模拟，震源位置位于模型内部(11000 m，3000 m)处，纵横向网格为25 m，时间采样间隔为2 ms.分别为伪谱法和伪解析法在相同参数下的Hess模型正演模拟波场快照如图 3.

通过插值组合的伪解析法很好得适应非均匀下的VTI介质数值模拟.对比可知，伪谱法出现了较为明显的频散(图 3a)，而伪解析法则仍能够保证较小的频散传播波场(图 3b).因此，在相同计算精度下，伪解析法可以用更大的时间步长，从而降低正演模拟过程的计算量.

其次则是VTI介质准P方程与VTI介质纯P波方程在Hess模型中的正演模拟结果对比.正演模拟采用的参数：时间采样间隔1 ms，采样时间长度为6 s，纵横向网格采样间隔为12.129 m，震源位于(10911 m，0 m)，采用30 Hz雷克子波，采用炮点右边500道接收，道间距12.129 m.吸收边界采用了Cerjan等提出的指数衰减吸收边界(Cerjan et al., 1985)，同时在边界吸收层中置入衰减速度场，经试验可获得较好的效果，正演结果如图 4.

其中，图 4a为VTI介质准P波方程正演模拟结果，图 4b则为纯P波方程正演模拟结果，正演方法均为伪解析法.可以看出，解耦合的VTI介质纯P波方程与传统准P波方程模拟波场的波形振幅与传播规律十分吻合.同时对比可以发现，相同的正演模拟参数情况下，VTI介质纯P波方程模拟结果的频散误差也较准P波方程更小.原因在于传统VTI介质准P波方程中速度值相对较低的伪横波分量导致了更大的频散误差，而VTI介质纯P波方程已经不包含伪横波分量，数值模拟过程更稳定，频散误差也更小.

推导了VTI介质解耦合声波近似方程，通过利用伪解析法解决数值求解问题.VTI介质解耦声波方程克服了传统VTI介质准P波方程的缺点，在保证准确描述qP波场特征的情况下，消除了伪SV波干扰，使得波场传播更稳定，数值频散也更小.而伪解析法则在保证了一定的精度情况下提高了计算效率.

对比伪谱法和伪解析法两种算法的数值模拟效果，可以看出伪解析法数值模拟结果的频散更小，故在相同计算精度下，可以用更大的时间步长，以减少计算量，该波场数值模拟的延拓方法也可以应用到逆时偏移中，或扩展到TTI介质等情况中去.

目前伪解析法同伪谱法一样，实际计算效率主要取决于二、三维傅里叶变换的计算效率，随着计算机技术的发展，利用现有的CPU/GPU并行计算(刘红伟等，2010)、或并行矢量化傅里叶变换优化程序，使得计算效率得以提高.