孫成禹 刁俊才 李文靜

(①中國石油大學(華東)地球科學與技術學院，山東青島 266580； ②青島海洋國家實驗室海洋礦產資源評價與探測技術功能實驗室，山東青島 266071； ③東方地球物理公司物探技術研究中心，河北涿州 072751)

0 引言

地震數據中噪聲嚴重影響資料處理和解釋結果[1]，壓制噪聲對后期地震資料的處理、反演和解釋等具重要作用[2]。地震資料常用的去噪方法一般分為空間域方法和變換域方法兩大類?？臻g域去噪方法包括均值濾波、中值濾波[3]和維納濾波[4]等，變換域去噪方法包括傅里葉變換、小波變換、曲波變換和Randon變換等閾值方法[5-6]。這些去噪算法本質都是利用了圖像本身的局部相關性，沒有充分挖掘、利用圖像的非局部相關性，因此在壓制噪聲的同時也會損害部分有效信號，不能精確保持斷層、裂縫等地質體的邊緣特征，甚至會出現偽影等現象[7]。

針對這一問題，Buades等[8]提出利用原始數據的結構相似性進行非局部去噪，與傳統(tǒng)的局部去噪算法相比，該算法能有效地保護圖像的邊緣細節(jié)信息。在非局部去噪算法中，Kostadin等[9]提出的三維塊匹配(Block-Matching 3D，BM3D)去噪算法利用塊匹配和三維變換域濾波技術進行串聯去噪，該方法綜合了非局部算法和變換域算法的優(yōu)勢[10-11]，對隨機噪聲具有很好的去噪效果。BM3D去噪是將圖像分割成許多個不同的小塊，根據不同小塊之間的相似性進行塊匹配，然后在三維變換域中去除噪聲，充分利用了數據的自相似性和冗余性信息，能較好地保留信號細節(jié)[12]。韓玉蘭等[13]針對BM3D算法運算量大、運行時間過長的問題，采用積分圖計算塊的相似性，代替相應的濾波過程，提高了計算效率；李文靜等[14]將改進后的BM3D算法引入地震數據處理領域，得到了較為理想的去噪效果。但該方法依然存在一些不足，濾波閾值參數需要根據經驗人為設定[15]，而閾值選取的不確定性嚴重影響基礎估計部分的降噪效果。由于在實際計算中，噪聲方差對閾值大小和塊匹配相似度的影響最大，而實際數據噪聲的方差卻是未知的[16]，極大地限制了該算法在地震資料去噪中的應用。

在變換域地震資料去噪方法中，Candès等[17]在脊波變換的基礎上提出了第一代曲波變換。張恒磊等[18]利用曲波尺度分解得到較為準確的噪聲方差估計值，使后續(xù)的曲波變換閾值選取更加精確。但是由于實際地震數據的復雜性，單一曲波變換方法的去噪精度不足，處理結果常常存在偽影和過度平滑等現象。針對這些問題，薛詩桂等[19]將曲波變換和循環(huán)平移技術結合，消除了曲波變換產生的偽吉布斯效應。姚振岸等[20]將各向異性擴散濾波和曲波變換結合保護地震數據的邊界特征；楊會等[21]將曲波變換和二維經驗模態(tài)分解結合，在去除噪聲的同時較好地保護了弱信號。曹靜杰等[22]將曲波變換作為稀疏變換自適應地進行稀疏反演去噪。

為了彌補曲波變換和常規(guī)BM3D去噪方法的缺陷，本文結合兩種方法的優(yōu)勢，提出了一種基于曲波噪聲估計的BM3D地震資料去噪方法。該方法無需人為設定相關參數，根據曲波變換得到的地震資料噪聲估計作為先驗信息，準確計算塊匹配相似度，自適應選取合適的濾波閾值進行去噪，可以得到更好的效果。數值模型和實際資料測試表明，與常規(guī)BM3D去噪方法、曲波變換去噪法相比，該方法去噪效果明顯、實用性較高。

1 基本原理

1.1 三維塊匹配去噪原理

BM3D去噪是圖像處理中的一種方法。其基本算法是：將含噪數據看作一個圖像，將其劃分為若干塊，并根據各塊與其他塊的相似度進行匹配分組；對組內數據進行三維線性變換，在變換域內對噪聲進行濾波后反變換回原來的域內，并將每個組內的塊像素返回到圖像原來的位置；最后對重疊的塊進行加權平均，得到去噪后的結果。為了提高匹配分組的正確性，通常對含噪數據進行前置硬約束閾值濾波作為地震信號的基礎估計，獲取維納濾波的收縮系數。因此，整個過程就包含基礎估計和重新估計兩個階段。每個階段也都包含塊匹配、變換域去噪和聚集三個環(huán)節(jié)。

在基礎估計的塊匹配階段，將含噪地震數據視為二維圖像，劃分成尺寸為N×N的多個塊Zi(i=1，2，…，n，n為塊總數)，依次選取各塊作為參考塊Zr(r=1，2，…，n)，定義塊與參考塊之間的距離為

(1)

(2)

(3)

(4)

(5)

(6)

將維納濾波應用于噪聲數據矩陣Yr，即將噪聲數據矩陣三維變換后的系數與維納濾波收縮系數相乘，再進行反變換后得到新的估計值為

(7)

最后將不同群組中的重疊塊進行加權平均得到重新估計值

(8)

1.2 基于曲波噪聲估計的BM3D地震資料去噪方法

常規(guī)的BM3D去噪方法需要根據噪聲方差σ2作為先驗信息求取濾波閾值和塊匹配參數。在實際地震資料的處理中，常用的噪聲方差估計方法是中值噪聲方差估計和基于小波分析的細節(jié)估計。當原始地震資料信噪比較高、噪聲方差較小時，這兩種方法的估計值與真值較為接近。但是，當原始地震資料的噪聲水平較高時，前者的估計值與真值偏差較大，而后者則不能完整地估計噪聲，導致估計值比實際噪聲方差偏小。

(9)

根據σ可以進一步確定BM3D去噪算法中的濾波閾值γ1及收縮系數φr。這樣濾波參數成為噪聲強度σ的自適應函數，弱反射和邊界細節(jié)信息可以得到有效的保護。

BM3D去噪算法結合了非局部算法和變換域算法的優(yōu)勢，可以有效去除高斯隨機噪聲。由于事先無法得到噪聲方差，需要人為設定濾波參數，因而去噪效果難以保證。本文將曲波變換估計出的噪聲方差作為先驗信息，求取適用于BM3D去噪算法的輸入噪聲強度σ，進而自適應地調整濾波閾值和塊匹配參數，完成去噪處理。當σ較小時，可直接應用式(1)計算塊匹配相似度；當σ為中等或較大時，直接塊匹配可能會導致錯誤的塊匹配分組，需要先對圖像塊進行預濾波處理，即先對兩個圖像塊進行正則二維線性變換，對其進行閾值濾波，然后再計算兩者之間的距離，即

(10)

式中：γ0=λ2Dσ為二維預濾波處理的閾值(λ2D為二維系數)；T2D為正則二維線性變換算子。去噪閾值參數通常都設為常量，為了有效地保護有效信號，將其設定為噪聲強度σ的自適應函數。當σ較小時，無需進行預濾波處理，圖像塊可以取值較小(如N取為8)； 當σ中等時，需要進行預濾波處理，濾波閾值γ0取為1.0，圖像塊應適當增大(如N取為9～10)；當原始資料噪聲強度σ較大時，二維預濾波閾值γ0可取為2.0， 圖像塊應取值較大(如N取為11～12)。

2 模型試算

為了測試自適應選取閾值的改進方法的優(yōu)越性，建立楔狀體模型，對其正演數據加入隨機噪聲，得到含噪記錄如圖1a所示。分別利用常規(guī)BM3D方法和改進的BM3D法進行去噪處理，結果如圖1b和圖1c所示。

圖1 楔狀體模型正演數據兩種方法去噪結果對比(a)原始含噪數據； (b)常規(guī)BM3D去噪結果； (c)本文方法去噪結果

選用峰值信噪比和信噪比兩個參數對去噪效果進行定量比較。二者分別定義為

(11)

(12)

對圖1a所示的模型正演數據，變化所加入噪聲的強度，分別使用常規(guī)BM3D和改進BM3D方法進行去噪處理，并計算其PSNR和SNR，結果如圖2

所示。從圖中曲線可以看出，改進BM3D去噪方法相對常規(guī)BM3D方法去噪后SNR和PSNR均更高，尤其是低信噪比情況下去噪效果更顯著。

為了進一步說明本文方法相對于常規(guī)方法的優(yōu)越型，使用圖3a所示的模型數據進行測試。對模型數據加入隨機噪聲，得到信噪比為3dB(即信號與噪聲的能量之比為2)的含噪記錄，如圖3b所示。分別使用常規(guī)BM3D方法、曲波變換濾波方法和本文方法進行去噪處理，結果如圖3c、圖3e和圖3g所示，去噪前后的差值剖面即去除的噪聲如圖3d、圖3f和圖3h所示。由圖可以看出，加入隨機噪聲后模型斷點處的反射基本淹沒在噪聲中，斷層邊緣位置模糊不清。 常規(guī)BM3D方法能夠在一定程度 上去除部分噪聲，但去噪后同相軸連續(xù)性降低(圖3c)；曲波變換濾波后的剖面(圖3e)同相軸彎曲部分連續(xù)性降低，斷層邊緣有效信號在殘差剖面中明顯(圖3f中紅框處)，對斷層細節(jié)信息保持較差。本文方法不僅能明顯去除隨機噪聲，且去噪后反射同相軸清晰，斷層邊界特征也得到了良好的保持(圖3g)。從圖3h所示的去除噪聲剖面上可以看出，本方法未將有效信號作為噪聲去除，很好地保留了有效信息，去噪效果相對最佳。圖4為三種方法去噪后的SNR和PSNR對比，可以看出本文方法的去噪效果最好。

圖2 兩種方法去噪效果的定量對比(a)峰值信噪比； (b)信噪比

圖3 不同方法去噪效果對比(a)原始模型數據； (b)加入隨機噪聲的模型數據； (c)常規(guī)BM3D方法去噪結果； (d)常規(guī)BM3D方法去除的噪聲； (e)曲波變換方法去噪結果； (f)曲波變換方法去除的噪聲； (g)本文方法去噪結果； (h)本文方法去除的噪聲

圖4 不同方法去噪結果的定量對比(a)SNR； (b)PSNR A：常規(guī)BM3D方法； B：曲波變換法； C：本文方法

3 實際資料試算

圖5a為實際疊后地震剖面，其中含有隨機噪聲，影響了資料品質。分別使用常規(guī)BM3D法、曲波變換法和本文方法進行去噪，結果如圖5b、圖5c和圖5d所示。

圖5 不同方法實際資料去噪結果對比

(a)原始剖面； (b)常規(guī)BM3D方法； (c)曲波變換法； (d)本文方法

為了說明本方法對斷層邊界信息的保持能力，從圖5中截取存在多個斷層的紅框區(qū)域進行放大顯示，如圖6a～圖6d所示?？梢姡撼Ｒ?guī)BM3D去噪結果邊界信息不夠清楚，弱反射層同相軸連續(xù)性較差(圖6b)；曲波變換方法保持邊界能力較差，斷層邊緣模糊(圖6c)；本文去噪方法在去噪后斷點清晰，斷層解釋更準確，較弱的反射同相軸連續(xù)性強(圖6d)。圖6e～圖6g為三種方法去除的噪聲剖面，可見常規(guī)BM3D法和本文方法去除的噪聲剖面均勻，表明未去除有效信號；而曲波變換去除的噪聲剖面中部分斷點處存在有效信號(紅框所示)。綜合去噪后剖面和殘差剖面可以看出，本文的方法去噪效果最好。

圖6 斷層區(qū)局部去噪效果對比放大顯示(a)實際資料去噪前； (b)常規(guī)BM3D法去噪結果； (c)曲波變換法去噪結果； (d)本文方法去 噪結果； (e)常規(guī)BM3D法去除的噪聲； (f)曲波變換法去除的噪聲； (g)本文方法去除的噪聲

4 結論

BM3D方法在地震資料去噪中具有良好的效果，但因無法準確估計噪聲方差并選取合適的處理參數，在實際應用中受到一定的限制。本文聯合曲波噪聲估計和BM3D去噪方法，通過曲波分析，實現了濾波閾值和塊匹配參數的自適應選取和去噪處理，提高了去噪精度和計算效率，具有較強的實際意義。理論模型及實際資料處理結果表明：

(1)基于曲波噪聲估計的BM3D去噪方法能準確選取濾波參數，去噪后剖面的信噪比和峰值信噪比更高；

(2)本文方法避免了參數選取的多次測試，提高了計算效率；

(3)由于實現了濾波閾值的自適應選取，本文方法較好地保護弱反射信號和邊界細節(jié)信息，為后續(xù)地震資料的精細解釋提供保障。