羅 娜,王利兵,王 靜,宋 昭,趙永紅,李細順,賈 華,陳凱男,趙志遠
(河北省地震局紅山基準地震臺,河北 邢臺 054000)
基于小波模極大值的地震信號去噪研究
羅娜,王利兵,王靜,宋昭,趙永紅,李細順,賈華,陳凱男,趙志遠
(河北省地震局紅山基準地震臺,河北邢臺054000)
摘要:小波分析在時域和頻域具有很好的局部化特性,是分析和處理數(shù)字信號強有力的工具。文章將基于小波變換的模極大值去噪算法應(yīng)用到地震信號的去噪研究中。首先依據(jù)相關(guān)理論驗證算法的有效性,并對紅山基準臺的地震數(shù)據(jù)進行去噪分析處理。 結(jié)果表明,去噪后的信號有效去除了大部分毛刺,去噪效果良好,噪聲得到很好的抑制。為實現(xiàn)數(shù)據(jù)處理的界面化及易操作,在基于小波變換模極大值去噪算法比較分析的基礎(chǔ)上,設(shè)計一個地震信號去噪系統(tǒng)分析界面,從而實現(xiàn)數(shù)據(jù)去噪的可視化應(yīng)用。
關(guān)鍵詞:地震信號;小波變換;模極大值
0引言
地震數(shù)字化記錄信號在采集、處理等過程中,不可避免地疊加了各種干擾信息,從而影響觀測質(zhì)量,增加震相分析的難度,為地震信號的精確識別造成一定的困難。因此,在地震信號的處理過程中,提取地震信號中的有用信息,提高地震資料的信噪比,對后續(xù)的信號分析研究是十分重要的[1-2]。
地震信號屬于非平穩(wěn)信號,傳統(tǒng)的高通濾波和低通濾波等方法處理非平穩(wěn)信號具有一定的局限性,信號各頻段的噪聲不能有效去除,且信號細節(jié)信息得不到有效地保留[3]。上世紀80年代中后期發(fā)展并成熟起來的小波理論,由于其多尺度分解的“顯微鏡”的特性,可以展示感興趣頻段和時間的精細信息,在地震、電磁、形變等學(xué)科的分析和研究中得到了越來越廣泛的應(yīng)用[4-6]。文章利用基于小波分析的模極大值去噪方法,對紅山基準地震臺(以下簡稱紅山臺)記錄的地震數(shù)據(jù)進行去噪研究,對原始記錄信號進行小波分解后,利用有用信號和噪聲的小波系數(shù)在各個尺度上具有的不同特性,去除屬于噪聲的小波系數(shù),然后對處理后的小波系數(shù)進行重構(gòu),得到降噪后的信號,并對去噪前后的數(shù)據(jù)進行對比分析,檢驗去噪效果,設(shè)計地震信號去噪系統(tǒng)分析界面,實現(xiàn)數(shù)據(jù)的可視化應(yīng)用。
1小波模極大值去噪算法
設(shè)ψ(t)為一平方可積函數(shù),即ψ(t)∈L2(R),若其傅里葉變換滿足條件[7-8]:
(1)
式中:Ψ(w)為ψ(t)的傅立葉變換,稱ψ(t)為一個基本小波或小波母函數(shù),稱式(1)為小波函數(shù)的可容許條件,且母小波是不唯一的,是可選擇的。
假設(shè)信號f(t)∈L2(R),則其連續(xù)小波變換定義為:
(2)
若ψa,τ(t)滿足容許條件(1),則連續(xù)小波存在逆變換,其相應(yīng)的反變換公式為:
(3)
信號f(t)的小波變換與李氏指數(shù)滿足如下關(guān)系[9-12]:
log2|W2jf(t)|≤log2k+jα,
(4)
式中:jα這一項把小波變換的尺度特征j與李氏指數(shù)(Lipschitz exponent)α聯(lián)系了起來。由式(4)可知,當函數(shù)f(t)的李氏指數(shù)α>0時,小波變換的模極大值將隨著尺度j的增大而增大;當α<0時,則隨著尺度j的增大而減小。對階躍情況α=0,則小波變換的模極大值不隨尺度j的改變而改變。地震記錄中反射波是連續(xù)可導(dǎo)的,也是有限的不連續(xù)點值,則該點的α值滿足α≥0,因而,有效信號的小波變換模極大值隨著尺度的增大而增大。然而,噪聲所對應(yīng)的李氏指數(shù)α往往是小于0的,其對應(yīng)的小波變換模極大值隨著尺度的增大而減小。因此,連續(xù)做若干次小波變換之后,由噪聲對應(yīng)的模極大值已經(jīng)基本去除或幅值很小,而所余模極大值點主要由信號控制?;诖嗽?,有如下去噪算法。
(1) 對含噪信號進行二進小波變換,求出每一尺度上小波變換細節(jié)系數(shù)的模極大值,一般分解尺度為3~5個尺度。
(2) 選取最大尺度上(設(shè)為J)的模極值點,選擇一閾值,若模極值點對應(yīng)的幅值小于此閾值,則去掉該極值點;否則,對應(yīng)的模極值點予以保留。
(3) 在尺度j-1上尋找尺度j(3≤j≤J)上小波變換模極大值點的傳播點,即保留j-1尺度上有效信號產(chǎn)生的極值點,去掉由噪聲產(chǎn)生的極值點。從理論上講,須在一個稠密的尺度序列上計算小波變換,從而找到那些模極值點傳播到下一尺度。但在實際的資料處理中,可以利用一個簡單ad hoc算法來搜索那些模極大值點傳播到下一尺度。具體算法為:在尺度上j的模極大值點位置,構(gòu)造一個鄰域O(nji,εj),其中nji為尺度j上的第i個極值點,εj為僅與尺度j有關(guān)的常數(shù)。在尺度j-1上的模極大值點中保留落在每一個鄰域O(nji,εj)上的模極大值點,而去除落在鄰域外面的模極值點,從而得到j(luò)-1尺度上新的極值點。然后令j=j-1,重復(fù)上述步驟,逐級搜索,到j(luò)=2為止。
(4) 對尺度j=1,首先將尺度j=1上的小波變換模極大值全部去掉,然后在尺度j=2上存在模極大值點的位置上,由尺度j=2上的模極大值進行非線性插值得到尺度j=1上的模極大值。每個模極大值點的Lipschitz指數(shù)α由式(5)決定。
W2jf(x)≤K(2j)α。
(5)
對于2個尺度j和i+1,對上式取對數(shù)并做相減可得:
(6)
當j≥2時,就可由上式近似求出2α。令式(6)中j=1可得:
(7)
從而在與尺度j=2上模極大值點位置相同的位置上,計算出尺度j=1上的模極大值。
(5) 根據(jù)每一尺度上保留下來的小波變換模極大值,對信號進行重構(gòu)。
2數(shù)據(jù)處理與分析
利用matlab編制相關(guān)程序,并對一測試含噪信號heavy sine作小波變換模極大值去噪實驗,驗證算法的有效性。在heavy sine信號中加入隨機高斯白噪聲,采用db4小波對信號進行分解與重構(gòu),原始信號加入噪聲后信噪比為15.694,分解尺度為4。
圖1為原始信號heavy sine、高斯白噪聲和加入高斯白噪聲后的heavy sine信號,圖2為經(jīng)過去噪處理后的信號,把去噪重構(gòu)后的信號和原始信號heavy sine做對比分析看出,原始信號heavy sine基本被恢復(fù),有效降低了噪聲的影響,去噪后信號的信噪比由SNR為15.694 5提升到23.915 9,信號信噪比有明顯的提高。
圖1 原始信號、高斯白噪聲和加噪信號圖Fig.1 The original signal, Gaussian white noise and signal with noise
圖2 經(jīng)過去噪處理后的信號Fig.2 Signal after denoising
2.2.1資料選取
研究資料選取河北紅山臺記錄的地震數(shù)據(jù),數(shù)據(jù)接收采用中國地震局武漢地震研究所生產(chǎn)的CTS-1寬頻帶地震計,數(shù)據(jù)采集系統(tǒng)采用中國地震局港震機電公司生產(chǎn)的EDAS-24地震數(shù)據(jù)采集器?!熬盼濉睌?shù)據(jù)采樣頻率為50 Hz,波形選取時查閱地震報告,根據(jù)報告調(diào)取地震波形作為研究數(shù)據(jù)。
2.2.2資料處理
(1) 選取紅山臺記錄的2012年10月2日13點26分發(fā)生在中國東南部4.0級的原始地震數(shù)據(jù)進行分析處理,長度為120 s左右,采用db4小波分解,尺度為4,處理結(jié)果如圖3所示。
圖3 原始信號小波4個尺度分解Fig.3 Original signal wavelet decomposed into 4 grades
圖3為原始信號201210021326.EVT 4個尺度分解,a1~a4為近似分解部分,d1~d4為細節(jié)分解部分,可以看出原始信號被噪聲干擾,噪聲基本覆蓋了整個有效信號區(qū)間,PG、SG波列不清晰,尤其是Pg波基本被淹沒在噪聲中,在細節(jié)尺度分解上存在著大部分噪聲。圖4為原始信號和去噪后信號的fft變換,把去噪重構(gòu)后信號和原始信號進行對比分析可以看出,去噪處理后的信號信噪比有明顯的提高,而經(jīng)過去噪后的fft變換高頻成分明顯減少,幅值大幅降低。去噪后的信號有效去除了大部分毛刺,噪聲得到很好的抑制,原始信號201210021326.EVT基本被恢復(fù),PG、SG波列變得清晰,到時可以得到準確的辨識。
(2) 2011年,由于安裝電擾動儀器需要,在臺站實施打井作業(yè),對數(shù)據(jù)造成一定的干擾。選取在此作業(yè)期間的波形數(shù)據(jù),進行小波變換模極大值去噪處理,采樣點數(shù)為4 096,分析結(jié)果如圖5、第15頁圖6所示。
從圖5可以看出,地震信號基本淹沒在打井噪聲信號中, fft變換圖也可以看出原始信號含有大量高頻信息量。經(jīng)過去噪處理后,打井噪聲基本被去除,波形質(zhì)量有明顯改善,地震信號基本被恢復(fù),原始fft變換的高頻成分幅值明顯降低,說明基本達到了去噪的效果。從圖中還可以看出,去噪后的信號有的細節(jié)地方噪聲處理的效果并不理想,所以,對于尺度和閾值的選取還需做進一步的研究。
圖4 原始信號和去噪后信號的fft變換Fig.4 FFT transforms of both original signal and signal after denoising
圖5 原始信號和去噪后的信號Fig.5 Original signal and signal after denoising
圖6 原始信號和去噪后信號的fft變換Fig.6 FFT transforms of both original signal and signal after denoising
3基于小波模極大值去噪系統(tǒng)的分析與設(shè)計
在基于小波變換模極大值去噪算法比較分析的基礎(chǔ)上,設(shè)計了一個地震信號去噪系統(tǒng)分析界面,實現(xiàn)數(shù)據(jù)去噪的可視化應(yīng)用。系統(tǒng)主窗口主要由數(shù)據(jù)文件名、數(shù)據(jù)載入、fft變換、模極大值去噪、退出5個模塊組成。
在系統(tǒng)的存放目錄下直接雙擊“waveFormAnalysis3.fig”文件運行系統(tǒng),或在matlab命令窗口中輸入“waveFormAnalysis3.fig”回車啟動系統(tǒng),運行系統(tǒng)后,點擊數(shù)據(jù)文件名,可以選擇文件路徑。數(shù)據(jù)載入之后,可以進行去噪處理(見圖7)。
圖7 小波模極大值去噪Fig.7 Denoising based on wavelet transform modulus maximum
如果系統(tǒng)計算量很大,數(shù)據(jù)處理時間太長,則運行速率就會降低,甚至模塊暫時出現(xiàn)運行錯誤,從而導(dǎo)致死機情況。因此,建議數(shù)據(jù)保存后,及時清理存儲空間的數(shù)據(jù),但是總體來看,基本達到預(yù)期效果,實現(xiàn)了可視化操作,具有一定的實用性。
4結(jié)論
小波具有良好的時頻局部化和多分辨率分析能力,由于有效信號的小波變換模極大值隨著尺度的增加而增加,而噪聲的小波變換模極大值隨著尺度的增加而減小,可以利用有效信號和噪聲模極大值隨尺度變化的不同特征來去除噪聲。文中將基于小波變換的模極大值算法用于測試含噪信號heavy sine和紅山臺地震數(shù)據(jù)的去噪處理中,噪聲得到了抑制,有效信號得到保留,可以基本保持信號的不失真,從而證明算法的有效性??傮w來看,基本達到預(yù)期去噪的目的,有實際的應(yīng)用效果。
參考文獻:
[1]魏紅梅,黃世源,許飛.小波包去噪在地震信號預(yù)處理中的應(yīng)用[J].東北地震研究,2008,24(2):45-49.
[2]李文,劉霞,段玉波,等.基于小波熵與相關(guān)性相結(jié)合的小波模極大值地震信號去噪[J].地震學(xué)報,2012,34(6):841-850.
[3]徐宏斌,李庶林,陳際經(jīng).基于小波變換的大尺度巖體結(jié)構(gòu)微震監(jiān)測信號去噪方法研究[J].地震學(xué)報,2012,34(1):85-96.
[4]孫志新.小波分析在地震數(shù)據(jù)噪聲處理中的應(yīng)用研究[D].中國地震局工程力學(xué)研究所,2008:25-35.
[5]李偉,馬欽忠,宋志平,等.小波變換在地電場數(shù)據(jù)分析中的應(yīng)用[J].地震學(xué)報,2013,35(1):26-35.
[6]李杰,劉希強,李紅,等.利用小波變換方法分析形變觀測資料的正常背景變化特征[J].地震學(xué)報,2005,27(1):33-41.
[7]彭玉華.小波變換與工程應(yīng)用[M].北京:科學(xué)出版社,1999.
[8]張德豐.MATLAB小波分析第2版[M].北京:機械工業(yè)出版社,2011.
[9]田金文,高謙,杜擁軍.基于小波變換模極大值的地震信號去噪處理方法[J].江漢石油學(xué)院學(xué)報,2001,23(1):22-25.
[10]姚勝利.地震信號的小波去噪方法研究[D].長沙:中南大學(xué),2007:22-37.
[11]孔祥茜,吳繼偉,岳繼光.地震信號小波變換的去噪方法[J].計算機輔助工程,2005,14(3):52-56.
[12]楊建國.小波分析及其工程應(yīng)用[M].北京:機械工業(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
作者簡介:第一羅娜(1981—),女,河北省保定人。2009年畢業(yè)于蘭州理工大學(xué),碩士研究生,工程師。
基金項目:河北省地震局碩(博)預(yù)研究項目(201201)。
收稿日期:2015-01-05
中圖分類號:P315.3
文獻標志碼:A
文章編號:1000-6265(2015)02-0012-04