王立輝，劉慶雅，許寧徽

（東南大學(xué) 儀器科學(xué)與工程學(xué)院 微慣性儀表與先進(jìn)導(dǎo)航技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室，南京 210096）

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

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

為了有效濾除常規(guī)隨機(jī)噪聲與高幅值磁脈沖干擾，本文首先分析了地磁信號(hào)特性與干擾信號(hào)特性，利用地磁信號(hào)的強(qiáng)相關(guān)性，計(jì)算各IMF的自相關(guān)指數(shù)，并自適應(yīng)地選取有效IMF。針對(duì)高幅值脈沖磁干擾，提出了一種多尺度的干擾檢測(cè)方法，引入排列熵計(jì)算有效IMF 的時(shí)間復(fù)雜度，并對(duì)復(fù)雜度高的部分進(jìn)行多尺度分析，檢測(cè)干擾成分，并濾除干擾部分后進(jìn)行信號(hào)重構(gòu)，獲得高信噪比的地磁信號(hào)。

1 地磁信號(hào)測(cè)量及特征分析

地磁場(chǎng)一般在幾nT 至幾萬(wàn)nT 之間，屬于弱磁場(chǎng)，一般選用磁通門、質(zhì)子或光泵磁力儀進(jìn)行測(cè)量。磁通門磁力儀能測(cè)量磁場(chǎng)矢量信息，具有動(dòng)態(tài)性能好、采樣頻率高、成本低、體積小和低功耗等優(yōu)點(diǎn)，在導(dǎo)航系統(tǒng)中應(yīng)用廣泛[14]。磁通門傳感器的測(cè)量原理導(dǎo)致其測(cè)量的磁信息并非絕對(duì)值，并且溫度效應(yīng)差，測(cè)量精度低。一般通過(guò)誤差補(bǔ)償與干擾濾波，提高磁通門傳感器的測(cè)量性能。

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

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

圖1 某地磁臺(tái)日測(cè)量值Fig.1 Daily measurement of a geomagnetic station

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

由于磁干擾頻譜特征與地磁圖測(cè)量信號(hào)頻譜接近，傳統(tǒng)的CEEMD 無(wú)法有效地將高幅值長(zhǎng)周期磁干擾的測(cè)量信號(hào)分解為獨(dú)立的IMF。為解決該問(wèn)題，提出了一種基于排列熵的多尺度干擾檢測(cè)方法，結(jié)合CEEMD 的正交性，對(duì)干擾段進(jìn)行二次濾波，取低頻殘余信號(hào)進(jìn)行信號(hào)重構(gòu)。

2.1 CEEMD

經(jīng)驗(yàn)?zāi)B(tài)分解（EMD）是一種針對(duì)非線性信號(hào)的自適應(yīng)濾波方法，它使用三次樣條插值法將信號(hào)分解為具有不同頻率的IMF。為了獲得目標(biāo)信號(hào)，可以通過(guò)分類規(guī)則將IMF 分類為以信號(hào)為主的IMF 和以噪聲為主的IMF，并通過(guò)以信號(hào)為主的IMF 對(duì)目標(biāo)信號(hào)進(jìn)行重構(gòu)。

CEEMD 是在EMD 分解的基礎(chǔ)上，通過(guò)對(duì)原信號(hào)多次疊加正負(fù)對(duì)應(yīng)的白噪聲，補(bǔ)充原始信號(hào)的頻率尺度，減少包絡(luò)線的擬合誤差。具體步驟如下[19]：

1）對(duì)原信號(hào)第i次疊加正負(fù)高斯白噪聲，得到兩組模態(tài)分量：

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

2）分解第i次添加白噪聲后的信號(hào)M1i,M2i，得到對(duì)應(yīng)模態(tài)分量對(duì)應(yīng)余量r1i、r2i。

3）對(duì)各個(gè)分量取平均，得到

4）對(duì)n次疊加白噪聲產(chǎn)生的分量對(duì)應(yīng)相加，求其平均值，最終得到CEEMD 分解結(jié)果：

其中，cj為CEEMD 分解后的第j個(gè)本征模態(tài)分量。

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

2.2 自適應(yīng)有效IMF 分離

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

其中，τ為信號(hào)偏移時(shí)間，T為信號(hào)長(zhǎng)度。在零偏移處，信號(hào)的自相關(guān)函數(shù)達(dá)到最大值。理想高斯白噪聲的歸一化自相關(guān)函數(shù)在零偏移處的值為1，其余偏移為0[20]。與干擾IMF 相比，信息IMF 的自相關(guān)函數(shù)具有更均勻的能量分布。為了比較多個(gè)信號(hào)的自相關(guān)性，需要對(duì)每個(gè)信號(hào)的自相關(guān)函數(shù)進(jìn)行歸一化。歸一化自相關(guān)函數(shù)為

歸一化自相關(guān)函數(shù)避免了信號(hào)能量對(duì)標(biāo)準(zhǔn)偏差的影響，并比較了信號(hào)之間的自相關(guān)。在本文中，將每個(gè)模態(tài)分量的歸一化的歸零自相關(guān)函數(shù)標(biāo)準(zhǔn)偏差用作IMF 的評(píng)估參數(shù)，可以計(jì)算如下：

其中，P(i)為第i個(gè)IMF 的標(biāo)準(zhǔn)差，ρi(τ)為第i個(gè)IMF的歸一化自相關(guān)函數(shù)，為ρ i(τ)的平均值，T為信號(hào)持續(xù)時(shí)間，閾值計(jì)算如下：

當(dāng)P(i)>μ時(shí)，對(duì)應(yīng)的IMF 為信號(hào)IMF，當(dāng)P(i)≤μ時(shí)，對(duì)應(yīng)的IMF 為噪聲IMF。

2.3 多尺度干擾檢測(cè)與濾波

測(cè)量信號(hào)被CEEMD 分解為多個(gè)IMF 后，自適應(yīng)篩選方法對(duì)IMF 進(jìn)行有效性分類，得到多個(gè)有效IMF，此時(shí)單個(gè)IMF 既包含基本測(cè)量信號(hào)又包含磁干擾值。對(duì)有效IMF 進(jìn)行干擾檢測(cè)，并對(duì)干擾段進(jìn)行濾波，可以有效濾除磁干擾。

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

圖2 排列熵滑動(dòng)窗示意圖Fig.2 The sliding windows schematic

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

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

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

根據(jù)Z-score 標(biāo)準(zhǔn)化，如下所示：

其中，μ是di的平均值，σ為di的標(biāo)準(zhǔn)差。

設(shè)置置換熵閾值為θ，將信號(hào)分為干擾段和正常段。如圖1所示，ds是干擾段的開(kāi)始，對(duì)應(yīng)到時(shí)間序列，干擾段的開(kāi)始為x s+1(m)，m= 1，de為干擾段的結(jié)束，對(duì)應(yīng)到時(shí)間序列，干擾段的結(jié)束時(shí)間點(diǎn)為xe(m),m=M。特殊情況下，當(dāng)只有ds大于θ時(shí)，干擾段對(duì)應(yīng)的時(shí)間序列為x s+1(m)，m= 1,2…M。

設(shè)ds是干擾段的開(kāi)始，d e為干擾段的結(jié)束，根據(jù)上述推斷，干擾段為x(p)，p=sφ+ 1，sφ+ 2…(e- 1)φ+M，對(duì)x(p)進(jìn)行CEEMD 分解為多個(gè)IMF 與一個(gè)殘余分量，將殘余分量加入原始信號(hào)，重構(gòu)IMF。

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

3 實(shí)驗(yàn)分析

圖4為實(shí)測(cè)地磁異常圖與隨機(jī)生成的運(yùn)動(dòng)軌跡，該軌跡地磁圖對(duì)應(yīng)的讀圖值時(shí)域及頻域特征如圖5所示，地磁測(cè)量周期為0.25 s。在低速行駛下，軌跡對(duì)應(yīng)地磁圖讀值大多為超低頻信號(hào)，與磁干擾信號(hào)頻譜重疊。

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

為了驗(yàn)證ID-CEEMD 對(duì)常規(guī)隨機(jī)噪聲的濾波效果，向軌跡地磁讀圖序列中添加均值為0，標(biāo)準(zhǔn)差為30 nT 的高斯白噪聲，與CEEMD、小波濾波、FIR 和中值濾波進(jìn)行特性對(duì)比，其中，CEEMD 疊加標(biāo)準(zhǔn)差為10 nT 的正負(fù)高斯白噪聲，正負(fù)噪聲疊加10 次，分解后得到11 個(gè)IMF 和殘余分量，小波濾波采用‘db3’5級(jí)基波，F(xiàn)IR 采用帶寬為0.001 Hz 的低通濾波器，中值濾波窗口為180 個(gè)采樣點(diǎn)，實(shí)驗(yàn)結(jié)果如圖6。表1列出了針對(duì)常規(guī)隨機(jī)噪聲，幾種傳統(tǒng)的濾波方法的濾波效果，通過(guò)信噪比（SNR）、均方根誤差（RMSE）與最大誤差作為評(píng)價(jià)參數(shù)。數(shù)據(jù)表示多種方法都可以對(duì)白噪聲進(jìn)行抑制，其中，中值濾波對(duì)白噪聲抑制效果最好，ID-CEEMD 對(duì)白噪聲也可以進(jìn)行抑制，抑制效果優(yōu)于傳統(tǒng)的CEEMD。

圖5 特征圖Fig.5 Characteristic map

圖6 常規(guī)隨機(jī)噪聲算法結(jié)果Fig.6 Experiment results of conventional random noise

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

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

為了驗(yàn)證ID-CEEMD 對(duì)磁通量門傳感器高振幅脈沖干擾的濾波效果，實(shí)驗(yàn)選取地磁臺(tái)某時(shí)段實(shí)測(cè)值作為干擾，其時(shí)域與頻域特征如圖7所示。從磁干擾頻域圖中可以看出，磁干擾從長(zhǎng)時(shí)段看來(lái)是類似脈沖干擾，由于脈沖持續(xù)時(shí)間長(zhǎng)，其頻率屬于低頻信號(hào)。

對(duì)測(cè)量信號(hào)進(jìn)行CEEMD 分解，其IMF 頻譜如圖8所示，其中1～6 IMF 頻率集中于10-2～1（Hz），7～11 IMF 和殘余分量頻率集中于10-2（Hz）以下。

圖8 地磁測(cè)量信號(hào)IMF 頻譜圖Fig.8 IMF spectrum of geomagnetic measurement signal

結(jié)合地磁圖讀值信號(hào)特征，從圖9中IMF 有效指數(shù)來(lái)看，7～11 IMF 包含主要有效信息，自適應(yīng)有效IMF 分離準(zhǔn)確地將7～11 IMF 分為有效分量，驗(yàn)證了該方法的可行性。

圖9 IMF 有效指數(shù)Fig.9 IMF effective index

實(shí)驗(yàn)對(duì)比了FIR、中值濾波、CEEMD、小波濾波與ID-CEEMD，濾波效果如圖10。從圖10 中可以看出，小波濾波無(wú)法制止該類信號(hào)干擾。經(jīng)典的CEEMD 和FIR 可以濾除脈沖干擾，但濾波信號(hào)失真嚴(yán)重，無(wú)法表現(xiàn)出地磁異常，中值濾波可以有效地濾除短時(shí)脈沖干擾，但當(dāng)干擾持續(xù)時(shí)間較長(zhǎng)時(shí)，濾波信號(hào)出現(xiàn)平值，無(wú)法完全濾除高幅值干擾。表2列出了測(cè)量信號(hào)和多種算法濾波結(jié)果數(shù)據(jù)，通過(guò)SNR、RMSE 和最大誤差的對(duì)比可以看出，ID-CEEMD 將信號(hào)信噪比提高了3 倍，最大誤差從151 nT 降低到23 nT，并保持了大部分地磁異常信息。

圖10 算法結(jié)果對(duì)比圖Fig.10 Comparison of experimental results

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

CEEMD 對(duì)高頻噪聲有著較好的抑制，但當(dāng)信號(hào)同時(shí)出現(xiàn)高頻和低頻干擾時(shí)，CEEMD 不能較好地抑制。為了驗(yàn)證ID-CEEMD 對(duì)高低頻噪聲同時(shí)存在時(shí)的抑制效果，將隨機(jī)高幅值磁脈沖噪聲與標(biāo)準(zhǔn)差為30 nT的高斯白噪聲疊加作為噪聲，算法對(duì)比結(jié)果如圖11。可以看出，ID-CEEMD 能有效抑制高低頻同時(shí)存在的噪聲。

圖11 算法結(jié)果對(duì)比圖Fig.11 Comparison of experimental results

表3列出了不同濾波方法統(tǒng)計(jì)數(shù)據(jù)，ID-CEEMD能有效提高信噪比，同時(shí)，最大誤差達(dá)到最小。

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

4 總 結(jié)

地磁測(cè)量信號(hào)的精度直接影響到地磁導(dǎo)航定位精度。本文重點(diǎn)研究了環(huán)境中高幅值磁脈沖干擾特性，分析了多種濾波方法的優(yōu)點(diǎn)與局限，提出了一種多尺度干擾檢測(cè)與CEEMD 濾波方法。根據(jù)地磁信號(hào)強(qiáng)相關(guān)特性，計(jì)算各IMF 的有效指數(shù)，并提出一種自適應(yīng)的有效IMF 分類方法；提出基于排列熵的干擾檢測(cè)方法提取各有效IMF 的干擾時(shí)段，并對(duì)干擾段進(jìn)行二次濾波，獲取低頻信息。將濾波信號(hào)SNR、RMSE 以及最大誤差作為評(píng)價(jià)濾波有效性的指標(biāo)，實(shí)驗(yàn)對(duì)比了多種典型的時(shí)域及頻域?yàn)V波算法的效果，結(jié)果表明ID-CEEMD 能有效濾除短時(shí)脈沖干擾，對(duì)高、低頻干擾有較好的抑制作用，并保持了有效的地磁異常信息。所提出的基于多尺度干擾檢測(cè)與完備集成經(jīng)驗(yàn)?zāi)Ｊ椒纸庀嘟Y(jié)合的降噪方法，運(yùn)算復(fù)雜度適中，實(shí)時(shí)性較好，能滿足中低動(dòng)態(tài)運(yùn)載器的地磁導(dǎo)航功能需求；隨著處理器技術(shù)的發(fā)展，運(yùn)算實(shí)時(shí)性將進(jìn)一步提升，有望滿足高速、高動(dòng)態(tài)運(yùn)載器的地磁導(dǎo)航需求。