VTI介质修正声学近似qP波波动方程与模拟

梁 锴,陈浩然,孙上饶

(中国石油大学(华东) 地球科学与技术学院,山东 青岛 266580)

各向异性介质中纵横波通常是耦合在一起传播的,因此各向异性弹性波波场解耦在地震成像和全波形反演中十分重要[1-2]。对应具有垂直对称轴的横向各向同性(transverse isotropy with a vertical axis of symmetry,VTI)介质是实际地质条件的有效简化,得到了广泛的关注。Alkhalifah[3]提出VTI介质声学近似方法,即通过设置垂向准横波(qSV波)速度(VS0)为0来解耦qP波场。然而,在声学假设下,准纵波(qP波)方程有两组共轭解,分别对应qP波和退化qSV波。后者是不需要的干扰波,并且各向异性参数ε<δ时,该方法随着迭代次数的增加会出现严重的数值溢出现象。

为了克服退化qSV波的干扰以及ε<δ时的不稳定现象,许多学者进行了各向异性介质改进或优化的纯qP波方程的研究[4-11]。然而,这些方程大多含有分数阶非椭圆项,导致模拟的计算成本较高。利用有限差分(finite difference,FD)方法求解这类纯波波动方程具有高效、适应复杂介质等特点[12],相比之下利用伪谱(pseudo-spectral,PS)法求解这类方程空间导数的精度更高[13],但对于复杂介质则会出现严重的数值发散。为了降低计算成本、提高计算精度,有限差分和伪谱混合的正演方法[14]以及低秩近似方法[15]等多种数值模拟方法相继被提出。

本研究将垂向qSV波速度VS0考虑成波数kx、kz和各向异性参数ε、δ的函数,推导了VTI介质修正声学近似qP波频散关系,通过傅里叶反变换得到修正声学近似的qP波波动方程,然后采用混合有限差分/伪谱算法实现了数值模拟,最后进行了精度和近似的频散关系分析,并给出了均匀介质模型和复杂介质模型的正演模拟数值示例。

1 VTI介质修正声学近似qP波波动方程

各向异性介质的Christoffel方程通常用于各向异性介质弹性波相速度或频散关系表征。根据VTI介质Christoffel方程可以确定VTI介质qP波和qSV波精确频散关系为:

(1)

(2)

(3)

其中:ωP、ωSV分别为qP波、qSV波的圆频率;kx=ksinθ和kz=kcosθ分别为沿x和z方向的波数,k为波数;θ为传播方向与z轴的夹角;VP0和VS0分别为qP波和qSV波沿VTI介质对称轴方向(z轴方向)的相速度;ε和δ为Thomsen参数[16]。

针对精确频散关系表达式复杂、难以应用的问题,Alkhalifah[3]提出了著名的声学假设近似,令VS0=0(qSV波沿VTI介质对称轴相速度为0),得到qP波声学近似频散关系方程:

(4)

其中,ωPa为声学假设近似qP波圆频率。传统声学近似频散关系形式简单、精度较高,但是对应的波动方程存在退化qSV波等问题。本研究对该方法进行了修正,不是令VS0直接为0,而是将VS0考虑成波数kx、kz和各向异性参数ε、δ的函数。首先根据式(2)得到VS0和ωSV的关系为:

(5)

将式(5)代入式(1)可得:

(6)

将VS0进行如下近似:

(7)

将式(7)代入式(5),此时ωSV可表示为:

(8)

将式(8)代入式(6),可以得到VTI修正声学近似qP波频散关系:

(9)

其中,ωPma为修正声学假设近似的qP波圆频率。

(10)

其中:P=P(x,z,t)为时空域的qP波波场;F、F-1分别表示傅里叶正变换和反变换。

VTI介质修正声学近似qP波波动方程(10)包含椭圆项和非椭圆项两部分。当ε=δ,即椭圆各向异性介质时,非椭圆项为0,该方程退化为椭圆波动方程,说明该方程在椭圆各向异性介质中是完全精确的。

2 混合有限差分/伪谱法正演模拟算法

对VTI介质修正声学近似qP波波动方程(10)进行数值模拟,由于该方程包含椭圆项和非椭圆项两部分,因此采用混合有限差分/伪谱法求解,即椭圆项采用有限差分算法求解,而非椭圆项包含波数kx、kz的高阶项,不适合直接采用有限差分算法求解,则采用伪谱法进行求解。

根据有限差分基本原理,将空间导数和时间导数用差分近似表征为:

(11)

对VTI介质修正声学近似qP波波动方程采用混合有限差分/伪谱法求解,最终的差分格式为:

(12)

其中,G表示用伪谱法计算的非椭圆项,即

(13)

3 数值示例

3.1 频散关系分析

为了验证VTI介质修正声学近似qP波波动方程的正确性和适用性,进行近似qP波频散关系的精度分析。设计了两个均匀VTI介质模型A和B,其相同参数为VP0=3 000 m/s,VS0=1 500 m/s,其各向异性参数ε、δ不同,如表1所示。其中,模型A中ε>δ,模型B中ε<δ。

表1 模型各向异性参数

利用式(1)、式(4)和式(9)分别计算精确qP波频散关系ωP、声学近似qP波频散关系ωPa和修正声学近似频散关系ωPma及其相对误差,结果如图1和图2所示。其中,图1(a)和图2(a)为频散关系曲线,图1(b)和图2(b)为相对误差曲线,蓝色实线为精确qP波频散关系,红色虚线为声学近似qP波频散关系及其误差,黑色虚线为修正声学近似频散关系及其误差。由图1和图2发现,声学近似和修正声学近似的频散关系与精确频散关系数值十分接近,本例中相对误差均小于0.25%,说明两种近似均具有很高的精度;并且对于ε>δ的VTI介质和ε<δ的VTI介质,修正声学近似频散关系的最大相对误差数值以及平均相对误差数值均小于声学近似的误差,说明修正声学近似的整体精度高于声学近似的精度。

图1 模型A精确和近似频散关系及其相对误差

3.2 正演模拟

为了进行波场对比,利用有限差分法求解VTI介质弹性波波动方程和声学近似qP波波动方程,利用有限差分和伪谱混合法求解修正声学近似qP波波动方程式(式(12)和式(13))。首先考虑表1中的两种均匀VTI介质模型,网格大小为301×301,横向和纵向网格间距均为10 m,采样间隔为1 ms,震源位于模型中心,采用主频20 Hz的雷克子波。

波场快照如图3和图4所示,其中图3(a)和图4(a)为弹性波波场,图3(b)和图4(b)为基于声学近似模拟的波场,图3(c)和图4(c)为基于修正声学近似模拟的qP波波场。在图3中,声学近似和修正声学近似模拟的qP波波场与弹性波模拟的qP波波场均吻合较好,说明两种方法均适用于ε>δ的模型。但是声学近似模拟的波场中存在退化的qSV波波场,说明声学近似模拟的波场不是纯qP波;修正声学近似模拟的波场中不存在退化的qSV波波场,说明修正声学近似模拟的波场是纯qP波。在图4中,声学近似模拟的波场存在数值溢出,使模拟的qP波波场淹没在数值溢出的干扰中,即声学近似无法适用于ε<δ的模型。而修正声学近似模拟的qP波波场与弹性波模拟的qP波波场也吻合较好,且不存在数值溢出现象。该数值示例表明修正声学近似qP波波动方程刻画了纯qP波波场,且既适用于ε≥δ的VTI介质,也适用于ε<δ的VTI介质。

图3 模型A的波场快照

图4 模型B的波场快照

网格点(130,120)处的地震记录如图5所示,其中图5(a)为模型A模拟结果,图5(b)为模型B模拟结果,蓝色实线为弹性波记录,红色点虚线为声学近似记录,黑色虚线为修正声学近似记录。由图5可见,弹性波记录中包含qP波和qSV波,声学近似记录中包含qP波和退化qSV波,而修正声学近似记录中只包含qP波,并且与弹性波模拟中的qP波记录吻合较好。

图5 网格点(130,120)处地震记录

将本研究方法应用于复杂介质模型(BP模型)。网格大小为601×601,横向和纵向网格间距均为10 m,采样间隔为1 ms。震源位于图6星号所示处,采用主频20 Hz的雷克子波。BP模型如图6所示,图6(a)为VP0模型,图6(b)为ε模型,图6(c)为δ模型。波场快照如图7所示,其中图7(a)为弹性波波场,图7(b)为声学近似模拟的波场,图7(c)为基于修正声学近似模拟的qP波波场。图7表明,复杂介质中声学近似模拟的qP波波场与弹性波模拟的qP波波场运动学规律一致,波前面吻合较好,但是声学近似模拟波场中可能存在退化qSV波(图7(b)中白色箭头所示),会增加纯qP波处理和解释的难度。而修正声学近似模拟的qP波波场的运动学规律与弹性波模拟的qP波波场一致,二者波前面吻合较好,且该方法模拟结果中不存在退化qSV波,说明该方法在复杂VTI介质中也能很好地实现纯qP波正演模拟。图8为震源深度处的单炮记录,其中图8(a)为弹性波记录,图8(b)为声学近似记录,图8(c)为修正声学近似记录。由图8可以看出,弹性波记录中qSV波比较强;声学近似记录中qP波比较明显,但是仍然存在退化qSV波的干扰;修正声学近似记录中只有qP波,并且与弹性波模拟中的qP波记录基本一致。

图6 BP模型

图7 波场快照

图8 单炮记录

4 结论

本研究对声学近似进行了修正,推导了VTI介质修正声学近似qP波频散关系和波动方程,采用混合有限差分/伪谱算法实现了VTI介质纯qP波正演模拟。通过理论分析和数值示例得到如下结论。

1) 基于修正声学近似的qP波波动方程是关于时间二阶导数的偏微分方程,其一组共轭解对应发散和汇聚的qP波,因此该方程的解不包含退化qSV波,而是纯qP波方程。

2) 基于修正声学近似的qP波波动方程包含椭圆项和非椭圆项两部分,非椭圆项在椭圆各向异性介质中等于零,因此该方程在椭圆各向异性介质中是完全精确的。

3) 因为退化qSV波在ε<δ的VTI介质中存在数值溢出,导致模拟结果不稳定,而基于修正声学近似的qP波波动方程不包含退化qSV波,因此该方程在ε≥δ和ε<δ的VTI介质中均是稳定的,能够适用于复杂VTI介质。

4) 直接用有限差分法求解波动方程中高阶导数是比较困难的,而混合有限差分/伪谱算法能够较好地解决波数kx、kz高阶项的计算问题。