羅　娜，王利兵，王　靜，宋　昭，趙永紅，李細(xì)順，賈　華，陳凱男，趙志遠(yuǎn)

(河北省地震局紅山基準(zhǔn)地震臺(tái)，河北　邢臺(tái)　054000)

基于小波模極大值的地震信號(hào)去噪研究

羅娜，王利兵，王靜，宋昭，趙永紅，李細(xì)順，賈華，陳凱男，趙志遠(yuǎn)

(河北省地震局紅山基準(zhǔn)地震臺(tái)，河北邢臺(tái)054000)

摘要：小波分析在時(shí)域和頻域具有很好的局部化特性，是分析和處理數(shù)字信號(hào)強(qiáng)有力的工具。文章將基于小波變換的模極大值去噪算法應(yīng)用到地震信號(hào)的去噪研究中。首先依據(jù)相關(guān)理論驗(yàn)證算法的有效性，并對(duì)紅山基準(zhǔn)臺(tái)的地震數(shù)據(jù)進(jìn)行去噪分析處理。 結(jié)果表明，去噪后的信號(hào)有效去除了大部分毛刺，去噪效果良好，噪聲得到很好的抑制。為實(shí)現(xiàn)數(shù)據(jù)處理的界面化及易操作，在基于小波變換模極大值去噪算法比較分析的基礎(chǔ)上，設(shè)計(jì)一個(gè)地震信號(hào)去噪系統(tǒng)分析界面，從而實(shí)現(xiàn)數(shù)據(jù)去噪的可視化應(yīng)用。

關(guān)鍵詞：地震信號(hào)；小波變換；模極大值

0引言

地震數(shù)字化記錄信號(hào)在采集、處理等過(guò)程中，不可避免地疊加了各種干擾信息，從而影響觀測(cè)質(zhì)量，增加震相分析的難度，為地震信號(hào)的精確識(shí)別造成一定的困難。因此，在地震信號(hào)的處理過(guò)程中，提取地震信號(hào)中的有用信息，提高地震資料的信噪比，對(duì)后續(xù)的信號(hào)分析研究是十分重要的[1-2]。

地震信號(hào)屬于非平穩(wěn)信號(hào)，傳統(tǒng)的高通濾波和低通濾波等方法處理非平穩(wěn)信號(hào)具有一定的局限性，信號(hào)各頻段的噪聲不能有效去除，且信號(hào)細(xì)節(jié)信息得不到有效地保留[3]。上世紀(jì)80年代中后期發(fā)展并成熟起來(lái)的小波理論，由于其多尺度分解的“顯微鏡”的特性，可以展示感興趣頻段和時(shí)間的精細(xì)信息，在地震、電磁、形變等學(xué)科的分析和研究中得到了越來(lái)越廣泛的應(yīng)用[4-6]。文章利用基于小波分析的模極大值去噪方法，對(duì)紅山基準(zhǔn)地震臺(tái)(以下簡(jiǎn)稱紅山臺(tái))記錄的地震數(shù)據(jù)進(jìn)行去噪研究，對(duì)原始記錄信號(hào)進(jìn)行小波分解后，利用有用信號(hào)和噪聲的小波系數(shù)在各個(gè)尺度上具有的不同特性，去除屬于噪聲的小波系數(shù)，然后對(duì)處理后的小波系數(shù)進(jìn)行重構(gòu)，得到降噪后的信號(hào)，并對(duì)去噪前后的數(shù)據(jù)進(jìn)行對(duì)比分析，檢驗(yàn)去噪效果，設(shè)計(jì)地震信號(hào)去噪系統(tǒng)分析界面，實(shí)現(xiàn)數(shù)據(jù)的可視化應(yīng)用。

1小波模極大值去噪算法

1.1　小波變換理論

設(shè)ψ(t)為一平方可積函數(shù)，即ψ(t)∈L2(R)，若其傅里葉變換滿足條件[7-8]：

(1)

式中：Ψ(w)為ψ(t)的傅立葉變換，稱ψ(t)為一個(gè)基本小波或小波母函數(shù)，稱式(1)為小波函數(shù)的可容許條件，且母小波是不唯一的，是可選擇的。

假設(shè)信號(hào)f(t)∈L2(R)，則其連續(xù)小波變換定義為：

(2)

若ψa,τ(t)滿足容許條件(1)，則連續(xù)小波存在逆變換，其相應(yīng)的反變換公式為：

(3)

1.2　基于小波變換模極大值的去噪算法

信號(hào)f(t)的小波變換與李氏指數(shù)滿足如下關(guān)系[9-12]：

log2|W2jf(t)|≤log2k+jα，

(4)

式中：jα這一項(xiàng)把小波變換的尺度特征j與李氏指數(shù)(Lipschitz exponent)α聯(lián)系了起來(lái)。由式(4)可知，當(dāng)函數(shù)f(t)的李氏指數(shù)α>0時(shí)，小波變換的模極大值將隨著尺度j的增大而增大；當(dāng)α<0時(shí)，則隨著尺度j的增大而減小。對(duì)階躍情況α=0，則小波變換的模極大值不隨尺度j的改變而改變。地震記錄中反射波是連續(xù)可導(dǎo)的，也是有限的不連續(xù)點(diǎn)值，則該點(diǎn)的α值滿足α≥0，因而，有效信號(hào)的小波變換模極大值隨著尺度的增大而增大。然而，噪聲所對(duì)應(yīng)的李氏指數(shù)α往往是小于0的，其對(duì)應(yīng)的小波變換模極大值隨著尺度的增大而減小。因此，連續(xù)做若干次小波變換之后，由噪聲對(duì)應(yīng)的模極大值已經(jīng)基本去除或幅值很小，而所余模極大值點(diǎn)主要由信號(hào)控制?；诖嗽?，有如下去噪算法。

(1) 對(duì)含噪信號(hào)進(jìn)行二進(jìn)小波變換，求出每一尺度上小波變換細(xì)節(jié)系數(shù)的模極大值，一般分解尺度為3~5個(gè)尺度。

(2) 選取最大尺度上(設(shè)為J)的模極值點(diǎn)，選擇一閾值，若模極值點(diǎn)對(duì)應(yīng)的幅值小于此閾值，則去掉該極值點(diǎn)；否則，對(duì)應(yīng)的模極值點(diǎn)予以保留。

(3) 在尺度j-1上尋找尺度j(3≤j≤J)上小波變換模極大值點(diǎn)的傳播點(diǎn)，即保留j-1尺度上有效信號(hào)產(chǎn)生的極值點(diǎn)，去掉由噪聲產(chǎn)生的極值點(diǎn)。從理論上講，須在一個(gè)稠密的尺度序列上計(jì)算小波變換，從而找到那些模極值點(diǎn)傳播到下一尺度。但在實(shí)際的資料處理中，可以利用一個(gè)簡(jiǎn)單ad hoc算法來(lái)搜索那些模極大值點(diǎn)傳播到下一尺度。具體算法為：在尺度上j的模極大值點(diǎn)位置，構(gòu)造一個(gè)鄰域O(nji,εj)，其中nji為尺度j上的第i個(gè)極值點(diǎn)，εj為僅與尺度j有關(guān)的常數(shù)。在尺度j-1上的模極大值點(diǎn)中保留落在每一個(gè)鄰域O(nji，εj)上的模極大值點(diǎn)，而去除落在鄰域外面的模極值點(diǎn)，從而得到j(luò)-1尺度上新的極值點(diǎn)。然后令j=j-1，重復(fù)上述步驟，逐級(jí)搜索，到j(luò)=2為止。

(4) 對(duì)尺度j=1，首先將尺度j=1上的小波變換模極大值全部去掉，然后在尺度j=2上存在模極大值點(diǎn)的位置上，由尺度j=2上的模極大值進(jìn)行非線性插值得到尺度j=1上的模極大值。每個(gè)模極大值點(diǎn)的Lipschitz指數(shù)α由式(5)決定。

W2jf(x)≤K(2j)α。

(5)

對(duì)于2個(gè)尺度j和i+1，對(duì)上式取對(duì)數(shù)并做相減可得：

(6)

當(dāng)j≥2時(shí)，就可由上式近似求出2α。令式(6)中j=1可得：

(7)

從而在與尺度j=2上模極大值點(diǎn)位置相同的位置上，計(jì)算出尺度j=1上的模極大值。

(5) 根據(jù)每一尺度上保留下來(lái)的小波變換模極大值，對(duì)信號(hào)進(jìn)行重構(gòu)。

2數(shù)據(jù)處理與分析

2.1　測(cè)試信號(hào)仿真實(shí)驗(yàn)與分析

利用matlab編制相關(guān)程序，并對(duì)一測(cè)試含噪信號(hào)heavy sine作小波變換模極大值去噪實(shí)驗(yàn)，驗(yàn)證算法的有效性。在heavy sine信號(hào)中加入隨機(jī)高斯白噪聲，采用db4小波對(duì)信號(hào)進(jìn)行分解與重構(gòu)，原始信號(hào)加入噪聲后信噪比為15.694，分解尺度為4。

圖1為原始信號(hào)heavy sine、高斯白噪聲和加入高斯白噪聲后的heavy sine信號(hào)，圖2為經(jīng)過(guò)去噪處理后的信號(hào)，把去噪重構(gòu)后的信號(hào)和原始信號(hào)heavy sine做對(duì)比分析看出，原始信號(hào)heavy sine基本被恢復(fù)，有效降低了噪聲的影響，去噪后信號(hào)的信噪比由SNR為15.694 5提升到23.915 9，信號(hào)信噪比有明顯的提高。

圖1　原始信號(hào)、高斯白噪聲和加噪信號(hào)圖Fig.1　The original signal, Gaussian white noise and signal with noise

圖2　經(jīng)過(guò)去噪處理后的信號(hào)Fig.2　Signal after denoising

2.2　實(shí)際地震資料的小波變換模極大值去噪處理

2.2.1資料選取

研究資料選取河北紅山臺(tái)記錄的地震數(shù)據(jù)，數(shù)據(jù)接收采用中國(guó)地震局武漢地震研究所生產(chǎn)的CTS-1寬頻帶地震計(jì)，數(shù)據(jù)采集系統(tǒng)采用中國(guó)地震局港震機(jī)電公司生產(chǎn)的EDAS-24地震數(shù)據(jù)采集器?！熬盼濉睌?shù)據(jù)采樣頻率為50 Hz，波形選取時(shí)查閱地震報(bào)告，根據(jù)報(bào)告調(diào)取地震波形作為研究數(shù)據(jù)。

2.2.2資料處理

(1) 選取紅山臺(tái)記錄的2012年10月2日13點(diǎn)26分發(fā)生在中國(guó)東南部4.0級(jí)的原始地震數(shù)據(jù)進(jìn)行分析處理，長(zhǎng)度為120 s左右，采用db4小波分解，尺度為4，處理結(jié)果如圖3所示。

圖3　原始信號(hào)小波4個(gè)尺度分解Fig.3　Original signal wavelet decomposed into 4 grades

圖3為原始信號(hào)201210021326.EVT 4個(gè)尺度分解，a1~a4為近似分解部分，d1~d4為細(xì)節(jié)分解部分，可以看出原始信號(hào)被噪聲干擾，噪聲基本覆蓋了整個(gè)有效信號(hào)區(qū)間，PG、SG波列不清晰，尤其是Pg波基本被淹沒(méi)在噪聲中，在細(xì)節(jié)尺度分解上存在著大部分噪聲。圖4為原始信號(hào)和去噪后信號(hào)的fft變換，把去噪重構(gòu)后信號(hào)和原始信號(hào)進(jìn)行對(duì)比分析可以看出，去噪處理后的信號(hào)信噪比有明顯的提高，而經(jīng)過(guò)去噪后的fft變換高頻成分明顯減少，幅值大幅降低。去噪后的信號(hào)有效去除了大部分毛刺，噪聲得到很好的抑制，原始信號(hào)201210021326.EVT基本被恢復(fù)，PG、SG波列變得清晰，到時(shí)可以得到準(zhǔn)確的辨識(shí)。

(2) 2011年，由于安裝電擾動(dòng)儀器需要，在臺(tái)站實(shí)施打井作業(yè)，對(duì)數(shù)據(jù)造成一定的干擾。選取在此作業(yè)期間的波形數(shù)據(jù)，進(jìn)行小波變換模極大值去噪處理，采樣點(diǎn)數(shù)為4 096，分析結(jié)果如圖5、第15頁(yè)圖6所示。

從圖5可以看出，地震信號(hào)基本淹沒(méi)在打井噪聲信號(hào)中， fft變換圖也可以看出原始信號(hào)含有大量高頻信息量。經(jīng)過(guò)去噪處理后，打井噪聲基本被去除，波形質(zhì)量有明顯改善，地震信號(hào)基本被恢復(fù)，原始fft變換的高頻成分幅值明顯降低，說(shuō)明基本達(dá)到了去噪的效果。從圖中還可以看出，去噪后的信號(hào)有的細(xì)節(jié)地方噪聲處理的效果并不理想，所以，對(duì)于尺度和閾值的選取還需做進(jìn)一步的研究。

圖4　原始信號(hào)和去噪后信號(hào)的fft變換Fig.4　FFT transforms of both original signal and signal after denoising

圖5　原始信號(hào)和去噪后的信號(hào)Fig.5　Original signal and signal after denoising

圖6　原始信號(hào)和去噪后信號(hào)的fft變換Fig.6　FFT transforms of both original signal and signal after denoising

3基于小波模極大值去噪系統(tǒng)的分析與設(shè)計(jì)

3.1　系統(tǒng)組成

在基于小波變換模極大值去噪算法比較分析的基礎(chǔ)上，設(shè)計(jì)了一個(gè)地震信號(hào)去噪系統(tǒng)分析界面，實(shí)現(xiàn)數(shù)據(jù)去噪的可視化應(yīng)用。系統(tǒng)主窗口主要由數(shù)據(jù)文件名、數(shù)據(jù)載入、fft變換、模極大值去噪、退出5個(gè)模塊組成。

3.2　系統(tǒng)界面設(shè)計(jì)及運(yùn)行

在系統(tǒng)的存放目錄下直接雙擊“waveFormAnalysis3.fig”文件運(yùn)行系統(tǒng)，或在matlab命令窗口中輸入“waveFormAnalysis3.fig”回車啟動(dòng)系統(tǒng)，運(yùn)行系統(tǒng)后，點(diǎn)擊數(shù)據(jù)文件名，可以選擇文件路徑。數(shù)據(jù)載入之后，可以進(jìn)行去噪處理(見圖7)。

圖7　小波模極大值去噪Fig.7　Denoising based on wavelet transform modulus maximum

如果系統(tǒng)計(jì)算量很大，數(shù)據(jù)處理時(shí)間太長(zhǎng)，則運(yùn)行速率就會(huì)降低，甚至模塊暫時(shí)出現(xiàn)運(yùn)行錯(cuò)誤，從而導(dǎo)致死機(jī)情況。因此，建議數(shù)據(jù)保存后，及時(shí)清理存儲(chǔ)空間的數(shù)據(jù)，但是總體來(lái)看，基本達(dá)到預(yù)期效果，實(shí)現(xiàn)了可視化操作，具有一定的實(shí)用性。

4結(jié)論

小波具有良好的時(shí)頻局部化和多分辨率分析能力，由于有效信號(hào)的小波變換模極大值隨著尺度的增加而增加，而噪聲的小波變換模極大值隨著尺度的增加而減小，可以利用有效信號(hào)和噪聲模極大值隨尺度變化的不同特征來(lái)去除噪聲。文中將基于小波變換的模極大值算法用于測(cè)試含噪信號(hào)heavy sine和紅山臺(tái)地震數(shù)據(jù)的去噪處理中，噪聲得到了抑制，有效信號(hào)得到保留，可以基本保持信號(hào)的不失真，從而證明算法的有效性。總體來(lái)看，基本達(dá)到預(yù)期去噪的目的，有實(shí)際的應(yīng)用效果。

參考文獻(xiàn)：

[1]魏紅梅，黃世源，許飛.小波包去噪在地震信號(hào)預(yù)處理中的應(yīng)用[J].東北地震研究,2008,24(2)：45-49.

[2]李文，劉霞，段玉波，等.基于小波熵與相關(guān)性相結(jié)合的小波模極大值地震信號(hào)去噪[J].地震學(xué)報(bào)，2012,34(6)：841-850.

[3]徐宏斌，李庶林，陳際經(jīng).基于小波變換的大尺度巖體結(jié)構(gòu)微震監(jiān)測(cè)信號(hào)去噪方法研究[J].地震學(xué)報(bào),2012,34(1)：85-96.

[4]孫志新.小波分析在地震數(shù)據(jù)噪聲處理中的應(yīng)用研究[D].中國(guó)地震局工程力學(xué)研究所,2008:25-35.

[5]李偉，馬欽忠，宋志平，等.小波變換在地電場(chǎng)數(shù)據(jù)分析中的應(yīng)用[J].地震學(xué)報(bào)，2013，35(1)：26-35.

[6]李杰，劉希強(qiáng)，李紅，等.利用小波變換方法分析形變觀測(cè)資料的正常背景變化特征[J].地震學(xué)報(bào),2005，27(1)：33-41.

[7]彭玉華.小波變換與工程應(yīng)用[M].北京：科學(xué)出版社,1999.

[8]張德豐.MATLAB小波分析第2版[M].北京：機(jī)械工業(yè)出版社,2011.

[9]田金文，高謙，杜擁軍.基于小波變換模極大值的地震信號(hào)去噪處理方法[J].江漢石油學(xué)院學(xué)報(bào)，2001,23(1)：22-25.

[10]姚勝利.地震信號(hào)的小波去噪方法研究[D].長(zhǎng)沙：中南大學(xué),2007:22-37.

[11]孔祥茜，吳繼偉，岳繼光.地震信號(hào)小波變換的去噪方法[J].計(jì)算機(jī)輔助工程，2005,14(3)：52-56.

[12]楊建國(guó).小波分析及其工程應(yīng)用[M].北京：機(jī)械工業(yè)出版社，2005.

Denoising of Seismic Signals Based on Wavelet Transform Modulus Maximum

LUO Na, WANG Li-bing, WANG Jing, SONG Zhao, ZHAO Yong-hong, LI Xi-shun, JIA Hua, CHEN Kai-nan, ZHAO Zhi-yuan

(Hongshan Reference Seismological Station of Earthquake Administration of Hebei Province, Xingtai, Hebei 054000, China)

Abstract:Wavelet transform analysis is of good localization properties both in time and frequency domain and is a powerful tool to analyze and process digital signal. Denoising algorithm based on wavelet transform modulus maximum is applied to denoising of seismic signals in the paper. The effectiveness of the algorithm is verified and then denoising is carried on seismic data of Hongshan Reference Station. The results show that the most of the burrs in signals are removed after denoising and noise is suppressed effectively. In order to make the data processing more easily, an analysis interface for seismic signals denoising system is designed to implement visualization of data denoising based on comparison of denoising algorithm of wavelet transform modulus maximum.

Key words:Seismic signals; Wavelet transform; Modulus maximum

作者簡(jiǎn)介：第一羅娜(1981—)，女，河北省保定人。2009年畢業(yè)于蘭州理工大學(xué)，碩士研究生，工程師。

基金項(xiàng)目：河北省地震局碩(博)預(yù)研究項(xiàng)目(201201)。

收稿日期：2015-01-05

中圖分類號(hào)：P315.3

文獻(xiàn)標(biāo)志碼：A

文章編號(hào)：1000-6265(2015)02-0012-04