基于曲波噪聲估計(jì)的三維塊匹配地震資料去噪

孫成禹 刁俊才 李文靜

(①中國(guó)石油大學(xué)(華東)地球科學(xué)與技術(shù)學(xué)院，山東青島 266580； ②青島海洋國(guó)家實(shí)驗(yàn)室海洋礦產(chǎn)資源評(píng)價(jià)與探測(cè)技術(shù)功能實(shí)驗(yàn)室，山東青島 266071； ③東方地球物理公司物探技術(shù)研究中心，河北涿州 072751)

0 引言

地震數(shù)據(jù)中噪聲嚴(yán)重影響資料處理和解釋結(jié)果[1]，壓制噪聲對(duì)后期地震資料的處理、反演和解釋等具重要作用[2]。地震資料常用的去噪方法一般分為空間域方法和變換域方法兩大類。空間域去噪方法包括均值濾波、中值濾波[3]和維納濾波[4]等，變換域去噪方法包括傅里葉變換、小波變換、曲波變換和Randon變換等閾值方法[5-6]。這些去噪算法本質(zhì)都是利用了圖像本身的局部相關(guān)性，沒有充分挖掘、利用圖像的非局部相關(guān)性，因此在壓制噪聲的同時(shí)也會(huì)損害部分有效信號(hào)，不能精確保持?jǐn)鄬印⒘芽p等地質(zhì)體的邊緣特征，甚至?xí)霈F(xiàn)偽影等現(xiàn)象[7]。

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

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

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

1 基本原理

1.1 三維塊匹配去噪原理

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

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

(1)

(2)

(3)

(4)

(5)

(6)

將維納濾波應(yīng)用于噪聲數(shù)據(jù)矩陣Yr，即將噪聲數(shù)據(jù)矩陣三維變換后的系數(shù)與維納濾波收縮系數(shù)相乘，再進(jìn)行反變換后得到新的估計(jì)值為

(7)

最后將不同群組中的重疊塊進(jìn)行加權(quán)平均得到重新估計(jì)值

(8)

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

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

(9)

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

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

(10)

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

2 模型試算

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

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

選用峰值信噪比和信噪比兩個(gè)參數(shù)對(duì)去噪效果進(jìn)行定量比較。二者分別定義為

(11)

(12)

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

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

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

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

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

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

3 實(shí)際資料試算

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

圖5 不同方法實(shí)際資料去噪結(jié)果對(duì)比

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

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

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

4 結(jié)論

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

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

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

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