地球物理学报  2013, Vol. 56 Issue (10): 3474-3486   PDF    
各向异性介质qP波传播描述I:伪纯模式波动方程
程玖兵 , 康玮 , 王腾飞     
同济大学海洋地质国家重点实验室, 上海 200092
摘要: 地球介质相对于地震波波长尺度的定向非均匀性会导致波速的各向异性,进而影响地震波场的运动学与动力学特征.各向异性弹性波动方程是描述该类介质波场传播的基本工具,在正演模拟、偏移成像与参数反演中起着关键作用.为了面向实际应用构建灵活、简便的各向异性波场传播算子,人们一直在寻求简化的各向异性波动方程.本文借鉴各向异性弹性波波型分离思想,通过对平面波形式的弹性波方程(即Christoffel方程)实施一种代表向波矢量方向投影的相似变换,推导出了一种适应任意各向异性介质、运动学上与原始弹性波方程完全等价,在动力学上突出qP波的新方程,即qP波伪纯模式波动方程.文中以横向各向同性(TI)介质为例,给出了相应的qP波伪纯模式波动方程及其声学与各向同性近似,并在此基础上开展了正演模拟和逆时偏移试验,展示了这种描述各向异性波场传播的新方程的特点与优势.
关键词: 各向异性      波动方程      相似变换      伪纯波      逆时偏移     
Description of qP-wave propagation in anisotropic media, Part Ⅰ:Pseudo-pure-mode wave equations
CHENG Jiu-Bing, KANG Wei, WANG Teng-Fei     
State Key Laboratory of Marine Geology, Tongji University, Shanghai 200092, China
Abstract: The ordered heterogeneity much smaller than the wave length will results in velocity anisotropy and thus affects the kinematic and dynamic properties of the seismic waves. As the tools to describe the wave propagation in earth′s media, anisotropic elastic wave equations play key roles in seismic modeling, migration and inversion. To generate flexible and efficient wave propagators in anisotropic media for practical us age, people never give up to construct simplified anisotropic wave equations in the past years.In this paper, based on wave-modes eparation theory, we propose so called pseudo-pure P-wave equations which are equivalent to the original elastic wave equation in kinematics but highlight the P-wave energy.We derive the pseudo-pure P-wave equations in TI media and their further simplified forms.The examples of synthesized wavefront snapshots and reverse time migration prove that, the pseudo-pure-mode wave equation is accurate anisotropic wave propagator with distinguishing features, and has great potential in seismic modeling, imaging and parameter in version..
Key words: Anisotropy      Waveequation      SimilarityTransform      Pseudo-purewavemode      Reverse     
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方程

设一般各向异性弹性介质的刚度张量为cijklijkl=1,2,3),则线性应力-应变关系(广义胡克定律)满足

(1)

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

(2)

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

(3)

其中uii=1,2,3)代表粒子振动位移分量.根据牛顿运动定律,存在如下运动方程:

(4)

其中ρ为介质的密度,fi为外力,üiui的二阶时间微分.(3)式代入(1)式,再把(1)式代入(4)式,可得位移表征的三维三分量各向异性弹性波动方程:

(5)

为了简化表示,记,于是(3)表示成矩阵形式:

(6)

若把偏导数矩阵记为L,波场位移矢量记为u=[uxuyuz],外力记为F,则波动方程(5)也写成矩阵形式

(7)

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

(8)

且有:

(9)

对无源的各向异性弹性波方程进行平面波分析,(7)式变换成Christoffel方程:

(10)

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

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

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

(14)

其中变换矩阵MP满足

(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)

其中C11C33C44C66C13为五个独立参数.为了保证刚度矩阵正定,这些独立刚度系数满足如下约束条件[36]C33>0,C44>0,C66>0,(C11-C66C33C132.于是前文Γ矩阵的元素简化成

(22)

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

(23)

将(23)式转入频率波数域,可知其Christoffel矩阵满足

(24)

根据(14)式对应的相似变换,新的Christoffel矩阵满足

(25)

将(25)式代入(17)式,反变换到时空域就得到VTI介质qP波伪纯模式波动方程:

(26)

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

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

(27a)

(27b)

(27c)

(27d)

(27e)

(27f)

其中vpzvsz代表qP与qS波垂向传播速度,vpn代表层内NMO速度,vpxvsx代表qP波与qSV波的横向传播速度,δεγ为其他Thomsen参数.于是方程(26)改写成

(28)

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

(29)

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

(30)

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

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

例如,在(29)式中取vsz=0,就得到声学近似qP波伪纯模式波动方程

(31)

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

(32)

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

(33)

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

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

(34)

其中

(35)

很容易推导出如下关系:

(36)

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

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

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

图 1 二维弱各向异性VTI介质波场模拟 (a)和(b)为原始弹性波方程位移场的水平与垂直分量;(c)和(d)为qP波伪纯模式波动方程模拟波场的水平与垂直分量,(e)为伪纯波场两个分量之和. Fig. 1 Wavefield modeling in a 2D weak VTI medium (a) horizontal and (b) vertical components synthesized by using the original elastic wave equation; (c) horizontal and (d) vertical components synthesized by using pseudo-pure qP-wave equation; (e) summation of (c) and (d).
图 2 二维强各向异性TTI介质波场模拟 (a)和(b)为原始弹性波方程位移场的水平与垂直分量;(c)和(d)为qP波伪纯模式波动方程模拟波场的水平与垂直分量,(e)为伪纯波场两个分量之和. Fig. 2 Wavefield modeling in a 2D strong TTI medium (a) horizontal and (b) vertical components synthesized by using the original elastic wave equation; (c) horizontal and (d) vertical components synthesized by using pseudo-pure qP-wave equation; (e) summation of (c) and (d).

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

图 3 二维VTI介质伪纯模式波动方程波前模拟与运动学分析 图中前两行对应伪纯波场的水平与垂直分量,第三行两个分量之和;第四行对应由群速度计算的qP波与qSV波的理论波前曲线;从左至右对应σ=∞,σ=6.0,σ=1.5和σ=0.75四种情况. Fig. 3 Wavefield modeling and kinematic analysis in a 2D VTI medium The firs ttwo rows represent horizontal and vertical components of pseudo-pure qP-wave fields; the third row represents the summation of the two components, namely the scalar pseudo-pure P-wave fields; the last row represents the theoretical curves of the qP and qSV wavefronts calculated based on group velocity.From left to right, the parameter σ changes from ∞ to 0.75.
5.2 三维波场模拟

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

图 4 三维均匀VTI介质伪纯模式波动方程正演模拟 前三行对应伪纯波场的三个分量,第四行为分量和波场.左、中、右三幅图分别对应过震源的iLine、xLine和水平切片. Fig. 4 Wavefiel dmodeling in a 3D VTI medium From left to right are iLine, xLine and depth slices of the snapshots.The first three row srepresent two horizontal and one vertical component of pseudo-pure qP-wave fields; the fourth row represent the summation of three components.
5.3 叠前逆时偏移

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

图 5 SEG/HESSVTI模型 (a)qP波垂直速度;(b)和(c)为Thomsen参数εδ Fig. 5 SEG、Hess VTI model (a) vertical qP-wave velocity; (b) and (c) are the two Thomsen parameters ε and δ.
图 6 伪纯qP波波前快照 (a)、(b)和(c)分别为伪纯qP波场的水平分量、垂直分量及其求和. Fig. 6 Snapshots of pseudo-pure qP-wave fields (a), (b) and (c) are the wavefield components and their summation, respectively.
图 7 基于qP波伪纯模式波动方程的逆时偏移结果 Fig. 7 Reverse-time migration result based on pseudo-pure qP-wave equation
6 结论

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

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

致谢

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

参考文献
[1] Tatham R H, M D McCormack.Multicomponent seismology in petroleum exploration.Investigations in Geophysics, vol. 6, SEG.1991.
[2] 何樵登, 张中杰. 横向各向同性介质中地震波及其数值模拟. 长春: 吉林大学出版社, 1996 . He Q D, Zhang Z J. Seismic Waves and Their Modeling in Transversely Isotropic Media (in Chinese). Changchun: Press of Jilin University, 1996 .
[3] 杜启振, 秦童. 横向各向同性介质弹性波多分量叠前逆时偏移. 地球物理学报 , 2009, 52(3): 801–807. Du Q Z, Qin T. Multicomponent prestack reverse-time migration of elastic waves intransverseisotropic medium. Chinese J. Geophys (in Chinese) , 2009, 52(3): 801-807.
[4] 张美根, 王妙月, 李小凡, 等. 各向异性弹性波场的有限元数值模拟. 地球物理学进展 , 2002, 17(3): 384–389. Zhang M G, Wang M Y, LiX F, et al. Finite element forward modeling of anisotropic elastic waves. Progress in Geophysics (in Chinese) , 2002, 17(3): 384-389.
[5] Zhang Q S, McMechan G A. 2D and 3D elastic wavefield vector decomposition in the wavenumber domain for VTI media. Geophysics , 2010, 75(3): D13-D26. DOI:10.1190/1.3431045
[6] Yan J, Sava P. Improving the efficiency of elastic wave-mode separation for heterogeneous tilted transverse isotropic media. Geophysics , 2011, 76(4): T65-T78. DOI:10.1190/1.3581360
[7] Yan J, Sava P. Elastic wave-mode separation for VTI media. Geophysics , 2009, 74(5): WB19-WB32. DOI:10.1190/1.3184014
[8] Alkhalifah T. An acoustic wave equation for anisotropic media. Geophysics , 2000, 65(4): 1239-1250. DOI:10.1190/1.1444815
[9] Zhang J F, Verschuur D J, Wapenaar C P A. Depth migration of shot recordsin heterogeneous, transversely isotropic media using optimum explicit operators. Geophysical Prospecting , 2001, 49(3): 287-299. DOI:10.1046/j.1365-2478.2001.00255.x
[10] Han Q Y, Wu R S. Aone-waydual-domain propagator for scalar qP-waves in VTI media. Geophysics , 2005, 70(2): D9-D17. DOI:10.1190/1.1884826
[11] Shan G J.Optimized implicit finite-difference migration for VTI media. ∥ 76th Ann.Internat. Mtg., Soc.Expl. Geophys., Expanded Abstracts, 2006:2367-2370.
[12] 刘礼农, 高红伟, 刘洪, 等. 三维VTI介质中波动方程深度偏移的最优分裂Fourier方法. 地球物理学报 , 2005, 48(2): 406–414. LiuL N, Gao H W, Liu H, et al. Wave equation depth migration with optimum split-step Fourier method in 3D VIT media. Chinese J. Geophys (in Chinese) , 2005, 48(2): 406-414. DOI:10.1002/cjg2.v48.2
[13] 吴国忱, 梁锴. VTI介质qP波方程高精度有限差分算子. 地球物理学进展 , 2007, 22(3): 896–904. Wu G C, Liang K. High precision finite-difference operators for qP wave equation in VTI media. Progress in Geophysics (in Chinese) , 2007, 22(3): 896-904.
[14] 朱兆林.倾斜横向各向同性介质中P-S波处理方法研究[博士论文].上海:同济大学, 2008. Zhu Z L.2008, The Processing Method of P-S Wave in Tilted Transversely Isotropic Media[Ph.D.Thesis] (in Chinese).Shanghai :Tongji University, 2008. http://www.geophy.cn/CN/abstract/abstract9787.shtml
[15] Zhou H, Zhang G, Bloor R.An anisotropic acoustic wave equation for VTI media. ∥ 68th EAGE Conference & Exhibition.Extended Abstracts, H033.2006. http://www.oalib.com/references/19003715
[16] Fletcher R P, Du X, Fowler P J.A new pseudo-acoustic wave equation forTI media.∥ 78th AnnualInternational Meeting.SEG, Expanded Abstracts, 2008:2082-2086. http://manu39.magtech.com.cn/Geophy/EN/abstract/abstract9787.shtml
[17] Zhou H, Zhang G Q, Bloor R.An anisotropic acoustic wave equation for modeling and migration in 2D TTI media.∥ 76th Annual International Meeting.SEG, Expanded Abstracts, 2006:194-198. http://www.oalib.com/references/18987281
[18] 黄翼坚, 朱光明, 刘池洋. 一个近似的VTI介质声波方程. 地球物理学报 , 2011, 54(8): 2117–2123. Huang Y J, Zhu G M, Liu C Y. An approximate acoustic wave equation for VTI media. Chinese J. Geophys (in Chinese) , 2011, 54(8): 2117-2123.
[19] Fowler P J, Du X, Fletcher R P. Coupled equations for reverse time migration in transversely isotropic media. Geophysics , 2010, 75(1): S11-S22. DOI:10.1190/1.3294572
[20] 康玮, 程玖兵. 横向各向同性介质拟声波方程及其在逆时偏移中的应用. 地球物理学报 , 2012, 55(3): 1033–1045. Kang W, Cheng J B. Pseudo-acoustic wave equations for reverse-timemigrationinTImedia. Chinese J. Geophys (in Chinese) , 2012, 55(3): 1033-1045.
[21] 段鹏飞, 程玖兵, 陈三平, 等. TI介质局部角度域射线追踪与叠前深度偏移成像. 地球物理学报 , 2013, 56(1): 269–279. Duan P F, Cheng J B, Chen S P, et al. Local angle-domain ray-tracing and prestack depth migration in TI medium. Chinese J. Geophys (in Chinese) , 2013, 56(1): 269-279.
[22] KlieH, W Toro.A new acoustic wave equation for modeling in anisotropic media.∥ 71stAnnualInternationalMeeting. SEG, ExpandedAbstracts, 2001:1171-1174. http://www.geophy.cn/CN/abstract/abstract9787.shtml
[23] Grechka V, Zhang L B, RectorJ W. Shear waves in acoustic anisotropic media. Geophysics , 2004, 69(2): 576-582. DOI:10.1190/1.1707077
[24] Duveneck E, Bakker P M. Stable P-wave modeling for reverse-time migration in tilted TI media. Geophysics , 2011, 76(2): S65-S75. DOI:10.1190/1.3533964
[25] Zhang Y, Zhang H Z, Zhang G Q. A stable TTI reverse time migration and it simple mentation. Geophysics , 2011, 76(3): WA3-WA11. DOI:10.1190/1.3554411
[26] 李博, 李敏, 刘红伟, 等. TTI介质有限差分逆时偏移的稳定性探讨. 地球物理学报 , 2012, 55(4): 1366–1375. Li B, Li M, Liu H W, et al. Stability of reverse time migration in TTI media. Chinese J. Geophys (in Chinese) , 2012, 55(4): 1366-1375.
[27] Liu Y, Sen M K. Acoustic VTI modeling with a time-space domain dispersion-relation-based finite-difference scheme. Geophysics , 2010, 75(3): A11-A17. DOI:10.1190/1.3374477
[28] Song X L, Fomel S. Fourier finite-difference wave propagation. Geophysics , 2011, 76(5): T123-T129. DOI:10.1190/geo2010-0287.1
[29] Kang W, Cheng J B.New coupled equation sofP-andSV-wave for RTM in TImedia.∥ 73rdEAGEAnnualMeeting, 2011, P267. http://manu39.magtech.com.cn/Geophy/EN/abstract/abstract9787.shtml
[30] Kang W, Cheng J B.Prestack scalar reverse time migration of elastic seismic datain TI media.∥ Expanded Abstracts of 81st Ann.Internat. Mtg., Soc.Expl.Geophys., 2011: 3367-3371. http://www.oalib.com/references/18992077
[31] Dellinger J, Etgen J. Wave-field separation in two-dimensional anisotropic media. Geophysics , 1990, 55(7): 914-919. DOI:10.1190/1.1442906
[32] Du X, Fletcher R P, Fowler P J.A new pseudo-acoustic wave equation for VTI media. ∥ 70th Conference and Exhibition.EAGE, ExtendedAbstracts, H033.2008. http://manu39.magtech.com.cn/Geophy/EN/abstract/abstract9787.shtml
[33] Thomsen L. Weak elastic anisotropy. Geophysics , 1986, 51(10): 1954-1966. DOI:10.1190/1.1442051
[34] Tsvankin I D, Chesnokov E M. Synthesis of body-wave seismograms from point sources in anisotropic media. Journal of Geophysics Research , 1990, 95(B7): 11317-11331. DOI:10.1029/JB095iB07p11317
[35] Psencik I, Gajewski D. Polarization, phase velocity, and NMO velocity of qP-waves in arbitrary weakly anisotropic media. Geophysics , 1998, 63(5): 1754-1766. DOI:10.1190/1.1444470
[36] Tsvankin I. Seismic Signatures and Analysisof Reflection Datain Anisotropic Media. Pergamon: Elsevier, 2001 .
[37] Fletcher R P, Du X, Fowler P J. Reverse time migration in tilted transversely isotropic (TTI) media. Geophysics , 2009, 74(60): WCA179-WCA187.