TTI介质弹性波近似解耦波动方程

梁锴， 孙上饶， 曹丹平， 印兴耀

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

0 引言

TTI介质弹性波频散关系方程可以通过求解Christoffel方程得到，利用不同的近似方法对其进行近似处理，并利用傅里叶逆变换将频率-波数域算子变换到时空域，可以得到对应的解耦波动方程.Alkhalifah和Tsvankin(1995)、Alkhalifah(1998, 2000)开创性地提出了声学假设近似方法，即将沿对称轴的横波速度近似视做零值，推导了一种简化的VTI介质弹性波频散关系方程，进而得到qP波方程.多种纯qP方程相继被提出(Klíe and Toro, 2001; Zhou et al., 2006; Hestholm, 2007; Du et al., 2008, 2010; Duveneck et al., 2008; Chu et al., 2011, 2013; Bloot et al., 2013; Schleicher and Costa, 2016; Li and Zhu, 2018; Xu et al., 2020)，但是声学假设条件下的波动方程存在两组共轭解，分别对应qP波和退化qSV波，而后者是不被期望的“干扰”波，在地震波数值模拟的过程中，随着迭代次数的增加，退化qSV波会产生严重的数值频散问题，特别在各向异性参数ε<δ的情况下，描述退化qSV波的解发散，误差随时间累计，导致数值解不稳定.Jin和Stovas(2018)定义一组新的参数用于描述VTI介质中退化qSV波运动学特征，并分析了其传播特征.为了能完全消除声学近似中的qSV波，Klíe和Toro(2001)在此基础上进行改进，消除了描述退化qSV波的一组解析解，得到纯qP波波动方程.此外也可以通过引入辅助变量，实现对VTI介质qP波方程的降阶处理，使其更容易实现(Zhou et al., 2006; Du et al., 2008; Chu et al., 2011).Liu等(2009)基于VTI介质声学近似方程，解耦得到了一种稳定的纯qP波和退化qSV波的波动方程，并应用于逆时偏移中；Fletcher等(2009)基于坐标旋转法推导了一般TTI介质qP波、qSV波的波动方程，并指出若人为设置TTI介质对称轴方向的横波速度足够大，则可以移除qSV的波面三角区，进而使波场能稳定传播；Chu等(2011)对Fletcher推导的TTI介质拟声波方程进行因式分解，利用一阶和高阶Taylor展开近似qP波方程，并指出采用一阶Taylor展开近似的qP波方程足以准确运用于大部分实际应用中；梁锴等(2009)从Thomsen弱各向异性近似和声学假设近似出发，对TTI介质弹性波波动方程进行了分解，导出TTI介质弱各向异性条件近似qP波和qSV波波动方程以及声学假设近似的qP波波动方程；黄翼坚等(2011)对VTI介质弹性波精确相速度表达式进行多项式近似，得到了VTI介质纯qP波方程；杨鹏等(2017)利用近似展开法，通过分解偏微分算子得到了TI介质纯qP波动方程；张庆朝等(2019)根据弱各向异性近似的假设，推导了任意空间取向TI介质弹性波相速度近似表达式，并进一步导出了qP波和qSV波的解耦波动方程，波场模拟结果显示该波动方程有着较好的稳定性；Chu等(2013)利用伪谱法对TTI介质拟声波方程(Chu et al., 2011)进行了数值模拟，误差分析表明一阶Taylor近似适合较弱的各向异性介质，强各向异性介质中需要采用更高阶的近似，以控制误差在允许范围内；杜启振等(2015)采用伪谱法和有限差分法相结合的混合法求解VTI介质纯qP方程，极大地提升了运算效率.

本文在TTI介质弹性波精确频散关系的基础上，利用近似的配方法(梁锴等，2018；孙上饶等，2021)，推导TTI介质qP波和qSV波近似频散关系，并将频率-波数域算子反变换到时-空域，导出了qP波和qSV波近似解耦的波动方程.使用两组模型参数对近似频散关系方程进行数值计算，得到了三维频散关系曲面和二维曲线，验证了近似频散关系的有效性.随后考察了近似式与精确式的差异项数值在不同各向异性强度下的变化情况，并分析了XOZ面内近似频散关系曲线的相对误差分布.最后，使用有限差分方法求解TTI介质纯qP波和纯qSV波近似波动方程，模拟了qP波和qSV波在均匀、层状及复杂TTI介质中的传播.数值模拟结果表明，在η<0和各向异性倾角变化较大的介质中，纯qP波和纯qSV波近似波动方程依然可以保持稳定.

1 方法原理

1.1 TTI介质弹性波近似频散关系方程

TTI介质弹性波相速度和频散关系可通过将平面波位移方程代入到Christoffel方程求解得到.TTI介质qP波和qSV波精确频散关系方程可表示为：

(1)

图1 TTI介质模型示意图Fig.1 Schematic diagram of TTI media model

依照近似配方法的思想(梁锴等, 2018)消除D项的一次根号，简化频散关系方程(1)，其中式(1)中D项可改写为：

(2)

(3)

(4)

令式(4)中右边第二项近似为零，即：

(5)

那么，D项可近似表示为D≈(a+b)2=Da.将D的近似值Da代入到精确qP波和qSV波频散关系(1)中，即可得到相对应的基于近似配方法的频散关系方程：

(6)

(7)

1.2 近似纯qP波和纯qSV波波动方程

(8)

(9)

解耦波动方程(8)和(9)既包含对时间和空间的混合偏导，也存在混合空间导数，这使得其计算量大于单一变量的空间导数(Fletcher et al., 2009)，不利于波动方程的数值求解.考虑波动方程(8)的二维形式，并引入辅助变量QP，变换得到等效的纯qP波波动方程为：

(10)

(11)

同理引入辅助变量QSV，对方程(9)进行变换，得到等效qSV波波动方程为：

(12)

(13)

可以利用方程(10)—(13)分别进行纯qP波和纯qSV波的数值模拟.

2 数值示例

2.1 数值试算和误差分析

Vernik和Liu(1997)对多组含油气地层进行了地震波相速度和各向异性参数的测量，本文从测量数据中选取两组中强各向异性页岩作为TTI介质模型，参数如表1所示.

表1 TTI介质模型参数Table 1 Parameters of TTI models

通过对比qP波和qSV波的精确和近似频散关系中的D项，定义差异项ΔD为：

(14)

利用表1的model 1的速度参数和倾角θ0对ΔD进行数值计算，在|σ|值固定不变的情况下，仅改变参数η和δ的值，结果如图2所示，图中空心箭头处为沿对称轴方向传播，实心箭头处为垂直对称轴方向传播.参数|η|越小，差异ΔD越小，并且在对称轴传播方向附近的误差相对垂直对称轴传播方向小；在沿对称轴和垂直对称轴的传播方向，近似值与精确值相等.特别地，当η=0时，介质呈现椭圆各向异性，近似频散关系方程(5)等于精确式(1).

图2 ΔD随相角θ的变化特征(a) η>0； (b) η<0.Fig.2 Variation characteristics of the difference ΔD with phase angle θ

对TTI介质qP波和qSV波频散关系进行三维数值计算，结果如图3所示.图3a、b为model 1的精确qP波和qSV波频散关系曲面，图3e、f为model 1的近似qP波和qSV波频散关系曲面，图3c、d为model 2的精确qP波和qSV波频散关系曲面，图3g、h为model 2的近似qP波和qSV波频散关系曲面.两组模型中的近似qP波和qSV波频散关系曲面(图3e—h)与精确频散关系曲面(图3a—d)方位特征较为相似，其中qP波的精度明显高于qSV波.

图3 TTI介质弹性波频散关系曲面(a) model 1 qP波精确值; (b) model 1 qSV波精确值; (c) model 2 qP波精确值; (d) model 2 qSV波精确值; (e) model 1 qP波近似值; (f) model 1 qSV波近似值; (g) model 2 qP波近似值; (h) model 2 qSV波近似值.Fig.3 Surfaces of elastic wave dispersion relation in TTI media(a) Exact value of qP wave in model 1； (b) Exact value of qSV wave in model 1； (c) Exact value of qP wave in model 2； (d) Exact value of qSV wave in model 2； (e) Approximate value of qP wave in model 1； (f) Approximate value of qSV wave in model 1； (g) Approximate value of qP wave in model 2； (h) Approximate value of qSV wave in model 2.

为考察近似值与精确值的误差，这里讨论频散曲面在XOZ平面内的频散关系曲线的相对误差分布情况.在观测坐标系XOZ平面内，基于近似配方法的频散关系方程(6)可表示为：

(15)

定义相对误差(RE)为：

(16)

根据表1中两组模型参数可得到二维TTI介质弹性波频散关系曲线，如图4所示，图中短实线表示TTI介质对称轴方向，qP波和qSV波近似频散关系曲线与精确值在对称轴和垂直对称轴方向误差为零，且误差关于上述两个传播方向呈现对称分布的形式.对比图4a、b可知，qP波近似频散关系精度整体上高于qSV波，原因在于qSV波近似频散关系的精度在较大程度上受组合参数σ的控制，该参数中纵横波速度比起到放大ε和δ之间差值的作用，从而导致了误差的增大.图5较为直观地呈现了qP波和qSV波近似频散关系曲线相对误差在θ∈[0,180°]范围内的分布情况，本例中沿对称轴(空心箭头处)和垂直对称轴(实心箭头处)的两个传播方向对应的相对误差为零，并且相对误差在这两个传播方向附近一定角度范围内很小，这同时也契合差异项ΔD的变化趋势(图2).

图4 TTI介质弹性波频散关系曲线(a) model 1 qP波； (b) model 1 qSV波; (c) model 2 qP波； (d) model 2 qSV波; 灰色实线表示精确频散关系曲线，黑色虚线表示近似频散关系曲线.Fig.4 Dispersion curves of elastic waves in TTI media(a) qP wave in model 1； (b) qSV wave in model 1； (c) qP wave in model 2； (d) qSV wave in model 2. Gray solid curve represents exact values, and black represents approximate values.

图5 qP波和qSV波近似频散关系曲线的相对误差(a) model 1相对误差； (b) model 2 相对误差.其中灰色实线表示qP波频散关系相对误差EP，黑色虚线表示qSV波频散关系相对误差ESV.Fig.5 Relative error of approximate dispersion relation curve for qP and qSV wave(a) Model 1； (b) Model 2. Gray solid curve represents relative error of dispersion relationship for qP wave (EP), and the black dashed curve represents relative error of dispersion relationship for qSV wave (ESV).

2.2 波场模拟

为验证文中导出的TTI介质qP波和qSV波近似波动方程(10)—(13)的有效性，本文利用有限差分法对其进行均匀、层状和复杂介质波场数值模拟.首先，设置均匀介质模型大小为401×401，空间间隔为Δx=Δz=10 m，时间采样间隔为Δt=1 ms.模型参数为：vP0=3.00 km·s-1，vS0=1.73 km·s-1，η=0.143，δ=0.2，σ=0.601，θ0=π/4；震源子波是主频为25 Hz的雷克子波，位于模型中心点.高阶空间有限差分近似有利于压制数值频散，但会增加运算成本，合理选择差分近似的阶数可以在适当增加运算时间的情况下得到更好的数值模拟结果.考虑到泊松方程的求解效率问题，本文对方程(10)—(13)取空间十阶、时间二阶精度的差分近似进行波场模拟.图6展示了TTI介质弹性波以及解耦的qP波和qSV波波场快照.对比图6a、c可知纯qP波的运动学特征基本与弹性波中qP波的特征相似，即使在大偏移距处，两者的波场特征也很相似；对比图6b、c可知，相比于弹性波qSV波场，纯qSV波的波前面三角区不明显，但其运动学特征整体保持较好的相似度.

图6 均匀TTI介质波场快照(a) 纯qP波； (b) 纯qSV波; (c) 弹性波.Fig.6 Snapshots of homogeneous TTI media(a) Pure qP wave; (b) Pure qSV wave; (c) Elastic wave.

其次，对层状各向异性介质模型进行数值模拟.模型大小为401×401，空间步长为Δx=Δz=10 m，时间采样间隔为Δt=1 ms，模型参数如表2所示.注意到表2中layer 1和layer 2为η>0的TTI介质模型，layer 3是η<0的TTI介质模型，这样设置是为了测试纯qP波和纯qSV波波动方程的稳定性.震源子波与均匀介质模型中相同，位于(2 km, 0.1 km)处，检波点设置于深度为0.1 km的水平面上，如图7a所示；图7b—d为0.8 s时刻的TTI介质地震波波场快照.通过对比可知，纯qP波(图7b)和纯qSV波(图7c)入射到界面时只产生反射、透射的同类波，不产生转换波.相对常规弹性波场(图7d)，前两者有着更加清晰的波场，利于研究同类型波的传播特征.此外，数值模拟结果表明，纯qP波方程在η<0和各向异性倾角变化较大情况下，仍然可以保持较好的稳定性.

图7 层状TTI介质模型(a)、纯qP波波场快照(b)、纯qSV波波场快照(c)以及弹性波波场快照(d)Fig.7 Layered TTI media model (a), snapshots of pure qP wave (b), pure qSV wave (c) and elastic wave (d)

表2 层状各向异性介质模型Table 2 Layered anisotropic medium model

最后，本文对BP盐丘模型进行数值模拟；模型大小为601×601，空间间隔和时间采样间隔与上述均匀介质模型一致，模型参数如图8所示.震源子波与均匀介质模型中相同，位于(3 km, 3 km)处.数值模拟结果表明，纯qP波波场中只可见透射和反射qP波的信息，其波前面与弹性波qP波的波前面位置几乎一致；相应地，纯qSV波中只可见透射和反射qSV波的信息；从相速度误差分析中可见，qSV波的误差大于qP波的误差，因此纯qSV波的波前面与弹性波SV波波前面存在一定的差异，但总体上在可接受的范围内(图9).

图8 BP盐丘模型(a) vP0; (b) vS0; (c) η; (d) δ; (e) θ0.Fig.8 Salt dome model of BP

图9 BP盐丘模型波场快照(a) 纯qP波; (b) 纯qSV波; (c) 弹性波.Fig.9 Snapshots of BP salt dome model(a) Pure qP wave; (b) Pure qSV wave; (c) Elastic wave.

3 结论和认识

本文在TTI介质弹性波精确频散关系的基础上，利用近似配方法推导了qP波和qSV波近似频散关系，通过傅里叶逆变换将其从频率-波数域变换到时-空域，进而导出了纯qP波和纯qSV波近似解耦波动方程.

(1) 理论分析和数值计算表明，近似和精确频散关系的差异项ΔD与|η|呈正相关，在椭圆各向异性介质中差异项为零，即近似式与精确式等效.

(2) 近似频散关系曲线相对误差整体较小，且在沿对称轴和垂直对称轴的传播方向上相对误差为零.

(3) 参数σ中纵横波速度比使得qSV波近似频散关系曲线相对误差普遍较qP波大.

(4) 数值结果表明，基于近似配方法解耦的纯qP波和纯qSV波波动方程在η<0和对称轴倾角变化较大的介质模型中依然可以保持稳定；相对于弹性波波场，两者均不产生转换波，波场更加清晰.

本文推导的TTI介质解耦波动方程可以被进一步推广到任意空间取向的TI介质中，此时需要从三维空间的角度讨论解耦纵横波的传播特征.本文的解耦波动方程可以用于各向异性介质正演模拟和逆时偏移算法中，但因用常规有限差分方法求解泊松方程的效率较低，可以考虑采用伪谱法、低秩分解法、有限差分法和伪谱法的混合法等效率较高的方法求解纯qP波和纯qSV波方程，以提高正演模拟和偏移成像的效率.

附录A

在qP波和qSV波近似频散关系(6)两边同时乘以波场PP(ω,kx,ky,kz)和PSV(ω,kx,ky,kz)，再进行傅里叶逆变换(ω→i∂/∂t,kx→i∂/∂x,ky→i∂/∂y,kz→i∂/∂z)，化简合并后得到三维TTI介质纯qP波和纯qSV波的近似波动方程：

(A1)

(A2)

(A3)

(A4)