结合泊松算法和波场分解成像条件的各向异性优化纯声波方程逆时偏移

徐世刚 包乾宗 任志明 刘 洋

(①长安大学地质工程与测绘学院地球物理系，陕西西安 710054； ②中国石油大学(北京)油气资源与探测国家重点实验室，北京 102249)

0 引言

地下介质普遍呈现各向异性，横向各向同性(TI)介质是地震勘探学界研究最广的一类各向异性模型，根据其对称轴方向可以分为VTI、HTI和TTI模型。逆时偏移技术受倾角限制小，能够适应复杂速度模型，获得高精度成像剖面，因此各向异性逆时偏移研究具有重要意义。

相比于各向异性弹性波方程，各向异性声波方程形式简单、计算量较小、波场特征明确。目前在各向异性波场模拟和偏移中主要采用伪声波方程，该类方程的一种推导思路是将耦合的各向异性P-SV波频散关系中的横波速度设置为0简化而来[1]。在此基础上，发展了多种形式的伪声波方程[2-8]。另一种伪声波方程的推导思路是基于声学近似的胡克定律[9-10]。但是当给定的各向异性参数不满足伪声波假设时，现有伪声波方程在进行P波模拟时，容易出现伪横波干扰及数值不稳定。为了彻底解决伪声波存在的问题，Liu等[11]解耦了VTI介质P-SV波频散关系，获得单独的P波和SV波方程。相比于显式表达的伪声波方程，各向异性纯声波方程包含复杂的拟微分算子，常规差分法难以直接求解。迄今为止，已经发展了多种数值算法求解纯声波方程，如伪谱法/伪解析法[12-13]、低秩分解法[14]、迭代求解法[15]、混合求解法[16]和快速泊松算法[17-18]等。目前关于纯声波方程的研究主要集中于平衡计算精度和效率方面，以及在偏移成像中的进一步应用等[19-23]。

传统逆时偏移技术对震源波场和检波点波场做互相关，结果中容易出现较强的低频噪声。Zhang等[24]采用拉普拉斯滤波压制低频噪声； Yoon等[25]、康玮等[26]利用Poynting矢量切除大角度成像结果，从而压制低频噪声； Liu等[27]提出波场分解构建逆时偏移成像条件，选择不同波场做互相关，能够较好地抑制成像噪声。在Liu等[27]的研究基础上，基于Hilbert变换构建复数波场的策略被用于提高计算效率[28-32]。波场分解成像条件已应用于各向同性声波和弹性波逆时偏移，考虑到该方法的优势，本文将其推广到各向异性介质逆时偏移。

作为对比，本文首先回顾了两种常用的TTI伪声波方程，分析了其存在的问题； 然后利用最小二乘方法求取了TTI优化纯声波频散关系，并采用泊松算法和有限差分相结合的数值方案进行求解； 最后结合介质各向异性特征对波场分解成像条件进行改进，将其推广至纯声波逆时偏移。相关理论和算例表明，将优化纯声波方程和波场分解成像条件联合应用于各向异性纯声波逆时偏移，能有效压制伪横波及低频噪声，获得高质量成像结果。

1 理论方法

1.1 基于限制性横波速度的TTI伪声波方程

通过将原始P-SV波频散关系中的横波速度设置为0，Alkhalifah[1]推导了各向异性伪声波方程，该方程能够有效简化计算，其多种改进版本被广泛应用于各向异性波场模拟和偏移成像[2-8]。但直接将垂向横波速度设置为0并不能完全满足各向异性介质中的物理学规律，当各向异性参数不满足近似条件时(即椭圆近似条件)，伪声波方程模拟容易出现SV波干扰及不稳定。为此，学者们提出了不同的解决方案，其中一种常用的基于限制性横波速度的TTI伪声波方程形式[4]为

(1)

(2)

其中θ为对称轴倾角。特别地，通过引入限制性横波速度能够有效提高稳定性。当vSz=0时，式(1)退化为传统伪声波方程，形式为

(3)

采用有限差分即可求解式(1)和式(3)。图1为两个方程计算的双层TTI模型(表1)的波场切片，其中上、下层介质厚度均为1500m。

表1 双层TTI模型参数

由图1可以看出：当介质具有椭圆各向异性特征(ε=δ)时，两种方程均可产生高精度P波响应，如图1左所示； 当各向异性参数不满足椭圆假设时(ε≠δ)，两种方程会产生较强的伪横波污染，如图1中和图1a右所示； 当ε<δ时，横波速度为0的伪声波方程(式(3))模拟出现不稳定，如图1b右所示。

1.2 TTI优化纯声波方程及其结合泊松算法和有限 差分的解法

多种改进版本的各向异性伪声波方程仅能在一定程度上压制伪横波干扰，提高稳定性，但不能从根本上解决问题。众多研究表明，发展纯声波方程能够彻底消除伪横波污染及不稳定。相比于显式表达的伪声波方程，纯声波方程包含复杂的偏微分算子，常规差分法难以直接求解。考虑到计算精度和效率，本文给出一种改进的、结合泊松算法和有限差分的各向异性优化纯声波方程求解方法[17-18]。

解耦的TTI纯声波频散关系[11]为

(4)

(5)

将式(5)代入式(4)，可得到目前最常用的各向异性纯声波频散关系。考虑到一阶泰勒展开的精度较低，高阶泰勒展开的计算量较大，为了平衡纯声波模拟时的精度和效率，本文借鉴优化的思想[17-18]获得一种高精度的纯声波频散关系。对式(4)中的根号项做如下近似[17-18]

图1 双层TTI模型式(1)(a)和式(3)方程模拟的波场切片(b)左：模型1； 中：模型2； 右：模型3

(6)

式(6)相对于待求近似系数(f1，f2，f3)为线性方程，因此可以通过极小化两端的误差获得近似系数值。构建的目标函数为

(7)

通过最小二乘方法求解上式，即可获得优化系数(f1，f2，f3)。

将式(6)代入式(4)，整理可得优化的纯声波频散关系

(8)

式中

(9)

特别地，当w1=1+2ε、w2=1、w3=-2(ε-δ)时，式(8)退化为常用的一阶泰勒展开纯声波频散式。式(8)在时间—空间域可写为[17]

(10)

引入辅助波场Q，将式(10)分解为[15，17]

(11)

(12)

式(11)为显式方程，采用常规差分即可求解。式(12)为隐式方程，其计算过程涉及线性方程组求解。Chu等[15]采用线性迭代求解隐式方程，但计算量较大。为了提高效率，本文给出一种基于泊松算法的实现策略[17]。

将式(11)和式(12)进一步整理为

(13)

∇2Q=P

(14)

式(13)和式(14)组成了新的TTI优化纯声波方程，其中式(13)可采用常规差分法求解。式(14)为泊松方程，基于前人研究成果，可通过泊松算法实现，计算步骤[17]如下。

(1)采用五点差分离散式(14)

Qi-1,j+Qi+1,j+Qi,j-1+Qi,j+1-4Qi,j=h2Pi,j

(15)

式中：h为网格间距；i、j分别为x、z方向网格点序号。

(16)

式中：Lx和Lz分别为x、z方向网格点数；n、l分别为波数域两个方向的网格点序号。

(3)计算辅助波场Q的傅里叶响应

(17)

(4)求解辅助波场Q

(18)

通过显式差分离散式(13)，通过式(15)～式(18)可实现式(14)的求解。这样能够保证纯声波方程式(13)和式(14)的高效求解。

图2对比了一阶泰勒展开方程(式(5))和优化展开方程(式(6))对根号部分的近似精度。由图可见，与精确值(图2a)相比，一阶泰勒展开存在较大误差(图2b)，而最小二乘优化展开能够获得更高的近似精度(图2c)。图3为优化纯声波方程模拟的双层TTI模型波场切片，与图1伪声波方程的计算结果相比，纯声波方程在不同参数条件下均能产生精确有效的P波响应，同时能够完全消除伪横波和数值不稳定。

1.3 改进波场分解成像条件

在逆时偏移成像中，互相关成像条件具有实现简单、计算高效等优点，得到广泛应用。常规互相关成像条件[34]为

(19)

式中：x为空间坐标；I(x)为成像结果；Ps(x,t)为震源波场；Pr(x,t)为检波点波场；tmax为记录长度。

常规互相关成像条件易引起较强低频噪声，原因是互相关过程中沿所有方向传播的震源波场和检波点波场均参与成像，同向传播的波场之间互相关叠加，易形成较强的低频噪声。为了压制传统互相关成像过程中的低频噪声，除了对低频噪声进行滤波处理外[24]，另一种有效策略是对参与成像的全波场沿不同方向进行分解，选取传播方向不同的分量波场参与成像，能够有效减弱低频噪声干扰[27-32]。

以二维情况为例，对式(19)中的震源波场Ps(x,t)和检波点波场Pr(x,t)进行分解

(20)

(21)

式中上标“↑”、“↓”、“←”、“→”分别表示上行、下行、左行和右行波场分量。

由式(19)～式(21)可知，震源波场分量与检波点波场分量之间两两互相关可得16种成像剖面。其中传播方向相反的波场分量互相关能够产生有效的成像结果，而同向传播的波场分量之间互相关主要产生低频噪声[27-32]。在前人关于波场分解成像条件的研究基础上，本文采用一种基于Hilbert变换构建复数波场的分解算法。该分解方法能够进行全波场的显式分离，目前被广泛用于各向同性介质的

图2 均匀各向异性模型的精确微分算子(a)及一阶泰勒展开误差(b)与最小二乘优化展开误差(c)的对比模型参数为ε=0.2，δ=0.1，θ=60°

图3 双层TTI模型优化纯声波方程模拟的波场切片(a)模型1； (b)模型2； (c)模型3

逆时偏移[27-32]。本文将其推广到各向异性介质的逆时偏移，实现步骤如下。

(1)基于式(13)和式(14)，震源波场的复数外推形式可以表示为

(22)

(2)基于步骤(1)构建的复数波场，以震源波场为例，上行波场可以表示为

(23)

(24)

图4 双层TTI模型4的波场切片(a)及分离的上(b)、下(c)、左(d)、右行波场(e)

2 模型算例

2.1 地堑模型

各向异性地堑模型的尺寸为4000m×3000m，相关参数如图5所示。网格间距为10m，时间步长为1ms，震源选用主频为28Hz的Ricker子波。

考虑到横波速度为0的伪声波方程存在稳定性问题，图6为基于限制性横波速度的伪声波和优化纯声波方程模拟的波场切片，可以看出，相比于伪声波方程的模拟结果(图6a)，纯声波方程能够有效压

图5 各向异性地堑模型

制伪横波干扰，产生高精度P波响应(图6b)。

图7为图6b的纯声波波场沿不同方向分解结果，可以发现，上、下、左、右行波得到了有效分离。传统互相关成像条件和波场分解成像条件的偏移结果如图8所示。由图8a可知，基于传统互相关条件的成像结果容易产生较强低频噪声干扰。对比部分分量波场的成像结果可见同向传播的波场之间互相关主要产生低频噪声，波场分解成像条件能够较好地将这部分干扰剔除，反向传播波场之间互相关能够产生有效的成像剖面。对图8a中的成像剖面进行滤波处理，结果如图8b所示。由图可知，不同剖面的低频噪音得到较好地去除。反向传播波场成像

图6 地堑模型不同方程模拟的波场切片(a)伪声波方程(式(1))； (b)优化纯声波方程(式(13)、式(14))

图7 地堑模型分离的上(a)、下(b)、左(c)、右行波场(d)

图8 地堑模型优化纯声波方程不同条件成像结果(a)及其拉普拉斯滤波结果(b)左：全波场互相关； 中：同向传播波场互相关之和； 右：反向传播波场互相关之和

结果之和略优于传统全波场成像结果。相关结果表明，不同方向行波之间的互相关结果可以作为传统方法的补充，成像结果精度更高。

2.2 修改的Marmousi TTI模型

修改的Marmousi TTI模型，如图9所示。将计算区域剖分为500×350个网格单元。用于逆时偏移的速度场和各向异性参数场均基于五点平滑获得。图10为各向同性声波、各向异性伪声波和优化纯声波的成像结果，相比之下，基于优化纯声波方程的逆时偏移能够获得更高精度的成像。图11为基于优化纯声波方程的不同条件成像结果的对比，可见，传统全波场参与的互相关成像，浅层存在严重的低频噪声污染。同向传播的波场之间互相关主要产生低频成像噪声，波场分解成像条件能够较好地将低频噪声去除，反向传播波场之间互相关能够产生高精度的成像结果。相比于经过拉普拉斯滤波的成像剖面而言，基于波场分解成像条件的偏移结果不需要滤波处理，能够较好地保证振幅信息。总体而言，基于波场分解的成像方法可以作为传统方法的有效补充，以获得不同展布方向地下构造的高精度成像。

图9 修改的Marmousi TTI模型(a)vPz； (b)ε； (c)δ； (d)θ

图10 修改的Marmousi TTI模型不同方程的成像结果(拉普拉斯滤波后)(a)各向同性声波方程； (b)伪声波方程(式(1))；(c)优化纯声波方程(式(13)、式(14))

图11 修改的Marmousi TTI模型的优化纯声波方程不同条件成像结果

3 结束语

本文首先回顾了TTI介质中两种常用的伪声波方程，在各向异性参数不满足声学近似时，利用该方程模拟容易产生伪横波干扰及不稳定现象。为了解决现有伪声波方程存在的问题，应用最小二乘方法获得了优化的纯声波频散关系，在此基础上采用一种泊松算法和有限差分相结合的高效精确各向异性纯声波求解算法，能够有效避免伪横波干扰及数值不稳定。将各向同性介质中基于Hilbert变换的复数波场分解成像条件推广到TTI介质优化纯声波逆时偏移成像，选择传播方向相反的分量波场参与最终成像。理论和数值算例验证了本文方法能够有效压制各向异性逆时偏移的伪横波干扰和低频噪声，获得高精度的成像结果。