地球物理学报  2016, Vol. 59 Issue (4): 1477-1490   PDF    
基于LS-RSGFD方法优化的横向各向同性(TI)介质一阶qP波高精度数值模拟
李振春, 杨富森, 王小丹    
中国石油大学(华东)地球科学与技术学院, 青岛 266580
摘要: 为克服各向异性介质弹性波数值模拟中存在着计算量大和波场分离困难等局限,研究了声学近似的VTI介质和TTI介质一阶qP波数值模拟方法.首先对VTI介质弹性波方程进行声学近似,推导了VTI介质一阶qP波方程;然后基于精确的TTI介质频散关系,引入一个包含各向异性控制参数σ的新辅助波场,推导了稳定的TTI介质二阶耦合qP波波动方程,并通过引入波场的伪速度分量,推导了等价的一阶应力-速度形式.结合旋转交错网格有限差分(RSGFD)和基于最小二乘优化的有限差分(LS-FD)两种各具优势的方法,研究了最小二乘旋转交错网格有限差分(LS-RSGFD)方法,并用其数值求解VTI和TTI介质一阶qP波方程,然后通过构造其LS-RSGFD格式,实现了高精度的各向异性介质qP波波场数值模拟.数值模拟结果表明:TI介质一阶qP波方程能够准确地模拟各向异性介质中qP波的运动学特征,引入控制参数σ能够有效地减弱不稳定性问题,保证非均匀TTI介质中qP波场的稳定传播;利用优化的LS-RSGFD方法可以得到高精度的合成地震记录,同时还可以相对地提高计算效率.
关键词: TI介质     LS-RSGFD方法     一阶qP波     数值模拟     PML边界条件     控制参数     不稳定性    
High-precision numerical simulation of first-order qP-waves in the transversely isotropic (TI) medium optimized by the LS-RSGFD method
LI Zhen-Chun, YANG Fu-Sen, WANG Xiao-Dan    
Department of Geoscience, China University of Petroleum (East China), Qingdao 266580, China
Abstract: In order to overcome the limitations of anisotropic elastic wave modeling, such as large computational costs, and the difficulty of wave-field separation, we have studied the numerical simulation methods of first-order qP-waves in the acoustic approximated VTI and TTI media. Starting from the acoustic approximated VTI elastic wave equation, we derived the VTI first-order qP-wave equation. Then based on the accurate TTI dispersion relation, the stable TTI second-order coupled qP-wave equation was derived by introducing a new auxiliary function including an anisotropic control parameter σ, and the equivalent first-order stress-velocity form was deduced. Combining the RSGFD with the LS-FD two advantageous method, we studied the LS-RSGFD method to solve VTI and TTI first-order qP-wave equation numerically, and then implemented high-precision numerical simulation through constructing the LS-RSGFD form of first-order qP-wave equation. Numerical results demonstrate that TI first-order qP-wave equation can describe the kinematics features of qP-waves in anisotropic media; introducing a control parameter can mitigate the instability problem and thus stabilize qP-wave propagation in heterogeneous TTI media; and applying the optimal LS-RSGFD method can acquire high-precision synthetic seismic recordings. Meanwhile it can enhance computational efficiency relatively.
Key words: TI media     LS-RSGFD method     First-order qP-wave     Numerical simulation     PML boundary condition     Control parameter     Instability    
1 引言

近年来,随着勘探目标的复杂化和大偏移距、宽方位角地震采集技术的发展,对资料处理的精度要求也越来越高,使得地球介质的各向异性成为勘探地球物理界的研究热点.横向各向同性(TI)介质是沉积岩层中最常见的一种各向异性介质,包括VTI、HTI和TTI介质.常规的基于各向同性假设的正演模拟方法已不足以准确地描述实际地球介质中地震波的传播特征,在以往的研究中,各向异性介质弹性波数值模拟方法取得一定的进展(裴正林,2004),但其存在弹性参数物理意义不明确、计算量大,以及波场分离困难等缺陷(韩令贺和何兵寿,2010).由于实际地震勘探中仍以纵波资料为主,因此,研究稳定的TI介质qP波波动方程及其高精度的数值模拟方法显得尤为重要.

Alkhalifah(2000)从精确的VTI介质qP-qSV耦合频散关系出发,利用声学假设近似,推导了VTI介质四阶qP波方程.Zhou等(2006a)Du等(2008)通过引入不同的辅助波场函数对其进行降阶处理,分别推导了两种二阶耦合的VTI介质qP波方程.Hestholm(2009)进一步推导了一阶偏微分形式的VTI介质qP波方程,并利用交错网格有限差分方法数值求解.Duveneck等(2008)从声学近似的VTI介质虎克定律出发,也推导了二阶和一阶应力-速度形式的VTI介质qP波方程.Zhang等(2005)Alkhalifah(2000)的VTI介质qP波方程推广到TTI介质;Zhou等(2006b)Fletcher等(2008)基于声学近似的TTI介质频散关系,分别推导了两种二阶耦合的TTI介质qP波方程.Fowler等(2010)研究了TI介质二阶耦合qP波方程的一般形式,并证明了它们的等价性.

采用声学近似的各向异性介质qP波方程会产生一种固有的qSV人为干扰波(称为伪横波),它的存在会干扰qP波正演模拟结果,进而还会使得偏移剖面中出现噪声;一些学者提出了不同的伪横波压制方法(Klie and Toro,2001Grechka et al.,2004Duveneck et al.,2008Zhang et al.,2009张岩和吴国忱,2013).同时,对于不均 匀的一般TTI介质,声学近似方法会造成地震波场在对称轴方向急剧变化的区域出现不稳定性现象.Fletcher等(2009)研究发现声学近似后的残余伪横波引起了数值不稳定性,因此引入一个非零的垂向剪切波速度来解决该问题;一些学者也给出了不稳定性问题的解决策略(Zhang H J and Zhang Y,2008Yoon et al.,2010),但这些方法 存在着修改模型参数的局限性;Tessmer(2010)认为声学近似不是真正的弹性各向异性,因而会造成不稳定性问题;因此,Duveneck和Bakker(2011)Zhang等(2011)从弹性各向异性理论出发,通过引入旋转坐标系下的自共轭旋转偏微分算子,分别推导了两种稳定形式的TTI介质二阶耦合qP波方程,但在数值实现过程中需用特殊方法处理.

Du等(2007)Pestana(2011)Zhan等(2012)在弱各向异性假设条件下,推导了时间-波数域的VTI和TTI介质近似解耦qP波方程,在一定程度上避免了伪横波干扰,但其仍需对方程中的波数算子进行带通滤波来保证稳定性,而且需要用伪谱法和REM方法数值求解,三维情况下不利于并行化实现.Liu等(2009)Chu等(2011)对声学近似的各向异性频散关系进行因式分解,推导了VTI和TTI介质的近似纯qP波方程,由于方程中含有四阶混合空间偏导数,计算量大且数值实现复杂.程玖兵等(2013)利用相似变换,推导了动力学上qP波能量占优的VTI和TTI介质伪纯模式qP波方程.Xu和Zhou(2014)将VTI介质声近似qP波方程中的伪偏微分算子进行分解,推导了适用于复杂各向异性介质的纯qP波方程,在此基础上实现了纯qP波正演模拟,并提高了计算效率.

以上方法主要都是基于二阶或四阶偏微分算子的qP波波动方程,采用常规矩形网格对其进行差分离散时,会引起较为严重的频散现象;而利用交错网格有限差分(SGFD)方法可以提高计算精度(李振春,2011),但其需要对弹性模量进行插值处理,对于参数变化较大的非均匀各向异性介质,会导致数值结果不准确.Saenger和Bohlen(2004)利用旋转交错网格有限差分(RSGFD)方法,有效地克服了插值处理带来的误差问题;但对于同样大小的网格单元,RSGFD法比常规SGFD需要更大的空间步长,这往往会增加空间频散,Liu(2013)研究了基于频散关系和最小二乘优化系数的有限差分(LS-FD)方法,可以有效地减弱数值频散,提高模拟精度.

为了提高各向异性介质波场数值模拟的精度和解决不稳定性问题,本文基于声学近似的TI介质弹性波方程与频散关系两种思路,首先推导了VTI介质一阶qP波方程和稳定的TTI介质一阶qP波方程,并给出了相应的PML边界控制方程.然后,结合RSGFD和LS-FD两种各具优势的方法,利用优化的LS-RSGFD方法数值求解TI介质一阶qP波方程,实现了高精度的波场数值模拟,同时引入控制参数σ来有效地减弱波场传播的不稳定性问题,以保证得到稳定的数值解.数值结果对比分析验证了本文方法的有效性和稳定性.

2 横向各向同性(TI)介质一阶qP波方程 2.1 VTI和TTI介质一阶应力-速度qP波方程

基于Alkhalifah(2000)提出的声学假设近似方法,二维VTI介质一阶应力-速度形式的弹性波方程和相应的声近似qP波方程可以分别表示为

式(1a)中,τij是应力分量,u、wxz方向上的速度分量,Cij是弹性常量,它可以表示为Thomsen参数形式(Tsvankin,2001).式(1b)中,p=τxxq=τzz分别是水平和垂直应力分量,vp0是垂向qP波速度;vpx是对称面内qP波速度,vpx2=vp02(1+2ε);vpn是qP波NMO速度,vpn2=vp02(1+2Δ);ε、δ是Thomsen参数(Thomsen,1986).

Alkhalifah(2000)提出的VTI介质四阶声近似qP波方程基础上,Hestholm(2009)引入几个中间变量进行降阶处理,推导了VTI介质一阶qP波方程;韩令贺和何兵寿(2010)印兴耀等(2014)分别利用交错网格和旋转交错网格有限差分方法求解Hestholm(2009)方程,实现了VTI介质一阶qP波数值模拟.杨富森等(2015)研究了几种等价形式的VTI介质一阶qP波方程及其正演模拟方法,并证明了方程(1b)是最有效的形式.因此,对于VTI介质情形,建议使用方程(1b)进行波场数值模拟.

由VTI介质推广到TTI介质只是计算更加复杂,而不会增加物理复杂性.二维情况下,引入对称轴倾角θ,可将方程(1b)推广得到TTI介质一阶应力-速度形式的qP波方程:

对于VTI介质,当各向异性参数满足ε≥δ时,利用方程(1b)进行数值模拟则能保证波场的稳定传播,研究表明,该参数限制条件对大多数实际沉积岩层是适用的.然而,对于非均匀TTI介质,利用方程(2)进行波场外推时会引起严重的数值不稳定性.这种不稳定性不是由速度场和各向异性参数场的差异造成,而是非均匀TTI介质中对称轴方向的剧烈变化造成伪横波相互作用引起的.

2.2 稳定的TTI介质一阶qP波方程

研究表明,各向异性参数σ=(vp0/vs0)2(ε-δ)在很大程度上控制着qSV波的强度和波前形态,同时qSV波的振幅随着σ值减小而减小;更重要的是,σ值的大小不影响qP波的运动学和动力学特征(Tsvankin,2001).因此,本文引入各向异性控制参数σ来保证qP波场的稳定传播,同时有效地衰减伪横波能量.首先从精确的二维TTI介质qP-qSV波耦合频散关系出发(Fowler et al.,2010):

式中,h1=kx2sin2θ+kz2cos2θ+kxkzsin2θ; h2=kx2cos2θ+kz2sin2θ-kxkzsin2θ.w是角频率,kx,kz是二维空间波数;θ是对称轴倾角;vp0vs0分别是沿对称面垂向的qP波和qSV波速度;vpxvpnε,δ的物理意义同方程(1b).然后,引入包含参数σ的新辅助波场q(w,kx,kz):

则可以得到一种新的二维时空域TTI介质二阶耦合qP波波动方程:

式中,;.p 是伪应力波场,q是为简化计算而引入的辅助波场.若令σ=∞(即vs0=0)时,可得到声学近似的特殊情形.

对于方程(5)中的二阶偏微分项,采用常规矩形网格进行差分离散时,会引起较为严重的频散现象,同时其边界反射问题也难以有效解决,因此可以研究等价的一阶应力-速度形式波动方程.为更好地结合旋转交错网格有限差分(RSGFD)技术,本文引入伪应力波场p和q的伪速度分量,并沿着x、z方向进行分解,令;. 则可以将方程(5)转化为等价的TTI介质一阶应力-速度qP波方程:

其中,一阶偏微分算子T1T2G1G2分别为

式中,upwpuqwq分别表示波场pqx、z方向上的伪速度分量;当对称轴倾角θ=0时,方程(6)退化为VTI介质情形.利用弱各向异性近似理论和介质分解理论,可将qSV波的反射系数表示为各向同性项RSViso(θ)和各向异性项RSVaniso(θ)之和(梁锴等,2011):

式中,σ1σ2分别表示反射界面上下层的各向异性控制参数.如果在整个计算空间选取一个足够小的恒定值σ,保证各处的qSV波各向异性项反射系数趋于0,则可以最大程度地衰减qSV波干扰能量,同时保证qP波场在非均匀TTI介质中的稳定传播.

3 TI介质一阶qP波数值模拟 3.1 最小二乘旋转交错网格高阶有限差分(LS-RSGFD)格式的构造

图 1a1b分别表示常规交错网格(SG)和旋转交错网格(RSG)的网格剖分示意图.在图 1a中,伪应力波场pq分布在不同的网格点上,在计算时需要对各向异性参数和倾角进行插值,由此会产生数值误差,甚至出现不稳定.而在图 1b中,伪应力波场分布在同一网格点上,而各向异性参数和倾角被定义在与应力相同的位置,这样可以避免对其进行插值所带来的误差.物理分量和参数的具体分布位置如表 1表 2所示.

图 1 二维网格剖分示意图
(a)常规交错网格;(b)旋转交错网格.
Fig. 1 Diagram of 2D grid discretization
(a)Ordinary staggered-grid(SG);(b)Rotated staggered-grid(RSG).

表 1 常规交错网格参数分布图 Table 1 Parameters distribution on ordinary SG

表 2 旋转交错网格参数分布图 Table 2 Parameters distribution on RSG

在计算空间导数时,RSGFD方法是先沿网格对角线方向进行计算,然后再将其计算结果进行线性叠加,得到沿坐标轴方向的空间导数.以方程(6)中波场p为例,其关于坐标轴方向的一阶空间偏导数的任意2N阶RSGFD格式可表示为(Saenger and Bohlen,2004)

其中,Δx和Δz分别是坐标轴方向的空间步长,Δr是对角线长度,Δr=(Δx2+Δz2)1/2.∂/∂x′和∂/∂z′分别表示沿两个对角线x′和z′方向的任意2N阶偏导数.设波场函数p(x,z,t)连续,则:

式中,Cn(N)为一阶空间偏导算子的差分系数,通常的求解方法是Taylor级数展开法.将式(9)代入式(8)中,即为任意2N阶精度的RSGFD格式,称为TE-RSGFD方法.

然而,对于同样大小的网格单元,RSGFD方法需要更大的空间步长,这往往会增加空间频散.研究表明,在同样网格离散程度和不增加计算量情况下,基于最小二乘优化系数的规则网格有限差分(LS-FD)法可以有效压制频散,提高模拟精度(Liu,2013).因此,将全局优化差分系数方法引入到RSGFD方法中.首先利用平面波理论,令:

式中,p0为常数,kxkz是空间波数.将式(10)代入式(8)中,可得到如下关系式:

式中,α=kxΔx=kzΔz,0≤α≤π.利用式(11)中α与α*之差可以估算出RSGFD方法引入的误差.设目标函数为

式中,Fobj是关于RSGFD系数C(N)n(n= 1,2,…,N)的二次函数,积分上限A≤π.对α在整个范围积分可求得误差,再结合最小二乘方法可以估算该误差的最小值.由极值的必要条件:∂Fobj/∂C(N)n=0,并结合式(12)进行求解,则可以得到RSGFD方法的全局优化新差分系数Cn(N).将新系数代入式(8)中,即为任意2N阶精度的优化LS-RSGFD格式.

以方程(6)中伪应力波场p及其伪速度分量upwp为例,给出TTI介质一阶qP波稳定方程的优化LS-RSGFD格式如下:

其中,i,j为空间离散点序号,k为时间离散点序号,Δt为时间步长.2N是差分阶数,Cn(N)为全局优化的新差分系数.同理,可以导出k时刻G1G2T2的差分格式.

3.2 PML边界条件和稳定性条件

压制人工边界反射是地震波数值模拟的一个关键问题,完全匹配层(PML)边界条件被认为是目前 吸收效果最好的边界条件而得到广泛应用(Berenger,1994). 本文依据PML边界条件的方程分裂思路,推导了适用于TI介质一阶qP波方程的时间域完全匹配层控制方程.以方程(6)中的伪应力场p及其速度分量upwp为例:

式中,pxpz为波场pxz方向的分裂算子;d(x)和d(z)为xz方向的阻尼因子.

RSGFD方法是利用差分算子和Taylor级数展开来近似求解偏导数,并进行时间外推,用离散数值解得到整个区域上的近似解,这样在取高阶差分时不可避免地产生截断误差.随着时间外推,这种误差积累越来越大,从而影响差分格式的精度.

时间二阶和空间2N阶精度的常规SGFD方法(董良国等,2000)和RSGFD方法(Saenger and Bohlen,2004)的稳定性条件分别为

式中,vmax=max{vp0(x,z),vpx(x,z),vpn(x,z)}表模型内的最大速度值,Δh是空间步长.对比两式可知,在相同的差分精度条件下,RSGFD方法比常规SGFD方法的稳定性更高.

4 模型试算与分析 4.1 均匀TTI介质模型

设计一个均匀TTI介质模型,其参数为:vp0= 2500 m·s-1ε=0.24,Δ=0.10,θ=45°,ρ=1 g·cm-3. 网格大小201×201,网格尺寸Δxz=10 m,采样间隔Δt=1 ms,选用主频为25 Hz的雷克子波震源,置于模型中心处.

采用空间十阶和时间二阶的优化LS-RSGFD方法进行正演模拟,分别得到VTI介质弹性波(θ=0°)、VTI介质一阶qP波(θ=0°)、TTI介质弹性波、TTI介质一阶qP波正演模拟在300 ms时刻的波场快照,如图 2所示.由图 2a2b可知,外层传播较快的均为qP波,两者的波阵面具有很好的一致性;图 2a中内层传播较慢的是qSV波,而图 2b中内层的菱状波形是声近似特殊情况下产生的伪横波.同样,图 2d2e中的TTI介质情形也能看到相同的波形现象.对比图 2b2c图 2e2f可以看出,在激发震源周围设置小范围的椭圆各向异性区域(即震源环)后,震源处产生的伪横波得到了很好的压制,则能保证VTI和TTI介质中仅含有qP波的波场传播.由此可知,应用VTI一阶qP波方程式(1b)和TTI一阶qP波稳定方程式(6)能够比较准确地模拟各向异性介质中qP波的运动学特征.需注意,下文中给出的各种波动方程的波场快照和地震记录均表示与图 2相同的波场分量.

图 2 VTI(a—c)和TTI(d—f)介质的波场快照对比(t=300 ms)
(a)和(d):弹性波正演的垂直应力分量;(b)和(e):一阶qP波正演的p波场(设置震源环前);(c)和(f):一阶qP波正演的p波场(设置震源环后).
Fig. 2 The wave-field snapshots comparison of VTI(a—c)and TTI(d—f)media(t=300 ms)
(a), (d): Vertical stress component of elastic wave modeling;(b), (e): Wave-field p of first-order qP-wave modeling without source box;(c), (f): Wave-field p of first-order qP-wave modeling with source box.

图 3为应用PML边界条件后的TTI介质qP波波场快照,从左至右的时间依次为300~700 ms,间隔为100 ms.由图可知,在进行各向异性TI介质qP波正演模拟时,利用本文推导的PML边界控制方程能获得理想的边界吸收效果.

图 3 应用PML边界条件后的qP波波场快照Fig. 3 The qP-wave wave-field snapshots with PML boundary condition

为更好地验证LS-RSGFD方法的精确性,首先将震源主频设置为30 Hz,其他模拟参数同图 2.然后利用不同空间差分阶数的常规TE-RSGFD和优化LS-RSGFD方法进行测试,得到如图 4所示的TTI介质qP波波场快照.对比图可知,对于相同的空间差分阶数,优化方法数值频散更弱,模拟精度更高(红色箭头所示).

图 4 常规TE-RSGFD(a—d)和优化LS-RSGFD(e—h)方法得到的波场快照对比(t=300 ms)
(a)和(e)8阶;(b)和(f)10阶;(c)和(g)12阶;(d)和(h)14阶.
Fig. 4 The wave-field snapshots comparison of conventional TE-RSGFD(a—d)and optimal LS-RSGFD(e—h)method(t=300 ms)
(a), (e)8-order; (b), (f)10-order; (c), (g)12-order; (d), (h)14-order.

图 5a是从图 4中地面位置600 m处抽取的单道波形对比结果,同样可以看出,对于相同的空间差分阶数,优化方法的精度更高.而且8阶优化方法与12阶常规方法的数值频散相当,10阶优化方法比14阶常规方法的数值频散更弱,表明用相对低阶的优化方法就能达到高阶常规方法的模拟精度,甚至更高.

图 5 不同差分阶数的TE-RSGFD和LS-RSGFD方法的(a)单道波形对比和(b)计算时间对比Fig. 5 The single trace wave-form comparison(a)and the computing time comparison(b)of TE-RSGFD and LS-RSGFD method with different difference orders

图 5b显示了不同空间差分阶数对应的TE-RSGFD方法(灰色)和LS-RSGFD方法(蓝色)的CPU平均计算时间对比结果(正演模拟参数同图 2,记录总时间为4.0 s).虽然计算机硬件不同可能会导致计算时间有所差异,但是在同等条件下得到的计算时间还是具有一定的参考价值.由图 5b可知,相同空间差分阶数条件下,优化方法的计算时间低于常规方法,前者相对于后者效率提高了约4%~15%.由于同一方法的计算时间随着差分阶数的增大而增加,因此,对于给定的模拟精度,利用低阶的优化方法还可以大大地提高计算效率.以8阶优化方法和12阶常规方法为例,两者的精度相当,但是前者相对于后者效率提高了约30%.

4.2 TTI介质层状模型及各向异性波场特征分析

图 6a为网格大小301×151的TTI介质层状模型及其各向异性参数,网格间距为10 m×10 m,采用空间十阶和时间二阶的LS-RSGFD方法进行数值模拟.模拟的观测系统为:选用主频为25 Hz的雷克子波震源,置于模型中(x=1500 m,z=10 m)位置处,采用中间放炮,两边接收的方式,共301道,道间距为10 m,采样间隔Δt=1.0 ms,记录时间为1.5 s.

图 6 TTI介质层状模型及波场快照(t=580 ms)
(a)各向异性参数;(b)VTI方法;(c)TTI方法.
Fig. 6 TTI layered model and snapshots(t=580 ms)
(a)Anisotropic parameters;(b)VTI method;(c)TTI method.

图 6b6c分别为VTI方法和TTI方法得到的t=580 ms时刻的波场快照.图 7分别显示了各向同性介质声波(图 7a)、VTI介质一阶qP波(图 7b)、TTI介质一阶qP波(图 7c)、VTI介质弹性波(图 7d)正演模拟得到的单炮地震记录.由图 6图 7可以简单看出,TTI介质的倾角因素会造成波前面形态发生扭曲,如图 6中深度1000~1500 m区域所示,同时各向异性因素会改变走时等运动学特征,如图 7中时间600~1100 ms区域所示.

图 7 合成地震记录
(a)各向同性声波正演;(b)VTI一阶qP波正演;(c)TTI一阶qP波正演;(d)VTI弹性波正演.
Fig. 7 Synthetic seismic recordings
(a)Isotropic acoustic wave modeling; (b)VTI first-order qP-wave modeling; (c)TTI first-order qP-wave modeling; (d)VTI elastic wave modeling.

为了更清晰地对比波场特征的变化,从图 7b7d中,以及图 7a7b7c中小偏移距处(Distance=1200 m)和大偏移距处(Distance=200 m)分别抽取一条地震道,得到如图 8图 9所示的单道波形对比图.对比图 8中蓝线和绿线可知,在两个反射层上,VTI弹性波与一阶qP波正演模拟的qP波反射走时吻合较好,表明本文的TI介质一阶qP波方程能够比较准确地模拟各向异性介质中qP的运动学(走时和速度等)特征,而动力学(振幅等)特征上也是TI介质弹性波方程的有效近似.

图 8 VTI介质弹性波与VTI介质qP波正演模拟的单道波形对比图
(a)小偏移距;(b)大偏移距.
Fig. 8 Single trace waveform comparison of VTI elastic and qP-wave modeling
(a)Small;(b)Large offset.

图 9 各向同性声波、VTI介质qP波、TTI介质qP波正演的单道波形对比图
(a)小偏移距;(b)大偏移距.
Fig. 9 Single trace waveform comparison of Isotropic, VTI, TTI qP-wave modeling
(a)Small ;(b)Large offset.

对比图 9中红线和绿线可知,相比于各向同性声波正演模拟,考虑各向异性因素的影响后,各层的反射走时和波形振幅均发生变化,特别是大偏移距处,这种影响更明显,且往往造成走时和振幅减小.同时对比图 9中绿线和蓝线可知,模型中第二层受到TTI介质的倾角影响,造成波前形态发生扭曲,使得振幅和走时等波场特征会进一步改变.因此,在针对实际资料中大倾角各向异性构造处理时,必须要考虑到TTI介质的倾角因素.

4.3 VTI介质Hess模型

为了进一步验证方法的有效性,采用空间十阶和时间二阶的LS-RSGFD方法对复杂VTI介质Hess模型进行一阶qP波数值模拟.模型网格大小为601×501,网格间距为10 m×10 m,其各向异性参数vp0ε、δ图 10所示,密度ρ=1 g·cm-3.模拟的观测系统为:选用主频为30 Hz的雷克子波震源,置于模型中(x=4000 m,z=10 m)位置处,采用中间放炮,两边接收的方式,共401道,道间距为10 m,采样间隔Δt=1.0 ms,记录时间为3.5 s.

图 10 VTI介质(a); HESS模型(b)及参数(c)Fig. 10 VTI-HESS model and anisotropic parameters

图 11是应用VTI一阶qP波方程式(1b)得到的qP波波场快照(t=1.3 s时刻)和单炮合成地震记录(局部放大结果),图 11a11c图 11b11d分别对应常规的TE-RSGFD方法和优化的LS-RSGFD方法.对比图 11a11b中的波场快照可以明显看出,常规TE-RSGFD方法得到的波场快照中存在明显的数值频散和波形扰动,而优化LS-RSGFD方法得到的波场快照中几乎没有频散,波前形态更清晰(黑色实线框区域所示).同时,对比图 11c11d中的炮记录也可以看出,常规TE-RSGFD方法得到的炮记录中频散严重,造成同相轴扰动,甚至导致部分区域的同相轴模糊不清(图 11c中黑色箭头所示);然而,优化LS-RSGFD方法得到的炮记录中频散现象得到很好地压制,其同相轴更清晰,连续性更好.因此,在相同的数值离散条件下,利用优化LS-RSGFD方法能够有效地压制数值频散,提高地震记录的模拟精度.

图 11 VTI介质Hess模型的qP波波场快照(t=1.3 s)及合成地震记录Fig. 11 The qP-wave wave-field snapshots(t=1.3 s)and seismic recordingds and of VTI Hess model
4.4 TTI介质逆掩断层模型

最后,利用空间十阶和时间二阶的LS-RSGFD方法对TTI介质逆掩断层模型进行正演模拟,模型的各向异性参数vp0ε、δ及倾角θ,如图 12a2b所示,密度ρ=1 g·cm-3;模型网格大小为401×401,网格间距为10 m×10 m.正演模拟的观测系统为:选用主频为30 Hz的雷克子波震源,置于模型中(x=2000 m,z=10 m)位置处,采用中间放炮,两边接收的方式,共401道,道间距为10 m,采样间隔Δt=1.0 ms,记录时间为2.5 s.

图 12 TTI逆掩断层模型及qP波波场快照(t=900 ms)
(a)vp0ε,δ;(b)θ;(c)不稳定方程;(d)稳定方程.
Fig. 12 TTI Over-thrust model and the qP-wave wave-field snapshots(t=900 ms)
(a)vp0, ε, δ;(b)θ;(c)Unstable equation;(d)Stable equation.

图 12c12d分别是应用稳定的TTI介质一阶qP波方程式(6)(设置控制参数σ=0.4)和不稳定的TTI介质一阶qP波方程式(2)得到的qP波波场快照(t=900 ms时刻).由图 12c可知,倾角变化引起的波场不稳定性严重干扰了有效的qP波能量,不利于得到准确的qP波合成地震记录;而图 12d中波场稳定传播,波场快照更加清晰.

图 13是应用TTI一阶qP波稳定方程式(6)得到的单炮地震记录,图 13a13c分别对应常规的TE-RSGFD方法和优化的LS-RSGFD方法,为了更清晰地对比,分别将两图中红色实线框区域进行 局部放大,得到图 13b13d的结果.由局部放大图可以明显看到,TE-RSGFD方法得到的地震记录中存在明显的数值频散和同相轴扰动现象,造成下方的同相轴变得模糊,连续性不好(图 13b中红色箭头所示);而LS-RSGFD方法得到的地震记录中同相轴更加清晰,连续性更好(图 13d中黑色箭头所示).因此,在相同的数值离散条件下,利用优化的LS-RSGFD方法能够得到更高精度的合成地震记录.

图 13 单炮地震记录
(a)常规的TE-RSGFD方法;(b)13a中虚线框区域局部放大;(c)优化的LS-RSGFD方法;(d)13c中虚线框区域局部放大.
Fig. 13 Single shot seismic recordings
(a)ConventionalTE-RSGFD method;(b)Local zoom of the dotted line box in 13a;(c)Optimal LS-RSGFD method;(d)Local zoom of the dotted line box in 13c.
5 结论与认识

本文基于声学近似的TI介质弹性波方程与频散关系两种思路,首先研究了VTI介质一阶qP波方程和稳定的TTI介质一阶qP波方程,并推导了相应的PML边界条件.然后,结合RSGFD法和LS-FD法的各自优势,研究了一种基于最小二乘优化差分系数的旋转交错网格有限差分(LS-RSGFD)方法,并通过构造TI介质一阶qP波方程的LS-RSGFD格式,实现了高精度的各向异性介质qP波波场数值模拟.

数值结果表明,TI介质一阶qP波方程数值模拟克服了弹性波数值模拟中波场分离困难的局限性,能够有效地模拟各向异性介质中qP波的运动学特征,便于研究qP波的传播规律,并为后续的各向异性纵波资料处理、解释和反演奠定理论基础.波场特征分析验证了各向异性因素会对地震波场特征(走时和振幅)产生一定影响,特别是远偏移距处这种影响更加明显,进而还会影响偏移和反演的精度.因此,在针对大偏移距、宽方位地震资料的实际应用中,各向异性因素越来越不能忽略.

在数值实现方面,相比于常规的TE-RSGFD方法,优化的LS-RSGFD方法压制数值频散效果更好,模拟精度更高;并且对于给定的模拟精度,利用LS-RSGFD方法还可以很大程度地提高计算效率.同时,在波场数值模拟过程中,利用PML边界条件可以实现对边界反射能量的有效吸收.复杂VTI介质Hess模型和TTI介质逆掩断层模型的数值测试结果也验证了LS-RSGFD方法能够得到高精度的qP波合成地震记录,引入各向异性控制参数σ可以有效地减弱不稳定性问题,并保证倾角突变的TTI介质中qP波场的稳定传播.利用稳定的TTI介质一阶qP波有限差分算子进行精确的偏移成像将是下一步的研究方向.

致谢 衷心感谢审稿专家和编辑部提出的宝贵意见和建议.感谢中国石油大学(华东)地震波传播与成像(SWPI)课题组的帮助,感谢国家863课题和国家自然科学基金的联合资助.
参考文献
[1] Alkhalifah T. 2000. An acoustic wave equation for anisotropic media. Geophysics, 65(4): 1239-1250.
[2] Berenger J P. 1994. A perfectly matched layer for the absorption of electromagnetic waves. Journal of Computational Physics, 114(2): 185-220.
[3] 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.
[4] Chu C L, Macy B K, Anno P D. 2011. Approximation of pure acoustic seismic wave propagation in TTI media. Geophysics, 76(5): WB97-WB107.
[5] Dong L G, Ma Z T, Cao J Z. 2000. A study on stability of the staggered-grid high-order difference method of one-order elastic wave equation. Chinese J. Geophys. (in Chinese), 43(6): 856-864.
[6] Du X, Bancroft J C, Lines L R. 2007. Anisotropic reverse-time migration for tilted TI media. Geophysical Prospecting, 55(6): 853-869.
[7] Du X, Fletcher R P, Fowler P J. 2008. A new pseudo-acoustic wave equation for VTI media.// 70th Annual International Conference and Exhibition, EAGE, Extended Abstracts, H033.
[8] Duveneck E, Milcik P, Bakker P M. 2008. Acoustic VTI wave equations and their application for anisotropic reverse-time migration.// SEG Technical Program Expanded Abstracts, 2186-2190.
[9] Duveneck E, Bakker P M. 2011. Stable P-wave modeling for reverse-time migration in tilted TI media. Geophysics, 76(2): S65-S75.
[10] Fletcher R, Du X, Fowler P J. 2008. A new pseudo-acoustic wave equation for TI media.// 78th Annual International Meeting, SEG, Expanded Abstracts, 2082-3086.
[11] Fletcher R P, Du X, Fowler P J. 2009. Reverse time migration in tilted transversely isotropic (TTI) media. Geophysics, 74(6): WCA179-WCA187.
[12] Fowler P J, Du X, Fletcher R P. 2010. Coupled equations for reverse time migration in transversely isotropic media. Geophysics, 75(1): S11-S22.
[13] Grechka V, Zhang L B, Rector J W. 2004. Shear waves in acoustic anisotropic media. Geophysics, 69(2): 576-582.
[14] Han L H, He B S. 2010. First-order quasi-P wave equation forward modeling and boundary conditions in VTI medium. Oil Geophysical Prospecting (in Chinese), 45(6): 819-825.
[15] Hestholm S. 2009. Acoustic VTI modeling using high-order finite differences. Geophysics, 74(5): T67-T73.
[16] Klíe H, Toro W. 2001. A new acoustic wave equation for modeling in anisotropic media.// 71st Annual International Meeting, SEG, Expanded Abstracts. 1171-1174.
[17] Li Z C. 2011. Seismic Pre-Stack Imaging Theory and Method (in Chinese). Dongying: China University of Petroleum Press, 24-39.
[18] Liang K, Yin X Y, Wu G C. 2011. Exact and approximate reflection and transmission coefficient for incident qP wave in TTI media. Chinese J. Geophys. (in Chinese), 54(1): 208-217, doi: 10.3969/j.issn.0001-5733.2011.01.022.
[19] 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.
[20] Liu Y. 2013. Globally optimal finite-difference schemes based on least squares. Geophysics, 78(4): T113-T132.
[21] Pei Z L. 2004. Three-dimensional numerical simulation of elastic wave propagation in 3-D anisotropic media with staggered-grid high-order difference method. Journal of the University of Petroleum, China (in Chinese), 28(5): 23-29.
[22] Pestana R C, Ursin B, Stoffa P L. 2011. Separate P- and SV-wave equations for VTI media.// 81st Annual International Meeting, SEG, Expanded Abstracts, 163-167.
[23] 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.
[24] Tessmer K. 2010. Reverse-time migration in TTI media.// 80th Annual International Meeting, SEG, Expanded Abstracts, 3193-3197.
[25] Thomsen L. 1986. Weak elastic anisotropy. Geophysics, 51(10): 1954-1966.
[26] Tsvankin I. 2001. Seismic Signatures and Analysis of Reflection Data in Anisotropic Media. Pergamon, 1-60.
[27] Xu S, Zhou H B. 2014. Accurate simulations of pure quasi-P-waves in complex anisotropic media. Geophysics, 79(6): T341-T348.
[28] Yang F S, Li Z C, Wang X D, et al. 2015. Equivalent anisotropic first-order pseudo-acoustic wave equations modeling and its efficiency analysis. Acta Seismologica Sinica (in Chinese), 37(4): 648-660.
[29] Yin X Y, Pu Y T, Liang K, et al. 2014. Rotated staggered finite-difference modeling for quasi-P wave in VTI medium. CT Theory and Applications (in Chinese), 23(5): 771-783.
[30] Yoon K, Suh S, Ji J, et al. 2010. Stability and speedup issues in TTI RTM implementation.// 80th Annual International Meeting, SEG, Expanded Abstracts, 3221-3225.
[31] 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.
[32] Zhang H J, Zhang Y. 2008. Reverse time migration in 3D heterogeneous TTI media.// SEG Technical Program Expanded Abstracts. Las Vegas, Nevada: SEG, 2196-2200.
[33] Zhang J H, Zhang G Q, Zhang Y. 2009. Removing S-wave noise in TTI reverse time migration. //79th Annual International Meeting. SEG, Expanded Abstracts. Houston, Texas: SEG, 2849-2853.
[34] Zhang L B, Rector J W III, Hoversten G M. 2005. Finite-difference modelling of wave propagation in acoustic tilted TI media. Geophysical Prospecting, 53(6): 843-852.
[35] Zhang Y, Zhang H Z, Zhang G Q. 2011. A Stable TTI Reverse time migration and its implementation. Geophysics, 76(3): WA3-WA11.
[36] Zhang Y, Wu G C. 2013. Methods of removing pseudo SV-wave artifacts in TTI media qP-wave reverse-time migration. Chinese Journal of Geophysics (in Chinese), 56(6): 2065-2076, doi: 10.6038/cjg20130627.
[37] Zhou H, Zhang G, Bloor R. 2006a. An anisotropic acoustic wave equation for VTI media.// 68th EAGE Conference and Exhibition. SPE, EAGE.
[38] Zhou H B, Bloor R, Zhang G Q. 2006b. An anisotropic acoustic wave equation for modeling and migration in 2D TTI media.// 76th Annual International Meeting. New Orleans, Louisiana: SEG, 194-198.
[39] 程玖兵, 康玮, 王腾飞. 2013. 各向异性介质qP波传播描述I: 伪纯模式波动方程. 地球物理学报, 56(10): 3474-3486, doi: 10.6038/cjg20131022.
[40] 董良国, 马在田, 曹景忠. 2000. 一阶弹性波方程交错网格高阶差分解法稳定性研究. 地球物理学报, 43(6): 856-864.
[41] 韩令贺, 何兵寿. 2010. VTI介质一阶准P波方程正演模拟及边界条件. 石油地球物理勘探, 45(6): 819-825.
[42] 李振春. 2011. 地震叠前成像理论与方法. 东营: 中国石油大学出版社, 24-39.
[43] 梁锴, 印兴耀, 吴国忱. 2011. TTI介质qP波入射精确和近似反射透射系数. 地球物理学报, 54(1): 208-217, doi: 10.3969/j.issn.0001-5733.2011.01.022.
[44] 裴正林. 2004. 三维各向异性介质中弹性波方程交错网格高阶有限差分法数值模拟. 石油大学学报(自然科学版), 28(5): 23-29.
[45] 杨富森, 李振春, 王小丹等. 2015. 等价各向异性介质一阶拟声波方程正演模拟及其效率分析. 地震学报, 37(4): 648-660.
[46] 印兴耀, 蒲义涛, 梁锴等. 2014. VTI介质准P波旋转交错有限差分数值模拟. CT理论与应用研究, 23(5): 771-783.
[47] 张岩, 吴国忱. 2013. TTI介质qP波逆时偏移中伪横波噪声压制方法. 地球物理学报, 56(6): 2065-2076, doi: 10.6038/cjg20130627.