地震各向异性在地球介质中是普遍存在的(Tsvankin，2001; Tsvankin et al．，2010)，如墨西哥湾、西非洲海岸和中国南海都有大量的各向异性地层存在，而横向各向同性(TI)介质是最常见的一种各向异性介质(Thomsen，1986).

Thomsen(1986)提出了著名的弱各向异性理论，Tsvankin(1996)基于Thomsen(1986)的弱各向异性理论推导得到P-SV波VTI介质相速度频散方程，这是后来很多学者研究TI介质正演模拟和偏移成像的出发点.

TI介质正演数值模拟是TI介质逆时偏移(张岩和吴国忱，2013)和TI介质全波形反演(Warner et al．，2013；Alkhalifah，2014)的基础.TI介质正演数值模拟一般在时间域进行，也可以在频率空间域进行，频率空间域正演的优势在于可以多炮同时模拟且不存在累积误差，此时频率空间域波动方程的求解问题变成了一个大型稀疏线性方程组的求解问题.近年来若干学者对TI介质频率空间域正演进行了详细的研究(Grini et al．，2007；Operto et al．，2007，2009；吴国忱等，2007；李桂花等，2011).吴国忱和梁锴(2005，2007)通过引入25点加权优化差分算子，实现了VTI介质的频率空间域正演模拟.吴国忱等(2007)将这种方法推广到TTI介质，推导出二维TTI介质频率空间域弹性波动方程，通过引入25点加权优化差分算子，采用高斯牛顿方法(Lines and Treitel，1984；Min et al．，2000)求取优化系数，从而在频率空间域模拟了弹性波在TTI介质中的传播过程.杜向东等(2009)将吴国忱和梁锴(2005，2007)的方法推广到三维，通过引入125点优化差分算子，实现了三维VTI介质的频率空间域正演模拟.李桂花等(2011)通过引入25点优化差分算子，实现了黏弹性VTI介质的频率空间域正演模拟.Grini等(2007)和Operto等(2007，2009)研究了黏声TTI介质一阶方程形式的混合网格频率空间域正演方法.Operto等(2014)将VTI介质四阶方程分解为一个椭圆型二阶方程和非椭圆校正项的形式，实现了三维VTI介质的频率空间域正演模拟.

以相速度误差1%为标准，VTI介质qP波波动方程常规9点频率空间域差分格式需要每个波长12个网格点才能正确模拟，计算精度差.Chen(2012)提出了一种平均导数优化方法的新思路，能有效提高计算精度，而且具有很好的灵活性，能适用于横纵向网格不等间距的情形.张衡等(2014)基于Chen(2012)利用平均导数优化方法构建差分格式的新思路发展了一种新的基于平均导数优化方法的声波方程频率空间域四阶25点格式，能将每个波长所需要的网格点数降低到2.78个，相比四阶常规9点方法明显提高了计算精度，也克服了Shin和Sohn(1998)提出的四阶25点混合格式只能适用于横纵向网格等间距情形的局限性.本文将平均导数优化方法进一步推广到VTI介质qP波波动方程，研究了基于平均导数优化方法的VTI介质频率空间域二阶9点算法.

本文首先推导了二维VTI介质频率空间域qP波波动方程，然后提出了基于平均导数优化方法(ADM)的二维VTI介质qP波波动方程频率空间域二阶9点格式.随后通过最小二乘优化方法求取优化系数并做了频散分析，表明VTI ADM 9点算法能够将每个波长所需的网格点数降低到3.57个，相比VTI常规9点格式每个波长需要12个网格点数明显提高了计算精度.最后进行数值模拟，测试表明VTI ADM 9点算法能实现高精度模拟.

利用VTI介质弹性张量矩阵，结合弹性动力学的本构方程(广义虎克定律)、运动微分方程(牛顿第二定律的微分形式)和几何方程(描述位移与应变关系)可得到VTI介质中的弹性波动方程.将平面波方程U=Ae^{i(k_xx+k_yy+k_zz－ωt)}代入VTI介质中的弹性波动方程并忽略体力项，推导得到Christoffel方程(吴国忱，2006；程玖兵等，2013).

由Christoffel方程出发，并利用TI介质中的Thomsen各向异性参数ε和Δ推导得到P-SV波VTI介质相速度频散方程(Tsvankin，1996；Alkhalifah，1998，2000)：

其中：.

式(1)中V_p(θ)表示相速度，θ表示相速度传播角度，f中V_Sz表示横波速度，V_Pz表示纵波速度，ε和Δ表示Thomsen各向异性参数(Thomsen，1986).

采用Alkhalifah(1998)提出的TI声波近似思想(即假设横波速度V_Sz为零)基于式(1)即可推导得到VTI介质频率空间域qP波波动方程.

令横波速度V_Sz=0，则f=1，代入式(1)得到：

进一步对式(2)进行代数变换并利用k_x²+k_z²=变换到频率波数域，即可得到二维VTI介质频率波数域qP波方程：

当Thomsen各向异性参数ε=δ时式(4)退化为二维椭圆各向异性介质方程：

因此二维椭圆各向异性介质方程和二维各向同性介质方程可看作二维VTI介质频率空间域qP波波动方程的特例.

式(4)表示的二维VTI介质频率空间域qP波波动方程可表示成VTI常规9点格式的形式(图 1a)，其差分格式为：

图 1

Fig. 1

图 1 (a) VTI常规9点差分格式示意图；(b) VTI ADM 9点差分格式示意图Fig. 1 (a) VTI conventional 9-point scheme；(b) VTI ADM 9-point scheme

但是VTI常规9点格式正演精度低，数值频散严重(图 2)，需要发展新的优化算法.本文基于Chen(2012)和张衡等(2014)采用平均导数优化方法构建声波方程频率空间域差分格式的思路，提出一种新的基于平均导数优化方法的二维VTI介质qP波波动方程频率空间域二阶9点有限差分格式(图 1b).

对式(4)二维VTI介质频率空间域qP波波动方程中的二维规则空间导数项和加速度项分别引入加权优化系数，二维规则空间导数项的二阶差分近似的波场值表示为正交方向上3个网格点波场值的加权平均形式，而加速度项ω⁴P则表示为ADM 9点差分格式中所有9点的加权平均.式(4)中规则空间导数项和引入平均导数优化的思想是为了充分利用全部9个网格点的信息，而混合空间导数项本身已经是一个致密的格式，必须利用全部9个网格点的信息才能表征(式(7))，因此在新格式中仍然保持与VTI常规9点一样的写法.据此构建得到的二阶VTI ADM 9点有限差分格式为：

其中：

采用经典的频散分析方法，引入一个平面波表示P(x，z，ω)=P₀e^－i(k_xx+k_zz)来进行频散分析的研究.

先考虑Δx≥Δz的情形，令r=Δx/Δz，将平面波P(x，z，ω)=P₀e^－i(k_xx+k_zz)代入VTI ADM 9点差分格式即式(8)可求取差分方程的相速度关系式，同时各向异性频率空间域波动方程频散关系的求取还需要特别考虑求取波动方程的相速度关系式(吴国忱，2006)，最后求得Δx≥Δz情形下的相速度频散关系式如下：

其中：

系数α，β，c，d由求E(α，β，c，d)的极小得到：

本文通过最小二乘优化方法来求取VTI ADM 9点格式的优化系数，求取的不同(Δx≥Δz)情况下的优化系数如表 1所示(以r=1，1.5，2，2.5，3，3.5，4为例).

表 1(Table 1)

表 1 当Δx≥Δz时对于不同求得的优化系数 Table 1 Optimization coefficients for different when Δx≥Δz 优化系数 r=1 r=1.5 r=2 r=2.5 r=3 r=3.5 r=4 α 0.805588317 0.817454805 0.820573933 0.812282362 0.797106514 0.781268301 0.766020269 β 0.772897959 0.805241484 0.814318890 0.810724448 0.806261295 0.804467906 0.804287094 c 0.628332546 0.647967013 0.641399953 0.620776768 0.596913581 0.572592034 0.546171078 d 0.094294523 0.087558376 0.092263981 0.101306788 0.111586150 0.122848081 0.135721278

表 1 当Δx≥Δz时对于不同求得的优化系数 Table 1 Optimization coefficients for different when Δx≥Δz

考虑到VTI ADM 9点格式系数在Δx≥Δz和Δx < Δz两种情况下的几何对称性，将Δx≥Δz情形下求取的优化系数α与β互换即可得到Δx<Δz情况下的优化系数，c与d保持不变(Chen，2012，2014).采用最小二乘优化方法求取的不同r=(Δx<Δz)情况下的优化系数如表 2所示(以r=1.5，2，2.5，3，3.5，4为例).

表 2(Table 2)

表 2 当Δx < Δz时对于不同求得的优化系数 Table 2 Optimization coefficients for different when Δx < Δz 优化系数 r=1.5 r=2 r=2.5 r=3 r=3.5 r=4 α 0.805241484 0.814318890 0.810724448 0.806261295 0.804467906 0.804287094 β 0.817454805 0.820573933 0.812282362 0.797106514 0.781268301 0.766020269 c 0.647967013 0.641399953 0.620776768 0.596913581 0.572592034 0.546171078 d 0.087558376 0.092263981 0.101306788 0.111586150 0.122848081 0.135721278

表 2 当Δx < Δz时对于不同求得的优化系数 Table 2 Optimization coefficients for different when Δx < Δz

图 2给出了VTI ADM 9点格式和VTI常规9点格式在不同r=ΔxΔz(Δx≥Δz)情况下根据表 1求取的优化系数计算得到的相速度频散曲线，其中令Thomsen各向异性参数ε=0.2，Δ=0.1.通过频散曲线分析得出，VTI常规9点格式当每个波长需要12个网格点的时候相速度误差才能控制在±1%以内，而达到同样的误差精度VTI ADM 9点格式只需要3.57个网格点数，VTI ADM 9点格式的频散误差明显小于VTI常规9点格式的频散误差，即VTI ADM 9点格式的计算精度明显提高.

图 2

Fig. 2

图 2 当Δx≥Δz时对于不同 VTI常规9点格式和VTI ADM 9点格式相速度频散曲线图Fig. 2 Phase velocity curves of the VTI conventional 9-point scheme and VTI ADM 9-point scheme for different when Δx≥Δz

引入频率空间域完全匹配层(Perfectly Matched Layer，PML)吸收边界条件(Berenger，1994；吴国忱等，2007；张衡等，2014；Moreira et al．，2014)，构建得到VTI ADM 9点格式PML波动方程(附录A)，并结合根据VTI模型不同的纵横向采样间隔求得的优化系数，即可实现VTI ADM 9点格式频率空间域数值模拟.下面通过一个复杂的BP2007 VTI海洋标准模型来分析VTI ADM 9点格式的正演精度和正演效率，VTI常规9点算法作为对比方法进行比较.

本文数值模拟时采用截取的部分BP2007 2D VTI标准模型(图 3).截取的BP2007 2D VTI模型模拟的主要是近海岸的地质构造，此模型特点是左部有一个被各向异性地层包围的高速各向同性盐丘，模型参数包括纵波速度V_P、Thomsen纵波各向异性参数ε、Thomsen变异系数Δ，模型区域大小为横向网格数N_x=401，纵向网格数N_z=361，模型横向网格间距为10 m，纵向网格间距为4 m，此时横纵向网格间距比达到了2.5，其中海水层为各向同性介质(图 3).

图 3

Fig. 3

图 3 截取的BP2007 2D VTI海洋标准模型参数 (a) VP； (b) ε； (c) δ.Fig. 3 Parameters of truncated BP2007 2D VTI ocean standard model

数值模拟时，雷克震源子波主频取15 Hz，震源位置设置在纵向第11层的横向300个网格点处，坐标为(3000 m，40 m)，检波点排列置于海平面以下40 m(即电缆沉放深度)处，为测线所有点接收，每个网格点设置一个检波器，全排列共401个检波器.采用频率空间域PML吸收边界条件，PML层数设为40层，传播时间为2 s，时间采样间隔为0.004 s.

值得一提的是，实际测试表明VTI频率空间域波场传播过程中经常会出现较强的伪横波(SV波)噪声，SV波干扰产生的原因是由于TI声波近似导致的.式(12)表示TI声波近似下SV波的相速度表达式，容易得出TI声波近似(直接令V_Sz=0)虽然使得沿对称轴方向(θ=0°)的横波速度为零，但是其他传播方向上的SV波相速度并不为零，从而产生了伪SV波噪声(Grechka et al．，2004).这种伪SV波噪声相比qP波而言虽然体现出低速低能量的特征，但是在实际偏移和反演中仍会产生比较严重的干扰，而且可能会导致算法不稳定，因此需要加以压制和消除.对式(12)的分析表明若把震源激发的位置设置为椭圆各向异性介质(ε=δ)或各向同性介质(ε=δ=0)即可消除伪SV波噪声.如通过采取在震源处添加一个各向同性薄层的策略即可有效实现频率空间域qP波波场传播过程中的伪SV波噪声消除.

因为本例震源置于各向同性的海水层中，VTI数值模拟时不会出现影响qP波场传播的伪横波(SV波)噪声，此种情况下不需要采取添加震源各向同性薄层消除伪SV波噪声的策略.

分别采用VTI常规9点格式和VTI ADM 9点格式进行频率空间域正演模拟，横纵向网格间距比为2.5情况下的VTI ADM 9点差分格式优化系数如表 1所示.从VTI ADM 9点格式频率空间域模拟35 Hz的单频波场切片(图 4a)可以看出频率空间域qP波场特征明显，其运动学特征得到了准确的刻画，在模型右侧速度渐变带，频率空间域波场变化不大，而在左侧陡峭盐丘处由于速度和各向异性参数变化复杂，入射波和反射波的相互干涉导致了反射界面上频率空间域波场的剧烈扰动.

图 4

Fig. 4

图 4 (a)VTI ADM 9点格式计算的35 Hz单频波场； (b)VTI常规9点格式计算得到的时间域地震记录； (c)VTI ADM 9点格式计算得到的时间域地震记录； (d)VTI介质时间域12阶高阶有限差分方法计算得到的时间域地震记录Fig. 4 (a) 35 Hz monochromatic wavefield computed by VTI ADM 9-point scheme； (b) Time-domain seismograms computed with the VTI conventional 9-point scheme； (c) Time-domain seismograms computed with the VTI ADM 9-point scheme； (d) Time-domain seismograms computed with the time domain 12-order high order finite difference scheme for VTI media

对所有频率的频率空间域波场进行反傅里叶变换到时间域就可得到所有时刻的时间域波场记录，提取每一时刻海面测线所有检波器采集得到的时间域波场记录即可得到时间域地震记录(图 4b，4c).从VTI常规9点格式的地震记录(图 4b)和VTI ADM 9点格式的地震记录(图 4c)对比可以清晰看出VTI ADM 9点的地震记录几乎没有频散，能精确模拟，而VTI常规9点的地震记录则出现了严重的数值频散，这是由于两种方法数值模拟精度差异较大而导致VTI常规9点格式在大网格间距下数值频散误差较大，而VTI ADM 9点格式仍能保持高精度.具体而言由于VTI常规9点格式达到无频散模拟所需的每个波长网格点数为12个，由于最大频率约为40 Hz且最小速度(即海水速度)为1492 m·s^-1，则无频散模拟最大允许网格间距大概为1492/40/12=3.1 m，而VTI ADM 9点格式达到无频散模拟所需的每个波长网格点数仅为3.57个，则无频散模拟最大允许网格间距可达到1492/40/3.57=10.4 m，因此对于本例横纵向网格间距分别为10 m和4 m的情况来说，横向网格间距10 m已经明显超过了VTI常规9点格式精确模拟所允许的3.1 m间距，因此产生了比较严重的数值频散，而VTI ADM 9格式仍能精确模拟.另外从地震记录上来看频率空间域PML吸收边界条件的引入很好地消除了人工边界反射(图 4b，4c).通过与图 4d中VTI介质高精度时间域12阶高阶有限差分方法(Duveneck and Bakker，2011)的结果进行对比，验证了本文方法与VTI介质时间域高阶方法精度相当.

对BP2007 2D VTI海洋标准模型相同计算区域同等计算精度下做了一个计算效率测试(本文同等计算精度即以恰恰不产生数值频散为限)，计算平台为联想ThinkCentre M8500t台式机电脑(酷睿i7八核，3.4 GHz).测试结果为VTI常规9点格式频率空间域单炮模拟耗时10262 s，而VTI ADM 9点格式频率空间域单炮模拟耗时仅为2730 s.因此在相同的精度要求下，VTI ADM 9点格式的计算效率要明显高于VTI常规9点格式，这是由于VTI ADM 9点格式相比VTI常规9点格式计算精度提升明显，因此可以通过采用更大的网格间距来提高计算效率.

就内存存储空间而言，如果VTI常规9点格式需要N_x×N_z(N_x和N_z分别为横向和纵向网格点数)个网格点来存储系数矩阵，则本文的VTI A

DM 9点格式只需要个网格点.因此VTI ADM 9点格式耗费的内存存储空间要明显少于VTI常规9点，从而减少了存储要求.

VTI介质频率空间域常规9点格式精度低，本文提出了基于平均导数优化方法的二维VTI介质qP波波动方程频率空间域二阶9点格式，通过优化方法使数值频散达到极小化，每个波长所需要的网格点数在1%的误差范围内仅为3.57个网格点数，而VTI常规9点差分格式在相同的误差范围内则需要约12个网格点数，因此VTI常规9点格式波场模拟精度差，相较本文方法更容易出现数值频散误差，本文方法体现出在克服数值频散误差方面的优越性.数值模拟结果表明本文的基于平均导数优化方法的二维VTI介质频率空间域正演方法不仅具有很好的计算精度和计算效率，而且平均导数优化方法的特性也使其同时具有很好的适用性和灵活性，能适用于不同的纵横向采样间隔情况.从新方法的应用前景上来看，由于各向异性频率空间域正演是各向异性频率空间域全波形反演的重要基础和关键环节，本文研究的基于平均导数优化方法的高精度VTI介质频率空间域正演方法可进一步应用于VTI介质频率空间域全波形反演的研究中，具有良好的应用前景.

引入PML技术(Berenger，1994；Moreira et al．，2014)，将频率域PML二维拉伸函数e_x和e_z与二维VTI介质频率空间域qP波波动方程(式(4))相结合推导得到二维VTI介质频率空间域qP波PML波动方程：

其中f_peak是震源主频，x_i和z_i分别为PML区域内部点与边界的x方向和z方向的距离，n_pml为PML层厚度，a₀取经验值1.79(吴国忱等，2007).

式(A1)结合二阶VTI ADM 9点有限差分格式(式(8))推导得到VTI ADM 9点格式的PML波动方程如下：

当波从内部区域过渡到PML边界层后波场逐渐衰减.