２０２０年１０
月第５５卷　第５期　
＊四川省成都市成华区二仙桥东三路１号成都理工大学地球物理学院，６１００５９。Ｅｍａｉｌ：ｌｉｙｏｎｇ
０７＠ｃｄｕｔ．ｃｎ本文于２０１９年１２月３０日收到，最终修改稿于２０２０年６月９
日收到。本项研究受国家科技重大专项“大型油气田及煤层气开发重大专项———时频聚集流体识别方法研究”（２０１６ＺＸ０５０２６００１－
００４）、四川省重点研发计划项目“基于人工智能去噪的页岩气储层微裂缝识别研究”（２０２０ＹＦＧ０１５７
）联合资助。·处理技术·
文章编号：１０００－７２１０（２０２０）０５－０９９７－
０８地震资料快速两步插值算法
马泽川①　李　勇＊①②　陈力鑫①　陈　杰①　王鹏飞①　李雪梅①
（①成都理工大学地球物理学院，
四川成都６１００５９；②油气藏地质及开发工程国家重点实验室（成都理工大学），四川成都６１００５９）马泽川，李勇，陈力鑫，陈杰，王鹏飞，李雪梅．地震资料快速两步插值算法．石油地球物理勘探，２０２０，５５（５）：９９７－１
００４．摘要　为了提高插值效率以及选择最优的插值方案，基于凸集投影（ＰＯＣＳ）算法和迭代阈值（ＩＳＴ）算法的分析公式，在前人的基础上发展了快速迭代收缩阈值（ＦＩＳＴ）算法和快速凸集投影（ＦＰＯＣＳ）算法。基本思想是：将前一步的插值结果与前两步的插值结果通过线性算子进行线性组合，得到迭代收缩算子；通过插值算法进行插值。同时引入质量控制新准则，提高了计算效率和精度。使用ＩＳＴ、ＰＯＣＳ、ＦＩＳＴ和ＦＰＯＣＳ等算法分别对由Ｓｅｉｓｍｉｃ　Ｌａｂ建立的四层地震模型、Ｍａｒｍｏｕｓｉ模型的不完整地震数据进行插值，筛选出最佳的阈值策略，并最终由实际地震资料进行验证。结果表明：阈值指数递减策略较恒阈值、阈值线性递减、数据驱动阈值等策略获
得的插值结果的信噪比更高；结合终止准则，最大迭代次数为３５～５０时，即可获得较好的插值效果。关键词　快速迭代收缩阈值算法　快速凸集投影算法　阈值策略　终止准则　地震资料插值中图分类号：Ｐ６３１  文献标识码：Ａ  ｄｏｉ：１０．１３８１０／ｊ
．ｃｎｋｉ．ｉｓｓｎ．１０００－７２１０．２０２０．０５．００８０　引言
实际地震数据在空间方向往往存在不规则的缺
失道［１－
２］，要想得到完整的地震数据，通常需要插值处理［
３］。地震资料插值方法大致可以分为两类［４－５］：①利用线性空间预测滤波器插值
［６－
７］；②借助数学
变换（如傅里叶变换）插值［８－
９］，其中凸集投影（ＰＯＣＳ）算法和迭代阈值（ＩＳＴ）
算法是较为经典的算法［１０］
，利用地震数据变换域的稀疏性对地震数据插值，由于原理直观、实现简单，并且具有稳定性，因此发展
迅速。Ｍｅｎｋｅ［１１］
认为ＰＯ
ＣＳ算法具有普适性，并将其引入地球物理领域。近年来，各种基于传统ＰＯＣＳ算法的改进方法被提出且被广泛应用。刘红喜
等［１２］提出了边缘保持的ＰＯＣＳ超分辨率重建方法；王本锋等［１３］引入了结合ＰＯＣＳ和Ｃｕｒｖｅｌｅｔ的超分辨率重建算法重建三维数字内核；Ａｂｍａ等［９］将
ＰＯ
ＣＳ算法引入地震插值领域，获得了很好的效果；王本锋等［１４］提出了ＰＯＣＳ插值重建地震数据的质
量控制新准则，在保证插值精度的同时提高了计算
效率。在ＰＯ
ＣＳ算法发展的同时，人们也对ＩＳＴ算法进行研究。Ｈｅｒｒｍａｎｎ等［１
５］
基于基追踪降噪模型，将ＩＳＴ算法用于非均匀记录地震道插值；
由于现有ＩＳＴ算法收敛速度较慢，Ｂｅｃｋ等［１
６］从一般的数学背景出发，提出了快速迭代收缩阈值（ＦＩＳＴ）
算法；Ｇａｏ等［１
７］使用指数递减阈值策略代替线性阈值策略。近年来，业界着重对ＩＳＴ和ＰＯＣＳ算法迭代慢的不足进行优化，提高了算法的计算效率和精度。
为了提高插值效率以及选择最优的插值方案，本文基于ＩＳＴ算法和ＰＯＣＳ算法的分析公式，发展了ＦＩＳＴ算法和快速凸集投影（ＦＰＯＣＳ）算法。基本思想是：将前一步的插值结果与前两步的插值结果通过线性算子进行线性组合，得到迭代收缩算子；通过插值算法进行插值。同时引入质量控制新准
则［１４］
，提高了计算效率和精度。通过理论模型筛选
出最优的阈值策略以及最大迭代次数，实际资料处理结果表明，文中方法在较少迭代次数下的插值效果较好，验证了算法的高效性。
　９９８　石油地球物理勘探２０２０年　
１　基本原理
１．１　地震资料插值原理
将采集到的不规则的随机地震资料记为ｄｏｂ
ｓ，期望得到的完整数据记为ｄ，
则ｄｏｂ
ｓ＝Ｍｄ（１）式中Ｍ为采样算子，是由元素０和１组成的对角矩阵。
为了方便插值，根据压缩感知原理［１８－２０］
，将地
震数据投影到变换域，即
ｄｏｂｓ＝Ｍ
ｄ＝ＭＡχ（２
）式中：Ａ∈Ｒ
ｎ×ｍ
为正交坐标系的合成算子，且是一个紧框架；χ为变换域的表示系数。
然而式（２）是一个欠定方程，具有多解性，因此需要添加一些限制性条件求解。由于稀疏性在变换
域内始终适用［２１］
，通常与χ相对应。因此，求解
式（２）就转化为Ｌ０范数的最小化问题
ｍｉｎχ∈Ｒ
ｎ
‖ｄｏｂｓ－Ｋ
χ‖０（３
）式中Ｋ＝ＭＡ。由于Ｌ０范数的最小化非凸复杂
性［２２］
，常常使用Ｌ１范数代替Ｌ０范数，并与Ｌ０范数的稀疏性相同
［１８］
。因此，式（３
）变为  ｍ
ｉｎχ
Ｊ（χ）＝１２‖ｄｏｂｓ－Ｋχ‖２
２＋
τ‖χ‖１（４）式中τ为正则化参数。由于Ａ为紧框架，根据Ａ＊Ａ
＝Ｉ（Ｉ为单位矩阵）
以及χ＝Ａ＊
Ａχ＝Ａ＊
ｄ（“上角＊”表示伴随），则
　ｍｉｎｄＪ（ｄ）＝１２‖ｄｏｂｓ－Ｋχ‖２
２＋
τ‖Ａ＊ｄ‖１（５）式（５
）为地震插值分析公式，因为该式直接将地震资料ｄ作为未知数进行分析［
２３］
。１．２　ＩＳＴ算法
有很多算法求解式（４）。Ｆｉｇｕｅｉｒｅｄｏ等［２４］
和Ｄａｕｂｅｃｈｉｅｓ等［２５］
提出了ＩＳＴ算法，
一般表示为χ
（ｋ＋
正则化是最小化策略的实现１）＝Ｔλ｛
χ（ｋ）＋Ｋ＊［ｄｏｂｓ－Ｋ
χ（ｋ
）
］｝ｋ＝１，２，…，Ｎ
（６
）式中：λ为正阈值参数；ｋ为迭代次数，Ｎ为最大迭代次数；Ｔλ为阈值函数，
与τ的取值有关，当λ＝τ时，Ｔλ为软阈值函数，
定义软阈值算子为Ｓλ（
ｕ）＝ｕ－λｕ｜ｕ｜｜
ｕ｜＞λ０｜ｕ｜≤烅烄烆λ
（７
）式中ｕ为单道地震数据。与式（５）相同的思路，将式（６）两边左乘Ａ，得到ＩＳＴ算法公式
ｄ（ｋ＋１）
＝Ａ
χ（ｋ＋
１）＝ＡＴλ｛χ（ｋ）＋Ｋ＊［ｄｏｂ
ｓ－Ｋχ（ｋ
）
］｝＝ＡＴλ｛Ａ＊ｄ（ｋ）
＋（ＭＡ）＊
［ｄｏｂ
ｓ－ＭＡχ（ｋ）
］｝＝ＡＴλ｛Ａ＊［ｄ（ｋ）＋Ｍ＊（ｄｏｂｓ－Ｍｄ（ｋ）
）］｝＝ＡＴλ｛Ａ＊［ｄｏｂｓ＋（
Ｉ－Ｍ）ｄ（ｋ）
］｝（８
）式中ｄ（ｋ）
、χ（ｋ）
分别为第ｋ次迭代时的ｄ与χ的值。
１．３　ＰＯＣＳ算法
定义
ｕ（ｋ）＝ｄｏｂ
ｓ＋（Ｉ－Ｍ）ｄ（ｋ）
（９）式中ｕ（
ｋ）
为第ｋ次迭代后原始地震道的插入数据。对于式（８
）
，有ｄ（ｋ＋１）＝ＡＴλ［Ａ＊ｕ（ｋ）
］
在第ｋ＋１次迭代时，结合式（９
），有  ｕ（ｋ＋１）＝ｄｏｂ
ｓ＋（Ｉ－Ｍ）ＡＴλ［Ａ＊ｕ（ｋ）
］（１０
）需要注意的是，ｌｉｍｋ→∞
ｄ（ｋ）
＝ｄ保证了ＩＳＴ算法的收
敛［１３］
，因此对于式（１０
）
，有  ｌ
ｉｍｋ→∞
ｕ（ｋ）
＝ｌｉｍｋ→∞
［ｄｏｂ
ｓ＋（Ｉ－Ｍ）ｄ（ｋ）
］＝ｄ（１１）由Ｍ（Ｉ－Ｍ）＝０，
得Ｍｕ（ｋ）
＝ｄｏｂｓ
（１２）在式（１０
）中用ｄ（ｋ）取代ｕ（ｋ）
，
得  ｄ（ｋ＋１）＝ｄｏｂｓ＋（Ｉ－Ｍ）ＡＴλ［Ａ＊ｄ（ｋ）］（１３
）这样，便得到ＰＯ
ＣＳ算法，该算法已广泛用于许多领域，如图像重建和修复等。１．４　两步插值：ＦＩ
ＳＴ和ＦＰＯＣＳ由于ＩＳＴ算法的收敛速度很慢。Ｂｅｃｋ等［１
６］
提出了快速迭代收缩阈值算法（ＦＩＳＴ）
，可以表示为珘ｘ（ｋ）
＝ｘ（ｋ）
＋ｔ（ｋ）
－１ｔ（ｋ＋
１）［ｘ（ｋ）－ｔ（ｋ－１）］ｘ（ｋ＋１）＝Ｔλ｛
珘ｘ（ｋ）＋Ｋ＊［ｄｏｂｓ－Ｋ珘ｘ（ｋ）
烅
烄烆］｝（１４）式中：珘ｘ（ｋ）
为迭代收缩算子；ｘ（ｋ）为第ｋ次迭代后的ｘ值；ｔ为迭代过程中的线性算子［１６］，ｔ（１）＝１，ｔ（
ｋ＋１）＝
１＋１＋４［ｔ（ｋ）］槡２
２。式（１４
）代表两步插值的过程：将前一步的插值结果ｘ（ｋ）与前两步的插值结果ｘ（ｋ－１）
通过线性算子ｔ进行线性组合，得到迭代收缩算子珘ｘ（ｋ）；通过插值算
法计算ｘ（ｋ＋１）
，提高了迭代速度。
与式（８）推导相同，式（１４）
两边左乘Ａ得到ＦＩＳＴ算法公式
　第５
５卷　第５期马泽川，
等：地震资料快速两步插值算法９
９９  珘ｄ（ｋ）＝ｄ（ｋ）＋ｔ（ｋ）
－１ｔ（ｋ＋
１）［ｄ（ｋ）－ｄ（ｋ－１）］ｄ（ｋ＋１）＝ＡＴλ｛Ａ＊［ｄｏｂｓ＋（Ｉ－Ｍ）珘ｄ（ｋ）烅
烄烆］｝（１５）式中：珘ｄ（ｋ）为地震数据插值的迭代收缩算子；ｄ（ｋ）为
第ｋ次迭代后的地震数据。
在式（１５）的基础上，按照式（１３）思路，可以推出ＦＰＯＣＳ算法公式
珘ｄ（ｋ）＝ｄ（ｋ）＋ｔ（ｋ）
－１ｔ（ｋ＋
１）［ｄ（ｋ）－ｄ（ｋ－１）］ｄ（ｋ＋１）＝ｄｏｂ
ｓ＋（Ｉ－Ｍ）ＡＴλ［Ａ＊珘ｄ（ｋ）烅
烄烆］（１６）１．
５　两步插值方法的阈值策略正则化参数τ具重要作用，并且与每次迭代的阈值λ密切相关。因此定义软阈值
λ＝τ
（１７）通过对比不同阈值策略的实际效果，筛选最佳的阈值策略。
（１
）恒阈值τｋ＝Ｃ  　
ｋ＝１，２，…，Ｎ（１８
）式中Ｃ为常数。
（２
）阈值线性递减τｋ＝τｍａ
ｘ－ｋ－１Ｎ－１
（τｍａｘ－τｍｉ
ｎ）  ｋ＝１，２，…，Ｎ（１９
）  （
３）阈值指数递减  τｋ＝τｍａ
ｘτ（）
ｍｉ
ｎｋ－
１Ｎ－１
τｍａｘ 
ｋ＝１，２，…，Ｎ（２０
）  （
４）数据驱动阈值τ１＝ｖ１τｋ＝ｖｊ
τＮ＝ｖＮ
ｊ＝「
（ｋ－１）Ｎｖ／（Ｎ－１）烅烄烆
?
（２１
）且
τｍａｘ＝ｐｍａｘ·ｍａｘ｛｜（Ａ＊ｄｏｂｓ）ｉ｜｝τｍｉｎ＝ｐｍｉｎ·ｍａｘ｛｜（Ａ＊ｄｏｂｓ）ｉ｜烅
烄
烆
｝（２２
）式中：τｍａｘ和τｍｉｎ分别为所选的最大和最小正则化参数；ｐｍａｘ和ｐｍｉｎ分别为所定义的最大和最小百分比；序列｛｜（Ａ＊
ｄｏｂｓ）ｉ｜｝∈［τｍｉｎ，τｍａ
ｘ］按其中元素振幅的下降顺序排列可以得到向量ｖ，Ｎｖ为ｖ的长度，ｖｊ为
ｖ中的第ｊ个元素；「
（ｋ－１）Ｎｖ／（Ｎ－１）?表示取不小于（ｋ－１）Ｎｖ／（Ｎ－１
）的最小整数。以上介绍了两步插值算法的原理。由于地震资
料在变换域具有稀疏性，根据压缩感知理论可知，将地震数据投影到变换域，选取阈值策略进行反复迭代，即可完成地震数据插值。
２　仿真实验和实际应用
２．１　两步插值法的终止准则
文中介绍的ＩＳＴ、ＰＯＣＳ、ＦＩＳＴ和ＦＰＯＣＳ等
４种插值算法都属于迭代方法，
这就意味着在达到所定义的最大迭代次数Ｎ之前，算法不会终止。因此，需要建立一个终止准则适时地结束迭代。Ｇａｏ
等［４－
５］等给出了两种终止准则，王本锋等［
１４］在此基础上提出了灵活、适用性更好的质量控制新准则
‖ｄ（ｋ）
－ｕ（ｋ）
‖２
２
‖ｄｏｂｓ‖２
２
＜η（２３
）进一步提高了计算效率。式中η为公差，需要在程
序运行前设置。
文中使用信噪比
ＳＮＲ＝１０ｌｇ‖ｄ‖２
２
‖ｄ－＾ｄ‖２
２
（２４
）评估插值质量。式中＾ｄ为重建后的地震记录。
综上所述，总结了地震资料插值流程（图１
）
。图１　地震资料插值流程
２．
２　算例模型使用ＩＳＴ、ＰＯＣＳ、ＦＩＳＴ和ＦＰＯＣＳ等算法分别对由Ｓｅｉｓｍｉｃ　Ｌａｂ
建立的四层地震模型（图２ａ）、Ｍａｒｍｏｕｓｉ模型（图２ｂ）的不完整地震数据（图２ｄ、图２ｅ）进行插值，筛选出最佳的阈值策略，并最终由实际地震资料（图２ｃ、图２ｆ
）进行验证。
　１０００　石油地球物理勘探２０２０年　
２．３　最优阈值策略
首先在由Ｓｅｉｓｍｉｃ　Ｌａｂ建立的四层地震模型上使用不同的阈值策略插值，以此筛选最优的阈值策略。应用ＩＳＴ、ＰＯＣＳ、ＦＩＳＴ和ＦＰＯＣＳ等算法对缺失数据（图２ｄ）插值，并绘制信噪比曲线（图３）。可见：在恒阈值时的信噪比为３８．５８６３ｄＢ，两步插值法在５０次迭代时便完成了图像重建，较一步插值法的图像重建速度更快（图３ａ）；在阈值线性递减时的信噪比为３５．１８７５ｄＢ，
在迭代次数较少时图像重建效果不明显，随着迭代次数增加，重构作用越来越大，但会引入噪声，导致信噪比降低（图３ｂ）；在阈值指数递减时的信噪比为４３．４５７１ｄＢ（图３ｃ）；数据驱动阈值时的信噪比为４１．１８１８ｄＢ（图３ｄ）。因此，阈值指数递减与数据驱动阈值的插值结果信噪比较高，数据重建效果较好。图４为使用阈值指数递减策略对图２ｄ的插值结果及其与图２ａ数据的差异。由图可见，经过插值后的地震数据恢复完好，且没有引入大量噪声
。
图２　实验算例
（ａ）由Ｓｅｉｓｍｉｃ　
Ｌａｂ建立的四层地震模型（３５Ｈｚ雷克子波），共５１２道，每道５１２个样点，采样间隔为０．００２ｓ；（ｂ）Ｍａｒｍｏｕｓｉ地震模型，共２０１道，每道６００个样点，采样间隔为０．００３ｓ；（ｃ）实际地震资料，共４３０道，每道５１２个样点；（ｄ）～（ｆ）对应图ａ～图ｃ缺失３
０％的数据２．
４　合理的最大迭代次数利用Ｍａｒｍｏｕｓｉ模型验证两步插值方法的收敛性。图５为应用阈值指数递减策略对图２ｅ插值的信噪比曲线。由图可见：在迭代次数较小时，重建效果较差，两步插值法的收敛速度略快于一步插值法，因此两步插值法可提高收敛速度；由于观测数据的复杂性，对Ｍａｒｍｏｕｓｉ地震模型的重构精度不高；
在５０次迭代后得到了较好的收敛解（图５ｃ），再增加迭代次数对插值精度的改善很小（图５ｄ），通过多次测试发现，当Ｎ∈［３５，５０］时可以取得满意的插值效果。图６为使用阈值指数递减策略迭代５０次对图２ｅ的插值结果及其与图２ｂ数据的差异。由图可见，四种算法具有相近的插值效果，在没有引入大量
噪声的情况下，地震数据恢复完好。
此外，在迭代次数较少时ＰＯＣＳ、ＦＰＯＣＳ算法的信噪比始终高于ＩＳＴ、ＦＩＳＴ算法，
这是由于所有算法在初始化时将ｄ（ｋ）
置０所致。由式（１２）可知，在ＰＯＣＳ、ＦＰＯＣＳ算法中，
有Ｍｄ（ｋ＋１）
＝ｄｏｂ
ｓ＝Ｍｄ    ｋ（２５
）因此
　‖ｄ（ｋ＋１）－ｄ‖＝‖（Ｉ－Ｍ）［ｄ（ｋ＋１）
－ｄ］‖２
（２６）这是ＰＯ
ＣＳ、ＦＰＯＣＳ算法的信噪比曲线的起点总是不为０的原因。由式（８）、式（２１）可知，在ＩＳＴ、ＦＩＳＴ算法中，
开始设计了较大的阈值，几乎忽略了所有的傅里叶系数，导致ｄ（ｋ＋１）
中的值过多，使
‖ｄ‖２
２≈‖ｄ－ｄ（ｋ＋１）‖２２、
ＳＮＲ≈０。
　第５５卷　第５期马泽川，
等：地震资料快速两步插值算法１
００１　
　图３　应用不同阈值策略对图２ｄ插值的信噪比曲线
（ａ）恒阈值（τ＝０．
００５，Ｎ＝２００）；（ｂ）阈值线性递减（Ｎ＝２００）；（ｃ）阈值指数递减（Ｎ＝２００）；（ｄ）数据驱动阈值（Ｎ＝２００
）图４　使用阈值指数递减策略对图２ｄ的插值结果及其与图２ａ数据的差异
（ａ）ＩＳＴ软阈值插值（ＳＮＲ＝３８．５８６３ｄＢ）；（ｂ）ＦＩＳＴ软阈值插值（ＳＮＲ＝４０．６５１４ｄＢ）；（ｃ）ＰＯＣＳ软阈值插值（ＳＮ
Ｒ＝３８．５８６４ｄＢ；（ｄ）ＦＰＯＣＳ软阈值插值（ＳＮＲ＝４０．６５１６ｄＢ）；（ｅ）～（ｈ）对应图ａ～图ｄ数据与图２
ａ数据的差值２．５　合理的终止准则以及实际地震数据插值效果
通过对不完整的理论模型数据（图２ｄ、图２ｅ）进行插值，筛选出最优的阈值策略和最大迭代次数，在此基础上验证实际地震数据（图２ｃ、图２ｆ
）的插值效果。图７为使用阈值指数递减策略迭代５０次对图２ｆ的插值结果及其与图２ｃ数据的差异，图８为应用阈值指数递减策略对图２ｆ插值的信噪比曲线。由图可见：①四种算法插值结果的信噪比约为１
４ｄＢ，

版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。