基于交疊組稀疏廣義全變分的地震信號隨機噪聲衰減

陳穎頻 彭真明 李美惠 喻 飛

(①電子科技大學信息與通信工程學院，四川成都 610054；②閩南師范大學物理與信息工程學院，福建漳州 363000；③電子科技大學信息地學研究中心，四川成都 610054)

0 引言

實際地震信號中存在大量隨機噪聲，這給地震資料的處理和解釋帶來較大的干擾。目前，去除地震信號隨機噪聲的方法主要有基于全變分模型的正則化去噪方法[1-4]、基于稀疏表示和字典學習的隨機噪聲去除方法[5-7]和基于變換域的濾波去噪方法[8-9]等。其中，基于稀疏表示和字典學習的去噪方法需觀察數據進行字典訓練，在獲得較好去噪效果的同時，要耗費大量時間在字典學習上，只要數據量稍大就會引起較大計算量。基于變換域的濾波去噪方法通常需將地震信號從空域映射到變換域中進行截斷濾波，易產生吉布斯效應。

全變分(Total variation，TV)模型被證明是一種可有效去除隨機噪聲的模型，自從它被Rudin等[1]提出之后，就引起學者們的廣泛關注[10-16]。TV模型充分挖掘了二維圖像的橫向縱向梯度信息，較好地契合了自然圖像的局部光滑和梯度稀疏等先驗知識。該模型被廣泛應用于地震圖像去噪[10-11]、地震波阻抗反演[12-13]、地震波全波形反演[14]、弱小目標檢測[15]和超分辨率分析[16]等眾多領域。

由于TV模型假設圖像是分片光滑常數，導致TV去噪模型存在較為嚴重的階梯效應。為抑制階梯效應和提高TV模型的去噪效果，Bredies等[17]提出廣義全變分(Total generalized variation，TGV)模型。該模型是全變分模型的推廣，具有凸性、下半連續性、旋轉不變性等眾多優秀的數學性質，并能逼近任意多項式[18]。TGV模型同時約束了圖像的一階梯度與二階梯度，因此有效緩解了全變分模型的階梯效應。Knoll等[18]將廣義全變分模型用于核磁共振成像，取得較好的應用效果。隨后，Guo等[19]提出一種基于TGV和剪切波變換的細節保留方法，用于MRI圖像重構中。Qin等[20]將這種多約束正則項用于圖像解卷積。Kong等[11]則將Shearlet和TGV正則項用于地震信號去噪。

近年來，Liu等[21]和Chen等[22-23]提出了交疊組稀疏正則項。該正則項是一種非分離正則項[24]。它不僅考慮了圖像差分域的稀疏性，還挖掘了每個點鄰域差分信息，因此可利用圖像梯度的結構化稀疏特性。通過交疊組合梯度可提高平滑區域與邊界區域的差異，從而抑制TV模型的階梯效應。Liu等[25]借鑒Chen等[22-23]的成果，將一維交疊組稀疏正則項推廣到二維領域，并將其引入各向異性全變分模型，用于椒鹽噪聲的去噪和解卷積問題。Liu等[26]將交疊組稀疏正則項用于Speckle噪聲的去除。

上述工作全部是基于各向異性全變分模型做的推廣，然而各向異性全變分模型沒有二階差分的約束，交疊組稀疏梯度的引入對階梯效應的抑制能力有限。需指出的是，TGV和交疊組稀疏收斂算子對階梯效應的抑制機理并不一樣，前者利用圖像一階、二階差分約束緩解階梯效應，后者則通過圖像梯度的結構特性抑制階梯效應。目前，TGV模型與交疊組稀疏收斂算子的交叉研究仍處于起步階段。

由于經典的TGV模型未考慮圖像梯度的鄰域結構特性，受Liu等[25]的啟發，本文提出一種基于交疊組稀疏收斂算子的改進廣義全變分模型(簡稱為TGV-OGS)，并將其應用于地震信號去噪中。鑒于該模型的復雜性，采用交替乘子迭代法(Alternating direction method of multipliers，ADMM)[27-28]對其求解。為提高算法效率，假定處理圖像滿足周期性邊界條件，并采用快速解卷積方法[29-31]，巧妙地避免了大型差分矩陣的相乘運算，將差分操作理解為卷積，再利用卷積定理，從頻域恢復圖像。

在后續實驗中，將對比常規各向異性全變分方法(Anisotropic total variation，ATV)[32]、各向異性組稀疏方法(Anisotropic total variation based on overlapping group sparsity，ATV-OGS)[25]、TGV去噪方法及本文去噪方法，并從峰值信噪比(Peak signal to noise ratio，PSNR)、結構相似性(Structural similarity，SSIM)[33]、信噪比(Signal to noise ratio，SNR)[7]及運算時間等指標客觀地對比這幾種算法。從實驗結果可見，本文提出的算法及所構建的模型能進一步改進TGV的去噪性能。

1 預備知識

經典的二階TGV模型[34]定義為

‖Dv*S-Vv‖1)+α1(‖Dh*Vh‖1+

‖Dv*Vv‖1+‖Dv*Vh+Dh*Vv‖1)]}

(1)

1.2 二維矩陣交疊組稀疏鄰近算子與MM算法

假定V0為待收斂矩陣，則其交疊組稀疏鄰近算子記為

P(V) =proxγφ(V0)

(2)

(3)

圖1 二維交疊組稀疏示意圖

利用優化—最小化方法(Majorization-minimization，MM)[36]可有效計算式(2)。根據MM算法，最小化P(V)需首先找到一個函數Q(V,U)(U為輔助變量)，該函數對所有V,U，都有Q(V,U)≥P(V)，并且當且僅當U=V時等號成立。據此，每次計算的Q(V,U)最小值為P(V)的優化值，式(2)的計算可轉化為如下迭代運算

(4)

考慮到組稀疏正則項的特殊性，并注意到下列不等式的存在

(5)

式中，當且僅當U=V時，等號成立。

(6)

將S(V,U)重寫為

(7)

式中：v是矩陣V的向量形式；C(U)與V無關，可視為關于V的常數項；D(U)∈RN2×N2，是一個對角矩陣，其對角元素定義為

[D(U)]m,m

(8)

結合式(4)、式(7)和式(8)，可將式(2)轉化為如下迭代優化問題

(9)

該式為二次規劃問題，其迭代最優解為

V(k+1)=mat{[I+γD2(V(k))]-1v0}

(10)

式中：I∈RN2×N2，表示單位矩陣；v0是V0的向量形式； mat表示向量矩陣化運算。

2 模型構建及其求解方法

2.1 基于交疊組稀疏的改進TGV模型

將TGV中的L1約束項改進為交疊組稀疏約束項，從而更好地挖掘圖像一階梯度與二階梯度的差分信息，建模如下

φ(Dv*S-Vv)]+α1[φ(Dh*Vh)+

φ(Dv*Vv)+φ(Dv*Vh+Dh*Vv)]}}

(11)

對比經典的TGV模型，本文構建的模型將一階梯度與二階梯度矩陣的每個像素點鄰域K2個信息點交疊組合，從而更充分挖掘一階梯度矩陣和二階梯度矩陣結構特性。顯然，每個像素點的鄰域梯度與二階鄰域梯度存在一定相似性，該結構稀疏先驗知識被交疊組稀疏模型較好地刻畫。

2.2 構建模型的ADMM求解

為求解式(11)定義的改進TGV模型，利用ADMM框架求解模型。令

則根據ADMM原理，可將式(11)的增廣拉格朗日函數表達為

μ0[φ(Z1)+φ(Z2)]+μ1[φ(Z3)+

φ(Z4)+φ(Z5)]-

〈β0Λ1,Z1-(Dh*S-Vh)〉-

〈β0Λ2,Z2-(Dv*S-Vv)〉-

〈β1Λ3,Z3-(Dh*Vh〉-〈β1Λ4,Z4-Dv*Vv〉-

〈β1Λ5,Z5-(Dv*Vh+Dh*Vv)〉+

(12)

式中：μ0=μα0；μ1=μα1；Λi(i=1,2,…,5)是拉格朗日乘子； 〈X,Y〉表示兩個矩陣X與Y的內積。

在ADMM框架下，分離變量及其對偶變量之間與三元組S、Vh、Vv是去耦合的。而三元組中三個變量是相互耦合的，因此需要建立關于三個變量的三元一次方程組。對于S子問題，其子目標函數為

(13)

本文假定處理數據滿足周期性邊界條件，因此可以利用快速傅里葉變換進行卷積計算。

根據卷積定律，兩個矩陣在空域卷積，對應到頻域，卷積結果的頻譜為兩個矩陣頻譜的點乘。將式(13)做傅里葉變換

(14)

(15)

整理可得

(16)

其中

式中：1是元素全為1的矩陣；上角“*”表示共軛運算。

同理，對Vh子目標函數關于Vh求導為零，可得Vh子問題的約束方程為

(17)

其中

詳細推導過程不再贅述。

同理，對Vv的子目標函數關于Vv求導為零，可得Vv子問題的約束方程為

(18)

其中

(19)

該式可用克萊姆法則和快速反傅里葉變換求解

(20)

式中：除法都為按元素點除運算；F-1表示二維逆快速傅里葉變換算子；定義矩陣行列式算子為

K12。K23。K31+K13。K21。K32-K13。K22。K31-

K12。K21。K33-K11。K32。K23

對于Z1子問題，其目標子函數為

(21)

根據式(10)，可得Z1的更新公式為

(22)

同理，Z2、Z3、Z4和Z5的更新公式為

(23)

且有

式中vec(·)表示將矩陣列向量化。

Λ1的目標子函數為

(24)

利用梯度上升法可得其更新公式

(25)

其中γ為學習率。

類似地，拉格朗日乘子變量Λ2、Λ3、Λ4和Λ5的更新公式對應為

(26)

至此，本文構建模型的所有子問題都得以解決。將該算法總結為表1。

表1 TGV-OGS算法偽代碼

其中： tol=10-4表示算法停止迭代的閾值；Max=30

3 數值實驗

將本文算法應用于地震信號，并通過峰值信噪比、結構相似性指示系數、信噪比及運行時間客觀評價該算法的去噪性能，并與ATV、ATV-OGS方法以及TGV去噪方法做全面對比。

3.1 評價指標

去噪領域最常用評價指標有：峰值信噪比、結構相似性指標、信噪比和運行時間。其中峰值信噪比、結構相似性指標[33]為

(27)

SSIM(X,Y)

(28)

為更好地表征去噪效果，還選用了信噪比參數，定義[7]如下

(29)

3.2 合成地震信號去噪性能對比測試

本節給出合成地震記錄，該記錄由20 Hz雷克子波與人工合成反射系數卷積獲得。圖2為截取的前50道合成記錄，可清晰觀察到去噪效果。

圖2 合成地震信號

3.2.1 交疊組稀疏正則項效果測試

為驗證交疊組稀疏正則項的有效性，在合成記錄中加入0均值、標準差為σ的高斯白噪聲，且設定K=3。ATV、ATV-OGS、TGV以及本文方法的測試結果如表2所示，最優指標用粗體標出。

表2 合成地震信號算法性能測試表

從表2中可看到，交疊組稀疏正則項不僅對經典的各向異性TV模型算法性能有所改進，而且對TGV去噪方法有較好的改進。為直觀反映其效果，以偽彩色形式(圖3)展示σ=30時上述四種算法的去噪效果。

觀察圖3中由矩形紅色框圈定區域(圖3d)，顯然是平滑區域，在放大區域中可以看到若干個被重噪聲污染的像素；ATV去噪方法在該區域有較明顯的塊效應。從放大圖(圖3f)中可見，仍然有兩個像素被重噪聲干擾。從圖3h可見，利用交疊組稀疏正則項則非常好地抑制了塊效應，且兩個重噪聲點被去除得較為干凈。值得指出的是，各向異性交疊組稀疏方法并未完全抑制ATV的塊效應。再觀察TGV的去噪結果(圖3i、圖3j)，可見TGV去噪方法也對ATV的塊效應有較好的抑制作用，在平滑區域獲得更加光滑的重構信號。但從圖3j可見，TGV去噪方法對于重噪聲污染點也是無能為力的。而本文提出的方法則將交疊組稀疏和TGV的優點結合起來，獲得比ATV-OGS和TGV更好的去噪效果。觀察圖3l可以看到，本文算法獲得的去噪結果在該區域更接近原圖。值得注意的是，設定K=3時，算法性能尚未達到最優。

為進一步觀察去噪效果，將圖3中六個子圖的第40道抽出進行觀察(圖4)。從圖4可見ATV存在明顯的階梯效應，而ATV-OGS去噪方法(圖4d)較好地緩解了ATV的階梯效應，但是去噪結果仍然不夠光滑。傳統TGV方法也存在類似問題。本文提出的方法在TGV理論框架上結合了組稀疏收斂方法，從而挖掘了數據鄰域梯度的結構信息，將兩者的優點充分結合，大幅提高了信號的重構質量。在圖4b中50～100ms區間，地震信號被一強幅度噪聲嚴重干擾。該強噪聲污染在圖4c和圖4e中依然存在，這反映傳統ATV和TGV方法的重大缺陷，即兩種算法都是以單一像素點為處理對象，孤立地迭代，這樣就無法充分考慮到信號鄰域的梯度相似信息。而ATV-OGS方法和本文提出方法則不存在這種問題，觀察圖4d和圖4f可知，基于交疊組稀疏收斂方法的ATV-OGS和TGV-OGS都能有效去除50～100ms存在的尖峰干擾。顯然，本文提出方法充分利用了地震信號的鄰域相似性，這對去除大幅度的噪聲污染尤其有效。通過上述實驗，得出如下結論，各向異性交疊組稀疏去噪方法能有效利用圖像一階梯度的結構特性，而廣義全變分模型中存在一階和二階圖像梯度，這些梯度同樣含有鄰似性，將交疊組稀疏正則項與TGV模型充分結合，進一步提升了TGV模型的去噪能力。

圖3 四種算法去噪效果對比

(a)原圖；(b)圖a的局部放大；(c)帶噪聲圖；(d)圖c的局部放大；(e)ATV去噪結果；(f)圖e的局部放大；(g)ATV-OGS去噪

結果；(h)圖g的局部放大圖；(i)TGV去噪結果；(j)圖i的局部放大圖；(k)TGV-OGV(K=3)去噪結果；(l)圖k的局部放大

3.2.2 參數敏感性分析

對本文算法的交疊組合數K進行測試對比，以評價其對算法的整體效果。采用PSNR和SSIM兩個指標對算法進行客觀評價，針對不同高斯噪聲水平，將K從1到13連續變化，并記錄PSNR和SSIM值，詳見圖5。從圖5可以看到，隨著K的增大，在不同噪聲下，K對PSNR和SSIM起到的作用是不一樣，顯然，在低噪聲的時候(例如σ=10)，K取11時，PSNR和SSIM都達到峰值，若K繼續加大，PSNR則下降。當σ=10時，本文算法的PSNR達到最高值37.7846 dB，高出TGV算法1.905 dB。從圖5可見，當噪聲較低的時候，圖像的鄰域信息對算法性能起到了正面作用。然而，K也不宜取得過大，否則可能會取到鄰域中圖像變化劇烈的區域，這種情況下，PSNR和SSIM就要下降。而當噪聲較大時，鄰域梯度的結構特性被破壞得比較嚴重，這種情況下，K稍大就會導致算法性能下降，例如當σ=40時，K=7時獲得了PSNR和SSIM的最高值。從圖5a中可見，當標準差為30時，設定K=9可以令PSNR值達到最高。圖6展示了在標準差為30時，

圖4 四種算法單道去噪效果對比

圖5 不同σ值的PSNR(a)和SSIM(b)隨K值變化曲線的對比

圖6 TGV-OGV(K=9)去噪效果圖

算法參數K=9的去噪效果，其PSNR為30.4824dB，SSIM為0.9789，SNR為15.8382dB，均遠超過經典TGV去噪模型。對比圖6b與圖4l可見，當K設置得足夠合理，恢復出來的地震信號將更加符合梯度的稀疏先驗知識。再對比圖6c與圖4f可知，當K取值合適時，本文算法的保護邊緣能力及對重噪聲的抗噪能力更強。

3.3 實際地震信號去噪性能測試

以二維地震信號作為測試信號(圖7a)。該地震信號為疊后地震信號，已經去除了噪聲。為客觀評價算法性能，在地震信號中加入σ=10，20，30，40的高斯白噪聲，對比ATV、ATV-OGS、TGV以及TGV-OGS等算法的去噪性能指標，各項指標被記錄在表3中。表3顯示，對于實際地震信號，本文提出算法的去噪效果最好。

圖7展示了σ=30時各種算法的二維地震信號去噪效果，可見本文算法有效去除了人為增加的高斯噪聲。對比原地震信號可以發現，本文方法將地震信號中的高頻干擾去除得較徹底，不存在ATV方法的塊效應問題，且同相軸具有較好的橫向連續性。

圖8展示了圖7中各個子圖中道號為70的單道地震信號。從圖8b可以觀察到，在1850～1900ms區間內，地震信號被一個大幅度噪聲污染，顯然，基于逐點迭代的ATV和TGV方法都無法有效去除

該位置的噪聲(圖8c和圖8e)。而基于交疊組稀疏收斂方法的ATV-OGS和TGV-OGS方法則非常好地去除了該位置的噪聲(圖8d和圖8f)。對比圖8d與圖8f可見，本文提出方法恢復的地震信號更接近原地震信號。從圖8還可知，交疊組稀疏正則項對于平滑區域的高噪聲污染有較為顯著的效果。

表3 實際地震信號算法性能測試表

圖7 四種算法去噪效果對比(a)二維地震信號；(b)帶噪數據；(c)ATV去噪結果；(d)ATV-OGS去噪結果；(e)TGV去噪結果；(f)TGV-OGV去噪結果

圖9給出圖8各個子圖的頻譜，可見加噪地震信號存在大量高頻干擾，而經過各種去噪方法去噪后的地震信號頻譜則一定程度壓制了高頻干擾。從圖9f可見，在大于100Hz范圍，本文方法壓制噪聲的效果明顯好于ATV、ATV-OGS、TGV方法。

圖10顯示圖7b中加入的噪聲(σ=30)和各種方法去除的噪聲剖面的對比。觀察圖10a的色標可發現，實際加入噪聲的幅度范圍是[-120,120]；而從圖10b可知，ATV方法估計的噪聲范圍是[-60,60]，顯然低估了噪聲幅度；從圖10c可見，ATV-OGS方法估計的噪聲范圍是[-70,70]，結果更準確些；觀察圖10d可知，TGV方法估計的噪聲范圍是[-50,50]；而本文提出方法(圖10e)估計的噪聲范圍是[-100,100]，最接近實際加入噪聲的幅度。因此，從噪聲分布來看，本文方法估計的噪聲范圍是最接近實際添加噪聲的。

圖8 四種算法單道去噪效果對比

圖9 圖8數據對應的頻譜

圖10 四種算法去除噪聲的對比

4 結論

本文從交疊組稀疏正則項出發，結合廣義全變分定義，在ADMM框架下提出一種改進廣義全變分去噪算法，并將其應用于地震信號去噪。該算法充分利用圖像一階、二階梯度的鄰域相似性，提高平滑區域與邊界區域的差異性，從而提高去噪的魯棒性，獲得相比于經典TGV更好的去噪性能。

將本文方法與ATV、ATV-OGS、TGV算法進行比較，從實驗結果得出如下結論和認識：

(1)將差分算子視為卷積算子，并結合ADMM能將復雜的模型轉換為一系列簡單的數學問題，并估計出較為準確的地震信號。

(2)本文算法去噪能力高于其他各類全變分去噪算法，尤其對平滑區域的大幅度噪聲污染有顯著的去噪效果，因此，提出算法對噪聲更加魯棒。

(3)取適當的K值可以有效提升整體去噪性能，在實際應用中需要調節參數取值，若K值取得過小，則鄰域信息利用不夠充分；反之，若取值過大，則會引入過多不相似的圖像塊，導致去噪性能下降。

值得注意的是，本文方法提出的通用正則項同樣適用于其他各類圖像重構問題，如地震波阻抗反演、自然圖像去噪、椒鹽噪聲去除、圖像解卷積，核磁共振圖像重構等反問題的求解中，針對此于問題有待于深入研究。另外，由于引入MM算法，導致涉及兩重循環，運算效率較低，因此算法的效率優化無疑是未來研究熱點。