0 引言

近年来，随着计算机性能和存储能力大幅提高，逆时偏移技术迅速发展成熟并得到广泛应用^[1-7]。目前最常见逆时偏移技术采用显式有限差分方法进行正反波场传播，利用互相关成像条件获取成像值，进而通过滤波方式去除低频噪声^[8-9]。这类方式在应用中存在两个局限：一是实际偏移时计算网格较大，显式有限差分法存在较严重的数值频散^[10-12]，且产生无规则噪声，不利于高频信息的保护与成像；二是成像后滤波在去除低频噪声的同时，会损伤有效低频信息，不利于低频信息的保护与成像。

空间差分数值频散由地震波的传播方向、空间差分精度和空间采样密度等三个因素决定。其中:因地震波向四周传播，故传播方向无法控制；提高空间差分精度，会使数值频散减弱；空间采样越密，一个波长范围内离散点数越多，数值频散越弱。因此，要想减弱有限差分算子的数值频散，通常可采用更高的差分阶数或减小空间网格尺寸^[13-14]，但这将导致计算量的显著增加。

伪谱法^[15-20]利用快速傅里叶变换将波场数据由空间域变换到波数域，避免了空间差分运算，对空间微分算子的逼近使其计算精度可达Nyquist谱精度，几乎没有空间数值频散。但因其时间方向二阶导数是采用有限差分算子近似，故仍存在时间频散，当时间采样间隔较大时，会出现不稳定现象。Etgen等^[21]提出拟解析法，在波数—空间域直接修正拉普拉斯算子，补偿因时间二阶差分近似而造成的解析误差，相对于伪谱法而言可有效消除时间方向数值频散。对于速度变化剧烈的介质，可由若干常速拟拉普拉斯算子插值得到变速拟拉普拉斯算子。为了更好处理变速介质情况下的波场传播问题，Song等^[22]提出傅里叶有限差分法，将速度分解为背景速度和扰动速度，分别利用快速傅里叶变换和有限差分进行波场传播计算，实现了对波场传播数值频散的有效压制。

逆时偏移互相关成像条件包含炮点与检波点传播方向相同波场的成像计算，而这部分成像与Claer-bout成像准则不相符，它是逆时偏移传播路径上低频、强振幅偏移噪声产生的根本原因^[4]。常规叠后滤波法虽简单快捷，但会损失成像剖面中的低频有效信息，影响高陡断面构造成像质量。

为了消除逆时偏移的低频噪声，Fei等^[23-24]提出在频率—波数域将波场分解为不同方向的传播波场，再做选择性互相关成像，从而实现对低频噪声和偏移假象的压制。然而，这种频率—波数域波场分解方式需保存所有时刻波场，并沿最慢维的时间方向做傅里叶变换，计算量和I/O量都很大，实际应用中难以实现。对此，Liu等^[25]提出仅在z方向波数域实现波场分解的隐式分解法，避免了时间域傅里叶变换导致的计算和存储问题，有利于逆时偏移中消除低频噪声。胡江涛等^[26]、Shen等^[27]提出基于解析时间波场的显式波场分解逆时偏移方法，消除低频噪声的同时压制了偏移构造假象，波场外推过程中需同时求解两次波动方程，利用原波场与构建的解析时间波场进行波场分解，避免了频率—波数域波场分解的存储问题。此类先波场分解再相关成像的逆时偏移方法，尽管成像条件与单程波偏移有些类似，不利于盐下构造回折波成像，但能满足一般复杂构造成像需求，因而受到广泛关注。

本文针对常规逆时偏移显式有限差分波场传播算法和拉普拉斯滤波方法存在的高低频信息损失问题，采用傅里叶有限差分算法计算波场传播，波场外推中直接在时间—波数域进行显式波场分解和互相关成像，从而实现对高、低频信息的充分保护与利用；针对傅里叶有限差分计算效率较低问题，研讨了依托GPU的算法，显著提高了傅里叶有限差分逆时偏移的计算速度。将该方法应用于实际地震资料处理，其成像效果优于常规逆时偏移。

1 方法原理 1.1 傅里叶有限差分波场传播算法

逆时偏移双程波方程为

$ \frac{{{\partial ^2}u\left( {\mathit{\boldsymbol{x}}, t} \right)}}{{\partial {t^2}}} = {v^2}\left( \mathit{\boldsymbol{x}} \right){\nabla ^2}u\left( {\mathit{\boldsymbol{x}}, t} \right) $

(1)

式中：u(x，t)是空间点x在时刻t的地震波场；v是地震波在介质中的传播速度。

据Soubaras等^[28]提出的时间匹配算法和傅里叶变换公式，假设地下介质为常速，则式(1)可写为

$ \begin{array}{l} u\left( {\mathit{\boldsymbol{x}}, t + \Delta t} \right) + u\left( {\mathit{\boldsymbol{x}}, t - \Delta t} \right) - 2u\left( {\mathit{\boldsymbol{x}}, t} \right)\\ = 2\int_{ - \infty }^{ + \infty } {\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over u} } \left( {\mathit{\boldsymbol{k}}, t} \right)\left[ {\cos \left( {\left| \mathit{\boldsymbol{k}} \right|{v_c}\Delta t} \right) - 1} \right]{{\rm{e}}^{ - {\rm{i}}\mathit{\boldsymbol{kx}}}}{\rm{d}}\mathit{\boldsymbol{k}} \end{array} $

(2)

式中：k是波数矢量；${\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over u} }$(k，t)表示时间—波数域的地震波场；Δt是时间步长；v_c是假设的常速度值。

实际地下介质是变速的，可用v(x)近似代替式(2)中的v_c，并引入背景速度v₀，则有

$ \begin{array}{l} \cos \left[ {\left| \mathit{\boldsymbol{k}} \right|v\left( \mathit{\boldsymbol{x}} \right)\Delta t} \right] - 1 = \left[ {\cos \left( {\left| \mathit{\boldsymbol{k}} \right|{v_0}\Delta t} \right) - 1} \right] \times \\ \frac{{\cos \left[ {\left| \mathit{\boldsymbol{k}} \right|v\left( \mathit{\boldsymbol{x}} \right)\Delta t} \right] - 1}}{{\cos \left( {\left| \mathit{\boldsymbol{k}} \right|{v_0}\Delta t} \right) - 1}} \end{array} $

(3)

上式右端的分式可展开为

$ \frac{{\cos \left[ {\left| \mathit{\boldsymbol{k}} \right|v\left( \mathit{\boldsymbol{x}} \right)\Delta t} \right] - 1}}{{\cos \left( {\left| \mathit{\boldsymbol{k}} \right|{v_0}\Delta t} \right) - 1}} \approx a + 2\sum\limits_{n = 1}^3 {{b_n}\cos \left( {{k_n}\Delta {x_n}} \right)} $

(4)

式中：n=1, 2, 3代表三个不同的方向；Δx_n、k_n分别是某一方向的网格间距和波数；a、b_n为系数，可分别表示为

$ a = \frac{{{v^2}\left( \mathit{\boldsymbol{x}} \right)}}{{v_0^2}}\left\{ {1 + \frac{{{{\left( {\Delta t} \right)}^2}\left[ {v_0^2 - {v^2}\left( \mathit{\boldsymbol{x}} \right)} \right]}}{6} \times \frac{1}{{{{\left( {\Delta {x_1}} \right)}^2} + {{\left( {\Delta {x_2}} \right)}^2} + {{\left( {\Delta {x_3}} \right)}^2}}}} \right\} $

(5)

$ {b_n}{\rm{ = }}\frac{{{{\left( {\Delta t} \right)}^2}{v^2}\left( \mathit{\boldsymbol{x}} \right)\left[ {{v^2}\left( \mathit{\boldsymbol{x}} \right) - v_0^2} \right]}}{{12v_0^2{{\left( {\Delta {x_n}} \right)}^2}}} $

(6)

式(2)等号右端可分为两步实现：首先基于背景速度利用傅里叶正反变换进行波场传播；然后再利用有限差分算法校正速度扰动影响。因而，双程波傅里叶有限差分波场传播算法实现步骤为：

(1) 利用FFT将时间—空间域波场u(x，t)变换成时间—波数域波场${\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over u} }$(k，t)；

(2) 将${\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over u} }$(k，t)乘以cos(|k|v₀Δt)-1)得到以背景速度传播的时间—波数域波场${\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over q} }$(k，t)；

(3) 利用反FFT将${\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over q} }$ (k，t)变换成其对应的时间—空间域波场q(x，t)；

(4) 利用式(5)、式(6)中差分系数，通过有限差分法由q(x，t)计算速度扰动下的传播波场q(x，t+Δt)；

(5) 基于时间差分格式，新时刻的时间—空间域波场u(x，t+Δt)表示为q(x，t+Δt)+2u(x，t)-u(x，t-Δt)；

(6) 重复步骤(1)~(5)，逐个计算各时刻波场。

利用上述步骤开展双程波波场传播试验。首先给定一个三维均匀介质模型，速度为3000m/s，模型尺寸为9000m×9000m×5000m，速度网格单元为20m×20m×20m，震源点位于模型中心(4500m，4500m，2500m)，地震子波是主频为40Hz的雷克子波，时间步长为1ms。

图 1是波场快照剖面，图 2是其水平切片。可见常规有限差分法波场传播过程中存在较为严重的频散噪声，而傅里叶有限差分法波场快照中，频散得到较彻底压制，有利于高频信息的保护。

进一步利用三维盐丘模型检测傅里叶有限差分法针对复杂介质的处理效果。速度网格单元为30m×30m×20m，考虑到海水速度为1500m/s，地震子波采用主频为20Hz雷克子波，时间步长为2ms。图 3是盐丘速度模型及对应的T1时刻常规有限差分法和傅里叶有限差分法波场快照，图 4是应用这两种方法得到的T2时刻波场快照的水平切片。对比亦能看出，常规有限差分法波场频散较严重，而傅里叶有限差分法波场信息更保真。

1.2 时间—波数域解析波场分解方法

解析波场表示为

$ q\left( {\mathit{\boldsymbol{x}}, t} \right) = u\left( {\mathit{\boldsymbol{x}}, t} \right) + {\rm{i}}p\left( {\mathit{\boldsymbol{x}}, t} \right) $

(7)

其实部是地震传播波场u(x, t)，虚部p(x, t)是波场u(x, t)的Hilbert变换。Hilbert变换是一个时间方向的卷积，直接进行Hilbert变换计算需将所有时刻的波场数据存盘并做卷积运算，数据量和计算量都很大，因而实用性较差。胡江涛等^[26]提出对波动方程的源项进行Hilbert变换的方法。考虑到波场传播算子的线性特征，传播波场即是解析时间波场的虚部，如此通过两个波场传播构建各时刻的解析波场。另外，傅里叶有限差分算法在波数域和空间域是交替进行波场传播，本文采用Zhang等^[29]给出的波数域Hilbert变换方法

$ \mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over p} \left( {\mathit{\boldsymbol{k}}, t} \right) = \frac{1}{{v\sqrt {k_x^2 + k_y^2 + k_z^2} }}\frac{{\partial \mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over u} \left( {\mathit{\boldsymbol{k}}, t} \right)}}{{\partial t}} $

(8)

式中：$\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over p} $(k, t)是${\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over u} }$(k, t)的Hilbert变换波场；k_x、k_y、k_z是三个方向的波数值。构造解析波场后，上、下行波分别表示为

$ {u_{\rm{U}}}\left( {\mathit{\boldsymbol{x}}, t} \right) = {\mathop{\rm Re}\nolimits} \left\{ {{{\rm{F}}^{ - 1}}\left[ {{u^ + }\left( {\mathit{\boldsymbol{k}}, t} \right)} \right]} \right\} $

(9)

$ {u_{\rm{D}}}\left( {\mathit{\boldsymbol{x}}, t} \right) = {\mathop{\rm Re}\nolimits} \left\{ {{{\rm{F}}^{ - 1}}\left[ {{u^ - }\left( {\mathit{\boldsymbol{k}}, t} \right)} \right]} \right\} $

(10)

式中：F^-1表示三维空间傅里叶反变换；Re表示对复数取实部。其传播波场为

$ {u^ + }\left( {\mathit{\boldsymbol{k}}, t} \right) = \left\{ {\begin{array}{*{20}{c}} {\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over q} \left( {\mathit{\boldsymbol{k}}, t} \right)\;\;\;{k_z} \ge 0}\\ {0\;\;\;\;\;\;\;\;\;\;\;{k_z} < 0} \end{array}} \right. $

(11)

$ {u^ - }\left( {\mathit{\boldsymbol{k}}, t} \right) = \left\{ {\begin{array}{*{20}{c}} {0\;\;\;\;\;\;\;\;\;\;\;{k_z} \ge 0}\\ {\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over q} \left( {\mathit{\boldsymbol{k}}, t} \right)\;\;\;{k_z} < 0} \end{array}} \right. $

(12)

式中${\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over q} }$(k, t)是时间—波数域的解析波场。

图 5是图 1b所示数据分离后的波场，可见上、下行波得到很好分离。此外，对三维盐丘模型的复杂波场也进行了波场分离试验，图 6a、图 6b对应为图 3c波场分离后的上行波和下行波，分离效果整体较好，但也存在一定残余。考虑到实际逆时偏移采用的速度模型通常是低频光滑模型，透射波能量明显大于反射波能量，且逆时偏移是震源顺时传播的下行波(透射波)与检波器逆时传播的下行波(透射波)相关成像，其上行波残余可忽略不计。

1.3 GPU加速方法

叠前偏移成像计算量大、效率低，通常采用GPU并行算法实现加速^[30-32]。傅里叶有限差分波场传播算法需做大量的正反三维傅里叶变换，计算效率低于常规高阶有限差分算法。本文利用GPU运算平台CUDA(Compute Unified Device Architecture)提供的用于傅里叶变换的函数库cuFFT(CUDA Fast Fourier Transform)进行加速。试算结果表明，cuFFT对信号长度较敏感，当信号长度为奇数时，运算速度较慢；当信号长度为偶数时，速度较快；当信号长度为64的整数倍时，速度最快。因此，执行cuFFT之前需先对三维空间网格进行镶边，使x、y、z方向空间网格数为64的整数倍，达到实现快速傅里叶变换的目的。

GPU加速傅里叶有限差分法有两种计算方案：一是仅将正反三维FFT置于GPU中，有限差分计算由CPU完成，该方案对内存需求小，但每一步波场传播都需做GPU数据I/O，影响计算效率；二是将整个计算步骤都置于GPU中，避免每步波场传播的数据I/O，显然会增加内存需求，特别是FD有限差分系数需占用大量内存空间。对此，在每步波场传播过程中重复计算FD有限差分系数，通过增加计算量换取减小内存空间，以降低内存需求。

同样采用上面的参数和硬件资源条件测试不同方法的计算效率，得到的对比结果如表 1。可见上述方案二效率提升较大，达到12倍加速比。GPU内存耗费方面，方案二需存储速度场、成像值和两个时刻的复数域波场，大致是常规有限差分法的1.5倍，目前常规GPU设备就能满足此内存需求。本文采用方案二进行GPU加速计算。

2 应用实例

M探区构造较简单，资料品质高，常规RTM成像分辨率较低(图 7a)。应用本文方法进行偏移处理得到成像剖面(图 7b)，可见本文方法成像结果分辨率更高，细节信息更丰富。图 8是两个结果由深度域转时间域后的频谱图，显然可见本文方法成像频带更宽。

N探区构造较复杂，资料品质一般。从常规RTM法(图 9a)和本文方法(图 9b)成像剖面看出，本文方法分辨率依然较高；从其成像频谱对比图(图 10)可见，本文方法具有更宽频带。

3 结论

针对常规逆时偏移中高阶有限差分算法、叠后滤波算法中的不足，本文研究了基于GPU的傅里叶有限差分波场传播算法和时间—波数域显式波场分解方法，形成了基于GPU的傅里叶有限差分逆时偏移成像技术，得到以下结论和认识。

(1) 逆时偏移的高频信息损失缘于高阶有限差分空间频散，傅里叶有限差分波场传播方法将速度分解为背景速度和扰动速度，分别利用快速傅里叶变换和有限差分计算波场传播，能较好压制数值频散，更好地保护高频数据信息。

(2) 基于解析波场的上、下行波分解是逆时偏移压制低频噪声的重要手段，在傅里叶有限差分波场传播中，直接对时间—波数域波场做波场分解，实现简单，额外计算量小。

(3) 傅里叶有限差分波场传播算法中包含大量的FFT计算，效率较低，通过GPU平台cuFFT函数加速并选择合适的I/O策略，能显著提高傅里叶有限差分法计算效率，使其实用性增强。

(4) 实际资料应用结果表明，该技术有利于高低频信息的充分保护与利用，提高了地震成像分辨率，成像效果优于常规逆时偏移。