应用GPU 的傅里叶有限差分逆时偏移

段 心 标

（①同济大学海洋与地球科学学院波现象与智能反演成像研究组，上海200092；②中国石化石油物探技术研究院，江苏南京211103）

0 引言

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

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

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

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

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

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

1 方法原理

1.1 傅里叶有限差分波场传播算法

逆时偏移双程波方程为

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

据Soubaras等［28］提出的时间匹配算法和傅里叶变换公式，假设地下介质为常速，则式（1）可写为式中：k是波数矢量；︵u（k，t）表示时间—波数域的地震波场；Δt是时间步长；vc是假设的常速度值。

实际地下介质是变速的，可用v（x）近似代替式（2）中的vc，并引入背景速度v0，则有

上式右端的分式可展开为

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

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

（1）利用FFT 将时间—空间域波场u（x，t）变换成时间—波数域波场（k，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是其水平切片。可见常规有限差分法波场传播过程中存在较为严重的频散噪声，而傅里叶有限差分法波场快照中，频散得到较彻底压制，有利于高频信息的保护。

图1 均匀速度模型的常规有限差分法（a）和傅里叶有限差分法（b）波场快照剖面

图2 均匀速度模型的常规有限差分法（a）和傅里叶有限差分法（b）波场快照水平切片

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

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

解析波场表示为

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

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

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

1.3 GPU 加速方法

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

图3 盐丘模型及其T1时刻波场快照剖面

图4 盐丘模型的常规有限差分法（a）和傅里叶有限差分法（b）T2时刻波场快照水平切片

图5 均匀速度模型的下行波（a）和上行波（b）波场分离结果

图6 盐丘模型的上行波（a）和下行波（b）波场分离结果

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

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

表1 不同算法效率对比

2 应用实例

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

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

图7 M 探区资料常规RTM（a）和本文方法（b）成像剖面

图8 M 探区资料常规RTM （蓝线）和本文方法（红线）成像剖面的频谱

3 结论

图9 N 探区资料常规RTM（a）和本文方法（b）成像剖面

图10 N 探区资料常规RTM （蓝线）和本文方法（红线）成像剖面的频谱

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

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

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

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

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