梁 躍董桂菊崔天時(shí)張東海

1）中國哈爾濱 150030 東北農(nóng)業(yè)大學(xué)

2）中國黑龍江 157009 黑龍江省牡丹江地震臺(tái)

3）中國黑龍江 158300 黑龍江省密山地震臺(tái)

密山地電場(chǎng)數(shù)據(jù)小波包去噪研究

梁 躍1），2）董桂菊1）崔天時(shí)1）張東海3）

1）中國哈爾濱 150030 東北農(nóng)業(yè)大學(xué)

2）中國黑龍江 157009 黑龍江省牡丹江地震臺(tái)

3）中國黑龍江 158300 黑龍江省密山地震臺(tái)

選取2009年11月21日密山市知一鎮(zhèn)地震發(fā)生當(dāng)天密山地電場(chǎng)數(shù)據(jù)，采用小波包閾值去噪算法，對(duì)地電場(chǎng)數(shù)據(jù)進(jìn)行去噪研究，還原地震發(fā)生后地電場(chǎng)數(shù)據(jù)變化信息，并結(jié)合去噪效果及誤差分析結(jié)果，選出適合密山地電場(chǎng)數(shù)據(jù)去噪的最優(yōu)小波函數(shù)。

地電場(chǎng)；小波包；小波；閾值去噪；小波函數(shù)

0 引言

小波變換自出現(xiàn)以來因其優(yōu)越的時(shí)頻分析特性，已經(jīng)廣泛應(yīng)用在信號(hào)分析、圖像處理、地震勘探數(shù)據(jù)處理等諸多領(lǐng)域。小波包變換是由小波變換進(jìn)一步發(fā)展而來的，同小波變換相比，小波包變換具有更高的時(shí)頻分辨率及更加優(yōu)越的去噪性能，尤其對(duì)于非平穩(wěn)信號(hào)，小波包去噪已經(jīng)成為信號(hào)處理的主流方向之一。魏紅梅等（2008）將小波包去噪方法應(yīng)用于地震信號(hào)的預(yù)處理中，削弱了地震信號(hào)中的高、低頻干擾信號(hào)。邱穎等（2009）應(yīng)用小波域閾值濾波方法對(duì)疊加在大地電場(chǎng)信號(hào)上的地電阻率人工供電干擾信號(hào)進(jìn)行濾波處理。王培茂等（2011）提出基于匹配算法的小波包去噪方法，很好地壓制了地震信號(hào)中的白噪聲。

地電觀測(cè)資料信息量大，各種高頻信息豐富、干擾因素復(fù)雜（張彩艷等，2011）。密山地電場(chǎng)觀測(cè)數(shù)據(jù)一直存在類似尖脈沖干擾信號(hào)，可能是由多種不同因素疊加到一起造成的（高研等，2012）。2009年11月21日密山市知一鎮(zhèn)發(fā)生ML3.7地震，密山地震臺(tái)測(cè)定震中距為27.8 km，地震發(fā)生前后密山地電場(chǎng)數(shù)據(jù)均未體現(xiàn)任何變化信息。本文以此次地震為例，采用小波包閾值去噪方法，還原地震發(fā)生后地電場(chǎng)數(shù)據(jù)中的變化信息，并選出適合密山地電場(chǎng)數(shù)據(jù)去噪的小波函數(shù)。

1 原理及方法

自1822年傅里葉變換發(fā)表以來，傅里葉變換一直是信號(hào)處理領(lǐng)域中應(yīng)用廣泛的分析手段。但是，由于傅里葉變換是一種全局變換，無法表述信號(hào)的時(shí)頻局部特性。盡管窗口傅里葉變換的提出實(shí)現(xiàn)了信號(hào)在時(shí)頻域內(nèi)局部化的聯(lián)合分析，但因其時(shí)頻窗的寬度是固定不變的，無法根據(jù)信號(hào)的頻率變化進(jìn)行調(diào)整。小波變換繼承和發(fā)展了窗口傅里葉變換的局部化思想，克服了窗口大小不隨頻率變化的缺點(diǎn)，在低頻部分具有較低的時(shí)間分辨率和較高的頻率分辨率，在高頻部分具有較高的時(shí)間分辨率和較低的頻率分辨率，是一種比較理想的信號(hào)處理工具。小波包變換是基于小波變換的進(jìn)一步發(fā)展，其優(yōu)一點(diǎn)在于，小波包分解對(duì)小波分解中沒有細(xì)分的高頻部分做進(jìn)一步分解，具有更高的時(shí)頻分辨率，圖1和圖2分別為小波及小波包的3層分解結(jié)構(gòu)圖。

圖1 小波3層分解結(jié)構(gòu)Fig.1 Three layer structure chart of wavelet decomposition

圖2 小波包3層分解結(jié)構(gòu)Fig.2 Three layer structure chart of wavelet packet decomposition

1.1 小波包變換定義

在多分辨率分析中，定義尺度函數(shù)空間Vj和小波函數(shù)空間Wj，尺度函數(shù)φ（t）和小波函數(shù)ψ（t）存在以下二尺度方程

將式（1）和式（2）推廣為

小波包變換的重構(gòu)算法為

1.2 閾值去噪法

閾值去噪法認(rèn)為，與有用信號(hào)對(duì)應(yīng)的小波包系數(shù)包含重要信息，所以其幅值大、數(shù)目少；而與噪聲對(duì)應(yīng)的小波包系數(shù)幅值小、數(shù)目多。通過選取適當(dāng)?shù)拈撝担瑢?duì)各尺度上的小波包系數(shù)進(jìn)行閾值處理，去掉與噪聲對(duì)應(yīng)的系數(shù)信息，從而達(dá)到信號(hào)去噪的目的。

小波包閾值去噪的步驟如下：①選擇小波函數(shù)并確定分解層數(shù)N，然后對(duì)數(shù)據(jù)進(jìn)行N層小波包分解；②提取閾值，并對(duì)小波包系數(shù)進(jìn)行閾值處理；③小波包重構(gòu)，重構(gòu)后得到去噪信號(hào)。在小波包去噪的過程中，核心的環(huán)節(jié)是在系數(shù)上作用閾值的過程，閾值選取的得當(dāng)與否直接影響去噪的質(zhì)量。

2 數(shù)據(jù)選取及去噪?yún)?shù)確定

2.1 地電場(chǎng)數(shù)據(jù)選取

北京時(shí)間2009年11月21日06時(shí)47分密山市知一鎮(zhèn)（131.81°E，45.53°N）發(fā)生ML3.7地震，密山地震臺(tái)測(cè)定震中距為27.8 km。選取密山地電場(chǎng)2009年11月21日00∶00至09∶59（北京時(shí)）東西（EW）向長極距和NS向長極距分鐘值數(shù)據(jù)進(jìn)行小波包去噪處理。

2.2 小波函數(shù)選擇及分解層數(shù)確定

小波函數(shù)具有緊支性、對(duì)稱性、正則性和消失矩階數(shù)等性質(zhì)。緊支性保證小波函數(shù)具有優(yōu)良的空間局部特性，函數(shù)的支撐寬度越小，局部化能力越強(qiáng)；對(duì)稱性又稱為線性相位特性，保證函數(shù)的濾波特性具有線性相移，不會(huì)造成濾波后信號(hào)的失真；正則性用來刻畫函數(shù)的光滑程度；消失矩階數(shù)反映了變換后能量的集中程度，消失矩階數(shù)越大正則性越好。幾種常用小波函數(shù)參數(shù)特性對(duì)比見表1。

表1 小波函數(shù)參數(shù)特性比較Table 1 Parameter characteristics of different wavelet functions

由于密山地電場(chǎng)數(shù)據(jù)中高頻干擾較多，所以選擇的小波函數(shù)必須具備良好的緊支性，以保證良好的局部性能；為了有利于分解后信號(hào)的精確重構(gòu)，小波函數(shù)應(yīng)具有較好的正交性；應(yīng)具備較好的對(duì)稱性及較高的消失矩階數(shù)。對(duì)比表1 中各小波函數(shù)的參數(shù)特性，Daubechies（dbN）函數(shù)符合本次研究對(duì)小波函數(shù)特性的要求。

因?yàn)樾〔ê瘮?shù)的各種特性不可能同時(shí)達(dá)到最優(yōu)，例如，函數(shù)的消失矩階數(shù)越大，正則性越好，函數(shù)越光滑；但函數(shù)的支撐長度隨之變大，緊支性降低。在dbN系列函數(shù)中，“N”代表消失矩階數(shù)，其支撐長度為“2N-1”?？紤]到消失矩階數(shù)與支撐長度之間的平衡，N不能太大也不能太小。因此，選擇Daubechies（dbN）系列中的 “db3”、“db4”、“db5”、“db6”函數(shù)進(jìn)行密山地電場(chǎng)數(shù)據(jù)的去噪處理，并根據(jù)去噪結(jié)果選出最佳小波函數(shù)。

正常地電場(chǎng)信號(hào)的能量分布在0—0.005 Hz，所以對(duì)含有干擾的地電場(chǎng)信號(hào)的分解層數(shù)可選擇6層或7層分解（邱影，2008），本次研究選擇7層小波包分解。

3 數(shù)據(jù)處理分析

首先采用ddencmp閾值獲取函數(shù)獲取EW向地電場(chǎng)數(shù)據(jù)在小波包去噪過程中的默認(rèn)閾值，然后使用“db3”、“db4”、“db5”及“db6”函數(shù)對(duì)數(shù)據(jù)進(jìn)行小波包閾值去噪處理，去噪方式為硬閾值去噪，熵標(biāo)準(zhǔn)默認(rèn)為sure熵標(biāo)準(zhǔn)，去噪結(jié)果見圖3，圖中箭頭所指時(shí)間為地震的發(fā)震時(shí)刻（下同）。

圖3 EW向數(shù)據(jù)小波包閾值去噪（a）EW向數(shù)據(jù)；（b）db3去噪；（c）db4去噪；（d）db5去噪；（e）db6去噪Fig.3 wavelet packet threshold de-noising on EW

與地震有關(guān)的自然電場(chǎng)異常是因?yàn)樵谠姓痼w及其周圍，受統(tǒng)一應(yīng)力場(chǎng)變化的影響，地下流體的沖流、過濾、定向作用使自然電場(chǎng)出現(xiàn)異常（馬欽忠等，2004），而自然電場(chǎng)的異常必然造成該區(qū)域內(nèi)大地電場(chǎng)發(fā)生變化。如圖3（a）所示，由于地電場(chǎng)數(shù)據(jù)干擾較多，EW向數(shù)據(jù)中未顯示任何震后變化信息。將圖3中各個(gè)函數(shù)的去噪結(jié)果進(jìn)行比較，隨著函數(shù)消失矩階數(shù)的增大，去噪后的信號(hào)越來越光滑，但函數(shù)的局部性能逐漸下降，去噪后細(xì)節(jié)信號(hào)丟失越嚴(yán)重。db3函數(shù)消失矩階數(shù)最小，去噪后信號(hào)不光滑，存在失真現(xiàn)象；db6函數(shù)消失矩階數(shù)最大，去噪后信號(hào)最光滑，但因局部性能較差，去噪后原始信號(hào)中的細(xì)節(jié)信息丟失較嚴(yán)重，未能還原地震發(fā)生后的任何變化信息。如圖3（c）所示，db4函數(shù)在去除大部分高頻干擾的同時(shí)，還原地震發(fā)生后地電場(chǎng)信號(hào)中的變化信息。由于地下流體的沖流、過濾等物理化學(xué)變化需要一定時(shí)間，所以地電場(chǎng)的同震變化信息略滯后且持續(xù)一段時(shí)間，圖3（c）中地震發(fā)生后地電場(chǎng)的波動(dòng)變化持續(xù)近半小時(shí)。

計(jì)算去噪后信號(hào)與原始信號(hào)之間的均方根誤差（RMSE）及相對(duì)誤差（error）值，見表2。

均方根誤差的定義為，原始信號(hào)與小波包去噪后重構(gòu)信號(hào)之差平方期望值的平方根；相對(duì)誤差的定義為原始信號(hào)與去噪后信號(hào)之間的差異。由表2可見，EW向數(shù)據(jù)，db4函數(shù)小波包去噪信號(hào)的兩項(xiàng)誤差值均小于其他函數(shù)，且圖3中db4函數(shù)的去噪效果最優(yōu)。所以，密山地電場(chǎng)EW向數(shù)據(jù)的最佳小波函數(shù)為db4函數(shù)。

同理，得出密山地電場(chǎng)NS向數(shù)據(jù)的最佳小波函數(shù)為db5函數(shù)，NS向數(shù)據(jù)小波包去噪信號(hào)的均方根誤差（RMSE）和相對(duì)誤差（error）值見表2。

表2 EW向數(shù)據(jù)誤差分析和NS向數(shù)據(jù)誤差分析Table 2 Error analysis of EW data

如表2所示，NS向數(shù)據(jù)db5函數(shù)小波包去噪信號(hào)的均方根誤差（RMSE）和相對(duì)誤差（error）值略大于db3函數(shù)，但小于db4和db6函數(shù)；此外，因db5函數(shù)的正則性和消失矩階數(shù)均優(yōu)于db3函數(shù)，所以，密山地電場(chǎng)NS數(shù)據(jù)的最佳小波函數(shù)為db5函數(shù)。

4 結(jié)論

通過對(duì)2009年11月21日密山地電場(chǎng)數(shù)據(jù)的小波包去噪結(jié)果進(jìn)行分析，得出以下結(jié)論。

（1）密山地電場(chǎng)數(shù)據(jù)確實(shí)存在大量干擾噪聲，通過小波包去噪方法去除大部分高頻干擾，還原地震發(fā)生后地電場(chǎng)數(shù)據(jù)的變化信息。

（2）綜合分析，密山地電場(chǎng)EW向數(shù)據(jù)去噪的最優(yōu)小波函數(shù)為db4函數(shù)，NS向數(shù)據(jù)去噪的最優(yōu)小波函數(shù)為db5函數(shù)。

（3）除地震發(fā)生時(shí)刻及其后一段時(shí)間以外，其他時(shí)刻的地電場(chǎng)數(shù)據(jù)仍存在類似突跳的干擾信息，需今后進(jìn)一步研究。

邊威.小波基的選取與構(gòu)造方法討論［D］.長春：東北師范大學(xué)，2007.

陳強(qiáng)，黃聲享，王韋.小波去噪效果評(píng)價(jià)的另一指標(biāo)［J］.測(cè)繪信息與工程，2008，33（5）：13-14.

高研，張亞江，李明忠，等.密山地震臺(tái)地電場(chǎng)干擾分析研究［J］.地震地磁觀測(cè)與研究，2012，33（1）：70-74.

馬欽忠，馮志生，宋治平，等.崇明與南京臺(tái)震前地電場(chǎng)變化異常分析［J］.地震學(xué)報(bào)，2004，26（3）：304-312.

邱影.地電場(chǎng)觀測(cè)中已知源干擾抑制研究［D］.北京：中國地震局地震預(yù)測(cè)研究所，2008.

邱穎，席繼樓.小波方法在地電場(chǎng)干擾處理中的分析研究［J］.地震，2009，29（2）：57-63.

陶珂，朱建軍.小波去噪質(zhì)量評(píng)價(jià)方法的對(duì)比研究［J］.大地測(cè)量與地球動(dòng)力學(xué)，2012，33（2）：128-133.

唐曉初.小波分析及其應(yīng)用［M］.重慶大學(xué)出版社，2006：96-98.

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

王培茂，楊冬紅，高茹，等.基于匹配小波包算法的地震信號(hào)去噪［J］.世界地質(zhì)，2011，30（2）：277-281.

席繼樓，趙家騮，王燕瓊.大地電場(chǎng)測(cè)量系統(tǒng)中的特殊抗干擾技術(shù)［J］.電子技術(shù)應(yīng)用，1999，125（12）：30-32.

張彩艷，趙潔，雷功明，等.甘肅嘉峪關(guān)地震臺(tái)大地電場(chǎng)觀測(cè)數(shù)據(jù)干擾分析［J］.山西地震，2011，4：25-28.

張德豐.Matlab小波分析［M］.北京：機(jī)械工業(yè)出版社，2010：158-196.

Research of wavelet packet de-noising on Mishan geoelectric fi eld data

Liang Yue1），2），Dong Guiju1），Cui Tianshi1）and Zhang Donghai3）
1） Northeast Agricultural University，Harbin 150030，China
2） Mudanjiang Seismic Station，Heilongjiang Province 157009，China
3） Mishan Seismic Station，Heilongjiang Province 158300，China

This paper，based on the data of the earthquake that took place in the Zhiyi town of Mishan on 21 November，2009 and was collected by the Mishan Seismic Station geoelectric fi eld（hereinafter referred to as Mishan geoelectric field），makes de-noising research of geoelectric field data by using wavelet packet threshold de-noising algorithm，and restores the change information in geoelectric fi eld data after the earthquake.In addition，it screens out the optimal wavelet functions which are appropriate for the Mishan geoelectric fi eld data de-noising according to the de-noising effect and the error analysis results.

geoelectric fi eld，wavelet packet，wavelet，threshold de-noising

10.3969/j.issn.1003-3246.2015.01.014

梁躍（1988—），男，黑龍江五常人，碩士，主要從事地震監(jiān)測(cè)工作，E-mail: 381665839@qq.com

中國地震局2014年度“地震監(jiān)測(cè)、預(yù)測(cè)、科研三結(jié)合”課題（140801）資助

本文收到日期：2014-09-29