基于干擾檢測和CEEMD的地磁信號降噪方法

王立輝，劉慶雅，許寧徽

（東南大學 儀器科學與工程學院 微慣性儀表與先進導航技術教育部重點實驗室，南京 210096）

地磁導航是運載器自主導航的重要研究方向之一[1]。地磁測量信息作為地磁導航的信號源，其測量精度與導航系統的精度和穩定性有關。由于地磁隨位置變化幅度較小，環境磁噪聲以及載體工作狀態變化引起的高幅值磁脈沖干擾將導致磁通門傳感器測量信號與先驗地磁圖存在較大的干擾誤差。在高幅值磁脈沖干擾信號中，載體工作狀態變化產生瞬態的磁脈沖干擾，類似椒鹽噪聲，可通過中值濾波去除；而不規則地磁脈動引起的環境磁噪聲脈動幅值范圍0.1 nT～1000 nT，頻率在1 mH～10 Hz 之間[2]，部分頻譜與地磁異常場重疊，常規的濾波算法無法將地磁信號從其中分離。

目前，多數算法針對地磁測量中的隨機誤差進行抑制[3-5]，以隨機誤差為高斯噪聲為基礎。當噪聲非高斯分布時，誤差抑制效果較差。中值濾波、均值濾波等針對干擾噪聲的濾波方法在去除圖片椒鹽噪聲、GPS 脈沖干擾等瞬態脈沖干擾上有著廣泛的應用[6]。由于磁場的磁化效應，磁脈沖干擾往往持續時間差異較大[7]，使用單一中值濾波難以抑制磁脈沖干擾對地磁測量的誤差。完備集成經驗模態分解（Complete ensemble empirical mode decomposition，CEEMD）具有自適應性、完備性和近似正交性的特點，在濾除異常光纖陀螺信號和異常GNSS 信號等方面有著廣泛的應用[8,9]。當CEEMD 分解含有高幅值脈沖干擾的測量信號，會出現模態混疊現象，干擾信號分布在各本征模態分量（IMF）中，需要對其進行干擾檢測。Bandt等[10]提出了一種檢測時間序列隨機性和動態變化的有效方法，來評價時間序列的自然復雜度，即排列熵。與傳統熵計算方法不同，排列熵考慮信號的時間順序，計算時間序列相鄰值的相關性，可以表現出信號隨時間變化的異常部分[11]。近年來，排列熵在抑制心電圖測量干擾、環境數據異常監測等領域有著廣泛的應用[12,13]。

為了有效濾除常規隨機噪聲與高幅值磁脈沖干擾，本文首先分析了地磁信號特性與干擾信號特性，利用地磁信號的強相關性，計算各IMF的自相關指數，并自適應地選取有效IMF。針對高幅值脈沖磁干擾，提出了一種多尺度的干擾檢測方法，引入排列熵計算有效IMF 的時間復雜度，并對復雜度高的部分進行多尺度分析，檢測干擾成分，并濾除干擾部分后進行信號重構，獲得高信噪比的地磁信號。

1 地磁信號測量及特征分析

地磁場一般在幾nT 至幾萬nT 之間，屬于弱磁場，一般選用磁通門、質子或光泵磁力儀進行測量。磁通門磁力儀能測量磁場矢量信息，具有動態性能好、采樣頻率高、成本低、體積小和低功耗等優點，在導航系統中應用廣泛[14]。磁通門傳感器的測量原理導致其測量的磁信息并非絕對值，并且溫度效應差，測量精度低。一般通過誤差補償與干擾濾波，提高磁通門傳感器的測量性能。

磁通門傳感器測量誤差主要包括三部分：傳感器固有誤差、航磁干擾與隨機磁干擾[15]，磁傳感器的磁測值可以表示為：

其中，Bm為磁傳感器測量值，Q為非正交補償矩陣，Bg為地磁矢量，Bh為航磁干擾，Be隨機磁干擾，E為傳感器零偏。針對傳感器非正交誤差和零偏，常采用預先標定進行校正[16]，針對穩定航磁干擾，常假設載體為均勻磁體，且剛性連接，載體機動補償所在區域地磁場恒定，使用T-L 方程進行穩定磁干擾補償。在沒有配備磁屏蔽設備的環境下校準，仍受到不可避免的外部磁場干擾。在實際應用中，除了上述的穩定磁場干擾，相對地磁基準圖，仍存在環境與載體的隨機高幅值磁干擾，如地質運動、地磁脈動等[17,18]。圖1為某地磁臺日測值曲線，地磁測量值存在較大波動，波動幅值在1000 nT左右，頻率在1 mHz～100 Hz之間，對導航系統穩定性有很大影響。

圖1 某地磁臺日測量值Fig.1 Daily measurement of a geomagnetic station

2 干擾檢查與完備集成模式分解（ID-CEEMD）

由于磁干擾頻譜特征與地磁圖測量信號頻譜接近，傳統的CEEMD 無法有效地將高幅值長周期磁干擾的測量信號分解為獨立的IMF。為解決該問題，提出了一種基于排列熵的多尺度干擾檢測方法，結合CEEMD 的正交性，對干擾段進行二次濾波，取低頻殘余信號進行信號重構。

2.1 CEEMD

經驗模態分解（EMD）是一種針對非線性信號的自適應濾波方法，它使用三次樣條插值法將信號分解為具有不同頻率的IMF。為了獲得目標信號，可以通過分類規則將IMF 分類為以信號為主的IMF 和以噪聲為主的IMF，并通過以信號為主的IMF 對目標信號進行重構。

CEEMD 是在EMD 分解的基礎上，通過對原信號多次疊加正負對應的白噪聲，補充原始信號的頻率尺度，減少包絡線的擬合誤差。具體步驟如下[19]：

1）對原信號第i次疊加正負高斯白噪聲，得到兩組模態分量：

其中，X為原始信號序列，Ni為第i次添加的白噪聲序列，

2）分解第i次添加白噪聲后的信號M1i,M2i，得到對應模態分量對應余量r1i、r2i。

3）對各個分量取平均，得到

4）對n次疊加白噪聲產生的分量對應相加，求其平均值，最終得到CEEMD 分解結果：

其中，cj為CEEMD 分解后的第j個本征模態分量。

在實際測量中，磁通門傳感器會受到很多不確定的電磁脈沖的干擾，例如局部異常電磁干擾，車輛工作狀態和姿態的變化。當EMD 分解這些高幅值脈沖信號時，會出現模態混疊現象。模態混疊是指地磁信號和干擾信號在IMF 中混合，而不能在時域中分離。盡管CEEMD 在某種程度上抑制了模態混疊，但是當原始信號的信噪比（SNR）低時，也會出現模態混疊現象。

2.2 自適應有效IMF 分離

經過CEEMD 分解后，原始信號被分解為多個IMF 和殘余分量。分解得到的IMF 需要進行分類以獲得有效IMF。自相關函數反映了函數在不同時間的相關程度，在地磁測量信號中，有效信號為超低頻分量，隨時間變換緩慢，具有較高的自相關性。設一個隨機信號為x(t)，其自相關函數為

其中，τ為信號偏移時間，T為信號長度。在零偏移處，信號的自相關函數達到最大值。理想高斯白噪聲的歸一化自相關函數在零偏移處的值為1，其余偏移為0[20]。與干擾IMF 相比，信息IMF 的自相關函數具有更均勻的能量分布。為了比較多個信號的自相關性，需要對每個信號的自相關函數進行歸一化。歸一化自相關函數為

歸一化自相關函數避免了信號能量對標準偏差的影響，并比較了信號之間的自相關。在本文中，將每個模態分量的歸一化的歸零自相關函數標準偏差用作IMF 的評估參數，可以計算如下：

其中，P(i)為第i個IMF 的標準差，ρi(τ)為第i個IMF的歸一化自相關函數，為ρ i(τ)的平均值，T為信號持續時間，閾值計算如下：

當P(i)>μ時，對應的IMF 為信號IMF，當P(i)≤μ時，對應的IMF 為噪聲IMF。

2.3 多尺度干擾檢測與濾波

測量信號被CEEMD 分解為多個IMF 后，自適應篩選方法對IMF 進行有效性分類，得到多個有效IMF，此時單個IMF 既包含基本測量信號又包含磁干擾值。對有效IMF 進行干擾檢測，并對干擾段進行濾波，可以有效濾除磁干擾。

排列熵僅利用值的順序，因此排列熵在信號的非線性失真下具有魯棒性，計算效率高。但是排列熵忽略了幅度的平均值和幅度值之間的差異[21]。結合地磁測量干擾信號的特點，對測量信號分解的有效IMF 進行多尺度干擾檢測。

圖2 排列熵滑動窗示意圖Fig.2 The sliding windows schematic

設IMF 的離散時間序列為x(n)，長度為N；窗口元素長度為M，滿足N >> M，窗口重疊率為ρ，x s和xs+1的離散時間距離為φ= (1 -ρ)M，xs與xs+1如圖1所示，為兩個相鄰的窗口。xs+1和xs的幅值變化率可以表示為：

隨著窗口的滑動，x(n)可以被分為xi(m)，i= 1,2,…,L，L= int((N-M)/φ)，int（）表示向下舍入。x(n)到xi(m)的映射可以表示為：

x(n)置換熵di可以表示為：

根據Z-score 標準化，如下所示：

其中，μ是di的平均值，σ為di的標準差。

設置置換熵閾值為θ，將信號分為干擾段和正常段。如圖1所示，ds是干擾段的開始，對應到時間序列，干擾段的開始為x s+1(m)，m= 1，de為干擾段的結束，對應到時間序列，干擾段的結束時間點為xe(m),m=M。特殊情況下，當只有ds大于θ時，干擾段對應的時間序列為x s+1(m)，m= 1,2…M。

設ds是干擾段的開始，d e為干擾段的結束，根據上述推斷，干擾段為x(p)，p=sφ+ 1，sφ+ 2…(e- 1)φ+M，對x(p)進行CEEMD 分解為多個IMF 與一個殘余分量，將殘余分量加入原始信號，重構IMF。

圖3 干擾段抑制示意圖Fig.3 The interference suppression schematic

3 實驗分析

圖4為實測地磁異常圖與隨機生成的運動軌跡，該軌跡地磁圖對應的讀圖值時域及頻域特征如圖5所示，地磁測量周期為0.25 s。在低速行駛下，軌跡對應地磁圖讀值大多為超低頻信號，與磁干擾信號頻譜重疊。

圖4 地磁圖與隨機軌跡Fig.4 Geomagnetic maps and random trajectories

為了驗證ID-CEEMD 對常規隨機噪聲的濾波效果，向軌跡地磁讀圖序列中添加均值為0，標準差為30 nT 的高斯白噪聲，與CEEMD、小波濾波、FIR 和中值濾波進行特性對比，其中，CEEMD 疊加標準差為10 nT 的正負高斯白噪聲，正負噪聲疊加10 次，分解后得到11 個IMF 和殘余分量，小波濾波采用‘db3’5級基波，FIR 采用帶寬為0.001 Hz 的低通濾波器，中值濾波窗口為180 個采樣點，實驗結果如圖6。表1列出了針對常規隨機噪聲，幾種傳統的濾波方法的濾波效果，通過信噪比（SNR）、均方根誤差（RMSE）與最大誤差作為評價參數。數據表示多種方法都可以對白噪聲進行抑制，其中，中值濾波對白噪聲抑制效果最好，ID-CEEMD 對白噪聲也可以進行抑制，抑制效果優于傳統的CEEMD。

圖5 特征圖Fig.5 Characteristic map

圖6 常規隨機噪聲算法結果Fig.6 Experiment results of conventional random noise

圖7 隨機高幅值磁脈沖噪聲特征圖Fig.7 Random high amplitude magnetic pulse interference

表1 常規隨機噪聲幾種算法的濾波效果Tab.1 Filtering effects of several algorithms for conventional random noise

為了驗證ID-CEEMD 對磁通量門傳感器高振幅脈沖干擾的濾波效果，實驗選取地磁臺某時段實測值作為干擾，其時域與頻域特征如圖7所示。從磁干擾頻域圖中可以看出，磁干擾從長時段看來是類似脈沖干擾，由于脈沖持續時間長，其頻率屬于低頻信號。

對測量信號進行CEEMD 分解，其IMF 頻譜如圖8所示，其中1～6 IMF 頻率集中于10-2～1（Hz），7～11 IMF 和殘余分量頻率集中于10-2（Hz）以下。

圖8 地磁測量信號IMF 頻譜圖Fig.8 IMF spectrum of geomagnetic measurement signal

結合地磁圖讀值信號特征，從圖9中IMF 有效指數來看，7～11 IMF 包含主要有效信息，自適應有效IMF 分離準確地將7～11 IMF 分為有效分量，驗證了該方法的可行性。

圖9 IMF 有效指數Fig.9 IMF effective index

實驗對比了FIR、中值濾波、CEEMD、小波濾波與ID-CEEMD，濾波效果如圖10。從圖10 中可以看出，小波濾波無法制止該類信號干擾。經典的CEEMD 和FIR 可以濾除脈沖干擾，但濾波信號失真嚴重，無法表現出地磁異常，中值濾波可以有效地濾除短時脈沖干擾，但當干擾持續時間較長時，濾波信號出現平值，無法完全濾除高幅值干擾。表2列出了測量信號和多種算法濾波結果數據，通過SNR、RMSE 和最大誤差的對比可以看出，ID-CEEMD 將信號信噪比提高了3 倍，最大誤差從151 nT 降低到23 nT，并保持了大部分地磁異常信息。

圖10 算法結果對比圖Fig.10 Comparison of experimental results

表2 隨機高幅值磁脈沖噪聲幾種算法的濾波效果Tab.2 Filtering effects of several algorithms for random high amplitude magnetic pulse interference

CEEMD 對高頻噪聲有著較好的抑制，但當信號同時出現高頻和低頻干擾時，CEEMD 不能較好地抑制。為了驗證ID-CEEMD 對高低頻噪聲同時存在時的抑制效果，將隨機高幅值磁脈沖噪聲與標準差為30 nT的高斯白噪聲疊加作為噪聲，算法對比結果如圖11?？梢钥闯?，ID-CEEMD 能有效抑制高低頻同時存在的噪聲。

圖11 算法結果對比圖Fig.11 Comparison of experimental results

表3列出了不同濾波方法統計數據，ID-CEEMD能有效提高信噪比，同時，最大誤差達到最小。

表3 高低頻噪聲幾種算法的濾波效果Tab.3 Filtering effects of several algorithms for random high and low interference

4 總 結

地磁測量信號的精度直接影響到地磁導航定位精度。本文重點研究了環境中高幅值磁脈沖干擾特性，分析了多種濾波方法的優點與局限，提出了一種多尺度干擾檢測與CEEMD 濾波方法。根據地磁信號強相關特性，計算各IMF 的有效指數，并提出一種自適應的有效IMF 分類方法；提出基于排列熵的干擾檢測方法提取各有效IMF 的干擾時段，并對干擾段進行二次濾波，獲取低頻信息。將濾波信號SNR、RMSE 以及最大誤差作為評價濾波有效性的指標，實驗對比了多種典型的時域及頻域濾波算法的效果，結果表明ID-CEEMD 能有效濾除短時脈沖干擾，對高、低頻干擾有較好的抑制作用，并保持了有效的地磁異常信息。所提出的基于多尺度干擾檢測與完備集成經驗模式分解相結合的降噪方法，運算復雜度適中，實時性較好，能滿足中低動態運載器的地磁導航功能需求；隨著處理器技術的發展，運算實時性將進一步提升，有望滿足高速、高動態運載器的地磁導航需求。