1 引言

磁梯度张量较传统磁探测量具有更为丰富的信息，更加有利于对磁性目标的定位和跟踪，为磁性目标探测提供了有效手段，相关理论也已成为国内外的研究热点(Clark, 2012, 2013；Schneider et al., 2013；于振涛等，2014a；张光等, 2015, 2016；Lee and Li, 2016；尹刚等，2016).在磁性目标探测领域中，磁偶极子是一种重要的磁源等效模型(张朝阳等，2010).因此，可以利用磁梯度张量测量值对等效为磁偶极子的目标进行在线定位.

海军舰艇研究与发展实验室(Naval Ship Research and Development Laboratory, NSRDL)的W. M. Wynn设计了包含五个三轴磁力计的十字形磁测系统，可以测量单点磁感应强度和磁梯度张量，在此基础上对C. P. Frahm提出的磁偶极子张量定位算法进行了实验验证(Wynn et al., 1975)；此后又证明了只有在已知目标磁场三分量或测量系统与目标的相对速度时，才可以唯一地确定目标真实位置(Wynn et al., 1975). Nara等(2006)以定位射频标签为应用背景，同样采用十字形张量测量系统，提出了一种更为简洁的通过单点目标磁梯度张量和目标磁场磁感应强度三分量定位磁偶极子的方法.尹刚等(2016)利用磁梯度张量中几何不变量的关系，结合目标磁感应强度三分量解算出了目标位置.

前述方法需要测量目标磁场磁感应强度三分量，而在多数实际应用场合中，地磁背景场的干扰使目标磁场磁感应强度三分量难以测量.而地磁背景场的梯度远远小于目标磁场的梯度，因此目标磁感应强度梯度及其组成的张量矩阵几乎不受地磁背景场的影响，基于这些特征量研究人员提出了多种不受地磁背景场干扰的方法.

李光等(2012)、于振涛等(2014b)、Yin等(2014)在Nara方法的基础上提出了基于磁感应强度一阶、二阶差分的磁偶极子定位方法.但由于二阶差分量本身较小，受仪器测量噪声影响较大，导致了这两种方法对仪器精度要求较高.Liu和Wang(2010)以反恐行动中探测路边炸弹为背景，设计了三轴磁力计测量阵列，通过测量多点磁梯度张量，建立了六维目标参数向量与磁梯度张量的关系方程，在用粒子群算法得到相关参数粗估计的基础上，利用最速下降法得到更精确的目标参数估计.但是由于需要估计的目标参数较多，且初始解随机，导致定位平均时间为2 s，无法应用在对实时性要求较高的场合.

美国海军水面作战中心的R. F. Wiegert等对张量约缩量的性质进行了研究，提出了利用张量约缩量的STAR(Scalar Triangulation And Ranging)定位方法(Wiegert et al., 2007, 2008)，采用了含有8个三轴磁力计的立方体张量测量阵列，通过同时测量6个点的磁梯度张量完成定位.但由于该方法建立在将张量约缩量近似为一个等值面为球面的不变量的基础上，所以存在较大系统误差.针对STAR方法存在的“非球面误差”，随阳轶等(Sui et al., 2012)通过对原始STAR方法的结果进行迭代处理，减小了定位误差.吕俊伟等(2015)则提出了新的等值面为球面的不变量，从根本上消除了“非球面误差”.经过分析与比较，吕俊伟改进的STAR方法是目前在相同的地磁背景场和仪器测量噪声的条件下性能最好的方法(刘继昊等，2016)，但是仍然存在张量测量系统复杂、需要较多传感器等问题.

为进一步提升地磁背景场和仪器测量噪声存在的情况下在线定位的精度，本文在Frahm方法解算张量方程组(Wynn et al., 1975)的基础上，结合Levenberg-Marquardt优化算法提出了一种通过六个三轴磁力计测量两点磁梯度张量进而对磁偶极子进行在线定位的方法.

2 Frahm方法原理

定义规范化磁矩为(Scaled Magnetic Moment)：

(1)

其中M为磁偶极子的磁矩矢量，r为张量测量点与磁偶极子的距离，m₁、m₂、m₃分别表示m在坐标系x、y、z轴上的投影.

方向向量为(Bearing Vector)：

(2)

其中r为张量测量点与磁偶极子的相对位置矢量，n₁、n₂、n₃分别表示n在坐标系x、y、z轴上的投影.

则张量矩阵的各个元素可以表示为(Wynn，1997)：

(3)

其中：

(4)

为Kronecker函数，μ_o为磁导率，G_ij表示G的第i行、第j列的分量.

上述由张量分量构成的9个方程中，有5个方程之间是独立的，再结合方向向量的模为1，可得到方程组(5) 为：

(5)

在得到张量测量值后，方程组(5) 中剩余的未知量显然可以解出，但是解法复杂.为了方便解算，利用特征值分解的方法，寻找一个旋转矩阵，使得旋转后的坐标系中，张量矩阵为原始张量矩阵的特征值矩阵，即：

(6)

其中：

(7)

它为特征值矩阵，λ₁、λ₂、λ₃为特征值，R为对应的特征向量矩阵，每一列为特征向量，R^-1即坐标系转换对应的旋转矩阵.在新的坐标系下，方程组化简为：

(8)

其中n₁′、n₂′、n₃′、m₁′、m₂′、m₃′是m和n在新的坐标系中对应坐标轴上的投影.

解方程组(8) 可以得到四组解，将得到的解旋转回原始坐标系为：

(9)

即得到原始坐标系下的四组解.之后，补充测量目标在测量点处的磁感应强度，便可从中筛选出真实的m和n并解算距离，完成定位.

3 两点磁梯度张量定位方法 3.1 方法原理

Frahm方法是获得相对位置向量的方向向量的有效方法，但为了从四组解中筛选出真实的m和n并解算距离r，必须测量目标磁感应强度，然而实际过程中磁力计所能测量的是目标磁感应强度与地磁场磁感应强度的叠加值，难以获得纯净的目标磁感应强度.而地磁场在小范围内可近似为均匀磁场，其差分量远远小于目标磁感应强度的差分量，目标磁感应强度差分量可近似认为不受地磁场干扰.因此，本文提出了仅依赖差分量的真实解筛选方法和r解算方法.

由Frahm方法可知，单凭一个点的磁梯度张量所蕴含的信息是不够的，需要补充新的观测量，本文补充观测量是另一点的磁梯度张量.设两个观测点分别为A点和B点，每个点对应的变量用下标A、B区分，例如G_A和G_B分别代表A点和B点的磁梯度张量.

(1) 真实解的筛选

真实解得筛选主要依据两个事实：首先，实际的应用场景中目标位于探测点的下方，据此对于每个点可以剔除两个包含虚假目标信息的解.其次，由于两个探测点观测的是同一个目标，所以真实解中两个点的规范化磁矩应该指向同一方向，即对应的单位向量应该相同，公式为：

(10)

而实际过程中由于干扰的存在，式(10) 不可能成立，所以笔者规定使式(10) 左侧绝对值达到最小的为真实解.

基于上述原则，真实解的筛选流程如图 1所示.

(2) 距离r_A的解算

为实现距离的计算，首先可以利用的信息是A、B两点以及目标之间的相对位置向量的空间关系，设A、B两点之间的相对位置向量为d，则有：

(11)

其次，还可以利用张量不变量与距离之间的比例关系.吕俊伟等(2015)提出的张量约缩量的改进量为：

(12)

本文对其进一步推导得到更为简洁的形式为：

(13)

对于两个点有：

(14)

变形得：

(15)

联立方程(11) 和(15) 可得方程组：

(16)

方程组(16) 是一个超定线性方程组，利用SVD方法可获得该方程组的最小二乘解.但是即使是最小二乘解，由于线性方程组对干扰的敏感性，此种方法解算的r_A只有在信噪比很高时才准确.

为了提高该方法的抗干扰能力，可补充利用磁力计测得的磁感应强度之差与磁偶极子参数的正演磁感应强度之差相等的条件，即：

(17)

其中：

(18)

下标i、j为磁力计标号，B_i代表第i号磁力计的测量值，Δr_i表示第i号传感器与探测点的相对位置向量，由测量阵列的设计方案确定，为已知量.实际应用中，由于干扰的存在，对r_A的解算转化为最优化问题，公式为:

(19)

以线性方程组解得的r_A作为初始解，结合Levenberg-Marquardt算法即可解出更为精确的r_A.

3.2 两点张量测量阵列

根据前文所述的原理，本着尽量减小磁力计个数、装配简单的原则，本文设计了两点张量测量阵列，如图 2所示，由6个坐标轴方向一致的三轴磁力计组成.

两点张量测量阵列中，磁力计1、2、3、4组成一个正方形张量测量阵列，完成G_A的测量，由于目标磁场的散度和旋度均为0，其磁梯度张量矩阵为对称矩阵且主对角线元素之和为0，因此只需要测量其磁感应强度分量在X、Y方向上的梯度即可得到磁梯度张量(Wynn et al., 1975)，具体计算方法为:

(20)

其中:

(21)

它表示第i号磁力计与第j号磁力计的磁场k分量均值，近似等于空间上两个磁力计连线中点处的磁场值.

磁力计3、4、5、6组成另一个正方形张量测量阵列，完成G_B的测量，3、4号磁力计的重复利用达到了减少系统磁力计个数的目的.

3.3 方法步骤

结合方法原理以及所设计的张量测量阵列，论文提出的具体算法为：

(1) 读取磁力计1到磁力计6的磁感应强度测量结果B_i, (i=1, 2, 3, 4, 5, 6)；

(2) 计算得到G_A、G_B；

(3) 基于方程组(8) 和图 1中方法，解算并筛选出真实的m_A与n_A，m_B与n_B；

(4) 解方程组(16) 得到r_A；

(5) 利用步骤(4) 得到r_A作为初始解，利用Levenberg-Marquardt算法进行上限为N次的迭代，求得方程(22) 的最优解，作为最终的r_A，方程(22) 为：

(22)

(6) 完成定位

(23)

4 仿真实验与结果分析 4.1 仿真实验方案与实验条件

取张量测量系统的坐标系为参考坐标系，为方便比较，每次定位时其张量测量系统姿态不变.以磁偶极子中心点为坐标系原点，定位点在(50, －50, 10)、(50, 50, 10)、(－50, 50, 10)、(－50, －50, 10) 确定的平面内通过网格法取得，网格大小为1 m×1 m，在每个点上对磁偶极子进行一次定位.

定义定位误差e的计算公式为：

(24)

其中为算法估计的位置向量，r为真实位置向量r的模值，即探测点与目标的距离.

考虑到实际应用，定位误差大于10%时认为定位失败，在下文的定位误差分布图中用纯白色网格表示定位失败的探测点，同时，定义定位成功的探测点所占所有探测点的比例为定位成功率，用以衡量算法的性能.

仿真实验过程中条件的选择过程如下：

(1) 针对不同的目标，同一方法的定位效果可能不同，所以目标的磁矩M要作为条件之一.

(2) 实际定位过程中，磁偶极子磁场往往不是单独存在的，而是与背景磁场叠加在一起.探测点所测得的磁感应强度表达式为：

(25)

其中B_m表示测得的磁感应强度，B_t表示目标的磁感应强度，N_B表示仪器测量噪声，B_e表示地磁背景场的磁感应强度.所以，仪器测量噪声与地磁背景场磁感应强度也应作为条件.

(3) 基线长度d是表示张量测量过程中，计算磁感应强度梯度的两个磁力计之间的距离，会对张量测量误差产生影响，故也应作为条件.

4.2 定位实验及对比分析

设定目标磁偶极子的磁矩M=[20000, 10000, 30000]^T(A·m²)，地磁背景磁场B_e=(49000, 100, 0)(nT)，本文方法Levenberg-Marquardt算法迭代次数上限设为10次，为更好对比，参考文献(吕俊伟等，2015)中的仿真条件，仪器测量噪声设定为均值为0 nT，标准差为0.01 nT的高斯噪声，张量测量系统的基线长度取为d=1 m.对本文方法和吕俊伟等(2015)改进的STAR方法同时进行仿真实验.实验结果如图 4所示.具体统计结果如表 1所示.

对比分析图 4a、b可知，本文方法的定位成功的范围要大于吕俊伟改进的STAR方法.从表 1可知，本文方法的定位成功率高，可定位范围大.但需要注意的是优化算法的采用导致了本文方法定位速度较慢，而直接以空间分布较为复杂的磁梯度张量分量解算定位信息导致定位成功区域内误差波动较大.

4.3 仪器测量噪声对定位效果的影响

保持4.2节中的其他条件不变，仪器测量噪声依次取均值为0 nT，标准差为10^-6、10^-5、10^-4 nT、10^-3、10^-2、10^-1、10⁰、10¹ nT的高斯噪声，对本文方法进行仿真实验，实验结果如图 5所示.

由图 5可知，随着噪声强度的增加，本文方法定位成功区域的半径逐渐缩小，这是因为随着半径的增大，距离目标越远，目标磁感应强度越小，越容易受到噪声的影响.分析图 5f的变化曲线可得，对于本仿真实验方案，当信号噪声小于10^-4 nT时，定位结果几乎不受噪声影响，而当噪声强度大于10^-4 nT时，定位成功率随着噪声强度的增加迅速下降.

4.4 优化算法迭代次数上限对定位效果的影响

保持4.2节中其他条件不变，Levenberg-Marquardt算法迭代次数上限依次取0、5、10、15、20、25、30、35次.对本文方法进行仿真实验，实验结果如图 6所示.

分析图 6可知，在Levenberg-Marquardt算法迭代次数上限小于10次时，定位的成功率随着迭代次数上限的增加而增加，且定位成功区域向外扩张，同时，定位成功范围内的定位误差空间分布特征也发生了改变，原来只采用线性方法(迭代次数为0) 时，可定位区域内包含的一个明显的白色定位失败区域，在使用优化算法修正后消失，以上现象均说明方法最后的优化算法确实提高了本文方法的精度.

然而，当迭代次数上限从10次再增加时，定位成功区域无明显变化，定位成功率也稳定在64%左右.这说明，迭代次数上限对方法性能的提升是有极限的，同时也证明了仅仅需要10次迭代，便能达到利用优化方法修正线性方法结果的目标，这保证了本文方法的实时性.

理论上分析，该种现象产生的原因是利用Levenberg-Marquardt算法解算距离只是整个定位方法中的一个步骤，其结果所能达到的精度受到之前步骤结果的制约，如方向向量解算精度，张量计算精度等.因此定位成功率并不能够随着迭代次数的增加一直增加到100%.

4.5 基线长度d对定位效果的影响

保持4.2节中其他条件不变，张量测量系统的基线长度d依次取0.1、1、2、3、4、5、6、7 m.对本文方法进行仿真实验，实验结果如图 7所示.

分析图 7可知，随着基线长度的增加，定位成功区域向外拓展，但同时内部的定位失败的“裂缝”也开始扩大，最终导致定位成功率随着基线长度的增加先增加后减小，在d=3 m时达到最大值.经过分析，笔者认为，定位成功区域向外拓展是因为随着基线长度的增加，非线性方法中求差的两个点之间距离越大，两点磁感应强度之差越大，而仪器测量噪声不变，所以仪器测量噪声对定位的影响减小，进而可定位区域的边界向外扩张；而“裂缝”的扩大是因为，磁梯度张量的测量原理是利用差分结果近似微分，随着基线长度的增加，差分与微分之间误差增大，导致张量测量误差增大，进而本文方法前半部分解算得到的方向向量误差增大，最终导致定位成功区域内“裂缝”扩大.

5 结论

仿真实验表明，本文方法利用6个三轴磁力计实现了对目标的定位，且不受地磁背景磁场的影响，在相同的条件下，仿真实验定位成功率优于吕俊伟等2015年提出的STAR改进方法.分析仿真实验结果，可以得出以下结论：

(1) 仪器测量噪声强度达到一定值后，定位成功率会随着噪声的增大迅速下降，降低仪器测量噪声是改善定位效果的有效途径.

(2) 兼顾定位效果和速度，本文方法中优化算法的迭代次数上限最好选择10次.

(3) 其他条件不变时，随着基线长度的增加，本文方法定位成功率先升高后降低，因此实际应用中应选择合适的基线长度.