地震資料快速兩步插值算法

馬澤川 李 勇*② 陳力鑫 陳 杰 王鵬飛 李雪梅

（①成都理工大學地球物理學院，四川成都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時，可以獲得較好的插值效果。