Research on Several Nonlinear Filtering Methods for Vehicle Magnetic Field Tracking*

YU Zhenhua，ZHOU Qianwei，XIE Dongfeng，LI Baoqing*

(Shanghai Institute of Microsystem and Information Technology，CAS，Shanghai 200050，China)

Research on Several Nonlinear Filtering Methods for Vehicle Magnetic Field Tracking*

YU Zhenhua1，ZHOU Qianwei1，XIE Dongfeng1，LI Baoqing1*

(Shanghai Institute of Microsystem and Information Technology，CAS，Shanghai 200050，China)

This article first summarizes the estimation and tracking method of vehicles magnetic field，and illustrates the advantages and existing problems to use non-linear filtering algorithms.Second，the principles and models of two methods，the“Two sensors”method and the"Gradiometer"method，in which nonlinear filtering algorithm is applied，are described.What’s more，a“Gradient tensor”method combined the Euler inversion with Kalman algorithm is proposed.The Cramer-Rao Bounds of nonlinear filtering algorithm is used to analyze and compare the performance of these methods mentioned above.At last，simulation results are carefully demonstrated and show that all these methods have the ability to estimate and track a vehicle.However，the“Gradient tensor”method performs much better than other methods in terms of the estimation of magnetic moment.

vehicle magnetic field tracking;nonlinear filtering;gradient tensor;ramer-Rao bounds

磁异常探测MAD(Magnetic Anomaly Detection)技术是当今潜艇探测的主要手段。潜艇探测领域中MAD装置一般有两种，一种是只测量异常磁场的磁力仪MAD装置;另一种是测量异常磁场及其磁梯度张量的磁梯度张量仪MAD装置［1］。和潜艇类似，车辆也可以等效为磁偶极子模型，将MAD技术运用到地面车辆上，我们就能够对车辆进行磁场探测、估计和跟踪［2－3］。很多UGS(Unattended Ground Sensors，无人值守地面传感器)系统中就使用了MAD技术，如美国海军陆战队战术远程传感系统TRSS(Tactical Remote Sensor System)的AN/GSQ－257UGS子系统中就使用了磁场探测节点MAGID和MAGID－Ⅱ。

车辆的磁场跟踪问题一般指的是估计车辆的位置和磁矩并对轨迹进行跟踪。可以根据估计出来的轨迹和磁矩将车辆的感应磁矩分离出来，从而使车辆分类成为可能［4］。该问题一般有3类方案:第1类是通过多个传感器测量不同点的场强来实时估计磁矩和位置，可以使用参数拟合算法。第2类是通过单个传感器测量场强来估计车辆的方向、速度(或最近点距离)，可以使用Anderson方程和子空间匹配［4］算法、非线性滤波［5］算法，文献［5－7］中说明了这种方法不能直接实现实时跟踪，需要依赖布设方案。第3类是利用复杂磁测仪测量单点的场强和梯度进行实时磁矩和位置估计，可以使用Wynn提出的算法［4］、Euler反演算法［8］，如美国的Quan-tum Magnetics公司生产的低功耗磁阻传感器系统。显然第2类方案与前面所说的磁力仪MAD装置对应，而第3类方案与磁梯度张量仪MAD装置对应。这种对应说明了车辆探测和潜艇探测的相似性。而第1类方法恰恰反映了在地面上使用传感器网络对车辆进行探测的特点和灵活性。

实际上还可以使用非线性滤波算法进行车辆跟踪，近些年才陆续有文献研究。非线性滤波算法能够实时进行轨迹跟踪，无需进行非线性方程求解或拟合运算，计算代价更小。而且由于利用了历史观测数据，提高了跟踪精度。文献［9－10］采用两个传感器方案，利用EKF、IEKF、UKF以及粒子滤波算法［11－12］对磁偶极子进行了轨迹跟踪。文献［13］采用梯度方案对磁偶极子进行了跟踪。但是多数文献都只是将非线性滤波算法用在两个传感器方案，少数针对梯度方案使用非线性滤波算法，基本上还没有文献针对梯度张量方案使用非线性滤波算法。而后两种都是单点测量方案，比起两个传感器等多点测量方案，在实际使用中更简单。同样对于两个传感器方案、梯度方案和梯度张量方案性能上的比较和分析也缺乏研究。但是这些内容涉及到在实际实现中选择哪种方案，比较很重要。本文将针对这两点进行研究。

1 几种方案分析

当磁场传感器与车辆的距离大于车辆自身最大几何尺寸的2.5倍时，可以将车辆的磁场模型等效为磁偶极子模型［6，14］，多数对车辆的磁场探测、估计和跟踪均是建立在此基础上的，这样可以简化计算。磁偶极子模型产生的磁场如下:

其中r=(x，y，z)为场源距离观测点的位置矢量，M=(Mx，My，Mz)为场源的磁矩矢量，B(r)为场源在观测点处产生的磁感应强度矢量，下文中这些定义保持不变。

一般地，可以用动态状态空间模型DSSM(Dynamic State Space Model)来描述磁偶极子跟踪问题。

其中，X为状态变量，A为状态转移矩阵，n为系统噪声，h(X)为观测函数，Y为观测变量，w为观测噪声，下文中的这些定义保持不变。状态变量由位置、速度及磁矩的组成，即X=(x，y，z，vx，vy，vz，Mx，My，Mz)。系统模型中的状态转移矩阵的具体形式由下式决定:

其中v=(vx，vy，vz)分别为场源的速度矢量，下文中定义不变，Ts为采样周期。考虑到车辆的机动性和磁矩的扰动，给速度和磁矩矢量一个随机噪声，这使得系统模型更符合实际情况。对于不同的方案，其系统模型相同，观测变量和观测模型不同，以下介绍几种常见方案。

1.1 “双传感器”方案

根据两个传感器同时刻观测的6个数据，如图1(a)，可以直接通过非线性方程组将位置和磁矩求解出来［5］。但是这种方法没有考虑历史信息，而且该非线性方程组求解比较复杂，难以保证实时性。文献［6］中证明了“双传感器”能够采用非线性滤波方法对目标进行跟踪。文献［6］中同时使用了Fisher矩阵方法证明了当两个传感器放在道路两侧对称位置且间距较大时，Fisher矩阵会更大，即估计误差会更小。两个传感器的位置必须已知为r1，r2，则观测变量为Y=(B1x，B1y，B1z，B2x，B2y，B2z)，其中B1和B2定义如下:

“双传感器”方案(简称方案1)是最简单的多点测量方案，没有硬件上的冗余。但“双传感器”方案要准确知道相对位置，这给布设带来了困难。而且每个传感器都需要知道自身的姿态，才能使得观测数据变换到同一个坐标系中，这要求每个传感器都要附带姿态角获取模块。数据的同步和传输也是个实际要考虑的问题。这是多点测量方案的共同缺点。

图1 磁场检测的两种传感器布设方案

1.2 梯度方案

文献［13，15］采用了一种梯度方案，并使用粒子滤波方法进行了跟踪，该文献结果说明了这种方案比直接Euler反演方法性能更好。其观测变量为:Y= (Bx，By，Bz，Gzx，Gzy，Gzz)，其中B和G定义如下:

其中，Gij的测量使用了两个传感器磁场数据进行一阶差分近似得到，传感器的基准距离(传感器的间距)为d，如图1(b)。此方案的观察噪声为均值为μ(Xk)的高斯噪声，μ(Xk)可以采用误差传递公式计算得到。

梯度计方案(简称方案2)基准距离可以较小，是一种单点测量方案，这显然弥补了“双传感器”方案实际使用中的缺点。

1.3 梯度张量方案

文献［8－9，15］中描述了磁梯度张量方法对磁偶极子的参数进行估计，文献［8］在此基础上使用了Euler反演法［16］(Euler反褶积法)，使得估计更加简单，避免了求解复杂的非线性方程组。磁梯度张量矩阵定义如下:

如图2是文献［17］中提出的优化的梯度张量测量系统方案。其中x轴分量的梯度由1号和5号传感器测量值做差分得到，y轴分量的梯度由2号和4号传感器测量值做差分得到，z轴分量的梯度通过计算得到［17］，3号传感器得到该点的磁感应强度矢量。

图2 梯度张量测量系统

Euler反演公式如下:

这里提出对Euler反演法的结果进行Kalman滤波(简称方案3)。Euler反演法的结果即为带噪声的位置和磁矩，其噪声是由梯度张量测量系统的测量噪声引起的，可以证明该噪声是零均值加性高斯噪声，其方差可以通过误差传递公式计算得到。很容易发现，此时观测模型和系统模型都是线性的，只要使用经典的Kalman滤波就能估计出系统状态。观测变量为:Y=(^rx，^ry，^rz，^Mx，^My，^Mz)，其中的分量由Euler反演公式计算得到。

与直接使用梯度张量来做观测变量相比，方案3使得计算得到大大简化。从3.1节图4中可以看出，方案3比Euler反演法的性能更好。方案3虽然在硬件上稍有冗余，但是可以实现单点测量，弥补了“双传感器”方案的缺点。另外这种冗余提高了估计的性能，使得基准距离可以做得很小。值得一提的是，梯度张量方案中的基准距离并不是随意确定的。基准距离太大，Euler反演公式的条件不能满足;基准距离太小，受到传感器准确度限制，可以通过仿真的方法确定合适的基准距离，这里略去说明。对于高精度的磁场传感器来说，可以使得基准距离做得很小，如JESSY STAR系统［1］中使用了超导量子干涉仪作为传感器，其基准距离只有3.5 cm［1］，可惜这种传感器还无法应用在车辆探测领域。相信随着磁传感器的发展，尤其是MEMS原子磁力仪的进步，在车辆探测领域中实现高精度、微型化的梯度张量方案是可以预期的。

2 Cramer-Rao下限分析

非线性滤波的Cramer-Rao下限决定了非线性滤波算法的性能下限，这里使用该方法来分析上面几种方案的性能，这样可以屏蔽具体非线性滤波算法的使用或者优化上的差异所导致的性能差距。文献［18］推导了非线性滤波Cramer-Rao下限。信息矩阵J的递推公式，如下:

式中变量定义如下，P(Xk+1|Xk)为转移概率密度，P (Yk+1|Xk+1)为似然概率密度。

上述所有方案的系统噪声和观测噪声都是加性的，所以有P(Xk+1|Xk)=P(Xk+1－AXk)，P(Yk+1| Xk)=P(Yk+1－h(Xk))。所有方案的系统噪声为均值为零的加性高斯白噪声，所以:D11k=AQ－1kA，

方案1和方案3中观测噪声是均值为零的加性高斯白噪声，有:

方案2中观察噪声是均值为μ(Xk)的加性高斯白噪声，有:

以下使用Matlab对3种方案的非线性滤波Cramer－Rao下限进行仿真计算。文献［4］中使用的武装车辆ATV的磁矩约为104A·m2～105A·m2数量级，轻型轿车约103A·m2数量级，SUV约104A·m2数量级，所以仿真中车辆的磁偶极子模型等效磁矩在103A·m2～105A·m2区间随机选择。目前磁传感器的准确度很容易达到nT级，如磁通门等。考虑到地磁场的漂变，在环境磁噪声较小的野外地区，假设单个单轴传感器的观测误差在5nT以内是比较合理的。仿真中取单个三轴传感器观测误差的协方差矩阵R如下:

仿真结果如图3所示。当方案1中的传感器位置距离为20 m，方案2的基准距离为0.5m，方案3的基准距离为0.1m时，其仿真结果如图3(a)所示。当方案1的传感器位置距离为0.5 m(此时布设在道路的同一侧)，方案2和方案3的基准距离为0.5 m时，其仿真结果如图3(b)所示。由仿真结果可以看出，方案2的性能较差。当方案1传感器间距较大时，距离和速度的估计上，方案1的性能最好，但随时间增加方案3和方案1趋于一致，而磁矩的估计上，方案3的性能最好。当所有方案的传感器间距均相同时，方案3性能最好，而方案1和方案2性能相差不大。由于方案1和方案2都是两个传感器组成的系统，没有增加新的信息量，所以当传感器间距相同时两者性能基本相同。方案3的性能优势是由于增加了传感器个数提供了更多的信息。实际上，梯度或者梯度张量能够对环境噪声有抑制作用，所以，在同样的条件下，方案2和方案3的性能应该要比仿真结果要稍微好一点。

图33 种方案的非线性滤波的Cramer-Rao下限比较

3 实际算法仿真

本节给出上述方案1和方案3的实际算法的性能仿真图，仿真中采用100次蒙特卡洛仿真得到每个时刻的估计误差的均值和均方根误差，并与Cramer-Rao下限进行比较，验证算法的性能。仿真中的参数设置与第3部分相同。仿真时非线性滤波代码大部分来自ReBEL工具箱。假定初始状态的值均未知，采用“两点法”计算得到初始状态。

3.1 梯度张量方案的算法仿真

根据前面方案3的描述，直接使用经典Kalman滤波就能够对车辆进行跟踪。基准距离选择0.1 m，算法其他参数和第3小节的一致。图4为Euler反演法结合Kalman算法(即方案3)与Euler反演法的比较。

图4Euler反演法与方案3比较(10次蒙特卡罗仿真的二维轨迹)

图5 为方案3在100次蒙特卡罗仿真后得到的性能。在当前参数设置下，使用信噪较高的地点做起始点，位置估计的均方根误差最小能达到0.3 m，速度估计的均方根误差最小能达到0.5 m/s，磁矩估计的均方根误差最小能达到200 A·m2，相对实际磁矩，误差约为0.7%。再次体现了方案3的良好性能。

3.2 “双传感器”方案的算法仿真

这里使用非线性Kalman滤波算法实现方案1，两个传感器距离为20 m，其他参数和第2小节的一致，如图6为EKF、srUKF和srCDKF算法的100次蒙特卡罗得到的性能。在目前参数设置下，使用信噪较高的地点做起始点，位置估计的均方根误差最小能达到0.2 m，速度估计的均方根误差最小能达到0.5 m/s，磁矩估计的均方根误差最小能达到1 000 A·m2，相对实际磁矩，误差约为3%。表1中列举了常见的非线性Kalman滤波算法在“双传感器”方案中的应用性能。可以看出UKF和srUKF算法性能较好。

表1 各种非线性Kalman滤波算法性能(100次蒙特卡罗仿真)

图6 “双传感器”方案的100次蒙特卡罗仿真的性能

4 总结

以上讨论了几种基于非线性滤波算法的磁偶极子跟踪的方案，“双传感器”方案、梯度方案、梯度张量方案，在梯度张量方案中提出来将Euler反演法和Kalman滤波结合进行磁偶极子跟踪。并使用了非线性滤波的Cramer-Rao下限对几种方案的性能进行分析和比较。几种方案均能够对车辆进行跟踪，梯度和梯度张量方案是单点测量方法，能够解决“双传感器”方案等多点测量方案的缺点。梯度张量方法在车辆跟踪问题上性能更好，尤其是磁矩估计更准确，计算代价也更小。

［1］张昌达.关于磁异常探测的若干问题［J］.工程地球物理学报，2008，4(6):549－553.

［2］刘石，李宝清，童官军，等.一种平面四轴向磁通门传感器的设计［J］.传感技术学报，2011，23(11):1565－1569.

［3］唐莉莉，宋勇，童官军，等.一种新的运动磁目标定位算法［J］.传感技术学报，2011，24(7):996－1000.

［4］Merlat L，Naz P.Magnetic Localization and Identification of Vehicles［J］.P Soc Photo-Opt Ins，2003，5090:174－185.

［5］任志良，黄玉盈，陈泽茂.磁偶极子源的被动测距与跟踪［J］.数据采集与处理，2001(3):380－383.

［6］Wahlstrom N，Callmer J，Gustafsson F.Magnetometers for Tracking Metallic Targets［C］//Proceedings of the Information Fusion(FUSION)，2010 13th Conference on，F，2010.IEEE.

［7］唐劲飞，龚沈光.基于磁偶极子模型的目标定位和参数估计［J］.电子学报，2002，30(4):614－616.

［8］Nara T，Suzuki S，Ando S.A Closed-Form Formula for Magnetic Dipole Localization by Measurement of Its Magnetic Field and Spatial Gradients［J］.IEEE Transactions on Magnetics，2006，42 (10):3291－3293.

［9］Czipott P V，Perry A R，Whitecotton B R，et al.Magnetic Detection and Tracking of Military Vehicles［J］.2002.

新型组合式直流断路器拓扑及其性能分析//郑晓铭,王钢,李海锋,郭彦勋,李健涛,洪潮//(24):105

［10］Wahlstrom N，Callmer J，Gustafsson F.Single Target Tracking U-sing Vector Magnetometers［C］//Proceedings of the Acoustics，Speech and Signal Processing(ICASSP)，2011 IEEE International Conference on，F，2011.IEEE.

［11］Birsan M.Unscented Particle Filter for Tracking a Magnetic Dipole Target［C］//Proceedings of the Oceans，2005 Proceedings of MTS/IEEE，F，2005.IEEE.

［12］Wu Z D，Zhou S H，Zhang H X.GMSPPF for Magnetic Ellipsoid Tracking with a Magnetic Gradiometer Tensor［J］.Advanced Materials Research，2013，779:1789－1792.

［13］Birsan M.Recursive Bayesian Method for Magnetic Dipole Tracking with a Tensor Gradiometer［J］.Ieee T Magn，2011，47(2):409－415.

［14］Caruso M J，Withanawasam L S.Vehicle Detection and Compass Applications Using AMR Magnetic Sensors［J］.Sensors Expo Proceedings，1999，477.

［15］Wynn W M.Magnetic Dipole Localization Using the Gradient Rate Tensor Measured by a 5-Axis Magnetic Gradiometer with Known Velocity［J］.Detection Technologies for Mines And Minelike Targets，1995，2496:357－367.

［17］Wynn W M，Frahm C P，Carroll P J，et al.Advanced Superconducting Gradiometer-Magnetometer Arrays And a Novel Signal-Processing Technique［J］.Ieee T Magn，1975，Ma11(2):701－707.

［18］Tichavsky P，Muravchik C H，Nehorai A.Posterior Cramér-Rao Bounds for Discrete-Time Nonlinear Filtering［J］.IEEE Transactions on Signal Processing，1998，46(5):1386－1396.

余振华(1972－)，男，硕士研究生。2011.6毕业于华中科技大学，专业为电子科学与技术专业。2011.9至今，于中科院上海微系统与信息技术研究所无线传感网实验室攻读硕士研究生，主要研究方向为无线传感网中的传感器信号处理，scholar.edward@gmail.com;

李宝清(1972－)，男，研究员。1999.9毕业于中科院上海冶金研究所获博士学位，专业为微电子。2001.5～2005.5美国新泽西理工学院微电子研究中心，访问学者，2006.4回国后在中科院上海微系统与信息技术研究所无线传感网实验室工作，主要研究方向为各类传感器的研制及其在无线传感网中的应用，libq@mail.sim.ac.cn。

车辆磁场跟踪问题的几种非线性滤波方案研究*

余振华1，周乾伟1，谢东峰1，李宝清1*
(中科院上海微系统与信息技术研究所无线传感网实验室，上海200050)

概述了车辆的磁场跟踪方案，并说明了在这些方案中运用非线性滤波算法的优点以及存在的问题。描述使用非线性滤波算法的“双传感器”方案和梯度计方案的原理和模型，并提出了一种Euler反演法与Kalman算法结合的梯度张量方案。使用了非线性滤波的Cramer-Rao下限分析并比较了这几种方案的性能。最后给出了方案的仿真结果。这几种方案均能实现对车辆的跟踪，其中梯度张量方案在磁矩估计上有明显的性能优势。

车辆磁场跟踪;非线性滤波;磁场梯度张量;Cramer-Rao下限

O441.5;TN850.6;TP274;TP212.9

A

1004－1699(2014)04－0505－07

2014－01－20修改日期:2014－04－02

C:7950;7310L

10.3969/j.issn.1004－1699.2014.04.016

项目来源:系统技术国防科技重点实验室基金项目(9140C18010213ZK34001)