随着地震勘探对深层地层的深入，对复杂地区的保幅成像越来越受到处理人员的重视了.鉴于地层各向异性的广泛性，如何对各向异性介质进行保幅成像也逐渐成为地震勘探的热门问题.由于地震波在各向异性介质中一般会发生横波分裂、不同方向传播的地震波速度不同、纵横波耦合强烈等现象；灰岩页岩相间的薄层交互层、有固定排列方向的垂直裂隙以及沉积岩中的旋回性砂泥岩薄互层都会引发地震波的各向异性.针对地球介质的实际情况进行地震勘探和高精度成像已经越来越重要和艰难.叠前逆时偏移技术利用全波长成像原理，对高陡构造、横向变速以及棱柱波等成像有较强的优势，具有更高的成像精确和更丰富的有效信息，故本文将其引入到各向异性介质的成像中优势是十分明显.

各向异性介质的逆时偏移成像一般有两种方式：一是直接利用弹性波的控制方程作为地震波场正向、反向传播的算子，进而进行成像.Dellinger和Etgen(1990)，Yan和Sava等(2010)通过求解各向异性介质弹性波的Christoffel方程得到P波和S波的偏振方向，从而求得空间域波场分离算子.逆时偏移中利用分离出纯P波、纯S波进行偏移成像.从某种意义上说，利用这种弹性波方程更加接近实际地质体，但实现的过程中存在许多问题：一是多参数初始模型的建立，这个对野外数据采集要求太高，各向异性介质参数将更多，很难实现；二是弹性波控制方程描述的地震波场是耦合的，这对计算机的计算效率和内存要求加大；三是在成像过程中，每个空间采样点的地震波场正、反传的每个成像时间步长都必须使用空间域波场分离算子实现波场分离，这太消耗时间和计算机的存储空间.考虑到以上困难，在实际生产中，主要还是使用标量波方程进行偏移成像.

另一中方式就是各向异性下的拟声波近似方程，即采用声学近似的方式导出控制方程.Alkhalifah(1998)最先提出qP波的控制方程：在频散方程中令对称轴上的SV波速度为零，简化频散关系，得到qP波方程；还有Zhou等(2006)、Du等(2008)等都对频散关系作了相应的优化简化，分别导出各自的qP波方程，但这样的简化也就决定其SV波并不能完全的消除，同时方程必须在介质满足特定的条件下才能保持稳定；Zhang等(2011)等人对Du的拟声波方程进行改进，通过保持能量守恒的方式证明其方程的稳定性；Pestana等(2011)、Zhan等(2011)等人则是通过求解出解耦的P波和SV波的频散关系从而得到纯P波和纯SV波波动方程，由于方程需要用隐格式求解，而隐格式求解效率较低，用在计算量较大的逆时偏移中并不是很合适.

本文根据笔者导出的拟声波一阶应力-速度方程，提出一套稳定的各向异性介质叠前逆时偏移技术.为了避免由各向异性参数引起的地震波传播的不稳定现象，文中根据Fletcher等(2009)理论，引入横波项，消除由残留的横波的不稳定现象.由本文导出的拟声波方程为一阶偏微分方程组，能很好的适用交错网格差分形式.通过对VTI_HESS模型和逆冲模型的数值模拟和逆时偏移，表明该拟声波方程能准确的描述TI介质中地震波运动学特征，并能很好的实现各向异性介质逆时偏移.

各向异性介质的拟声波方程是研究该介质下地震波正演、偏移成像、反演等的重要依据，传统的拟声波方程都是由频散关系通过不同方式推导而来，其实它们之间可以相互转化，只是由于引入的辅助函数不同，才使得求解出来的波动方程有所差别.本文将基于弹性波的基本理论，即本构方程、运动微分方程还有几何方程，推导TI介质中的拟声波一阶速度-应力方程.

我们知道TTI介质可由VTI介质通过坐标变换而来，所以VTI介质可认为是TTI介质的一种特例.下面以VTI介质为例，在观测坐标系下，其微分运动方程可写为

由几何方程我们可得到应变和速度分量的关系为

(5)式是VTI介质中弹性波一阶速度-应力方程，下面对其进行声学近似，刚度系数与Thomsen参数关系为

(8)式就是所要求的拟声波一阶速度-应力方程(为了更清楚表示物理量，将下标由1、2、3改为x、y、z).较其传统的拟声波方程，该方程组为一阶形式，更容易求解；另一方面该方程组以应力表示更具有物理意义.

在三维观测系统内，设TI介质的对称轴和OXZ平面内的观测系统Z轴的夹角为θ，这个角度称为极化角；在OXY平面内与X轴的夹角为φ，称为方位角.这样只要对(8)进行相应的坐标变换，就能将VTI介质中的一阶方程转换到TTI介质中.设转换矩阵为C(详细的推导过程参考附件B)

将(9)式带入到方程(8)式中就可得到TTI介质中拟声波一阶速度-应力方程为

对二维均匀各向异性介质试验，各向异性参数：vp0=3000 m/s，ε=0.24，δ=0.1，其中TTI介质的XOZ平面内的旋转角度：θ=45°.本文新推导的qP波数值模拟效果如图 1所示，为了能够对比本文推导的一阶方程与基于频散关系推导的二阶方程差异，笔者选取近年来较为常用的Zhou等(2006)所推导的二阶方程做对比(其方程是由Alkhalifah的四阶方程降阶转变而来)，图 2为Zhou等(2006)qP波方程波场快照.

由模型一的波场快照图 1、2可知，当各向异性参数ε≥δ时，Zhou的方程和新推导的方程的qP波波前面呈椭圆状，符合各向异性介质中qP波的运动学特征；而且都没有出现频散等不稳定现象，只是横波干扰不一样，Zhou的方程的横波干扰明显强于新推导一阶方程，主要是由于二阶方程放大了横波干扰项.

各向异性介质中的拟声波方程主要由两种方式推导：一是由P-SV波的频散关系简化变换得到，另一就是像本文由弹性波的基本理论推导而来.针对以上两种方式的拟声波方程，笔者将从能量守恒的角度来分析其控制方程的稳定性.

在物理学中我们常把动能表示为，势能表示为Ep(t)=∫Fdh.其中势能即保守力做的功，在弹性波动力学中，其内力就是应力，也就是说其势能为质点在应力作用下产生位移所做的功，可以表示为Ep(t)=∫σdu(其中u是位移矢量，σ是应力矢量)，再利用位移和应变关系：，将该关系式带入势能表达式中得到势能密度为

由上述(11)、(12)式我们就能定义弹性波的能流密度为(单位质量的质点所携带的能量)

将本文所推导的拟声波一阶速度-应力方程组带入动能和势能表达式中，求取能流密度时间偏导(详细的推导过程请参考附件A，正文中笔者将直接写出计算结果)，经过简化计算可以得到动能和势能对时间导数的表达式为(VTI介质)

为了对比新推导的qP方程和传统的qP波方程数值模拟的稳定性及各向异性参数的适应性，对较为复杂的逆冲模型进行试算，将模型改为二维TTI介质模型，如图 3分别是逆冲模型的Vp0、ε、δ和θ模型，其中θ是指ZOX平面内的旋转角度.

图 4和图 5分别为新推导的qP方程和Zhou等(2006)推导的qP方程正演结果，由上图可知新推导的方程对逆冲模型的数值稳定性要好于Zhou等(2006)方程，特别是随着波传播的时间增大，差异更明显.这是由于Zhou、X.Du等人推导的TTI介质中的qP波方程为二阶方程，这使得横波残留波场也成二次方增加，特别是在TTI介质的各向异性参数变换较大的区域更为明显.

如图 6分别为Zhou方程与新qP方程RTM成像结果，对比可见，新方程时的成像结果都比Zhou方程成像结果更加清晰，对各向异性介质的刻画更加符合真实的模型，在高陡构造成像保幅效果更好，由此，体现了新方程的在TI介质逆时偏移成像方面的优越性，同时验证了该方程数值求解的稳定性.

为了验证逆时偏移中考虑各向异性可行性和必要性，利用本文推导一阶拟声波方程对HESS模型进行试验，如图 7所示，分别是HESS模型的纵波速度场，ε模型，δ模型，其中P波速度范围为1500 m/s到4500 m/s，ε和δ的范围分别为0到0.28和0到0.16，震源为主频25赫兹的雷克子波，总共300炮，炮间距：30 m，每炮排列数：600，检波器间隔：10 m，时间采样长度：4.1 s，时间采样间隔：0.4 ms，空间采样间隔：10 m.采用震源归一化互相关成像条件以及Laplace滤波.

如图 8为Zhou推导的qP方程和新导的qP方程正演结果(其中震源在模型中心位置)，由图可知在HESS模型中新导的方程和Zhou方程的qP波波前都能很好的保持运动学特征，但在各向异性参数ε、δ变化加大区域新导的方程数值稳定性要好于Zhou方程.

图 9分别为常规声波逆时偏移、拟声波逆时偏移成像结果.由图可知常规逆时偏移构造形态发生错乱，深部能量弱，成像效果差.VTI逆时偏移对该模型断层构造进行了很好的刻画，高陡倾角断层面处同相轴清晰，地下盐丘轮廓明显，尤其是大倾角地区，保幅成像效果更加明显.

本文基于弹性波的基本理论，推导出各向异性介质拟声波一阶速度-应力方程.文中从能量守恒的角度出发，证明了新推导的方程比基于频散关系推导出的拟声波方程更加稳定且能量保持守恒.通过对均匀二维模型和逆冲模型进行正演模拟，表明新推导的方程不仅能很好的维持qP波的运动学特征，而且能够很好的适应于任意旋转角度的TTI模型，而其他控制方程在都会出现不同程度的数值不稳定现象.本文将该各向异性介质拟声波一阶速度-应力方程引用到叠前逆时偏移技术中，引入归一化互相关成像条件.模型试算表明，各向异性介质拟声波一阶速度-应力方程逆时偏移能够校正介质的各向异性影响，成像剖面更加清晰，深部能量加强，整体能量更均衡，成像效果更加符合真实地层的构造特征.

由文中的定义，我们定义弹性波的能流密度为(单位质量的质点所携带的能量)

将(8)式应力应变关系改写为

对(A-4)求时间导数：

将(A-6)带入(A-5)中：

对上式求分部积分：

由(A-10)和(A-8)可知，现在我们证明了本文推导的一阶速度-应力方程组在传播过程中能量是守恒.

基于频散关系推导拟声波方程主要有：Alkhalifah(1998)、Zhou等(2006)、Du等(2008)和Fowler(2009)的方程.下面于较为常用的Du等(2008)方程为例.

同理可得：

对上式求时间导数：