基于“方程解耦”算法的各向异性介质纵横波分离方法

孔丽云， 胡志方*， 王一博

1 中国地质调查局油气资源调查中心， 北京 100083 2 中国科学院地质与地球物理研究所， 北京 100029

0 引言

地球介质中广泛存在着各向异性，其纵横波耦合的特性，使得基于各向同性假设的正演模拟方法不足以准确地描述实际地球介质中地震波的传播规律.在以往的研究中，考虑各向异性的弹性波动方程正演模拟取得了一定的进展.该方法是用三分量矢量场来描述弹性波波场，其正演结果能够获得比标量波场更丰富的波场信息.但其存在一定的缺陷(何兵寿和张会星，2006)：①各向异性参数复杂，物理意义不明确，在实际生产中难以获得；②多分量数据占用资源大，计算效率低；③由于多分量地震勘探成本高，目前野外采集仍然以纵波地震资料为主，多分量弹性波理论的应用受到资料不足的限制.因此，需要对各向异性介质的纵波分量进行提取，以适应目前阶段地震数据采集和处理的需求.目前主要有两种思路，一种是对介质的弹性波场进行纵横波波场的分离，另外一种就是通过频散关系式对qP波近似方程进行推导.

纵横波分离方面，常用的方法是Helmholtz分解(Yan and Sava，2008).其基本原理是基于纵横波偏振方向的不同(纵波偏振方向平行于传播方向，横波偏振方向垂直于传播方向)，利用波场的散度和旋度来分别对纵波和横波进行求取.该算法在各向同性情况下适用，但在各向异性介质中，纵、横波的传播方向发生了偏离，与极化方向不再是平行和垂直的关系，所以Helmholtz分解对于各向异性介质不适用.Yan和Sava(2009)提出了基于空间域波场分离算子的各向异性介质纵横波分离算法，但该算法的计算效率太低，在3D情况下几乎无法实现.为了提高计算效率，Yan 和Sava(2011)提出利用不同权重的均匀反射模型对非均质模型进行近似，这种方法实质上是在计算效率和精度之间进行了一个权衡，但反射模型的选择和计算代价都依赖于模型本身，难以应用于实际情况.Cheng 和Fomel(2014)转换思路，将波场传播方向和极化方向之间的投影关系重新改写为空间波数域的积分算子，通过低阶近似的方法对其进行数值求解.计算效率得到了部分提高，但仍然会随着模型的尺寸和复杂度而降低.Zhou 和Wang(2017)通过将波场传播方向局部旋转到qP波的极化方向，发展了VTI介质波场分离的新算子，有效提高了计算效率.另外，由于在分离过程中采用散度场代表纵波，旋度场代表横波，这就意味着，纵波仅代表整体位移的输入输出，而横波则代表旋转变换的角速度，其物理意义与原波场并不相同，且其振幅和相位信息都发生了改变.在各向同性介质中，Duan和Sava(2015)提出，可以通过在时间上对震源子波积分后乘以纵、横波相速度的外推波场对分离后的波场进行矫正.在各向异性介质中，由于波形改变变得更加复杂，没有跟各向同性相似的矫正方法.Sun等(2001)通过Hilbert变换矫正相位的改变；Zhang和McMechan(2010)对VTI介质下波数域二维和三维矢量波波场分离方法进行了研究；Sun等(2011)用纵横波速度恢复 P波和S波的相对振幅比，通过地表P波速度和频率恢复P和S波的真正振幅值.该方法虽然能恢复振幅与相位，但存储空间耗费量巨大.

qP波近似方程方面，Alkhalifah从VTI介质的精确qP-qSV波频散关系出发，假设垂向剪切波速度Vs0=0，推导了VTI介质四阶拟声波方程.为降低四阶偏微分方程计算的复杂性(Alkhalifah，2000)，周熙焱等(2018)、Zhou等(2019)和Du等(2008)从声学近似的频散关系出发，通过引入不同的辅助波场函数，分别推导了两种VTI介质二阶耦合拟声波方程.为结合交错网格技术实现精确的有限差分数值模拟，Hestholm推导了一阶偏微分形式的VTI介质拟声波方程(Hestholm，2009).Duveneck等(2008)从声学近似的胡克定律出发，也推导了二阶形式和一阶应力-速度形式的VTI介质拟声波方程.与Alkhalifah方程一样，几种拟声波方程均存在低速度、低振幅的qSV人为干扰波问题(Grechka et al.，2004)，而这种干扰波的存在会影响正演模拟和偏移成像的最终结果.同时，由于无法考虑到转换波的信息，利用qP波近似方程得到的各向异性介质纵波数据并不完整.

为得到振幅和相位信息相对准确且波场信息相对完整的各向异性介质纵波数据，本文采用波场分离的算法进行各向异性介质的纵横波分离，其基本思想是源于Zhang等(2007)提出的“方程解耦”的方法.该方法是直接将弹性刚度系数分量分解为纵波和横波相关的系数，进而达到纵横波方程分离的效果，因此不需要额外的计算.在各向同性介质中，Zhang等(2007)得到精确的纵横波分离的一阶速度-应力方程；Zhou等(2016)将这种方法进一步拓展到三维各向同性模型的情况.在各向异性介质中，纵横波是耦合在一起的，这就使得其弹性波场的纵横波分离需要在一定的近似条件下才能完成.Bloot等(2013)从Thomsen(1986)弱各向异性参数出发，通过假设δ=0(近垂直P波传播中关键的各向异性参数)，从弹性刚度系数矩阵中的耦合项C13中分离出纵波模量.

本文借助于Bloot等(2013)的近似方法，将弹性刚度系数分量中的耦合项C13近似为纵波模量和横波模量的线性函数，推导得到了各向异性介质的纵横波分离的一阶速度-应力方程，从而能够实现各向异性介质的纵横波波场分离.

1 传统的波场分离算法

1.1 Helmholtz分解

对于纵横波波场分离技术，本文所使用的是基于旋度场和散度场的交错网格高阶有限差分的纵横波波场分离方法.首先，根据场的理论，任何一个矢量场都可定义为一个标量场的散度和一个矢量场的旋度相加的结果，因此，质点速度v=(vx,vz)可以写为如下表达形式：

(1)

其中，ψp和ψs分别代表速度场的标量位和矢量位，即纵波和横波的位移位.其次，对于上述速度场分别求其散度和旋度，可以得到速度矢量的无旋场与无散场，也即纵波波场vp与横波波场vs.其表达式分别为：

(2)

(3)

在二维情况下，上述方程可以简化为：

(4)

(5)

从式(4)和式(5)中可以看出，Helmholtz分解得到的纵波是标量，在二维情况下，横波也只有平行于y轴的一个分量，可视为标量.

1.2 声波近似方程

首先，对VTI介质的虎克定律进行经典的声波近似假设(Alkhalifah，2000)，令垂向剪切波速度vs0=0，得到用Thomsen参数表示的VTI介质声波近似胡克定律：

(6)

式中，σij(i,j=x,y,z)为应力分量，eij(i,j=x,y,z)为应变分量，ρ为密度，Vp0为VTI介质对称轴方向上的纵波速度，ε和δ为Thomsen各向异性参数.

表达式(6)两端对时间求导，并结合运动方程

ρ∂vi/∂t=∂σij/∂xii,j=x,y,z

(7)

可以得到一阶应力-速度形式的VTI介质qP波方程

(8)

式中，vx,vy,vz分别为沿x,y,z方向的速度分量，p=σ11=σ22和q=σ33为水平应力分量和垂直应力分量.

在二维情况下，方程(8)简化为：

(9)

2 纵横波方程解耦算法

2.1 基本思想

声波和弹性波的区别在于其传播介质的不同，前者是在声学介质中传播，只能传播纵波；而后者则是在弹性空间中传播，除了纵波还会产生横波.而传播介质的不同主要体现在刚度系数矩阵上，因此，如果能将弹性刚度系数矩阵分离为分别只与纵波和横波相关的两部分，那么就可以将整个弹性波方程分解为纵波方程和横波方程，从而达到纵横波分离的目的.基于这一思想，给出适用于交错网格高阶有限差分算法的各向异性纵横波分离方程的具体推导过程.

2.1.1 二阶波动方程

三维弹性波动方程通常表示为：

(10)

式中，ρ为介质密度；u=(ux,uy,uz)T为位移矢量；C为刚度系数矩阵；G为偏微分算子矩阵.

根据Zhang等(2007)的思想，假设完整的弹性系数矩阵C可以分为与纵波和横波分别相关的弹性系数矩阵Cp和Cs两部分，那么弹性波动方程就可以拆分为：

(11)

其中，up和us分别为纵波波场和横波波场，全波场u=up+us.

2.1.2 一阶速度-应力方程

胡克定律表达式为：

T=CE=CGTu,(12)

其中，T=(σxx,σyy,σzz,σyz,σxz,σxy)T为应力向量，E=(exx,eyy,ezz,eyz,exz,exy)T为应变向量.将上述表达式代入方程(11)中，可以得到纵横波分离的一阶速度-应力方程表达式：

(13)

从上述推导过程中可以看出，如何将完整的弹性系数矩阵C分解为分别只与纵波和横波相关的Cp和Cs两部分，是该方法的核心问题.对于各向同性介质，由于纵波和横波能够完全解耦，所以能够得到完全精确的纵波和横波方程(Zhang et al.，2007).但是在各向异性介质中，由于纵横波无法解耦，我们需要对相关系数进行一定的假设，才能将弹性系数矩阵C进行分解，下面给出各向异性VTI介质弹性系数矩阵分解的过程.

2.2 弹性系数矩阵的分解

2.2.1C13的线性近似表达

对于VTI介质，刚度系数不为零分量用Thomsen参数可以表达为：

其中，Vp0、Vs0分别为qP波和qSV波沿VTI介质对称轴方向的相速度；ε、δ和γ是表示VTI介质各向异性强度的三个无量纲因子：ε是度量qP波各向异性强度参数，ε越大，介质的纵波各向异性越大，ε=0，纵波无各向异性；δ是连接Vp0和Vp90之间的一种过渡性参数；γ可以看成是度量qS波各向异性强度或横波分裂强度的参数，γ越大，介质的横波各向异性越大，γ=0时，横波无各向异性.

(15)

从(15)式可以看出，系数A、B的求解是C13和C23分解的关键.由于C13和C23的表达式相同，以C13为例，给出其系数A和B的近似求解过程.

2.2.2 系数的求取

对式(14)中C13的表达式进行重新整理得：

(16)

为求取横波分量系数B，首先假设δ=0，得到

进一步，假设Vp0=0，可得：

B=-2.

(17)

将(17)式代入C13近似表达式(15)中，就可以计算得到系数A为：

(18)

2.2.3Cp和Cs近似表达式

将系数表达式(17)和(18)代入C13的近似表达式(15)中，结合式(14)中其他分量的表达式，可以得到VTI介质中的纵波弹性模量Cp和横波模量Cs的分量表达式分别为：

(19)

(20)

将表达式(19)、(20)代入方程(13)中，可以得到VTI介质在三维情况下的纵横波分离的一阶速度-应力方程.

2.3 二维分解公式

在二维情况下，VTI介质中的一阶速度-应力方程表达式为：

(21)

其中，

D11m=(C13C13m-C33C11m)M-1,

D12m=(C13C11m-C11C13m)M-1,

D13m=C55m/C55,

D21m=(C13C33m-C33C13m)M-1(m=p,s)

D22m=(C13C13m-C11C33m)M-1,

D23m=C55m/C55,

(22)

3 理论模型测试

为说明“方程解耦”算法的精确性和可靠性，基于x和z两个方向的地震波场，利用二维公式(22)对三层水平层状介质模型进行测试，并将波场分离结果与同样考虑x、z两个方向地震波场的Helmholtz分解(式(4)、(5))和拟声波近似(式(9))等传统算法进行比较.

3.1 理论模型及测试

设计三层水平层状介质模型如图1所示，中间层为VTI介质，其顶层和底层为各向同性弹性介质模型，介质参数如表1 所示.模型大小为704 m×704 m×528 m，空间网格间距为8 m；模拟的物理时间为0.32 s，时间步长为0.4 ms.震源采用爆炸震源，主频25 Hz，波源函数表达式为：

(23)

其中，h(t)为雷克子波，(x0,z0)为震源中心点，α决定了震源子波的宽度，f0为主频，t0为延迟时间.

表1 三层模型物理参数表Table 1 Physical parameters for three-layer model

图1 三层模型示意图Fig.1 Diagram of three-layer model

图2 爆炸震源激发的纵波路径示意图(去除直达波)，包含入射纵波-反射纵波(P-P，r1，橙色)，入射纵波-透射转换横波-反射横波-透射转换纵波(P-S-S-P，r2，绿色)，入射纵波-透射转换横波-反射转换纵波-透射纵波(P-S-P-P，r3，红色)，入射纵波-透射纵波-反射转换横波-透射转换纵波(P-P-S-P，r3，红色)，入射纵波-透射纵波-反射纵波-透射纵波(P-P-P-P，r4，蓝色)Fig.2 Diagram of P-wave path excited by explosion source (removing direct wave), including incident P-refeccted P (P-P, r1, orange), incident P-transmitted converted S-reflected S-transmitted converted P(P-S-S-P, r2, green), incident P-transmitted converted S-reflected converted P-transmitted P (P-S-P-P, r3, red), incident P-transmitted P-reflected converted S-transmitted converted P (P-P-S-P, r3, red), incident P-transmitted P-reflected P-transmitted P (P-P-P-P, r4, blue) (P is longitudinal wave and S is transverse wave)

3.2 波场的完整性比较

图2和图3给出了爆炸震源激发后，在接收点接收到的来自I、II两个地层界面的纵横波射线路径示意图.从图中可以看出，纵波和横波都存在四条反射路径(其中红色路径是两条路径的重合，由于到达时间相同，在这里用一条路径代替).所以，在分离后的地震剖面上，应该存在四条同相轴.图4和图5分别为基于方程解耦算法在x方向和z方向的纵横波结果，为了说明该算法的有效性，将其与Helmholtz分解(图6)和各向异性拟声波方程(图7)的纵横波分离结果进行对比.由于Helmholtz分解算法在二维情况下得到的纵波和横波都是标量，图6中只显示了纵波和横波的综合分离结果，没有方向性.另外，由于拟声波近似算法的结果中只有纵波信息，所以图9中只给出了纵波分别在x方向和z方向的分离结果.

图3 爆炸震源激发的横波路径示意图(去除直达波)，包含入射纵波-反射转换横波(P-S，r′1，橙色)，入射纵波-透射转换横波-反射横波-透射横波(P-S-S-S，r′2，绿色)，入射纵波-透射转换横波-反射转换纵波-透射转换横波(P-S-P-S，r′3，红色)，入射纵波-透射纵波-反射转换横波-透射横波(P-P-S-S，r′3，红色)，入射纵波-透射纵波-反射纵波-透射转换横波(P-P-P-S，r′4，蓝色)Fig.3 Diagram of S-wave path excited by explosion source (removing direct wave), including incident P-reflected converted S (P-S, r′1, orange), incident P-transmitted converted S-reflected S-transmitted S (P-S-S-S, r′2, green), incident P-transmitted converted S-reflected converted P-transmitted converted S (P-S-P-S, r′3, red), incident P-transmitted P-reflected converted S-transmitted S (P-P-S-S, r′3, red), incident P-transmitted P-reflected P-transmitted convered S (P-P-P-S, r′4, blue)(P is longitudinal wave and S is transverse wave)

图4 基于方程解耦算法的数值模拟结果(x分量)(a) 完整的波场； (b) 纵波波场； (c) 横波波场.Fig.4 Numerical simulation results based on equation decoupling algorithm (x component)(a) Complete wave field; (b) P-wave; (c) S-wave.

图5 基于方程解耦算法的数值模拟结果(z分量)(a) 完整的波场； (b) 纵波波场； (c) 横波波场.Fig.5 Numerical simulation results based on equation decoupling algorithm (z component)(a) Complete wave field; (b) P-wave; (c) S-wave.

图6 基于Helmholtz分解的数值模拟结果(a) 纵波波场； (b) 横波波场.Fig.6 Numerical simulation results based on Helmholtz decomposition(a) P-wave; (b) S-wave.

图7 基于拟声波近似的纵波模拟结果(a) x分量； (b) z分量.Fig.7 Numerical simulation results of longitudinal wave based on acoustic approximation(a) x component; (b) z component.

通过对比可以看出，基于方程解耦算法(图4和图5)和基于Helmholtz分解算法得到的波场(图6)，都存在四条反射波同相轴，而基于拟声波近似方法得到的波场(图7)中，只存在两条反射波同相轴.因此，基于方程解耦算法和Helmholtz分解算法得到的纵波数据，都能够得到完整的反射波场.

3.3 振幅和相位的精确性比较

为了说明方程解耦算法在地震波振幅和相位方面的精确性，分别取方程解耦算法、Helmholtz分解、拟声波近似等三种方法数值模拟结果中的第100道地震数据进行对比.图8—11分别是x方向纵波分量、z方向纵波分量、x方向横波分量、z方向横波分量与相对应全波场信息的比较.为方便对比，将Helmholtz分解的标量结果同样显示在x方向和z方向.

对比纵波(图8和图9)可以看出，三种方法在地震波走时上与全波场都保持一致.方程解耦算法(红色实线)和拟声波近似方法(紫色实线)的振幅和相位信息与全波场(黑色虚线)基本保持一致，但是Helmholtz分解算法(蓝色实线)得到的纵波信息，无论在振幅还是相位方面，都与全波场(黑色虚线)存在很大误差.在横波的对比结果中(图10和图11)同样可以看出，相对于Helmholtz分解算法(蓝色实线)，方程解耦算法在横波振幅和相位上与原始波场保持高度一致性.因此，基于方程解耦算法得到的纵波数据，能够得到准确的振幅和相位信息.

图8 不同方法(蓝实线：Helmholtz分解；紫实线：拟声波方程；红实线：方程解耦)得到的纵波信息与全波场信息(黑色虚线)的对比(x分量)(a) 0.1～0.38 s； (b) 0.38～0.64 s.Fig.8 Comparison of full wavefield (black dotted line) and P-wave obtained by different methods (blue solid line: Helmholtz decomposition; purple solid line: acoustic approximation; red solid line: equations decoupling)(x component)

图9 不同方法(蓝实线：Helmholtz分解；紫实线：拟声波方程；红实线：介质分离)得到的纵波信息与全波场信息(黑色虚线)的对比(z分量)(a) 0.1～0.38 s； (b) 0.38～0.64 s.Fig.9 Comparison of full wavefield (black dotted line) and P-wave obtained by different methods (blue solid line: Helmholtz decomposition; purple solid line: acoustic approximation; red solid line: equations decoupling) (z component)

图10 不同方法(蓝实线：Helmholtz分解；红实线：介质分离) 得到的横波信息与全波场信息(黑色虚线)的对比(x分量)(a) 0.1～0.38 s； (b) 0.38～0.64 s.Fig.10 Comparison of full wavefield (black dotted line) and S-wave obtained by different methods (blue solid line: Helmholtz decomposition; red solid line: equations decoupling) (x component)

图11 不同方法(蓝实线：Helmholtz分解；红实线：介质分离) 得到的横波信息与全波场信息 (黑色虚线)的对比(z分量)(a) 0.1～0.38 s； (b) 0.38～0.64 s.Fig.11 Comparison of full wavefield (black dotted line) and S-wave obtained by different methods (blue solid line: Helmholtz decomposition; red solid line: equations decoupling)(z component)

4 结论

基于方程解耦的思想，通过将弹性刚度系数分量中的耦合项C13近似表示为纵波模量和横波模量的线性函数，推导得到了各向异性介质的纵横波分离的一阶速度-应力方程，从而能够实现各向异性介质的纵横波波场分离.与拟声波方程算法相比，方程解耦算法得到的纵波波场信息较为完整；与Helmhotlz分解算法相比，方程解耦算法得到的纵波波场的振幅和相位信息都相对准确.因此，对于各向异性介质而言，方程解耦算法既能保证分离后波场的完整性，又能保持波场振幅、相位信息的准确性，为各向异性非常规储层地震数据处理、解释的可靠性提供了重要的数据保障.