基于EMD與1D全變分的地震信號(hào)去噪

王均榮，張 濤，安素珍，周仲禮，王茂芝

(1.成都理工大學(xué)數(shù)學(xué)地質(zhì)四川省重點(diǎn)實(shí)驗(yàn)室，成都 610059;2.電子科技大學(xué)神經(jīng)信息教育部重點(diǎn)實(shí)驗(yàn)室，成都 610054)

引 言

地震信號(hào)中噪聲的壓制是地震數(shù)據(jù)處理中重要的預(yù)處理任務(wù)，關(guān)系到后續(xù)的地震信號(hào)屬性提取和數(shù)據(jù)解釋等相關(guān)工作。嚴(yán)格意義上地震信號(hào)屬于非平穩(wěn)信號(hào)［1］，隨著現(xiàn)代地震勘探對(duì)地震資料信噪比的要求越來越高，各種現(xiàn)代數(shù)字信號(hào)處理方法被廣泛應(yīng)用，如ICA、PCA、小波閾值、EMD、TV等去噪方法。但是，目前很多方法難以有效地應(yīng)付更為復(fù)雜的地震信號(hào)，往往在去噪的同時(shí)，也會(huì)造成嚴(yán)重的細(xì)節(jié)性有效信號(hào)的丟失。因此，在有效消除地震信號(hào)噪聲的同時(shí)，如何有效保持地震信號(hào)中的構(gòu)造邊緣信息，已成為人們更為關(guān)心的問題。經(jīng)驗(yàn)?zāi)B(tài)分解時(shí)頻分析法是最近發(fā)展起來的一種處理非線性非平穩(wěn)信號(hào)的新方法，由于EMD算法能自主分解輸入信號(hào)，因此EMD去噪是自適應(yīng)的，不需要先驗(yàn)知識(shí)，與傳統(tǒng)的傅里葉變換濾波相比有著明顯的優(yōu)勢(shì)［2］。另外，1D-TV算法能較好地實(shí)現(xiàn)地震信號(hào)與噪聲信號(hào)的分離，并能有效的保護(hù)地震信號(hào)邊緣信息［3-4］。

將EMD技術(shù)與1D-TV方法相結(jié)合，應(yīng)用于地震信號(hào)噪聲的壓制。首先通過EMD算法將地震信號(hào)分解為一系列表征信號(hào)特征時(shí)間尺度的固有模態(tài)函數(shù)(Intrinsic Mode Function，IMF)，然后應(yīng)用1D-TV方法對(duì)所選擇的高頻成分IMF分量進(jìn)行處理，最后利用處理后的IMF分量重構(gòu)信號(hào)，從而達(dá)到去噪的目的。該方法能充分發(fā)揮兩者的優(yōu)點(diǎn)，在更精細(xì)的尺度上對(duì)噪聲進(jìn)行壓制，并使信號(hào)邊緣信息得到有效的保護(hù)。實(shí)際的地震數(shù)據(jù)處理結(jié)果表明該新算法是可行的，且去噪效果優(yōu)于單獨(dú)使用EMD方法。

1 EMD/TV(全變分)方法原理

1.1 經(jīng)驗(yàn)?zāi)B(tài)分解(EMD)去噪原理

EMD是對(duì)非平穩(wěn)信號(hào)進(jìn)行時(shí)頻處理的方法，該方法將信號(hào)中不同尺度的波動(dòng)逐級(jí)分解開來，產(chǎn)生一系列具有不同特征尺度的數(shù)據(jù)信號(hào)，每一個(gè)數(shù)據(jù)信號(hào)就是一個(gè)固有模態(tài)函數(shù)［5］IMF，IMF分量按從高頻到低頻依次排列。其中，IMF必須滿足以下兩個(gè)條件:

(1)信號(hào)的極值點(diǎn)數(shù)目和過零點(diǎn)數(shù)目相等或最多相差一個(gè)。

(2)且由局部極大值構(gòu)成的上包絡(luò)線和局部極小值構(gòu)成的下包絡(luò)線的平均值為零［6］。

以實(shí)際地震信號(hào)分解(圖1)為例，EMD算法的實(shí)施步驟:

(1)找到信號(hào)y(t)所有的局部極大值并且用三次樣條函數(shù)插值連接，獲得上包絡(luò)線［同理，連接局部極小值點(diǎn)作為下包絡(luò)線。

(2)上下包絡(luò)的平均值組成均值信號(hào)m，并定義信號(hào)h1=y(t)-m。

(3)如果信號(hào)h1滿足IMF的條件，那么h1就是求得的第一個(gè)IMF分量;否則將h1作為原始信號(hào)進(jìn)行(1)～(2)的步驟，直到第k次迭代后的差值成為第一個(gè)IMF，記為 c1。

(4)從原信號(hào)中減去c1，得到第一階剩余信號(hào)r1=y(t)-c1，將r1作為新數(shù)據(jù)序列重復(fù)上述過程，可以遞推得到其他IMF，直到趨勢(shì)分量rn單調(diào)或只有一個(gè)極值。

圖1 原始地震信號(hào)通過EMD分解成不同的IMF成分

假設(shè) y=(y［1］，···，y［N］)∈RN，N ＞ 1，采集到的地震信號(hào)，x=(x［1］，···，x［N］)∈ RN為干凈的地震信號(hào)，n為采集過程中混入的隨機(jī)噪聲，則

其中，x(t)為干凈信號(hào)，n(t)為噪聲信號(hào)，對(duì)y(t)進(jìn)行EMD分解后滿足下式:

其中，ci為第i個(gè)IMF分量，rn+1為剩余分量。

EMD去噪的基本原理是從分解的IMF分量中去除高頻噪聲引起的IMF分量。由于EMD分解的IMF分量從高到低，因此第一個(gè)IMF分量的頻率最高，是由高頻隨機(jī)噪聲引起的。去除高頻分量后，用剩余的IMF分量重構(gòu)信號(hào)，即達(dá)到去噪的目的［5，7］。

1.2 1D－TV模型去噪原理

全變分去噪模型首先由Rudin、Osher and Fatemi于1992年提出［8］，主要用于圖像的去噪、修復(fù)、分割等方面，該方法基于噪聲圖像的全變分明顯高于無噪聲圖像的全變分，把圖像去噪問題轉(zhuǎn)化為求解圖像模型的能量泛函最小化問題，其實(shí)質(zhì)就是各向異性擴(kuò)散，從而實(shí)現(xiàn)去除圖像的噪聲，其優(yōu)點(diǎn)在于平滑噪聲的同時(shí)，又能很好的保持圖像邊緣特性［9-10］。

對(duì)于式(1)給出的噪聲信號(hào)，其全變分的離散形式定義為［3］:

由于噪聲圖像的全變分明顯要大于無噪聲圖像的全變分，因此全變分的信號(hào)去噪問題就轉(zhuǎn)化為最小化全變分泛函問題［12］:

其可以等價(jià)于式(5)，即為求解非線性最優(yōu)化問題［3-4］:

式(5)中，第一項(xiàng)為保真項(xiàng)，它主要起降低原信號(hào)失真度的作用;第二項(xiàng)為正則化項(xiàng)，λ為正則化參數(shù)，λ的取值與信號(hào)的噪聲水平有關(guān)，對(duì)平衡去噪與平滑起重要作用。

本文采用Little［11］提出的原始-內(nèi)點(diǎn)對(duì)偶算法來求解，通過在原始信號(hào)上加上高斯白噪聲進(jìn)行仿真實(shí)驗(yàn)(圖2)。從圖2可以看出，經(jīng)1D-TV處理后，噪聲得到很好的壓制。

2 新算法的構(gòu)造

為了更好地去除地震信號(hào)中的噪聲，同時(shí)保護(hù)地震構(gòu)造邊緣信息，基于1D-TV算法和EMD算法各自的優(yōu)點(diǎn)，提出了一種新的混合去噪算法，其基本思想:首先，通過EMD算法將地震信號(hào)分解為一系列表征信號(hào)特征時(shí)間尺度的固有模態(tài)函數(shù)IMF;然后，應(yīng)用1D-TV方法對(duì)所選擇的高頻成分IMF分量進(jìn)行處理;最后，利用處理后的IMF分量重構(gòu)信號(hào)，從而達(dá)到去噪的目的。

2.1 EMD分解高頻部分確定原則

圖2 1D－TV去噪效果

信號(hào)經(jīng)EMD分解后得到有限個(gè)頻率從高到低的IMF，其中階數(shù)小的IMF對(duì)應(yīng)于信號(hào)的高頻成分，一般包含的是信號(hào)尖銳部分或噪聲，階數(shù)大的IMF對(duì)應(yīng)于信號(hào)的低頻成分，一般認(rèn)為低頻成分中的噪聲影響很小。基于EMD去噪的主要思想就在于，對(duì)大多數(shù)被噪聲污染的信號(hào)，其主要能量集中在低階IMF段，越往高階 IMF段，其含的能量就越?。?，7］。因此，一定存在某個(gè)IMF分量，使得對(duì)于該分量之后的IMF(k+1)中的信號(hào)為主導(dǎo)模態(tài)，而其前k個(gè)IMF中噪聲為主導(dǎo)模態(tài)。

對(duì)于含有噪聲的地震信號(hào)表達(dá)式可以寫成:

式中s(t)為干凈地震信號(hào)，m(t)為噪聲信號(hào)?，F(xiàn)采用噪聲能量最小的準(zhǔn)則重構(gòu)地震信號(hào)，用?(t)表示EMD分解信號(hào)得到的本征模態(tài)函數(shù)重構(gòu)得到的地震信號(hào)。去掉前k-1個(gè)本征模態(tài)函數(shù)后的信號(hào):

利用均方差［13］使?(t)盡量接近s(t)。均方差表示為:

依據(jù)公式(8)可得IMF組中各個(gè)分量的能量:

由此求得能量最小分量對(duì)應(yīng)的k值，則有:

式(11)中，n是重構(gòu)地震信號(hào)的起始IMF索引，從第n個(gè)本征模態(tài)函數(shù)開始，以后的都為信號(hào)的主要成分，從而確定k值。

2.2 新算法的描述

對(duì)地震信號(hào)y(t)去噪混合算法描述如下:

Step1:將地震信號(hào)y(t)進(jìn)行多尺度EMD分解，得到一組從高頻到低頻排列的本征函數(shù)IMF。

Step2:根據(jù)上述能量法則選取前k階IMF作為高頻信號(hào)。

Step3:對(duì)選取的高頻信號(hào)進(jìn)行選取不同的正則化參數(shù)λ進(jìn)行1D-TV去噪。

Step4:對(duì)處理后的高頻IMF成分和剩余IMF成分進(jìn)行重構(gòu)得到去噪后的地震信號(hào)。

其流程如圖3所示。

圖3 EMD－1DTV去噪算法流程圖

3 實(shí)際地震信號(hào)處理

利用新提出的算法，對(duì)某一金屬礦床地震信號(hào)行了處理。圖4為含有噪聲的原始地震剖面圖，此地震數(shù)據(jù)是二進(jìn)制格式，選取了100道，道間距50 m，每道500個(gè)采樣點(diǎn)，采樣間隔為1 ms。圖5為利用新提出的基于EMD+1DTV方法去噪后的地震剖面圖，圖6為傳統(tǒng)EMD方法去噪后的地震剖面圖。

為了評(píng)價(jià)提出的新算法信號(hào)去噪的效果，利用峰值性噪比(PSNR)和邊緣保持度兩個(gè)評(píng)價(jià)標(biāo)準(zhǔn)對(duì)實(shí)際的處理結(jié)果進(jìn)行評(píng)價(jià)［14］:

(1)峰值信噪比:

圖4 原始地震信號(hào)波形圖

圖5 EMD+1DTV方法去噪后的地震信號(hào)波形圖

圖6 EMD方法去噪后的地震波形圖

其中，q(i，j)為原始的地震數(shù)據(jù)，p(i，j)為經(jīng)過處理后的地震數(shù)據(jù);(i，j)為地震數(shù)據(jù)中的某一點(diǎn);M，N代表地震數(shù)據(jù)中點(diǎn)的集合。

(2)邊緣保持度［15］:

其中，EPI為邊緣保持度，p(i，j)，q(i，j)分別為處理后的地震數(shù)據(jù)和原始地震數(shù)據(jù)的某一點(diǎn)。EPI為去噪后的地震數(shù)據(jù)和未經(jīng)去噪的地震數(shù)據(jù)的邊緣對(duì)比度之比，EPI越接近于1，去噪后的地震信號(hào)邊緣保持度就越好。

通過表1可以看出，基于EMD與1DTV的新算法對(duì)金屬地震資料處理結(jié)果要比單獨(dú)的運(yùn)用EMD效果要好，在峰值信噪比和邊緣保持度上都有著很大的提高。

表1 峰值信噪比與邊緣保持度處理結(jié)果對(duì)比表

4 結(jié)束語

EMD和全變分方法都是新興的信號(hào)處理方法，EMD分解對(duì)信號(hào)有良好的自適應(yīng)性，而全變分方法在信號(hào)降噪的過程中能保護(hù)好信號(hào)的邊緣信息，有著廣泛的應(yīng)用，但是二者也有各自不同的缺陷，并不完美，因此把二者有效的結(jié)合具有十分重要的理論與實(shí)踐意義。關(guān)于一維全變分在處理一維信號(hào)上還有很多的限制，因此，對(duì)一維全變分在地震信號(hào)上的處理還值得進(jìn)一步的研究。

［1］(美)L科恩.時(shí)頻分析理論及應(yīng)用［M］.白居憲，譯.西安:西安交通大學(xué)出版社，1998.

［2］楊凱，劉偉.基于改進(jìn)EMD的地震信號(hào)去噪［J］.西南石油大學(xué)學(xué)報(bào):自然科學(xué)版，2012，34(4):75-81.

［3］ Karahanoglu F I，Bayram I，Van D V D.A signal processing approach to generalized 1-D total variation［J］.Signal Processing，IEEE Transactions on，2011，59(11):5265-5274.

［4］ Condat L.A direct algorithm for 1D total variation denoising［J］.IEEE Signal Proc.Letters，2013，20(11):1054-1057.

［5］史恒，李桂林，王偉，等.基于總體經(jīng)驗(yàn)?zāi)Ｊ椒纸獾牡卣鹦盘?hào)隨機(jī)噪聲消除［J］.地球物理學(xué)進(jìn)展，2011，26(1):71-78.

［6］ Huang N E，Shen Z，Long S R，et al.The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis［J］.Proceedings of The Royal Society A Mathematical Physical and Engineering Sciences，1998，454(1971):903-995.

［7］王 婷.EMD算法研究及其在信號(hào)去噪中的應(yīng)用［D］.哈爾濱:哈爾濱工程大學(xué)，2010.

［8］ Rudin L I，Osher S，F(xiàn)atemi E.Nonlinear total variation based noise removal algorithms［J］.Physica D:Nonlinear Phenomena，1992，60(1):259-268.

［9］ Chan T F，Golub G，Mulet P.A nonlinear primal-dual method for total variation-based image restoration［J］.SIAM J.Sci.Comp.，1999，20(6):1964-1977.

［10］ Chambolle A，Lions P L.Image recovery via total variation minimization and related problems［J］.Numer-ische Mathematik，1997，76(2):167-188.

［11］Little M A，Jones N S.Sparse Bayesian step-filtering for high-throughput analysis of molecular machine dynamics［C］//Proceedings of the IEEE International Conference on Acoustics，Speech，and Signal Processing，Dallas，Texas，USA，March 14-19，2010:4162-4165.

［12］ Tang G，Ma J W.Application of Total-Variation-Based curvelet shrinkage for Three-Dimensional seismic data denoising［J］.IEEE geoscience and remote sensing letters，2011，8(1):103-105.

［13］ Boudraa A O，Cexus J C.EMD-Based signal filtering［J］.IEEE Transactions on Instrumentation and easurement，2007，56(6):2196-2202.

［14］屈 勇，曹俊興，朱海東，等.一種改進(jìn)的全變分地震圖像去噪技術(shù)［J］.石油學(xué)報(bào)，2011，32(5):815-819.

［15］鐘勇，孫東，陳磊.圖像增強(qiáng)技術(shù)在地震資料處理中的應(yīng)用［J］.物探與化探，2010，34(3):392-341.