复杂TTI介质稳定的纯qP波波场模拟方法

胡书华　王宇超　刘文卿　周　阳

(①中国石油勘探开发研究院西北分院，甘肃兰州 730020；②同济大学海洋与地球科学学院波现象与反演成像研究组，上海 210092)

1　引言

随着勘探目标的复杂化与多样化，地震勘探对地下介质精确成像的需求不断提高。Baysal等[1]于1983年提出的逆时偏移(Reverse time migration，RTM)技术，因其成像精度高、适应复杂介质等优点，已得到广泛的研究与应用[2-5]； 而作为纵波勘探中逆时偏移技术的基础，高效而精确的准纵波(quasi-P wave，qP波)模拟算法也处于持续发展中。同时，由于地下介质普遍具有各向异性特性，因此当前精确成像算法对地下各向异性高精度近似的要求越来越高[6]。

自1998年Alkhalifah[7]提出将各向异性物性对称轴方向横波速度设置为零的声学近似思想以来，基于该思想的各向异性介质qP波模拟算法得到快速发展： Alkhalifah[8]导出声学近似下TI介质四阶控制方程， Zhou等[9]将该方程分解为两个二阶微分方程； Du等[10]和Fowler等[11]也导出了不同形式的二阶qP波控制方程。

近年，Fletcher等[12]、Duveneck等[13]、Zhang等[14]发现这些方程在实际应用中会遇到“声学近似导致的横波干扰”和“物性对称轴变化剧烈处波场传播不稳定”两个问题，分析其原因并提出了不同的解决方案： Flecher等[12]采用添加有限横波速度的方法，导出有限横波方程，消除伪横波并保证稳定性； Duveneck等[13]则从一般介质本构方程出发，导出了更加符合物理意义的稳定的qP波传播控制方程； Zhang等[14]在该方程中添加自共轭二阶微分算子来保证其稳定性。

近期关于纯qP波模拟的研究主要集中在有效各向同性模型近似算法和拟微分算子分解方法。如Alkhalifah等[15]的各向同性近似算法在很大程度上改善了RTM效率； Xu等[16，17]提出拟微分算子分解方法并对其进行了改进，将伪微分算子分解为一个标量算子和一个微分算子，该方法得到的波场不受伪横波干扰，且波场模拟精度较高。

本文基于Xu等[17]的拟微分算子分解思想，分析了分解后标量算子和椭圆微分算子的频散特性，推导出对应的TTI介质一阶qP波控制方程，并用交错网格高阶有限差分算法对方程实施求解，得到复杂TTI介质稳定的纯qP波模拟的结果。

2　新的纯qP波控制方程

2.1　算子分解思路

从一般介质本构关系出发，求解Christoffel方程得到qP-qSV耦合相速度公式，进而可导出各向异性介质频散关系。据Alkhalifah[7]声学近似思想，令VSz=0，有VTI介质qP-qSV耦合频散关系

(1)

式中：ω为角频率；VPz为P波垂向速度；ε和δ为Thomsen各向异性参数；kx、ky和kz分别为x、y和z方向的波数。

考虑Alkhalifah[8]给出的VTI介质标量拟微分方程

(2)

该式即为VTI介质qP波传播控制方程，注意方程中的所有参数都是空间变化的。式(2)的频散关系与式(1)相同，即与弹性P波的频散关系一致，表明其能够准确地模拟qP波传播。然而，式(2)是一个伪微分方程，在数值上不易求解。因此，Xu等[17]提出将式中伪微分算子分解为一个椭圆微分算子和一个标量算子的思路。式(2)对应的频散关系可容易地表达成

(3)

定义单位波数矢量

(4)

将其代入式(3)，同时对式(3)中根式项进行变换，得到

(5)

记式(5)等号右侧含有单位波数矢量分量的乘积项为标量算子Se，剩余部分变换到时空域，则为椭圆微分算子

(6)

(7)

式(6)即为椭圆微分算子； 式(7)为标量算子。

式(3)～式(7)即为算子分解下的椭圆微分算子分解方法。算子分解的主要目的在于改造原伪微分方程，将其中不易求解的伪微分算子分解为容易求解的椭圆微分算子和标量算子，在保持波场模拟精度的同时实现方程高效求解。

2.2　算子频散特性分析

考虑对新算子的频散特性进行分析，同时与精确qP波频散关系进行对比，以研究新算子对于波场的作用。

式(6)对应的频散关系可写成

(8)

考虑二维情况，给定参数VPz=3000m/s、ε=0.24、δ=0.1，绘制精确频散关系曲线(式(1))、椭圆微分算子频散关系曲线(式(8))及本文推导的式(6)和式(7)对应的频散关系曲线。由图1可见，椭圆微分算子接近于精确qP波频散关系，然而两者存在明显的数值差； 分析式(8)可知两者相差的量来自于标量算子Se的作用。同时，采用本文新方程所得频散关系曲线与精确频散关系曲线相一致。因此可知，该伪微分算子分解方法实际上相当于用标量算子Se对椭圆微分算子进行校正。

进一步分析校正项Se的作用，记标量算子Se=1+ΔSe，则有

(9)

二维情况下，将频散关系式转换到极坐标系(R，θ)下，绘制的精确频散关系曲线、椭圆微分算子关系曲线、标量算子ΔSe项频散曲线以及本文新算子的频散关系曲线如图2所示。可见ΔSe项的数值在-0.05～0之间，是对椭圆微分算子的一种微小补偿，目的在于将其校正到精确频散关系。同时从该图可知，式(6)与式(7)组成的新算子可保持精确的qP波频散关系， 因此具有较高的qP波模拟精度。

图1　　　　精确的qP波频散关系、椭圆微分算子及本文频散关系对比

图2　　　　极坐标下精确频散关系、标量算子、椭圆微分算子及本文频散关系对比

2.3　TTI介质纯qP波控制方程

为了数值实施方便，引入辅助波场p、q和r，同时将式(6)写为一阶形式

(10)

Se=

(11)

式(11)与式(10)共同组成了VTI介质纯qP波控制方程。

对于TTI介质，仿照Zhang等[14]和Bube等[18]的确保稳定性做法，在旋转坐标系下引入自共轭微分算子

(12)

式中： “·”为u、p、q、r中任一种量；θ表示TTI介质物性对称轴倾角；φ为对称轴方位角； 上标“T”表示转置运算。用新的算子依次替换式(10)和式(11)中的微分算子，可得旋转坐标系下的TTI介质纯qP波传播控制方程

(13)

(14)

式(12)～式(14)共同构成TTI介质纯qP波传播控制方程。该控制方程来自于伪微分算子分解得到的新算子，与精确的qP波频散关系相同，且方程推导过程中未做任何近似，将此用于TTI介质qP波模拟具有较高的精度。

3　数值算法

本文采用交错网格高阶有限差分数值解法求解式(12)～式(14)，在t时刻计算质点位移u，在t-Δt/2时刻计算辅助波场p、q和r。

考察式(12)～式(14)发现，t+Δt时刻波场u的计算不仅需要t时刻的波场u，还需t+Δt/2时刻辅助波场p、q和r分别关于x、y和z的微分； 同样地，t+Δt/2时刻辅助波场p、q和r的计算不仅需t-Δt/2时刻的波场p、q和r，还需t时刻波场u关于x、y和z的微分。因此，为了避免插值，提高计算精度，本文采用Lebedev交错网格[19，20](图3)。其中：在方点位置处计算波场u和标量算子Se； 在圆点位置处计算辅助波场p、q和r。

(15)

(16)

式中：f(t)为震源函数；F为f的离散变量。

图3　式(12)～式(14)的Lebedev交错差分网格格式

这里需要注意的是，式(13)和式(14)中虽然都有Gxu、Gyu和Gzu项，但前者在p、q和r的位置计算，而后者在u的位置计算。两者计算的位置不同，虽然可以通过插值得到不同位置的数值，但为了保证计算的精度，对两处位置分别计算波场u的梯度。两者差分公式略有区别，可根据网格位置类似写出。计算复杂度虽然增加了波场u梯度的1次计算，但保证了波场模拟的精度。

至此，给出了用Lebedev交错网格高阶有限差分求解式(12)～式(14)的数值算法。为了保证波场模拟的精度和稳定性，引入了Lebedev交错网格，同时采用波场梯度2次计算的数值算法策略。

4　模型数据测试

4.1　均匀TTI介质模型

给定均匀TTI介质参数：VPz=3000m/s，ε=0.1，δ=0.024，网格数为400×400，网格尺寸10m×10m，震源子波主频为16Hz，震源位于网格中心位置，时间步长1ms，时间方向采用2阶差分，空间方向采用10阶差分。其0.6s的波场快照如图4所示。可见Fletcher有限横波法和本文方法得到的波场中qP波成分与弹性波场中P波成分一致，表明两者均能较好地保持P波传播特征。但是有限横波法得到的qP波场中存在明显的残余横波干扰，而本文方法得到的纯qP波场已经完全消除了残余横波干扰，表明本文方法不受伪横波干扰的影响。

4.2　双层TTI介质模型

设计双层TTI介质模型如图5所示，模型在2500m深度位置存在反射界面。对其分别进行常规网格有限差分弹性波模拟、有限横波法模拟以及本文提出的纯qP波方法模拟，所得结果如图6所示。模拟参数与上述均匀TTI介质一致。

图4　均匀TTI介质模型600ms波场快照对比

图5　水平双层TTI介质模型

图6　水平双层TTI介质模型550ms波场快照对比

分析图6可知，从全弹性波波场快照可以分辨出直达纵波、直达横波以及各种转换波，如反射纵波、反射横波、透射纵波、透射横波；而有限横波法得到的波场中不仅存在残余直达横波干扰，也存在残余的透射和反射横波干扰；本文方法得到的纯qP波波场中只有纵波成分，不存在残余横波干扰，且纵波波前特征与全弹性波场中的纵波特征一致，证明本文方法的有效性。

4.3　复杂TTI介质模型

本文从BP-TTI模型中截取一块各向异性参数复杂的区域，如图7所示。模型网格数为1801×1801，网格尺寸为6.25m×6.25m。可看到模型中存在大量陡倾断层，各向异性对称轴倾角参数θ在某些位置有突变现象，这对于波场的稳定模拟带来挑战。

将震源放置于模型中心位置处，采用常规网格有限差分弹性波模拟方法、有限横波方法、本文方法分别进行波场模拟，得到1000ms和1500ms波场快照如图8和图9所示。

图9　BP TTI模型1500ms波场快照对比

由波场模拟结果可知，有限横波法和本文方法所得到的波场，其P波波前特征与全弹性波波场一致，表明两者对于波场中的纵波成分能较好地进行模拟。但是，有限横波法波场中存在残余伪横波干扰，且残余伪横波在传播过程中又产生残余的反射、透射横波干扰；而本文方法能得到干净的纯P波波场，不存在残余横波干扰。表明本文提出的方法能够适应各向异性物性对称轴参数剧烈变化情况，得到稳定的无横波干扰的纯qP波波场。

5　结论

本文基于Xu的伪微分算子分解思路，结合旋转坐标系下的自共轭微分算子导出稳定的一阶纯qP波控制方程，同时在Lebedev交错网格框架下给出其高阶有限差分算法，在复杂各向异性模型中得到了稳定传播的qP波场。由算子频散分析以及均匀TTI模型、双层TTI模型以及复杂TTI介质模型纯qP波波场模拟试验可知，相对于其他qP波模拟算法，本文提出的纯qP波模拟算法不受伪横波的干扰，且对各向异性物性对称轴参数变化有更好的适应性，在物性对称轴倾角参数变化剧烈时也能得到稳定的纯qP波波场。

[1]Baysal E,Kosloff D D,Sherwood J W C.Reverse time migration.Geophysics,1983,48(11)：1514-1524.

[2]刘文卿,王宇超,雍学善等.基于GPU/CPU叠前逆时偏移研究及应用.石油地球物理勘探,2012,47(5)： 712-716.

Liu Wenqing,Wang Yuchao,Yong Xueshan et al.Prestack reverse time migration on GPU/CPU co-parallel computation.OGP,2012,47(5)：712-716.

[3]刘文卿,王西文,刘洪等．盐下构造速度建模与逆时偏移成像研究及应用．地球物理学报,2013,56(2)：612-625．

Liu Wenqing,Wang Xiwen,Liu Hong et al.Application of velocity modeling and reverse time migration to subsalt structure.Chinese Journal of Geophysics,2013,56(2)：612-625.

[4]刘定进,蒋波,李博等.起伏地表逆时偏移在复杂山前带地震成像中的应用.石油地球物理勘探,2016,51(2)：319-324.

Liu Dingjin,Jiang Bo,Li Bo et al.Rugged-topography reverse time migration in complex piedmont zone.OGP,2016,51(2)：219-324.

[5]姜国博,曾庆芹,李虹等.逆时偏移技术在土库曼斯坦地区的应用.石油地球物理勘探,2017,52(增刊2): 81-85,103.

Jiang Guobo,Zeng Qingqin,Li Hong et al.Application of reverse time migration in Turkmenistan area.OGP,2017,52(S2): 81-85,103.

[6]韩令贺,何兵寿.VTI介质一阶准P波方程正演模拟及边界条件.石油地球物理勘探,2010,45(6)：819-825.

Han Linghe,He Bingshou.First-order Quasi-P wave equation forward modeling and boundary conditions in VTI medium.OGP,2010,45(6)：819-825.

[7]Alkhalifah T.Acoustic approximation for processing in transversely isotropic media.Geophysics,1998,63(2)： 623-631.

[8]Alkhalifah T.An acoustic wave equation for aniso-tropic media.Geophysics,2000,65(4)：1239-1250.

[9]Zhou H,Zhang G,Bloor R.An anisotropic acoustic wave equation for modeling in 2D TTI Media.SEG Technical Program Expanded Abstracts,2006,25：194-198.

[10]Du X,Fletcher R P,Fowler P J.A new pseudo-acoustic wave equation for VTI media.70th EAGE Conference & Exhibition,Extended Abstracts,2008：H033.

[11]Fowler P J,Du X,Fletcher R P.Coupled equations for reverse time migration in transversely isotropic media.Geophysics,2010,75(1)：S11-S22.

[12]Fletcher R P,Du X,Fowler P J.Reverse time migration in tilted transversely isotropic (TTI) media.Geo-physics,2009,74(6)：WCA179-WCA187.

[13]Duveneck E,Bakker P M.Stable P-wave modeling for reverse time migration in tilted TI media.Geophy-sics,2011,76(2)：S65-S75.

[14]Zhang Y,Zhang H Z,Zhang G Q.A stable TTI re-verse time migration and its implementation.Geophysics,2011,76(3)：WA3-WA11.

[15]Alkhalifah T,Ma X,Waheed U B et al.Efficient anisotropic wavefield extrapolation using effective isotropic models.75th EAGE Conference & Exhibition,2013：Tu-01-16.

[16]Xu S,Zhou H.Accurate simulations of pure quasi-P-waves in complex anisotropic media.Geophysics,2014,79(6)：1-8.

[17]Xu S,Tang B,Mu J et al.Elliptic decomposition of quasi-P wave equation.77th EAGE Conference & Exhibition,2015：We P6 12.

[18]Bube K P,Nemeth T,Stefani J P et al.On the instability in second-order systems for acoustic VTI and TTI media.Geophysics,2012,77(5)：T171-T186.

[19]Hu S H,Wang X W,Sun J Q et al.Stable simulating algorithm of pure quasi-P wavefield in complex anisotropic media.78th EAGE Conference & Exhibition,2016：We SRS3 04.

[20]李娜,李振春,黄建平等.Lebedev网格与标准交错网格耦合机制下的复杂各向异性正演模拟.石油地球物理勘探,2014,49(1)：126-131.

Li Na,Li Zhenchun,Huang Jianping et al.Numerical simulation with coupling Lebedev and standard staggered grid schemes for complex anisotropic media.OGP,2014,49(1)：126-131.

[21]董良国,马在田,曹景忠等.一阶弹性波方程交错网格高阶差分解法.地球物理学报,2000,43(3)：411-419.

Dong Liangguo,Ma Zaitian,Cao Jingzhong et al.A staggered-grid high-order difference method of one-order elastic wave equation.Chinese Journal of Geophysics,2000,43(3)：411-419.