劉淑聰,高爾根,陳遜,劉春俠
(防災(zāi)科技學(xué)院,河北三河065201)
小波包多閾值法在地震信號(hào)去噪中的應(yīng)用研究
劉淑聰,高爾根,陳遜,劉春俠
(防災(zāi)科技學(xué)院,河北三河065201)
地震信號(hào)中通常含有各種干擾噪聲,嚴(yán)重影響了地震資料的信噪比和分辨率,小波包變換是地震資料去噪的有效方法之一。針對(duì)傳統(tǒng)小波包閾值去噪不明顯和存在失真的問(wèn)題,提出一種基于多閾值函數(shù)的小波包地震信號(hào)去噪方法。對(duì)地震波信號(hào)進(jìn)行小波包分解,并對(duì)小波包分解系數(shù)按照頻率大小的順序進(jìn)行排列,根據(jù)分解的系數(shù)處于不同頻帶選取不同的閾值準(zhǔn)則進(jìn)行去噪處理,對(duì)得到的系數(shù)進(jìn)行重構(gòu),可有效地去除地震信號(hào)中的噪聲。對(duì)仿真地震信號(hào)以及實(shí)際地震信號(hào)進(jìn)行小波包多閾值去噪處理,實(shí)驗(yàn)結(jié)果表明,該方法較好地去除了干擾噪聲保留了有用信號(hào),去噪效果明顯且失真小,有效地提高了地震資料的分辨率。
地震信號(hào);小波包變換;多閾值函數(shù);信號(hào)去噪
在地震勘探中,由于地震信號(hào)的復(fù)雜性,不可避免地要受到各種干擾噪聲的影響,使得采集到的原始地震信息剖面模糊不清,地震資料信噪比和分辨率降低,地震資料解釋與處理更為困難,因此降噪是地震資料處理的關(guān)鍵。很多學(xué)者根據(jù)信號(hào)和噪聲的各種差異,設(shè)計(jì)了多種去除噪聲、提高信噪比的方法[1-2]。小波變換作為一種新方法技術(shù),由于其具有時(shí)頻分析、多分辨率和去相關(guān)性等特點(diǎn),在信號(hào)處理方面得到了廣泛應(yīng)用[3-4]。小波包變換是小波變換的推廣,不僅對(duì)信號(hào)的低頻部分進(jìn)行分解,而且可以對(duì)小波分析未分解的高頻部分提供更精細(xì)的分解,能夠?qū)Π小⒏哳l信息的信號(hào)進(jìn)行更好地時(shí)頻局部化分析。小波包變換是一種更加精細(xì)的分析和重構(gòu)方法,對(duì)地震信號(hào)中的高頻噪聲有一定的去噪效果[5-6],但傳統(tǒng)的小波包去噪方法沒(méi)有對(duì)信號(hào)和噪聲的特征進(jìn)行充分的研究,去噪效果不明顯,去噪精度有待進(jìn)一步提高。
本文提出一種對(duì)不同頻率成分選擇不同閾值函數(shù)的小波包地震信號(hào)去噪方法,對(duì)小波包分解系數(shù)按照頻率大小的順序進(jìn)行排列,根據(jù)信號(hào)和噪聲的分解系數(shù)處于不同頻帶選取不同的閾值準(zhǔn)則進(jìn)行去噪處理,對(duì)得到的系數(shù)進(jìn)行重構(gòu),可有效地去除信號(hào)中的噪聲。通過(guò)本文小波包多閾值去噪方法對(duì)人工合成地震記錄以及實(shí)際地震信號(hào)進(jìn)行降噪處理,結(jié)果表明小波包多閾值地震信號(hào)去噪方法能夠有效提取信號(hào)中各頻段的有用信息,消除干擾信號(hào),去噪能力相對(duì)于傳統(tǒng)方法有所提高,可獲得較理想的去噪效果。
1.1 小波包分解
在小波分析去噪中實(shí)際信號(hào)往往被認(rèn)為由有用信號(hào)和噪聲組成,經(jīng)過(guò)小波變換后,有用信號(hào)的能量主要集中在低頻的小波系數(shù)上,噪聲信號(hào)的能量主要集中在高頻的小波系數(shù)上,因此通常把高頻部分去掉(置零),低頻部分用于進(jìn)一步分解和重構(gòu),這樣就損失了包含在高頻段的有用信號(hào)。小波包分解是一種比小波分析更加精細(xì)的時(shí)頻分析方法,它可以對(duì)小波空間進(jìn)行進(jìn)一步分解,能同時(shí)對(duì)信號(hào)的高頻部分和低頻部分進(jìn)行多層次劃分,具有更為精確的時(shí)頻局部分析能力,更加有利于信號(hào)和噪聲的分離。
通常情況下,小波包降噪步驟為[7]:
(1)信號(hào)的小波包分解:選擇合適的小波基函數(shù)并確定小波包分解的層次N,對(duì)地震波記錄信號(hào)S進(jìn)行N層小波包分解。A表示低頻,D表示高頻,末尾序號(hào)數(shù)表示小波包分解的層數(shù),即尺度數(shù)。原始信號(hào)S經(jīng)過(guò)3層分解后等價(jià)于:
(2)計(jì)算最佳小波基:根據(jù)一個(gè)給定的熵標(biāo)準(zhǔn)計(jì)算最佳樹(shù),常用熵標(biāo)準(zhǔn)有shannon,threshold,norm,log energy,sure和user等幾種類(lèi)型。
(3)小波包分解系數(shù)的閾值量化:為有效區(qū)分信號(hào)和噪聲,需對(duì)分解的每一個(gè)小波包系數(shù)選擇合適的閾值進(jìn)行閾值處理。常用的閾值處理方法有強(qiáng)制消噪處理、默認(rèn)閾值消噪處理、軟(或硬)閾值消噪處理。軟(或硬)閾值消噪最為常用,但是仍然存在問(wèn)題。
(4)信號(hào)重構(gòu):利用經(jīng)過(guò)閾值處理后的小波包系數(shù)進(jìn)行重構(gòu)即可獲得去噪處理后的信號(hào)。
1.2 小波包的頻率順序
以N=3為例,信號(hào)S經(jīng)過(guò)3層小波包分解后,原始信號(hào)S=AAA3+DAA3+ADA3+DDA3+AAD3+DAD3+ADD3+ DDD3。
其中A表示低頻,D表示高頻,末尾序號(hào)數(shù)3表示小波包分解層數(shù)(尺度數(shù)為3),分解結(jié)構(gòu)如圖1所示。
圖1 小波包3層分解結(jié)構(gòu)圖
以noischir信號(hào)進(jìn)行sym6小波3層分解,分解后的各小波包系數(shù)如圖2所示。通過(guò)圖2可以發(fā)現(xiàn),小波包系數(shù)的自然順序與其頻率順序不一致,最低頻部分對(duì)應(yīng)的是AAA3,其次是DAA3,最高頻部分對(duì)應(yīng)的是AAD3。
圖2 noischir信號(hào)3層小波包分解系數(shù)
由文獻(xiàn)[8]可知小波包分解后的結(jié)果并不是頻率由低到高排列,任何小波包都會(huì)產(chǎn)生頻率順序錯(cuò)位,由于小波包分解時(shí),高通濾波器會(huì)進(jìn)行一次“翻轉(zhuǎn)”操作,故小波包分解后的低頻部分按頻率從小到大排列,高頻部分按頻率從大到小排列??梢酝ㄟ^(guò)小波包的性質(zhì)證明[8],任何小波包的分解都會(huì)產(chǎn)生自然順序與頻率順序不一致的現(xiàn)象,并且不一致現(xiàn)象的情況是一樣的。
2.1 最佳小波包基的確定(計(jì)算最佳樹(shù))
一個(gè)長(zhǎng)為L(zhǎng)=2N的信號(hào)最多有2N種不同的分解方法,同時(shí)一個(gè)深度為N的完全二叉子樹(shù)的個(gè)數(shù)為2N。這個(gè)數(shù)字很大,不可能對(duì)每一種情況一一列舉,而通過(guò)最小熵標(biāo)準(zhǔn)可以得到一種最優(yōu)的信號(hào)分解方法。傳統(tǒng)的基于熵的標(biāo)準(zhǔn)有shannon熵、threshold熵、norm熵、log energy熵、sure熵和user熵等幾種類(lèi)型。shannon熵定義如下:
對(duì)信號(hào)進(jìn)行逐層分解,分解的每個(gè)節(jié)點(diǎn)計(jì)算熵值,通過(guò)比較某一節(jié)點(diǎn)與其子節(jié)點(diǎn)的熵值,獲取熵值最小的基即為最優(yōu)小波包基。
2.2 多閾值選取
在小波閾值去噪中通常認(rèn)為噪聲表現(xiàn)為高頻信號(hào),對(duì)小波分解的高頻系數(shù)進(jìn)行門(mén)限閾值處理,小于閾值認(rèn)為是由噪聲產(chǎn)生的并置零,大于閾值相應(yīng)于有用信號(hào)將其保留,從而實(shí)現(xiàn)信號(hào)和噪聲分離。小波包分解可以對(duì)小波分解產(chǎn)生的高頻信號(hào)進(jìn)一步分解從而獲得高頻部分更詳細(xì)的信息,有助于將高頻噪聲和高頻信號(hào)區(qū)分開(kāi)來(lái),從而獲得更為理想的去噪效果。經(jīng)過(guò)小波包分解后的小波包系數(shù)包含著不同的信息,有用信號(hào)和噪聲信號(hào),如果對(duì)所有的小波包系數(shù)采取同一種閾值方法處理,便容易出現(xiàn)去噪效果不明顯或是去噪過(guò)度等現(xiàn)象影響去噪精度,若對(duì)小波包系數(shù)分情況進(jìn)行多閾值函數(shù)處理就很好地克服了上述缺點(diǎn)。
多閾值處理就是對(duì)每個(gè)小波包分解系數(shù)靈活地選用不同的閾值準(zhǔn)則,從而最大程度地保留有用信號(hào)去除噪聲信號(hào)。在了解了小波包分解系數(shù)的頻率順序之后,就可以將小波包分解系數(shù)按照頻率順序從小到大排列,從而進(jìn)行頻帶分割。根據(jù)噪聲與信號(hào)在各頻帶上的小波包分解系數(shù)的不同表現(xiàn),采用多個(gè)閾值進(jìn)行去噪。在不同的頻帶采取不同的閾值處理方法,最后進(jìn)行信號(hào)重構(gòu),實(shí)現(xiàn)不同頻帶多閾值去噪處理。
小波包分析中常用的四種閾值準(zhǔn)則有固定形式閾值準(zhǔn)則(sqtwolog準(zhǔn)則)、自適應(yīng)閾值準(zhǔn)則(rigrsure準(zhǔn)則)、啟發(fā)式閾值準(zhǔn)則(heursure準(zhǔn)則)以及極小化極大閾值準(zhǔn)則(minimaxi準(zhǔn)則),各自的選取規(guī)則不同,適用范圍不同,去噪效果也不相同[9]。
(1)固定形式閾值準(zhǔn)則(sqtwolog),采用的是固定形式的閾值,產(chǎn)生的閾值大小為sqrt(2?log(length(X)))。
(2)自適應(yīng)閾值準(zhǔn)則(rigrsure),基于Stein的無(wú)偏似然估計(jì)原理的自適應(yīng)閾值選擇。
(3)啟發(fā)式閾值準(zhǔn)則(heursure),是無(wú)偏似然估計(jì)和固定閾值估計(jì)原則的折中。如果按無(wú)偏似然估計(jì)原則處理的信號(hào)噪聲很大,在這種情況下就采用這種固定的閾值。
(4)極小化極大閾值準(zhǔn)則(minimaxi),采用極大極小原理選擇閾值,它產(chǎn)生一個(gè)最小均方差的極值,而不是沒(méi)有誤差。
minimaxi閾值準(zhǔn)則和rigrsure準(zhǔn)則更加保守,不容易丟失信號(hào)中的有用成份,但只除去較少的噪聲。sqtwolog閾值準(zhǔn)則以及heursure閾值準(zhǔn)則兩種方法類(lèi)似,都是將全部系數(shù)進(jìn)行處理,因此可以較強(qiáng)地去除噪聲,但是也容易引起過(guò)度去噪使信號(hào)失真。因此,本文采用小波包多閾值處理,將分解后的系數(shù)分為不同的頻帶,中低頻部分采用minimaxi閾值準(zhǔn)則和rigrsure閾值準(zhǔn)則,而在高頻部分采用sqtwolog閾值準(zhǔn)則以及heursure閾值準(zhǔn)則。閾值準(zhǔn)則選取之后可對(duì)每層的系數(shù)進(jìn)行噪聲層的估計(jì)來(lái)調(diào)整閾值。
2.3 算法步驟
(1)根據(jù)所給信號(hào)選定一個(gè)合適的小波基函數(shù),先對(duì)信號(hào)進(jìn)行N層小波包分解,獲得與樹(shù)T節(jié)點(diǎn)N對(duì)應(yīng)的系數(shù)cfi(i=0,1,2,…),并按頻帶進(jìn)行排序。
(2)優(yōu)化小波包樹(shù),確定熵標(biāo)準(zhǔn)從而選定最優(yōu)小波基,得到最優(yōu)樹(shù)。
(3)依據(jù)信號(hào)獲取sqtwolog準(zhǔn)則對(duì)應(yīng)的閾值th1,rigrsure準(zhǔn)則對(duì)應(yīng)的閾值th2,heursure準(zhǔn)則對(duì)應(yīng)的閾值th3,minimaxi閾值準(zhǔn)則對(duì)應(yīng)的閾值th4。
(4)利用不同的閾值對(duì)小波包分解的系數(shù)cfi(i=0,1,2,…)進(jìn)行閾值處理,低頻段采用閾值th2,中頻段采用閾值th4或th1,高頻段采用閾值th3??稍O(shè)定各頻段比例系數(shù)。
(5)對(duì)處理后的小波包系數(shù)進(jìn)行重構(gòu)。
3.1 模擬地震記錄去噪
為了驗(yàn)證新閾值函數(shù)小波包消噪的有效性,對(duì)一段加入隨機(jī)噪聲的模擬地震信號(hào)進(jìn)行試驗(yàn),信號(hào)的長(zhǎng)度為1 000,圖3是對(duì)含噪聲的模擬地震記錄處理對(duì)比圖。其中圖3(a)為模擬地震信號(hào),信噪比為7 dB,圖3(b)為小波包硬閾值去噪后的信號(hào),圖3(c)為小波包軟閾值去噪后的信號(hào),采用小波包分解時(shí),選取“dB2”小波基,最大分解尺度為3。圖3(d)為小波包多閾值方法去噪后的信號(hào)。
為了比較不同閾值降噪方法的降噪效果,選用均方根誤差(RMSE)、信噪比(SNR)為評(píng)價(jià)指標(biāo)用來(lái)定量比較去噪效果,如表1所示。將原始信號(hào)記為f(n),含噪信號(hào)為s(n),信號(hào)長(zhǎng)度為L(zhǎng),信噪比(SNR)公式定義為:
原始信號(hào)與估計(jì)信號(hào)之間的均方根誤差(RMSE)定義為:
其中均方根誤差越小、信噪比越大則去噪效果越好。
圖3 模擬地震信號(hào)去噪
表1 模擬記錄各種方法去噪效果比較
從圖3和表1中可以看出,以上的幾種方法去噪后信噪比都有所提高,大部分噪聲都得到了抑制,但是使用小波包多閾值法去噪后,失真最小,信噪比(SNR)和均方根誤差(RMSE)最大,去噪效果最好。
3.2 合成地震剖面去噪
在Matlab下可仿真合成地震剖面,通常采用一個(gè)低頻ricker子波和一個(gè)高頻ticker子波疊加合成共炮點(diǎn)道集[10],設(shè)置參數(shù)深度為100 m,最小偏移距為10 m,道數(shù)為50。圖4為合成的共炮點(diǎn)道集,圖5為加入20 dB隨機(jī)噪聲后的共炮點(diǎn)道集。采用sym4小波進(jìn)行小波包去噪,圖6是小波包硬閾值去噪后的共炮點(diǎn)道集,圖7是小波包軟閾值去噪后的共炮點(diǎn)道集,圖8是小波包多閾值法去噪后的共炮點(diǎn)道集。從圖6~圖8可以看出,采用小波包多閾值方法去噪后噪聲有效地減小,相比常規(guī)小波包閾值去噪方法,失真小,去噪后的地震信號(hào)最接近于原始信號(hào)。從表2中可以看出,經(jīng)過(guò)小波包硬閾值去噪、小波包軟閾值去噪和小波包多閾值去噪后的信噪比分別為36.700 6 dB,35.430 9 dB和39.974 6 dB,采用小波包多閾值方法去噪的信噪比最高,效果最好。
圖4 合成的共炮點(diǎn)道集
圖5 加入噪聲后的共炮點(diǎn)道集
圖6 小波包硬閾值去噪后的共炮點(diǎn)道集
實(shí)際采集的原始地震數(shù)據(jù)有大量的噪聲,如隨機(jī)干擾,線性干擾等噪聲,這些噪聲相互混合影響原始地震數(shù)據(jù),使整個(gè)地震剖面分辨率較低,必須對(duì)采集的數(shù)據(jù)進(jìn)行消噪處理。實(shí)驗(yàn)在某地區(qū)觀測(cè)點(diǎn)采用分布式地震儀進(jìn)行觀測(cè),26通道,道間距為1 m,偏移距為10 m,炮點(diǎn)距(兩炮點(diǎn)間距離)為2 m,檢波器陣列與震源陣列共線,地震儀采樣率為8 kHz,記錄時(shí)間為6 s。對(duì)采集到的實(shí)際地震信號(hào)選用sym4小波進(jìn)行3次分解。
圖7 小波包軟閾值去噪后的共炮點(diǎn)道集
圖8 小波包多閾值去噪后的共炮點(diǎn)道集
表2 地震剖面各種方法去噪效果比較
4.1 單通道地震信號(hào)去噪
圖9(a)為其中的一個(gè)通道實(shí)際信號(hào)中的一部分?jǐn)?shù)據(jù),圖9(b)為小波包軟閾值去噪后的信號(hào),圖9(c)為小波包硬閾值去噪后的信號(hào),圖9(d)為小波包多閾值去噪后的信號(hào)。
從圖9(b)和圖9(c)可以看出硬閾值函數(shù)去噪和軟閾值函數(shù)去噪后的地震信號(hào)都去除了一定的噪聲,但去噪效果不明顯。經(jīng)過(guò)本文小波包多閾值函數(shù)去噪后的地震信號(hào)如圖9(d)所示,濾除了大部分噪聲,有用信號(hào)分量卻很好地保留下來(lái),失真較小,去噪效果最好,經(jīng)過(guò)小波包多閾值函數(shù)去噪后的地震信號(hào)也最接近于原始信號(hào)。
圖10分別為實(shí)際信號(hào)及三種方法去噪后的信號(hào)頻譜圖,從圖中可以看出信號(hào)的主要頻率在20 Hz和60 Hz左右,經(jīng)過(guò)去噪后噪聲都有所壓制,但小波包多閾值法去噪后的有用信號(hào)得到增強(qiáng),噪聲得到了更好地壓制,可以看出其去噪效果最好。
圖9 單通道地震信號(hào)去噪處理
圖10 單通道地震信號(hào)去噪頻譜圖
4.2 多通道地震信號(hào)去噪
圖11(a)為26通道部分實(shí)際地震信號(hào),經(jīng)過(guò)小波包多閾值函數(shù)去噪后如圖11(b)所示,從圖中可以看出信號(hào)中的噪聲大部分已被消除,有用信號(hào)很好地保留,同相軸更加清晰。圖11(c)和圖11(d)分別為從圖11(a)和圖11(b)中取0~300 ms區(qū)間內(nèi)的部分波形,可發(fā)現(xiàn)26通道實(shí)際地震信號(hào)主要集中在該區(qū)間,從圖11(d)中也可以發(fā)現(xiàn)噪聲被較好地去除,有用信號(hào)增強(qiáng),同相軸更加清晰。
圖11 多通道地震信號(hào)去噪
本文提出一種基于小波包多閾值函數(shù)的地震信號(hào)去噪方法,將信號(hào)經(jīng)過(guò)小波包分解后的系數(shù)按照頻率大小的順序排列,在不同頻段范圍內(nèi)選取合適的閾值方法進(jìn)行自適應(yīng)閾值處理。理論分析和實(shí)際地震信號(hào)應(yīng)用表明該方法能夠?qū)π盘?hào)按不同頻段進(jìn)行去噪處理,既能有效地去除低頻中的噪聲,又很好地保留了高頻中的有用信號(hào),有效地解決了其他方法去噪效果不明顯或過(guò)度去噪的問(wèn)題,可提高地震資料的信噪比和分辨率,可應(yīng)用于實(shí)際地震資料的處理中。
[1]張軍華,呂寧,田連玉,等.地震資料去噪方法技術(shù)綜合評(píng)述[J].地球物理學(xué)進(jìn)展,2006,21(2):546-553.
[2]劉霞,潘洪屏,高曉春,等.基于小波閾值的地震信號(hào)去噪處理[J].科學(xué)技術(shù)與工程,2010,10(29):7251-7254.
[3]裴正林.小波理論及其在地震數(shù)據(jù)處理中的應(yīng)用[J].地球物理學(xué)進(jìn)展,2002,17(3):486-490.
[4]吳招才,劉天佑.地震數(shù)據(jù)去噪中的小波方法[J].地球物理學(xué)進(jìn)展,2008,23(2):493-498.
[5]張瑞紅,林大超,喬蘭.最優(yōu)小波包變換在地震信號(hào)去噪中的應(yīng)用[J].地震研究,2011,34(3):358-363.
[6]周懷來(lái),李錄明,李枚,等.基于小波包分析的地震資料聯(lián)合去噪方法研究[J].物探化探計(jì)算技術(shù),2007,29(5):391-395.
[7]胡昌華,李國(guó)華,周濤.基于Matlab 7.x的系統(tǒng)分析與設(shè)計(jì):小波分析[M].3版.西安:西安電子科技大學(xué)出版社,2008.
[8]紀(jì)躍波.小波包的頻率順序[J].振動(dòng)與沖擊,2005,24(3):96-98.
[9]章浙濤,朱建軍,匡翠林,等.小波包多閾值去噪法及其在形變分析中的應(yīng)用[J].測(cè)繪學(xué)報(bào),2014,43(1):14-20.
[10]劉霞,潘洪屏,高曉春,等.基于小波閾值的地震信號(hào)去噪處理[J].科學(xué)技術(shù)與工程,2010,10(29):7251-7254.
[11]趙志剛,萬(wàn)嬌娜,管聰慧.基于小波包變換與自適應(yīng)閾值的圖像去噪[J].中國(guó)圖象圖形學(xué)報(bào),2007,12(6):977-980.
[12]DUHAMEL P.Algorithms meeting the lower bounds on the multiplicative complexity of length-2nDFTs and their connection with practical algorithms[J].IEEE Transactions on Acoustics,Speech and Signal Processing,1990,38(9):1504-1511.
Application of wavelet packet multi-threshold method in seismic signal denoising
LIU Shucong,GAO Ergen,CHEN Xun,LIU Chunxia
(Institute of Disaster Prevention,Sanhe 065201,China)
Since the seismic signal generally contains a variety of interference noises,which affects the signal-to-noise ratio and resolution of seismic data seriously,and because the wavelet packet transform is an effective way of the seismic signal denoising,a seismic signal denoising method based on wavelet packet of multi-threshold function is proposed for the traditional wavelet packet threshold denoising method has unobvious denoising effect and distortion.The wavelet packet decomposition is conducted for seismic wave signal,and the wavelet packet decomposition coefficients are arranged according to the frequency real values.Different threshold criterions are selected to conduct denoising processing according to the decomposed coefficients in different bands to reconstruct the obtained coefficients,which can eliminate the noise in seismic signal effectively.The denoising processing of wavelet packet multi-threshold for the simulation seismic signal and practical seismic signal were conducted.The experimental results show this method can eliminate the interference noise availably while reserving the useful signal,has good denoising effect and low distortion,and can improve the resolution of seismic data efficiently.
seismic signal;wavelet packet transform;multi-threshold function;signal denoising
TN911.7-34;TH763
A
1004-373X(2015)23-0054-06
10.16652/j.issn.1004-373x.2015.23.016
劉淑聰(1983—),女,河北燕郊人,講師,研究生。主要研究方向?yàn)樾盘?hào)處理。
2015-07-20
中國(guó)地震局教師科研基金項(xiàng)目(20140105);國(guó)家自然科學(xué)基金2011年度面上項(xiàng)目(41174043);中央高校基本科研業(yè)務(wù)經(jīng)費(fèi)項(xiàng)目(ZY20110101)