基于磁偶极子磁场分布特征的磁矩方向估算方法

，， ，,2

(1.海军工程大学，湖北 武汉 430033；2.海军海洋测绘研究所，天津 300061)

0 引言

磁矩是描述物体磁特性的一个重要参数，磁矩测量对于卫星的姿态控制[1-2]、磁性目标的消磁及磁防护[3]工作都具有重要的意义。

为解决磁性目标磁矩的反演问题，陈进明等[4]利用神经网络优化算法对磁矩反演进行了优化计算，并对反演过程中的鲁棒性较差的现象进行了自适应修正；洪咸东等[5]对水雷磁矩测量过程中的测量误差进行了建模，分析了测量距离、传感器长度以及传感器直径等因素对测量误差的影响；亓亮等[6]利用磁传感器阵列系统测量辐射体的空间磁场分布，然后通过计算解算出辐射体的磁矩，但该方法同样需要进行迭代收敛计算，求解的精度与选用的方法有关。这些磁矩反演的问题都转化为了非线性优化问题，直接根据磁场数据优化求解磁矩。

对于非线性问题的求解，其计算方法的选取具有一定的随机性，求解结果受所选用的方法以及初值的设定影响较大。为避免磁矩反演时的非线性优化问题，本文提出了基于磁偶极子磁场分布特征的磁矩方向反演方法。这里的磁场分布特征主要包含磁梯度张量[7-8]分布和磁场矢量分布两部分。

1 特征平面上磁梯度张量奇异性分析

设源点磁偶极子所在位置为坐标原点O，磁偶

极子磁矩为M=(mx,my,mz)，场点A(x,y,z)处的空间磁场为：

(1)

对式(1)中Bx、By、Bz分别在x、y、z三个方向上求偏导，得到磁梯度张量G中的9个量，磁梯度张量可表示为：

(2)

式中：

计算得到磁梯度张量G的3个特征值分别为

(3)

根据式(3)中特征值的表达式可知特征值λ2<0、λ3>0，若特征值λ1=0，即xmx+ymy+zmz=0，磁梯度张量为奇异阵，则此时场点位于磁偶极子的特征平面[9]上；反过来若场点位于特征平面上，则有xmx+ymy+zmz=0，此时特征值λ1=0，磁梯度张量为奇异阵。可知，所有满足磁梯度张量矩阵奇异的场点均位于特征平面P上，满足这一条件的场点位置矢量r与磁矩M是垂直的。

2 磁矩方向估算方法

2.1 特征平面法向量的计算

如图1所示，O为磁偶极子源点，M是磁矩，平面P是磁偶极子的特征平面，位于平面P上的两条直线L1和L2称为特征线。特征线L1、L2与磁矩M是相互垂直的。

任意与特征平面不平行的面上的磁梯度张量奇异点都位于同一条直线上，这条直线是该平面与特征平面的相交线。因此若能计算得到不重合的两个平面上的特征线L1和L2，就可以唯一确定特征平面P的位置，进一步可求得特征平面P的法向量Im，法向量Im与磁矩M是平行关系，若要进一步使Im的方向与磁矩M的方向一致，需要结合磁偶极子磁场分布特征对Im方向进一步进行标定。

2.2 磁矩方向的计算

图2为磁偶极子磁场沿磁矩的剖视图，图中给出了磁力线及其指向。坐标原点O为磁源位置，以磁矩的指向为Z轴的正向，P是特征平面，A是场点，B为磁场矢量。

图中磁偶极子的特征平面P将磁场空间划分成为上下两部分，上半部分磁力线的走势呈现从磁源向四周发散的趋势，记为磁场发散区，在发散区内场点A位置失量OA与磁矩M的夹角小于90°；下半部分磁力线的走势呈现从四周向磁源汇聚的趋势，记为磁场汇聚区，在汇聚区内场点A位置矢量OA与磁矩M的夹角大于90°。失量OA与磁矩M的夹角计算公式：

(4)

场点A处的磁场B与场点位置矢量OA的夹角：

(5)

根据式(4)、式(5)可知，cosθ和cosβ的符号相同，仅与xm1+ym2+zm3的正负有关。因此，若cosβ>0，则有cosθ>0，即场点位于发散区，此时场点坐标矢量OA与磁矩M的夹角θ<90°；若cosβ<0，则有cosθ<0，即场点位于汇聚区，此时场点坐标矢量OA与磁矩M的夹角θ>90°。即在某一指定的场点A处，可以根据OA与B的夹角关系确定OA与磁矩M的夹角关系，进而对向量Im的方向进行判定，最后实现磁矩方向的估算。

3 仿真及误差分析

在实际的测量计算中，由于存在多种因素的影响使得测量数据与理论值之间存在差异，根据矩阵特征值λ1=0来判断某一矩阵是否为奇异阵往往是比较困难的，较为可行的做法是根据条件数的大小来判断矩阵是否奇异或接近于奇异，条件数大的为奇异矩阵或接近于奇异矩阵，若奇异，则该场点位于特征平面上。因此，可以在某一区域范围内取多个场点，根据场点处的磁梯度张量矩阵的条件数来判断该场点是否位于特征平面上。

设偶极子的磁矩M=(1 000，-1 430，3 600)A·m2，在空间Z=-3 m的平面上，取定矩形测量区域X、Y的范围为-19

结合图3中的三维视图分析可知，在计算空间中出现了磁张量矩阵条件数异常的现象，大部分场点处磁张量矩阵的条件数都比较小，不超过10，小部分场点处的条件数明显超过其他场点很多，达到几百甚至上千。从俯视图上可以看出，条件数较大的点基本分布在一条直线附近，根据上节的分析可以判断，这条直线是特征平面P与Z=-3平面的相交线，在此称为特征线。拟合得到特征线L1(图中斜线)的曲线方程：

同样，在Z=3的平面上选择同样大小的测量区域进行场点磁张量矩阵条件数的计算，将场点对应的条件数绘制于图4中。

与图3中的结果相同，图4显示出在Z=3的平面区域上也存在磁张量矩阵条件数异常的现象，部分场点的条件数明显超过其他场点很多，最大达到1 000以上。从俯视图上可以看出，这些条件数较大的场点近似分布在一条直线附近，这条直线是特征平面P与Z=3平面的相交线。拟合得到特征线L2(图中斜线)的曲线方程：

如图5所示，特征线L1和L2位于平面P上，磁矩M与L1和L2是相互垂直的关系，因此可根据特征线L1和L2的空间位置关系初步计算出与磁矩M平行的空间向量Im，即特征平面的法向量。

直线L1、L2在Z=0平面上投影的位置关系：

1)在Y轴截距相差：

Δb=21.598 8-6.463 6=15.135 2

2)投影直线的倾斜角(与X轴夹角)：

α=arctan(0.701 6)=0.611 8

3)两投影直线间距：

Δb′=Δbcosα=12.389 9

直线L1、L2所确定的特征平面P与Z=0平面的夹角：

β=arctan(Δh/Δb′)=
arctan((3-(-3))/13.810 3)=0.451 0

设Iy=±(0,1,0)是平行于Y轴的单位向量，将Iy首先绕Z轴逆时针旋转α角，然后再绕X轴逆时针旋转β角，即是特征平面P的一个法向量，同时也是一个与磁矩M平行的单位向量。根据欧拉旋转关系确定出磁矩M的方向：

(6)

下面结合磁场分布特征对Im的正负号进行判定。

选取坐标为(-19,-20,3)的场点A，对应的磁场B=(-11.747 9,-0.747 0,-15.732 7)nT，根据式(5)可以求得：

计算可知，当式(6)中Im取正号时可满足cosθ> 0，最终确定磁矩方向为Im=(1,-1.425 5,3.595 7)。

计算磁矩方向Im与设定磁矩M=(1 000,-1 430,3 600)A·m2之间的夹角为

4 结论

本文提出了基于磁偶极子磁场分布特征的磁矩方向估算方法，通过分析磁偶极子的磁场分布特性，利用磁梯度张量奇异点确定特征平面的位置，得到与磁矩相平行的特征平面的法向量，再根据位置矢量与磁矩、磁场矢量所成的两个夹角之间的关系确定法向量的具体指向，从而得到磁矩的方向。仿真结果表明，该方法有效避免了以往磁矩求解工作中需进行非线性优化的难题，实现了目标磁矩的反演。估算得到的磁矩方向与设定的磁矩方向基本一致的，两者之间的角度偏差仅为0.14°。

参考文献:

[1]吴卫权.地磁场中卫星不侧置态下磁矩测试方法[J].上海航天，2014，31(1)：50-55.

[2]易忠, 史尧宜.地磁场中测量卫星磁矩的一种新的方法[J].环模技术，1996(4)：14-24.

[3]张朝阳，虞伟乔，陆鹏飞.基于远场等效磁矩的潜艇磁防护技术[J].舰船科学技术，2015，37(2)：97-100.

[4]陈进明，隗燕琳.基于神经网络的磁性目标磁矩反演方法[J].船电技术，2012，32(9)：57-60.

[5] 洪咸东，任志良，俞伟，等.鱼雷辐射磁矩测量误差研究[J].鱼雷技术，2007，15(1)：42-44.

[6]亓亮，陆志宏，李迪.基于磁传感器阵列的电磁辐射体辐射磁矩测量[J].舰船电子对抗，2015，38(6)：109-112.

[7]李金朋，张英堂，范红波，等.基于磁梯度张量的地下小目标相关成像方法[J].探测与控制学报，2016，38(3)：75-78.

[8]李光，随阳轶，刘丽敏，等.基于差分的磁偶极子单点张量定位方法[J].探测与控制学报，2012，35(5)：50-54.

[9]任来平，欧阳永忠，陆秀平，等.水下铁磁体磁场特征平面[J].海洋测绘，2005，25(4)：1-4.