1 引言

大量的研究与观察表明，地球介质普遍存在弹性各向异性，主要表现就是地震波传播速度随方向变化.在地震勘探中，人们主要关注的各向异性类型有^[1]：1）具有垂直对称轴的横向各向同性（VTI），其诱导因素是岩石的层状沉积结构（如由近似定向排列粘土板构成的页岩以及周期性薄层沉积）.2）具有倾斜对称轴的横向各向同性（TTI），如倾斜的页岩地层或薄互层沉积.3）具有水平对称轴的横向各向同性（HTI），诱导因素包括近似定向排列的垂直裂缝系统以及非均衡水平应力等.4）具有三个相互正交的镜像对称平面的正交各向异性（ORT），适合描述定向排列垂直裂缝系统的层状沉积岩石，或者是各向同性背景岩石中含两套相互正交的垂直裂缝系统的情况.如果出现倾斜沉积地层或倾斜裂缝系统，就表现为更复杂的各向异性类型.在弹性各向同性介质中，存在两种模式或类型的波，即纵波与横波.而在弹性各向异性介质中，除了纵波之外还存在快、慢横波.通常情况下，各向异性介质中纵波的偏振（或极化）方向并不与传播方向平行，横波的偏振方向也并不与传播方向垂直，因此人们把其中的纵波称为qP波，把两种横波分别称为qS1波、qS2波，或在横向各向同性（TI）介质中以qSV波和SH波来代表.

在描述各向异性介质中地震波的传播过程时，若采用原始的各向异性弹性波动方程，依赖的参数多，算法复杂，计算成本目前很难承受^[2-4].而且，各向异性弹性波波型分离在实际应用中面临诸多困难^[5-7].因此，近十年人们针对P波信号为主的垂直分量地震资料的正演模拟与偏移成像主要采用简化的各向异性介质qP波方程，而且主要考虑VTI和TTI介质两种情况.推导VTI介质qP波方程最常用的方法就是采用“声学近似”，即将Thomsen参数表征的耦合的qP-qSV波频散关系或波动方程中沿对称轴的横波速度设为零.Alkhalifah^[8]首先从声学近似频散关系推出了简化的时空域qP波（四阶）波动方程以及分裂后二阶方程（组）.紧接着许多学者基于声学近似频散关系推导出各种VTI和TTI介质的单程波传播算子^[9-14].而后一些学者仍然从声学近似频散关系出发推导出了更容易求解的耦合二阶拟声波方程^[15-18].理论上，将一个高阶线性偏微分方程分解成低阶偏微分方程组是不唯一的.Fowler等^[19]指出，VTI介质四阶频散方程可分解成一系列等价的耦合二阶频散方程，而且这些矩阵形式表示的耦合方程组之间还可以借助相似变换相互转换.他们还在文中列举了几种所谓“最小的”耦合方程.尽管这些不同的耦合方程在运动学上等价的，但在非均匀介质中动力学特性存在差别，波场变量的量纲及其能量强弱也各有不同^[19-20].为了考虑倾斜层状各向异性，通过坐标系旋转可将VTI介质qP波拟声波方程扩展到TTI介质^[16-17].此外，段鹏飞等^[21]基于qP波拟声波方程推导出的射线追踪方程，提出了高效的TI介质局部角度域叠前深度偏移成像方法.

声学近似是对各向异性介质波场传播的一种很有效的简化，极大地促进了TI介质逆时偏移（RTM）技术的发展与应用.不过，声学近似同时也引起了一些棘手的问题.首先，声学近似只能消除沿对称轴传播的qSV波，其他方向残余的横波干扰会给纵波波场造成污染^[22-23].其次，声学近似改变了VTI介质原有的弹性性质，只有在Thomsen参数满足一定条件（即ε≥δ）时才能保证拟声波方程的数值稳定性^[23].对于TTI介质，当对称轴方向空间变化剧烈时，声学近似引起的数值不稳定性问题非常突出.近年来，Duveneck等^[24]、Zhang和Zhang ^[25]等从第一原理（即运动方程与广义胡克定律）出发，借助自共轭约束等数学技巧构建出了相对稳定的qP波拟声波方程.康玮和程玖兵^[20]从VTI介质弹性波方程出发推导了波场分量为正应力的qP-qSV波耦合方程及其在TTI介质中的扩展形式，讨论了耦合方程及其拟声波方程在计算成本、稳定性等方面的差异.李博等^[26]也通过数值试验讨论了TTI介质拟声波方程RTM算法的稳定性问题.TI介质拟声波方程的数值解法也是许多学者关注的课题^[27-28].

为了推导不含声学近似的各向异性波动方程，Kang和Cheng ^[29-30]从VTI介质二维三分量弹性波方程出发，借助具有波场方向投影意义的相似变换，分别推导出了动力学上qP波占优的或qSV波占优的耦合波动方程.这种新的各向异性波动方程具有三方面特点：其一，它们在形式上与Fowler^[19]提到的“最小的”二阶耦合方程相似，也不存在混合偏导数，数值计算比原始各向异性波动方程方便.其二，它们所依赖的相似变换具有明确物理含义，在波数域，变换后的矢量波场所有分量之和与原始弹性波场向各种波模式偏振方向的各向同性参考方向的标量投影是等价的.例如，向qP波各向同性参考方向投影获得的耦合波动方程运动学上与原始波动方程完全等价，波场分量求和得到的标量波场在动力学上突显qP波能量，但残余较弱的qS波.根据上述特点，我们将这些单一波模式占优的耦合方程称为伪纯模式波动方程.

下文从一般各向异性介质三维三分量弹性波方程出发，基于Christoffel方程的特殊相似变换，推导qP波伪纯模式波动方程，给出在VTI和TTI介质中的具体形式及其拟声波近似，然后通过波场正演模拟与RTM试验，展示伪纯模式波动方程的特点与应用潜力.

2 各向异性介质弹性波动方程与Christoffel方程

设一般各向异性弹性介质的刚度张量为c_ijkl（i，j，k，l=1，2，3），则线性应力-应变关系（广义胡克定律）满足

(1)

其中σ_ij为应力张量，ε_kl为应变张量.上式中采用了爱因斯坦求和符号.采用Voigt符号记法，上式还可写成矩阵形式：

(2)

应变与位移之间满足如下柯西方程：

(3)

其中u_i（i=1，2，3）代表粒子振动位移分量.根据牛顿运动定律，存在如下运动方程：

(4)

其中ρ为介质的密度，f_i为外力，ü_i为u_i的二阶时间微分.（3）式代入（1）式，再把（1）式代入（4）式，可得位移表征的三维三分量各向异性弹性波动方程：

(5)

为了简化表示，记和，于是（3）表示成矩阵形式：

(6)

若把偏导数矩阵记为L，波场位移矢量记为u=[u_x，u_y，u_z]，外力记为F，则波动方程（5）也写成矩阵形式

(7)

其中Γ为3×3的微分算子对称矩阵，满足

(8)

且有：

(9)

对无源的各向异性弹性波方程进行平面波分析，（7）式变换成Christoffel方程：

(10)

其中ω代表频率，是位移波场的傅里叶域形式，就是所谓的Christoffel矩阵，代表空间偏导数矩阵L的波数域等价形式，它由波数k_x，k_y和k_z确定.

Christoffel方程对应的频散关系描述了波场的运动学特征.对于固定的平面波成分，（10）式与经典的特征值问题对应.Christoffel矩阵的三个特征值分别对应三种波模式或波型（qP、qS1和qS2）的相速度，与每个特征值对应的特征向量代表该波型的粒子偏振矢量.由于Cristoffel矩阵是实对称的，它一定存在三个相互正交的特征向量，因此三种波型的偏振方向相互正交^[31].将（10）式反映的频散关系按矩阵特征值系统进行分析，允许我们通过相似变换获得许多与原始各向异性弹性波动方程等价的耦合二阶线性方程（组）^[19].事实上，这些年针对VTI介质出现了很多种不同的耦合波动方程^{[15, 16, 32]}. Fowler等^[19]从相似变换的角度指出这些耦合方程在运动学上是等价的，但对偏移成像而言只有那些形式简洁的“最小的”耦合方程才有利于数值计算.由于这些方程都是从数学角度推导出来的，波场变量的物理含义往往不直观，所模拟的波场既不代表qP波，也不代表qS波，而是二者的某种耦合形式.对于常规纵波资料RTM，人们常常取qP波能量相对较强的分量波场施加成像条件.为了消除qSV波假象，通常还会在这些耦合方程基础上进一步施加声学近似条件，或者直接从声学近似频散关系出发进行相似变换获得声学近似耦合方程.当然，这就必然遇到引言提到的声学近似引起的问题.

3 基于相似变换构建qP波伪纯模式波动方程

通常情况下，各向异性介质中qP波的偏振方向并不与波场传播方向平行，qS波的偏振方向也不和波场传播方向垂直，只有严格区分偏振方向才能实现qP与qSV波的解耦.准确计算偏振方向需要求解Christoffel方程^[31].然而，在非均匀各向异性介质中，偏振方向随局部弹性参数变化而改变，qP与qSV波分离涉及非常复杂的数值计算.也就是说，基于偏振方向差异直接推导只传播单一波模式的波动方程是不现实的.本文的思路源自于对弹性波模式分离的经典算法的观察.

在各向同性介质中，弹性波场的散度对应标量纵波数据，其旋度对应矢量形式的横波数据.例如，空间域弹性波场的散度可写成

(11)

它在傅立叶域或波数域的等价形式满足

(12)

其中代表散度算子∇·在傅立叶域的对应算子.上式还可表示成矢量内积形式：

(13)

其中k=[k_x，k_y，k_z]为波矢量， 代表与波矢量k方向一致的单位矢量.傅立叶域每点（k_x，k_y，k_z）处的波场代表沿方向以空间频率传播的平面波.给定平面波传播方向，通常对应两种传播速度和偏振方向不同的波型或波模式，即P波与S波，其中前者的偏振方向与传播方向平行，后者的偏振方向与传播方向正交.我们知道，两个矢量的内积物理上代表一个矢量向另一个矢量方向的标量投影.因此Yan和Sava^[6]将（13）式代表的弹性波模式分离解释为弹性波场在波矢量方向的标量投影.注意，该式还存在一个与空间频率有关的因子ik.严格地讲散度算子或其傅立叶域等价算子不单纯是一个波模式分离算子，它还会放大波场的高波数成分并引起90°相位移^[5-6].尽管如此，它们在实现各向同性弹性波中P波和S波数据分离方面是非常有效的，故一直在多分量地震数据处理中被广泛采用.

根据上述观察，我们将波矢量方向视为各向异性qP波偏振方向的各向同性参考方向.为了把傅立叶域的原始弹性波场向该方向进行标量投影，本文引入如下针对Christoffel矩阵的相似变换：

(14)

其中变换矩阵M_P满足

(15)

于是，与此变换对应的新波场为

(16)

满足新的Christoffel方程，即：

(17)

将其反变换到时空域就得到一种新的二阶耦合波动方程：

(18)

其中与相似变换前的微分算子矩阵Γ相对应.

在运动学上，方程（7）、（10）、（17）与（18）完全等价，都准确描述了各向异性介质三种波模式的运动学特征.不难发现，新的波场在傅立叶域满足：

(19)

与（12）式对比可知，本文引入的相似变换与新波场的分量求和，共同达到将傅立叶域原始弹性波场向波矢量方向标量投影的效果.相应地，空间域的新波场分量求和对应原始弹性波场u的散度，即：

(20)

对于各向同性介质，这样获得的标量波场u仅含P波，不含S波.但是，就新的矢量波场而言，它同时包含P波与S波两种成分.在各向异性介质中，qP波的偏振方向与波矢量方向不一致，由相似变换与波场分量求和引起的标量投影只达到部分波型分离效果，最终的标量波场u中会同时包含qP与qS波成分.

对大多数地下岩石，即使存在弹性各向异性，各向异性都不是很强^[33].在绝大多数平面波方向上，qP波的偏振方向与波矢量的方向偏差并不大^[34-35].由于qS波偏振矢量总是与qP波偏振矢量垂直，故通常比qP波更偏离波矢量方向.因此可以推断，本文这种新的标量波场u中qP波能量是占优的，qS波能量相对弱一些（后文数值试验也证实了这一推断）.也就是说，（18）式的运动学完全正确，具有部分波型分离功能（即动力学上突出标量qP波的成分），故本文称之为qP波伪纯模式波动方程.

针对典型的各向异性介质，只需将（18）式中Christoffel矩阵代入相应的刚度系数即可.下文以最常见的横向各向同性介质为例，先给出VTI介质中的qP波伪纯模式波动方程及其各种近似，然后讨论如何扩展到对称轴倾斜的TTI介质.

4 VTI介质qP波伪纯模式波动方程

代表VTI介质弹性性质的刚度系数矩阵可写成

(21)

其中C₁₁、C₃₃、C₄₄、C₆₆和C₁₃为五个独立参数.为了保证刚度矩阵正定，这些独立刚度系数满足如下约束条件^[36]：C₃₃＞0，C₄₄＞0，C₆₆＞0，（C₁₁-C₆₆）C₃₃＞C₁₃².于是前文Γ矩阵的元素简化成

(22)

进而写出VTI介质三维三分量弹性波动方程：

(23)

将（23）式转入频率波数域，可知其Christoffel矩阵满足

(24)

根据（14）式对应的相似变换，新的Christoffel矩阵满足

(25)

将（25）式代入（17）式，反变换到时空域就得到VTI介质qP波伪纯模式波动方程：

(26)

它在运动学上与原始弹性波方程（23）是完全等价的，且不含有混合偏导数项（故而有利于数值计算）.

若用Thomsen参数来表征VTI介质的刚度系数，即：

(27a)

(27b)

(27c)

(27d)

(27e)

(27f)

其中v_pz和v_sz代表qP与qS波垂向传播速度，v_pn代表层内NMO速度，v_px与v_sx代表qP波与qSV波的横向传播速度，δ、ε和γ为其他Thomsen参数.于是方程（26）改写成

(28)

对于构造成像，可以利用VTI介质x-y平面各向同性偏振特征，将（28）式中前两个方程合并成一个方程，从而得到数值计算更经济的波动方程

(29)

其中代表伪纯波两个水平分量之和.在二维情况下，上式简化为

(30)

注意，方程（29）形式上与Fowler^[19]推导的“最小的”耦合方程很相似，计算代价也相当.不过，本文采用的相似变换具有明确的物理含义，数值上也能够达到突显标量qP波能量的目的.由于常规地震勘探所接收的信号以纵波能量为主，因此本文qP波伪纯模式波动方程的应用价值是不言而喻的.

对于以纵波成分为主的常规地震资料偏移成像处理而言，为了压制残余横波造成的假象，文中借鉴前人常用的两种方法.一是把震源附近小区域视为各向同性或椭圆各向异性介质.二是对伪纯模式波动方程进一步施加经典的声学近似，既可简化方程，降低计算成本，也可大幅度压制残余横波假象.

例如，在（29）式中取v_sz=0，就得到声学近似qP波伪纯模式波动方程

(31)

注意，方程（31）与Fletcher等^[16]推导的VTI介质拟声波方程形式相近，计算成本相当.不过，这里伪纯波场分量之和主要代表qP波能量.二维情况下，上式简化成

(32)

对于各向同性情况（ε=δ=0），（31）式简化后进一步合并为一个方程，即

(33)

其中v_P代表各向同性P波的传播速度，代表三个波场分量之和，与原始弹性波场的散度对应，（34）式恰好就是经典的常密度声波方程.这也进一步证实本文采用的相似变换获得的耦合方程具有明确的物理含义.

对于TTI介质，对称轴方向不再是垂直方向.但从物理上讲，TTI介质和VTI介质没有本质差异，只不过观测坐标系与描述介质对称性的物理坐标系不一致而已.通过坐标旋转，可由VTI介质波动方程推导出TTI介质中的对应形式^[20].设对称轴方位角为φ，倾角为θ，则物理坐标系x′=（x′，y′，z′）与实际观测坐标x=（x，y，z）系满足如下变换关系：

(34)

其中

(35)

很容易推导出如下关系：

(36)

将（26）或（28）式及其近似方程中空间偏导数项用（36）式中的关系代换，就能够推导出TTI介质qP波伪纯模式波动方程.注意，坐标旋转引入了一些混合偏导数项，增加了计算的复杂度.对于背景TI介质中含有定向垂直裂缝系统的ORT介质而言，坐标旋转还要考虑裂缝走向.

5 数值试验 5.1 二维波场模拟

首先，通过原始各向异性弹性波方程、伪纯模式波动方程正演模拟的波前快照，观察伪纯模式波动方程的特点.图 1显示了弱VTI介质（v_pz=3000.0 m/s，v_sz=1500.0 m/s，ε=0.1和δ=0.05）中弹性波场的水平与垂直分量（如图 1a、b），以及qP波伪纯波场的水平与垂直分量及其求和（如图 1c、d和e）.图 2显示了强TTI介质（v_pz=3000.0 m/s，v_sz=1500.0 m/s，ε=0.25，δ=-0.25和θ=30°）中的波场快照.可见，伪纯波波场的水平与垂直分量之和以qP波能量占优.随着各向异性增强，残余qSV波波前几何形状越来越复杂，其能量也变强，但仍然比qP波弱得多.

然后，利用本文伪纯模式波动方程正演模拟四种各向异性参数组合情况下的波前快照，观察运动学特征与残余波型情况（图 3）.根据Tsvankin^[36]的分析，TI介质qSV波运动学特征受参数σ=（ε-δ）v_pz²/v_sz²控制.对于均匀VTI介质，利用本文方程（30）模拟计算σ=∞，σ=6.0，σ=1.5和σ=0.75四种情况下的波前快照.图 3第一行分别为四种参数情况下波前快照的水平分量，第二行为相应波前快照的垂直分量，第三行为相应波前快照水平分量与垂直分量之和.可见，在分量求和之后的波场中以qP波能量为主，但残余较弱的qSV波能量.图 3第四行对应四种参数情况下根据群速度计算的几何波前曲线.从模拟波前运动学特征与理论曲线的一致性看出，本文伪纯模式波动方程在运动学上能够准确描述qP波与qSV波.此外，可以观察到当σ比较小时，qSV波波前形状比较规则；当σ较大时，qSV波波前存在多重路径现象.

5.2 三维波场模拟

本例进一步通过三维波前快照展示VTI介质伪纯模式波动方程的特点.模型参数为v_pz=3300 m/s，v_sz=1770 m/s，ε=0.2，δ=0.1和γ=-0.2.为了突出横波能量，这里γ的绝对值设得比较大，表示具有很强的横波各向异性.从正演模拟的某时刻的三维波场中，截取显示其中一些过震源（图片中心点）的波前快照（如图 4）.根据波速判定，最外面的波前对应qP波，最里面的波前对应SH波，中间的波前对应qSV波.从伪纯波的三个分量波场来看，当横波各向异性较强时，伪纯波波场中qSV波与SH波能量也比较强.由于VTI介质中SH波是在各向同性平面（即水平面）内偏振的，因此垂直分量中没有SH波能量（如图 4g、h、i）.不过，与二维情况类似，将伪纯波波场三个分量求和后，qP波能量得到突显，SH波完全消除，但还是残余一定的qSV波能量（如图 4j、k、l）.

5.3 叠前逆时偏移

利用从SEG官方网站共享数据服务器下载的Hess二维VTI模型与合成炮记录，测试本文伪纯模式波动方程在RTM中的应用效果.图 5显示了该模型qP波垂直速度与两个Thomsen参数.注意，该数据体中缺少qSV波垂直速度信息.这里采用Fletcher等^[37]提出的有限（非零）qSV波速度方法确定伪纯模式波动方程需要的横波速度信息.图 6为qP波伪纯模式波动方程正演模拟的1.1秒时波场的水平、垂直分量，以及二者之和（震源位于计算区域中心）.从分量求和得到的标量qP波伪纯波波前快照看，波场中残余qSV波较弱，而且传播距离也较短（因为横波速度较小）.图 7显示了多炮偏移叠加得到的成像剖面，所有反射界面包括断层和盐丘边界都获得了清晰的图像，也合理地描述了与盐体右侧斜交的楔状结构.此例表明，虽然基于qP波伪纯模式波动方程正向与逆时传播的波场中含有少量残余qSV波，但是多炮偏移叠加过程可以有效地干涉抵消这部分能量造成的假象.

6 结论

根据各向异性介质地震波传播模拟与偏移成像的需要，本文借鉴弹性波波型分离思想，对原始弹性各向异性波动方程对应的Christoffel矩阵施加特殊的相似变换，推导出了适合任意各向异性介质的标量qP波能量占优的伪纯模式波动方程.虽然该方程形式上与前人推导的“最小”耦合方程形式非常接近（计算代价相当），但其波场物理含义更明确.伪纯波场在运动学上与原始弹性波场完全等价，其分量求和获得的标量波场在动力学上以qP波能量为主.这对常规垂直分量为主的单分量地震资料偏移成像来说是非常有利的.与前人提出的拟声波近似方程不同，伪纯模式波动方程能够正确反映各种波型的运动学特征，故在多分量勘探中转换波正演、成像与速度反演等应用中同样具有一些价值.针对伪纯模式波动方程的声学近似可进一步简化方程、降低计算成本、压制横波干扰，但在TTI介质对称轴方向剧烈变化时可能会存在数值计算不稳定问题.对于构造成像而言，偏移速度模型以低频分量为主，相对光滑的各向异性速度模型有利于提高数值计算稳定性.与前人研究主要关注TI介质和纵波不同，本文推导伪纯模式波动方程的思路既适合任意各向异性介质，也适合横波，因此在今后多波地震勘探中也具有应用潜力.

本文qP波伪纯模式波动方程隐含向其波矢量方向（即各向同性P波的偏振方向）投影波场的特点，但由于各向同性与各向异性偏振方向的差异，导致qP波伪纯波波场中尚残余qS波能量.依据各向异性波型分离理论，可以根据投影方向偏差对伪纯波波场进行校正，进而提取出完全分离的qP波数据.这方面的工作将在本文姊妹篇中详细论述.

致谢

感谢HESS公司与SEG提供二维VTI模型及其合成地震数据.