基于交疊組稀疏非凸Lp偽范數高階全變分的地震隨機噪聲壓制

梁上林 胡天躍* 崔 棟 隋京坤

(①北京大學地球與空間科學學院，北京 100871； ②中國石油勘探開發研究院，北京 100083)

0 引言

地震記錄中常常混雜著大量的隨機噪聲，嚴重降低了資料的信噪比。有效壓制地震數據中的隨機噪聲是地震資料處理的重要環節。然而，由于其產生具有隨機性且無固定頻率和視速度，隨機噪聲壓制一直是勘探地球物理界的難點之一。目前，壓制隨機噪聲的方法可歸納為時間域濾波[1-2]、頻率域濾波[3]、空間域濾波[4]和各種變換域濾波[5-9]等。Rudin等[4]于1992年首次提出空間域的全變分(Total variation，TV)模型，利用含噪數據總變分比無噪信號總變分大的性質，構造能量泛函以達到去噪的目的。鑒于其有效性，該模型被廣泛應用于噪聲壓制[10-12]、波阻抗反演[13]和全波形反演[14]等領域。

盡管TV模型具有保護圖像邊緣信息的優勢，但存在嚴重的階梯效應，導致信號分片光滑。針對此問題，Selesnick等[15]率先將交疊組稀疏(Overlapping group sparsity，OGS)引入TV模型，指出信號的一階差分不僅具有稀疏特性，還具有結構稀疏特性。之后Selesnick等[16]通過正則化技術有效識別出平滑區域內被強噪聲污染的信號。為挖掘圖像一階梯度和二階梯度的結構稀疏性，陳穎頻等[17]將OGS與TV模型相結合，構建了一種改進的去噪模型。

上述方法[15-17]從本質上都是在尋找不同的正則項及懲罰函數，以達到最佳去噪效果。考慮到L1正則項是一種凸近似，往往使重構的非零信號比真實值小而導致過擬合，Ahlad等[18]提出組內使用L2范數和組間使用L1范數的加權方式。Lp偽范數本質是一種非凸正則項，它具有比凸核更接近秩函數的優點，能在恢復信號細節與給定殘差稀疏解之間達到平衡，被廣泛應用于圖形修復、地震數據重建和重力反演等領域[19-21]。Adam等[22-23]基于高階TV模型，引入非凸正則項技術，在圖形去模糊方面取得一定效果； 同時指出OGS與非凸高階TV模型在去除斑塊虛假信息上是互補的。然而，目前尚未見到將Lp偽范數、OGS和TV模型三者相結合應用于地震資料去噪的實例。

鑒于一階TV模型不能有效保護信號的邊緣及紋理特征，存在嚴重的階梯效應，本文將高階TV模型與OGS相結合，目的是充分挖掘和利用兩者的優點，在減弱階梯效應的同時較好地保護局部信息； 同時，為解決信號重構過程中產生的斑點和塊狀問題，保護非零有效信號，引入非凸Lp偽范數； 最后通過模擬數據和實際資料對本文方法進行驗證，并與常規五種去噪方法進行對比。

1 去噪原理

1.1 傳統TV去噪模型

實際地震數據可表示為

y=s+n

(1)

式中：y表示檢波器接收到的地震數據；s表示有效信號；n為隨機噪聲。TV去噪模型對應如下數學極值問題[4]

(2)

1.2 二維OGS正則項

Selesnick等[15]定義了窗尺度為K×K且中心點為(i,j)的二維矩陣

(3)

m1=K-12、m2=K2,

式中表示向下取整，K表示組稀疏交疊程度。二維OGS正則項定義為

(4)

可見，OGS將鄰域信息作為參考形成組合梯度，能充分挖掘信號的局部相似性。

1.3 基于OGS非凸Lp偽范數高階TV模型

考慮到OGS能有效恢復全局較大輪廓而高階TV模型可較好地保護局部弱信號，將OGS與高階TV模型結合，充分利用兩者優點形成一個混合去噪模型； 同時，引入非凸Lp偽范數以更準確地重構噪聲數據中非零信號。該模型所對應的數學極值問題[22]表示為

(5)

改進的方法是以信號周圍點的信息作為參考形成組合梯度，使孤立的噪聲污染點與鄰域的組合梯度下降； 同時保持圖像邊緣點與鄰域的組合梯度強度，通過設定合理的閾值，可在有效去除噪聲的同時，減弱階梯效應，較好地保護有效信號。

采用分離變量法將式(5)轉化成如下受約束的極值問題

s.t.a=s,b=2s,c=s

(6)

該式對應的無約束增廣拉格朗日極值為

L(s,a,b,c;ξ1,ξ2,ξ3)=

(7)

為求解式(7)所示的復雜耦合問題，本文采用交替方向乘子法(Alternating direction method of multipliers，ADMM)[25]，按照下式交替迭代更新

(8)

可見，式(7)所示耦合問題被分解成四個子問題。下面利用式(9)～式(13)求解相應子問題。

sk+1是一個最小平方問題，其解析解為

sk+1=[λ1I+βT+β(2)T2+βI]-1[λ1y-Tξ1+

βTak-(2)Tξ2+β(2)Tbk-ξ3+βck]

(9)

ak+1的約束如下

(10)

式中包含了式(4)所示的OGS正則項，當K=1時式(10)退化成傳統的高階TV模型； 當K>1時，式(10)表示求解OGS正則項。鑒于最大最小值(Majorization minimization，MM)算法具有高效利用函數凹凸性搜尋原函數最值的優點，當原目標函數難優化時，不直接對其求解，而是求解一個易于優化且逼近于原目標函數的替代函數[26]。本文通過二次函數替代φ(a)實現MM算法。

同理，bk+1的約束為

(11)

式(11)是一個非凸優化問題。非凸優化算法和閾值的選取是決定去噪效果的兩個重要因素。鑒于經典的加權方向迭代L1算法能有效平衡信號的正負二階微分[27-28]，因此本文采用該算法求解

(12)

式中：ζ和ε都是較小量，用以避免分母為零； max表示求二者較大值； sign為符號函數。

類似地，ck+1的約束為

(13)

另外，求解式(7)所需的三個拉格朗日乘子可表示如下

(14)

至此，四個子問題已求解完成。

將整個算法流程總結如下：

(1)輸入二維地震數據y，給定參數λ1>0，λ2>0，K，p；

(2)初始化模型：s0=y，k=0，β=0.004，si=1,2,3=0，模型更新率εmin=1×10-6；

(3)更新式(9)中的sk+1；

(4)采用MM迭代算法更新式(10)中的ak+1；

(5)更新式(12)中的bk+1；

(6)更新式(13)中的ck+1；

2 模擬數據驗證

2.1 評價指標

信噪比通常用信號總能量與噪聲總能量的比值SNR(Signal to noise ratio)表征，是評價整體去噪效果的常用指標；結構相似性參數SSIM(Structural similarity)從亮度、對比度、結構三方面度量兩幅圖像的相似性，常作為衡量圖像失真或噪聲壓制情況的客觀標準[29]； 同時，為表征高斯噪聲強弱，定義零均值噪聲的標準差為δ。三者的具體定義如下

(15)

(16)

(17)

式中：M、N分別表示二維地震數據的行數和列數；μs、μy對應表示s、y的均值；σs、σy分別表示s、y的方差；σsy表示s與y的協方差；C1和C2是維持穩定的兩個常量，一般分別取2.252和6.752。為了更全面地評判去噪效果，本文將SNR、SSIM和運算耗時作為三個客觀評價指標。

2.2 模擬數據

構建一個水平層狀模型，通過有限差分聲波正演模擬得到圖1所示的共炮記錄。模擬采用主頻為20Hz的雷克子波，采樣間隔為4ms，采樣長度為4s，接收道數為500。為了模擬地震數據中不同強度的隨機噪聲，向該記錄加入標準差為δ的高斯白噪聲，得到圖2所示的含噪數據。由圖2可知，隨著δ增大，背景噪聲加強，SNR降低，SSIM變差，弱反射信號逐漸被強噪聲淹沒。

圖1 正演模擬共炮記錄

圖2 不同δ的含噪數據(a)δ=10； (b)δ=20； (c)δ=30； (d)δ=40

2.3 參數優選

OGS從一個點向四周延伸K×K個點，可有效挖掘信號的局部相似性。為考察K值對去噪效果的影響，通過不斷調整其大小，分別對圖2所示的四種不同噪聲背景下的地震數據做去噪處理。

圖3展示不同K值得到的去噪結果所對應的三種指標的變化曲線。從圖3a可見，隨著K值的增大，SNR先從小到大逐漸增加，達到最大值后減小。當SNR取最大值時，K值分別為3、5、7和10，說明K值的優選與背景噪聲強度有關。當背景噪聲較弱時，較小K值就能達到滿意的去噪效果；當背景噪聲增強時，應適當增大K值以提高SNR。圖3b中的SSIM表現出與SNR類似形態，區別在于當SSIM達到最大值后緩慢減小。這說明K值超過最優值后，對SSIM的影響不顯著。圖3c顯示，隨著K值的增大運算耗時也不斷增加，這是由于鄰域信息的引入，導致計算量增大。

圖3 去噪結果的三指標隨K值的變化曲線(a)SNR； (b)SSIM； (c)耗時

因此，K值對本文方法去噪效果有顯著影響。當K值過小時，鄰域信息不足導致SNR和SSIM達不到最優，去噪效果欠佳； 而當K值過大時，又會引入過多不相似的數據點，反而導致評價指標減小，去噪能力下降，這不僅會產生平滑效應，還會因使用過多鄰域信息增加不必要的計算量。

邊界與局部細節的恢復主要受非凸正則化階數p控制，對圖2含噪數據進行處理并得到圖4所示三種指標隨p的變化曲線。圖4a表明，當δ=10時，SNR隨p值增大而快速增至最大值，然后快速減小；當δ=20、30和40時，SNR曲線都具有平緩段、快速上升段和快速下降段。在平緩段SNR對p值不敏感，經快速上升后達到最大值。SNR取最大值時對應p值分別為0.23、0.46、0.62和0.75，這說明p值的優選與噪聲強度有關。圖4b表明，弱噪聲情況下SSIM曲線先上升，然后保持平緩；當背景噪聲變強時，SSIM曲線呈現平緩段—快速上升段—平緩段分布。該指標取最大值時對應p值分別為0.16、0.47、0.62和0.83。圖4c表明，隨著p值增大，運算耗時先增至最大值，然后逐漸減小，顯然該指標與背景噪聲強度無關。

圖4 去噪結果的三種指標隨p值的變化曲線(a)SNR； (b)SSIM； (c)耗時

綜上，p值對本文去噪方法有同樣顯著的影響。高SNR資料應取較小p值；反之，低SNR資料應適當增大p值，以達到最佳去噪效果。

圖5展示K與p最優情形下對圖2所示四種不同強度噪聲背景下地震數據的去噪效果。不同強度隨機噪聲均得到有效壓制，剖面整體干凈，反射波同相軸清晰，深層弱能量信號得到有效保護。

圖5 含有不同δ背景噪聲數據的去噪結果(a)δ=10； (b)δ=20； (c)δ=30； (d)δ=40

以上去噪試驗表明，本文方法能壓制階梯效應，去除不同強度隨機噪聲，提高信噪比，并能有效保護弱能量信號。

2.4 方法對比

為進一步驗證本文方法的有效性，分別與各向異性全變分(Anisotropic total variation，ATV)、交疊組稀疏全變分(Overlapping group sparsity total variation，OGSTV)、中值濾波(Median filter，MF)、三維塊匹配(Block matching 3D，BM3D)和K奇異值分解(K-singular value decomposition，KSVD)等五種方法進行對比。其中，ATV和OGSTV與本文方法屬于同類，后三者為常用去噪方法。

表1展示這些方法去噪后各項指標對比結果。從SNR和SSIM指標看，本文方法明顯優于其他五種方法。從運算耗時指標看，本文方法用時雖然大于MF、ATV和OGSTV，但都未超過8s，在可接受范圍內。對比后發現，用時明顯小于BM3D和KSVD，這是由于最后兩種方法涉及到耗時的非局部塊匹配和字典訓練。

表1 針對四種含噪數據的六種方法去噪結果評價指標對比

圖6展示δ=40時上述六種方法去噪結果。強噪聲背景下，ATV(圖6a)去噪并不理想，殘留明顯斑點和小塊，階梯效應嚴重。OGSTV(圖6b)引入鄰域信息，階梯效應得到一定壓制，但仍存部分斑點偽影。MF(圖6d)去噪效果最差，存在大量殘余噪聲。BM3D(圖6e)產生了平滑效應，對深層弱信號的保護能力差。KSVD(圖6f)存在一定斑點，影響了最終去噪效果。圖6c所示剖面干凈，階梯效應弱，深層反射波同相軸連續性好，說明本文方法具有保護弱信號能力。

圖6 不同去噪方法對比(a)ATV； (b)OGSTV； (c)本文方法； (d)MF； (e)BM3D； (f)KSVD

從上述不同方法縱橫向對比可知，本文方法能有效去除隨機噪聲，壓制階梯效應，更好地保護弱信號，在效率與精度上能達到良好的平衡。

3 實際資料應用

將本文去噪方法應用于M工區實際地震資料。所選單炮記錄(圖7a)共有240道，采樣點數為700，采樣間隔為2ms。由于存在大量隨機噪聲，反射波同相軸連續性差，弱能量信號被掩蓋，信噪比較低。其疊后剖面(圖7b)共有800道，采樣點數為800，采樣間隔為2ms。該剖面深部地層有一定起伏，且發育微斷層。隨機噪聲的存在致使剖面不清晰，同相軸錯斷，給地質解釋帶來諸多不便。

分別用上述六種方法進行去噪處理并得到圖8(單炮)和圖9(剖面)所示結果。歷經多次試驗并優選，針對圖7所示實際資料，本文方法參數設置為K=4、p=0.18和K=3、p=0.21。對比所得單炮記錄去噪結果可知，ATV(圖8a)和MF(圖8d)去噪效果差，信號失真嚴重，尤其是1.2～1.5s和2.1～2.5s區間的弱信號； 相對而言，OGSTV(圖8b)階梯效應明顯減少，去噪效果有一定提升，但仍存殘余噪聲； BM3D(圖8e)能保留強能量信號，但弱信號恢復能力差，產生了嚴重的平滑效應； KSVD(圖8f)和本文方法(圖8c)去噪效果相對較好，但KSVD在一些細節(弱信號)上不如本文方法。

圖7 原始地震資料(a)單炮資料； (b)疊后剖面

圖8 不同單炮記錄去噪結果對比(a)ATV； (b)OGSTV； (c)本文方法； (d)MF； (e)BM3D； (f)KSVD

對比、分析圖9及表2可知，MF去噪效率高，但剖面(圖9d)上仍有大量噪聲殘留； OGSTV(圖9b)去噪效果有改善，但仍不理想； ATV(圖9a)去噪后產生的階梯效應使部分細節丟失； 同樣地，BM3D(圖9e)出現平滑模糊區塊，不利于微斷裂等信息的恢復； KSVD(圖9f)和本文方法(圖9c)都有效壓制了隨機噪聲，但KSVD計算效率低。

表2 不同去噪方法的計算耗時對比

圖9 疊后資料去噪結果對比(a)ATV； (b)OGSTV； (c)本文方法； (d)MF； (e)BM3D； (f)KSVD

針對實際資料的應用結果表明，本文方法能壓制不同強度隨機噪聲，減弱階梯效應，有效保護信號細節特征，去噪后整個剖面同相軸具有更好清晰度和連續性，因此具有良好應用潛力。

4 結論

(1)本文去噪方法適用于不同強度的背景噪聲，在有效壓制隨機噪聲的同時兼顧計算效率，能顯著提高地震資料的品質，具有廣闊的應用前景。

(2)OGS考慮了信號的鄰域信息，能充分挖掘局部相似性； 將它與高階TV模型和Lp偽范數結合，不僅能壓制階梯效應，提高算法去噪能力，而且能有效保護圖像邊緣信息，具有較高保真度。

(3)K值決定鄰域信息的多少，若K值過小，則不能充分挖掘鄰域信息；反之，若引入非相似信息，則不僅增大了計算量，而且導致平滑效應。p值關系到局部非零值信號的恢復，對于弱信號局部細節的保護有重要作用。