磁梯度张量测量及解释应用被视为磁法勘探工作的一次突破(张昌达，2006； 2007)，其在资源勘 探、军事、环境等领域都有着重要的应用价值(Schmidt and Clark，2000； 2006).磁梯度张量系统测量的对象是磁场矢量的梯度场，测量的结果受磁化方向影响小，能够反映目标体的矢量磁矩信息，且能更好地反演场源参数(方位、磁矩等)并对场源进行定位和追踪(张凤旭等，2007)，提高了磁源体的分辨率(马国庆等，2012).

磁梯度张量具有良好的数学性质，但其搭载在运动载体上使用时，载体姿态的变化将影响其测量和数据解释的精度(张光等，2013a)，由张量元素计算得到的张量不变量因其不受姿态变化的影响而得到了较多的关注，被广泛应用于磁性目标的定位及磁异常源的边界识别等多种场合(张光等，2013b).当探测距离大于2.5倍磁性物体长度时，可将磁性目标简化为磁偶极子模型(张朝阳等，2010)，因此，本文在张量不变量的基础上，以磁偶极子产生的磁梯度张量为对象，研究磁梯度张量矩阵、磁偶极子位置及磁矩间的本质关系，推导出了蕴含在其中的、不受姿态变化影响的固定几何关系，本文称之为磁梯度张量的几何不变量.

磁梯度张量是磁场矢量的3个分量在相互正交的3个方向的空间变化率.若 B 为磁场矢量场，则磁梯度张量可以表示为两个三元素矢量的矩阵乘法运算：

由式(1)和(2)可知，磁梯度张量矩阵 G 为对称矩阵，其中9个元素中，只有5个是独立的，即只需得到其中的5个元素就可测得磁梯度张量.在磁法探测中，磁性目标简化为磁偶极子模型时可由6个未知量描述，即磁性目标的三维位置和三维磁矩.因此，磁矩矢量为 m =(m_x，m_y，m_z)的磁偶极子，在距离其 r =(x，y，z)处产生的磁矢量场(卞光浪等，2011)及磁梯度张量矩阵中的5个独立分量(Rim et al.，2012)可分别由下式计算得到

其中，μ₀为真空磁导率，R= |r| 为磁偶极子位置与测量点之间的距离.

易知，磁梯度张量矩阵的特征值与坐标的选择无关，为张量的基本不变量(余天庆和毛为民，2006)，且张量矩阵的其余不变量可由特征值表示.设矩阵 G 的3个特征值为λ₁，λ₂和λ₃，则由 G 的特征方程可得磁梯度张量的另外4个不变量，分别表示为

为完成磁梯度张量的几何不变量的推导，首先给出矩阵 G 的特征值和其对应特征向量的求解公式.对于任意形状的磁性目标产生的磁梯度张量矩阵，令：，则磁梯度张量矩阵的3个特征值可表示为

将式(4)及式(5)代入式(6)可得磁偶极子产生的磁梯度张量矩阵的特征值：

令：δ=ec²+bcd+abc+d²e+ade+e³，ε=a²c+ adc+abe+c³+ce²+bde，Φ=bc²+cde－be²－ace，则磁偶极子产生的梯度张量的三个特征值对应的特征向量 V ₁，V ₂，V ₃分别为

将式(8)、式(9)及式(10)进行归一化处理即可得到3个单位特征向量：₁，₂和₃.

由式(7)中磁偶极子张量矩阵的特征值公式可知，假设3个特征值已由磁梯度张量系统测得的实际信号求解得到，则磁偶极子位置矢量与磁矩矢量间的夹角可由下式求解得到

其中，和分别为位置矢量及磁矩矢量的单位向量.

由式(11)可知，磁梯度张量的特征值为张量不变量，因此，在磁梯度张量系统发生姿态变化时，由该公式计算得到的磁偶极子位置矢量与磁矩矢量间的夹角保持不变.本文将其定义为磁梯度张量的几何不变量1，其几何示意如图 1所示.

图 1

Fig. 1

图 1 磁偶极子梯度张量场的几何不变量示意图Fig. 1 Sketch map of magnetic gradient tensor′s geometric invariants

由磁偶极子产生的磁矢量、磁梯度张量场公式及几何不变量1可得：

因此，基于本文推导的几何不变量1，若已知磁偶极子产生的总磁场强度及磁梯度张量信息，可通过特征值分析及公式(12)求得磁偶极子相对于测量点的距离R和磁矩模M.

将式(4)代入式(10)化简可得如下等式：

由式(13)可知，磁梯度张量矩阵 G 的绝对值最小的特征值对应的特征向量垂直于磁矩矢量和位置矢量，且不随着张量系统姿态的变化而变化.本文将此固定的垂直关系定义为磁梯度张量的几何不变量2.

由对称矩阵的特征向量之间的几何关系可知，特征向量 V ₃垂直于 V ₁和V ₂；由几何不变量2可知，特征向量 V ₃垂直于磁矩矢量 m和位置矢量 r ，因此，假设向量 V ₁，V ₂，m和 r 共面，则4个向量之间可互相表示，令

又磁偶极子产生的梯度张量场与矢量场满足欧拉反褶积公式(Reid et al.，1990；Nara et al.，2006；Nara and Ito，2014)，且构造指数为3，即

由式(18)及式(20)可知，位置矢量及磁矩矢量与最大及最小特征值对应的特征向量共面，且不随着系统的姿态变化而变化，本文将其定义为几何不变量3.因此，几何不变量3可用于磁性目标的定位和磁矩反演.

在已知磁偶极子产生的磁梯度张量矩阵的情况下，可通过特征值分析，利用式(12)求解磁偶极子相对于测量点的距离和磁矩模，然后利用式(18)及式(20)求解磁偶极子的位置及磁矩矢量的单位向量，最终实现磁偶极子位置及磁矩矢量的正确反演，其几何示意如图 2所示.

图 2

Fig. 2

图 2 几何不变量3表示的磁矩及位置矢量示意图Fig. 2 Inverse solution of source-dipole vector and magnetic moment based on geometric invariant 3

由图 2可知，该方法得到真实位置及磁矩矢量的同时，也得到了3个虚假解，因此，在实际应用中，可根据所测目标在磁梯度张量系统的上方或下方，直接去除两个虚假解，然后根据磁梯度张量系统测得的多点总磁场强度及磁矢量场进行另一个虚假解的去除，进而实现真实磁偶极子位置及磁矩矢量的正确反演.

利用本文提出的几何不变量进行磁偶极子位置及磁矩矢量反演的三组仿真实验，以验证所提几何不变量的准确性.其中磁梯度张量系统位于坐标原点，磁偶极子运动轨迹的具体参数设置如表 1所示.

表 1(Table 1)

表 1 仿真参数设置 Table 1 The simulation parameters 序号 初始位置/m 最终位置/m 运动速度/(m·s -1) 采样频率/Hz 磁矩/Am 2 1 (-10,-50,-5) (-10,50,-5) (0,10,0) 20 (2000,2000,1000) 2 (-55,-40,5) (45,60,5) (10,10,0) 20 (2000,2000,1000) 3 (50,-50,3) (-50,0,3) (-10,5,0) 20 (2000,2000,1000)

表 1 仿真参数设置 Table 1 The simulation parameters

仿真实验流程如下：

(1)根据磁偶极子数学模型，计算磁偶极子位于不同位置时在坐标原点产生的磁梯度张量值.

(2)对计算得到的磁梯度张量数据进行特征值分析，求解其对应的特征值和特征向量.

(3)利用公式(12)求解磁偶极子相对于测量点的距离及磁矩模.

(4)将计算得到的特征值及特征向量代入式(19)及式(20)，求解磁偶极子位置及磁矩矢量.

利用本文所提几何不变量反演得到的磁偶极子在x-y平面上的运动轨迹如图 3—5所示.由图可知，几何不变量估计得到磁偶极子真实运动轨迹的同时也求解得到3个虚假解，必须采用有效方法准确去除虚假解.虚假解的去除流程如图 6所示，去除虚假解后的真实轨迹在图 3—5中以深色表示.

图 3

Fig. 3

图 3 估计得到的磁偶极子运动轨迹(仿真实验1)Fig. 3 Estimated track of magnetic dipole coming from the first simulation

图 4

Fig. 4

图 4 估计得到的磁偶极子运动轨迹(仿真实验2)Fig. 4 Estimated track of magnetic dipole coming from the second simulation

图 5

Fig. 5

图 5 估计得到的磁偶极子运动轨迹(仿真实验3)Fig. 5 Estimated track of magnetic dipole coming from the third simulation

图 6

Fig. 6

图 6 虚假解去除流程Fig. 6 Flow chart of elimination of ghost solutions

由于在实际应用中所探测目标在测量平面的上方或下方为已知信息，且单点测量的数据反演得到的四组磁偶极子的三维坐标关于原点对称，因此，根据已知信息中的Z坐标信息可直接删除四组解中的两个虚假解.将剩下的两组解代入式(3)和式(4)中可计算得到两组不同的解在测量点产生的磁梯度张量值.若为真实解，则该计算值与实际测量值相吻合；若为虚假解，则计算得到的张量数据与测量数据

存在较大的偏差，因此，根据计算值与测量值之间的差异可进一步去除第3个虚假解，进而得到真实的磁偶极子位置及磁矩矢量.值得注意的是，虚假解的去除是实时的，即每个测量点反演得到的四组解都需要进行虚假解的去除.

由于在实际应用中，磁梯度张量系统存在测量误差将在一定程度上影响磁偶极子探测和位置反演的距离，因此，将上述理论模型仿真中加入随机噪声进行研究.考虑真实磁梯度张量系统多用来测量5 个独立的磁梯度张量分量及张量矩阵对称的原则，在5个独立的磁梯度张量分量信号中加入均值为0 nT·m^-1，方差为0.5 nT·m^-1的高斯白噪声.图 7为表 1中第一组仿真实验加入噪声后的运动轨迹反演结果.由图可知，加入随机噪声后，磁偶极子的可跟踪范围为y轴方向上的[-20 m，20 m].在可跟踪范围内，本文所提算法可较为准确地跟踪磁偶极子，而超出此跟踪范围后，反演得到的磁偶极子的位置坐标存在较大的偏差和不确定性，反演结果无法为磁偶极子的运动轨迹提供参考.同样地对表 1中的第二及第三组仿真实验加入随机噪声，则第二组仿真实验的可跟踪的坐标范围为[-20，-5]至[8，23]，第三组仿真实验的可跟踪的坐标范围为[0，-25]至[-22，-14].由式(3)及式(4)可知，磁矢量信号及磁梯度张量信号随着距离的增加衰减极为迅速，因此，距离磁偶极子目标较远处的磁场信号较弱，有效信号可能被测量噪声所淹没，进而反演的得到的磁偶极子位置及磁矩矢量存在较大的偏差.即在实际应用中，磁偶极子的跟踪距离不仅仅受到其自身磁矩大小的影响，还受到测量系统中测量噪声的影响.

图 7

Fig. 7

图 7 存在测量噪声时估计得到的磁偶极子运动轨迹Fig. 7 Estimated track of magnetic dipole existing measurement noises

为验证本文所提几何不变量在磁性目标定位及跟踪中的可行性，采用4个磁通门传感器构建了平面十字形磁梯度张量系统进行实验研究.实验时某铁磁性运动载体(在测量范围内可近似等效为磁偶极子)以60 km·h^-1的速度在道路上直线运动，磁梯度张量系统放置于道路一侧的地面上，同步采集各传感器的输出数据，采样频率为2048 Hz.选取两条行驶路线进行目标轨迹的估计，假设运动载体沿着y轴正方向前进，在x方向上距离磁梯度张量系统的距离(Closest Proximity Approach，CPA)分别为1.4 m和2 m，z轴方向竖直向下构建右手坐标系进行数据计算和分析.

由于磁梯度张量系统静止不动，因此可在无磁性目标的情况下采集背景地磁场，实时采集得到的磁传感器数据减去背景地磁场即为磁性目标产生的磁场数据.在此，利用总磁场强度的变化判断是否有磁性目标的出现，然后截取采集得到的实时信号进行运动轨迹的估计.由于运动载体等效的磁偶极子位于测量系统的上方，可直接去除两个虚假解，则反演得到的载体在x-y平面的运动轨迹如图 8、图 9所示.利用图 6所示的方法去除剩余的虚假解，则图 8、图 9所示右侧深色曲线为估计得到的运动轨迹的真实解.

图 8

Fig. 8

图 8 估计得到的载体运动轨迹(CPA=1.4 m)Fig. 8 Estimated track of mobile vehicle when CPA equal to 1.4 m

图 9

Fig. 9

图 9 估计得到的载体运动轨迹(CPA=2 m)Fig. 9 Estimated track of mobile vehicle when CPA equal to 2 m

由图 8和图 9可知，虚假解表示的铁磁体运动轨迹与真实运动轨迹相差较大，利用图 6所示的虚假解去除流程有效地识别出了由磁梯度张量数据计算得到的真实解.在可探测范围内，利用几何不变量估计得到的铁磁性载体运动轨迹的真实解可较好地跟踪运动载体的实际轨迹，可实现基于磁梯度张量数据的运动载体实时跟踪.由于磁传感器自身的噪声以及实验场地中地磁梯度的存在，导致估计得到的目标位置与真实值存在一定的偏差，可通过进一步提高磁传感器的精度实现目标的远距离跟踪.

(1)基于数学解析方法对磁偶极子产生的磁梯度张量矩阵的特征值及特征向量进行了理论推导，得出了磁偶极子位置矢量与磁矩矢量间的夹角关于特征值的数学表达式，并在此基础上给出了磁矩模及磁偶极子与测量点间的距离的求解公式.

(2)基于磁偶极子产生的矢量场、磁梯度张量场及其欧拉反褶积公式，验证了磁偶极子位置矢量、磁矩矢量、最大特征值对应的特征向量及最小特征值对应的特征向量共面的假设，并给出了位置矢量及磁矩矢量关于特征值及特征向量的数学表达式.

(3)对本文提出的几何不变量在磁性目标跟踪上的应用进行了仿真和实测分析，实验结果证明了公式推导的正确性及其在目标跟踪中的实用性.