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

（①成都理工大學地球物理學院，四川成都610059；②油氣藏地質(zhì)及開發(fā)工程國家重點實驗室（成都理工大學），四川成都610059）

0 引言

實際地震數(shù)據(jù)在空間方向往往存在不規(guī)則的缺失道［1-2］，要想得到完整的地震數(shù)據(jù)，通常需要插值處理［3］。地震資料插值方法大致可以分為兩類［4-5］：①利用線性空間預(yù)測濾波器插值［6-7］；②借助數(shù)學變換（如傅里葉變換）插值［8-9］，其中凸集投影（POCS）算法和迭代閾值（IST）算法是較為經(jīng)典的算法［10］，利用地震數(shù)據(jù)變換域的稀疏性對地震數(shù)據(jù)插值，由于原理直觀、實現(xiàn)簡單，并且具有穩(wěn)定性，因此發(fā)展迅速。Menke［11］認為POCS算法具有普適性，并將其引入地球物理領(lǐng)域。近年來，各種基于傳統(tǒng)POCS算法的改進方法被提出且被廣泛應(yīng)用。劉紅喜等［12］提出了邊緣保持的POCS超分辨率重建方法；王本鋒等［13］引入了結(jié)合POCS和Curvelet的超分辨率重建算法重建三維數(shù)字內(nèi)核；Abma 等［9］將POCS算法引入地震插值領(lǐng)域，獲得了很好的效果；王本鋒等［14］提出了POCS插值重建地震數(shù)據(jù)的質(zhì)量控制新準則，在保證插值精度的同時提高了計算效率。在POCS算法發(fā)展的同時，人們也對IST 算法進行研究。Herrmann 等［15］基于基追蹤降噪模型，將IST 算法用于非均勻記錄地震道插值；由于現(xiàn)有IST 算法收斂速度較慢，Beck 等［16］從一般的數(shù)學背景出發(fā)，提出了快速迭代收縮閾值（FIST）算法；Gao等［17］使用指數(shù)遞減閾值策略代替線性閾值策略。近年來，業(yè)界著重對IST和POCS算法迭代慢的不足進行優(yōu)化，提高了算法的計算效率和精度。

為了提高插值效率以及選擇最優(yōu)的插值方案，本文基于IST 算法和POCS算法的分析公式，發(fā)展了FIST 算法和快速凸集投影（FPOCS）算法。基本思想是：將前一步的插值結(jié)果與前兩步的插值結(jié)果通過線性算子進行線性組合，得到迭代收縮算子；通過插值算法進行插值。同時引入質(zhì)量控制新準則［14］，提高了計算效率和精度。通過理論模型篩選出最優(yōu)的閾值策略以及最大迭代次數(shù)，實際資料處理結(jié)果表明，文中方法在較少迭代次數(shù)下的插值效果較好，驗證了算法的高效性。

1 基本原理

1.1 地震資料插值原理

將采集到的不規(guī)則的隨機地震資料記為dobs，期望得到的完整數(shù)據(jù)記為d，則

式中M 為采樣算子，是由元素0 和1 組成的對角矩陣。

為了方便插值，根據(jù)壓縮感知原理［18-20］，將地震數(shù)據(jù)投影到變換域，即

式中：A∈Rn×m為正交坐標系的合成算子，且是一個緊框架；χ 為變換域的表示系數(shù)。

然而式（2）是一個欠定方程，具有多解性，因此需要添加一些限制性條件求解。由于稀疏性在變換域內(nèi)始終適用［21］，通常與χ 相對應(yīng)。因此，求解式（2）就轉(zhuǎn)化為L0范數(shù)的最小化問題

式中K＝MA。由于L0范數(shù)的最小化非凸復(fù)雜性［22］，常常使用L1范數(shù)代替L0范數(shù)，并與L0范數(shù)的稀疏性相同［18］。因此，式（3）變?yōu)?/p>

式中τ為正則化參數(shù)。由于A 為緊框架，根據(jù)A*A＝I（I為單位矩陣）以及χ＝A*Aχ＝A*d（“上角*”表示伴隨），則

式（5）為地震插值分析公式，因為該式直接將地震資料d作為未知數(shù)進行分析［23］。

1.2 IST算法

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

式中：λ為正閾值參數(shù)；k 為迭代次數(shù)，N 為最大迭代次數(shù)；Tλ為閾值函數(shù)，與τ 的取值有關(guān)，當λ＝τ時，Tλ為軟閾值函數(shù)，定義軟閾值算子為

式中u 為單道地震數(shù)據(jù)。與式（5）相同的思路，將式（6）兩邊左乘A，得到IST 算法公式

式中d（k）、χ（k）分別為第k次迭代時的d與χ 的值。

1.3 POCS算法

定義

式中u（k）為第k 次迭代后原始地震道的插入數(shù)據(jù)。對于式（8），有

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

在第k＋1次迭代時，結(jié)合式（9），有

由M（I－M）＝0，得

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

這樣，便得到POCS算法，該算法已廣泛用于許多領(lǐng)域，如圖像重建和修復(fù)等。

1.4 兩步插值：FIST和FPOCS

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

式（14）代表兩步插值的過程：將前一步的插值結(jié)果x（k）與前兩步的插值結(jié)果x（k－1）通過線性算子t進行線性組合，得到迭代收縮算子～x（k）；通過插值算法計算x（k＋1），提高了迭代速度。

與式（8）推導(dǎo)相同，式（14）兩邊左乘A 得到FIST 算法公式

式中：～d（k）為地震數(shù)據(jù)插值的迭代收縮算子；d（k）為第k次迭代后的地震數(shù)據(jù)。

在式（15）的基礎(chǔ)上，按照式（13）思路，可以推出FPOCS算法公式

1.5 兩步插值方法的閾值策略

正則化參數(shù)τ具重要作用，并且與每次迭代的閾值λ密切相關(guān)。因此定義軟閾值

通過對比不同閾值策略的實際效果，篩選最佳的閾值策略。

（1）恒閾值

式中C 為常數(shù)。

（2）閾值線性遞減

（3）閾值指數(shù)遞減

（4）數(shù)據(jù)驅(qū)動閾值

且

式中：τmax和τmin分別為所選的最大和最小正則化參數(shù)；pmax和pmin分別為所定義的最大和最小百分比；序列｛|（A*dobs）i|｝∈［τmin，τmax］按其中元素振幅的下降順序排列可以得到向量v，Nv為v 的長度，vj為v 中的第j 個元素；「（k－1）Nv／（N－1）?表示取不小于（k－1）Nv／（N－1）的最小整數(shù)。

以上介紹了兩步插值算法的原理。由于地震資料在變換域具有稀疏性，根據(jù)壓縮感知理論可知，將地震數(shù)據(jù)投影到變換域，選取閾值策略進行反復(fù)迭代，即可完成地震數(shù)據(jù)插值。

2 仿真實驗和實際應(yīng)用

2.1 兩步插值法的終止準則

文中 介 紹 的IST、POCS、FIST 和FPOCS 等4種插值算法都屬于迭代方法，這就意味著在達到所定義的最大迭代次數(shù)N 之前，算法不會終止。因此，需要建立一個終止準則適時地結(jié)束迭代。Gao等［4-5］等給出了兩種終止準則，王本鋒等［14］在此基礎(chǔ)上提出了靈活、適用性更好的質(zhì)量控制新準則

進一步提高了計算效率。式中η 為公差，需要在程序運行前設(shè)置。

文中使用信噪比

評估插值質(zhì)量。式中^d 為重建后的地震記錄。

綜上所述，總結(jié)了地震資料插值流程（圖1）。

圖1 地震資料插值流程

2.2 算例模型

使用IST、POCS、FIST 和FPOCS等算法分別對由Seismic Lab 建立的四層地震模型（圖2a）、Marmousi模型（圖2b）的不完整地震數(shù)據(jù)（圖2d、圖2e）進行插值，篩選出最佳的閾值策略，并最終由實際地震資料（圖2c、圖2f）進行驗證。

2.3 最優(yōu)閾值策略

首先在由Seismic Lab建立的四層地震模型上使用不同的閾值策略插值，以此篩選最優(yōu)的閾值策略。應(yīng)用IST、POCS、FIST 和FPOCS等算法對缺失數(shù)據(jù)（圖2d）插值，并繪制信噪比曲線（圖3）?？梢姡涸诤汩撝禃r的信噪比為38.5863 d B，兩步插值法在50次迭代時便完成了圖像重建，較一步插值法的圖像重建速度更快（圖3a）；在閾值線性遞減時的信噪比為35.1875dB，在迭代次數(shù)較少時圖像重建效果不明顯，隨著迭代次數(shù)增加，重構(gòu)作用越來越大，但會引入噪聲，導(dǎo)致信噪比降低（圖3b）；在閾值指數(shù)遞減時的信噪比為43.4571dB（圖3c）；數(shù)據(jù)驅(qū)動閾值時的信噪比為41.1818d B（圖3d）。因此，閾值指數(shù)遞減與數(shù)據(jù)驅(qū)動閾值的插值結(jié)果信噪比較高，數(shù)據(jù)重建效果較好。圖4為使用閾值指數(shù)遞減策略對圖2d的插值結(jié)果及其與圖2a數(shù)據(jù)的差異。由圖可見，經(jīng)過插值后的地震數(shù)據(jù)恢復(fù)完好，且沒有引入大量噪聲。

圖2 實驗算例

2.4 合理的最大迭代次數(shù)

利用Marmousi模型驗證兩步插值方法的收斂性。圖5為應(yīng)用閾值指數(shù)遞減策略對圖2e插值的信噪比曲線。由圖可見：在迭代次數(shù)較小時，重建效果較差，兩步插值法的收斂速度略快于一步插值法，因此兩步插值法可提高收斂速度；由于觀測數(shù)據(jù)的復(fù)雜性，對Marmousi地震模型的重構(gòu)精度不高；在50次迭代后得到了較好的收斂解（圖5c），再增加迭代次數(shù)對插值精度的改善很?。▓D5d），通過多次測試發(fā)現(xiàn)，當N∈［35，50］時可以取得滿意的插值效果。圖6 為使用閾值指數(shù)遞減策略迭代50 次對圖2e的插值結(jié)果及其與圖2b數(shù)據(jù)的差異。由圖可見，四種算法具有相近的插值效果，在沒有引入大量噪聲的情況下，地震數(shù)據(jù)恢復(fù)完好。

此外，在迭代次數(shù)較少時POCS、FPOCS 算法的信噪比始終高于IST、FIST 算法，這是由于所有算法在初始化時將d（k）置0所致。由式（12）可知，在POCS、FPOCS算法中，有

因此

圖4 使用閾值指數(shù)遞減策略對圖2d的插值結(jié)果及其與圖2a數(shù)據(jù)的差異

2.5 合理的終止準則以及實際地震數(shù)據(jù)插值效果

通過對不完整的理論模型數(shù)據(jù)（圖2d、圖2e）進行插值，篩選出最優(yōu)的閾值策略和最大迭代次數(shù)，在此基礎(chǔ)上驗證實際地震數(shù)據(jù)（圖2c、圖2f）的插值效果。圖7 為使用閾值指數(shù)遞減策略迭代50 次對圖2f的插值結(jié)果及其與圖2c數(shù)據(jù)的差異，圖8 為應(yīng)用閾值指數(shù)遞減策略對圖2f插值的信噪比曲線。由圖可見：①四種算法插值結(jié)果的信噪比約為14dB，均能較好地重建缺失數(shù)據(jù)（圖7）。②兩步插值算法相對于單步插值算法具有更快的收斂速度（圖8）；與沒有加入終止準則的結(jié)果（圖8a）相比，應(yīng)用終止準則后，算法能更快地結(jié)束迭代，從而提高計算效率（圖8b），如分別在35、30次迭代時IST、FIST 算法得到了收斂解。

圖5 應(yīng)用閾值指數(shù)遞減策略對圖2e插值的信噪比曲線

圖6 使用閾值指數(shù)遞減策略迭代50次對圖2e的插值結(jié)果及其與圖2b數(shù)據(jù)的差異

圖7 使用閾值指數(shù)遞減策略迭代50次對圖2f的插值結(jié)果及其與圖2c數(shù)據(jù)的差異

圖8 應(yīng)用閾值指數(shù)遞減策略對圖2f插值的信噪比曲線

3 結(jié)論

本文研究了加速迭代的快速凸集投影算法和快速迭代閾值算法，基本思想是：將前一步的插值結(jié)果與前兩步的插值結(jié)果通過線性算子進行線性組合，得到迭代收縮算子；通過插值算法進行插值。同時引入質(zhì)量控制新準則，提高了計算效率和精度。仿真實驗和實際應(yīng)用證明：閾值指數(shù)遞減策略較恒閾值、閾值線性遞減、數(shù)據(jù)驅(qū)動閾值等策略獲得的插值結(jié)果的信噪比更高；結(jié)合終止準則，最大迭代次數(shù)為35～50時，可以獲得較好的插值效果。