徐斌， 顧偉

(上海海事大學(xué) 物流工程學(xué)院，上海 201306)

0 引 言

在現(xiàn)代海軍發(fā)展中，測(cè)磁與消磁已經(jīng)成為關(guān)系到艦艇生存的重要技術(shù).港口消磁站的建設(shè)也已成為海軍艦艇基地的重要項(xiàng)目.在港口消磁站測(cè)量的信號(hào)中除艦艇的磁場(chǎng)外還包含地磁場(chǎng)和其他環(huán)境或鐵磁物質(zhì)的信號(hào)，因此要進(jìn)行信號(hào)分析降噪，才能準(zhǔn)確地得到所需的磁場(chǎng)信號(hào).

在現(xiàn)代地磁測(cè)量中，有不少可行的降噪方法在地震電磁檢測(cè)站得到實(shí)際應(yīng)用[1]，其中比較有效的為小波降噪算法.然而在現(xiàn)有的港口消磁站測(cè)磁環(huán)節(jié)研究中，還沒(méi)有如何降噪的應(yīng)用算法.

本文在研究小波算法與小波包算法[2]的異同和優(yōu)點(diǎn)后，對(duì)消磁監(jiān)測(cè)站測(cè)量的艦船磁場(chǎng)原始信號(hào)運(yùn)用小波包降噪算法進(jìn)行分析，分離出地磁日變干擾磁場(chǎng)及港口消磁站附近日常作業(yè)產(chǎn)生的干擾磁場(chǎng)，從而得到平穩(wěn)的地磁測(cè)量信號(hào).

1 測(cè)量實(shí)驗(yàn)平臺(tái)

在港口測(cè)磁過(guò)程中采用磁通門(mén)傳感器[3]陣列采集磁場(chǎng)信號(hào).磁通門(mén)傳感器利用材料的B-H飽和特性進(jìn)行弱磁場(chǎng)測(cè)量，其測(cè)量范圍在10-11～10-2T，能夠有效捕捉極弱磁場(chǎng)的特征信號(hào).在港口測(cè)磁過(guò)程中，磁通門(mén)采用多傳感器(皆為三分量傳感器)陣列方式，Z軸平行于地磁Z軸方向.

實(shí)際測(cè)量磁場(chǎng)存在干擾，因此需要進(jìn)行數(shù)據(jù)處理濾除干擾磁場(chǎng).

磁場(chǎng)實(shí)測(cè)值=被測(cè)對(duì)象磁場(chǎng)+人為干擾磁場(chǎng)+背景地磁場(chǎng)

人為干擾磁場(chǎng)=港口日常作業(yè)的干擾磁場(chǎng)+港口船只出入的信號(hào)瞬變干擾磁場(chǎng)

背景地磁場(chǎng)=本地地磁磁場(chǎng)+地磁日變磁場(chǎng)

地磁日變磁場(chǎng)基本可分為地磁每日都存在的地磁靜日變化和地磁擾日變化,其中：地磁靜日變化為具有周期性的連續(xù)磁場(chǎng)變化;地磁擾日變化為無(wú)周期性的隨機(jī)磁場(chǎng)變化，不是恒定值.因此在實(shí)際測(cè)量中，將背景地磁場(chǎng)中的地磁日變磁場(chǎng)作為干擾量.本地地磁磁場(chǎng)為由地殼以下地核場(chǎng)、地殼場(chǎng)與本地周?chē)潭ńㄖ旱拇艌?chǎng)之和，為測(cè)量時(shí)基本固定的本地環(huán)境磁場(chǎng).因此，實(shí)際測(cè)量中的干擾磁場(chǎng)為人為干擾磁場(chǎng)與地磁日變磁場(chǎng)之和.

在實(shí)際測(cè)量中干擾信號(hào)一直存在，因此先對(duì)無(wú)艦船時(shí)的環(huán)境磁場(chǎng)進(jìn)行測(cè)量和分析.在實(shí)際測(cè)量中對(duì)環(huán)境中的地磁日變磁場(chǎng)以及正常運(yùn)行的港口附近高壓電、電機(jī)等干擾磁場(chǎng)進(jìn)行25 h不間斷監(jiān)測(cè).數(shù)據(jù)由磁通門(mén)傳感器陣列采集，采樣時(shí)間自第一天10:40起至次日11:40結(jié)束，每分鐘采樣6次，得到磁通門(mén)傳感器陣列磁場(chǎng)數(shù)據(jù).

2 干擾磁場(chǎng)分析

2.1 地磁日變磁場(chǎng)

地磁日變磁場(chǎng)為一個(gè)太陽(yáng)日中主要由固體地球外部原因引起的、疊加在地球基本磁場(chǎng)之上的各種短期的地磁變化.按照成因不同，變化磁場(chǎng)可分為平靜變化和干擾變化兩大類(lèi).[4]前者的基本成因是電離層在地磁場(chǎng)中運(yùn)動(dòng)產(chǎn)生較為穩(wěn)定的電流體系;后者的基本成因是太陽(yáng)風(fēng)與地磁場(chǎng)相互作用，在磁層和電離層中形成各種短暫的電流體系.地球外部的各種電流體系，都能在固體地球內(nèi)部感應(yīng)出相應(yīng)的內(nèi)部電流體系.每一種地磁變化都是外部電流體系與內(nèi)部電流體系產(chǎn)生的磁場(chǎng)之和,一般前者約占總和的70%，后者約占30%.靜日變化Sq雖然依賴(lài)于當(dāng)?shù)靥?yáng)時(shí), 并以一個(gè)太陽(yáng)日為周期, 但它在時(shí)空分布特征上是有差異的, 而且緯度不同的地區(qū)靜日變化規(guī)律不同.靜日變化的日變曲線在形態(tài)上基本一致,白晝起伏大,夜間相對(duì)穩(wěn)定,極大值和極小值主要出現(xiàn)在白天.擾日變化為無(wú)規(guī)律的隨機(jī)磁場(chǎng)變化，雖然在地磁日變磁場(chǎng)中僅占小部分，即使在地磁活動(dòng)十分平靜期間，也或多或少地存在擾日變化.圖1為測(cè)點(diǎn)附近地磁日變曲線示例，數(shù)據(jù)來(lái)源于測(cè)點(diǎn)附近地磁監(jiān)測(cè)站公示的地磁日變采樣數(shù)據(jù)，橫坐標(biāo)為測(cè)試地點(diǎn)的時(shí)間(當(dāng)?shù)貢r(shí))軸，且對(duì)應(yīng)的每3 h的Kp指數(shù)值均≤2.圖中所示的磁場(chǎng)為：在測(cè)試地點(diǎn)周?chē)鷧^(qū)域存在鐵磁性物質(zhì)，但測(cè)試期間無(wú)鐵磁性物質(zhì)在附近移動(dòng)及所在區(qū)域無(wú)強(qiáng)電供應(yīng)的情況下測(cè)量得到的磁場(chǎng).由Kp指數(shù)可以看出，圖中所示測(cè)試時(shí)段內(nèi)，測(cè)點(diǎn)附近的地磁變化的主要影響因素為靜日變化,擾日變化雖然存在但相對(duì)較為平靜.圖中的磁場(chǎng)隨時(shí)間進(jìn)行著緩慢、平穩(wěn)的變化，總體變化趨勢(shì)具有一定的規(guī)律,短時(shí)變化為無(wú)規(guī)律的、幅值≤±10 nT的、隨機(jī)變化的白噪聲.

圖1 測(cè)點(diǎn)附近地磁日變曲線示例

2.2 港口測(cè)磁點(diǎn)附近港口作業(yè)產(chǎn)生的干擾磁場(chǎng)

港口測(cè)磁點(diǎn)附近各種設(shè)施的運(yùn)作所用的高壓電由岸電供應(yīng).由港口作業(yè)大型機(jī)械、岸電供電站、岸電供電線路、港口行車(chē)軌道、大地形成電流回路.在港口行車(chē)軌道與大地不完全絕緣時(shí)，會(huì)造成電流部分泄漏.這部分泄漏電流本身會(huì)產(chǎn)生磁場(chǎng)且供電電路中電流不平衡會(huì)形成磁場(chǎng).而磁傳感器的安裝位置離港口只有幾公里，因此測(cè)量時(shí)會(huì)受到港口磁場(chǎng)的影響.圖2為一組從早上10:40開(kāi)始測(cè)量的25 h三分量磁場(chǎng)數(shù)據(jù)圖，橫坐標(biāo)為采樣序列，每分鐘采樣6次.從圖中可以看出，在三分量磁傳感器中，在上午10:40至午夜12:00磁場(chǎng)幅值波動(dòng)強(qiáng)烈，雖然Z分量磁場(chǎng)還是依照日變曲線的變化規(guī)律浮動(dòng)，但是其幅值振蕩無(wú)規(guī)律，其振幅比午夜12:00至早上6:00擴(kuò)大十幾倍，X與Y分量受到附近電磁干擾影響較小.磁場(chǎng)變化特征符合日常港口作業(yè)時(shí)間特點(diǎn).因此可以斷定，消磁站附近的日常作業(yè)產(chǎn)生的電磁干擾是存在且明顯的.

(a)X分量磁場(chǎng)

(b)Y分量磁場(chǎng)

(c)Z分量磁場(chǎng)

2.3 消磁站附近港口船只出入瞬變干擾

港口日常作業(yè)時(shí)，會(huì)不時(shí)有船只進(jìn)出港.船只靠近磁傳感器時(shí)勢(shì)必會(huì)引起磁場(chǎng)的異變，其特征為瞬變、幅值大、時(shí)間短，在數(shù)據(jù)圖上表現(xiàn)為一個(gè)突變峰值，很容易識(shí)別和抑制處理.在本次實(shí)驗(yàn)中，臨時(shí)禁止任何船只進(jìn)出港口，使港口船只出入引起的瞬變干擾得到完全抑制.

3 小波包降噪

3.1 小波包分析重構(gòu)

小波包變換由小波變換發(fā)展而來(lái)，可以視為普通小波函數(shù)的線性組合.小波變換從空間(時(shí)間)和頻率的局部變換中提取信號(hào)中有用的信息.[5]小波變換的本質(zhì)是低通濾波,然而其直接濾除高頻信號(hào)的特點(diǎn)會(huì)使一些有用的高頻信號(hào)細(xì)節(jié)丟失.小波包變換在處理信號(hào)時(shí)同時(shí)分解高頻與低頻部分，因此更具有靈活性.因此，對(duì)含噪信號(hào)進(jìn)行小波包變換，對(duì)其小波包系數(shù)進(jìn)行閾值操作，然后進(jìn)行重構(gòu)，得到的消噪后信號(hào)優(yōu)于小波變換的處理結(jié)果[6].

小波包降噪步驟為：信號(hào)的小波包分解、確定最佳小波包基、小波包分解系數(shù)的閾值量化、小波包信號(hào)重構(gòu).圖3為小波包分解樹(shù)，其中:S為原信號(hào)，f為信號(hào)頻率，A代表低頻段，D代表高頻段.由圖3可以看出，原信號(hào)S可以進(jìn)行多層分解，直到其分解的小波包能夠滿(mǎn)足分析需求.第1層分解為尺度參數(shù)為1的互不重疊的高頻D((fl+fh)/2,fh)與低頻A(fl,(fl+fh)/2)兩部分，原信號(hào)S的完整尺度(fl,fh)被分割.第2層分解得到尺度參數(shù)為2的4個(gè)頻率范圍AA(fl, (fl+fh)/4)，AD((fl+fh)/4, (fl+fh)/2)，DA((fl+fh)/2, 3(fl+fh)/4)，DD(3(fl+fh)/4,fh).如果進(jìn)行n層分解，可以得到尺度參數(shù)為n的、2n個(gè)節(jié)點(diǎn)的小波包分解樹(shù).圖3中進(jìn)行的2層分解僅用于演示說(shuō)明，實(shí)際分解層數(shù)需按信號(hào)計(jì)算得到最佳分解樹(shù)才能確定.

圖3 小波包分解樹(shù)狀圖

信號(hào)以f(c)正交小波分解的公式[7]為

Pj-1f(t)=Pjf(t)+Djf(t)

(1)

(2)

(3)

小波包的重構(gòu)公式為

(4)

3.2 閾值的選取

閾值如何選擇、小波包系數(shù)如何進(jìn)行閾值量化一直是小波包降噪算法核心中的核心，直接關(guān)系到信號(hào)去噪的品質(zhì).一般數(shù)據(jù)的數(shù)組序列中,xi=fi+ei(i=1,2,3,…,n),其中,fi為信號(hào)序列部分，ei為高斯噪聲序列，且ei與fi互不相關(guān).xi的性質(zhì)可用它的小波包系數(shù)描述.小波包系數(shù)反映信號(hào)的能量，因此通常將小波包系數(shù)看作“能量元”，其系數(shù)大小代表其攜帶能量的多少.[8]噪聲中信號(hào)信息的小波包系數(shù)值一般小于fi小波包系數(shù)值.噪聲小波“能量元”在經(jīng)過(guò)正交鏡像濾波分解后，變成在整個(gè)小波包系數(shù)軸的均勻分布.[9]因此,可以選擇一個(gè)閾值，濾除所有等于或小于閾值的小波包系數(shù)，保留所有大于閾值的小波包系數(shù)，用以重構(gòu)信號(hào).[10]在選取閾值的過(guò)程中，閾值太大不僅會(huì)濾除噪聲信息而且會(huì)濾除有效的原有信號(hào)細(xì)節(jié)部分，使重構(gòu)得到的信號(hào)損失部分特征甚至完全失真;閾值太小會(huì)使噪聲信號(hào)得不到濾除，無(wú)法達(dá)到預(yù)期的降噪效果.[11]由此可見(jiàn)，不同閾值所形成的信噪比截然不同，因此閾值的選取至關(guān)重要，需在實(shí)驗(yàn)中不斷嘗試，得到最優(yōu)閾值后再進(jìn)行降噪處理.

小波包信號(hào)閾值去噪方法通常有：固定閾值、自適應(yīng)閾值、混合型閾值和最小最大準(zhǔn)則閾值.實(shí)際工程中，軟閾值和硬閾值應(yīng)用最多.[12]軟閾值模型為

(5)

硬閾值模型為

(6)

式中：λ為閾值；ω為原信號(hào)節(jié)點(diǎn)能量元.

本文采用一種新的閾值函數(shù)，并用MATLAB進(jìn)行新閾值函數(shù)下的小波包信號(hào)去噪運(yùn)算[13].

(7)

為衡量上述3種閾值消噪的效果，通常采用信噪比和均方差作為衡量標(biāo)準(zhǔn).信噪比

(8)

均方差

(9)

式中:I和Id分別為原始信號(hào)和消噪后的信號(hào).

表1 3種閾值消噪的θ和ρ

表1列出上述實(shí)驗(yàn)中降噪后信號(hào)的信噪比θ以及原始信號(hào)與消噪后的信號(hào)之間的均方差ρ.由表1可以看出，改進(jìn)方案的降噪效果明顯比單純的軟閾值或硬閾值方法好，從而驗(yàn)證改進(jìn)方案的有效性.

3.3 小波包MATLAB實(shí)現(xiàn)

MATLAB小波包程序[14]如下：

>>save mag.mat S；%將地磁數(shù)據(jù)信號(hào)S存儲(chǔ)在mag.mat文件中；

>>load mag；%將文件裝載到matlab環(huán)境中；

>>S=A(1：8860)；%將信號(hào)的8860個(gè)數(shù)據(jù)給變量S；

>>wpt = wpdec(S,3,‘db1’, ‘user’) ；%使用db1小波包對(duì)S進(jìn)行3層分解；

>>wpt = wpsplt(wpt,[3 0]); %分解小波包節(jié)點(diǎn)(3,0);

>>cfs = wpcoef(wpt,[3 0]); %讀取小波包(3,0)的系數(shù);

>> s130=wprcoef(wpt,[3,0]); %重構(gòu)節(jié)點(diǎn)(3,0);

…

>> s137=wprcoef(wpt,[3,7]); %重構(gòu)節(jié)點(diǎn)(3,7);

>> subtlot(311);tlot(a1);

上述程序中的S為對(duì)每個(gè)節(jié)點(diǎn)系數(shù)設(shè)置閾值后得到保留的有用系數(shù).

4 數(shù)據(jù)處理結(jié)果

圖4為采用軟、硬閾值方法去噪后的Z分量磁場(chǎng)數(shù)據(jù).圖5為采用新閾值方法去噪后的三分量磁場(chǎng)數(shù)據(jù).比較圖4和5的Z分量數(shù)據(jù)可知，新閾值方法比普通軟閾值和硬閾值方法降噪效率更好，能得到更好的信噪比.從圖5可以看出，濾波后的信號(hào)在保持地磁的日變形態(tài)基本不變的同時(shí)，大幅降低由港口正常作業(yè)引起的電磁干擾幅度.圖6為理想無(wú)干擾Z分量磁場(chǎng)數(shù)據(jù)圖，數(shù)據(jù)來(lái)源于當(dāng)日地磁站公示的日變磁場(chǎng)查表數(shù)據(jù).從圖5和6中可以得出，降噪處理后的數(shù)據(jù)基本符合理想實(shí)測(cè)數(shù)據(jù)，說(shuō)明小波包降噪能在有效降噪的同時(shí)，比較完整地保留原始數(shù)據(jù)中的有效數(shù)據(jù).圖7為從處理后數(shù)據(jù)中任意截取的50 min地磁數(shù)據(jù)圖.從圖7中可以看出,地磁日變的白噪聲可以控制在5 nT之內(nèi).因此可以得出結(jié)論：小波包降噪算法能夠較為有效地抑制人為干擾,并且能對(duì)地磁日變干擾磁場(chǎng)起到一定的抑制效果.

(a) 軟閾值方法

(b) 硬閾值方法

(a) X分量磁場(chǎng)

(b) Y分量磁場(chǎng)

(c) Z分量磁場(chǎng)

圖6 理想狀態(tài)下的Z分量磁場(chǎng)數(shù)據(jù)

圖7 50 min磁場(chǎng)數(shù)據(jù)

5 結(jié) 論

小波包分解重構(gòu)算法在實(shí)際應(yīng)用中對(duì)港口正常作業(yè)所引起的干擾磁場(chǎng)的抑制效果比較明顯，對(duì)測(cè)量艦艇信號(hào)，降低由地磁日變磁場(chǎng)、港口正常作業(yè)所引起的地磁干擾起到積極作用.不過(guò)小波包算法對(duì)信號(hào)進(jìn)行降噪處理時(shí)計(jì)算量較大，具有一定的延遲，如何實(shí)現(xiàn)高度實(shí)時(shí)性還需在以后的研究中繼續(xù)探索.

參考文獻(xiàn)：

[1]劉宗堯. 地磁信號(hào)檢測(cè)系統(tǒng)設(shè)計(jì)及誤差補(bǔ)償研究[D]. 南京: 南京理工大學(xué), 2011.

[2]李鳳閣, 李鵬. 基于小波和小波包的信號(hào)消噪方法分析[J]. 山西科技, 2010, 25(4): 70-73.

[3]郭愛(ài)煌. 磁通門(mén)技術(shù)及其應(yīng)用[J]. 傳感器技術(shù), 2000, 19(4): 1-4.

[4]卞光浪, 邊剛, 于波, 等. 基于多分辨率分析地磁日變信號(hào)去噪方法[C]//中國(guó)測(cè)繪學(xué)會(huì)第十八屆海洋測(cè)繪綜合性學(xué)術(shù)研討會(huì)論文集. 2010: 66-68.

[5]韓朝暉. 利用小波包分析提取信號(hào)分量[C]//2006北京地區(qū)高校研究生學(xué)術(shù)交流會(huì)——通信與信息技術(shù)會(huì)議論文集(上). 2006: 582-589.

[6]MALLAT S G. A theory for multiresolution signal decomposition: the wavelet representation[J]. Pattern Anal & Machine Intelligence, IEEE Trans, 1989, 11(7): 674-693.

[7]ZHANG Xuming, XIONG Youlun. Impulse noise removal using directional difference based noise detector and adaptive weighted mean filter[J]. Signal Processing Letters, IEEE, 2009, 16(4): 295-298.

[8]CHUI C K, LIAN J. A study of orthonormal multi-wavelets[J]. Appl Numerical Math, 1996, 20(3): 273-298.

[9]任大男, 任韌. 自適應(yīng)小波包正交性系統(tǒng)研究[J]. 應(yīng)用光學(xué), 2010, 31(S): 61-65.

[10]YUAN Jing, HE Zhengjia, ZI Yanyang,etal. Adaptive multiwavelets via two-scale similarity transforms for rotating machinery fault diagnosis[J]. Mech Systems & Signal Processing, 2009, 23(5): 1490-1508.

[11]孫潔娣, 靳世久. 基于小波包能量及高階譜的特征提取方法[J]. 天津大學(xué)學(xué)報(bào), 2010, 43(6): 562-566.

[12]魯懷偉, 杜三山. 一種小波包去噪自適應(yīng)閾值算法[J]. 蘭州鐵道學(xué)院學(xué)報(bào): 自然科學(xué)版, 2001, 20(6): 11-15.

[13]曲國(guó)慶, 黨亞民, 章傳銀, 等. 小波包消噪方法分析及改進(jìn)[J]. 大地測(cè)量與地球動(dòng)力學(xué), 2008, 28(4): 102-106.

[14]王嘉梅. 基于MATLAB的小波包信號(hào)消噪處理[J]. 計(jì)算機(jī)與網(wǎng)絡(luò), 2001(1): 52-53.