石油物探  2018, Vol. 57 Issue (4): 555-569  DOI: 10.3969/j.issn.1000-1441.2018.04.009
0
文章快速检索     高级检索

引用本文 

王美霞, 张晓慧, 张钋, 等. 波场分解算法与逆时偏移角道集[J]. 石油物探, 2018, 57(4): 555-569. DOI: 10.3969/j.issn.1000-1441.2018.04.009.
WANG Meixia, ZHANG Xiaohui, ZHANG Po, et al. Wavefield decomposition and RTM angle gathers[J]. Geophysical Prospecting for Petroleum, 2018, 57(4): 555-569. DOI: 10.3969/j.issn.1000-1441.2018.04.009.

作者简介

王美霞(1987—), 女, 硕士, 主要从事波传播和地震成像研究。Email:meixia66.wang@gmail.com

文章历史

收稿日期:2018-04-27
改回日期:2018-05-18
波场分解算法与逆时偏移角道集
王美霞1 , 张晓慧2 , 张钋1 , 唐冰1 , 徐昇1     
1. Statoil Gulf Services, Houston 77042;
2. 前斯塔托尔北京技术服务有限公司, 北京 100022
摘要:逆时偏移生成的角度域共成像点道集(ADCIGs)有利于提高成像质量, 可用于偏移速度分析和AVA/AVO研究等。使用Poynting矢量估计波传播方向可以快速生成角道集, 但当模型和波场比较复杂时, 不同方向传播的波场分量相互重叠, 无法准确计算波的传播方向。波场分解可以提高Poynting矢量的准确性。为此, 研究了波场分解算法, 在计算Poynting矢量之前先将波场分解成不同方向传播的分量, 用于生成角道集。为解决传统的傅里叶变换在频率-波数域进行波场分解时计算量大的问题, 对波场分解算法进行了优化, 并提出了SSE(单指令多数据流扩展)向量化和多线程并行计算实现算法, 从而提高了计算效率。理论和数值计算结果表明:使用波场分解可以显著提高逆时偏移角道集的质量, 此优化算法能够得到与传统方法一致的结果并且提高了计算效率。二维Sigsbee2A模型和三维SEAM TTI模型的数值计算也证明了该算法用于生成高质量角道集和逆时偏移图像的有效性和稳定性。优化后的波场分解算法并不局限于上、下行波分解, 可容易地拓展到其它方向的波场分解, 也可应用于最小二乘偏移和全波形反演, 以消除成像噪声, 提高成像质量。
关键词波场分解    Poynting矢量    角道集    希尔伯特变换    逆时偏移    计算效率    并行计算    
Wavefield decomposition and RTM angle gathers
WANG Meixia1, ZHANG Xiaohui2, ZHANG Po1, TANG Bing1, XU Sheng1     
1. Statoil Gulf Services, Houston 77042, USA;
2. Formerly Statoil (Beijing) Technology Service, Beijing 100022, China
Abstract: Angle domain common image gathers (ADCIGs) from reverse time migration (RTM) can help with AVO/AVA study, improve image quality, and perform velocity analysis.Using the Poynting vectors to estimate wave propagation directions can efficiently generate ADCIGs, but it loses accuracy when the model is complex and different events in the wavefields overlap.Wavefield decomposition can improve the accuracy of Poynting vectors.Therefore, we investigate the wavefield decomposition method, separating the wavefield into different directional components prior to the calculation of Poynting vectors, and apply it to generate ADCIGs.The conventional algorithm is computationally intensive as it is implemented straightforwardly in frequency-wavenumber (f-k) domain by Fourier transforms.To solve this problem, we investigate the optimization of the wavefield decomposition algorithm, and propose to implement it with vectorization by Streaming SIMD Extensions (SSE) and parallelization with multiple threads.The efficiency improvement is about five times.Theoretical and numerical examples showed that using wavefield decomposition can largely improve the quality of ADICGs, and the optimized algorithm can provide equivalent results to the traditional algorithm and improve the computational efficiency.The numerical examples from the 2D Sigsbee2A and 3D SEAM TTI models validate the effectiveness and robustness of the proposed algorithm for high-quality ADCIGs and RTM images.The proposed algorithm is not limited to up-down wavefield decomposition and is easy to extend to wavefield decomposition in other directions.It can be applied to least-squares RTM and FWI to remove artifacts and improve the imaging quality as well.
Key words: wavefield decomposition    Poynting vector    angle gather    Hilbert transform    reverse time migration    computational efficiency    parallel computing    

复杂构造区的地震勘探, 如大角度地层或盐丘等需要高质量的速度模型和高精度的成像方法。逆时偏移[1-4]利用双程波传播方程, 能够准确模拟反射波、折射波、多次波等各种地震波, 是目前最为精确的成像技术之一。

逆时偏移算法在20世纪80年代就已经提出, 但由于其巨大的计算量和超量的计算资源需求, 直到21世纪初才有实际应用。逆时偏移应用面临着理论和计算两方面的挑战。理论方面主要有两点:①逆时偏移低频噪声问题[4]。逆时偏移成像是通过对震源波场和检测器波场应用传统的零延迟互相关成像条件[5]得到的, 如果存在强反射, 这个过程不仅在实际反射界面处产生图像, 也会沿着波路径产生严重的低频噪声。这些噪声可能掩盖真实的构造, 给解释工作造成困难。将这个传统的成像条件和真振幅偏移方法[6-11]做比较, 我们发现主要区别在于依赖于震源和检波器波场夹角的权重函数。实际上, 真振幅偏移中关于反射角的函数可以压制很多反向散射噪声, 关键在于如何计算逆时偏移中的反射角。②快速生成偏移共成像点道集问题。共成像点道集是质量控制、速度建模、属性分析等的重要输入数据, 被广泛应用于复杂地区数据处理。常用的共成像点道集是共偏移距道集, 也就是使用地面偏移距作为索引的道集。然而, 在复杂情况下, 共偏移距道集具有很多由于多路径传播而引入的偏移假像, 而且很难直接由逆时偏移得到高质量的共偏移距道集, 需要进行特殊处理[12]。相比之下, 由逆时偏移产生角道集[13-15]就比较自然。此时角道集由地下反射角作为索引, 偏移噪声也相对较弱。如果使用真振幅逆时偏移, 角道集还可以为AVA分析提供输入数据。角道集可以使用拓展成像条件的方法得到[16-19], 也可以直接由波场产生[11-20]。这些方法在二维情况下计算效率都比较高, 但是在三维情况下计算量都很大, 所以在实现逆时偏移算法时首先要考虑算法效率以及算法的并行计算和优化等。计算效率方面, 引入Poynting矢量方法提高了生成角道集的计算效率, 但是由于检波器波场比较复杂, 一般Poynting矢量仅被用于震源波场, 使用逆时偏移成像或反射界面倾角而不是检波器端波传播方向信息[21-22]。在三维道集逆时偏移应用中, 面对存储量大、计算量高等挑战, Poynting矢量方法更具吸引力。YOON等[23]利用Poynting矢量方法估计波传播方向, 该方法的缺点是在一个时空点处, Poynting矢量只能给出一个方向, 因此, 当波场复杂时, 不同传播方向的波混叠在一起, Poynting矢量就不能准确计算方向。若用于角道集计算, 就会使得有些能量被成像在错误的角度上, 进而降低角道集质量或导致角道集存在部分能量缺失, 这严重影响Poynting矢量方法在检波器端波场的应用。

本文进一步分析了Poynting矢量方法。我们认为只需要在使用Poynting矢量计算波传播方向之前对波场进行分解, 便可以将其同时用于震源和检波器波场。波场分解越精细, Poynting矢量就越精确, 但计算量也会随之增加。对于一般的地面观测系统, 初至波主要是震源的下行波, 反射波主要是检波器的下行波。如果模型中存在强速度差界面, 波场中初至波和反射波就会叠加, 因此进行简单的上、下行波场分解[2-4]就可以显著提高Poynting矢量的精度。

传统的波场分解算法[2-4]一般是通过傅里叶变换[25]在频率-波数域进行的。分解计算效率取决于快速傅里叶变换的实现效率。然而, 傅里叶变换是全局运算, 计算耗时。另一种方法是使用希尔伯特变换(HT)[24, 26-28]进行波场分解。下行波和上行波可以通过原始波场和在时间-深度方向进行希尔伯特变换后的波场组合得到[2-7]。这和利用傅里叶变换得到的结果一致。优点是希尔伯特变换的实现比较灵活; 可以通过快速傅里叶变换实现; 也可以在时间-空间域使用有限长度卷积算子实现, 以便进一步进行算法优化以提高计算效率。本文使用SSE和多线程优化, 提高利用希尔伯特变换进行波场分解的计算效率。

利用希尔伯特变换进行波场分解时, 需要考虑时间和空间方向的希尔伯特变换。有两种方法可以实现波场在时间方向的希尔伯特变换。一种是利用希尔伯特变换后的震源函数进行波传播得到一个新的波场[27], 该方法会加倍波传播计算量和存储量, 因此读写时间也会增加。另一种办法是只使用原始波场, 在计算道集的过程中边计算边对原始波场进行希尔伯特变换。若效率和算法不是瓶颈, 后者需要的数据读写较少。因此, 本文采用后者。

本文首先简单回顾真振幅逆时偏移成像和角道集理论以及Poynting矢量方法。其次, 我们证明通过傅里叶变换或希尔伯特变换都能够进行波场分解, 且理论上等价, 进一步使用单指令多数据流(SIMD)优化希尔伯特变换算法, 并使用数值例子验证了优化后算法的精度和效率。然后, 将优化的希尔伯特算法用于波场分解, 进一步使用多线程优化并行算法。利用数值计算例子说明这种方法在得到和传统傅里叶方法一致结果同时, 能够显著提高数值计算效率。最后, 将这种加速的波场分解算法应用于逆时偏移产生角道集, 利用二维Sigsbee2A和三维SEAM TTI计算实例说明本文方法的有效性和稳定性。

1 方法理论

利用双程波方程进行波传播的数值计算[29-35], 可以得到震源波场us(x, t)和检波器波场ur(x, t)。这里t表示时间, 向量x表示空间位置。在二维介质中, x=(z, x)。在三维介质中, x=(z, x, y)。下面我们首先回顾成像条件和角道集理论, 之后引出一种快速的波场分解算法。

1.1 成像条件

偏移成像的传统互相关成像条件[5]可以表示为:

$ \begin{array}{l} I\left( x \right) = \int {\int {\int {u_{\rm{s}}^ * \left( {\mathit{\boldsymbol{x}},\omega ;{x_{\rm{s}}},{y_{\rm{s}}}} \right) \cdot {u_{\rm{r}}}\left( {\mathit{\boldsymbol{x}},\omega ;{x_{\rm{s}}},{y_{\rm{s}}}} \right) \cdot } } } \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{\rm{d}}\omega {\rm{d}}{x_{\rm{s}}}{\rm{d}}{y_{\rm{s}}} \end{array} $ (1)

式中:xsys表示震源位置; ω为角频率。此成像条件在反射点处产生像, 但同时也在满足成像条件的路径上产生假像。假像和真像的区别在于在假像处震源波场和检波器波场的传播方向相异。YOON等[23]提出只有当震源和检波器波传播方向所成角度在一定范围内时才使用互相关条件成像, 且波传播方向由Poynting矢量计算得到, 该方法的一个缺点是很难选取用于成像的最大角度。COSTA等[36]提出了依赖于角度的光滑函数作为偏移权重, 在一定程度上克服了最大成像角度选取的困难。与真振幅偏移[6-7, 10-11]相比, 其和传统成像条件的区别在于依赖于震源和检波器波场夹角的权重函数。因此, 一个更加适当的成像条件可以表达为[6, 11, 36-38]:

$ \begin{array}{l} I\left( \mathit{\boldsymbol{x}} \right) = \int {\int {\int {u_{\rm{s}}^ * \left( {\mathit{\boldsymbol{x}},\omega ;{x_{\rm{s}}},{y_{\rm{s}}}} \right) \cdot {u_{\rm{r}}}\left( {x,\omega ;{x_{\rm{s}}},{y_{\rm{s}}}} \right) \cdot } } } \\ \;\;\;\;\;\;\;\;\;\;\;{\cos ^2}\theta {\rm{d}}\omega {\rm{d}}{x_{\rm{s}}}{\rm{d}}{y_{\rm{s}}} \end{array} $ (2)

式中:θ表示反射角。由于偏移噪声是当震源和检波器波传播方向成较大角度时产生的, 所以权重函数cos2θ可以有效压制偏移中的低频噪声。

1.2 使用Poynting矢量产生真振幅逆时偏移角道集

共炮逆时偏移的真振幅三维角度域共成像点道集[11, 39, 40]可以按下式计算:

$ \begin{array}{l} R\left( {\mathit{\boldsymbol{x}},{\theta _0},{\varphi _0}} \right) = \int {\int {\int {\delta \left( {\theta - {\theta _0}} \right)\delta \left( {\varphi - {\varphi _0}} \right)v\left( \mathit{\boldsymbol{x}} \right)\frac{{{{\cos }^2}\theta }}{{\sin \theta }} \cdot } } } \\ \;\;\;\;u_{\rm{s}}^ * \left( {\mathit{\boldsymbol{x}},\omega ;{x_{\rm{s}}},{y_{\rm{s}}}} \right) \cdot {u_{\rm{r}}}\left( {\mathit{\boldsymbol{x}},\omega ;{x_{\rm{s}}},{y_{\rm{s}}}} \right){\rm{d}}\omega {\rm{d}}{x_{\rm{s}}}{\rm{d}}{y_{\rm{s}}} \end{array} $ (3)

式中:R(x, θ0, φ0)表示在空间位置x处方位角φ0和反射角θ0方向上的反射系数; δ(θ-θ0)和δ(φ-φ0)指狄拉克函数, 它们将积分限制在要求的角度θ0φ0, v(x)为速度。

反射角和方位角计算如下:

$ \cos \left( \theta \right) = \frac{{\mathit{\boldsymbol{p}} \cdot {\mathit{\boldsymbol{p}}_{\rm{r}}}}}{{\left| \mathit{\boldsymbol{p}} \right| \cdot \left| {{\mathit{\boldsymbol{p}}_{\rm{r}}}} \right|}} $ (4)
$ \cos \left( \varphi \right) = \frac{{\left[ {\left( {{\mathit{\boldsymbol{p}}_{\rm{s}}} \times {\mathit{\boldsymbol{p}}_{\rm{r}}}} \right) \times \mathit{\boldsymbol{z}}} \right] \cdot \mathit{\boldsymbol{x}}}}{{\left| {\left( {{\mathit{\boldsymbol{p}}_{\rm{s}}} \times {\mathit{\boldsymbol{p}}_{\rm{r}}}} \right) \times \mathit{\boldsymbol{z}}} \right|}} $ (5)

其中, p=1/2(ps+pr), x, z为坐标轴的单位向量; pspr分别表示震源和检波器的波传播方向。震源波传播方向可以按下式利用Poynting矢量进行计算:

$ {\mathit{\boldsymbol{p}}_{\rm{s}}}\left( {\mathit{\boldsymbol{x}},t} \right) = - \frac{{\partial {u_{\rm{s}}}\left( {\mathit{\boldsymbol{x}},t} \right)}}{{\partial t}} \cdot \nabla {u_{\rm{s}}}\left( {\mathit{\boldsymbol{x}},t} \right) $ (6)

类似地, 检波器波传播方向pr也可以使用Poynting矢量计算得到。

Poynting矢量计算是在时间-空间域进行的, 当波场比较复杂时, Poynting矢量并不精确。图 1a展示了从Sigsbee2A模型得到的一个波场快照, 可以看到波场很复杂, 几种不同类型的波交织在一起, 这使得Poynting矢量很难给出一个正确的方向。然而, 如果我们仅关注下行波(图 1b), 波场就变得简单很多, Poynting矢量可以较准确地计算波传播方向。因此, 波场分解, 特别是上、下行波场分解(对于地表观测系统), 能显著提高Poynting矢量计算方向的精度, 这也可以进一步提高角道集的质量。

图 1 由Sigsbee2A得到的波场快照(a)及其下行波场分量(b)
1.3 波场分解 1.3.1 利用傅里叶变换进行波场分解

传统方法利用傅里叶变换在频率-波数域进行波场分解。约定时间-深度方向的二维傅里叶变换如下:

$ U\left( {{k_z},\omega } \right) = \int {\int {u\left( {z,t} \right){{\rm{e}}^{ - {\rm{i}}\omega t - {k_z}z}}{\rm{d}}t{\rm{d}}z} } $ (7)

于是下行波在频率-波数域可表示为:

$ {U_{\rm{D}}}\left( {{k_z},\omega } \right) = \left\{ {\begin{array}{*{20}{c}} {U\left( {{k_z},\omega } \right)}&{\omega \cdot {k_z} < 0}\\ {\frac{1}{2}U\left( {{k_z},\omega } \right)}&{\omega \cdot {k_z} = 0}\\ 0&{\omega \cdot {k_z} > 0} \end{array}} \right. $ (8)

上行波在频率-波数域可表示为:

$ \begin{array}{l} {U_{\rm{U}}}\left( {{k_z},\omega } \right) = U\left( {{k_z},\omega } \right) - {U_{\rm{D}}}\left( {{k_z},\omega } \right)\\ \;\;\;\;\;\; = \left\{ {\begin{array}{*{20}{c}} {U\left( {{k_z},\omega } \right)}&{\omega \cdot {k_z} > 0}\\ {\frac{1}{2}U\left( {{k_z},\omega } \right)}&{\omega \cdot {k_z} = 0}\\ 0&{\omega \cdot {k_z} < 0} \end{array}} \right. \end{array} $ (9)

在公式(8)和公式(9)中, 对于满足ω·kz=0的U(kz, ω), 我们使用了系数1/2以便使得UD(kz, ω)+UU(kz, ω)=U(kz, ω)。这是合理的, 因为满足ω·kz=0的波既非上行波也非下行波。

1.3.2 利用希尔伯特变换进行波场分解

SHEN等[27]指出上、下行波可以使用希尔伯特变换进行波场分离。他们使用一对波场实现分解:震源函数传播得到的原始波场和进行希尔伯特变换后的震源函数传播得到的波场。这种方法很有效, 但计算效率不高, 它需要两倍的波传播计算量。本文提出另一种实现方法, 仅使用原始波场, 在计算道集过程中进行希尔伯特变换。重点以上、下行波分解来描述算法, 并将该方法拓展到其它方向的波场分解。

通过傅里叶变换和希尔伯特变换之间的关系, 最终下行波可以表示为:

$ {u_{\rm{d}}}\left( {\mathit{\boldsymbol{x}},t} \right) = \frac{1}{2}u\left( {\mathit{\boldsymbol{x}},t} \right) + \frac{1}{2}{H_t}{H_z}\left[ {u\left( {\mathit{\boldsymbol{x}},t} \right)} \right] $ (10)

上行波可以表示为:

$ {u_{\rm{u}}}\left( {\mathit{\boldsymbol{x}},t} \right) = \frac{1}{2}u\left( {\mathit{\boldsymbol{x}},t} \right) - \frac{1}{2}{H_t}{H_z}\left[ {u\left( {\mathit{\boldsymbol{x}},t} \right)} \right] $ (11)

式中:Ht[u(x, t)]表示波场u(x, t)在时间方向的希尔伯特变换; Hz[u(x, t)]表示深度方向的希尔伯特变换。公式(10)和公式(11)与SHEN等[27]提出的公式一致, 具体推导过程见附录A。

公式(10)和公式(11)提供了另一种波场分解的办法。希尔伯特变换可以简单通过傅里叶变换实现, 但是这样计算效率与传统频率-波数域的波场分解一样。为了提高计算效率, 我们提出直接在时间-空间域通过卷积进行希尔伯特变换, 并进一步使用SSE优化该算法。

对于一维输入数组aj, j=0, 1, …, n-1, 使用(2m+1)长度卷积算子进行的离散希尔伯特变换如下:

$ {b_j} = \sum\limits_{k = 0}^{2m} {{h_k}{a_{j + m - k}}} ,\;\;\;\;\;\;\;j = 0,1, \cdots ,n - 1 $ (12)

其中, ${h_k} = \left\{ \begin{array}{l} 0,\left( {m - k} \right){为偶数}\\ \frac{2}{{{\rm{ \mathsf{ π} }}\left( {m - k} \right)}},\left( {m - k} \right){为奇数} \end{array} \right.$。数组{bj}为数组{aj}的希尔伯特变换。在公式(12)中, 如果j+m-k < 0或j+m-k>n-1, 则aj+m-k=0。考虑到计算效率, 可选取m为一个相对小的整数。公式(12)中数据的截断会引入数值噪声, 可以使用Tukey窗进行平滑处理, 以消除这些边界效应。

我们注意到进行希尔伯特变换的卷积算子是反对称的, 即:hk=-h2m-k, k=0, 1, …, m。而且大约一半的系数hk都是0。基于这两条性质, 公式(12)可以进一步改写为如下形式以减少计算量:

$ \begin{array}{l} {b_j} = \sum\limits_{k = m + 1,k + = 2}^{2m} {{h_{2m - k}}\left( {{a_{j - m + k}} - {a_{j + m - k}}} \right)} ,\\ \;\;\;\;\;\;j = 0,1, \cdots ,n - 1 \end{array} $ (13)

式中:k+=2表示k的步长为2。

以下的实现和数值计算实例都是基于公式(13)。

1.3.3 希尔伯特变换波场分解算法的优化

若将波场分解用于逆时偏移中进行上、下行波分解, 我们需要时间和深度方向的希尔伯特变换(如公式(10)和公式(11))。约定波场在内存中的存储规律为最快维为深度, 最慢维为时间。因此, 提高波场分解算法速度的关键是分别优化快维和慢维方向的希尔伯特变换。下面先说明如何优化二维数组的快维和慢维的希尔伯特变换, 然后, 将其用于波场分解。

给定二维数组ai, j, i=0, 1, …, n2-1, 且j=0, 1, …, n1-1(下角标i表示慢维索引, j表示快维索引), 记数组{ai, j}的希尔伯特变换为数组{bi, j}。传统方法在快维方向的希尔伯特变换可直接表示为:

$ {b_{i,j}} = \sum\limits_{k = m + 1,k + = 2}^{2m} {{h_{2m - k}}\left( {{a_{i,j - m + k}} - {a_{i,j + m - k}}} \right)} $ (14)

若要计算慢维方向的希尔伯特变换, 传统方法是先通过转置将慢维变为快维, 再由公式(14)计算, 最后再做一次数组转置。对数据做两次转置增加了计算量。我们提出一种不用转置的算法, 即按照公式(15)直接对慢维进行数值运算:

$ {b_{i,j}} = \sum\limits_{k = m + 1,k + = 2}^{2m} {{h_{2m - k}}\left( {{a_{i - m + k,j}} - {a_{i + m - k,j}}} \right)} $ (15)

由于我们使用相对较短的卷积算子, k相对于慢维长度n2较小, 就可以使用SIMD沿快维j方向加速算法, 即将每一列ai-m+k, *-ai+m-k, *看作一个元素, 同一个系数h2m-k与该元素相乘, 这样就可以使用SSE同时处理多个数据。我们将这种优化的算法记为慢维方向的向量化希尔伯特变换(vectorized Hilbert transform, VHT)。

通过使用SSE优化, 慢维方向的希尔伯特变换不再需要对数组进行转置。理论上, 数值计算效率可以提高到4倍(SSE)或8倍(AVX, Advanced Vetor Extensions)。但是考虑到从缓存中获取数据也需要时间, 实际提高水平可能略低于理论值。

对于快维方向的希尔伯特变换, 我们发现如下使用SSE优化的算法比直接利用公式(14)计算更加高效:

1) 将输入数组ai, j做转置得到aj, i, i=0, 1, …, n2-1, j=0, 1, …, n1-1;

2) 使用VHT计算aj, i在慢维方向的希尔伯特变换:

$ {{\bar b}_{j,i}} = \sum\limits_{k = m + 1,k + = 2}^{2m} {{h_{2m - k}}\left( {{{\bar a}_{j - m + k,i}} - {{\bar a}_{j + m - k,i}}} \right)} $ (16)

3) 将bj, i进行转置得到bi, j

在以上过程中, 二维转置可以使用优化的SSE转置函数(_MM_TRANSPOSE4_PS)实现。尽管使用了转置, SSE向量化仍然可以大幅提高计算效率。我们称此为快维方向的VHT。

由于重点是上、下行波分解, 所以进行波场分解算法的核心是二维波场的分解。对于三维情况, y方向是最慢维, 只需对每一个y切面波场进行二维波场分解。在使用SSE向量化的同时, 还可以使用OpenMP进行多线程并行计算以提高效率。对于地震成像应用, 目前的缓存一般足够容纳二维子波场。因此, 基于公式(10)和公式(11), 我们提出以下流程进行高效波场分解:

1) 对每一个ix, 将二维子波场u(z, ix, t)从内存加载到缓存中, 并使用慢维方向的向量化希尔伯特变换计算Ht[u(z, ix, t)];

2) 对每一个it, 将二维子波场Ht[u(z, x, it)]从内存加载到缓存中, 并使用快维方向的VHT计算HzHt[u(z, x, it)];

3) 使用公式(10)(或公式(11))得到下行(或上行)波场。

对于三维波场的分解, 我们只需要对每一个iy切片进行以上二维波场分解。

1.3.4 数值算例

先以一个二维波场快照u(z, x)为例说明VHT的计算效率及其在波场分解中的应用。输入波场是由常速度v=1790m/s模型通过伪谱法[41]模拟产生, 震源函数是主频为14Hz的Ricker子波, 我们对子波进行了3Hz的高通滤波。波传播网格为dx=dz=15m, nz=nx=701, 且深度z为快维。震源在区域中心。我们使用一个相对较长的卷积算子, m=40, 即81个点, 进行希尔伯特变换以便得到较好的精度。

首先, 我们验证与传统方法相比, 慢维方向的VHT能够提高计算效率。图 2t=2s时刻的波场快照。图 3a是使用向量化希尔伯特变换计算得到的二维波场在慢维x方向的希尔伯特变换波场。图 3b比较了慢维方向的VHT和传统希尔伯特变换的计算效率。为了得到可靠的计算时间统计, 我们重复了5000次计算。从图 3b可以看到, 对于慢维方向的希尔伯特变换, 向量化算法比传统算法效率高三倍多。为检验VHT的精度, 我们取图 3a中一条深度方向的记录与傅里叶方法得到的结果进行比较, 见图 3c。由于傅里叶方法对于所有波数都是精确的, 可以认为其得到的结果为理想结果。图 3c中VHT结果与傅里叶方法的结果十分吻合, 说明本文的VHT方法准确有效。在进行VHT时, 卷积算子越长, 精度越高。在实际应用中, 如果不需要特别高的精度, 可以使用一个较短的卷积算子, 这样可以进一步提高计算效率。

图 2 t=2s时刻的二维波场快照
图 3 使用VHT得到的波场快照(a)、VHT和传统希尔伯特变换(HT)计算效率对比(b)以及深度z=5250m处VHT结果和傅里叶变换结果(FT)比较(c)

其次, 我们再验证快维方向的VHT。图 4a图 3a所示波场快照在深度方向的希尔伯特变换。图 4b比较了VHT方法和传统HT方法用于计算快维方向希尔伯特变换的计算效率。我们同样进行了5000次重复计算。可以看到向量化方法可以减少大概一半的计算时间。

图 4 快维(z方向)的希尔伯特变换 a 使用VHT得到的波场快照在z方向的希尔伯特变换; b z方向的VHT和传统希尔伯特变换(HT)的效率比较(计算重复了5000次)

从以上两个实例(图 3b图 4b), 我们可以看到:传统HT方法在计算慢维方向希尔伯特变换时比计算快维方向希尔伯特变换时需要时间更多; 相反地, 向量化方法在计算快维方向希尔伯特变换时比计算慢维方向希尔伯特变换时需要更多时间。这均是由于是否需要数据转置所致。

最后, 我们进一步检验VHT用于波场分解的综合计算效率, 并将其和传统的利用傅里叶变换的波场分解进行比较。为了合理地进行比较, 在进行傅里叶变换时使用了FFTW(the Faster Fourier Transform in the West)[42]。波场时间方向采样为:nt=660, dt=4ms。图 5a图 5b分别是由本文算法和傅里叶方法得到的下行波场在t=2s时的波场快照。图 5c比较了本文方法和傅里叶变换方法的计算效率。两种方法都是用16核计算机进行计算。为得到可靠的计算时间统计, 运算重复了100次。可以看到本文提出的使用VHT的分解方法能够很大地提高计算效率, 比传统使用傅里叶变换进行分解的方法快了6倍多。为进行更详细的精度比较, 图 6比较了从图 5a图 5b中取出的x=4950m处一条深度方向的记录, 可以看到两种方法得到的结果非常一致。

图 5 上、下行波分解示例(采用t=2s时刻的下行波分量) a VHT结果; b FFTW计算结果; c 两者效率比较(计算重复100次)
图 6 下行波场中x= 4950m处深度方向记录的比较
1.3.5 波场分解拓展

以上主要阐述了如何进行上、下行波分解, 我们也可以将该方法拓展到其他方向的波场分解。为了理论的完整性, 附录B列出了其他方向波分量的公式, 包括左行波(ul)、右行波(ur)、左上行波(uul)、右上行波(uur)、左下行波(udl)、右下行波(udr)、左上后行波(uulb)、左上前行波(uulf)、右上后行波(uurb)、右上前行波(uurf)、左下后行波(udlb)、左下前行波(udlf)、右下后行波(udrb)、右下前行波(udrf)。

同理, 用于上、下行波分解的优化方法也可以用于其它方向波场的分解。这里给出一个左、右行波场分解的例子来说明其有效性。我们同样使用图 2所使用的模型。图 7a图 7b分别是t=2s时刻的左行波和右行波。图 7c比较了两种方法的计算效率。可以看到, 本文提出的VHT方法能够显著提高左、右行波分解的计算效率。值得注意的是本文方法用于左、右行波分解时不需要任何数据转置, 这一点也优于传统方法。

图 7 向量化希尔伯特变换得到的t=2s时刻的左、右行波分解结果 a 左行波分量; b 右行波分量; c 两者效率比较(计算重复100次)
2 应用实例

将波场分解算法用于逆时偏移成像来检验其在逆时偏移和角道集生成中的作用。偏移像可以通过叠加角道集得到, 这使得成像更加灵活, 例如可以只叠加某一角度范围内的道集或使用合适的关于角度的权重函数。下面介绍了二维Sigsbee2A模型和三维SEAM TTI模型例子。

2.1 二维Sigsbee2A例子

首先使用Sigsbee2A模型和数据。偏移中共使用500炮的合成地震数据。最大频率为36Hz。为了准确模拟波场, 我们使用伪谱法[41]进行波传播, 并对时间频散进行了校正[34]图 8a是不使用波场分解时得到的角道集, 可以观察到盐丘下边界之上道集并不连续也不清晰。图 8b为使用VHT进行波场分解得到的结果, 道集连续性变好, 噪声也得到压制。这说明波场分解能够提高Poynting矢量方向计算的精度。图 9是使用互相关成像条件(方程(1))和未进行分解的波场得到的叠加图像, 从图上可以看到严重的低频噪声, 以致掩盖了真实的反射界面。图 10a是使用包含角度权重的成像条件(方程(2))和未进行分解的波场得到的叠加图像, 图 10b是使用包含角度权重的成像条件(方程(2))和下行波得到的叠加图像。比较图 9图 10a, 可以看到依赖于角度的权重函数可以压制很多噪声, 盐丘以下的反射层也变得清楚。但是图 10a中还存在一些噪声, 使盐丘下边界之上的图像比较模糊, 如图中红色箭头所示区域。若进一步使用波场分解, 叠加图像就变得比较清晰, 质量得到提高, 如图 10b所示。这也说明波场分解之后使用包含角度权重的成像条件(方程(2))可以有效降低偏移噪声。

图 8 Sigsbee2A模型的角道集(道集抽取位置6004.56~25816.56m, 间隔1524m) a 不使用波场分解得到的角道集; b 使用VHT进行波场分解得到的角道集
图 9 利用互相关成像条件得到的二维Sigsbee2A模型叠加图像(没有使用波场分解)
图 10 利用0~60°角道集叠加得到二维Sigsbee2A模型叠加图像(成像条件为依赖于角度的权重函数) a 未使用波场分解; b 采用VHT对下行波场进行波场分解
2.2 三维SEAM TTI例子

实际上, 我们更加重视三维应用的计算效率。为此, 我们使用三维SEAM (SEG Advanced Modeling) TTI模型来说明本文算法在三维逆时偏移应用中的计算效率。SEAM模型区域为40km×35km×15km。包含一个复杂的盐丘[43], 这在墨西哥湾地区非常典型。不仅模型大小符合实际情况, 合成地震数据也比较符合真实情形。由于计算资源有限, 我们选取了215炮数据进行测试。图 11展示了炮的位置。在算法中, 我们使用了XU等[33]提出的拟P波方程, 这样不仅效率高而且没有S波产生的噪声。测试中, 偏移最高频率为18Hz。以下我们选取位置y=32250m处观测偏移结果。图 12a图 12b分别是不使用和使用波场分解得到的结果。可以看到:图 12a中, 特别是红色椭圆标记区域, 丢失了很多成像细节; 图 12b相对清晰很多, 也有更多细节信息, 特别是在盐丘之上。这个比较说明即便对于非常复杂的三维模型, 波场分解也能够提高Poynting矢量的精度, 从而提高角道集的质量。图 13是不使用波场分解得到的叠加偏移图像(使用成像条件(1)), 可以看到严重的噪声, 真实的构造被掩盖。图 14是使用成像条件(2)得到的图像, 图 14a图 14b分别使用波场分解前后得到的结果。比较图 14a图 13, 可见关于角度的权重函数能够压制一部分噪声, 但是图像质量还不是很高。进一步使用波场分解, 可以得到更加清晰、细节更丰富的图像(图 14b)。因此依赖于角度的权重函数可以压制一些噪声, 波场分解可以进一步提高成像质量, 去除偏移噪声。

图 11 三维SEAM TTI模型实例(选取215炮的位置)
图 12 三维SEAM TTI模型角道集(角道集位置11~25km, 间隔为1km) a 未进行波场分解得到的角道集; b 利用VHT波场分解后得到的角道集
图 13 三维SEAM TTI模型叠加偏移图像(采用互相关成像条件(方程(1)), 未进行波场分解)
图 14 利用0~60°的角道集叠加得到三维SEAM TTI模型叠加偏移图像(采用依赖于角度的权重函数(公式(2)成像) a未使用波场分解; b 采用VHT进行波场分解
3 结论

角度域共成像点道集对复杂构造的地震成像十分重要, 可以使用Poynting矢量高效计算波传播方向, 但是当波场复杂时, 传统的Poynting方法精度较低。为提高Poynting矢量的精度, 可以在计算Poynting矢量之前将波场进行分解, 再进行Poynting矢量计算, 得到高精度的传播方向, 从而得到高质量的角道集和偏移图像。我们提出了VHT和多核并行的波场分解算法。经过此优化, 波场分解算法加快大约5倍。二维Sigsbee2A和三维SEAM TTI的数值计算实例均说明本文的波场分解算法能够有效去除逆时偏移图像和角道集中的噪声, 提高精度和计算效率。该方法并不局限于上、下行波场分解, 可以很容易拓展到其它方向的波场分解。本文算法也可以进一步应用于最小二乘偏移和全波形反演以消除成像噪声和提高成像质量。

附录A 公式(10)和公式(11)的推导

我们首先定义如下两个辅助函数:

$ {g_1}\left( \xi \right) = \left\{ {\begin{array}{*{20}{c}} 1&{\xi > 0}\\ {\frac{1}{2}}&{\xi = 0}\\ 0&{\xi < 0} \end{array}} \right. $ (A1)
$ {g_2}\left( \xi \right) = \left\{ {\begin{array}{*{20}{c}} 0&{\xi > 0}\\ {\frac{1}{2}}&{\xi = 0}\\ 1&{\xi < 0} \end{array}} \right. $ (A2)

通过使用辅助函数(A1)式和(A2)式, 可以将频率-波数域下行波的定义公式(8)改写成如下形式:

$ {U_{\rm{D}}}\left( {{k_z},\omega } \right) = U\left( {{k_z},\omega } \right) \cdot \left[ {{g_1}\left( \omega \right) \cdot {g_2}\left( {{k_z}} \right) + {g_2}\left( \omega \right) \cdot {g_1}\left( {{k_z}} \right)} \right] $ (A3)

同理, 使用辅助函数(A1)式和(A2)式, 可以将频率-波数域上行波的定义公式(9)改写成如下形式:

$ {U_{\rm{U}}}\left( {{k_z},\omega } \right) = U\left( {{k_z},\omega } \right) \cdot \left[ {{g_1}\left( \omega \right) \cdot {g_1}\left( {{k_z}} \right) + {g_2}\left( \omega \right) \cdot {g_2}\left( {{k_z}} \right)} \right] $ (A4)

对于任意信号, 其希尔伯特变换和傅里叶变换之间存在如下关系:

$ {F_{\rm{h}}}\left( \xi \right) = F\left( \xi \right) \cdot \sigma \left( \xi \right) $ (A5)

式中:ξ表示频率或波数, F(ξ)表示信号的傅里叶变换, Fh(ξ)表示信号的希尔伯特变换。根据我们关于希尔伯特变换和傅里叶变换的定义, σ(ξ)是如下函数:

$ \sigma \left( \xi \right) = \left\{ {\begin{array}{*{20}{c}} i&{\xi > 0}\\ 0&{\xi = 0}\\ { - i}&{\xi < 0} \end{array}} \right. $ (A6)

根据函数定义(A1)式、(A2)式和(A6)式, 我们可得到以下g1(ξ), g2(ξ)关于σ(ξ)的表示:

$ {g_1}\left( \xi \right) = \frac{1}{2} + \frac{1}{{2i}}\sigma \left( \xi \right) $ (A7)
$ {g_2}\left( \xi \right) = \frac{1}{2} - \frac{1}{{2i}}\sigma \left( \xi \right) $ (A8)

将(A7)式和(A8)式代入(A3)式可以得到下行波的表达式:

$ {U_{\rm{D}}}\left( {{k_z},\omega } \right) = \frac{1}{2}U\left( {{k_z},\omega } \right) + \frac{1}{2}\sigma \left( \omega \right) \cdot \sigma \left( {{k_z}} \right) \cdot U\left( {{k_z},\omega } \right) $ (A9)

同理, 将(A7)式和(A8)式代入(A4)式可以得到上行波的表达式:

$ {U_{\rm{U}}}\left( {{k_z},\omega } \right) = \frac{1}{2}U\left( {{k_z},\omega } \right) - \frac{1}{2}\sigma \left( \omega \right) \cdot \sigma \left( {{k_z}} \right) \cdot U\left( {{k_z},\omega } \right) $ (A10)

(A9)式和(A10)式是下行波和上行波在频率-波数域的表达式。再使用傅里叶变换和希尔伯特变换的关系(A5)式, 可以观察到下行波和上行波都可表示为原始波场及其在时间和深度方向进行希尔伯特变换后波场的组合。因此, 在时间-空间域中, 下行波场和上行波场可以分别表示为公式(10)和公式(11)。

附录B 波场分解拓展

如公式(10)和公式(11)所示, 上、下行波场都可以通过希尔伯特变换得到。在此附录中, 我们将此理论加以延伸。首先推导左行波、右行波和左上行波的表达式。其它方向波分量可以同理得到, 我们仅给出最终表达式。

在频率-波数域中, 右行波可以定义为:

$ {U_{\rm{R}}}\left( {{k_x},\omega } \right) = \left\{ {\begin{array}{*{20}{c}} {U\left( {{k_x},\omega } \right)}&{\omega \cdot {k_x} < 0}\\ {\frac{1}{2}U\left( {{k_x},\omega } \right)}&{\omega \cdot {k_x} = 0}\\ 0&{\omega \cdot {k_x} > 0} \end{array}} \right. $ (B1)

左行波可以定义为:

$ {U_{\rm{L}}}\left( {{k_x},\omega } \right) = \left\{ {\begin{array}{*{20}{c}} {U\left( {{k_x},\omega } \right)}&{\omega \cdot {k_x} > 0}\\ {\frac{1}{2}U\left( {{k_x},\omega } \right)}&{\omega \cdot {k_x} = 0}\\ 0&{\omega \cdot {k_x} < 0} \end{array}} \right. $ (B2)

类似附录A中上、下行波的推导, 可以得到以下由希尔伯特变换表示的左、右行波表达式:

$ {u_{\rm{r}}}\left( {\mathit{\boldsymbol{x}},t} \right) = \frac{1}{2}u\left( {\mathit{\boldsymbol{x}},t} \right) + \frac{1}{2}{H_t}{H_x}\left[ {u\left( {\mathit{\boldsymbol{x}},t} \right)} \right] $ (B3)
$ {u_{\rm{l}}}\left( {\mathit{\boldsymbol{x}},t} \right) = \frac{1}{2}u\left( {\mathit{\boldsymbol{x}},t} \right) - \frac{1}{2}{H_t}{H_x}\left[ {u\left( {\mathit{\boldsymbol{x}},t} \right)} \right] $ (B4)

左上行波可以通过进一步将上行波uu(x, t)进行左、右行波分解得到:

$ {u_{{\rm{ul}}}}\left( {\mathit{\boldsymbol{x}},t} \right) = \frac{1}{2}{u_u}\left( {\mathit{\boldsymbol{x}},t} \right) - \frac{1}{2}{H_t}{H_x}\left[ {{u_u}\left( {\mathit{\boldsymbol{x}},t} \right)} \right] $ (B5)

将上行波表达式(公式(11))代入此方程可以得到:

$ {u_{{\rm{ul}}}}\left( {\mathit{\boldsymbol{x}},t} \right) = \frac{1}{4}\left\{ {u\left( {\mathit{\boldsymbol{x}},t} \right) - {H_t}{H_z}\left[ {u\left( {\mathit{\boldsymbol{x}},t} \right)} \right] - {H_t}{H_x}\left[ {u\left( {\mathit{\boldsymbol{x}},t} \right)} \right] - {H_x}{H_z}\left[ {u\left( {\mathit{\boldsymbol{x}},t} \right)} \right]} \right\} $ (B6)

同样, 右上行波可以表达为:

$ {u_{{\rm{ur}}}}\left( {\mathit{\boldsymbol{x}},t} \right) = \frac{1}{4}\left\{ {u\left( {\mathit{\boldsymbol{x}},t} \right) - {H_t}{H_z}\left[ {u\left( {\mathit{\boldsymbol{x}},t} \right)} \right] + {H_t}{H_x}\left[ {u\left( {\mathit{\boldsymbol{x}},t} \right)} \right] + {H_x}{H_z}\left[ {u\left( {\mathit{\boldsymbol{x}},t} \right)} \right]} \right\} $ (B7)

左下行波可以表达为:

$ {u_{{\rm{dl}}}}\left( {\mathit{\boldsymbol{x}},t} \right) = \frac{1}{4}\left\{ {u\left( {\mathit{\boldsymbol{x}},t} \right) + {H_t}{H_z}\left[ {u\left( {\mathit{\boldsymbol{x}},t} \right)} \right] - {H_t}{H_x}\left[ {u\left( {\mathit{\boldsymbol{x}},t} \right)} \right] + {H_x}{H_z}\left[ {u\left( {\mathit{\boldsymbol{x}},t} \right)} \right]} \right\} $ (B8)

右下行波可以表达为:

$ {u_{{\rm{dr}}}}\left( {\mathit{\boldsymbol{x}},t} \right) = \frac{1}{4}\left\{ {u\left( {\mathit{\boldsymbol{x}},t} \right) + {H_t}{H_z}\left[ {u\left( {\mathit{\boldsymbol{x}},t} \right)} \right] + {H_t}{H_x}\left[ {u\left( {\mathit{\boldsymbol{x}},t} \right)} \right] - {H_x}{H_z}\left[ {u\left( {\mathit{\boldsymbol{x}},t} \right)} \right]} \right\} $ (B9)

左上后行波可以表达为:

$ \begin{array}{l} {u_{{\rm{ulb}}}}\left( {\mathit{\boldsymbol{x}},t} \right) = \\ \frac{1}{8}\left\{ {u\left( {\mathit{\boldsymbol{x}},t} \right) - {H_t}{H_z}\left[ {u\left( {\mathit{\boldsymbol{x}},t} \right)} \right] - {H_t}{H_x}\left[ {u\left( {\mathit{\boldsymbol{x}},t} \right)} \right] - {H_t}{H_y}\left[ {u\left( {\mathit{\boldsymbol{x}},t} \right)} \right] - {H_x}{H_z}\left[ {u\left( {\mathit{\boldsymbol{x}},t} \right)} \right] - } \right.\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left. {{H_x}{H_y}\left[ {u\left( {\mathit{\boldsymbol{x}},t} \right)} \right] - {H_y}{H_z}\left[ {u\left( {\mathit{\boldsymbol{x}},t} \right)} \right] + {H_t}{H_x}{H_y}{H_z}\left[ {u\left( {\mathit{\boldsymbol{x}},t} \right)} \right]} \right\} \end{array} $ (B10)

左上前行波可以表达为:

$ \begin{array}{l} {u_{{\rm{ulf}}}}\left( {\mathit{\boldsymbol{x}},t} \right) = \\ \frac{1}{8}\left\{ {u\left( {\mathit{\boldsymbol{x}},t} \right) - {H_t}{H_z}\left[ {u\left( {\mathit{\boldsymbol{x}},t} \right)} \right] - {H_t}{H_x}\left[ {u\left( {\mathit{\boldsymbol{x}},t} \right)} \right] + {H_t}{H_y}\left[ {u\left( {\mathit{\boldsymbol{x}},t} \right)} \right] - {H_x}{H_z}\left[ {u\left( {\mathit{\boldsymbol{x}},t} \right)} \right] + } \right.\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left. {{H_x}{H_y}\left[ {u\left( {\mathit{\boldsymbol{x}},t} \right)} \right] + {H_y}{H_z}\left[ {u\left( {\mathit{\boldsymbol{x}},t} \right)} \right] - {H_t}{H_x}{H_y}{H_z}\left[ {u\left( {\mathit{\boldsymbol{x}},t} \right)} \right]} \right\} \end{array} $ (B11)

右上后行波可以表达为:

$ \begin{array}{l} {u_{{\rm{urb}}}}\left( {\mathit{\boldsymbol{x}},t} \right) = \\ \frac{1}{8}\left\{ {u\left( {\mathit{\boldsymbol{x}},t} \right) - {H_t}{H_z}\left[ {u\left( {\mathit{\boldsymbol{x}},t} \right)} \right] + {H_t}{H_x}\left[ {u\left( {\mathit{\boldsymbol{x}},t} \right)} \right] - {H_t}{H_y}\left[ {u\left( {\mathit{\boldsymbol{x}},t} \right)} \right] + {H_x}{H_z}\left[ {u\left( {\mathit{\boldsymbol{x}},t} \right)} \right] + } \right.\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left. {{H_x}{H_y}\left[ {u\left( {\mathit{\boldsymbol{x}},t} \right)} \right] - {H_y}{H_z}\left[ {u\left( {\mathit{\boldsymbol{x}},t} \right)} \right] - {H_t}{H_x}{H_y}{H_z}\left[ {u\left( {\mathit{\boldsymbol{x}},t} \right)} \right]} \right\} \end{array} $ (B12)

右上前行波可以表达为:

$ \begin{array}{l} {u_{{\rm{urf}}}}\left( {\mathit{\boldsymbol{x}},t} \right) = \\ \frac{1}{8}\left\{ {u\left( {\mathit{\boldsymbol{x}},t} \right) - {H_t}{H_z}\left[ {u\left( {\mathit{\boldsymbol{x}},t} \right)} \right] + {H_t}{H_x}\left[ {u\left( {\mathit{\boldsymbol{x}},t} \right)} \right] + {H_t}{H_y}\left[ {u\left( {\mathit{\boldsymbol{x}},t} \right)} \right] + {H_x}{H_z}\left[ {u\left( {\mathit{\boldsymbol{x}},t} \right)} \right] - } \right.\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left. {{H_x}{H_y}\left[ {u\left( {\mathit{\boldsymbol{x}},t} \right)} \right] + {H_y}{H_z}\left[ {u\left( {\mathit{\boldsymbol{x}},t} \right)} \right] + {H_t}{H_x}{H_y}{H_z}\left[ {u\left( {\mathit{\boldsymbol{x}},t} \right)} \right]} \right\} \end{array} $ (B13)

左下后行波可以表达为:

$ \begin{array}{l} {u_{{\rm{dlb}}}}\left( {\mathit{\boldsymbol{x}},t} \right) = \\ \frac{1}{8}\left\{ {u\left( {\mathit{\boldsymbol{x}},t} \right) + {H_t}{H_z}\left[ {u\left( {\mathit{\boldsymbol{x}},t} \right)} \right] - {H_t}{H_x}\left[ {u\left( {\mathit{\boldsymbol{x}},t} \right)} \right] - {H_t}{H_y}\left[ {u\left( {\mathit{\boldsymbol{x}},t} \right)} \right] + {H_x}{H_z}\left[ {u\left( {\mathit{\boldsymbol{x}},t} \right)} \right] - } \right.\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left. {{H_x}{H_y}\left[ {u\left( {\mathit{\boldsymbol{x}},t} \right)} \right] + {H_y}{H_z}\left[ {u\left( {\mathit{\boldsymbol{x}},t} \right)} \right] - {H_t}{H_x}{H_y}{H_z}\left[ {u\left( {\mathit{\boldsymbol{x}},t} \right)} \right]} \right\} \end{array} $ (B14)

左下前行波可以表达为:

$ \begin{array}{l} {u_{{\rm{dlf}}}}\left( {\mathit{\boldsymbol{x}},t} \right) = \\ \frac{1}{8}\left\{ {u\left( {\mathit{\boldsymbol{x}},t} \right) + {H_t}{H_z}\left[ {u\left( {\mathit{\boldsymbol{x}},t} \right)} \right] - {H_t}{H_x}\left[ {u\left( {\mathit{\boldsymbol{x}},t} \right)} \right] + {H_t}{H_y}\left[ {u\left( {\mathit{\boldsymbol{x}},t} \right)} \right] + {H_x}{H_z}\left[ {u\left( {\mathit{\boldsymbol{x}},t} \right)} \right] + } \right.\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left. {{H_x}{H_y}\left[ {u\left( {\mathit{\boldsymbol{x}},t} \right)} \right] - {H_y}{H_z}\left[ {u\left( {\mathit{\boldsymbol{x}},t} \right)} \right] + {H_t}{H_x}{H_y}{H_z}\left[ {u\left( {\mathit{\boldsymbol{x}},t} \right)} \right]} \right\} \end{array} $ (B15)

右下后行波可以表达为:

$ \begin{array}{l} {u_{{\rm{drb}}}}\left( {\mathit{\boldsymbol{x}},t} \right) = \\ \frac{1}{8}\left\{ {u\left( {\mathit{\boldsymbol{x}},t} \right) + {H_t}{H_z}\left[ {u\left( {\mathit{\boldsymbol{x}},t} \right)} \right] + {H_t}{H_x}\left[ {u\left( {\mathit{\boldsymbol{x}},t} \right)} \right] - {H_t}{H_y}\left[ {u\left( {\mathit{\boldsymbol{x}},t} \right)} \right] - {H_x}{H_z}\left[ {u\left( {\mathit{\boldsymbol{x}},t} \right)} \right] + } \right.\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left. {{H_x}{H_y}\left[ {u\left( {\mathit{\boldsymbol{x}},t} \right)} \right] + {H_y}{H_z}\left[ {u\left( {\mathit{\boldsymbol{x}},t} \right)} \right] + {H_t}{H_x}{H_y}{H_z}\left[ {u\left( {\mathit{\boldsymbol{x}},t} \right)} \right]} \right\} \end{array} $ (B16)

右下前行波可以表达为:

$ \begin{array}{l} {u_{drf}}\left( {\mathit{\boldsymbol{x}},t} \right) = \\ \frac{1}{8}\left\{ {u\left( {\mathit{\boldsymbol{x}},t} \right) + {H_t}{H_z}\left[ {u\left( {\mathit{\boldsymbol{x}},t} \right)} \right] + {H_t}{H_x}\left[ {u\left( {\mathit{\boldsymbol{x}},t} \right)} \right] + {H_t}{H_y}\left[ {u\left( {\mathit{\boldsymbol{x}},t} \right)} \right] - {H_x}{H_z}\left[ {u\left( {\mathit{\boldsymbol{x}},t} \right)} \right] - } \right.\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left. {{H_x}{H_y}\left[ {u\left( {\mathit{\boldsymbol{x}},t} \right)} \right] - {H_y}{H_z}\left[ {u\left( {\mathit{\boldsymbol{x}},t} \right)} \right] - {H_t}{H_x}{H_y}{H_z}\left[ {u\left( {\mathit{\boldsymbol{x}},t} \right)} \right]} \right\} \end{array} $ (B17)
致谢: 感谢Hongbo Zhou和Jun Mu的帮助和讨论。感谢MIT的FFTW, SMAART JV的Sigsbee2A模型和数据, 以及SEG的SEAM模型和数据。
参考文献
[1]
BAYSAL E, KOSLOFF D D, SHERWOOD J W C. Reverse time migration[J]. Geophysics, 1983, 48(11): 1514-1524. DOI:10.1190/1.1441434
[2]
MCMECHAN G A. Migration by extrapolation of time-dependent boundary values[J]. Geophysics Prospecting, 1983, 31(3): 413-420. DOI:10.1111/gpr.1983.31.issue-3
[3]
ZHOU H, ZHANG G, BLOOR R. An anisotropic acoustic wave equation for modeling and migration in 2D TTI media[J]. Expanded Abstracts of 76th Annual Internat SEG Mtg, 2006, 194-198.
[4]
ZHANG Y, SUN J. Practical issues of reverse time migration:true amplitude gathers, noise removal and harmonic source encoding[J]. First Break, 2009, 26(1): 29-35.
[5]
CLAERBOUT J F. Toward a unified theory of reflector mapping[J]. Geophysics, 1971, 36(3): 467-481. DOI:10.1190/1.1440185
[6]
MILLER D, ORISTAGLIO M, BEYLKIN G. A new slant on seismic imaging:migration and integral geometry[J]. Geophysics, 1987, 52(7): 943-964. DOI:10.1190/1.1442364
[7]
LAMBARÉ G, VIRIEUX J, MADARIAGA R, et al. Iterative asymptotic inversion in the acoustic approximation[J]. Geophysics, 1992, 57(9): 1138-1154. DOI:10.1190/1.1443328
[8]
THIERRY P, LAMBARÉ G, PODVIN P, et al. 3-D preserved amplitude prestack depth migration on a workstation[J]. Geophysics, 1999, 64(1): 222-229. DOI:10.1190/1.1444518
[9]
ZHANG Y, ZHANG G, BLEISTEIN N. True amplitude wave equation migration arising from true amplitude one-way wave equations[J]. Inverse Problems, 2003, 19(5): 1113-1138. DOI:10.1088/0266-5611/19/5/307
[10]
DENG F, MCMECHAN G A. 3D true-amplitude prestack depth migration[J]. Expanded Abstracts of 77th Annual Internat SEG Mtg, 2007, 2160-2164.
[11]
XU S, ZHANG Y, TANG B. 3D angle gathers from reverse time migration[J]. Geophysics, 2011, 76(2): S77-S92. DOI:10.1190/1.3536527
[12]
GIBOLI M, BAINA R, DUQUET B. Reverse time migration surface offset gathers part 1:a new method to produce 'classical' common image gathers[J]. Expanded Abstracts of 82nd Annual Internat SEG Mtg, 2012, 3334-3338.
[13]
DE BRUIN C G M, WAPENAAR C P A, BERKHOUT A J. Angle-dependent reflectivity by means of prestack migration[J]. Geophysics, 1990, 55(9): 1223-1234. DOI:10.1190/1.1442938
[14]
XU S, CHAURIS H, LAMBARÉ G, et al. Common angle image gather:a strategy for imaging complex media[J]. Expanded Abstracts of 68th Annual Internat SEG Mtg, 1998, 1538-1541.
[15]
XU S, CHAURIS H, LAMBAR G, et al. Common angle migration:a strategy for imaging complex media[J]. Geophysics, 2001, 66(6): 1877-1894. DOI:10.1190/1.1487131
[16]
SAVA P, FOMEL S. Angle-domain common-image gathers by wavefield continuation methods[J]. Geophysics, 2003, 68(3): 1065-1074. DOI:10.1190/1.1581078
[17]
SAVA P, FOMEL S. Time-shift imaging condition in seismic migration[J]. Geophysics, 2006, 71(6): S209-S217. DOI:10.1190/1.2338824
[18]
BIONDI B, SYMES W W. Angle-domain common-image gathers for migration velocity analysis by wavefield-continuation imaging[J]. Geophysics, 2004, 69(5): 1283-1298. DOI:10.1190/1.1801945
[19]
SHAN G, BIONDI B. Angle-domain common-image gathers for steep reflectors[J]. Expanded Abstracts of 78th Annual Internat SEG Mtg, 2008, 3068-3072.
[20]
ZHANG Q, MCMECHAN G A. Direct vector-field method to obtain angle-domain common-image gathers from isotropic acoustic and elastic reverse time migration[J]. Geophysics, 2011, 76(5): WB135-WB149. DOI:10.1190/geo2010-0314.1
[21]
VYAS M, NICHOLS D, MOBLEY E. Efficient RTM angle gathers using source directions[J]. Expanded Abstracts of 81st Annual Internat SEG Mtg, 2011, 3104-3108.
[22]
YOON K, GUO M, CAI J, WANG B. 3D RTM angle gathers from source wave propagation direction and dip of reflector[J]. Expanded Abstracts of 81st Annual Internat SEG Mtg, 2011, 3136-3140.
[23]
YOON K, MARFURT K J. Reverse-time migration using the Poynting vector[J]. Exploration Geophysics, 2006, 37(1): 102-107. DOI:10.1071/EG06102
[24]
LIU F, ZHANG G, MORTON S A, et al. An effective imaging condition for reverse-time migration using wavefield decomposition[J]. Geophysics, 2011, 76(1): S29-S39. DOI:10.1190/1.3533914
[25]
COOLEY J W, TUKEY J W. An algorithm for the machine calculation of complex Fourier series[J]. Mathematics of Computation, 1965, 19(90): 297-301. DOI:10.1090/S0025-5718-1965-0178586-1
[26]
HARDY G H, LITTLEWOOD J E, PÓLYA G. Inequalities[M]. Cambridge: Cambridge University Press, 1952: 1-336.
[27]
SHEN P, ALBERTIN U. Up-down separation using Hilbert transformed source for causal imaging condition[J]. Expanded Abstracts of 85th Annual Internat SEG Mtg, 2015, 4175-4179.
[28]
FEI T W, LUO Y, YANG J, LIU H, QIN F. Removing false images in reverse time migration:the concept of de-primary[J]. Geophysics, 2015, 80(6): S237-S244. DOI:10.1190/geo2015-0289.1
[29]
AKI K, RICHARDS P G. Quantitative seismology:theory and methods[M]. San Franscisco: W H Freeman and Co, 1980: 1-700.
[30]
YANG D, PENG J, LIU M, TERLAKY T. Optimal nearly analytic discrete approximation to the scalar wave equation[J]. Bulletin of the Seismological Society of America, 2006, 96(3): 1114-1130. DOI:10.1785/0120050080
[31]
ZHOU H, ZHANG G. Prefactored optimized compact finite-difference schemes for second spatial derivatives[J]. Geophysics, 2011, 76(5): WB87-WB95. DOI:10.1190/geo2011-0048.1
[32]
CHEN F, XU S. Pyramid-shaped grid for elastic wave propagation[J]. Expanded Abstracts of 82nd Annual Internat SEG Mtg, 2012, 1-5.
[33]
XU S, ZHOU H. Accurate simulations of pure quasi-P-waves in complex anisotropic media[J]. Geophysics, 2014, 79(6): T341-T348. DOI:10.1190/geo2014-0242.1
[34]
WANG M, XU S. Finite-difference time dispersion transforms for wave propagation[J]. Geophysics, 2015, 80(6): WD19-WD25. DOI:10.1190/geo2015-0059.1
[35]
ZHANG Y, ZHANG G. One-step extrapolation method for reverse time migration[J]. Geophysics, 2009, 74(4): A29-A33. DOI:10.1190/1.3123476
[36]
COSTA J C, SILVA NETO F A, ALCÂNTARA M R M, et al. Obliquity-correction imaging condition for reverse time migration[J]. Geophysics, 2009, 74(3): S57-S66. DOI:10.1190/1.3110589
[37]
JIN S, MADARIAGA R, VIRIEUX J, et al. Two-dimensional asymptotic iterative elastic inversion[J]. Geophysical Journal International, 1992, 108(2): 575-588. DOI:10.1111/gji.1992.108.issue-2
[38]
OPERTO M S, XU S, LAMBARÉ G. Can we quantitatively image complex structures with rays?[J]. Geophysics, 2000, 65(4): 1223-1238. DOI:10.1190/1.1444814
[39]
ZHANG Y, XU S, TANG B, et al. Angle gathers from reverse time migration[J]. The Leading Edge, 2010, 29(11): 1364-1371. DOI:10.1190/1.3517308
[40]
TANG B, XU S, ZHANG Y. 3D angle gathers with plane-wave reverse-time migration[J]. Geophysics, 2013, 78(2): S117-S123. DOI:10.1190/geo2012-0313.1
[41]
KOSLOFF D D, BAYSAL E. Forward modeling by a Fourier method[J]. Geophysics, 1982, 47(10): 1402-1412. DOI:10.1190/1.1441288
[42]
FRIGO M, JOHNSON S G. The design and implementation of FFTW3[J]. Proceedings of the IEEE, 2005, 93(2): 216-231. DOI:10.1109/JPROC.2004.840301
[43]
STORK C, COMPTON S, HEUERMANN P. RTM images from SEAM data show interesting features[J]. Expanded Abstracts of 81st Annual Internat SEG Mtg, 2011, 3196-3200.