基于干扰检测和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 能有效滤除短时脉冲干扰，对高、低频干扰有较好的抑制作用，并保持了有效的地磁异常信息。所提出的基于多尺度干扰检测与完备集成经验模式分解相结合的降噪方法，运算复杂度适中，实时性较好，能满足中低动态运载器的地磁导航功能需求；随着处理器技术的发展，运算实时性将进一步提升，有望满足高速、高动态运载器的地磁导航需求。