地磁背景下基于传感器阵列的磁偶极子目标跟踪方法

1 引言

随着浅层陆地资源的日益匮乏与开发难度的增大，世界各国对海洋资源的开发与海域权益的争夺日趋激烈。面对部分邻国蓄意挑起的海洋争端及来自境外军事侦查、非法入侵和资源掠夺的敌方潜艇及无人潜航器等对我国的临海主权，尤其是对港口、码头及军事设施等的安全构成严重威胁[1,2]。磁探测技术是一种被动探测技术，因其能够有效弥补声学近海、浅海探测的不足受到国内外的广泛关注[3]。

磁性目标的特性识别主要包括优化定位与实时滤波跟踪两大类。其中，优化定位算法通常依赖大量观测数据，比如传感器网络、磁梯度系统或航空测线数据等，通过最优化的方法获取目标位置、磁矩等参量。优化定位算法研究的典型代表有，以色列的Sheinker等人[4,5]在2005～2008年提出依次基于遗传算法、模拟退火算法及种群增量学习算法实现磁偶极子目标的定位。2010年Yang等人[6]将粒子群与(Levenberg-Marquardt, LM)优化算法结合实现了3个磁偶极子目标的定位与跟踪，2018年Gao等人[7]将LM算法用于运动磁性目标的定位。此外，利用磁梯度张量系统对UXO目标的定位与识别UXO同样吸引了大量研究人员的关注[8–10]，其通常将搭载在移动平台上，如无人机等。在目标的实时滤波跟踪方面，通常利用静止平台的磁传感器阵列测量目标磁异常，实现对目标体的实时跟踪。滤波跟踪算法研究的典型代表有，Birsan等人[11,12]在2003～2006年利用扩展/无迹卡尔曼滤波方法和无迹粒子滤波算法等迭代贝叶斯算法对磁偶极子目标进行实时跟踪；2008年，Kozick等人[13]提出磁传感器阵列目标跟踪的矢量与标量模型，并进行仿真实验研究；2012年日本学者Alhmiedat等人[14]提出将无线磁传感器网络用于实现军事目标的识别与跟踪应用。国内研究学者在近几年针对磁性目标的跟踪已开展大量工作，主要的研究成果有：2010年，张朝阳等人[15]对磁性目标跟踪中磁偶极子模型适用性进行研究；2012年，于振涛等人[16]利用海底磁力仪阵列对水中磁目标定位与跟踪方法进行研究；2013年，吴志东等人[17,18]提出用于实现磁偶极子目标的非线性跟踪算法；2014年，张宏欣等人[19]提出改进的粗糙化粒子滤波的算法改善磁偶极子目标跟踪中的粒子贫化；2015年，周穗华等人[20]提出混合卡尔曼滤波算法解决经典粒子滤波在磁偶极子跟踪中高维模型估计时维数灾问题；2017年吴垣甫等人[21]提出一种递推更新的卡尔曼滤波方法对磁偶极子目标进行跟踪，有效抑制磁偶极子跟踪中初始误差过大导致的滤波发散问题；张宏欣等人[22]提出渐近贝叶斯估计算法和多模型自适应滤波算法对磁偶极子目标进行跟踪，有效抑制目标跟踪中初始误差导致的性能下降和发散等问题；在2018年张宏欣[23]提出多模型自适应滤波算法对磁偶极子目标进行跟踪，算法可在目标源及位置等先验信息完全未知条件下准确估计出真实目标状态，提升真实目标跟踪中的适应力。

分析当前国内外研究现状，磁目标跟踪的主要研究对象在于非线性跟踪模型的优化与不同体制跟踪系统的研究，其中很少考虑到复杂跟踪环境对目标跟踪性能的影响，比如海水中目标在运动过程中存在的涡流磁场、尾流磁场[18]，以及目标跟踪过程中存在地磁背景干扰等，都会对磁目标的跟踪性能产生影响。本文主要针对磁目标跟踪过程中存在的地磁背景干扰的问题开展相关的研究，并在传统磁跟踪模型的基础上提出利用差量磁场对磁偶极子进行跟踪的模型，并考虑到差量磁场模型的高度非线性度，提出基于蒙特卡洛卡尔曼滤波(Monte Carlo Kalman Filter, MCKF)跟踪方法提升跟踪性能。

2 磁偶极子运动模型

铁磁性目标周围产生磁场叠加原背景磁场导致周围的磁场发生畸变，称为磁异常。当目标探测距离较近时，磁性目标可等效为若干个偶极子源叠加；当探测距离为目标尺寸3倍以上时，可将目标视为磁偶极子，其矢量磁场如式(1)，式中 µ0为自由空间磁导率， m 为目标磁矩矢量，r 为距离矢量，其幅度为r。

确立直角坐标系，式(1)可以写为如式(2)形式。式中(x0,y0,z0)代表了在3维坐标系中矢量磁力仪的坐标，(x,y,z)代表坐标系中磁偶极子的坐标；(Bx,By,Bz)代表了磁场矢量在各坐标轴上的投影分量，而(mx,my,mz)代表了磁矩矢量在3维坐标系统的投影分量。

由式(2)可见，任意磁偶极子目标产生的磁场仅依赖于目标与传感器的相对位置及磁偶极子的磁矩。因此，通过对目标位置与磁矩参量的估计，即可完成对运动目标的实时跟踪。

磁偶极子目标跟踪过程中状态变量由两部分组成：在t 时刻，描述目标状态的位置参数[x(t),y(t),z(t)]和速度参数[vx(t),vy(t),vz(t)]以及目标的磁矩参数[mx(t),my(t),mz(t)]。因此，目标状态向量x t如式(3)。

考虑到水面舰船或航行器等磁性目标运动状态相对稳定，目标运动状态可以用恒速模型和恒磁矩模型描述，因此目标运动状态方程如式(4)。

根据目标的位置参量与磁矩参量，通过式(2)，可以获得不同位置传感器的磁异常，当使用3个矢量传感器阵列时，各传感器测量场表示为因此，t时刻的观测场z t可表示如式(5)。

综上，偶极子目标跟踪过程中的状态方程式(4)为线性的；而量测方程式(2)、式(5)为非线性方程。

3 基于差量磁场的MCKF跟踪算法

3.1 差量磁场信号

传感器测量磁场中同时包含目标磁异常与背景地磁场。通常，采用远参考传感器来消除背景地磁场的影响，但基于参考传感器的方法存在两个明显缺陷：其一，参考传感器与传感器阵列的距离很难把握，距离过近则容易受目标磁异常影响；距离过远，则由于相关性的下降致使背景场消除不彻底。其二，参考传感器设置处需要可靠的背景环境，要保证工作中不受外界磁干扰影响。本文针对上述缺陷，提出传感器阵列测量磁场差量实现磁偶极子目标跟踪，在消除背景地磁场干扰的同时，又避免参考传感器的设置难题。

通过仿真示例，阐述磁偶极子差量磁场与测量场及梯度磁场的特性。仿真磁偶极子运动轨迹如图1所示，其轨迹中类似波浪式的运动用于模拟真实目标运动过程中的随机性，“*”代表测量磁传感器，仿真磁偶极子磁矩为[10, 20, 30]A·m2。

传感器测量磁场如图2所示，其中图2(a)代表传感器1测量的磁场3分量；图2(b)代表传感器1处的磁梯度场中3个典型量，依次是x分量场在x方向梯度、y分量场在y方向梯度及z分量场在z方向的梯度；图2(c)代表传感器1, 2的差量场。

从图2可知，差量场与传感器原始测量磁场在磁异常幅度值上相当，而梯度磁场最小。因此，与磁梯度场跟踪相比，采用差量磁场进行目标跟踪时可以获得更远的探测距离。

3.2 基于MCKF的跟踪算法

常用的非线性滤波跟踪算法有扩展卡尔曼滤波(Extended Kalman Filter, EKF)和无迹卡尔曼滤波(Unscented Kalman Filter, UKF)等。EKF算法利用泰勒级数展开，对非线性的系统方程作1阶线性化近似处理，仅适用于弱非线性系统。UKF以确定性的方式选择一组近似高斯分布的样本点，通过非线性方程变换，获得均值和方差的估计，是最常用的非线性卡尔曼滤波算法。本文提出的蒙特卡洛卡尔曼滤波方法，与UKF原理相似，区别在于通过蒙特卡洛的方法生成随机sigma样本点，通过非线性方程变换，获取非线性系统均值与方差估计，完成卡尔曼滤波增益计算与状态的更新，算法流程如图3所示。

图1 仿真磁偶极子运动轨迹与传感器分布位置

图2 传感器测量结果

图3 基于MCKF跟踪算法的流程图

基于图3的算法流程，目标状态变量x及 其状态方程与观测方程分别用f(x)和h(x)表示，则磁偶极子的蒙特卡洛卡尔曼滤波跟踪算法的步骤可归结如下：

(1) 滤波器初始化：定义初始状态变量为9×1维列向量，记为x0；初始化协方差矩阵为9×9矩阵，标记为P0xx。

(2) 状态变量预测：根据蒙特卡洛规则随机生成N 个样本点， 记x(i)n−1|n−1 =xn−1|n−1+Dn−1|n−1c(i),i=1,2,···,N， 样本点权重记为wi =1/N,i=1,2,···,N。根据随机样本点获取状态量均值预测记为式(6)，其协方差矩阵预测如式(7)，式中c (i)～ N(0,I),i=1,2,···,N ; Dn−1|n−1 =Q为变量预测协方差矩阵。

(3) 观测值及其协方差矩阵估计：根据蒙特卡洛规则生成样本点i=1,2,···,N，通过量测方程估计观测均值如式(8)。利用输出样本点及输入样本点获取观测值协方差及观测量与状态变量互协方差矩阵如式(9)，式中R为观测噪声协方差矩阵。

(4) 增益计算与状态变量更新：利用第(3)步得到的协方差矩阵求得滤波的增益估计，如式(10)所示。利用卡尔曼增益与实时观测值，更新状态变量及其协方差矩阵，如式(11)和式(12)。

(5) 存储与迭代：将该时刻状态变量及其协方差存储至跟踪文件，返回第(2)步进行下一时刻迭代。

4 实验与算法验证

4.1 仿真实验

仿真实验如3.1节中描述，磁偶极子目标运动轨迹如图1。用于实现目标跟踪的传感器阵列由3个矢量磁力仪构成，坐标分别为(–2 m, 0, 0), (0, 0,0), (2 m, 0, 0)，依次标记为传感器1、传感器2、传感器3。通过传感器1、传感器2和传感器2、传感器3分别获得差量磁场，并在原始差量场上叠加高斯白噪声模拟工作的噪声。利用观测磁场，分别采用EKF, UKF与本文提出的基于MCKF跟踪算法，实现目标实时跟踪。为评估不同算法的跟踪性能，采用100次随机实验，将跟踪误差的统计平均值作为评判指标。不同仿真算法误差结果如图4所示，其中图4(a)为x轴投影轨迹的误差，图4(b)为y轴投影轨迹的误差，图4(c)为z轴投影轨迹的误差。

从图4可以看出，EKF与UKF具有类似的跟踪效果，MCKF的初期跟踪效果较差，而后期跟踪效果明显优于EKF/UKF的跟踪结果。定性来讲，MCKF在收敛速度上较EKF/UKF更慢，但稳定后在跟踪效果上跟踪精度最高。为定量分析不同算法跟踪性能，计算3种跟踪方法跟踪误差结果如表1。

表1中列举不同时间区间，x, y, z等3个方向上投影轨迹跟踪误差。从表1中看出MCKF算法与EKF/UKF算法相比，在1～40点内跟踪误差较大，随后在40～80点内跟踪误差相当，而在稳定的80～200点内，MCKF的跟踪精度明显优于其余两种跟踪算法；表1结果表明，在z方向3种算法跟踪误差相当。

依据MCKF算法原理，蒙特卡洛样本点数会影响算法跟踪性能。为评估不同蒙特卡洛样本点数的影响，依次采用50, 100, 200, 400样本点进行仿真，结果如表2。从表2中得出随粒子数由50增至200，算法的跟踪误差逐步减小，说明增加粒子数可有效提升算法跟踪性能；而当粒子数由200增至400时，跟踪性能的提升并不明显。由于粒子数多少决定计算量，因此目标跟踪过程中应在满足性能时尽量减少粒子数。

4.2 地磁背景跟踪实验

地磁背景跟踪实验在张家口张北县某地草原进行，选用Mag-03磁通门传感器，传感器分布与磁性目标运动轨迹示意如图5，其中模拟跟踪磁性目标体为通电螺旋圈如图6，磁矩约为60 A·m2。根据图5(b)中运动轨迹，磁性目标在运动过程中的拐点位置由于磁性源的姿态发生变化，此时目标分量磁矩会产生突变。

传感器1、传感器2、传感器3测量的磁异常如图7(a)、图7(b)、图7(c)所示，传感器1、传感器2的差量磁场如图7(d)，传感器2、传感器3的差量磁场如图7(e)。差量磁异常信号与原始磁异常信号在幅度上相近，甚至大于原始测量磁异常。利用本节中提出的基于MCKF跟踪方法与传统的基于EKF/UKF的滤波跟踪方法对目标进行跟踪，目标位置跟踪结果如图8所示，其中图8(a)、图8(b)、图8(c)分别代表x, y, z的跟踪结果，黑色实线为参考轨迹；目标磁矩的跟踪结果如图9所示，其中图9(a)、图9(b)、图9(c)分别代表磁矩x, y, z分量磁矩的跟踪结果，黑色实线为参考磁矩。

由图8和图9的跟踪结果可以看出：(1)在跟踪初始阶段及拐点位置，EKF与UKF算法跟踪响应速度要优于MCKF跟踪方法，与仿真结果一致；(2)当目标运动状态稳定后(图中3500点后)跟踪过程中，MCKF的跟踪性能明显优于其余两种算法。其中，EKF与UKF算法在4300点后，由于信噪比下降算法的跟踪轨迹发生明显偏离，而MCKF直到4800点后轨迹才逐渐发生偏移。综上所述，真实目标的跟踪结果与仿真结果表现一致，说明本文提出的算法具有较好的抗地磁干扰性能，在低信噪比时跟踪性能优于传统滤波方法。

图4 仿真磁偶极子跟踪结果对比

表1 不同跟踪算法各方向投影轨迹跟踪误差(m)

图5 真实目标跟踪实验

表2 不同蒙特卡洛样本点数的算法各方向投影轨迹跟踪误差(m)

图6 跟踪实验模拟磁性目标体

5 结束语

本文分析了地磁场背景下磁偶极子目标跟踪过程中存在背景干扰问题，提出利用传感器差量磁场信号实现目标跟踪的模型。新模型的优势有：(1)传感器阵列测量磁场差量与原始测量场相比磁场强度相似，因此可以保证足够远的跟踪距离；(2)测量磁场差量不易受到背景地磁场的干扰。同时，考虑差量磁场模型的高度非线性，本文提出基于MCKF的滤波跟踪算法，从而有效解决减小跟踪模型的非线性误差。通过仿真实验，结果表明基于新模型的跟踪算法较经典的UKF、EKF算法具有更高的跟踪精度；地磁背景跟踪实验，同样证实算法的有效性。

图7 阵列传感器测量磁场与差量磁场

图8 模拟磁性目标的位置跟踪结果

图9 模拟磁性目标的磁矩跟踪结果