孫思亮 ，劉懷山,2

(1.中國(guó)海洋大學(xué) 海底科學(xué)與探測(cè)技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室，山東 青島 266100；2.海洋國(guó)家實(shí)驗(yàn)室 海洋礦產(chǎn)資源評(píng)價(jià)與探測(cè)技術(shù)功能實(shí)驗(yàn)室，山東 青島 266071)

1 引 言

實(shí)際采集到的地震數(shù)據(jù)常?；祀s著隨機(jī)噪聲，地震資料的質(zhì)量較低[1]，后續(xù)的地震數(shù)據(jù)處理分析難度大，因此地震數(shù)據(jù)去噪是必不可少的環(huán)節(jié)。

目前常用的壓制隨機(jī)噪聲的方法有f-x域預(yù)測(cè)濾波[2]、中值濾波[3]、小波分析[4,5]、奇異值分解[6,7]等。近年來(lái)，隨著小波理論的深入研究，Ridgelet[8]、Curvelet[9]、Seislet[10]、Shearlet[11]等變換的多尺度幾何分析方法相繼問(wèn)世，地震數(shù)據(jù)處理更加的多元化。

Curvelet變換因其多尺度、多方向性的特點(diǎn)已在圖像處理方面取得了一些成果[12-14]。在地震勘探領(lǐng)域，楊凱等[15]給出了曲波域地震數(shù)據(jù)隨機(jī)噪聲壓制的軟硬閾值折中方法。孫成禹等[16]利用曲波變換估計(jì)地震數(shù)據(jù)噪聲結(jié)合改進(jìn)的自適應(yīng)三維塊匹配算法對(duì)地震數(shù)據(jù)進(jìn)行去噪，取得了較好的效果。路鵬飛等[17]將Curvelet變換和傅里葉變換聯(lián)合，進(jìn)而壓制面波。楊會(huì)等[18]將Curvelet變換與二維經(jīng)驗(yàn)?zāi)B(tài)分解結(jié)合，對(duì)隨機(jī)噪聲進(jìn)行壓制，取得了一定的效果。

非局部均值(Non-Local Means)算法提出以后，常常被用于圖像降噪處理[19，20]。后來(lái)，Peng等[21]和Jacques等[22]改進(jìn)了該算法，提出了快速非局部均值濾波(Fast Non-Local Means)，大大減小了NLM算法的時(shí)間復(fù)雜度。FNLM算法原理簡(jiǎn)單、容易實(shí)現(xiàn)并且濾波效果優(yōu)良，在地震數(shù)據(jù)處理方面具有廣大的應(yīng)用前景。

實(shí)際地震數(shù)據(jù)通常較為復(fù)雜，常規(guī)曲波變換硬閾值方法的去噪精度往往不夠，并且處理結(jié)果常常存在偽影和過(guò)度平滑等現(xiàn)象，不能達(dá)到理想的去噪效果。針對(duì)這些問(wèn)題，本文引入循環(huán)平移技術(shù)[23]和塊狀復(fù)數(shù)閾值估計(jì)技術(shù)[24]與曲波變換進(jìn)行結(jié)合，對(duì)傳統(tǒng)的曲波閾值去噪進(jìn)行改進(jìn)，進(jìn)而提出一種改進(jìn)的曲波變換閾值去噪和快速非局部均值濾波相結(jié)合的去噪方法。為地震數(shù)據(jù)隨機(jī)噪聲壓制提供一種新的思路。

2 方法原理

2.1 第二代曲波變換

第二代曲波變換易于實(shí)現(xiàn)，運(yùn)算效率高，適用于地震數(shù)據(jù)的處理分析。

對(duì)于信號(hào)f[t1,t2],0≤t1≤t2≤n，它的離散Curvelet變換為[9]：

(1)

將信號(hào)頻譜進(jìn)行分割，得到的局部窗函數(shù)如下[9]：

(2)

定義φ為一維低通窗口內(nèi)積，其中φj(ω1,ω2)=φ(2-jω1)φ(2-jω2)。引入一組等間隔的斜率tanθl=l×2?j/2?,l=-2?j/2?,...,2?j/2?-1，則有:

(3)

(4)

2.2 快速非局部均值算法

非局部均值(NLM)算法去噪的基本思想是：通過(guò)搜索窗口中像素值的加權(quán)平均來(lái)估計(jì)目標(biāo)像素值。其具體去噪過(guò)程為：

假設(shè)輸入的噪聲圖像v為經(jīng)典的噪聲模型：

v(p)=u(p)+n(p)

(5)

(6)

(7)

式(6)和式(7)中，p為目標(biāo)像素；Ω為搜索窗口；q為位于搜索窗口的像素；Z(p)為權(quán)重歸一化參數(shù)；w(p,q)為像素點(diǎn)p，q之間的權(quán)重，一般用歐式距離衡量其相似度[22]。

(8)

式(8)中，σ為噪聲的標(biāo)準(zhǔn)差;h表示和σ相關(guān)的濾波參數(shù);G為兩個(gè)像素點(diǎn)的鄰域的歐式距離;max(·)為取最大值的函數(shù)。

由于非局部均值算法時(shí)間復(fù)雜度過(guò)高，基于積分圖的快速非局部均值算法(FNLM)進(jìn)而出現(xiàn)，其基本思路如下：

與NLM一樣，設(shè)置搜索窗大小為2D×2D,領(lǐng)域窗大小為2d×2d。像素點(diǎn)z和z+r之間的灰度值距離為

(9)

式(9)中，z=(z1,z2)∈Ω，r=(r1,r2)∈[-D,D]×[-D,D]；積分圖定義為：

(10)

圖像的積分圖遞推公式如下[22]：

?(x,y)∈Ω,x≤1,y≤1;

Sum(x,y)=s(x,y)+Sum(x-1,y)

+Sum(x,y-1)-Sum(x-1,y-1)

(11)

則以像素點(diǎn)(x,y)，(x+r1,x+r2)為中心的鄰域窗區(qū)域的灰度值歐氏距離為

(12)

綜上可知，基于積分圖的FNLM算法在計(jì)算兩個(gè)鄰域窗相似度時(shí)，只需要求差值，大幅度降低了NLM算法的時(shí)間復(fù)雜度。

2.3 Curvelet變換和FNLM聯(lián)合去噪方法

當(dāng)?shù)卣饠?shù)據(jù)的信噪比較低，弱有效信息被隨機(jī)噪聲淹沒(méi)時(shí)，直接采用Curvelet變換硬閾值降噪會(huì)在濾除大部分噪聲的同時(shí)，也會(huì)濾除弱有效信號(hào)，并且處理結(jié)果會(huì)出現(xiàn)偽影和過(guò)度平滑等現(xiàn)象。而直接采用快速非均值濾波又會(huì)造成隨機(jī)噪聲壓制不完全的問(wèn)題。因此本文提出一種改進(jìn)的Curvelet變換閾值方法與FNLM方法相結(jié)合去噪方法，以解決上述兩種方法在地震數(shù)據(jù)去噪過(guò)程中的出現(xiàn)問(wèn)題。

本文去噪方法的主要思路如下：①將循環(huán)平移和塊狀復(fù)數(shù)域閾值引入Curvelet變換去噪過(guò)程。循環(huán)平移方法能夠消除傳統(tǒng)曲波變換閾值去噪方法存在的偽吉布斯現(xiàn)象的“振鈴”效應(yīng)[25，26]。塊狀復(fù)數(shù)域閾值能夠較為準(zhǔn)確地獲取地震數(shù)據(jù)的過(guò)濾范圍，盡可能地濾除隨機(jī)噪聲，保留有效信號(hào)，提升地震數(shù)據(jù)的質(zhì)量[27,28]。②利用改進(jìn)的曲波變換閾值方法對(duì)含噪地震數(shù)據(jù)進(jìn)行去噪處理，分離有效信號(hào)和噪聲。③利用FNLM算法對(duì)分離出的噪聲進(jìn)行處理，提取殘留在噪聲中的有效信號(hào)。④將上面兩步的處理結(jié)果整合得到去噪后的地震數(shù)據(jù)。

3 理論模型實(shí)驗(yàn)及分析

為驗(yàn)證改進(jìn)曲波變換閾值去噪結(jié)合FNLM算法在去除地震信號(hào)隨機(jī)噪聲中的優(yōu)越性，首先在層狀模型地震記錄中加入隨機(jī)噪聲，再分別采用該方法和小波閾值去噪、FNLM算法、曲波變換閾值去噪對(duì)層狀模型含噪地震記錄進(jìn)行去噪處理。計(jì)算去噪前后地震數(shù)據(jù)的信噪比(Signal Noise Ratio,SNR)和均方根誤差(Root Mean Square Error,RMSE)，評(píng)價(jià)上述四種方法的去噪效果。SNR值越大，RMSE值越小，表示去噪效果越好。

信噪比參數(shù)：

(13)

均方根誤差參數(shù)：

(14)

式(13)和(14)中，s0表示無(wú)噪聲的地震數(shù)據(jù);s表示去除噪聲后的地震數(shù)據(jù);N表示地震數(shù)據(jù)的大小。

圖1(a)為一個(gè)共100道，每道512個(gè)采樣點(diǎn)的模擬地震記錄。對(duì)模擬地震記錄中加入的隨機(jī)噪聲，信噪比為-9.59 dB，有效信號(hào)被淹沒(méi)(圖1b)。

圖1 層狀模型合成地震記錄及其加噪后記錄Fig.1 Layered model synthetic seismogram and seismic records with added noise

圖2 四種去噪方法結(jié)果對(duì)比Fig.2 Comparison of four de noising methods

圖3 合成地震記錄去噪前后剖面頻譜對(duì)比Fig.3 Spectrum correlation of synthetic seismogram before and after denoising

分別采用FNLM算法、小波閾值去噪、Curvelet變換閾值去噪和本文方法對(duì)層狀模型含噪地震記錄進(jìn)行隨機(jī)噪聲壓制處理。

圖2(a)～圖2(d)為四種方法對(duì)于層狀模型含噪地震記錄的去噪結(jié)果。可以看到，這四種方法都能夠在一定程度上壓制隨機(jī)噪聲。但FNLM算法、小波閾值方法和Curvelet變換閾值方法去噪都不夠徹底，去噪結(jié)果中仍然有殘余的噪聲，并且弱有效信號(hào)的同相軸沒(méi)能夠較好恢復(fù)(圖2(a)、圖2(b)和圖2(c)矩形方框區(qū))。小波閾值去噪結(jié)果同相軸都發(fā)生了畸變，有效信號(hào)損失較為嚴(yán)重(圖2(b))。圖2(d)為本文方法去噪結(jié)果，與圖1(a)相比，隨機(jī)噪聲幾乎被完全壓制，弱有效信號(hào)也得到了恢復(fù)，保護(hù)了地震信號(hào)的完整性。

表1分別給出了層狀模型含噪地震記錄和應(yīng)用四種去噪方法去噪后地震記錄的信噪比和均方根誤差。對(duì)比可知本文方法得到的地震數(shù)據(jù)的信噪比最高，并且均方根誤差最小，說(shuō)明利用本文方法壓制地震數(shù)據(jù)隨機(jī)噪聲最優(yōu)。

為了進(jìn)一步說(shuō)明本文方法壓制隨機(jī)噪聲的優(yōu)勢(shì)，分別從原始記錄、加噪記錄和四種去噪方法去噪后的記錄中抽取第50道記錄進(jìn)行頻譜分析。比較原始記錄頻譜(圖3(a))和加噪記錄的頻譜(圖3(b))可知，加噪以后的地震記錄高頻部分明顯增加，隨機(jī)噪聲主要表現(xiàn)為高頻。比較圖3(c)、圖3(d)與圖3(a)可以發(fā)現(xiàn)(見(jiàn)圖中的紅色方框)，F(xiàn)NLM算法和小波閾值方法壓制了大部分高頻噪聲，但去噪結(jié)果中仍然存留高頻噪聲。比較圖3(e)、圖3(f)與圖3(a)可以發(fā)現(xiàn)，曲波變換閾值方法和本文方法去噪結(jié)果中幾乎不含高頻噪聲，但圖3(e)相較于圖3(a)中的最大振幅有所降低，說(shuō)明曲波變換閾值方法對(duì)有效信號(hào)有損傷，部分有效信號(hào)被濾除。圖3(f)和圖3(a)的形態(tài)基本一致，說(shuō)明本文方法既壓制了隨機(jī)噪聲也保護(hù)了有效信號(hào)。

表1 四種方法去噪結(jié)果評(píng)價(jià)指標(biāo)對(duì)比

4 實(shí)際地震資料應(yīng)用

將本文方法應(yīng)用到實(shí)際的單炮和疊后地震數(shù)據(jù)處理中，驗(yàn)證本文方法壓制隨機(jī)噪聲的優(yōu)勢(shì)。圖4(a)為疊前地震單炮記錄，隨機(jī)噪聲較多，數(shù)據(jù)信噪比較低，弱同相軸幾乎被淹沒(méi)(見(jiàn)圖4(a)箭頭處)，有效信號(hào)難以識(shí)別。

圖4 疊前單炮地震數(shù)據(jù)及其去噪結(jié)果Fig.4 Pre—stack single shot seismic data and its denoising results by four methods

同樣利用FNLM算法、小波閾值方法、Curvelet變換閾值方法和本文方法對(duì)疊前單炮地震記錄進(jìn)行隨機(jī)噪聲壓制處理。圖4(b)～圖4(e)為四種方法的去噪結(jié)果。從圖中可以發(fā)現(xiàn)，F(xiàn)NLM算法和小波閾值方法去噪結(jié)果中隨機(jī)噪聲殘留較多，部分有效信號(hào)被濾除(圖4(b)、圖4(c))；Curvelet變換閾值去噪效果較好，但仍有少量有效信號(hào)被濾除，對(duì)于弱有效信號(hào)的保護(hù)不夠理想(圖4(d))；本文方法去噪結(jié)果弱有效信號(hào)豐富，隨機(jī)噪聲去除地較為徹底(圖4(e))。

圖5為含有大量的隨機(jī)噪聲，信噪比較低的實(shí)際疊后地震記錄，層間弱信號(hào)被淹沒(méi)，剖面同相軸的連續(xù)性較差(見(jiàn)圖5箭頭處)，對(duì)地質(zhì)解釋帶來(lái)了嚴(yán)重影響。

為了消除隨機(jī)噪聲影響，為地震資料解釋奠定基礎(chǔ)，同樣分別利用上述四種方法對(duì)實(shí)際疊后地震數(shù)據(jù)進(jìn)行隨機(jī)噪聲壓制處理。圖6(a)～圖6(d)為4種方法的去噪結(jié)果，由圖可知，4種方法均具有一定的去噪能力。但FNLM算法在去除隨機(jī)噪聲的同時(shí)，部分弱幅度有效信號(hào)遭到破壞，導(dǎo)致部分同相軸不連續(xù)(圖6(a))。小波閾值的去噪結(jié)果中還有隨機(jī)噪聲存在，并且有效信號(hào)的同相軸發(fā)生了畸變，同相軸的連續(xù)性差(圖6(b))；Curvelet閾值去噪結(jié)果隨機(jī)噪聲去除地較為干凈，有效信號(hào)較為豐富(圖6(c))；本文方法去噪結(jié)果地震剖面的質(zhì)量得到明顯提升，絕大部分隨機(jī)噪聲被壓制，弱幅度有效信息豐富、同相軸清晰連續(xù)(圖6(d))。

圖5 某工區(qū)疊后地震剖面Fig.5 Post-stack seismic profile of a work area

圖6 實(shí)際疊后地震數(shù)據(jù)去噪結(jié)果Fig.6 Denoising result of post stack seismic data

為了更好地比較FNLM算法(圖6(a))、小波閾值方法(圖6(b))、Curvelet變換閾值方法(圖6(c))和本文方法(圖6(d))的去噪效果，對(duì)圖中紅框標(biāo)示區(qū)域進(jìn)行放大，如圖7所示。從圖中的箭頭所指位置可以發(fā)現(xiàn)，F(xiàn)NLM算法在壓制隨機(jī)噪聲的同時(shí)，有效信號(hào)也受到了衰減(圖7(b))；小波閾值方法去噪結(jié)果中仍殘留較多隨機(jī)噪聲，同相軸發(fā)生畸變，去噪效果欠佳(圖7c)；Curvelet變換閾值方法和本文方法去噪效果較好，在壓制隨機(jī)噪聲的同時(shí)，弱地震信號(hào)也得到了一定恢復(fù)，同相軸的連續(xù)性較好(圖7(d)、圖7(e))。但通過(guò)對(duì)比圖7(d)、圖7(e)和圖7(a)中的矩形方框內(nèi)的地震記錄，可以發(fā)現(xiàn)，相較于Curvelet閾值去噪方法，本文方法更加能夠保護(hù)弱有效信號(hào)，能夠?qū)⑽⑷跤行盘?hào)從強(qiáng)隨機(jī)噪聲中提取出來(lái)。

圖7 實(shí)際疊后地震數(shù)據(jù)局部放大圖比較Fig.7 Comparison of local enlarged images of actual post-stack seismic data

5 結(jié) 論

本文提出了一種基于改進(jìn)Curvelet閾值去噪與FNLM算法相結(jié)合的去噪方法，充分發(fā)揮了兩者的優(yōu)點(diǎn)，克服了傳統(tǒng)曲波變換閾值去噪方法和FNLM算法的缺點(diǎn)。理論模型和實(shí)際地震資料的應(yīng)用結(jié)果表明：與小波閾值方法、FNLM算法和Curvelet變換閾值方法相比，本文方法明顯優(yōu)于其他方法，對(duì)隨機(jī)噪聲的壓制更加徹底，去除隨機(jī)噪聲的能力更強(qiáng)；去噪后的地震數(shù)據(jù)信噪比得到了極大的提升，同相軸更加清晰、連續(xù)性更好，并且能夠最大限度地保留弱有效信號(hào)。