基于水下运动平台的磁性目标跟踪与特征反演

马剑飞， 丁凯， 颜冰， 林春生

(1.海军工程大学 兵器工程学院， 湖北 武汉 430033；2.近地面探测技术重点实验室， 江苏 无锡 214035)

0 引言

相比于直接测量目标的三分量磁场，磁梯度张量信息包含5个独立变量，具备不受地磁场波动的影响以及可以更加全面地反映磁场细节等突出的优势。在浅海环境下，声纳探测易受水文条件以及地形等因素的影响，磁梯度张量被动探测相对而言具有抗干扰能力强、定位精度高以及隐蔽性好的优势[1]。在军事应用中，特别是在港口地区，可以有效探测敌方的潜艇、水雷目标或者未爆弹(UXO)。在民用领域，对于沉船打捞或者海洋矿物勘探等应用都有一定的应用价值[2-3]。

目前，基于磁梯度张量的探测技术逐步成为磁探测领域的研究热点，德国、美国、澳大利亚等国家已开发出磁梯度张量探测系统，并开展了对应的航空航海试验[4-6]。相比于航空磁探，基于水下运动平台的磁探测能够更加接近磁性目标，因此其测量的信噪比较高，特别是对于微弱的磁性目标检测而言意义尤为重要，而且利用有限的几个测量点能够很快反演出磁性目标的位置，进一步可以反演出目标的形状以区分目标种类。水下运动平台受风浪的影响较小，比在海面上运行更加稳定[1]，因此水下运动平台的运动噪声也会比较低。

现有的磁梯度定位方法大多都是将磁性目标视为偶极子模型[7-9]，在远场情况下这种简化模型是可靠的，但在近场情况下这种忽视目标尺度的模型导致定位结果出现较大的误差。另外，远场的磁梯度信号较微弱，定位结果受测量噪声影响较大，而近场情况下信噪比较高，目标尺度因素(即模型不匹配)则是定位误差的主要来源。Clark[10]与Beiki等[11]研究了常见几何体的构造指数、几何项以及其位置的反演方法，但是一般情况下目标的构造指数和几何项是未知的，因此很难应用于实际。

在磁性目标模型选择方面，吴志东等[12]研究了基于磁性椭球体模型的磁梯度张量水面舰船目标跟踪问题，利用提出的高斯混合采样粒子滤波算法取得了较好的跟踪结果，在拟合目标磁梯度时假定椭球体的焦距是已知的，然而焦距参数是一个未知的参数，即使给定一个初值，由于焦距观测性比较低[13]，滤波器更新过程中也很难趋近于真值，因此利用该模型去反演目标尺度特征是不合适的。戴忠华等[14]利用偶极子模型与椭球体的混合模型实现对舰船磁场模拟，已知目标与测量节点的距离，利用大量节点的磁场测量数据去反演混合模型的参数，实现了与舰船磁场很好的拟合，得到了相对较好的反演效果。张宏欣[15]以线列阵偶极子模型模拟舰船的磁场，并利用水下双节点测量的三轴磁场对实现了舰船目标的准确定位跟踪。事实上，利用多个偶极子完全可以对尺度目标的磁场特性进行精确逼近，磁偶极子的分布方式也可以有效反演目标的尺度特征。

在磁性目标分类的研究中，葛健等[16]研究了不同姿态航弹对应的垂直磁梯度空间分布特征，证明了利用磁梯度特征识别掩埋磁性目标的可行性。张敬东[17]利用小波多尺度分解数据处理方法有效提高了弱磁梯度异常信息的分辨率。郑建拥等[18]以磁总量模值、磁梯度5个独立变量和3个特征值不变量作为特征，利用改进的支持向量机方法提升了目标分类的准确率。Zhang等[19]基于多面体磁模型采用最陡下降法对目标的形状进行了反演，证明了目标磁梯度能够弥补磁场信息的不足，提升反演的精度。然而，上述两种方式只局限于在特定距离条件下对特定目标的分类，在水下目标分类应用中不具有实际意义。

磁性目标跟踪滤波器的观测信号一般为磁矢量信号，目标磁矩、尺寸以及速度等信息隐含于观测信号之中，而且滤波跟踪算法先验信息缺失，设定的初值范围不可避免地引入了人为经验，大量的未知参数也导致其跟踪精度有限。相比较而言，以磁梯度张量作为估值滤波器观测信号的优势在于：1)磁梯度张量拥有5个独立变量，可以提高观测矢量的维数；2)磁梯度张量相对于磁矢量可以更好地抑制环境噪声的干扰；3)磁梯度定位算法能够为跟踪滤波器提供较为准确的滤波器磁矩初值以及位置初值。值得注意的是，随着估值滤波器观测值的不断引入，偶极子的空间分布特性也逐渐显现，其分布特性能够在一定程度上反映磁性目标的尺度和结构特征。

针对磁性尺度目标模型不匹配问题，本文旨在结合磁梯度张量与滤波估值算法的优势，以多磁偶极子为模型研究磁性尺度目标的跟踪方法，并根据多偶极子的空间分布状态研究目标尺度和结构特征的反演算法。

1 滤波估值算法

目前以目标磁场作为信号源对运动目标进行跟踪的研究已经十分广泛，其估值滤波算法主要分为粒子滤波框架与卡尔曼滤波框架，然而这两类滤波器的估值精度对滤波器初值比较敏感，不能保证任意初值条件下的滤波收敛性[15]，但利用磁梯度定位算法可以为滤波器提供较准确的距离、速度以及磁矩初值，从而提升磁性目标跟踪的准确性。

鲁棒估计是指在出现异常量测情况下，利用适当方法对目标状态进行估计，减小异常量测值对状态估计带来的误差。当出现观测野值时，传统卡尔曼滤波算法在对状态进行估计时会产生较大的误差。为了对目标状态进行稳健的估计，人们将鲁棒统计学应用于滤波算法中，提出了基于Huber方法的线性卡尔曼滤波算法[20]，实质是一种广义极大似然估计的卡尔曼滤波算法。文献[21]在线性鲁棒卡尔曼滤波的基础上，揭示了Huber方法在滤波算法中对卡尔曼增益和状态估计的影响，并基于此推导基于无迹卡尔曼滤波(UKF)器的非线性鲁棒滤波算法，以避免线性化带来的额外误差。

考虑如下运动目标：

xk+1=f(xk)+wk,

(1)

其量测方程为

zk+1=h(xk+1)+vk+1,

(2)

式中：xk为k时刻的运动状态矢量；f为状态转移函数；wk为k时刻的运动噪声；zk+1为k+1时刻的观测矢量；h为观测函数；vk+1为k+1时刻的测量矢量。

构造如下非线性回归方程：

(3)

定义如下变量：

式中：Rk+1为观测噪声的协方差矩阵。则回归变量

yk+1=g(xk+1)+ξk+1.

(4)

据此可定义Huber代价函数定义为

(5)

式中：ek+1,i为yk+1-g(xk+1)的第i个分量；m为回归变量的维数；ρ函数为

(6)

u为门限值，e为估计误差。

对ρ函数关于e求导，可得φ函数为

(7)

由(6)式可知，当误差|e|

对(5)式进一步求导，可得

(8)

(9)

(10)

(11)

由于真实的目标状态xk+1未知，令残差δxk+1|k=0，则Ψx为单位矩阵I，即

(12)

可得Huber方法改进后的UKF算法等效量测噪声为

(13)

由上述推导过程可知，非线性Huber无迹卡尔曼滤波(NHUKF)算法通过最小化目标代价函数来抑制异常观测值对目标状态估计带来的影响，而且没有对非线性方程的线性化过程，从而避免了传统算法的线性化误差[20]。

2 建立跟踪估值模型

图1所示为一个水下运动平台磁梯度跟踪定位的典型示意图，其中平台的运动参数可以通过平台的惯性导航准确获取，因此平台的速度参数不需要作为状态估计量。考虑到需要通过偶极子的分布特征来判断目标形状，在三维空间中最少需要4个点才能表达立体结构特征，而过多的偶极子会导致状态估计量的维数激增从而大幅增加运算量，故选用4偶极子模型作为NHUKF算法的状态估计模型。

图1 载体平台运动示意图Fig.1 Schematic diagram of moving carrier platform

目标的磁性来源主要是铁磁性材料，而且同一物体经过的外磁场磁化过程也大致相同，同时也是出于降低观测量的维数和利用磁偶极子分布特征反映物体形状特征的考虑，假设4个偶极子的磁矩是相同的。综上所述，设定滤波器状态估计量为

(14)

由于磁梯度张量仅包含5个独立变量，状态变量是15维的，而观测量维数过低导致卡尔曼增益产生的状态量更新方向不一定是真实解的方向，容易陷入歧义解的范围。在水下运动平台磁定位的应用中，平台速度是已知的，可以利用前面时刻的测量磁梯度扩充观测量的维数，相当于变相地进行阵列跟踪定位，因此设置观测量为

zk=[Gk-2n,Gk-n,Gk],

(15)

式中：Gk、Gk-n、Gk-2n分别表示k时刻、k-n时刻和k-2n时刻5个独立磁梯度变量组成的向量。不难看出，此时观测量和状态量都是15维的向量，可以进行准确地磁梯度跟踪定位。

3 目标尺度反演方法

NHUKF估值滤波算法的偶极子位置估计值是4个离散的三维空间点坐标，对于三维离散点要提取其分布特征，直观的想法是通过线性拟合的方式提取其主要的分布尺度与方向。然而使用多元线性回归在三维空间拟合的结果是一条曲线，很难反映目标的主尺度特征，而且现有拟合方式大多通过计算预测值和观测值的均方误差来解算拟合函数，对于求解主尺度特征而言，更希望以各点到直线的距离为代价函数。

主成分分析(PCA)法是一种通过正交变换将一组可能存在线性相关性的变量转换为一组线性不相关变量的统计方法，常用于特征降维以更好地发掘数据模式。下面结合PCA法的计算过程和目标主尺度分析的需要进行简单说明。假设偶极子的坐标分别为(xi,yi,zi)，首先计算这些偶极子坐标的协方差矩阵，即

(16)

对C进行矩阵正交分解，可以得到3个特征值λ1>λ2>λ3，分别对应3个正交的单位特征向量u1、u2、u3.显然，特征值越大表示这些离散点在该特征向量方向投影的方差就越大，最大特征值对应的特征向量为第1主成分。对于目标尺度分析的应用而言，第1主成分表示目标的最大尺度在此方向上通过各离散点在第1主成分上的投影就可以估计出目标的主尺度，各离散点在第1主成分方向的投影为

pi=[xi,yi,zi]u1.

(17)

常见的掩埋目标形状主要为圆柱体以及圆盘两种典型结构，这里约定当圆柱体的高H大于底面半径时称之为圆柱，当圆柱体的高H小于底面半径时称之为圆盘。对于圆柱体而言，第一主成分方向与圆柱体的轴线一致，其主尺度对应圆柱体的高，偶极子分布接近于一条直线，此时λ1应该远大于λ2和λ3，因此目标的主尺度可估计为

(18)

式中：N为状态估计模型的偶极子数。

对于圆盘而言，λ1和λ2的大小相差不多，且远大于λ3，其第1主尺度和第2主尺度反映的是圆盘的直径特征。鉴于本文利用4个偶极子模拟目标的磁场特征，在分析主尺度和圆盘直径的代数关系时采用4个小圆来模拟大圆的磁场。根据等面积法，小圆盘的半径r等于大圆盘半径的一半，显然，若小圆覆盖的区域与大圆覆盖的区域重合部分越多，则拟合的效果越好。图2所示为两种典型的覆盖策略。其中第1种策略是小圆之间相交且与大圆相切，此时所估计的小圆半径为

图2 偶极子拟合策略Fig.2 Fitting strategy of multiple dipoles

r1=L/2,

(19)

式中：L为4个小圆圆心组成的正方形的对角线长度。第2种策略是小圆之间相切且与大圆相交，此时所估计的小圆半径为

(20)

对于正方形顶点组成的离散点集合，利用PCA法可以推导出其协方差矩阵的最大特征值λ1与其对角线长度L的关系为

(21)

综合以上两种策略，可以推导得到圆盘目标的主尺度为

(22)

4 仿真实验

为了验证算法的可行性，开展利用搭载磁梯度张量系统的运动载体平台探测静止磁性目标的仿真实验，实验内容主要分为两部分：第一部分检验NHUKF算法的定位性能；第二部分检验PCA法在目标尺度反演时的性能。

4.1 基于NHUKF算法的目标定位

设磁性目标的中心为空间直角坐标系的原点，x轴与水下运动平台的轴线重合并指向自主水下机器人(AUV)平台的艏部，y轴指向右舷，z轴竖直向下。仿真条件设置为：初始时刻目标的位置为(-40 m， 5 m， -5 m)；磁梯度张量测量系统基线Ls为0.5 m；采样频率fs为5 Hz；AUV平台运动速度为(5 m/s，0 m/s，0 m/s)；磁梯度系统的测量精度为0.01 nT. 如图3所示，目标是高为3 m、半径为0.5 m的圆柱体，目标横滚角、俯仰角以及方位角分别为π/6 rad、π/4 rad和π/3 rad，磁梯度正演基于密集的偶极子网格，网格边长为0.1 m，每个磁偶极子的磁矩M为(-3 A·m2， 1 A·m2， 5 A·m2)。跟踪滤波器的初值由3阶张量法确定[8]。

图3 磁性目标正演模型Fig.3 Forward modeling of magnetic target

如图4所示，在无异常观测的情况下，NHUKF算法与UKF算法的定位性能相近，随着观测值的不断引入，其定位性能要优于3阶张量法，定位误差小于0.1 m. 在接近目标时，多数时刻3阶张量法都有较好的定位结果，但有的时刻其定位矩阵条件数很大(与磁偶极子的特征平面有关)，稍微的扰动或者观测精度不足就会导致定位结果发生较大的偏差。如图5所示，异常观测的引入在距目标-20 m、-10 m以及0 m的时刻，不难看出NHUKF算法相对于UKF算法而言，能够有效消除异常观测对定位性能带来的影响。另外可以发现3阶张量法的定位性能受异常观测影响十分大，在距目标-20 m、-10 m以及0 m的时刻，其定位误差都出现了比较大的峰值。

图4 无异常观测时的定位性能比较Fig.4 Comparison of localization performances without abnormal observation

图5 存在异常观测时的定位性能比较Fig.5 Comparison of localization performances with abnormal observation

由上述分析可知，利用磁梯度张量的NHUKF算法够有效定位尺度目标，其定位性能要优于传统磁梯度张量定位算法，而且能够有效消除异常观测的影响。

4.2 基于PCA算法的目标尺度反演

对于掩埋或者半掩埋的磁性目标而言，大多数目标都具有规则的几何形状，通过反演的偶极子分布特征、磁矩特征和目标尺度特征，在很大程度上对目标的种类进行区分，也可以排除环境的干扰。本节基本仿真条件与4.1节一致，分别检验PCA法在不同磁矩、不同姿态、不同尺度以及不同形状情况下对目标尺度的反演性能。

4.2.1 不同磁矩的目标尺度反演

表1所示为不同目标磁矩[Mx,My,Mz]对应的反演结果，其中Mx、My和Mz分别表示单个偶极子在x轴、y轴和z轴方向的磁矩，主尺度方向表示偶极子分布协方差矩阵的最大特征值λ1对应的单位特征向量。从表1中可以看出，反演主尺度方向误差在1°的范围内，反演主尺度模值的误差小于0.3 m，正演偶极子的磁矩方向不会影响算法反演的性能，而且算法可以准确估计出目标的整体磁矩特征。

表1 不同磁矩的目标反演

图6所示为不同偶极子磁矩对应的反演结果。由图6可以看出，反演的偶极子位置虽然略有差别，但基本都是沿圆柱体的轴线分布，反演的主尺度参数能够很好地反映出其轴线的尺度和方向的特征。在滤波器状更新过程中，有限的几个偶极子通过不断调整位置更好地拟合目标磁场特征，显然当偶极子沿主尺度分布时能够更好地减小滤波器的预测值和观测值的误差。

图6 不同磁矩的目标反演Fig.6 Inversion of targets with different magnetic moments

4.2.2 不同姿态的目标尺度反演

表2所示为在不同目标姿态条件下对目标尺度的反演结果，其中[θ,β,φ]中θ为横滚角、β为俯仰角、φ为方位角。从表2中可以看出，不同姿态目标的反演结果都比较准确，反演的主尺度方向误差在2°范围内，反演的主尺度模值误差小于0.3 m.

表2 不同姿态的目标反演

图7所示为不同姿态的目标反演结果，反演的偶极子沿圆柱体轴线分布，反演的主尺度也对应于圆柱体轴线，反演的姿态角也与目标真实姿态角十分接近，表明本文的主尺度估计算法不受姿态角影响，而且该算法可以准确估计目标的姿态。

图7 不同姿态的目标反演Fig.7 Inversion of targets with different attitudes

4.2.3 不同尺度的目标尺度反演

表3所示为不同目标尺度对应的反演结果，其中R、H分别表示圆柱的底面半径和高。从表3中可以看出，不同尺度目标的反演结果都比较准确，反演的主尺度方向误差在2°的范围内，反演的主尺度模值误差小于0.3 m.

表3 不同尺度的目标反演

图8所示为不同尺度目标的反演结果，其中反演的偶极子沿圆柱体轴线分布，在不同目标尺度情况下，反演的主尺度模值和方向与圆柱体轴线的模值和方向相近，表明反演算法可以准确估计不同大小磁性圆柱体的尺度特征。

图8 不同尺度的目标反演Fig.8 Inversion of targets with different scales

表4 不同的圆柱目标反演

4.2.4 不同形状的目标尺度反演

表4对比了实心以及空心圆柱体的反演结果，表5对比了实心以及空心圆盘的反演结果，其中实心表示无空腔的目标；外壳表示空心目标的外部壳体；空腔表示目标内部的中空部分。由表4和表5可以看出，空腔对目标反演的结果影响不大，算法仍然可以准确估计目标主尺度的方向以及模值，圆柱体偶极子分布协方差矩阵的特征值满足λ1≫λ2>λ3的关系，圆盘偶极子分布协方差矩阵的特征值满足λ1≈λ2≥λ3的关系，以此特征值关系可以区分目标为圆柱或者圆盘，联合目标的磁矩特征以及尺度特征就可以推断目标的类别。

表5 不同的圆盘目标反演

如图9(a)、图9(b)、图9(c)、图9(d)所示，圆柱状结构目标的反演偶极子基本沿轴线分布，圆盘状结构目标的反演偶极子沿中心切面对称分布。对于圆柱状结构目标，第1主尺度反映的是其轴线的尺度；对于圆盘状结构目标，第1主尺度和第2主尺度反映的其横截面的尺度。

图9 不同形状的目标反演Fig.9 Inversion of targets with different structures

5 结论

本文利用多个偶极子组合模型，结合具有鲁棒性的NHUKF算法研究一种针对尺度磁性目标的磁梯度定位方式，该方法能够有效提高磁性尺度目标的定位精度，特别是当出现异常观测时仍然可以准确定位目标，同时也可以准确反演偶极子的分布特征和磁矩特征，具有良好的鲁棒性。进一步提出一种基于偶极子分布模式反演目标尺度特征的PCA算法，该方法对不同磁矩、不同姿态、不同尺度以及不同形状的目标，主尺度反演误差小于0.3 m，对目标主尺度方向的反演误差小于2°，亦可识别目标主尺度为线结构或者面结构，为水下目标种类判别提供了一定的参考。