涇陽臺(tái)數(shù)字地磁信號(hào)的干擾分析及去噪處理

魯 權(quán)， 邢西淳， 李西京， 毛 娟

(1.戶縣地震辦公室，陜西 戶縣 710300；2.陜西省地震局涇陽地震臺(tái)，陜西 咸陽 713704；3.陜西省地震局乾陵地震臺(tái)，陜西 咸陽 713300）

0 引言

隨著經(jīng)濟(jì)的發(fā)展，地震臺(tái)附近開山、爆破、破石等經(jīng)濟(jì)活動(dòng)以及大型電器設(shè)備突然啟動(dòng)、關(guān)閉或負(fù)荷突然改變和重型機(jī)動(dòng)車輛來往頻繁等[1、2]，致使涇陽臺(tái)GM4、FHD-1數(shù)字地磁信號(hào)的噪聲干擾日益嚴(yán)重。這些電磁噪聲有的離觀測(cè)點(diǎn)很近，在某些情況下噪聲干擾甚至能大過真實(shí)的地磁信號(hào)本身，觀測(cè)到的地磁信號(hào)形成尖脈沖、毛刺、臺(tái)階干擾。另外地磁觀測(cè)儀器產(chǎn)生的儀器噪聲也會(huì)對(duì)地磁數(shù)據(jù)產(chǎn)生干擾[3、4]。涇陽GM4磁通門磁力儀和FHD-1核旋質(zhì)子磁力儀觀測(cè)數(shù)據(jù)，很明顯的可以看到兩類干擾，一類是儀器供電設(shè)備對(duì)儀器觀測(cè)數(shù)據(jù)的干擾，另一類是上述各種環(huán)境干擾對(duì)觀測(cè)數(shù)據(jù)的影響。

傳播介質(zhì)的影響及接收器本身和周圍環(huán)境的干擾[5、6]，給信號(hào)檢測(cè)與識(shí)別帶來一些困難。如何識(shí)別與消除干擾是應(yīng)用數(shù)字地磁前兆資料進(jìn)行地震預(yù)測(cè)的基礎(chǔ)問題，而如何分離數(shù)字地磁資料中的高頻與低頻信息以及如何識(shí)別地震的趨勢(shì)與短期異常，也是地磁數(shù)字化前兆資料分析應(yīng)用必須解決的問題。因此利用信號(hào)處理技術(shù)，抑制和剔除干擾，提高信噪比減小各種干擾對(duì)地磁數(shù)據(jù)質(zhì)量的影響，最大限度保留原始信號(hào)，是十分必要的。

小波變換是一種信號(hào)的時(shí)頻分析，它具有多分辨率的特點(diǎn)，可以方便地從混有強(qiáng)噪聲的信號(hào)中提取原始信號(hào)，被譽(yù)為分析信號(hào)的顯微鏡[5～7]。噪聲通常包含在高頻中，因此對(duì)實(shí)際信號(hào)進(jìn)行小波分解，并對(duì)小波分解的高頻系數(shù)進(jìn)行門限閾值量化處理后進(jìn)行小波重構(gòu)，達(dá)到消除噪聲的目的，即抑制信號(hào)的噪聲，在實(shí)際信號(hào)中恢復(fù)真實(shí)信號(hào)。

本文利用小波分析方法，對(duì)涇陽臺(tái)數(shù)字化地磁觀測(cè)資料消除高頻干擾，提取真實(shí)地磁變化進(jìn)行了初步分析探討。

1 小波變換及資料處理方法

設(shè) ψ（t）為一在 （-∞,+∞）平方可積函數(shù)， 即 ψ（t）∈L2（R）， 并滿足如下條件：

函數(shù) f（t）∈L2（R）的連續(xù)小波變換 DWTa，b定義為

在連續(xù)小波變換式（2）中,DWTa，b中尺度因子 (伸縮因子)a和平移因子b(a，b∈R，a≠0)的連續(xù)變化的值[9，10]。實(shí)際應(yīng)用中，信號(hào)f（t）是離散序列， a與b也需要離散化，成為離散小波變換 DWTa，b。 即

一般采用a=2k。隨著k的增加，信號(hào)從最高頻向低頻分解。當(dāng)時(shí)k=0時(shí)，信號(hào)為采樣頻率，k=1時(shí)將頻率二等分，依此類推[3]。

當(dāng)信號(hào)的采樣率滿足Nyquist frequency（奈奎斯特頻率）要求時(shí)，用理想低通和高通濾波器將信號(hào)分解成頻帶在0～π/2的低頻部分與頻帶在π/2～π的高頻部分。分別得到信號(hào)的概貌與細(xì)節(jié)。對(duì)低頻部分重復(fù)進(jìn)行下去，直至將原始信號(hào)作多分辨分解。

一個(gè)含噪聲的一維時(shí)間序列信號(hào)模型的數(shù)學(xué)表達(dá)式[3]：

其中， f（i）為真實(shí)信號(hào)，e（i）為噪聲， s（i）為含噪聲的信號(hào)。

一般情況下，e（i）近似看作為高斯噪聲，通常表現(xiàn)為高頻性質(zhì)，其噪聲信號(hào)包含在較高頻率的細(xì)節(jié)之中。所以消噪過程可先對(duì)信號(hào)作小波分解，以門限閾值等形式對(duì)小波系數(shù)進(jìn)行處理，然后對(duì)信號(hào)進(jìn)行重構(gòu)即可達(dá)到消噪的目的。在 s（i）中恢復(fù)出現(xiàn)真實(shí)信號(hào)f（i）。

采用小波自適應(yīng)閾值方法確定噪聲閾值，通過在不同尺度上選取合適的閾值，并將小于該閾值的小波系數(shù)置零或弱化，有效去除隨小波分解尺度增大而減小的小波模極大值。而保留大于閾值的小波系數(shù)，保留隨尺度增大而增大的模極大值，然后用剩余的小波變換模極大值去重構(gòu)信號(hào)，實(shí)現(xiàn)噪聲信號(hào)與地磁信號(hào)的有效分離。

2 電源對(duì)地磁觀測(cè)數(shù)據(jù)干擾的抑制

在觀測(cè)實(shí)踐中電源連續(xù)穩(wěn)定的工作是保證正常觀測(cè)的基礎(chǔ)，數(shù)字地磁儀器電源工作不正常，可能導(dǎo)致觀測(cè)數(shù)據(jù)的斷記和數(shù)據(jù)質(zhì)量的下降，從而影響觀測(cè)數(shù)據(jù)的使用。涇陽臺(tái)FHD-1核旋磁力儀供電采用交直流雙電源模式，市電網(wǎng)停電后直流電瓶自動(dòng)給儀器供電[4]。圖1為2011年3月21日涇陽臺(tái)FHD-1記錄曲線，可以看出11時(shí)30分左右開始出現(xiàn)明顯的高頻干擾，觀測(cè)室周圍環(huán)境沒有干擾現(xiàn)象，儀器接地良好，探頭聯(lián)接正常。15時(shí)30分左右改用直流電瓶給儀器供電，記錄曲線高頻毛刺消失，因此認(rèn)為高頻干擾與直流穩(wěn)壓電源有關(guān)，可能由于FHD-1核旋磁力儀原配的直流穩(wěn)壓電源由于濾波性能下降，對(duì)觀測(cè)記錄產(chǎn)生了干擾。其后用東芝筆記本電源適配器 (15 V，電流2A)更換了原直流穩(wěn)壓電源給儀器供電，干擾現(xiàn)象排除。

圖1 FHD-1電源干擾（各分量顯示存在干擾）Fig.1 FHD-1 power supply interference

3 小波消噪在涇陽臺(tái)地磁觀測(cè)數(shù)據(jù)應(yīng)用實(shí)例

涇陽臺(tái)地磁記錄室東側(cè)約200 m、西側(cè)約400 m，有兩條鄉(xiāng)間道路，受重型機(jī)動(dòng)車運(yùn)輸所產(chǎn)生的脈沖狀干擾的影響，產(chǎn)生脈沖狀干擾。一般的脈沖狀干擾周期在4～20 s左右。脈沖干擾的幅度與鐵磁性物體的磁化率、體積、運(yùn)行速度以及與觀測(cè)儀器探頭之間距離等參數(shù)有關(guān)，涇陽臺(tái)記錄到的干擾幅度從幾個(gè)nT到幾十nT都有。圖2（a）為涇陽臺(tái)2011年10月8日地磁觀測(cè)的原始曲線，圖2（b）為采用小波降噪后的結(jié)果，圖2（c）為原始資料和降噪后結(jié)果的殘差。資料作小波變換消噪處理。采用此法的理由是該小波有較好的信號(hào)自適應(yīng)性，結(jié)果顯示能有效的去除部分干擾，信號(hào)經(jīng)重構(gòu)后，得到相對(duì)理想的曲線,與其它正常時(shí)段的觀測(cè)結(jié)果保持相似的形態(tài)，地磁真實(shí)信號(hào)保留的很完整。中誤差是估算這一波動(dòng)幅度的可行指標(biāo)之一[11～13]，計(jì)算誤差結(jié)果顯示，降噪后結(jié)果的中誤差為0.360 5，明顯優(yōu)于原始觀測(cè)資料的中誤差為6.366 1。同樣圖3給出了2011年4月30日汽車干擾的處理結(jié)果，觀測(cè)數(shù)據(jù)誤差由0.652 6降為降噪后的0.304 6，說明在抑制效果和誤差上小波方法是令人滿意的。

為進(jìn)一步分析結(jié)果的有效性，圖4給出了FHD-1 Z分量2012年7月1日至6日小波抑制噪聲對(duì)比圖，可以看出仍顯示了比較好的處理效果。圖5給出了距離涇陽地震臺(tái)約50 km的乾陵地震臺(tái)觀測(cè)資料的對(duì)比結(jié)果，通過小波抑制方法抑制后的數(shù)據(jù)既不失高頻細(xì)節(jié)信號(hào)，在地磁日變形態(tài)等長周期信號(hào)上也與乾陵臺(tái)站保持一致。

在處理過程中根據(jù)地磁信號(hào)每個(gè)細(xì)節(jié)尺度上的波動(dòng)的大小進(jìn)行估計(jì)是選取該尺度上的閾值的關(guān)鍵。由于小波基函數(shù)的選擇單一，不能涵蓋所有干擾特征形態(tài)，小波系數(shù)的物理含義也較為模糊[14]，在選擇要抑制的小波系數(shù)方面存在一定經(jīng)驗(yàn)性。

圖2 FHD-1 2011-10-08原始圖、小波消噪圖及其殘差值圖Fig.2 Original graph,wavelet denoising and their residuals value of FHD-1 on October 8th,2011

圖3 涇陽臺(tái)GM4 2011年4月30日世界時(shí)秒值Z分量原始及抑制汽車干擾對(duì)比圖Fig.3 Comparison of the original and suppression automotive interference of Z component of the seconds value(World Time)at Jingyang station GM4 on April 30th,2011

圖4 涇陽臺(tái)FHD-1 Z分量2012-07-01～06小波抑制噪聲對(duì)比Fig.4 Comparison of noise suppression of Z component at Jingyang station from July 1th to 6th,2012

圖5 乾陵臺(tái)FHD-1 Z分量原始數(shù)據(jù)圖與涇陽臺(tái)FHD-1 Z分量小波消噪數(shù)據(jù)對(duì)比Fig.5 Comparison of Z component of the original data at Qianling station FHD-1 and wavelet denoising at Jingyang station FHD-1

4 結(jié)論與討論

本文針對(duì)涇陽臺(tái)地磁觀測(cè)的情況對(duì)儀器供電設(shè)備導(dǎo)致的干擾可采用高性能的直流穩(wěn)壓電源供電 （或采用直流電瓶供電）消除干擾，對(duì)環(huán)境噪聲干擾，采用小波變換模極大值方法對(duì)受干擾的地磁信號(hào)進(jìn)行處理，結(jié)果表明：含有噪聲干擾的地磁觀測(cè)數(shù)據(jù)經(jīng)過小波變換將低頻趨勢(shì)信號(hào)與高頻噪聲及突變信號(hào)很好地分離，Z分量的細(xì)節(jié)保留充分。因此小波變換對(duì)非平穩(wěn)信號(hào)的消噪效果顯著,可以用小波變換模極大值抑制方法對(duì)地磁觀測(cè)數(shù)據(jù)進(jìn)行預(yù)處理，預(yù)處理后的結(jié)果提高了資料的質(zhì)量，可更好的研究地磁信號(hào)對(duì)震磁關(guān)系進(jìn)行探討。同時(shí)對(duì)不同的干擾特征，通過具體調(diào)試，合理的選取小波去噪閾值是有效的抑制信號(hào)噪聲的關(guān)鍵。

[1]Jankowski J,C Suksdorff.Guide For Magnetic Measurements and Observatory Practice[M].IAGA.1996.

[2]全國地震標(biāo)準(zhǔn)化技術(shù)委員會(huì).GB/T 19531.2-2004地震臺(tái)站觀測(cè)環(huán)境技術(shù)要求第二部分：電磁觀測(cè)[S].北京：中國標(biāo)準(zhǔn)出版社.

[3]胡昌華，張軍波，夏 軍，等.基于Matlab的系統(tǒng)分析與設(shè)計(jì)-小波分析[M].西安：西安電子科技大學(xué)出版社，1999.

[4]邢西淳，毛 娟.FHD-1數(shù)字化地磁觀測(cè)系統(tǒng)常見故障維修及電源改進(jìn)[J].地震地磁觀測(cè)與研究，2008， 29(6)： 131-138.

[5]Daubechies I.Ten lectures on wavelets[M].CBMS-NSF Regional Conferences Series in Applied Mathematics,1992.

[6]謝凡，滕云田，胡星星，等.地磁臺(tái)的城市軌道交通干擾的小波抑制方法研究[J].地球物理學(xué)報(bào)，2011， 54(10)： 2698-2707.

[7]Asimopolos L,Pestina A M,Asimopolos N S.Considerations on Geomagnetic Data Analysis[J].Chinese J.Geophys， 2010， 53(3)： 765-772.

[8]潘泉，孟晉麗，張磊，等.小波濾波方法及應(yīng)用[J].電子與信息學(xué)報(bào)，2007，29(1)：236-242.

[9]邵輝成，杜興信，金學(xué)申，等.小波分析在地震趨勢(shì)預(yù)測(cè)中的應(yīng)用[J].中國地震，2000，16(1)：48～52.

[10]邵輝成，杜常娥.中國大陸地震活動(dòng)的多尺度分析[J].地震學(xué)報(bào)，2004，26（1）：102-105.

[11]Bruce A G,Gao H Y.Understanding Wave Shrink:Variance and bias estimation[J].Biometrika.1996，83(4)：727-745.

[12]Gao H Y,Gao H,Bruce A G.WaveShrink with semisoft shrinkage[J].StatSci Division of MathSoft，Inc.1995.

[13]Gao H Y.Wavelet shrinkage denoising using the non-negative garrote[J.Journalof Computational and Graphical Statistics.1998， 7(4)： 469-488.

[14]高靜懷，朱光明.地震資料處理中小波函數(shù)的選取研究[J].地球物理學(xué)報(bào)，1996，39(3)：392～400.