溫 睿 劉國昌* 冉 揚

(①中國石油大學(北京)油氣資源與探測國家重點實驗室，北京 102249；②中國石油大慶油田測試技術服務分公司，黑龍江大慶 163000)

1 引言

在地震數(shù)據(jù)采集過程中，受物理因素或經(jīng)濟條件限制，所采集數(shù)據(jù)往往會出現(xiàn)采樣不足的現(xiàn)象，比如采集時受障礙物、風化帶等影響，造成地震道缺失。另外有時出于經(jīng)濟成本的考慮，不得不降低采樣率，這些因素都會造成觀測地震數(shù)據(jù)的不規(guī)則化，導致假頻的產(chǎn)生，從而對后續(xù)的地震數(shù)據(jù)處理成像造成不良影響。為了解決這個問題，需對地震數(shù)據(jù)進行插值處理，對缺失的地震數(shù)據(jù)進行重建，地震數(shù)據(jù)插值重建對后續(xù)的SRME多次波去除，波動方程偏移成像等環(huán)節(jié)具有重要的意義。

地震數(shù)據(jù)重建方法大致可分為兩大類[1]，第一大類是基于波動方程的方法[2]，此類方法利用地下信息進行數(shù)據(jù)重建，但是計算成本較高，可細分為基于動校正以及基于成像算子的方法。另一大類是基于信號分析和統(tǒng)計的方法，可細分為三個分支： ①基于預測濾波的方法[3-6]，通過低頻部分的濾波算子預測高頻部分的信息，實現(xiàn)數(shù)據(jù)重建； ②基于降秩的方法[7，8]，將地震數(shù)據(jù)看作矩陣或者張量，數(shù)據(jù)的缺失和不規(guī)則會導致矩陣秩的增加，通過SVD分解等矩陣分解算法達到降秩的目的，從而重建缺失數(shù)據(jù)； ③基于稀疏變換的壓縮感知方法，壓縮感知理論認為：如果信號是稀疏的或者通過稀疏變換可以在某個變換域得到稀疏表示，設計一個與變換基不相關的測量矩陣進行觀測并得到少量觀測數(shù)據(jù)，那么通過各種求解最優(yōu)化問題的重構算法可大概率地將觀測數(shù)據(jù)恢復為原來的數(shù)據(jù)[9，10]。這些重構算法常用的稀疏變換包括Fourier變換[11]、Radon變換[12，13]、Curvelet變換[14，15]、Seislet變換[16，17]和Shearlet變換[18]等。近年來，隨著壓縮感知理論的發(fā)展，基于稀疏變換的重建算法在地震數(shù)據(jù)重建領域流行起來，其中凸集投影算法(POCS)以其簡便、高效的特點而被廣泛應用。該算法由Bregman[19]提出，被Stark等[20]應用于圖像處理領域， Abma等[21]進一步將其應用于地震數(shù)據(jù)重建領域。以硬閾值迭代算法(IHT)為代表的稀疏促進迭代閾值類算法與POCS原理相近，同樣被應用到地震數(shù)據(jù)插值重建。Herrmann等[22]提出了基于曲波變換的促稀疏算法CRSI(Curvelet Recovery by Sparsity-promoting Inversion)。近幾年在圖像處理領域開始流行的線性Bregman迭代算法也開始應用于地震數(shù)據(jù)重建[23-25]。線性Bregman迭代算法雖然收斂速率快，但是最終的重建質量并不理想。Bai等[26]考慮到線性Bregman迭代算法與CRSI算法的優(yōu)缺點，提出了基于曲波變換聯(lián)合促稀疏反演算法(JRSI)，通過加權的方法將Bregman迭代算法與CRSI算法相結合，其收斂速度比CRSI快，且重建質量優(yōu)于傳統(tǒng)的Bregman迭代算法。

Abma等[21]的研究表明，當利用壓縮感知地震數(shù)據(jù)重建方法時，閾值模型的選取會影響重建質量，在迭代過程中使用了線性遞減的閾值模型。這之后，許多學者研究并提出了其他閾值模型進行地震數(shù)據(jù)的插值重建： Galloway等[27]提出了將指數(shù)閾值模型應用到POCS算法中，提高了收斂速率； 張華等[28]和王本鋒等[29]對指數(shù)模型做出了改進，取得了更好的重建效果； Gao等[30]提出了基于數(shù)據(jù)驅動的閾值模型，僅用少量的迭代即可達到較高的重建質量； 葛子建等[31]提出了反比例閾值模型，通過調節(jié)參數(shù)得到任意下降速率的閾值曲線，因此具有非常好的適用性。

壓縮感知地震數(shù)據(jù)重建算法受稀疏變換、迭代算法和閾值模型的影響，實際應用中選擇哪種方法較困難，因此本文結合理論推導和算例分析，對該類算法的三個關鍵因素進行對比研究，為在實際地震數(shù)據(jù)重建中選擇稀疏變換、迭代算法及閾值模型提供理論依據(jù)和建議。

本文論述了基于稀疏變換的壓縮感知地震數(shù)據(jù)重建算法： 首先對比了Fourier、Curvelet、Seislet三種稀疏變換，分析了不同稀疏變換的優(yōu)勢與局限性； 然后比較了POCS、IHT、Bregman和JRSI四種不同迭代算法的差異； 最后研究了線性、指數(shù)和數(shù)據(jù)驅動三種不同閾值模型的特點，并通過模擬和實際算例對比分析了壓縮感知地震數(shù)據(jù)重建過程中影響重建效果和效率的三方面關鍵因素。

2 基于稀疏變換的壓縮感知地震數(shù)據(jù)重建

不規(guī)則地震數(shù)據(jù)重建可描述為如下反演問題

dobs=Rm+n

(1)

式中：dobs為觀測到的不規(guī)則缺失地震數(shù)據(jù)；n為隨機噪聲；m是待求的重建之后的完整數(shù)據(jù)；R是由元素1和0構成的矩陣。R在壓縮感知領域被稱為采樣矩陣，其形式以如下例子作簡要說明

(2)

在不考慮噪聲且假設d1=m1，d3=m3，d5=m5的情況下，有

(3)

基于模擬計算方便考慮，本文假設dobs與m兩個向量的長度相等，將R視為mask矩陣[32]，即R為由元素1和0組成的對角陣，那么式(2)可改寫為

(4)

矩陣R是高度稀疏的，使式(1)成為一個欠定方程組。為了求解這個欠定問題，往往需要增加約束條件?？紤]到地震數(shù)據(jù)在各種變換域所表現(xiàn)的稀疏性，通過增加稀疏約束條件，將地震數(shù)據(jù)重建問題轉化為L0或者L1范數(shù)優(yōu)化問題

Φ(x)=‖dobs-RCTx‖22+λP(x)

(5)

式中：x為稀疏域中數(shù)據(jù)；CT為反稀疏變換；λ為正則化參數(shù)；P(x)為稀疏約束項，如L0或L1范數(shù)。地震數(shù)據(jù)插值重建是求取使目標函數(shù)(式(5))達到最小時的x。式(5)中稀疏變換C有多種可供選擇的類型，如Fourier變換、Radon變換、Curvelet變換、Seislet變換等，但不同變換的重建效果和計算效率存在差別。

當P(x)選為L1或L0范數(shù)時，優(yōu)化問題(式(5))是非線性的，可通過迭代類算法求解該非線性優(yōu)化問題。常用的解決L1范數(shù)約束問題的迭代算法有POCS、Bregman和JRSI等，解決L0范數(shù)約束問題的迭代算法有IHT算法等。POCS算法是求解非線性優(yōu)化問題(式(5))的常用迭代算法，其迭代公式為

mk+1=dobs+(I-R)CTTτkCmk

(6)

式中：mk為第k次迭代后數(shù)據(jù)空間域的解；C和CT分別為正、反稀疏變換；τk為閾值；Tτk為硬閾值算子，可由下式表示

(7)

另一種與POCS算法相似的IHT迭代算法可寫為[33]

xk+1=Tτk[xk+(RCT)T(dobs-RCTxk)]

(8)

式中xk為第k次迭代在稀疏域中的解。式(8)可轉化為

(9)

比較IHT迭代式(式(9))與POCS迭代式(式(6))可發(fā)現(xiàn)，理論上IHT迭代算法等價于在POCS迭代結束后又做了一次稀疏變換閾值濾波。若將最后一次迭代的閾值設為零，則兩種算法重建結果完全相同； 若最后一次迭代的閾值不為零，則兩種算法最后的重建結果會有差異。

Bregman迭代算法是Osher等[34]在研究全變分圖像去噪時提出的一種新型迭代正則化方法，該方法被拓展到基于小波的圖像去噪、非線性逆尺度空間、醫(yī)療核磁共振圖像等[35-37]問題。隨后又被應用到壓縮感知領域，成為求解L1范數(shù)以及L1范數(shù)相關最優(yōu)化問題最有效的方法之一[38]。近年來，Bregman迭代算法[23]開始應用到缺失地震數(shù)據(jù)重建問題。其算法如下

(10)

Bregman迭代算法還可發(fā)展為下面的迭代形式

(11)

式中：dk為每一次迭代的輸入數(shù)據(jù)；T為閾值算子。經(jīng)過N次迭代之后，得到稀疏域的解xN，反變換回時間域就得到數(shù)據(jù)域的解，也就是最終的重建結果mN=CTxN。

Bai等[26]提出一種新的促稀疏迭代反演算法，通過引入加權因子將CRSI迭代算法與Bregman迭代算法(式(10))結合起來，取得了不錯的重建效果。CRSI迭代算法原理與IHT迭代算法相似，仍然可以用式(8)進行迭代，CRSI迭代算法與IHT算法的區(qū)別在于CRSI算法的閾值不是每次迭代都會變化，而是每隔幾次迭代閾值變化一次。CRSI算法雖然穩(wěn)定，但收斂速率不快。Bai等[26]針對CRSI穩(wěn)定的特點以及Bregman算法收斂速率較快的特點，將二者結合，得到了聯(lián)合促稀疏迭代算法(JRSI)，其迭代公式如下

(12)

式中：vk為第k次迭代閾值之前的Curvelet系數(shù)；β為權重因子。為了實現(xiàn)快速高效的重建，權重因子β應該在前幾次迭代中取較小值，而在最后幾次迭代中取較大值，β可通過階梯函數(shù)或指數(shù)函數(shù)定義，即

(13)

(14)

式中Nd為小于迭代次數(shù)N的實數(shù)。以階梯函數(shù)定義權重因子，意味著在Nd+1次迭代之前，只有Bregman加速項對迭代結果有貢獻； 在Nd次迭代之后，只有穩(wěn)定項對迭代結果有貢獻。若以指數(shù)函數(shù)定義權重因子，意味著迭代前期vk權重大于xk，從而加快迭代速率。隨著迭代次數(shù)增多，vk權重減少，直到xk占全部比重，這時迭代的穩(wěn)定性可得到保證，重建之后能達到較高信噪比。

在壓縮感知地震數(shù)據(jù)重建方法中，稀疏變換、迭代算法和閾值模型是影響插值重建效果和計算效率的三個關鍵因素，下面分別對其進行討論。

3 稀疏變換

對于基于稀疏變換的地震數(shù)據(jù)重建方法，有多種數(shù)學變換可供選擇，如Fourier變換、Radon變換、Curvelet變換、Seislet變換、dreamlet變換和shearlet變換等。在最常用的幾種變換中： Fourier變換最易實現(xiàn)，但它主要適用于線性同相軸居多的情況，而對于彎曲同相軸居多的情況，其重建效果略差； Curvelet變換具有多尺度和多方向性[39]，可對非平穩(wěn)地震數(shù)據(jù)進行最優(yōu)局部分解，將地震數(shù)據(jù)進行稀疏表示； Seislet變換是基于小波提升算法而開發(fā)的針對地震數(shù)據(jù)的變換方法[40]，可提供比經(jīng)典小波變換更有效的地震數(shù)據(jù)壓縮能力。

圖1分析、對比了Fourier、Curvelet和Seislet三種稀疏變換的系數(shù)衰減速率： Seislet變換的系數(shù)衰減速率最快，Curvelet變換次之，F(xiàn)ourier變換最慢。一般地，稀疏變換的系數(shù)衰減速率越快，地震數(shù)據(jù)在該變換域就越可得到更稀疏的表示，也就更適用于基于稀疏變換的數(shù)據(jù)重建中?；赟eislet變換的數(shù)據(jù)重建通常能取得比Curvelet變換及Fourier變換更好的重建效果，但對于過于復雜的地震數(shù)據(jù)，比如多同相軸交叉情況，Seislet變換可能出現(xiàn)傾角求取不準確的情況而惡化重建效果。

圖1 不同稀疏變換系數(shù)衰減速率對比

圖2 三種稀疏變換重建結果對比

圖3 三種稀疏變換重建的誤差對比

對于圖2a所示的不規(guī)則缺失模擬地震數(shù)據(jù)，分別選擇Fourier、Curvelet及Seislet變換作為稀疏變換，選用POCS算法進行重建。從圖2重建后剖面可見，F(xiàn)ourier變換的重建效果明顯不如Curvelet和Seislet變換； 從圖3重建誤差剖面不難發(fā)現(xiàn)，F(xiàn)ourier變換重建誤差相對較大，而另外兩種變換重建誤差則相對較小。由于Curvelet和Seislet變換具有局部性質，對于彎曲同相軸居多的復雜地震剖面，筆者傾向于選擇Curvelet和Seislet等容易處理傾角的變換作為迭代所需的稀疏變換。但從計算效率考慮，F(xiàn)ourier變換比Curvelet和Seislet變換更具優(yōu)勢，若是進行高維地震數(shù)據(jù)重建[41]，Curvelet和Seislet變換計算量巨大，而Fourier變換則可顯著節(jié)省計算成本，此時選擇Fourier變換作為稀疏變換更合適。

4 迭代算法

本文首先對POCS算法與IHT算法進行對比分析，這兩種方法理論上具有相似性，對Sigsbee模型數(shù)值模擬結果(圖4a、圖4b)，分別用POCS算法與IHT算法來對缺失的地震道進行重建，閾值模型選取線性硬閾值模型，稀疏變換選擇Curvelet變換，迭代次數(shù)設為50次。從重建結果(圖4c～圖4f)及其誤差剖面(圖5)來看，兩種迭代算法的差異很小。為了更進一步對兩種迭代方法進行對比，計算每一次迭代的信噪比并繪制信噪比恢復曲線以更直觀地考察兩種重建迭代算法的表現(xiàn)。信噪比據(jù)下式計算[42]

圖4 缺失地震數(shù)據(jù)及POCS與IHT算法重建結果對比

圖5 POCS與IHT算法重建結果的誤差對比

(15)

式中：O為原始無噪聲完整數(shù)據(jù)；mn為第n次迭代后的重建結果。從圖6的信噪比恢復曲線對比圖可看出： 如果第n次迭代的閾值等于零，那么兩條曲線(虛線)在最后一次迭代會相交， 也就是得到相同的重建結果； 如果第n次迭代的閾值不為零，那么兩條曲線(實線)不會相交，也就是說迭代結束后會得到不同的重建結果，而且POCS算法的重建結果會略好于IHT算法。該算例表明在最后一次迭代的閾值為零的情況下，POCS算法和IHT算法會得到一樣的重建結果。如果原始數(shù)據(jù)中不含有噪聲，最后一次迭代閾值可選為零，這樣能得到較好效果； 但如果原始數(shù)據(jù)中含有噪聲，最后一次迭代閾值應該大于零，其大小與原始數(shù)據(jù)的信噪比有關。

圖6 POCS算法與IHT算法信噪比恢復曲線對比

下面對線性Bregman算法和JRSI算法進行對比。從圖7可看出，JRSI算法的重建質量要明顯好于Bregman算法。繪制上述4種迭代算法的信噪比恢復曲線(圖8)，可見POCS算法與IHT算法最為穩(wěn)定，最終達到的信噪比最高，在最后一次迭代閾值為零的情況下，二者最終達到相同的信噪比。Bregman算法雖然收斂速度比較快，但最終能達到的信噪比較低。JRSI算法則結合了Bregman算法和IHT算法的優(yōu)點，相對IHT算法來說前期可達到更高信噪比，最終的重建質量與IHT相比差異不大。

圖7 Bregman算法與JRSI算法重建結果及其誤差對比

采用圖9a所示的二維疊前實際數(shù)據(jù)來測試各個迭代算法，稀疏變換固定為Curvelet變換，閾值模型固定為線性閾值模型。重建道集以及重建誤差如圖9及圖10所示。從迭代50次之后的重建結果可以看出，Bregman迭代算法的重建結果稍差，而另外三種迭代算法都有不錯的重建效果。從圖11的信噪比恢復曲線可以看出，經(jīng)過足夠多的迭代，POCS算法、IHT算法以及JRSI算法都能達到一個較高的信噪比，JRSI算法最終達到的信噪比略高于前二者，而且用相對較少的迭代次數(shù)就可達到很高的信噪比。Bregman算法雖然收斂速率快，但是最終達到的信噪比不高。

圖8 不同迭代算法信噪比恢復曲線對比

圖9 實際缺失地震數(shù)據(jù)及其幾種迭代算法重建結果

圖10 四種迭代算法重建結果的誤差對比

5 閾值模型

Abma等[21]在應用POCS算法進行地震數(shù)據(jù)重建時指出閾值模型的選取是影響重建質量的關鍵因素，合理的閾值選擇能顯著提高重建速率和重建質量。閾值模型可分為兩個大類：硬閾值與軟閾值。

圖11 四種迭代算法信噪比恢復曲線對比

其公式分別為

(16)

(17)

其中τ為閾值，至于具體的閾值計算公式，Abma等[21]在應用POCS算法進行地震數(shù)據(jù)重建時使用線性閾值計算公式

τk=τi+(τf-τi)(k-1)/(N-1)

(18)

式中：τi與τf分別為初始閾值和終止閾值，初始閾值通常選取最大的稀疏變換系數(shù)，終止閾值一般為一個接近或者等于零的小值，其大小與原始數(shù)據(jù)的信噪比有關；k=1,2,…,N，為迭代次數(shù)，隨著迭代

次數(shù)增加，閾值從初始閾值線性降至終止閾值。

為了提高收斂速率，Galloway等[27]提出了指數(shù)閾值計算公式

(19)

與線性閾值計算相似，指數(shù)閾值計算仍然是使閾值隨迭代次數(shù)增加從初始的大閾值降至終止的小閾值，區(qū)別在于閾值按指數(shù)規(guī)律下降。

Gao等[30]提出了數(shù)據(jù)驅動閾值模型，數(shù)據(jù)驅動閾值可按照以下流程選?。?/p>

(1)計算缺失數(shù)據(jù)的頻譜，選定初始閾值與終止閾值，將大于初始閾值及小于終止閾值的振幅移除；

(2)將剩余的振幅值降序排列，得到一條閾值曲線；

(3)假設總迭代次數(shù)為N次，則將曲線等分為N-1份，同時得到N個節(jié)點；

(4)第一次迭代，將第一個節(jié)點所對應的縱坐標的值作為此次迭代的閾值，以此類推，可得到每一次迭代的閾值。

仍選取圖4中的Sigsbee模型數(shù)據(jù)進行測試，迭代算法均采用相對比較穩(wěn)定的POCS算法，分別使用線性、指數(shù)、數(shù)據(jù)驅動閾值模型進行數(shù)據(jù)重建。各閾值模型重建結果、誤差剖面以及信噪比恢復曲線如圖12、圖13及圖14所示。從信噪比恢復曲線可以看出，線性閾值模型雖然最終能達到較高的信噪比，但是收斂速度較慢，而指數(shù)閾值模型不僅具有較快的收斂速率，同時能保證很好的重建質量，數(shù)據(jù)驅動閾值模型相對于前兩種能夠更快地達到一個較高的信噪比，但是由于閾值隨迭代次數(shù)下降太快，最終能夠達到的信噪比略低，而且對于不同的數(shù)據(jù)，重建效果也可能有所不同。相對于硬閾值，軟閾值模型的重建質量則普遍較低。

圖12 幾種閾值模型重建結果

圖13 幾種閾值模型重建結果誤差

圖14 幾種閾值模型重建信噪比恢復曲線對比圖

采用圖9的實際數(shù)據(jù)測試不同閾值模型，迭代算法均采用相對比較穩(wěn)定的POCS算法，分別應用線性、指數(shù)及數(shù)據(jù)驅動模型進行數(shù)據(jù)重建，重建結果、重建誤差及信噪比恢復曲線見圖15、圖16及圖17。與模擬數(shù)據(jù)的測試結果相似，指數(shù)閾值模型無論是收斂速率還是重建質量都要優(yōu)于線性閾值模型，而數(shù)據(jù)驅動模型雖然收斂速率快，但是最終能達到的信噪比不高，也說明了其受不同數(shù)據(jù)的影響，具有不穩(wěn)定性。軟閾值模型與硬閾值模型相比，前者的重建結果稍差。

圖15 幾種閾值模型實際數(shù)據(jù)重建結果剖面

圖16 幾種閾值模型實際數(shù)據(jù)重建誤差剖面

圖17 幾種閾值模型信噪比恢復曲線對比圖

6 結論

本文對壓縮感知地震數(shù)據(jù)重建中稀疏變換、迭代算法和閾值模型三個關鍵因素進行了研究，為在實際地震數(shù)據(jù)重建中選擇稀疏變換、迭代算法以及閾值模型提供了理論依據(jù)和建議。結論如下：

(1)FFT在處理具有彎曲同相軸的非平穩(wěn)數(shù)據(jù)時具有一定局限性，而Curvelet變換和Seislet變換等因其具有多方向性和局部性質，在進行復雜地震數(shù)據(jù)重建時更具優(yōu)勢，但是Curvelet變換和Seislet變換等局部性稀疏變換計算量大，尤其對于高維地震數(shù)據(jù)重建計算成本過高，計算效率不如FFT。

(2)POCS算法與IHT算法較為穩(wěn)定，經(jīng)過足夠的迭代次數(shù)能夠達到很高的信噪比，如果最后一次迭代的閾值為零，那么經(jīng)過同樣的迭代次數(shù)，兩種算法會得到相同的重建結果； Bregman迭代算法雖然前期迭代速率快，但經(jīng)多次迭代后，并不能達到一個很高的信噪比； JRSI算法結合了Bregman算法與IHT算法的優(yōu)點，可更快地達到較高的信噪比，同時也能保證迭代結束后的重建質量。

(3)線性閾值模型最終所能達到的信噪比雖然很高，但是迭代速率較慢；指數(shù)閾值模型具有很高的收斂速率并且能保證重建質量；數(shù)據(jù)驅動模型往往具有較快的迭代速率，但是重建效果受到數(shù)據(jù)影響，穩(wěn)定性較差。

感謝吉林大學劉洋教授、中國地質大學(北京)高建軍副教授、同濟大學王本鋒博士和中國地震局白蘭淑博士的有益討論。