基于改进匹配追踪傅里叶插值的地震数据规则化重构

张春枫 杜启振，2 张富源 符力耘，2 魏国华

（1. 中国石油大学（华东）山东省深层油气重点实验室，山东 青岛 266580；2. 青岛海洋科学与技术试点国家实验室海洋矿产资源评价与探测技术功能实验室，山东 青岛 266237；3. 中国石化胜利油田分公司物探研究院，山东 东营 257022）

0 引 言

随着中国油气资源勘探开发逐步深化，对野外采集的地震资料质量提出了更高的要求。理论上常常假定，地震数据沿着空间方向上按等距离均匀采样。在野外地震数据采集过程中，地震资料在时间方向上可以做到均匀、密集采样，受制于实际施工环境等因素的影响，在空间方向上的采样往往是不规则甚至是稀疏的，这会对后续的地震数据处理流程（地震数据去噪、CMP 道集叠加、偏移成像等）产生不良影响［1］。为使得基于规则采样假设条件的常规地震数据处理流程得以顺利进行，不规则地震数据的插值与规则化重构具有十分重要的意义。

地震数据规则化的方法主要分为以下6 类：第1 类是采用线性预测滤波的方法进行重构［2-5］；第2类是以波动方程为基础的重构方法［6］；第3 类是矩阵降秩的重构方法［7-8］；第4 类是相干倾角的重构方法［9］；第5 类是基于数学变换域的重构方法；第6 类是基于人工智能的重构方法［10-11］。目前，基于傅里叶变换域的重构方法受到颇多关注，其优点在于不受空间采样影响，可以很好地拓展到高维领域。A.J.W.Duijndam 等［12］将傅里叶变换应用于二维不规则地震数据的恢复与重构；K.Hindriks 等［13］又将其拓展到三维领域；B.Liu 等［14］采用最小加权范数对不规则地震数据进行插值重构；S.Xu 等［15-16］结合了非均匀离散傅里叶变换，提出了基于抗泄露傅里叶变换（Anti-leakage Fourier Transform，ALFT）的不规则地震数据重构，并将其应用于三维地震数据重构；在S.Xu 等研究的基础上，ALFT 算法中傅里叶系数计算部分也被改进与优化［17-19］，包括引入非均匀快速傅里叶变换、GPU 加速算法等。傅里叶变换方法的重构效率受傅里叶系数计算时长的影响，加快傅里叶系数的计算速度是提升不规则地震数据重构效率的关键之一。

为解决傅里叶系数计算耗时高这一问题，本文采用与ALFT 原理相近的匹配追踪傅里叶插值（Matching Pursuit Fourier Interpolation，MPFI） 算法，在MPFI 算法的基础上对其迭代结构进行改进，在不规则地震数据重构前计算其傅里叶系数矩阵。和传统MPFI 算法相比，改进后的算法优化了循环中大量复杂的计算，减少迭代所需时间，从而提高不规则地震数据的重构效率。

1 地震数据重构原理

采用傅里叶变换将不规则数据规则化符合信号稀疏分解理论。信号稀疏分解是指具有稀疏性的信号可以在重构数据最少、重构误差最小的情况下，由有限个最优重构数据的线性加权组合表示，即使用较少的数据重构出与原信号最接近或者较为接近的解。对于规则缺失的地震数据，傅里叶变换重构十分有效；但是对于不规则采样地震数据，其傅里叶基函数不再互相正交，导致各个频率之间发生频谱泄露，即某个傅里叶系数的能量泄漏到其他的傅里叶系数上，能量最强的频率造成的泄漏最大。此时，常规的傅里叶重构方法不再适用，亟待发展新的方法。S.G.Mallat 等［20］提出了一种基于稀疏计算的匹配追踪算法（Matching Pursuit，MP）来解决这一问题。任何信号的波形都可以通过一系列正弦波累加得到，可用傅里叶级数表示为

式中：f(x)——信号；an(n=0，±1，±2，…)——权重系数；i——虚数单位；ω——角频率；x——空间坐标，m。

地震数据规则化与匹配追踪算法所解决的问题类似，采用匹配追踪与傅里叶变换相结合的方式对地震数据作规则化与重构，称其为匹配追踪傅里叶插值。MPFI 算法以傅里叶变换为基础，其基本思想是先计算不规则数据的傅里叶系数，将其保存到系数谱中，循环多次后再把重构的傅里叶系数谱经过反傅里叶变换，输出到期望的空间位置上。

实现MPFI 算法的前提条件是不规则数据具有稀疏性。实际上野外采集得到的t-x域地震数据并不具有稀疏性。野外采集地震数据在时间方向上是规则采样，采样间隔相同，采用FFT 或者DFT就能将其从时间域变换到频率域；受各种因素的影响，空间方向上的检波点排布通常是不规则的，需要在地震数据的空间方向上作非均匀离散傅里叶变换，这样才能在f-k域中通过迭代计算、更新傅里叶系数谱。变换后，地震数据中大多数波数接近于0，满足不规则数据的稀疏性条件。而且变换后傅里叶系数谱上的能量更加集中，有利于后续傅里叶系数的计算与选取。将最终选取的一系列傅里叶系数组成新的系数谱，再将其经过反傅里叶变换后重构成规则采样数据。

常规的MPFI 算法针对不规则数据进行多次迭代，能够将其泄漏的能量回归，从而恢复原始数据。为了压制能量泄露，计算地震数据中其他频率成分的系数之前需要对输入数据进行更新。绝对值最大的傅里叶系数对频谱泄露的影响最强。故在更新不规则数据时应从绝对值最大的频率成分开始，在原数据中减去该频率在f-x域的贡献，以此类推，逐渐递归到绝对值较小的频率。

MPFI 与ALFT 算法在计算和寻找傅里叶系数方面基本相同，T.Nguyen 等［21］开发了一种快速MPFI 算法，其中NDFT 在每个时间频率下仅计算一次，从而大大降低了MPFI 算法的计算时间。

对于变换后的f-x域地震数据f(xl)，可以用离散和的形式表示为

其逆变换为

式中：f s(xl)——第s次迭代中的f-x域地震数据；s——迭代次数；xl(l= 1，2，…，N)——检波器坐标，m；N——检波器数量；q——空间波数索引；Nf——空间频率数量；(kq)——第s次迭代中的傅里叶系数；kq——空间波数。

MPFI 算法中需要不断地对地震数据进行更新，更新量由所选取傅里叶系数对应的贡献值决定，系数越大，泄露在频谱上的能量就越多。每次迭代中最强傅里叶系数对应空间波数ks可由公式计算得到

式中ks——第s次迭代中所选取的空间波数。

初始化迭代次数s= 1，s+ 1 次迭代后，不规则地震数据的剩余量表示为

式中：f s+1(xl)——第s+ 1 次迭代中的f-x域地震数据；̂(ks)——第s次迭代中的最大傅里叶系数。

常规MPFI 算法的重构步骤［22］：

（1）在地震数据的时间方向上作FFT 变换，得到f-x域地震数据；

（2）在空间方向上，对输入的数据作非均匀离散傅里叶变换，计算其傅里叶系数谱；

（3）选择能量最强的傅里叶系数，保存在新系数谱中；

（4）将选中的系数按照采样点的实际空间位置作非均匀反傅里叶变换；

（5）原始不规则f-x域数据减去该次反傅里叶变换的结果，更新不规则数据，再次循环；

（6）不断重复循环步骤（2）到（5），直到满足迭代误差或者预定义的迭代次数时为止；

（7）将保存的傅里叶系数谱按照期望空间位置进行反傅里叶变换，得到重构后规则地震数据。

2 重构算法的改进

MPFI 算法的主要不足是数据重构计算时间长、效率低，这是因为循环中涉及到大量非均匀离散正反傅里叶变换的计算，特别是循环中的第一步，傅里叶系数的计算量巨大。循环过程中不规则数据每更新一次，其剩余量中的傅里叶系数都要重新计算一遍。如果能够在迭代前提前估算出这些傅里叶系数，那么在循环中就可以省略大量复杂运算，不规则地震数据的重构效率将会大大提高。计算傅里叶系数的方法多种多样。本文采用直接计算傅里叶系数的方式，对公式（5）两边同时作非均匀离散傅里叶逆变换，并结合公式（3）得到

整理后，可得到不规则采样地震数据的剩余量在f-k域中的更新迭代公式

式中(kq)——第s+ 1 次迭代中的傅里叶系数。

公式（7）右边指数项的计算量仍然很大，但是其中的空间波数kq、偏移距xl只与地震数据有关，而每次选中的最强傅里叶系数所对应波数ks是空间波数集合kq中的元素，因此改进后迭代公式中的复杂指数项可以提前计算得到，并且只计算一次。与迭代公式（5）相比，迭代公式（7）将不规则地震数据剩余量的计算从f-x域转换到f-k域，地震数据只在f-k域中更新迭代。此优势是能够避免迭代中不规则地震数据在f-x域与f-k域之间频繁转换，减少了循环中大量非均匀离散傅里叶变换，从而达到提高计算效率的目的。

改进后使用MPFI 算法重构地震数据的过程变为：

（1）计算地震数据的傅里叶分量矩阵；

（2）在地震数据的时间方向和空间方向分别作傅里叶变换；

（3）选择能量最强的傅里叶系数，并保存在新傅里叶系数谱中；

（4）使用改进后的地震数据剩余量迭代公式在f-k域中更新不规则地震数据；

（5）重复进行步骤（3）到（4），直到剩余量满足误差或者达到预定义的迭代次数为止；

（6）将保存的傅里叶系数谱按照规则的空间采样点位置进行反傅里叶变换，得到规则地震数据。

3 模型测试

3.1 凹陷模型

为了检验改进后MPFI 算法的可行性，采用二维凹陷速度模型作时间二阶、空间十阶声波正演模拟测试。

二维凹陷速度模型的横向采样点为800，纵向采样点为440，横向网格宽度为5 m，纵向网格宽度为5 m，上层速度为2 500 m/s，下层速度为2 800 m/s（图1）。

图1 凹陷速度模型Fig. 1 Graben velocity model

采用完全匹配层（Perfectly Matched Layer，PML）吸收边界条件，炮点位于横向2 000 m、纵向25 m 处，空间方向共160 道，道间距为25 m，时间方向共1 700 个采样点，采样间隔时间为1 ms，雷克子波的主频为25 Hz。正演模拟后得到规则采样单炮地震记录（图2（a））。

图2 凹陷模型正演与重构结果Fig. 2 Graben model forward modeling and reconstruction result

一般来说，对于单炮地震记录，直达波的能量强于地下分界面所反射的有效波的能量。重构地震数据时，如果直达波与反射波位置较为接近或者叠合在一起，应先切除强直达波，否则在选择傅里叶系数时往往选取的是直达波的信息而不是有效波的信息，会严重影响到后续的地震数据重构。对规则正演地震记录中的检波器作随机缺失（缺失数量为32），得到欠采样单炮地震记录（图2（b））。相对于规则采样数据，随机欠采样地震数据的同相轴出现了中断、不连续现象。使用MPFI 算法与改进MPFI 算法重构后，随机缺失的检波点被重构到规则的期望空间位置上，且重构后正演结果中同相轴变得光滑、连续，随机欠采样地震数据被重构为规则采样地震数据（图2（c）、（d））。

相比于规则采样数据的频波谱（图3（a）），不规则采样数据的频波谱上还出现了能量泄露，其能量峰值降低，频率泄漏到整个频波谱内（图3（b）），使用MPFI算法与改进MPFI算法重构后，频波谱中泄露的能量回归到有效频率中（图3（c）、（d））。单道地震数据对比表明（图4），使用2种算法重构后的单道地震数据与原始单道地震数据相比，误差基本相同，2 种方法重构后的单道数据十分接近，可认为2 种方法重构效果一致。上述结果表明在二维凹陷速度模型中，改进后的重构方法具有可行性。

图5 Marmousi声波速度模型Fig. 5 Marmousi acoustic velocity model

3.2 复杂模型

为了检验改进后MPFI 算法在复杂模型中的可行性，采用抽稀Marmousi 速度模型作时间二阶、空间十阶有限差分声波正演模拟测试。为了减弱直达波对反射波的干扰，将模型上方水层加厚。如图5 所示，该速度模型的横向采样点数为530，纵向采样点数为227，采用PML 吸收边界条件，其色标表示介质中声波传播速度。

正演模拟中设置的参数为：地震数据的空间方向共有106 道，横向网格宽度为5 m，纵向网格宽度为5 m，道间距为25 m，时间方向共有1 601 个采样点，时间采样间隔为1 ms，炮点坐标位于横向1 325 m、纵向25 m 处，雷克子波的主频设为25 Hz。正演模拟后得到规则采样单炮地震记录（图6（a））。将正演模拟中的检波器按照随机缺失的方式进行数据采样（缺失数量为24），以模拟实际采集得到的不规则地震数据，可以看出，不规则采样地震数据中的反射波同相轴出现了随机中断，有明显的不连续现像（图6（b）），经过重构，同相轴缺失信息得到补偿，地下波场得到恢复（图6（c）、（d））。

图6 Marmousi模型正演与重构结果Fig. 6 Marmousi model forward modeling and reconstruction results

本文采用的2 种重构方法均是依靠计算傅里叶系数谱重构不规则地震数据，改进后的MPFI 算法是在常规算法的基础上对傅里叶系数的计算方式作出优化，对迭代公式两端同时应用傅里叶变换，因此2 种方法在重构效果上基本相同，但改进后的重构方法计算效率更加高效，这是因为使用MPFI 方法重构不规则地震数据的效率还与重建过程中傅里叶变换次数有关。

对于相同的地震数据，其数据量越大、频率范围越宽，循环中所涉及的非均匀离散傅里叶变换就越多。这些非均匀离散傅里叶变换在计算机中需要使用欧拉公式计算，因此重构时间就越长。采用MPFI 算法重构凹陷不规则采样与抽稀Marmousi 不规则采样数据，在计算傅里叶系数谱时分别需要进行大量空间方向上的非均匀傅里叶变换，在所占内存相同的情况下，使用改进后的MPFI 算法在计算傅里叶系数谱的过程中不需进行空间方向上的非均匀傅里叶变换，计算量相比于优化前大大减少，而且迭代公式（7）中涉及的复杂运算在循环前就已经一次完成，不需要在每次循环中去重复计算傅里叶系数。此外2 种方法在时间方向上的快速傅里叶变换次数相同，因此地震数据重构效率得到提高。

表1 说明了程序运行的环境配置，以及改进前后所需的总运行时间，可以看出，改进后的MPFI算法的重构时间明显缩短，运行效率显著提升。

表1 改进前、后重构计算时间对比Table 1 Comparison of reconstruction calculation before and after improvement

4 应用实例

为了验证改进后MPFI 算法在实际地震资料中的应用，采用某陆上工区叠前CMP 道集记录进行测试。由于数据量过大，现截取原地震数据中某一部分。原始CMP 道集道间距为25 m，共70 道，时间方向共500 个采样点，时间采样间隔为2 ms，其波形如图7（a）所示。不规则CMP 道集中空间方向为56 道，其采样方式为随机欠采样，占原始道集采样的80%（图7（b））。使用改进后MPFI 算法对随机欠采样CMP 道集数据进行重构，重构后其同相轴变密，缺失信息得到补偿，有效信号明显得到改善（图7（c））。

图7 实际地震资料应用Fig. 7 Application of real seismic data

5 结 论

（1）改进后的重构算法在计算、更新傅里叶系数时，不依赖大量的正反非均匀离散傅里叶变换，大大减少了迭代中的计算量，重构效率相比常规MPFI 算法得到提高。

（2）使用二维凹陷速度模型与抽稀Marmousi速度模型的正演单炮地震记录对改进后的重构算法进行测试，测试结果表明改进的MPFI 算法可以将不规则采样地震数据规则化，并且在相同条件下改进后算法的重构效率更高；将改进后的算法应用到实际地震资料中，取得了良好的重构效果。