能量歸一塊匹配三維協(xié)同濾波及地震去噪應(yīng)用

梁晨曦,張固瀾,2,李 磊,羅 帆,羅一梁,李 勇,段 景,杜 皓

(1.西南石油大學(xué)地球科學(xué)與技術(shù)學(xué)院,四川成都610500;2.油氣藏地質(zhì)與開(kāi)發(fā)國(guó)家重點(diǎn)實(shí)驗(yàn)室,四川成都610500;3.中國(guó)石油集團(tuán)東方地球物理勘探有限責(zé)任公司物探技術(shù)研究中心,河北涿州072751)

在地震勘探過(guò)程中,隨機(jī)噪聲的存在影響地震數(shù)據(jù)的后續(xù)處理精度,因此壓制噪聲對(duì)后期地震資料的處理、反演和解釋等具重要作用。地震資料去噪方法一般分為原始數(shù)據(jù)域方法[1-2]和變換域方法[3-6]。原始數(shù)據(jù)域方法主要包括中值濾波[7-8]、維納濾波[9]、各向異性擴(kuò)散濾波等[10];變換域方法主要包括離散余弦變換[11-12]、小波變換[13-14]、曲波變換等[15]。這些算法由于缺少對(duì)非局部信息的利用,因此在去除噪聲的同時(shí)會(huì)損害部分有效信號(hào)。

DABOV等[16]提出的塊匹配三維協(xié)同濾波(Block-matching and 3D collaborative Filtering,BM3D)方法將數(shù)據(jù)域與變換域結(jié)合,充分利用數(shù)據(jù)自相似性和冗余性信息對(duì)隨機(jī)噪聲進(jìn)行壓制,其核心步驟如下:①將原始數(shù)據(jù)分解成多個(gè)相同大小的分塊數(shù)據(jù)[17-18];②對(duì)各分塊數(shù)據(jù)進(jìn)行二維變換,即先在列(行)方向進(jìn)行一維變換,直至所有列(行)都處理結(jié)束,在此基礎(chǔ)上,再在行(列)方向進(jìn)行一維變換,直至所有行(列)都處理結(jié)束,得到二維變換域分塊數(shù)據(jù);③在二維變換域,按目標(biāo)塊和各分塊數(shù)據(jù)的相似程度和給定的閾值進(jìn)行塊匹配[19-20](按照匹配后相似塊的塊順序號(hào)進(jìn)行),得到分組后的三維數(shù)據(jù)體;④對(duì)分組后的三維數(shù)據(jù)體在塊序號(hào)方向進(jìn)行一維變換得到三維變換域數(shù)據(jù),并進(jìn)行三維協(xié)同濾波(硬閾值濾波和經(jīng)驗(yàn)維納濾波)去噪處理;⑤對(duì)分塊數(shù)據(jù)塊聚集重構(gòu)圖像,得到最終的去噪結(jié)果。

BM3D方法在去噪[21-22]的同時(shí)可較好地保留原始數(shù)據(jù)的有效信息[23],因此被廣泛用于圖像處理和地震數(shù)據(jù)處理等[24-25]領(lǐng)域。劉向樂(lè)等[26]提出的變換域方法為小波變換的BM3D方法,在圖像處理中取得了較好的應(yīng)用效果。任婷婷等[27]提出的變換域方法為離散余弦變換(discrete cosine transform,DCT)的BM3D方法,可用于地震數(shù)據(jù)隨機(jī)噪聲壓制。孫成禹等[28]將曲波變換與BM3D相結(jié)合,提出了基于曲波估計(jì)噪聲的BM3D方法,并用于實(shí)際地震數(shù)據(jù)的隨機(jī)噪聲壓制,該方法利用曲波變換估計(jì)地震數(shù)據(jù)噪聲強(qiáng)度,并基于估計(jì)的噪聲強(qiáng)度自適應(yīng)調(diào)整BM3D中的濾波閾值和塊匹配參數(shù)。張歡等[29]提出結(jié)合主成分分析的BM3D方法,該方法通過(guò)對(duì)地震數(shù)據(jù)進(jìn)行降維,得到地震數(shù)據(jù)噪聲估計(jì)值,并基于估計(jì)噪聲值規(guī)避傳統(tǒng)噪聲預(yù)估值的局限性,優(yōu)化了BM3D對(duì)地震數(shù)據(jù)的噪聲估計(jì)。

BM3D方法的塊匹配及分組結(jié)果受塊間能量差異影響較大,去噪精度有待提升;另外,BM3D方法中的參數(shù)選取(如濾波閾值等)需經(jīng)過(guò)多次試驗(yàn)才能確定,參數(shù)選取效率較低。為提高BM3D方法的去噪效果和參數(shù)選取效率,本文提出了基于分塊能量歸一化的BM3D,介紹了該方法的基本原理及具體實(shí)現(xiàn)步驟,并將其與傳統(tǒng)方法(中值濾波、DCT濾波、BM3D)進(jìn)行對(duì)比。最后通過(guò)模型數(shù)據(jù)和實(shí)際數(shù)據(jù)測(cè)試與應(yīng)用,驗(yàn)證本文方法的有效性。

1 塊匹配三維協(xié)同濾波基本原理

BM3D由基礎(chǔ)估計(jì)和最終估計(jì)組成,其核心步驟如圖1所示,主要分為以下3個(gè)方面:①硬閾值塊匹配及分組;②三維協(xié)同濾波,即基礎(chǔ)估計(jì)階段的硬閾值濾波和最終估計(jì)階段的經(jīng)驗(yàn)維納濾波;③塊聚集。

圖1 傳統(tǒng)BM3D算法基本流程

1.1 硬閾值塊匹配及分組

(1)

且

(2)

(3)

其中,{ }表示集合。

假設(shè)Di[m][n],Dj[m][n]和Gj[m][n]分別為分塊數(shù)據(jù)Di,Dj和相似塊匹配結(jié)果Gj的離散表示,則:

(4)

且

(5)

(6)

且

(7)

由(1)式、(2)式、(4)式和(5)式可知,τ1和τ2越小,搜索到的相似塊越少;τ1和τ2越大,搜索到的相似塊越多。由(3)式和(5)式可知,歐氏距離d(Bi,Bj)和d(Di,Dj)受分塊數(shù)據(jù)總能量差異影響較大,并不能完全反映分塊數(shù)據(jù)間的相似度。因此,在分塊數(shù)據(jù)總能量差異較大的情況下,τ1和τ2的選擇會(huì)影響B(tài)M3D的計(jì)算效率,需多次試驗(yàn);且搜索到的相似塊不能真正反映搜索塊與目標(biāo)塊間的相似度,影響了三維協(xié)同濾波結(jié)果,最終影響B(tài)M3D去噪效果。

1.2 三維協(xié)同濾波

(8)

且

(9)

其中,|·|表示取模。

(10)

且

(11)

其中,δ∈[0,+∞)。

由(8)式至(11)式可知,λ1越大,噪聲壓制效果越明顯;δ越大,噪聲壓制效果越明顯。在分塊數(shù)據(jù)總能量差異較大情況下,λ1和δ會(huì)影響B(tài)M3D計(jì)算效率和三維協(xié)同濾波結(jié)果,最終影響B(tài)M3D去噪效果,其值的選擇需多次試驗(yàn)。

1.3 塊聚集

H0[x][y]=

(12)

H1[x][y]=

(13)

2 塊能量歸一匹配三維協(xié)同濾波基本原理

BM3D方法易受塊間能量差異影響且參數(shù)選取效率不高,相似塊匹配精度、處理參數(shù)選取效率及濾波效果都有待提升。為消除上述影響,本文提出塊能量歸一BM3D方法,主要包括以下3個(gè)方面:①能量歸一軟/硬閾值塊匹配及分組;②能量歸一軟/硬閾值三維協(xié)同濾波;③能量歸一塊聚集。其核心步驟如圖2 所示(改進(jìn)部分為圖2中紅色模塊),與BM3D方法的不同之處主要體現(xiàn)在①和②兩方面。

圖2 改進(jìn)的BM3D算法基本流程

2.1 能量歸一軟/硬閾值塊匹配及分組

(14)

(15)

式中:max(·)表示取最大值;μ1為實(shí)數(shù),且μ1∈[0,1];Ei表示分塊數(shù)據(jù)Fi的能量。

變換域軟/硬閾值相似塊匹配可表示為:

(16)

且

(17)

式中:軟閾值τ2為常數(shù),且τ2∈[0,1];硬閾值τ3為常數(shù),且τ3∈[0,1]。

(18)

2.2 能量歸一軟/硬閾值三維協(xié)同濾波

(19)

式中:μ2為實(shí)數(shù),且μ2∈[0,1]。

2.3 能量歸一塊聚集

最終估計(jì)塊聚集過(guò)程如下:①濾波數(shù)據(jù)分塊;②原始圖像重構(gòu)。實(shí)現(xiàn)方式與傳統(tǒng)方法一致。

3 模型數(shù)據(jù)試算

通過(guò)模型試算來(lái)對(duì)比分析傳統(tǒng)BM3D方法與本文方法的相似塊匹配精度和濾波效果。本文方法中,變換域方法均為二維DCT,圖中色標(biāo)數(shù)值均代表振幅值。

3.1 模型試算一

圖3a為12×12個(gè)像素點(diǎn)的無(wú)噪理論模型數(shù)據(jù)的數(shù)字顯示,其中,圖3a右上角、左下角和右下角藍(lán)框區(qū)域像素值分別為紅框區(qū)域像素值的2倍、5倍和10倍;圖3b為圖3a的圖像顯示;圖3c為對(duì)圖3a加入隨機(jī)噪聲后的圖像顯示。

圖3c中黃色方塊為目標(biāo)塊,利用傳統(tǒng)的塊匹配方法和本文改進(jìn)的塊匹配方法,進(jìn)行相似塊匹配處理。圖4為不同τ1的傳統(tǒng)塊能量歸一硬閾值塊匹配結(jié)果。圖5為不同τ2的塊能量歸一軟閾值塊匹配結(jié)果。圖6為不同τ3的塊能量歸一硬閾值塊匹配結(jié)果。對(duì)比圖3至圖6可知,傳統(tǒng)塊匹配結(jié)果受塊間能量差異影響,精度不高;且隨著硬閾值的不斷增大,其相似塊匹配精度提升效果不明顯;基于塊能量歸一的軟/硬閾值塊匹配均消除了塊間能量差異影響,其相似塊匹配精度較高。

圖3 無(wú)噪理論模型的數(shù)字顯示(a)、無(wú)噪理論模型的圖像顯示(b)和含噪理論模型的圖像顯示(c)

圖4 不同τ1的傳統(tǒng)塊能量歸一硬閾值塊匹配結(jié)果a τ1=16;b τ1=200;c τ1=1000

圖5 不同τ2的塊能量歸一軟閾值塊匹配結(jié)果a τ2=0.08;b τ2=0.20;c τ2=0.50

圖6 不同τ3的塊能量歸一硬閾值塊匹配結(jié)果a τ3=0.02;b τ3=0.05;c τ3=0.20

3.2 模型試算二

圖7a為含有復(fù)雜斷層的無(wú)噪合成地震數(shù)據(jù),圖7b 為對(duì)圖7a加入隨機(jī)噪聲后的合成地震數(shù)據(jù)。圖8a、圖8b、圖8c和圖8d分別為對(duì)圖7b使用二維中值濾波、二維DCT濾波、傳統(tǒng)BM3D方法和本文方法后的濾波結(jié)果,濾除的噪聲分別如圖9a、圖9b、圖9c和圖9d所示。圖10為4種方法去噪后的峰值信噪比(peak signal to noise ratio,PSNR)[30]的對(duì)比結(jié)果。

由圖7至圖10可知:二維中值濾波和二維DCT濾波后的剖面斷點(diǎn)不清晰,邊緣保持能力不強(qiáng)且對(duì)有效信息傷害較大(如圖8a、圖8b、圖9a和圖9b中黑色箭頭所示);傳統(tǒng)的BM3D方法濾后的剖面斷點(diǎn)較清晰,邊緣保持能力較強(qiáng),且對(duì)有效信號(hào)傷害較小(如圖8c 中黑色箭頭所示);本文方法較傳統(tǒng)的BM3D方法的濾波效果有所提升,且濾除的噪聲剖面中幾乎不含有效信號(hào)(如圖9c、圖9d中黑色箭頭所示)。對(duì)比4種方法的PSNR值可知,本文方法峰值信噪比最大,較傳統(tǒng)方法去噪效果有明顯提升。

圖7 原始無(wú)噪合成地震數(shù)據(jù)(a)和含噪合成地震數(shù)據(jù)(b)

圖8 采用不同濾波方法的濾波結(jié)果a 二維中值濾波;b 二維DCT濾波;c 傳統(tǒng)BM3D方法;d 本文方法

圖9 采用不同濾波方法濾除的噪聲a 二維中值濾波;b 二維DCT濾波;c 傳統(tǒng)的BM3D方法;d 本文方法

圖10 4種方法去噪后的峰值信噪比對(duì)比

4 實(shí)際地震數(shù)據(jù)應(yīng)用

采用傳統(tǒng)BM3D濾波方法和本文改進(jìn)的BM3D濾波方法對(duì)實(shí)際三維地震數(shù)據(jù)進(jìn)行去噪處理,具體步驟如下:①對(duì)三維地震數(shù)據(jù)體按Inline(或Crossline)線讀取數(shù)據(jù),并對(duì)讀取的Inline線剖面數(shù)據(jù)進(jìn)行去噪處理;②重復(fù)步驟①,直至所有的Inline線(或Crossline線)地震數(shù)據(jù)去噪處理都結(jié)束,得到去噪后的三維地震數(shù)據(jù)體。

圖11a、圖11b分別為某三維地震數(shù)據(jù)中Inline100和Crossline100原始數(shù)據(jù)。圖12a和圖12b分別為對(duì)圖11a采用傳統(tǒng)BM3D方法和能量歸一塊匹配軟閾值BM3D方法進(jìn)行濾波處理后得到的剖面,濾除的噪聲分別如圖12c和圖12d所示。圖13a和圖13b 分別為對(duì)圖11b采用傳統(tǒng)BM3D方法和能量歸一塊匹配軟閾值BM3D方法進(jìn)行濾波處理后得到的剖面,濾除的噪聲分別如圖13c和圖13d所示。圖14a 是某三維地震數(shù)據(jù)中1s等時(shí)切片,圖14b和圖14c 分別為對(duì)圖14a采用傳統(tǒng)BM3D方法和能量歸一塊匹配軟閾值BM3D方法進(jìn)行濾波處理后得到結(jié)果,濾除的噪聲分別如圖14d和圖14e所示。

圖11 Inline100(a)和Crossline100(b)原始數(shù)據(jù)

由圖12、圖13和圖14可見(jiàn):兩種方法濾波處理后的剖面上同相軸連續(xù)性和信噪比均得到有效提升;改進(jìn)的BM3D方法的濾波結(jié)果對(duì)構(gòu)造刻畫(huà)精度較傳統(tǒng)BM3D方法的濾波結(jié)果有所提升;另外,傳統(tǒng)BM3D方法濾除的噪聲中含有較多有效信號(hào),而改進(jìn)的BM3D方法濾除的噪聲中幾乎不含有效信號(hào),去噪后反射同相軸更清晰且參數(shù)選取效率較高。

圖12 對(duì)Inline100原始數(shù)據(jù)采用傳統(tǒng)BM3D方法和本文方法進(jìn)行濾波處理后得到的剖面a 傳統(tǒng)BM3D方法;b 本文方法;c 傳統(tǒng)BM3D方法濾除的噪聲;d 本文方法濾除的噪聲

圖13 對(duì)Crossline100原始數(shù)據(jù)采用傳統(tǒng)BM3D方法和本文方法進(jìn)行濾波處理后得到的剖面a 傳統(tǒng)BM3D方法;b 本文方法;c 傳統(tǒng)BM3D方法濾除的噪聲;d 本文方法濾除的噪聲

圖14 1s處等時(shí)切片a 原始數(shù)據(jù);b 傳統(tǒng)BM3D方法濾波處理后的結(jié)果;c 本文方法濾波處理后的結(jié)果;d 傳統(tǒng)BM3D方法濾除的噪聲;e 本文方法濾除的噪聲

5 結(jié)論

1)本文提出的基于分塊數(shù)據(jù)總能量歸一的塊匹配及分組算法,可有效提升相似塊匹配精度,從而提升BM3D的應(yīng)用效果;

2)本文提出的基于分塊數(shù)據(jù)總能量歸一的軟/硬閾值濾波方法,通過(guò)軟閾值濾波減少參數(shù)試驗(yàn)次數(shù),從而提高了參數(shù)選取效率;

3)本文方法提升了去噪效果和參數(shù)選取效率,可被廣泛用于地震資料的去噪處理。