基于小波包降噪的地磁場干擾數據分析

徐斌， 顧偉

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

0 引 言

在現代海軍發展中，測磁與消磁已經成為關系到艦艇生存的重要技術.港口消磁站的建設也已成為海軍艦艇基地的重要項目.在港口消磁站測量的信號中除艦艇的磁場外還包含地磁場和其他環境或鐵磁物質的信號，因此要進行信號分析降噪，才能準確地得到所需的磁場信號.

在現代地磁測量中，有不少可行的降噪方法在地震電磁檢測站得到實際應用[1]，其中比較有效的為小波降噪算法.然而在現有的港口消磁站測磁環節研究中，還沒有如何降噪的應用算法.

本文在研究小波算法與小波包算法[2]的異同和優點后，對消磁監測站測量的艦船磁場原始信號運用小波包降噪算法進行分析，分離出地磁日變干擾磁場及港口消磁站附近日常作業產生的干擾磁場，從而得到平穩的地磁測量信號.

1 測量實驗平臺

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

實際測量磁場存在干擾，因此需要進行數據處理濾除干擾磁場.

磁場實測值=被測對象磁場+人為干擾磁場+背景地磁場

人為干擾磁場=港口日常作業的干擾磁場+港口船只出入的信號瞬變干擾磁場

背景地磁場=本地地磁磁場+地磁日變磁場

地磁日變磁場基本可分為地磁每日都存在的地磁靜日變化和地磁擾日變化,其中：地磁靜日變化為具有周期性的連續磁場變化;地磁擾日變化為無周期性的隨機磁場變化，不是恒定值.因此在實際測量中，將背景地磁場中的地磁日變磁場作為干擾量.本地地磁磁場為由地殼以下地核場、地殼場與本地周圍固定建筑群的磁場之和，為測量時基本固定的本地環境磁場.因此，實際測量中的干擾磁場為人為干擾磁場與地磁日變磁場之和.

在實際測量中干擾信號一直存在，因此先對無艦船時的環境磁場進行測量和分析.在實際測量中對環境中的地磁日變磁場以及正常運行的港口附近高壓電、電機等干擾磁場進行25 h不間斷監測.數據由磁通門傳感器陣列采集，采樣時間自第一天10:40起至次日11:40結束，每分鐘采樣6次，得到磁通門傳感器陣列磁場數據.

2 干擾磁場分析

2.1 地磁日變磁場

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

圖1 測點附近地磁日變曲線示例

2.2 港口測磁點附近港口作業產生的干擾磁場

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

(a)X分量磁場

(b)Y分量磁場

(c)Z分量磁場

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

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

3 小波包降噪

3.1 小波包分析重構

小波包變換由小波變換發展而來，可以視為普通小波函數的線性組合.小波變換從空間(時間)和頻率的局部變換中提取信號中有用的信息.[5]小波變換的本質是低通濾波,然而其直接濾除高頻信號的特點會使一些有用的高頻信號細節丟失.小波包變換在處理信號時同時分解高頻與低頻部分，因此更具有靈活性.因此，對含噪信號進行小波包變換，對其小波包系數進行閾值操作，然后進行重構，得到的消噪后信號優于小波變換的處理結果[6].

小波包降噪步驟為：信號的小波包分解、確定最佳小波包基、小波包分解系數的閾值量化、小波包信號重構.圖3為小波包分解樹，其中:S為原信號，f為信號頻率，A代表低頻段，D代表高頻段.由圖3可以看出，原信號S可以進行多層分解，直到其分解的小波包能夠滿足分析需求.第1層分解為尺度參數為1的互不重疊的高頻D((fl+fh)/2,fh)與低頻A(fl,(fl+fh)/2)兩部分，原信號S的完整尺度(fl,fh)被分割.第2層分解得到尺度參數為2的4個頻率范圍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).如果進行n層分解，可以得到尺度參數為n的、2n個節點的小波包分解樹.圖3中進行的2層分解僅用于演示說明，實際分解層數需按信號計算得到最佳分解樹才能確定.

圖3 小波包分解樹狀圖

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

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

(1)

(2)

(3)

小波包的重構公式為

(4)

3.2 閾值的選取

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

小波包信號閾值去噪方法通常有：固定閾值、自適應閾值、混合型閾值和最小最大準則閾值.實際工程中，軟閾值和硬閾值應用最多.[12]軟閾值模型為

(5)

硬閾值模型為

(6)

式中：λ為閾值；ω為原信號節點能量元.

本文采用一種新的閾值函數，并用MATLAB進行新閾值函數下的小波包信號去噪運算[13].

(7)

為衡量上述3種閾值消噪的效果，通常采用信噪比和均方差作為衡量標準.信噪比

(8)

均方差

(9)

式中:I和Id分別為原始信號和消噪后的信號.

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

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

3.3 小波包MATLAB實現

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

>>save mag.mat S；%將地磁數據信號S存儲在mag.mat文件中；

>>load mag；%將文件裝載到matlab環境中；

>>S=A(1：8860)；%將信號的8860個數據給變量S；

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

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

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

>> s130=wprcoef(wpt,[3,0]); %重構節點(3,0);

…

>> s137=wprcoef(wpt,[3,7]); %重構節點(3,7);

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

上述程序中的S為對每個節點系數設置閾值后得到保留的有用系數.

4 數據處理結果

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

(a) 軟閾值方法

(b) 硬閾值方法

(a) X分量磁場

(b) Y分量磁場

(c) Z分量磁場

圖6 理想狀態下的Z分量磁場數據

圖7 50 min磁場數據

5 結 論

小波包分解重構算法在實際應用中對港口正常作業所引起的干擾磁場的抑制效果比較明顯，對測量艦艇信號，降低由地磁日變磁場、港口正常作業所引起的地磁干擾起到積極作用.不過小波包算法對信號進行降噪處理時計算量較大，具有一定的延遲，如何實現高度實時性還需在以后的研究中繼續探索.

參考文獻：

[1]劉宗堯. 地磁信號檢測系統設計及誤差補償研究[D]. 南京: 南京理工大學, 2011.

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

[3]郭愛煌. 磁通門技術及其應用[J]. 傳感器技術, 2000, 19(4): 1-4.

[4]卞光浪, 邊剛, 于波, 等. 基于多分辨率分析地磁日變信號去噪方法[C]//中國測繪學會第十八屆海洋測繪綜合性學術研討會論文集. 2010: 66-68.

[5]韓朝暉. 利用小波包分析提取信號分量[C]//2006北京地區高校研究生學術交流會——通信與信息技術會議論文集(上). 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]任大男, 任韌. 自適應小波包正交性系統研究[J]. 應用光學, 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]. 天津大學學報, 2010, 43(6): 562-566.

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

[13]曲國慶, 黨亞民, 章傳銀, 等. 小波包消噪方法分析及改進[J]. 大地測量與地球動力學, 2008, 28(4): 102-106.

[14]王嘉梅. 基于MATLAB的小波包信號消噪處理[J]. 計算機與網絡, 2001(1): 52-53.