王維紅，郭雪豹，石 穎

(東北石油大學(xué)地球科學(xué)學(xué)院，黑龍江大慶163318)

非局部平均濾波噪聲壓制方法及其在VSP資料逆時(shí)偏移中的應(yīng)用

王維紅，郭雪豹，石 穎

(東北石油大學(xué)地球科學(xué)學(xué)院，黑龍江大慶163318)

基于雙程波的逆時(shí)偏移會(huì)產(chǎn)生低頻成像噪聲，在成像后運(yùn)用拉普拉斯濾波法可以取得較好的壓噪效果，但是，該方法嚴(yán)重依賴于角度參數(shù)，使得濾波后的成像剖面上常常存在同相軸不光滑和噪聲壓制不完全的情況?；诏B后地震資料同相軸結(jié)構(gòu)的特點(diǎn)，為進(jìn)一步提高逆時(shí)偏移疊加數(shù)據(jù)的信噪比，引入了非局部平均濾波法，針對拉普拉斯濾波后的成像數(shù)據(jù)進(jìn)行進(jìn)一步的壓噪處理。非局部平均濾波法將輸入的地震數(shù)據(jù)分解為不含噪聲的地震數(shù)據(jù)和噪聲數(shù)據(jù)兩部分，利用不同成像點(diǎn)與其它成像點(diǎn)間的相似系數(shù)，實(shí)現(xiàn)濾波處理。二維復(fù)雜模型VSP正演模擬數(shù)據(jù)逆時(shí)偏移結(jié)果的試算表明，應(yīng)用非局部平均濾波后的剖面信噪比得到進(jìn)一步提高，同相軸的連續(xù)性明顯增強(qiáng)。實(shí)際VSP資料逆時(shí)偏移低頻噪聲壓制試處理也表明非局部平均濾波方法具有計(jì)算精度高、算法穩(wěn)定性好和易于實(shí)現(xiàn)的特點(diǎn)。

非局部平均;噪聲壓制;逆時(shí)偏移;垂直地震剖面;拉普拉斯濾波

隨著計(jì)算機(jī)硬件技術(shù)的發(fā)展，近年來基于雙程波方程的逆時(shí)偏移方法得到了快速發(fā)展和應(yīng)用[1-2]。除地面地震資料的逆時(shí)偏移成像外，井中地震資料，特別是VSP資料的逆時(shí)偏移技術(shù)也得到一定程度的發(fā)展，為井中地震高精度成像奠定了基礎(chǔ)[3-5]。在逆時(shí)偏移方法研究中，一般來說選擇互相關(guān)成像條件進(jìn)行成像計(jì)算，互相關(guān)成像條件可為成像提供準(zhǔn)確的動(dòng)力學(xué)特性[6]，而且實(shí)現(xiàn)方法簡單，但是該成像條件會(huì)使逆時(shí)偏移結(jié)果產(chǎn)生大量強(qiáng)振幅的低頻噪聲[7]。對于低頻噪聲的壓制，主要有波場延拓噪聲壓制方法、應(yīng)用有效成像條件的噪聲壓制方法和成像后濾波的噪聲壓制方法等。

Baysal等[8]提出了無反射波動(dòng)方程來壓制低頻噪聲，但由于該方法與疊前逆時(shí)偏移噪聲的產(chǎn)生機(jī)理有本質(zhì)上的區(qū)別，所以不適用于疊前處理;Leowenthal等[9]提出通過平滑速度模型來減少反射波。對入射角較大的情況，上述兩種方法依舊會(huì)產(chǎn)生反射，因此噪聲壓制效果不很明顯。Fletcher等[10]提出了雙程無反射方程來壓制低頻噪聲，該方法在方程中引入密度項(xiàng)，使得波阻抗為常數(shù)，雖然解決了逆向散射的問題，但只在垂直入射時(shí)應(yīng)用效果好，隨著入射角度的增大，反射波的能量逐漸增強(qiáng)，該方法的實(shí)用性就變得很差。應(yīng)用有效成像條件壓制噪聲的實(shí)質(zhì)就是對互相關(guān)成像條件進(jìn)行改進(jìn)，以達(dá)到噪聲壓制的效果。Kaelin等[11]通過實(shí)驗(yàn)得出結(jié)論，可以借助檢波照明來改善成像條件，分別試驗(yàn)了波場分解成像條件、歸一化成像條件、坡印廷矢量成像條件和立體成像條件，這些方法存在的共同缺陷是對噪聲不能有效完全壓制，使得地震成像數(shù)據(jù)體的解釋和應(yīng)用受到很大影響。成像后濾波法中最簡單的是高通濾波，但這種方法對噪聲壓制不徹底，并且會(huì)破壞有效波場的低頻信息[12]。Guitton等[13]提出在波數(shù)域?qū)υ肼曔M(jìn)行壓制，分別對微分濾波、拉普拉斯濾波、最小平方濾波進(jìn)行了試驗(yàn)。Zhang等[14]指出在疊加偏移的角道集中，通過對大反射角的噪聲抑制來消除噪聲是逆時(shí)偏移低頻噪聲壓制的有效方法，并提出在成像過程中引入拉普拉斯算子進(jìn)行噪聲壓制，運(yùn)用拉普拉斯濾波，相當(dāng)于給角道集加權(quán)。拉普拉斯濾波法在逆時(shí)偏移的研究中的應(yīng)用較為廣泛且有效。

非局部平均(Non-local Means，NLM)濾波最早由Buades等[15]提出的，該算法利用成像點(diǎn)窗口與臨近窗口間的相似性來加強(qiáng)構(gòu)造信息，從而實(shí)現(xiàn)隨機(jī)噪聲的壓制。非局部平均方法最初應(yīng)用于圖像處理，由于非局部平均算法對每一個(gè)成像點(diǎn)降噪時(shí)需要計(jì)算全部成像點(diǎn)的相似系數(shù)，因此具有較大的計(jì)算量[16]。許多學(xué)者為減少計(jì)算時(shí)間做了不同的研究，Sheng等[17]將算法在GPU上運(yùn)行，取得了較好的效果。Mahmoudi等[18]、Brox等[19]為減少每個(gè)成像點(diǎn)的計(jì)算時(shí)間，將求取相似系數(shù)限制在成像點(diǎn)為中心的一個(gè)范圍內(nèi)。這些改進(jìn)算法較原本的算法節(jié)省了大量的計(jì)算時(shí)間。目前非局部平均算法已成功應(yīng)用于醫(yī)學(xué)數(shù)據(jù)、雷達(dá)數(shù)據(jù)、音頻數(shù)據(jù)、顯微成像等領(lǐng)域[20-23]。Bonar等[24]將該方法引入地震資料處理中，進(jìn)行了地震資料隨機(jī)噪聲壓制處理。對應(yīng)用拉普拉斯濾波壓制低頻噪聲后的逆時(shí)偏移地震數(shù)據(jù)體，再應(yīng)用非局部平均濾波方法進(jìn)行噪聲的壓制處理，可進(jìn)一步提高地震數(shù)據(jù)的信噪比，為地震資料的精細(xì)解釋和有效應(yīng)用提供高質(zhì)量的基礎(chǔ)數(shù)據(jù)。

1 非局部平均算法原理

非局部平均濾波依據(jù)圖像結(jié)構(gòu)或地震數(shù)據(jù)同相軸結(jié)構(gòu)的特點(diǎn)，在不同像素點(diǎn)通過對圖像結(jié)構(gòu)增強(qiáng)，實(shí)現(xiàn)噪聲壓制的目的[24]。其基本原理如下。

對一含有噪聲的成像結(jié)果I，可應(yīng)用下式表示：

(1)

式中：U為不含噪聲的成像數(shù)據(jù);N表示噪聲。

(2)

式中：w(m,k)為m和k之間成像值的相似性，其滿足0≤w(m,k)≤1且∑kw(m,k)=1。其中，每個(gè)m都存在與k相關(guān)的相似系數(shù)。為定量表示相似性，定義Nm為以m為中心的窗口，為計(jì)算Nm和Nk之間的相似性，采用高斯加權(quán)歐幾里得距離表示：

(3)

(4)

式中：x0，y0為高斯函數(shù)中心，其中坐標(biāo)x，y與上式中的l項(xiàng)對應(yīng)，u為偏離參數(shù)。

Nm和Nk之間的相似系數(shù)可以表示為：

(5)

圖1說明了成像點(diǎn)間相似系數(shù)計(jì)算原理，圖中虛線框內(nèi)列出了4個(gè)圖像框。計(jì)算以框1中圓點(diǎn)為中心的窗口與以周圍點(diǎn)為中心的窗口之間的相似系數(shù)，如計(jì)算框1中心點(diǎn)與框2中心點(diǎn)之間的相似系數(shù)。首先確定窗口大小，即以點(diǎn)為中心的框所包括的數(shù)據(jù)大小;然后對框1與框2之間的數(shù)據(jù)計(jì)算高斯歐幾里得距離;最后算出兩框中心點(diǎn)之間的相似系數(shù)。Budas等[15]給出的非局部平均濾波是對全局計(jì)算相似系數(shù)，也就是說對于圖1，需要計(jì)算框1中心點(diǎn)與成像體內(nèi)其余所有點(diǎn)之間的相似系數(shù)，這樣會(huì)產(chǎn)生較大的計(jì)算量。所以，僅用以框1圓點(diǎn)為中心的虛線邊框內(nèi)的所有點(diǎn)來計(jì)算相似系數(shù)，可大大減小計(jì)算量，提高計(jì)算效率。

圖1中僅就框1和框2內(nèi)點(diǎn)的順序計(jì)算相似系數(shù)可能較為片面，實(shí)際上為使濾波效果更好，應(yīng)充分考慮輸入數(shù)據(jù)的不同對稱性。圖2給出了以成像點(diǎn)為中心的窗口的選取，在計(jì)算成像點(diǎn)之間相似系數(shù)的時(shí)候，如計(jì)算框1中心點(diǎn)和框2中心點(diǎn)相似系數(shù)時(shí)，對框2做如圖2所示的翻轉(zhuǎn)，并分別與框1求取相似系數(shù)，取其中最大值作為該點(diǎn)的相似系數(shù)，可使濾波效果更好。

圖1 成像點(diǎn)間相似系數(shù)計(jì)算原理

圖2 成像點(diǎn)窗口對稱性示意

2 VSP資料逆時(shí)偏移中的低頻噪聲

眾所周知，單程波偏移結(jié)果沒有強(qiáng)能量噪聲或者說不存在成像噪聲，是由于單程波偏移方法限制了波場的傳播方向，炮點(diǎn)正向延拓的波場遇到反射面時(shí)只考慮向下傳播，因此成像只發(fā)生在炮點(diǎn)與檢波點(diǎn)波場傳播方向相反的位置。

與單程波偏移方法不同，基于雙程波的逆時(shí)偏移會(huì)產(chǎn)生低頻噪聲。根據(jù)Claerbout[6]提出的時(shí)間一致性成像原理及互相關(guān)成像條件，成像點(diǎn)位于震源波場外推與接收波場外推時(shí)間相一致之處。對于雙程延拓的逆時(shí)偏移來說，由于不限制傳播方向，在射線路徑上的每個(gè)點(diǎn)都滿足ts=tr，也就是說在整個(gè)射線傳播的路徑上，成像條件都是滿足的，因此會(huì)在整條射線傳播路徑上成像。如果直接利用公式作為成像條件，會(huì)產(chǎn)生非常強(qiáng)的低頻噪聲，特別是在淺層強(qiáng)反射界面上方，產(chǎn)生的低頻噪聲足以將淺層構(gòu)造完全淹沒。圖3給出了疊前逆時(shí)偏移低頻噪聲的產(chǎn)生機(jī)制[7]。

逆時(shí)偏移低頻噪聲壓制應(yīng)用最為廣泛的方法是拉普拉斯濾波法。Zhang等[14]采用對大反射角的噪聲抑制法來壓制低頻噪聲，指出拉普拉斯濾波相當(dāng)于給角道集加權(quán)。拉普拉斯算子噪聲壓制方法簡單易操作，一般將其寫成二階差分的形式用于逆時(shí)偏移成像后濾波，對噪聲壓制效果較為明顯。

二維拉普拉斯算子可表示為：

(6)

在波數(shù)域可以表示為

(7)

其中，kx和kz分別為x和z方向的波數(shù)。kI為成像域中的波數(shù)矢量kI=kr-ks,其中ks和kr分別為震源處波場和檢波點(diǎn)處波場的波數(shù)矢量。

圖3 逆時(shí)偏移低頻噪聲產(chǎn)生機(jī)制(據(jù)劉紅偉等[7])

由余弦定理得到：

(8)

由(7)式和(8)式可得:

(9)

從(9)式可以看出，采用拉普拉斯算子進(jìn)行噪聲壓制相當(dāng)于給逆時(shí)偏移成像結(jié)果乘以一個(gè)權(quán)系數(shù)，這個(gè)系數(shù)與入射角θ有關(guān)，θ越大，權(quán)系數(shù)越小;θ越小，權(quán)系數(shù)越大。當(dāng)θ=90°時(shí)，噪聲可完全消除。為滿足保幅的處理要求，在將拉普拉斯算子用于逆時(shí)偏移數(shù)據(jù)體的噪聲壓制時(shí)，(9)式中的速度v可以直接從速度模型提取，頻率ω需要在成像前完成補(bǔ)償。

為試驗(yàn)VSP逆時(shí)偏移低頻噪聲壓制效果，設(shè)計(jì)了如圖4所示的復(fù)雜二維理論模型，縱向、橫向網(wǎng)格點(diǎn)數(shù)分別為800和600，網(wǎng)格間距均為5m。采用地面放炮井中接收的觀測方式。震源從地表左端100m開始，向右每隔100m模擬一炮VSP記錄，共39炮。設(shè)計(jì)兩口VSP井接收，分別放于模型的最左端和最右端，分別記為1號井(左端)和2號井(右端)，井中自地表開始部署600個(gè)檢波器，檢波器的間隔為5m。震源為主頻30Hz的Ricker子波，時(shí)間采樣間隔(計(jì)算步長)為0.5ms，采樣點(diǎn)數(shù)為5000。圖5給出的是左端1號井的正演炮記錄，圖5a為第1炮的炮記錄(震源位于模型地表100m處)，圖5b為第20炮的記錄，也就是震源位于模型地表中間位置(偏移距2000m)的正演單炮記錄?？梢钥闯觯瑘D5a的上行波和下行波同相軸很清晰，易于識(shí)別;圖5b震源遠(yuǎn)離井口，波場特征變得比較復(fù)雜，同相軸很難識(shí)別，同時(shí)也表現(xiàn)出了復(fù)雜構(gòu)造模型的特點(diǎn)。

圖4 雙井變偏VSP的井間二維速度模型

圖5 正演模擬的左端1號井VSP單炮記錄

基于二維聲波方程，采用優(yōu)化十二階差分，對以上模擬VSP數(shù)據(jù)進(jìn)行了逆時(shí)偏移計(jì)算。圖6a為變偏VSP逆時(shí)偏移互相關(guān)成像的結(jié)果，具有明顯的低頻噪聲；圖6b為經(jīng)拉普拉斯濾波后的成像數(shù)據(jù)，可以看出低頻噪聲得到較好的壓制，表明拉普拉斯濾波方法可以取得不錯(cuò)的低頻噪聲壓制效果。但從圖6b中也可以看出，在高速地質(zhì)體的兩側(cè)，同相軸的光滑性和連續(xù)性都不令人滿意，對后續(xù)地質(zhì)構(gòu)造解釋造成較大困難。

圖6 模擬VSP數(shù)據(jù)逆時(shí)偏移結(jié)果(a)和經(jīng)拉普拉斯濾波法壓制低頻噪聲后的結(jié)果(b)

3 非局部平均濾波法噪聲壓制

基于非局部平均濾波基本理論，編程實(shí)現(xiàn)了非局部平均濾波噪聲壓制方法。對圖6給出的經(jīng)拉普拉斯濾波后的成像數(shù)據(jù)體進(jìn)行非局部平均濾波處理。在實(shí)際處理時(shí)，所需輸入的速度文件可直接從逆時(shí)偏移的速度模型讀出，理論上講，輸入的模型越精確，濾波后噪聲壓制效果越好。為檢驗(yàn)方法的實(shí)用性，將圖4所示的二維速度模型進(jìn)行平滑，平滑之后的速度模型如圖7所示，用該模型作為非局部平均濾波求取相似系數(shù)的輸入速度文件。將以上模擬VSP數(shù)據(jù)的逆時(shí)偏移結(jié)果(圖6a)進(jìn)行拉普拉斯濾波后(圖6b)，再進(jìn)行非局部平均濾波獲得的結(jié)果如圖8所示。對比圖6b，可以看出圖8中噪聲明顯減少，同相軸的光滑性和連續(xù)性均得到增強(qiáng)。為更清晰地對比兩種方法的噪聲壓制效果，進(jìn)行局部放大顯示。圖9a是圖6b的局部放大圖，圖9b是圖8對應(yīng)位置的局部放大，對比可知，經(jīng)非局部平均濾波后，成像數(shù)據(jù)體的噪聲得到明顯壓制，地震資料的信噪比得到有效提高，為后續(xù)地震資料的精細(xì)解釋和應(yīng)用提供了高質(zhì)量的基礎(chǔ)數(shù)據(jù)。

圖7 圖4所示速度模型平滑后的結(jié)果

圖8 圖6b剖面再經(jīng)非局部平均濾波后的結(jié)果

圖9 圖6b剖面(a)和圖8剖面(b)去噪效果的局部放大顯示

為進(jìn)一步測試非局部平均濾波方法的實(shí)際應(yīng)用效果，對某地區(qū)實(shí)際變偏VSP數(shù)據(jù)進(jìn)行了試處理。該數(shù)據(jù)共411炮，40級檢波器接收。速度模型橫向和縱向網(wǎng)格間距均為10m，時(shí)間步長為0.5ms。圖10給出了實(shí)際VSP數(shù)據(jù)逆時(shí)偏移結(jié)果的低頻噪聲壓制效果對比，圖10a為逆時(shí)偏移結(jié)果經(jīng)拉普拉斯濾波后的成像剖面，圖10b為進(jìn)一步進(jìn)行非局部平均濾波后的成像剖面。對比圖10a與圖10b可以看出，在沒有改變同相軸形態(tài)的情況下，應(yīng)用非局部平均濾波后偏移噪聲進(jìn)一步得到壓制，剖面整體信噪比得到了有效提高。

需要說明的是，非局部平均濾波方法是對疊后的數(shù)據(jù)體進(jìn)行濾波操作，因此其計(jì)算時(shí)間與疊前偏移的炮數(shù)無關(guān)，僅與偏移數(shù)據(jù)體的網(wǎng)格大小和非局部平均濾波所采用的計(jì)算參數(shù)有關(guān)，以圖10b所給出的實(shí)際數(shù)據(jù)為例，非局部平均濾波計(jì)算時(shí)間為105s，可見該濾波方法的計(jì)算效率很高，具有很強(qiáng)的實(shí)用性。

圖10 實(shí)際變偏VSP數(shù)據(jù)逆時(shí)偏移結(jié)果的低頻噪聲壓制效果對比

4 結(jié)論與建議

1) 基于雙程波動(dòng)方程的逆時(shí)偏移是當(dāng)前復(fù)雜構(gòu)造地震成像中精度高、應(yīng)用廣泛的方法，但是該方法會(huì)產(chǎn)生能量很強(qiáng)的低頻噪聲。在當(dāng)前的低頻噪聲壓制方法中，拉普拉斯濾波法應(yīng)用效果較好。

2) 非局部平均濾波噪聲壓制方法應(yīng)用地震數(shù)據(jù)中不同成像點(diǎn)與其它成像點(diǎn)之間的相似性進(jìn)行噪聲壓制，對逆時(shí)偏移后的成像數(shù)據(jù)體幾乎沒有任何假設(shè)，具有較強(qiáng)的適應(yīng)性。另外，該方法在疊后數(shù)據(jù)體上進(jìn)行噪聲壓制處理，計(jì)算效率較高，易于實(shí)現(xiàn)，可根據(jù)成像數(shù)據(jù)體的情況選擇不同的計(jì)算參數(shù)。

3) 在應(yīng)用拉普拉斯濾波對理論模擬的復(fù)雜構(gòu)造變偏VSP數(shù)據(jù)逆時(shí)偏移成像結(jié)果進(jìn)行低頻噪聲壓制后，再應(yīng)用非局部平均濾波噪聲壓制方法可使數(shù)據(jù)的信噪比得到進(jìn)一步提高，同相軸的連續(xù)性得到增強(qiáng)。實(shí)際數(shù)據(jù)試算結(jié)果表明該方法可有效應(yīng)用于VSP地震資料逆時(shí)偏移中。當(dāng)然，該方法也可以應(yīng)用于地面地震資料逆時(shí)偏移成果數(shù)據(jù)的低頻噪聲壓制，具有廣泛的適用性。

[1] Zhang Y,Duan L,Xie Y. A stable and practical implementation of least-squares reverse time migration[J]. Geophysics,2015,80(1):23-31

[2] 劉文卿,王西文,劉洪，等. 鹽下構(gòu)造速度建模與逆時(shí)偏移成像研究及應(yīng)用[J].地球物理學(xué)報(bào),2013,56(2):616-625 Liu W Q,Wang X W,Liu H,et al. Application of velocity modeling and reverse time migration to subsalt structure[J]. Chinese Journal of Geophysics (in Chinese),2013,56(2):612-625

[3] 孫贊東. 三維三分量VSP方法原理及應(yīng)用[M]. 北京:石油工業(yè)出版社,2011:1-4 Sun Z D. The principle and application of 3D three-component VSP method[M]. Beijing:Petroleum industry press,2011:1-4

[4] 朱金明,嚴(yán)俊華. VSP逆時(shí)偏移處理[J]. 石油地球物理勘探,1991,26(5):564-570 Zhu J M,Yan J H. VSP reverse-time migration[J]. Oil Geophysical Prospecting,1991,26(5):564-570

[5] 朱金明,董敏煜,李承楚. VSP的雙程無反射波動(dòng)方程逆時(shí)偏移[J]. 石油地球物理勘探,1989,24(3):256-270 Zhu J M,Dong M Y,Li C C. VSP reverse-time migration using two-way nonreflection wave equation[J]. Oil Geophysical Prospecting,1989,24(3):256-270

[6] Claerbout J F. Toward a unified theory of reflector mapping[J]. Geophysics,1971,36(3):467-481

[7] 劉紅偉,劉洪,鄒振,等. 地震疊前逆時(shí)偏移中的去噪與存儲(chǔ)[J]. 地球物理學(xué)報(bào),2010,53(9):2171-2180 Liu H W,Liu H,Zou Z,et al. The problems of denoise and storage in seismic reverse time migration[J]. Chinese Journal of Geophysics (in Chinese),2010,53(9):2171-2180

[8] Baysal E,Kosloff D,Sherwood J. A two-way nonreflecting wave equation[J]. Geophysics,1984,49(2):132-141

[9] Loewenthal D,Stoffa P A,Faria E L. Suppressing the unwanted reflections of the full wave equation[J]. Geophysics,1987,52(7):1007-1012

[10] Fletcher R F,Fowler,Kitchenside P. Suppressing artifacts in prestack reverse-time migration[J]. Expanded Abstracts of 75thAnnual Internat SEG Mtg,2005,2049-2051

[11] Kaelin B,Guitton A. Imaging condition for reverse time migration[J]. Expanded Abstracts of 76thAnnual Internat SEG Mtg,2006,2594-2598

[12] Mulder W,Plessix R. One-way and two-way wave-equation migration[J]. Expanded Abstracts of 73rdAnnual Internat SEG Mtg,2003,881-884

[13] Antoine Guitton,Bruno Kaelin,Biondo Biond. Least-square attenuation of reverse time migration artifacts[J]. Expanded Abstracts of 76thAnnual Internat SEG Mtg,2006,2348-2352

[14] Zhang Y,Sun J. Practical issues of reverse time migration:true amplitude gathers,noise removal and harmonic-source encoding[J]. First Break,2009,26(1):29-35

[15] Buades A,Coll B,Morel J M. A review of image denoising algorithms,with a new one[J]. Multiscale Modeling and Simulation,2005(4):490-530

[16] Buades A,Coll B,Morel J M. Image denoising methods:a new nonlocal principle[J].SIAM Review,2010,52(1):113-147

[17] Sheng B,Li P,Sun H. Image-based material restyling with fast non-local means filtering[J]. Proceedings of the 2009 Fifth International Conference on Image and Graphics,2009,841-846

[18] Mahmoudi M,Sapiro G. Fast image and video denoising via nonlocal means of similar neighborhoods[J]. IEEE Signal Processing Letters,2005,12,839-842

[19] Brox T,Kleinschmidt O,Cremers D. Efficient nonlocal means for denoising of textural patterns[J]. IEEE Transactions on Image Processing,2008,17(7)：1083-1092

[20] Coupé P,Yger P,Prima S,et al. An optimized blockwise nonlocal means denoising filter for 3-D magnetic resonance images[J]. IEEE Transactions on Medical Imaging,2008,27(4):425-441

[21] Deledalle C A,Denis L,Tupin F. NL-InSAR:nonlocal interferogram estimation[J]. IEEE Transactions on Geoscience and Remote Sensing,2011,49(4):1441-1452

[22] Zoican S. Speech de-noising system with nonlocal means algorithm[J]. 9thInternational Symposium on Electronics and Telecommunications(ISETC),2010,315-318

[23] Wei D Y,Yin C C. An optimized locally adaptive non-local means denoising filter for cryo-electron microscopy data[J]. Journal of Structural Biology,2010,172(1):211-218

[24] Bonar D,Sacchi M. Denoising seismic data using the nonlocal means algorithm[J]. Geophysics,2012,77(1):A5-A8

(編輯：朱文杰)

Non-local means filtering denoising approach and its application in VSP reverse time migration seismic data

Wang Weihong,Guo Xuebao,Shi Ying

(CollegeofEarthSciences,NortheastPetroleumUniversity,Daqing163318,China)

Low-frequency noise is the inherent property of reverse-time migration (RTM) based on two-way wave equation.Filtering after imaging is widely used to suppress this kind of noise,especially Laplacian filtering method.But the result of Laplacian method depends on angle parameters,the case of rough events and residual noises always exist.In order to improve signal to noise ratio of imaging data of reverse-time migration,Non-local Means (NLM) filtering approach is introduced in this paper,which is used to suppress noise based on the characteristics of event structure of post-stack seismic data,and is applied to the seismic data after Laplacian filtering.By using NLM method,the input seismic data can be divided into noise and data,and the similarity coefficient of different imaging points is applied to realize filtering.For 2D complicated theoretical model,we carried out forward modeling and RTM on VSP data.The results demonstrate that the combination of Laplacian and NLM filtering is better than the result of Laplacian filtering.Moreover,the application of actual VSP seismic data also suggests NLM filtering has advantages in computational accuracy,stability and practicability.

non-local means,noise suppression,reverse time migration,vertical seismic profiling (VSP),Laplacian filtering

2014-09-25;改回日期：2014-12-08。

王維紅(1975—)，男，博士，副教授，現(xiàn)主要從事地震資料數(shù)字處理方面的研究。

國家自然科學(xué)基金項(xiàng)目(41474118)和國家高技術(shù)研究發(fā)展計(jì)劃(863)(2012AA061202)資助。

P631

A

1000-1441(2015)02-0165-07

10.3969/j.issn.1000-1441.2015.02.007