地球物理学报  2017, Vol. 60 Issue (1): 258-270   PDF    
改进的TTI介质纯P波方程正演模拟与逆时偏移
郭成锋1 , 杜启振1 , 张明强2 , 郑华灿3 , 韩冬1     
1. 中国石油大学(华东)地球科学与技术学院, 青岛 266580;
2. 中海油田服务股份有限公司, 天津 300451;
3. 东方地球物理公司新兴物探开发处, 河北涿州 072751
摘要: 逆时偏移作为一种高精度偏移方法已成为复杂构造成像的重要技术,描述纵波独立传播的延拓方程是各向异性介质逆时偏移的一个关键问题.在对VTI介质几个经典相速度近似公式回顾的基础上,针对常用于描述纯P波的Harlan近似公式在各向异性参数ε较大情况下近似精度较低的问题,本文对Harlan公式中的非椭圆项进行了修正,在非椭圆项前添加了一个与各向异性参数ε有关的修正系数,得到了三种改进型Harlan公式,并以近似精度最高的改进式为基础,推导了TTI介质纯P波方程.针对该伪微分方程,本文利用伪谱法和有限差分法联合实现波场延拓,对于常密度二阶方程,基于中心网格实现;对于一阶应力-速度方程则基于旋转交错网格实现.通过数值试验分析了TTI介质纯P波一阶应力-速度方程的近似精度,并以一阶纯P波方程为基础进行了TTI介质逆时偏移数值模拟试验.结果表明,本文给出的方法能够较准确地描述TTI介质纯P波波场特征,可以应用至各向异性介质逆时偏移.
关键词: TTI介质      纯P波      混合法      正演模拟      逆时偏移     
Numerical simulation and reverse time migration using an improved pure P-wave equation in tilted transversely isotropic media
GUO Cheng-Feng1, DU Qi-Zhen1, ZHANG Ming-Qiang2, ZHENG Hua-Can3, HAN Dong1     
1. School of Geoscience, China University of Petroleum(East China), Qingdao 266580, China;
2. China Oilfield Services Limited, Tianjin 300451, China;
3. Development and Production Seismic Division, Hebei Zhuozhou 072751, China
Abstract: Seismic anisotropy has been observed in many exploration areas, especially in sedimentary rocks. The most commonly used anisotropic model is the transverse isotropy (TI) with a vertical (VTI) or a titled (TTI) symmetric axis. For practical application, transverse isotropy is commonly considered under the acoustic approximation. In general, such acoustic approximation can be divided into two categories:coupled acoustic anisotropic wave equations and decoupled pure P-wave equations. One of the problems of the coupled acoustic equations is that they are not really free of shear waves. The artificial SV component is usually considered as numerical artifacts and may cause numerical instabilities in TTI media. However, decoupled pure P-wave equations are completely free from shear wave. In this paper, we firstly evaluate some commonly used phase velocity approximations and conclude that the Harlan's approximation, widely used to derive pure P-wave equations, has low accuracy in a medium with large ε. Then we present a method to improve its accuracy by modifying the anelliptic (fractional) term and obtain three modifications. Based on the most accurate one, we derive an improved TTI pure P-wave equation. The resulting wave equation can be formulated either as a second-order equation or as a system of first-order equations. The proposed TTI pure P-wave equation involves pseudo-differential operators which are difficult to be calculated by the finite-difference method. We use a hybrid pseudo-spectral/finite-difference (PS/FD) scheme to solve the TTI pure P-wave equation. The pseudo-differential operators are calculated by the PS method and the differential operators are calculated by the high-order rotated-staggered-grid FD method (RSGFD). We test the hybrid PS/RSGFD scheme with various kinds of TTI models. In the first numerical example, we compare the first-order TTI pure P-wave equation with the counterpart elastic equation. We conclude that pure P-wave equation can provide good kinematic and dynamic approximations to the elastic wave equation for homogeneous media with weak-to-moderate anisotropy. The second test model is a heterogeneous wedge TTI model which contains sharp contrasts. This numerical example demonstrates such a hybrid PS/RSGFD scheme can provide a stable P-wave propagator. To further verify the proposed algorithm, we implement an anisotropic reverse-time migration (RTM) example on a small region of BP2007 TTI model. These numerical examples demonstrate the potential of the proposed method..
Key words: TTI media      Pure P wave      Hybrid method      Modeling      Reverse time migration     
1 引言

由于具有成像精度高且无倾角限制、适应任意复杂速度模型等诸多优点,逆时偏移(Whitmore,1983McMechan,1983Baysal et al.,1983)已经逐步成为复杂构造区域成像的重要方法(Zhang et al.,2011).对于各向异性介质,若采用基于各向同性假设的逆时偏移算法进行成像,则会产生构造假象,降低成像精度(Huang et al.,2009杜启振和秦童,2009),而各向异性介质全弹性波逆时偏移则存在计算量大、波场分离算法复杂、各向异性参数估计困难等问题,实施起来仍然比较困难.考虑到实际资料多以纵波数据为主,目前针对各向异性介质逆时偏移的研究主要是基于声学近似开展的.对地震勘探来说,最具有实际意义的是横向各向同性介质(TI).如果地层的对称轴沿着垂直方向,一般采用具有垂直对称的横向各向同性(VTI)模型来近似地球介质.当地层受到构造运动或者应力作用,其对称轴将发生倾斜时,具有任意倾斜对称轴的横向各向同性(TTI)模型则更符合实际情况.

Tsvankin和Thomsen(1994)Alkhalifah(1998)研究发现,横波速度对纵波的运动学特征的影响甚微.随后,Alkhalifah(2000)以VTI介质为例,通过将对称轴方向上的横波速度设置为零,提出了著名的各向异性声学近似假设.尽管声学近似在物理上不可实现,但是它却能够较准确地描述纵波的运动学特征.随后许多学者对各向异性声学近似进行了深入研究,得到了一系列形式不同的VTI和TTI介质声学方程.Zhou等(2006a2006b)通过引入辅助变量将Alkhalifah(2000)提出的四阶偏微分方程降阶为更容易计算的二阶方程组.Duveneck等(2008)则从胡克定律和运动方程出发,推导了另一种形式的各向异性声学近似方程,该组方程的优点是波场变量具有明确的物理意义,并且十分便于数值计算.Fowler等(2010)对VTI介质二阶声波方程的一般形式进行了深入探讨,证明了这类偏微分方程在运动学上的等价性.程玖兵等(2013)推导了一种适应任意各向异性介质的伪纯模式波动方程,该方程在运动学上与原弹性波方程完全等价,在动力学上突出P波.

尽管声学近似极大地促进了各向异性介质地震波传播与成像技术的发展和应用,然而它也有突出的问题.首先,声学近似并不能把横波完全去除,残留的横波波前面呈现为菱形,这种低速度、低振幅的菱形横波会给纵波波场造成干扰(Grechka et al.,2004).尽管一些滤波方法(Zhang et al.,2009; 张岩和吴国忱,2013)在一定程度上可以压制菱形横波干扰,但是仍无法将其彻底消除.其次,声学近似对介质的各向异性是有一定要求的,只有在各向异性参数ε≥δ时才能保证声学方程的数值稳定性(Alkhalifah,2000; 李博等,2012).另外,对于TTI介质,当对称轴的方向发生剧烈变化时,声学近似也会导致数值不稳定.针对这个问题,Fletcher等(2009)通过引入一个有限的横波速度,在一定程度上缓解了TTI介质中声学近似的数值不稳定.Zhang等(2011)则通过引入自共轭偏微分算子构建了一种相对稳定的TTI介质声学近似波动方程.Duveneck和Bakker(2011)直接从胡克定律出发,通过对应力和应变旋转亦推导出了一种稳定的TTI介质声学近似方程,但是计算过于繁琐.

由于声学假设存在菱形横波干扰、稳定性条件要求高的问题,近些年来,各向异性纯P波方程受到广泛关注.基于Harlan(1995)给出的一种VTI介质纯P波近似频散关系,Etgen和Brandsberg-Dahl(2009)采用伪解析法进行了正演模拟;Crawley等(2010)进而将伪解析法拓展到TTI介质,成功实现了TTI介质逆时偏移技术.Liu等(2009)Alkhalifah(2000)给出的VTI介质四阶方程出发得到了分离的纯P和SV波方程.Pestana等(2011)则利用泰勒近似对频散公式进行解耦,进而推导出VTI介质纯P波波动方程;Zhan等(2012)进而利用旋转波数的概念将其拓展到TTI介质中,采用快速扩展法(REM,Pestana和Stoffa,2010)对方程进行了求解.Zhan等(2013)为了提高计算效率,减少伪谱法中的傅里叶变换的次数,采用有限差分和伪谱法混合求解纯P波方程.杜启振等(2015)提出了VTI介质的纯P波一阶应力-速度方程交错网格混合法.

多数VTI和TTI介质纯P波方程逆时偏移方法(例如Crawley et al.,2010; Pestana et al.,2011; Zhan et al.,20122013)是基于Harlan(1995)近似公式实现的,但是Harlan相速度近似式在各向异性参数ε较大的情况下精度比较低.针对这个问题,本文利用与参数ε有关的替换速度对Harlan公式中的非椭圆项进行了修正,得到了三种改进型Harlan公式,并以近似精度最高的改进式为基础,推导了一种 TTI介质纯P波方程.

尽管纯P波方程无横波干扰,但是其形式上要比声学近似波方程复杂,纯P波方程涉及到分数形式的伪微分算子,而这类算子无法直接利用显式有限差分进行求解.虽然伪谱法可以方便地处理这类算子,但是所需的计算量较大.本文采用伪谱法和有限差分混合法(Zhan et al.,2013; 杜启振等,2015)实现波场延拓,对于常密度二阶方程,基于中心网格实现;对于一阶应力-速度方程则基于旋转交错网格实现.数值试验证明了旋转交错网格有限差分混合法在求解纯P波一阶应力-速度方程中的可行性,利用具有强分界面的非均匀楔形模型验证了混合法纯P波一阶方程的数值稳定性,并实现了BP2007 TTI模型(部分)的正演模拟和逆时偏移.

2 TTI介质纯P波方程

从精确的VTI介质相速度公式出发,利用一阶泰勒近似得到P波近似频散公式,再通过坐标旋转的方式得到TTI介质纯P波方程.在这个过程中,对描述纯P波相速度的Harlan近似公式进行改进,提高其近似精度.VTI介质P波相速度精确公式如下所示(Fowler,2003):

(1)

其中,vPzvSz分别是P波和SV波沿对称轴方向的相速度,θ是相角,vPx2=(1+2)vPz2vPn2=(1+2)vPz2εδ是各向异性参数(Thomsen,1986).根据声学近似假设(Alkhalifah,2000),设沿着对称轴方向的横波速度vSz=0,并记vPe2(θ)=vPx2sin2θ+vPz2cos2θ,可以得到如下的相速度近似式:

(2)

公式(2)与Alkhalifah(2000)提出的频散公式是密切相关的,两者的不同之处在于,除了P波,Alkhalifah(2000)给出的频散公式还包含了一个相速度为的菱形横波.为了下文表述方便,我们仍旧称公式(2)为Alkhalifah近似式.利用泰勒一阶近似消除公式(2)中的根号项,可以得到如下的相速度近似式:

(3)

公式(3)是Muir和Dellinger(1985)给出的相速度公式的一种特殊情况,这里称公式(3)为Muir-Dellinger近似式.将关系式和cosθ=代入公式(3),其中ω为角频率,kxkz分别为xz方向上的波数,则可以得到如下的频散公式:

(4)

通过公式(4)我们可以推导得到多种不同形式的波动方程,例如:(1)将公式(4)两侧同时乘上vPx2kx2+vPz2kz2,并进行时间和空间上的Fourier反变换,则可以得到一个四阶的偏微分方程(Klié and Toro,2001),或者通过引入辅助变量,将其转化为一个二阶偏微分方程组(Chu et al.,2011; 杜启振等,2015),但是不论是四阶方程还是二阶方程组都需要求解大型线性方程组,求解过程复杂,计算效率低;(2)对公式(4)进行时间方向的Fourier反变换,得到时间-波数域方程,利用谱类方法进行求解,比如伪谱法(Kosloff and Baysal,1982).但是公式(4)中的分数项耦合着介质速度,而伪谱法一般要求波数算子不依赖于介质参数(Chu et al.,2013).除此之外,在VTI介质频散公式向TTI介质拓展的过程中,我们还希望利用下面的关系式来简化计算:

(5)

其中${\hat k_x}$和${\hat k_z}$为旋转坐标系下的波数,因此公式(4)中分数项(又称非椭圆项,Fowler,2003)的分母必须不包含与介质有关的参数.Chu等(2013)利用无穷几何级数的方式将分母处的参数ε分离出来,进而采用伪谱法进行求解,虽然他们指出仅需要前几项就可以得到较为满意的精度,但是对于TTI介质,其计算量仍然会显著增加.更为简单的做法是直接将公式(3)中分母处的各向异性参数ε舍掉(例如Harlan,1995; Etgen and Brandsberg-Dahl,2009; Crawley et al.,2010; Pestana et al.,2011; Zhan et al.,2012),即令非椭圆项中的vPe2(θ)=vPz2.若将公式(3)中分母处的各向异性参数ε舍掉,我们可以得到如下的相速度近似式:

(6)

公式(6)与Harlan(1995)所提出的频散公式是相等价的,这里称公式(6)为Harlan近似式.由于利用了关系式vPe2(θ)≈vPz2因此Harlan近似式仅适用于各向异性参数ε比较小的情况,当参数ε比较大时,Harlan近似式便会产生较大的误差.

事实上,我们可以把舍弃参数ε的做法看作是利用vPz2来代替vPe2(θ).基于这种认识,那么我们便可以利用其他的试探速度vPt2代替vPe2(θ),进而有可能寻找到精度更高的近似公式.理论上试探速度vPt有无限种形式,只要其满足vPz2vPt2vPx2即可.除了vPz2之外,本文利用(vPx2+vPz2)/2,vPxvPz和2/(vPx-2+vPx-2)构建了三个改进型Harlan相速度近似公式:

(7)

(8)

(9)

为了对比不同近似公式的精度,本文以六组页岩样品(具体参数见表 1)计算了不同近似公式在0°~90°内的均方根误差(见表 2),并绘制成图(见图 1),然后以North Sea页岩为例,绘制了这6个近似公式相对应的相速度曲线(见图 2a)和相对误差曲线(见图 2b).从图 1图 2中可以看出:(1)Alkhalifah近似式和Muir-Dellinger近似式的精度大体一致,都具有比较高的近似精度,与精确相速度吻合得比较好;(2)常用于推导VTI和TTI纯P波方程的Harlan近似式的精度最低;(3)本文给出的三个改进型Harlan近似式,即公式(7)、(8)和(9)的精度基本一致,虽然低于Alkhalifah近似式和Muir-Dellinger近似式的精度,但是高于Harlan近似式的精度.综合比较,本文采用公式(9)作为进一步推导TTI介质纯P波波动方程的基础公式.

表 1 页岩样品 Table 1 Shale samples
表 2 相速度近似公式均方根误差(%) Table 2 RMS relative errors(%)of phase-velocity approximations
图 1 由不同近似公式得到的相速度均方根误差 其中灰色点线表示公式(2),灰色虚线表示公式(3),灰色实线表示公式(6),黑色点线表示公式(7),黑色虚线表示公式(8),黑色实线表示公式(9). Fig. 1 RMS relative errors of phase velocity approximations The gray dotted line represents Eq.(2). Gray dash line represents Eq.(3). Gray solid line represents Eq.(6). Black dotted line represents Eq.(7). Black dashed line represents Eq.(8). And black solid line represents Eq.(9).
图 2 P波相速度近似曲线(a)及相对误差曲线(b) Fig. 2 P-wave phase velocity approximations;(b)Relative errors compared to the exact P-wave phase velocity

将关系式代入相速度近似公式(9),可以得到如下的频散公式:

(10)

相较与Harlan频散公式,近似式(10)只在其非椭圆项前添加了一个与各向异性参数ε有关的系数(1+ε)/(1+2ε),所以不增加任何计算量.为了得到TTI介质纯P波方程,需要对频散公式应用坐标旋转.设φ表示对称轴相对于垂直方向上的倾角,借助如下的旋转矩阵:

(11)

那么有旋转波数

(12)

将旋转波数代入公式(10)可以得到TTI介质纯P波近似频散公式:

(13)

若将(13)式中具有伪微分算子的乘号项展开,并采用伪谱法进行波场延拓,那么在每个时间步长内需要进行8次Fourier变换.为了尽可能的减少Fourier变换的次数,本文采用混合法进行数值求解,即所有常规微分算子采用有限差分方法计算,而分数形式的伪微分算子采用伪谱法进行计算.应用倍角公式cos2φ=1-2sin2φ=2cos2φ-1进一步合并,有

(14)

将波场变量p(kx,kz,ω)代入公式(14)的两边,并引入辅助变量x,由Fourier变换的性质,i是复数单位,将公式(14)反变换到时间-空间域有

(15)

其中L1=(kz2-kkx2)/(kx2+kz2),L2=(kxkz)/(kx2+kz2),F[]和F-1[]分别表示Fourier变换和逆变换.对于上面的伪微分方程,可以采用伪谱法和中心有限差分进行联合求解,该混合法在每个时间步长内只需要三次Fourier变换和三次有限差分.中心网格有限差分的精度比较低,数值频散严重,尤其方程(15)还涉及到交叉偏导数,因此我们希望能够利用到高精度的交错网格有限差分方法.引入介质密度,将方程改写为如下的一阶应力-速度方程:

(16)

其中,ψ是辅助变量,ρ为介质密度,类比于声波方程,uw分别表示xz方向上的质点振动速度.当ε=δ时,公式(16)退化为椭圆各向异性,当各向异性参数为零时,(16)式退化为声波方程,其解即为声波压力场.

3 旋转交错网格有限差分格式

仔细观察,公式(16)中的变量p和变量u同时对xz方向进行微分,因此公式(16)是无法应用标准交错网格混合法求解,鉴于此,我们采用旋转交错网格有限差分混合法.旋转交错网格是在交错网格基础上发展起来的.在地震波模拟方面,Saenger等(2000)首先系统地研究了旋转交错网格有限差分算法,现今已经成功应用于弹性和黏弹性、各向同性和各向异性介质中的波场模拟(Saenger and Bohlen,2004).旋转交错网格有限差分将同一物理量的不同分量定义在同一网格点上,其中速度分量和介质密度定义在整网格点上,而应力分量(本文为声波压力场)和弹性常量(本文为各向异性参数)定义在半网格点上,如图 3所示.在计算空间导数时,不再沿xz方向进行计算,而是先沿对角线进行差分,然后再将这些对角线方向的计算结果进行线性叠加来得到沿坐标轴方向的空间导数.在二维空间的情况下,两坐标系的变换关系如下(Saenger et al.,2000):

(17)

图 3 旋转交错网格的二维网格单元示意图 Fig. 3 Schematic diagram of 2-D rotated staggered grid

其中,$\tilde x$和$\tilde z$分别代表沿对角线的旋转坐标系方向,Δx和Δz为网格间距,,两坐标系下的求导算子的关系为

(18)

类比标准交错网格有限差分格式,函数f(x,z,t)在$\tilde x$和$\tilde z$方向上的2N阶偏导数分别表示为(Saenger et al.,2000)

(19)

对于地震波场数值模拟来说,数值计算的稳定性是一个重要问题.直接分析TTI介质一阶纯P波方程的旋转交错网格有限差分混合解法的稳定性是相当繁琐的,实际应用也不方便.Saenger和Bohlen(2004)也曾指出对于一般各向异性介质,旋转交错网格的稳定性估计是比较困难的.由前面的推导可知,TTI介质纯P波方程可以认为是VTI介质纯P波方程经过φ角旋转得到的,因而我们只需要考察φ=0°的情况即可.下面是TTI介质一阶纯P波方程的旋转交错网格混合法的时间二阶精度、空间2N阶精度的稳定性条件表达式(推导过程较为繁琐,见附录):

(20)

Δh为空间采样间隔,Δt是时间采样间隔,Vmax是垂直方向的最大速度,am(m=1,2,…,M)是2M阶旋转交错网格差分系数.

4 数值试验与分析

本小节我们对旋转交错网格有限差分混合法一阶纯P波方程进行了三组数值模拟试验,第一组是为了分析TTI介质纯P波方程的运动学和动力学精度,第二组是为了测试TTI纯P方程在强速度变化、强倾角变化介质中的数值稳定性,第三组进行了BP2007 TTI(部分)模型逆时偏移试算.

4.1 均匀模型

首先,在各向异性参数为ε=0.2,δ=-0.05的均匀介质中进行数值模拟,同时与各向异性介质弹性波应力-速度方程进行对比.本数值测例是为了分析TTI介质纯P波方程的运动学和动力学精度,图 4所示为纯P波波场和弹性波波场快照对比图.第一组依次是纯P波方程u变量波场(图a1),弹性波方程u变量波场(图b1)以及两者的波形对比(图c1);第二组是w波场变量.在数值模拟时,时间采用2阶有限差分格式,空间采用10阶差分格式.模型计算区域大小为400×400的正方形,网格间距为Δxz=10 m,沿垂直对称轴P波速度为3000 m·s-1,纵横波速度比为1.73,对称轴倾角为45°,时间步长dt=1 ms.震源子波采用Ricker子波,子波中心频率fm=25 Hz,震源置于模型中间.从图 4可以看出,除了没有SV波之外,纯P波方程模拟的uw波场与弹性波方程模拟的波场传播特征十分吻合,具有很高的近似精度.

图 4 波场变量uw在0.5 s时刻的波场快照 (a)混合法纯P波;(b)旋转交错网格有限差分弹性波;(c)波形对比.(c)中黑色实线代表弹性波,灰色虚线代表纯P波. Fig. 4 Wavefields u and w computed by two equations at the moment of 0.5 s (a)Hybrid PS/RSGFD pure P wave;(b)RSGFD elastic wave;(c)Waveform comparison. In(c),the black solid line represents elastic wave,and the gray dashed line represents pure P wave.

纯P波方程中的声波压力场与弹性波方程中的应力场没有直接的对应关系,不过我们可以认为-2pσ11+σ33(σ11为弹性波场应力水平分量,σ33为应力垂直分量),图 5所示的是由方程(16)计算p和由弹性波方程计算得到的应力组合场-(σ11+σ33)/2,可以明显看到两者在振幅上并不能匹配起来,但是其走时信息是基本一致的.值得注意的是,不论伪声学近似方程还是解耦的纯P方程,它们都不能处理转换波,也就是说,在非均匀介质中,由近似方程控制的波场振幅信息与弹性波方程是不匹配的,这方面的研究仍旧是今后的重要课题.

图 5 混合法纯P波方程声波压力场p,(b)旋转交错网格有限差分弹性波组合应力场-(σ11+σ33)/2, (c)波形对比.(c)中黑色实线代表弹性波,灰色虚线代表纯P波. Fig. 5 Pressure p snapshot computed by hybrid PS/RSGFD pure P-wave equation;(b)Combined stress fields of elastic waves -(σ11+σ33)/2 computed by RSGFD;(c)Waveform comparison. In(c),the black solid line represents elastic wave, and the gray dashed line represents pure P wave.
4.2 楔形模型

图 6a所示是一个具有强分界面的楔形模型.Duveneck和Bakker(2011)利用该楔形模型测试了不同的TTI声波近似方程的稳定性,Zhan等(2012)亦利用该楔形模型测试了基于快速扩展算法的TTI纯P方程的稳定性.我们对这个模型做了修改,使其中的各向异性参数存在ε<δ的情况,模型从区域1到区域2物性参数发生了很大的变化,对称轴的倾角变化接近90°,具体物性参数见表 3.模型计算区域为400×400的正方形,网格间距为Δxz=10 m,时间步长dt=1 ms.震源子波采用Ricker子波,子波中心频率fm=25 Hz,震源置于模型中间.图 6b—6c所示分别为0.8 s、1.2 s和4 s时的波场快照,地震波在整个延拓过程十分稳定,并未出现数值溢出的不稳定现象.

图 6 非均匀楔形模型(a)和0.8 s(b)、1.2 s(c)及4 s(d)时的波场快照 Fig. 6 Heterogeneous wedge model;(b)—(d)Snapshots of wavefield at time 0.8 s,1.2 s and 4 s
表 3 2D非均匀楔形模型参数表 Table 3 Parameters of the 2D heterogeneous wedge model
4.3 逆时偏移数值试验

BP 2007 TTI模型是典型的各向异性介质模型,我们选取了其右侧的一部分模型体(图 7),进行了重新采样,模型的网格大小设置为10 m×10 m,垂直方向包含900个采样点最深达9 km,水平方向上包含1100个采样点.震源子波采用主频为30 Hz的Ricker子波,时间步长dt=1 ms,置于地下12.5 m处.图 8是单炮地震波场快照.图 9是基于声波逆时偏移(a)、VTI逆时偏移(b)和TTI逆时偏移(c)的成像结果.从图 9的对比中明显可见,TTI逆时偏移的成像效果要优于各向同性逆时偏移的结果,构造更加清晰,成像位置也更准确.表 4是这3种逆时偏移算法的单炮平均耗时,计算机处理器为24核Intel(R)Xeon(R)E5-2695型CPU,各向同性声波逆时偏移基于一阶声波应力-速度方程,VTI逆时偏移则基于杜启振等(2015)所采用的VTI纯P波一阶应力-速度方程.

图 7 BP2007 TTI模型(部分) Fig. 7 BP2007 TTI model(part)
图 8 BP2007 TTI模型1.5 s和3 s时的波场快照 Fig. 8 Wavefield snapshots of BP2007 TTI model at 1.5 s and 3 s
图 9 BP2007 TTI模型逆时偏移对比.(a)各向同性逆时偏移剖面;(b)VTI逆时偏移剖面;(c)TTI逆时偏移剖面 Fig. 9 RTM comparison on BP2007 TTI model with(a)isotropic algorithm,(b)VTI algorithm and(c)TTI algorithm
表 4 2D各向同性、VTI及TTI逆时偏移算法的单炮平均耗时 Table 4 Average computing time of 2D isotropic,VTI and TTI RTM
5 结论

为了得到TTI介质高精度的纯P波延拓算子,本文对常用于推导TTI介质纯P波方程的Harlan相速度公式进行了分析与改进,利用不同的替换速度代替近似式中的非椭圆项中的vPe2(θ),得到了三个改进型Harlan近似公式,在非椭圆项前添加了一个与各向异性参数ε有关的修正系数.相速度曲线对比结果表明,本文给出的改进型相速度近似式比Harlan近似式的精度更高.以其中近似精度最高的改进式为基础,本文推导了TTI介质纯P波方程,并以此为基础构建波场正、反向延拓算子,实现了TTI介质纯P波方程正演模拟与逆时偏移.在数值计算的过程中采用旋转交错网格混合求解此纯P波方程,减少求解过程中Fourier变换的数量,提高计算效率.数值试验证明本文给出的方法能够较准确地描述TTI介质P波波场的运动学特征,同时具有稳定性好、数值频散低、模拟精度高等特点,可以应用到TTI介质逆时偏移.

附录 TTI介质一阶纯P波方程旋转交错网格混合解法稳定性分析

假设波动方程有平面波解ξ(x,z,t)=ξ(0,0,0)ei(kxx+kzz)e-iωt,那么有

(A1)

有限差分时间二阶递推格式可以认为是对方程(A1)中三角函数cos(ωΔt)的Taylor级数二阶截断,即

(A2)

为了数值稳定,必须限制在[-1,1],那么必须有

(A3)

结合前面公式(10),有

(A4)

(A5)

首先考察不等式(A4).当ε≤δ时,式(A4)成立;当ε>δ时,只需要2(εδ)kx2≤(1+2ε)kx2成立即可,那么

(A6)

接下来考察不等式(A5).由于

因此只需要即可.

将波动方程平面波解代入离散格式(19)有

(A7)

将式(A7)代入式(18)有

(A8)

由关系式

(A9)

令Δxzh,那么有

(A10)

因此只要就成立,约束条件(A5)必成立,那么有

(A11)

参考文献
Alkhalifah T. 1998. Acoustic approximations for processing in transversely isotropic media. Geophysics, 63(2): 623-631. DOI:10.1190/1.1444361
Alkhalifah T. 2000. An acoustic wave equation for anisotropic media. Geophysics, 65(4): 1239-1250.
Baysal E, Kosloff D D, Sherwood J W C. 1983. Reverse time migration. Geophysics, 48(11): 1514-1524. DOI:10.1190/1.1441434
Cheng J B, Kang W, Wang T F. 2013. Description of qP-wave propagation in anisotropic media, Part I:Pseudo-pure-mode wave equations. Chinese J. Geophys. (in Chinese), 56(10): 3474-3486. DOI:10.6038/cjg20131022
Chu C, Macy B K, Anno P D. 2011. An accurate and stable wave equation for pure acoustic TTI modeling.//81st Annual Meeting, SEG, Expanded Abstracts, 179-184.
Chu C L, Macy B K, Anno P D. 2013. Pure acoustic wave propagation in transversely isotropic media by the pseudospectral method. Geophysical Prospecting, 61(3): 556-567. DOI:10.1111/gpr.2013.61.issue-3
Crawley S, Brandsberg-Dahl S, McClean J, et al. 2010. TTI reverse time migration using the pseudo-analytic method. The Leading Edge, 29(11): 1378-1384. DOI:10.1190/1.3517310
Du Q Z, Qin T. 2009. Multicomponent prestack reverse-time migration of elastic waves in transverse isotropic medium. Chinese J. Geophys. (in Chinese), 52(3): 801-807.
Du Q Z, Guo C F, Gong X F. 2015. Hybrid PS/FD numerical simulation and stability analysis of pure P-wave propagation in VTI media. Chinese J. Geophys. (in Chinese), 58(4): 1290-1304. DOI:10.6038/cjg20150417
Duveneck E, Milcik P, Bakker P M, et al. 2008. Acoustic VTI wave equations and their application for anisotropic reverse-time migration.//78th Annual International Meeting, SEG, Expanded Abstracts, 2186-2190.
Duveneck E, Bakker P M. 2011. Stable P-wave modeling for reverse-time migration in tilted TI media. Geophysics, 76(2): S65-S75. DOI:10.1190/1.3533964
Etgen J T, Brandsberg-Dahl S. 2009. The pseudo-analytical method:Application of pseudo-Laplacians to acoustic and acoustic anisotropic wave propagation.//79th Annual International Meeting. SEG, Expanded Abstracts, 2552-2556.
Fletcher R P, Du X, Fowler P J. 2009. Reverse time migration in tilted transversely isotropic (TTI) media. Geophysics, 74(6): WCA179-WCA187. DOI:10.1190/1.3269902
Fowler P J. 2003. Practical VTI approximations:a systematic anatomy. Journal of Applied Geophysics, 54(3-4): 347-367. DOI:10.1016/j.jappgeo.2002.12.002
Fowler P J, Du X, Fletcher R P. 2010. Coupled equations for reverse time migration in transversely isotropic media. Geophysics, 75(1): S11-S22. DOI:10.1190/1.3294572
Grechka V, Zhang L B, Rector J W Ⅲ. 2004. Shear waves in acoustic anisotropic media. Geophysics, 69(2): 576-582. DOI:10.1190/1.1707077
Harlan W S. 1995. Flexible seismic traveltime tomography applied to diving waves. Stanford Exploration Project Report, 89: 145-164.
Huang T, Zhang Y, Zhang H Z, et al. 2009. Subsalt imaging using TTI reverse time migration. The Leading Edge, 28(4): 448-452. DOI:10.1190/1.3112763
Jones L E A, Wang H F. 1981. Ultrasonic velocities in Cretaceous shales from the Williston basin. Geophysics, 46(3): 288-297. DOI:10.1190/1.1441199
Klié H, Toro W. 2001. A new acoustic wave equation for modeling in anisotropic media.//71st Annual Meeting, SEG, Expanded Abstracts, 1171-1174.
Kosloff D D, Baysal E. 1982. Forward modeling by a Fourier method. Geophysics, 47(10): 1402-1412. DOI:10.1190/1.1441288
Li B, Li M, Liu H W, et al. 2012. Stability of reverse time migration in TTI media. Chinese J. Geophys. (in Chinese), 55(4): 1366-1375. DOI:10.6038/j.issn.0001-5733.2012.04.032
Liu F Q, Morton S A, Jiang S S, et al. 2009. Decoupled wave equations for P and SV waves in an acoustic VTI media.//79th Annual International Meeting. SEG, Expanded Abstracts, 2844-2848.
McMechan G A. 1983. Migration by extrapolation of time-dependent boundary values. Geophysical Prospecting, 31(3): 413-420. DOI:10.1111/gpr.1983.31.issue-3
Muir F, Dellinger J. 1985. A practical anisotropic system. Stanford Exploration Project Report, 44: 55-58.
Pestana R C, Stoffa P L. 2010. Time evolution of the wave equation using rapid expansion method. Geophysics, 75(4): T121-T131. DOI:10.1190/1.3449091
Pestana R C, Ursin B, Stoffa P L. 2011. Separate P- and SV-wave equations for VTI media.//SEG Technical Program Expanded Abstracts, 163-167.
Saenger E H, Gold N, Shapiro S A. 2000. Modeling the propagation of elastic waves using a modified finite-difference grid. Wave Motion, 31(1): 77-92. DOI:10.1016/S0165-2125(99)00023-2
Saenger E H, Bohlen T. 2004. Finite-difference modeling of viscoelastic and anisotropic wave propagation using the rotated staggered grid. Geophysics, 69(2): 583-591. DOI:10.1190/1.1707078
Thomsen L. 1986. Weak elastic anisotropy. Geophysics, 51(10): 1954-1966. DOI:10.1190/1.1442051
Tsvankin I, Thomsen L. 1994. Nonhyperbolic reflection moveout in anisotropic media. Geophysics, 59(8): 1290-1304. DOI:10.1190/1.1443686
Vernik L, Liu X Z. 1997. Velocity anisotropy in shales:A petrophysical study. Geophysics, 62(2): 521-532. DOI:10.1190/1.1444162
Wang Z J. 2002. Seismic anisotropy in sedimentary rocks, part 2:Laboratory data. Geophysics, 67(5): 1423-1440. DOI:10.1190/1.1512743
Whitmore N D. 1983. Iterative depth migration by backward time propagation.//53rd Annual International Meeting, SEG, Expanded Abstracts, 382-385.
Zhan G, Pestana R C, Stoffa P L. 2012. Decoupled equations for reverse time migration in tilted transversely isotropic media. Geophysics, 77(2): T37-T45. DOI:10.1190/geo2011-0175.1
Zhan G, Pestana R C, Stoffa P L. 2013. An efficient hybrid pseudospectral/finite-difference scheme for solving the TTI pure P-wave equation. J. Geophys. Eng., 10(2): 025004. DOI:10.1088/1742-2132/10/2/025004
Zhang H Z, Zhang G Q, Zhang Y. 2009. Removing S-wave noise in TTI reverse time migration.//79th Annual International Meeting, SEG, Expanded Abstracts, 2849-2853.
Zhang Y, Zhang H Z, Zhang G Q. 2011. A stable TTI reverse time migration and its implementation. Geophysics, 76(3): WA3-WA11. DOI:10.1190/1.3554411
Zhang Y, Wu G Z. 2013. Methods of removing pseudo SV-wave artifacts in TTI media qP-wave reverse-time migration. Chinese J. Geophys. (in Chinese), 56(6): 2065-2076. DOI:10.6038/cjg20130627
Zhou H B, Zhang G Q, Bloor R. 2006a. An anisotropic acoustic wave equation for modeling and migration in 2D TTI media.//76th Annual International Meeting, SEG, Expanded Abstracts, 194-198.
Zhou H B, Zhang G Q, Bloor R. 2006b. An anisotropic acoustic wave equation for VTI media.//68th EAGE Conference & Exhibition. SPE, EAGE.
程玖兵, 康玮, 王腾飞. 2013. 各向异性介质qP波传播描述I:伪纯模式波动方程. 地球物理学报, 56(10): 3474–3486. DOI:10.6038/cjg20131022
杜启振, 秦童. 2009. 横向各向同性介质弹性波多分量叠前逆时偏移. 地球物理学报, 52(3): 801–807.
杜启振, 郭成锋, 公绪飞. 2015. VTI介质纯P波混合法正演模拟及稳定性分析. 地球物理学报, 58(4): 1290–1304. DOI:10.6038/cjg20150417
李博, 李敏, 刘红伟, 等. 2012. TTI 介质有限差分逆时偏移的稳定性探讨. 地球物理学报, 55(4): 1366–1375. DOI:10.6038/j.issn.0001-5733.2012.04.032
张岩, 吴国忱. 2013. TTI介质qP波逆时偏移中伪横波噪声压制方法. 地球物理学报, 56(6): 2065–2076. DOI:10.6038/cjg20130627