地震资料快速两步插值算法

马泽川 李 勇*② 陈力鑫 陈 杰 王鹏飞 李雪梅

（①成都理工大学地球物理学院，四川成都610059；②油气藏地质及开发工程国家重点实验室（成都理工大学），四川成都610059）

0 引言

实际地震数据在空间方向往往存在不规则的缺失道［1-2］，要想得到完整的地震数据，通常需要插值处理［3］。地震资料插值方法大致可以分为两类［4-5］：①利用线性空间预测滤波器插值［6-7］；②借助数学变换（如傅里叶变换）插值［8-9］，其中凸集投影（POCS）算法和迭代阈值（IST）算法是较为经典的算法［10］，利用地震数据变换域的稀疏性对地震数据插值，由于原理直观、实现简单，并且具有稳定性，因此发展迅速。Menke［11］认为POCS算法具有普适性，并将其引入地球物理领域。近年来，各种基于传统POCS算法的改进方法被提出且被广泛应用。刘红喜等［12］提出了边缘保持的POCS超分辨率重建方法；王本锋等［13］引入了结合POCS和Curvelet的超分辨率重建算法重建三维数字内核；Abma 等［9］将POCS算法引入地震插值领域，获得了很好的效果；王本锋等［14］提出了POCS插值重建地震数据的质量控制新准则，在保证插值精度的同时提高了计算效率。在POCS算法发展的同时，人们也对IST 算法进行研究。Herrmann 等［15］基于基追踪降噪模型，将IST 算法用于非均匀记录地震道插值；由于现有IST 算法收敛速度较慢，Beck 等［16］从一般的数学背景出发，提出了快速迭代收缩阈值（FIST）算法；Gao等［17］使用指数递减阈值策略代替线性阈值策略。近年来，业界着重对IST和POCS算法迭代慢的不足进行优化，提高了算法的计算效率和精度。

为了提高插值效率以及选择最优的插值方案，本文基于IST 算法和POCS算法的分析公式，发展了FIST 算法和快速凸集投影（FPOCS）算法。基本思想是：将前一步的插值结果与前两步的插值结果通过线性算子进行线性组合，得到迭代收缩算子；通过插值算法进行插值。同时引入质量控制新准则［14］，提高了计算效率和精度。通过理论模型筛选出最优的阈值策略以及最大迭代次数，实际资料处理结果表明，文中方法在较少迭代次数下的插值效果较好，验证了算法的高效性。

1 基本原理

1.1 地震资料插值原理

将采集到的不规则的随机地震资料记为dobs，期望得到的完整数据记为d，则

式中M 为采样算子，是由元素0 和1 组成的对角矩阵。

为了方便插值，根据压缩感知原理［18-20］，将地震数据投影到变换域，即

式中：A∈Rn×m为正交坐标系的合成算子，且是一个紧框架；χ 为变换域的表示系数。

然而式（2）是一个欠定方程，具有多解性，因此需要添加一些限制性条件求解。由于稀疏性在变换域内始终适用［21］，通常与χ 相对应。因此，求解式（2）就转化为L0范数的最小化问题

式中K＝MA。由于L0范数的最小化非凸复杂性［22］，常常使用L1范数代替L0范数，并与L0范数的稀疏性相同［18］。因此，式（3）变为

式中τ为正则化参数。由于A 为紧框架，根据A*A＝I（I为单位矩阵）以及χ＝A*Aχ＝A*d（“上角*”表示伴随），则

式（5）为地震插值分析公式，因为该式直接将地震资料d作为未知数进行分析［23］。

1.2 IST算法

有很 多 算 法 求 解 式（4）。Figueiredo 等［24］和Daubechies等［25］提出了IST 算法，一般表示为

式中：λ为正阈值参数；k 为迭代次数，N 为最大迭代次数；Tλ为阈值函数，与τ 的取值有关，当λ＝τ时，Tλ为软阈值函数，定义软阈值算子为

式中u 为单道地震数据。与式（5）相同的思路，将式（6）两边左乘A，得到IST 算法公式

式中d（k）、χ（k）分别为第k次迭代时的d与χ 的值。

1.3 POCS算法

定义

式中u（k）为第k 次迭代后原始地震道的插入数据。对于式（8），有

d（k＋1）＝A Tλ［A*u（k）］

在第k＋1次迭代时，结合式（9），有

由M（I－M）＝0，得

在式（10）中用d（k）取代u（k），得

这样，便得到POCS算法，该算法已广泛用于许多领域，如图像重建和修复等。

1.4 两步插值：FIST和FPOCS

由于IST 算法的收敛速度很慢。Beck等［16］提出了快速迭代收缩阈值算法（FIST），可以表示为

式（14）代表两步插值的过程：将前一步的插值结果x（k）与前两步的插值结果x（k－1）通过线性算子t进行线性组合，得到迭代收缩算子～x（k）；通过插值算法计算x（k＋1），提高了迭代速度。

与式（8）推导相同，式（14）两边左乘A 得到FIST 算法公式

式中：～d（k）为地震数据插值的迭代收缩算子；d（k）为第k次迭代后的地震数据。

在式（15）的基础上，按照式（13）思路，可以推出FPOCS算法公式

1.5 两步插值方法的阈值策略

正则化参数τ具重要作用，并且与每次迭代的阈值λ密切相关。因此定义软阈值

通过对比不同阈值策略的实际效果，筛选最佳的阈值策略。

（1）恒阈值

式中C 为常数。

（2）阈值线性递减

（3）阈值指数递减

（4）数据驱动阈值

且

式中：τmax和τmin分别为所选的最大和最小正则化参数；pmax和pmin分别为所定义的最大和最小百分比；序列｛|（A*dobs）i|｝∈［τmin，τmax］按其中元素振幅的下降顺序排列可以得到向量v，Nv为v 的长度，vj为v 中的第j 个元素；「（k－1）Nv／（N－1）⏋表示取不小于（k－1）Nv／（N－1）的最小整数。

以上介绍了两步插值算法的原理。由于地震资料在变换域具有稀疏性，根据压缩感知理论可知，将地震数据投影到变换域，选取阈值策略进行反复迭代，即可完成地震数据插值。

2 仿真实验和实际应用

2.1 两步插值法的终止准则

文中 介 绍 的IST、POCS、FIST 和FPOCS 等4种插值算法都属于迭代方法，这就意味着在达到所定义的最大迭代次数N 之前，算法不会终止。因此，需要建立一个终止准则适时地结束迭代。Gao等［4-5］等给出了两种终止准则，王本锋等［14］在此基础上提出了灵活、适用性更好的质量控制新准则

进一步提高了计算效率。式中η 为公差，需要在程序运行前设置。

文中使用信噪比

评估插值质量。式中^d 为重建后的地震记录。

综上所述，总结了地震资料插值流程（图1）。

图1 地震资料插值流程

2.2 算例模型

使用IST、POCS、FIST 和FPOCS等算法分别对由Seismic Lab 建立的四层地震模型（图2a）、Marmousi模型（图2b）的不完整地震数据（图2d、图2e）进行插值，筛选出最佳的阈值策略，并最终由实际地震资料（图2c、图2f）进行验证。

2.3 最优阈值策略

首先在由Seismic Lab建立的四层地震模型上使用不同的阈值策略插值，以此筛选最优的阈值策略。应用IST、POCS、FIST 和FPOCS等算法对缺失数据（图2d）插值，并绘制信噪比曲线（图3）。可见：在恒阈值时的信噪比为38.5863 d B，两步插值法在50次迭代时便完成了图像重建，较一步插值法的图像重建速度更快（图3a）；在阈值线性递减时的信噪比为35.1875dB，在迭代次数较少时图像重建效果不明显，随着迭代次数增加，重构作用越来越大，但会引入噪声，导致信噪比降低（图3b）；在阈值指数递减时的信噪比为43.4571dB（图3c）；数据驱动阈值时的信噪比为41.1818d B（图3d）。因此，阈值指数递减与数据驱动阈值的插值结果信噪比较高，数据重建效果较好。图4为使用阈值指数递减策略对图2d的插值结果及其与图2a数据的差异。由图可见，经过插值后的地震数据恢复完好，且没有引入大量噪声。

图2 实验算例

2.4 合理的最大迭代次数

利用Marmousi模型验证两步插值方法的收敛性。图5为应用阈值指数递减策略对图2e插值的信噪比曲线。由图可见：在迭代次数较小时，重建效果较差，两步插值法的收敛速度略快于一步插值法，因此两步插值法可提高收敛速度；由于观测数据的复杂性，对Marmousi地震模型的重构精度不高；在50次迭代后得到了较好的收敛解（图5c），再增加迭代次数对插值精度的改善很小（图5d），通过多次测试发现，当N∈［35，50］时可以取得满意的插值效果。图6 为使用阈值指数递减策略迭代50 次对图2e的插值结果及其与图2b数据的差异。由图可见，四种算法具有相近的插值效果，在没有引入大量噪声的情况下，地震数据恢复完好。

此外，在迭代次数较少时POCS、FPOCS 算法的信噪比始终高于IST、FIST 算法，这是由于所有算法在初始化时将d（k）置0所致。由式（12）可知，在POCS、FPOCS算法中，有

因此

图4 使用阈值指数递减策略对图2d的插值结果及其与图2a数据的差异

2.5 合理的终止准则以及实际地震数据插值效果

通过对不完整的理论模型数据（图2d、图2e）进行插值，筛选出最优的阈值策略和最大迭代次数，在此基础上验证实际地震数据（图2c、图2f）的插值效果。图7 为使用阈值指数递减策略迭代50 次对图2f的插值结果及其与图2c数据的差异，图8 为应用阈值指数递减策略对图2f插值的信噪比曲线。由图可见：①四种算法插值结果的信噪比约为14dB，均能较好地重建缺失数据（图7）。②两步插值算法相对于单步插值算法具有更快的收敛速度（图8）；与没有加入终止准则的结果（图8a）相比，应用终止准则后，算法能更快地结束迭代，从而提高计算效率（图8b），如分别在35、30次迭代时IST、FIST 算法得到了收敛解。

图5 应用阈值指数递减策略对图2e插值的信噪比曲线

图6 使用阈值指数递减策略迭代50次对图2e的插值结果及其与图2b数据的差异

图7 使用阈值指数递减策略迭代50次对图2f的插值结果及其与图2c数据的差异

图8 应用阈值指数递减策略对图2f插值的信噪比曲线

3 结论

本文研究了加速迭代的快速凸集投影算法和快速迭代阈值算法，基本思想是：将前一步的插值结果与前两步的插值结果通过线性算子进行线性组合，得到迭代收缩算子；通过插值算法进行插值。同时引入质量控制新准则，提高了计算效率和精度。仿真实验和实际应用证明：阈值指数递减策略较恒阈值、阈值线性递减、数据驱动阈值等策略获得的插值结果的信噪比更高；结合终止准则，最大迭代次数为35～50时，可以获得较好的插值效果。