随着宽方位地震勘探技术的发展及勘探精度的不断提高，各向异性地震成像受到业界更多的关注.逆时偏移方法基于精确的双程波动理论，在时间方向上对震源波场进行正向延拓、检波点波场反向延拓(Baysal et al.，1983；McMechan，1983；Whitmore，1983)，正反向延拓的波场零相位相关产生成像结果.该方法具有无反射倾角限制、能适应强横向变速，能够对回转波及多次反射波进行正确成像的特点.在各向同性介质中，可以对标量声波方程进行求解，得到逆时偏移成像结果.严格来说实际接收的地震波场更接近于弹性波方程所描述的波场特征.因此，一些学者考虑能否直接用弹性波方程外推成像(Sun et al.，1986；Chang and McMechan，1986).但利用弹性波方程进行波场外推成像存在很多困难，首先是多参数建模困难，其次是直接外推耦合的方程组计算量巨大，另外多分量数据波场分解及相关成像难度很大.因此，目前开展弹性波逆时偏移成像实用化的条件还不成熟.但并不能否认弹性波逆时偏移成像的重要性.

对于各向异性介质，严格意义上的标量波是不存在的.但对于一般的偏移成像，与各向同性介质相似，可通过声学近似的方式，能够找到一个能代表 TI介质中P波运动学特征的控制方程.利用这样的方程就能进行P波波场的外推和偏移成像.Alkhalifah(1998)首先提出了拟声波近似思想，通过声学近似简化耦合的频散关系，引入了TI介质拟声波近似频散关系并导出方便数值求解的qP波控制方程.尽管这个标量波场的频散关系具有和矢量波场中P波相似的运动学特征，但它会产生伪横波.对于VTI介质，基于Alkhalifah的频散关系式，许多学者已经推导出了拟声波双程波动方程并用于VTI介质逆时深度偏移成像中，Alkhalifah(2000)，Zhou等(2006)，Du等(2008)，Fletcher等(2008)，Duveneck等(2008)利用不同辅助变量简化了拟声波方程的求解.同时在二维情况下通过解耦的P-SV关系导出‘纯P’波方程，但‘纯P波’方程的快速数值算法仍不成熟(Liu et al.，2009；Chu et al.，2011).

从Hooke定律和运动方程出发，通过声学近似也可推导出VTI介质拟声波方程.从平面波传播的角度，VTI介质向TTI介质转化可以看成在旋转后的坐标系下的平面波传播.坐标系旋转后，仅增加了数学运算的复杂度，并不增加物理量纲.通过设置对称轴方向SV波速度为零，可推导出TTI介质下的拟声波方程(Fletcher et al.，2009；Duveneck and Bakker，2011).

本文将着重介绍一个新的TTI介质拟声波方程推导及实现策略，该方程从精确P-SV波频散关系推导而来，并且使用Thomsen(1986)定义的各向异性参数ε和δ进行参数化.对于对称轴方向不变且ε≥δ的介质(VTI或HTI介质)来说，通过声学近似，可得到稳定的数值解，可以直接将该方程应用于正演和偏移中去.然而，对于非均匀TTI介质中的正演和逆时偏移，各向异性对称轴方向的变化会引起数值计算的不稳定问题.为了得到更加稳定的TTI介质中qP波方程，Flecter等(2009)通过在声学近似方程中加入非零横波速度得到有限横波方程来增加方程求解稳定性.更自然的解决办法由Duveneck和Bakker(2011)首先提出，即对弹性波方程直接做声学近似，避开从频散关系推导qP波方程时对称轴方向缓变的假设，导出了TTI介质中与弹性波方程对应的qP波方程，该方程可以稳定模拟P波分量的传播.Zhang Y和Zhang H Z(2009)从构造自共轭微分算子角度导出了一种TTI介质中的qP波方程.该方程可以看作是Duveneck和Bakker(2011)导出的qP波方程的近似方程，其忽略了部分对称轴方向空间导数项以提高计算效率，但保留了Duveneck方程中微分算子的自共轭性来保证稳定性(Zhang Y and Zhang H Z，2009；Zhang H Z and Zhang Y，2011).本文采用一种新的方式增强qP波传播时的稳定性，采用的是有限横波正则化的方法，这样就可以消除由倾斜对称轴变化以及由ε<δ引起的传播不稳定问题，同时最大程度地减弱波场传播中的伪横波影响.在数值试算及实际资料试算部分，展示了采用所导出方程的正演和逆时偏移结果，证明了方法理论的正确性.

从平面波传播的角度，从VTI介质向TTI介质转化可看作是在旋转坐标系下的平面波传播.首先，基于Christoffel方程从相速度/频散关系出发，可导出VTI介质耦合频散关系：

式中(k_x，k_y，k_z)表示波数向量，v_Pz表示垂向P波速度，表示P波动校正速度，v_Px=表示横向P波速度，v_Sz表示横向SV波速度，ε和δ是Thomsen各向异性参数.(1)式的时间域四阶偏微分方程很难直接求解，由于涉及到的时间层及储存的波场较多，不利于数值求解(Alkhalifah，2000).因此，对频散关系(1)式做了一系列的变量替换得到更利于求解的二阶偏微分方程.在(1)式基础上两边同时乘上p(ω，k_x，k_y，k_z)，并引入辅助变量：

则(1)式可以重新写为：

联立(2)、(3)式，并将其反变换到时空域，得到VTI介质中的有限横波方程(Zhou et al.，2006)：

同理，通过求解非均匀TTI介质的Christoffel方程，可推导出旋转坐标系下TTI介质P-SV波耦合频散关系:

式中(k_x′，k_y′，k_z′)表示与观测坐标对应的旋转坐标系下的波数，ω表示角频率，v_Pz表示垂直于对称面方向的P波速度，表示相对于垂直对称面方向的P波动校正速度，v_Px=表示对称面内的P波速度，v_Sz表示垂直于对称面的SV波速度，ε和δ是Thomsen(1986)定义的各向异性参数.

由观测坐标与旋转坐标系下波数矢量的关系可得：

在常角度假设前提下，将(6)式反变换到空间域，可得到旋转坐标系下的二阶微分算子：

将(4)式中的二阶微分算子换成(7)式旋转坐标系下的形式，即可得TTI介质中的有限横波方程：

将(7)式对应的空间微分算子用H₁、H₂重写可得：

则有限横波方程(8)可表示为：

式中，θ是倾角，Φ是对称轴的方位角.注意对于VTI介质，微分算子H₁只作用于垂直方向，而H₂只作用于水平方向.两个算子都不包含空间交叉导数项.在各向同性介质中，式(10)中的两个二阶偏微分方程是等价的，波场p等价于波场q.

方程(10)可以采用伪谱法求解，求解中要用到波数域空间导数计算(Fletcher，2009).然而，由于伪谱法需要多次傅里叶变换，不利于大型数据的波传播并行计算.因此，采用时空域高阶有限差分法求解(10)式.尽管波场p和q都是四阶偏微分方程组的解，但是我们将波场p看作压力场，将q看作辅助波场时，式(10)可看作是两个二阶偏微分方程，与各向同性声波方程类似，可通过高阶差分方法求解(董良国，2000；刘红伟，2010；李博，2010；刘文卿，2012，2013)，本文不再讨论其解法.

方程(10)即为Fletcher等(2009)提出的所谓的有限横波方程计算(即剪切波速度不为0).从方程(10)出发，根据Alkhalifah的思想，将对称轴方向剪切波速度设为0，可得到一个更为简化的TTI介质拟声波方程，这个方程更容易求解，但数值求解也更加不稳定(Fletcher et al.，2009)，Zhang Y和Zhang H Z(2009)采用伪谱法计算时，也发现了类似的不稳定现象.当然，我们可以对模型做平滑处理能使波场稳定传播，但这种方式同时会显著地改变波的运动学特征，影响成像精度.值得注意的是，有限横波方程也并非完全稳定，这种不稳定性通常开始于模型中对称轴倾角存在剧烈变化的地方.本文提出了一种带正则化项的有限横波方法，数值试验表明相较于有限横波方程，这种方程能使波场更加精确稳定地传播下去.下面讨论带正则化项的有限横波方程实现过程.

所谓带正则化项方程，即是设计一个滤波器，对波场传播中不稳定的高波数成分做适当衰减，以达到减少频散或使波场计算稳定的目的.首先来具体推导各向同性介质中带正则化项平面波解的形式，以形象说明正则化项的作用.对于构建的带正则化项的方程：

方程(11)两边同时对时间/空间做傅氏变化，经过整理并带入平面波解通解形式，可得到各向同性介质下带正则化项声波方程平面波解形式：

其中 k =(k_x，k_y，k_z)， x =(x，y，z)分别为波数以及空间坐标向量，σ为正则化系数(通常取值(0~1)). 与不带正则化项的声波方程平面波解相比，(12)式多出了一项振幅衰减项，其物理意义为，对平面波波场中传播时间较长的高波数成分做衰减.利用这一思想，假设波场中不稳定成分主要为高波数情况下，构造TTI介质中带正则化项的有限横波方程，使得波场数值计算更加稳定.对于方程(10)可以构建成带正则化项的有限横波方程：

带正则化项的方程(13)看起来似乎要比原始的有限横波方程(10)复杂许多，但在实际计算中，其增加的计算量和储存量并不多.原因在于正则化项中旋转坐标系下波场微分无需重复计算，仅需多储存两个波场，记录上一时刻所有正则化项波场之和对应的两个波场，然后对时间求微分即可.实验表明选择合理的σ值，即使倾斜同相轴的方位及倾角变化较大，依然能够保证波场稳定地传播.这也是这种方式实际资料计算中比较可行的原因.

利用带正则化校正的TTI介质qP波方程进行逆时偏移成像主要采用以下步骤实现：首先求解炮波场的沿时间的正向外推定解问题，如式(14)—(15).

其初始条件为：

其中 s 代表震源位置，f(t)为模拟震源函数.在此RTM实现方案中，首先将边界层最里侧的波场记录下来.二维情况下是四个边界面处的波场记录，三维情况下是六个边界面处的三维波场.视机器内存大小，可存储在内存中，也可存储在硬盘上.

其次，在求解炮波场正向外推基础上再求解波场沿时间逆向外推定解问题，如式(16)—(17).

其边值条件为：

其中Bp(x，t)、Bq(x，t)分别是两个波场在边界处的值.利用每个时间层上的边界波场，解上述边值问题，沿时间逆时外推可实现炮波场逆向外推.此时，与检波点波场的逆向外推相同步.因此在每个时间步，利用成像条件进行炮波场与检波点波场零延迟相关成像.

为了测试采用带正则化校正的TTI介质qP波方程波场传播的稳定性，我们选择了Foothill TTI模型中参数变化最为剧烈的部分作为测试模型.图 1展示了该模型的参数数据，模型中包含了速度、倾角及各向异性参数.图 2展示的是不同控制方程不同时刻的波场快照，图 2a在正演方程中采取不带正则化项的有限横波方程，在模型中存在对称轴倾角剧烈变化的地方，波场传播不稳定.模型正演中采用25 Hz的雷克子波作为震源，震源放在(3550，0)m处，计算网格在两个方向上均为10 m.利用空间10阶、时间2阶有限差分去逼近控制方程中相应导数项.图 2c、2d分别为采用带正则化校正的TTI介质qP波方程2000 ms、3000 ms的波场快照，与图 2a、2b同一时刻的波场快照对比可以看出，此时波的传播是稳定的.通过数值试算可看到，直接由频散关系导出的TTI介质中的有限横波qP波方程在忽略参数的空间的导数时，由于受对称轴倾角变化的影响仍然会产生数值求解不稳定的现象，我们通过正则化的方式衰减高波数成分，最终使得整个波场计算稳定.

图 1

Fig. 1

图 1 Foothill模型参数 (a) vP速度； (b) 倾角模型； (c) epsilon模型； (d) delta模型.Fig. 1 Model parameters of Foothill (a) P-wave velocity vP; (b) Angle θ defining symmetry axis direction; (c) Thomsen parameter ε; (d) Thomsen parameter δ.

图 2

Fig. 2

图 2 不同方程传播的波场快照 (a) 有限横波方程2 s波场； (b) 有限横波方程3 s波场； (c) 正则化有限横波方程2 s波场； (d) 正则化有限横波方程3 s波场.Fig. 2 Wave field snapshots from (a) Finite shear wave equation at t=2 s; (b) Finite shear wave equation at t=3 s; (c) Regularized finite shear wave equation at t=2 s; (d) Regularized finite shear wave equation at t=3 s.

我们将带正则化校正的TTI介质qP波方程方法应用到TTI逆时偏移中，Foothill模型作为该测试所用模型，模型参数如图 1所示，波场传播的计算均是采用高阶有限差分法来求解，我们采用的是时间方向2阶差分，空间方向10阶差分的有限差分算法.图 3展示了逆时偏移成像结果，反射层都正确归位，陡倾角地层都能成像.但在地层倾角较大、各向异性较强的区域，其下覆地层影响较大，各向同性及VTI各向异性逆时偏移成像效果都不理想，而TTI各向异性逆时偏移能使陡倾角地层及下覆地层准确成像.因此，对于各向异性介质的地震数据，采用TTI各向异性算法进行成像可得到更为精确的结果.

图 3

Fig. 3

图 3 Foothill 模型不同方法逆时偏移结果对比 (a) 各向同性逆时偏移结果； (b) VTI各向异性逆时偏移结果； (c) TTI各向异性逆时偏移结果.Fig. 3 RTM results of Foothill model using different methods (a) Isotropic migration result; (b) VTI migration result; (c) TTI migration result.

本文所采用的实际数据是我国西部某油田的资料.该区发育大规模岩溶缝洞型储层，非均质性及各向异性都比较明显.以往由于各向异性的影响，人们认为窄方位角资料效果好，如今方位各向异性现象逐步被人们认识，它可以提供更多的裂缝信息和储层信息.因此，常规各向同性地震成像方法由于没有考虑地层各向异性特性，很难准确刻画缝洞的空间位置.传统的偏移距道集是反射点所有信息的综合，并没有准确包含倾角信息和方位角信息，无法将角度信息从道集中提取出来、也无法精确描述与方位有关的各向异性(王西文，2010).因此，只有全方位的角度道集才能获得准确的全方位速度信息、倾角信息及与方位有关的各向异性参数信息，可以基于全方位共反射角道集进行各向异性速度建模，建模流程如图 4.

图 4

Fig. 4

图 4 全方位角度域网格层析各向异性速度建模流程Fig. 4 Work flow of full-azimuth angle-domain grid tomography method and anisotropic velocity model building

图 5为传统Kirchhoff方法获得的共反射点(CRP)道集与全方位反射角道集对比.从道集特征可看出，常规的CRP道集是反射点所有信息的综合，掩盖了不同方位的速度变化特征及各向异性特征.全方位反射角道集定义了一个新的数据结构(Koren and Ravve，2011)，其角度为5°、10°、15°、20°、25°、30°、35°，每一个角度包含全方位的信息.从道集特征可看出，随反射角度增大，相同反射角，不同方位道集不平直，具有不同的速度，其各向异性特征更加明显.利用全方位反射角道集不同方位的剩余时差(图 5中绿线)可获得有效的各向异性参数，通过网格层析成像方法建立精确的各向异性速度模型，如图 6.

图 5

Fig. 5

图 5 全方位反射角道集(a)与传统CRP道集(b)对比Fig. 5 Comparison of full-azimuth reflection angle gathers (a) and traditional CRP gathers (b)

图 6

Fig. 6

图 6 各向异性参数剖面 (a) epslion参数剖面； (b) delta参数剖面； (c)层速度参数剖面； (d) 地层倾角参数剖面； (e)方位角参数剖面.Fig. 6 Profiles of anisotropic parameters (a) Thomsen parameter ε; (b) Thomsen parameter δ; (c) Interval velocity; (d) Dip parameter; (e) Azimuth angle.

图 7a为各向同性逆时偏移成像结果，图 7c为TTI各向异性逆时偏移成像结果，由于TTI各向异性逆时偏移成像考虑整个波场信息、地层的各向异性参数及倾角和方位角等空间属性参数，所以串珠成像更聚焦、深度误差更小，成像精度更高.从成像结果及钻井误差对比分析得出，各向同性逆时偏移结果深度误差为190 m(图 7a)，而各向异性逆时偏移结果深度误差为30 m(图 7b，图 7c)，对于TTI各向异性逆时偏移成像由于上覆地层倾角不大，垂向消除深度误差的效果与VTI介质逆时偏移基本相同，但在横向上串珠位置相差12.5 m，与钻井更加吻合.同时深层成像TTI介质各向异性逆时偏移成像相对VTI介质及各向同性介质的成像有较大的改善.

图 7

Fig. 7

图 7 各向同性逆时偏移(a)、VTI各向异性逆时偏移(b)与TTI各向异性逆时偏移(c)对比Fig. 7 Comparison of isotropic RTM result (a), VTI RTM result (b), and (c) TTI RTM result

对于复杂构造成像，叠前逆时偏移是一个有力的工具.宽方位数据采集的普及使得各向异性逆时偏移变得必不可少.本文着重讨论了TI介质下的声学近似、基于频散关系的控制方程推导及高阶有限差分数值解法.基于Christoffel方程从相速度/频散关系出发推导出的qP方程，通过强制方程中沿对称轴方向的剪切波速度为0，可以得到一个简化的波动方程.这个退化的方程可以提高求解的效率，并且在ε－δ≥0的VTI介质中是稳定的.若将其推广到TTI介质，在速度模型复杂时存在数值计算的不稳定问题.我们采用一个新的TTI介质拟声波方程及其解决方案，在有限横波基础上提出了一种带正则化项的有限横波方法，保证波场能精确稳定地传播下去，同时最大程度地减弱SV反射系数的影响.

Foothill TTI模型的正演和逆时偏移试验表明，在波场传播步数较多时，常规的有限横波方程也开始出现计算不稳定的现象，通过正则化的方式衰减高波数成分，最终使得整个波场计算稳定，同时成像效果也最好.我们提出的方案在3D TTI介质实际资料中也能得到非常好的效果.在模型参数剧烈变化的区域，上述方法仍然会产生SV波.为了进行P波偏移，这些SV波被看作是噪声.但在实际资料偏移中，我们并没有在最终偏移的剖面上发现来自这些SV能量的污染，可能是SV能量相对于P波能量比较弱的缘故，这些SV波能量是否在RTM角道集上有反映，还有待于进一步研究.从具有实际物理含义及数值计算稳定角度出发，我们应在弹性波框架下做近似，即TI介质中qP控制方程的推导应在时空域进行，这些是下一步研究的方向.