CN112611310A - 一种磁偶极子目标测距测向方法 - Google Patents

Abstract

本发明公开了一种磁偶极子目标测距测向方法，使用测量轴对称配置的三轴磁强计十字阵列测量阵列中心处的磁偶极子磁梯度张量，依据磁矩大小与磁性目标姿态无关的特点可以得到被探测目标的磁矩大小范围，由该点磁梯度张量的比值以及目标磁矩大小范围反演出被探测目标与十字阵列中心的距离和十字阵列测量坐标系下的磁偶极子目标的方位角。本发明所提的方法无需磁梯度张量测量***对磁性目标的磁场矢量的测量，也无需磁梯度张量测量***移动以测量空间多个点的磁梯度张量值；本发明所提的方法能使用三轴磁强计十字阵列作为磁梯度张量测量***，所使用的三轴磁强计相对较少，降低了测量***的复杂性和重量，便于轻型平台搭载使用。

Description

一种磁偶极子目标测距测向方法

技术领域

本发明涉及一种磁偶极子目标测距测向方法，特别是一种基于磁梯度张量分量比值的磁偶极子目标测距测向方法，属于磁性目标探测及其参数反演技术领域

背景技术

基于磁异常场的被动式目标探测技术具有隐蔽性好、能连续侦测、效率高、使用简单可靠、反应迅速等优点，还不受天文、气象和水文等的影响，可作为声学探测手段的有益补充，在水下目标探测与追踪、矿产勘探、体内微型诊疗装置的无创定位等方面有着重要的应用价值。

作为共轭抑制背景磁场的差分测量方式，磁梯度张量的测量结果受地磁场的倾角和偏角的影响小，所构成的张量不变量无需特殊处理就可很好地描述磁场源。因此，基于磁梯度张量的磁性目标参数反演技术得到了关注与研究。

Wiegert等提出了张量约缩量STAR的定位方法，张量约缩量为球等值面的不变量(Roy Wiegert，J.Oeschger.Generalized magnetic gradient contraction basedmethod for detection,localization and discrimination of underwater mines andunexploded ordnance[J].Mine Warfare&Ship Self-Defence,2010,2:1325-1332)。Clark利用磁梯度张量特征值计算归一化磁源强度，结合梯度张量确定磁源位置的唯一解(DavidA Clark.New methods for interpretation of magnetic vector and gradient tensordata I:Eigenvector analysis and the normalized source strength[J].Explorationgeophysics,2012,43(4):267-282)。Teixeira等利用张量矩阵最小特征值对应的特征向量与相对位置矢量之间的正交关系，提出一种张量欧拉反褶积和磁场梯度张量特征分析的联合反演方法，增强了反演算法的收敛性(Francisco Curado Teixeira,AntónioPascoal.Magnetic navigation and tracking of underwater vehicles[J].IFACProceedings Volumes,2013,46(33):239-244)。Lee等用磁梯度张量导出的两个标量参数来表征目标的磁矩和位置，再由梯度法求解位置并实验验证这种方法(Kok-Meng Lee,MinLi.Magnetic tensor sensor for gradient-based localization of ferrous objectin geomagnetic field[J].IEEE Transactions on Magnetics,2016,52(8):1-10)。

海军工程大学提出了一种基于等效衰减磁矩的定位方法，分析了地磁场起伏、三轴磁强计噪声及测量基线长等对定位结果的影响(贾文抖,林春生,孙玉绘,翟国君.地磁场中磁目标的等效衰减磁矩定位方法[J].兵工学报,2018,39(2):283-289)。海军航空工程学院提出了一种基于正六面体磁梯度张量***的单点定位改进方法，消除了张量约缩量STAR法的***误差(吕俊伟,迟铖,于振涛,毕波,宋庆善.磁梯度张量不变量的椭圆误差消除方法研究[J].物理学报,2015,64(9):52-59)。陆军工程大学研究了一种两点磁梯度张量定位算法并进行了铁磁体定位实验(Gang Yin,Yingtang Zhang,Zhining Li,Hongbo Fan andGuoquan Ren.Detection of ferromagnetic target based on mobile magneticgradient tensor system[J].Journal of Magnetism and Magnetic Material,2016,402(6):1-7)。

通常的磁梯度张量定位方法需要测量磁性目标的磁场矢量或者测量磁性目标在空间多个点的磁梯度张量。一方面，在地球背景场存在的条件下不易测量磁性目标的磁场矢量；另一方面，要想测量空间多个点的磁梯度张量需要移动磁梯度张量测量***或者用更多的三轴磁强计构建合适的测量阵列，这增加了测量***的复杂性。

发明内容

针对上述现有技术，本发明要解决的技术问题是提供一种基于磁梯度张量分量比值的磁偶极子目标测距测向方法，所使用的三轴磁强计相对较少，降低测量***的复杂性和重量，便于轻型平台搭载使用。

为解决上述技术问题，本发明的一种磁偶极子目标测距测向方法，包括以下步骤：

步骤1、将四个同规格的三轴磁强计按照测量轴对称配置方式安装在由非磁性材料制作而成的十字基座上，构成三轴磁强计十字阵列；四个三轴磁强计位置点分别为A、B、C和D点，其中A点和C点关于y轴对称，B点和D点关于x轴对称；

步骤2、读取四个三轴磁强计的测量输出，计算十字阵列中心处的磁梯度张量的独立分量值a、b、d、e和f，具体为：

式中，δ_x和δ_y分别为三轴磁强计的x_s轴和y_s轴偏离z_s轴的位置误差，L_x为十字阵列的点A和点C之间的长度，L_y为点B和点D之间的长度，B_Ax、B_Ay和B_Az为点A处的磁性体磁场的三分量，B_Bx、B_By和B_Bz为点B处的磁性体磁场的三分量，B_Cx、B_Cy和B_Cz为点C处的磁性体磁场的三分量，B_Dx、B_Dy和B_Dz为点D处的磁性体磁场的三分量；

步骤3、设位置坐标向量X＝[x,y,z]^T∈R³，构造关于X的方程式，建立优化目标函数F(X)，具体为：

F(X)＝|f₁(X)|+|f₂(X)|+|f₃(X)|+|f₄(X)|

式中，

k_ax＝(3r²-5x²)x，k_ay＝(r²-5x²)y，k_az＝(r²-5x²)z，k_bx＝(r²-5y²)x，k_by＝(3r²-5y²)y，k_bz＝(r²-5y²)z，k_dx＝k_ay，k_dy＝k_bx，k_dz＝-5xyz，k_ex＝k_dz，k_ey＝k_bz，k_ez＝y(r²-5z²)，k_fx＝k_az，k_fy＝k_dz，k_fz＝(r²-5z²)x，η_ab＝a/b，η_ad＝a/d，η_ae＝a/e和η_af＝a/f，x、y和z为磁梯度张量的测量点相对于磁偶极子的空间坐标，r为磁梯度张量的测量点与磁偶极子之间的距离；

步骤4、根据探测目标的类型确定磁矩大小的范围[m_min,m_max]，求解约束最优化问题计算磁偶极子目标与三轴磁强计十字阵列中心之间的距离

约束最优化问题具体为：

式中，h₁(X)＝xv₃₁+yv₃₂+zv₃₃，v₃₁、v₃₂和v₃₃是磁梯度张量矩阵的第三个特征值对应的特征向量元素；

式中，μ为磁导率，

λ₁、λ₂和λ₃为磁梯度张量矩阵的三个特征值；

步骤5、计算得到磁偶极子磁矩的三分量m_x、m_y和m_z，具体为：

式中，矩阵[K_ij]的元素分别为K₁₁＝(3r²-5x²)x、K₁₂＝(r²-5x²)y、K₁₃＝(r²-5x²)z、K₂₁＝(r²-5y²)x、K₂₂＝(3r²-5y²)y、K₂₃＝(r²-5y²)z、K₃₁＝K₁₂、K₃₂＝K₂₁、K₃₃＝-5xyz、K₄₁＝K₃₃、K₄₂＝K₂₃、K₄₃＝(r²-5z²)y、K₅₁＝K₁₃、K₅₂＝K₃₃和K₅₃＝(r²-5z²)x；

步骤6、求解关于k的多项式，得到k的六个根，具体为：

A₆k⁶+A₅k⁵+A₄k⁴+A₃k³+A₂k²+A₁k+A₀＝0

式中，A₆＝d²(a+2b)-e²(a-b)+2def，A₅＝-2d[(a-b)(a+2b)+(d²+e²+f²)]，A₄＝(a-b)²(a+2b)+d²(4a-7b)+(f²-2e²)(a-b)+6def，A₃＝-4d[(a-b)²+(-d²+e²+f²)]，A₂＝(a-b)²(2a+b)+d²(4b-7a)+(2f²-e²)(a-b)+6def，A₁＝2d[(a-b)(2a+b)-(d²+e²+f²)]，A₀＝d²(2a+b)+f²(a-b)+2def。

步骤7、舍弃k的复数根，根据k的剩余实数根计算q值，具体为：

步骤8、根据十字对称阵列的中心点位置处于偶极子磁性目标的上方或下方以及建立的测量坐标系确定z的正负，并根据k、q和t的值计算z值，具体为：

其中t是磁偶极子磁场大小，当十字对称阵列的中心点位置处于偶极子磁性目标的上方，z为负，当十字对称阵列的中心点位置处于偶极子磁性目标的下方，z为正。

分别计算出x和y，具体为：

x＝kqz

y＝qz

得到位置坐标的实数集合{x_l,y_l,z_l|l＝1,2,…,N_r}，N_r为k的实数根的个数；

步骤9、利用{x_l,y_l,z_l|l＝1,2,…,N_r}计算出误差集合{Δ_l}，具体为：

步骤10、从{Δ_l}寻找最小值对应的l，即

则k的值选取为：

k＝k_L

然后得到q和z的值；

步骤11、确定x和y的符号，具体为：

式中，sign(·)表示求符号值，磁偶极子目标相当于阵列中心点P的位置坐标x_Q和y_Q的符号分别为sign(x_Q)＝-sign(x)和sign(y_Q)＝-sign(y)；

计算磁偶极子目标的方位角ψ∈[0,2π)，具体为：

当sign(x_Q)为“+”且sign(y_Q)为“+”时，ψ＝arctan(1/k)；

当sign(x_Q)为“+”且sign(y_Q)为“-”时，ψ＝2π+arctan(1/k)；

当sign(x_Q)为“-”且sign(y_Q)为“+”或“-”时，ψ＝π+arctan(1/k)；

其中，当|x|-0≤A时且sign(y_Q)为“+”时，ψ＝π/2，；当|x|-0≤A时且sign(y_Q)为“-”时，ψ＝3π/2，A为给定值，A≤10^-4；

步骤12、利用距离值r计算得到磁偶极子磁场和磁矩的大小，具体为：

式中，

本发明的有益效果：本发明提出一种基于磁梯度张量分量比值的磁偶极子目标测距测向方法，使用测量轴对称配置的三轴磁强计十字阵列测量阵列中心处的磁偶极子磁梯度张量，依据磁矩大小与磁性目标姿态无关的特点可以得到被探测目标的磁矩大小范围，由该点磁梯度张量的比值以及目标磁矩大小范围反演出被探测目标与十字阵列中心的距离和十字阵列测量坐标系下的磁偶极子目标的方位角。本发明所提的方法无需磁梯度张量测量***对磁性目标的磁场矢量的测量，也无需磁梯度张量测量***移动以测量空间多个点的磁梯度张量值；本发明所提的方法能使用三轴磁强计十字阵列作为磁梯度张量测量***，所使用的三轴磁强计相对较少，降低了测量***的复杂性和重量，便于轻型平台搭载使用。

相比于现有的基于磁梯度张量的磁偶极子参数反演方法，本发明提出一种基于磁梯度张量分量比值的磁偶极子目标测距测向方法，无需磁梯度张量测量***测量磁性目标的磁场矢量，也无需磁梯度张量测量***移动以测量空间多个点的磁梯度张量值，拓展了方法的应用对象；本发明所提的方法能使用三轴磁强计十字阵列作为磁梯度张量测量***，所使用的三轴磁强计相对较少，降低了测量***的复杂性和重量，便于轻型平台搭载使用，具有很好的实际应用价值；所提的方法通过磁梯度张量分量之间的比值运算消减了介质磁导率的影响。

附图说明

图1是三轴磁强计敏感轴对称配置的十字阵列；

图2(a)是测距误差随三轴磁强计噪声标准差的变化曲线；

图2(b)是测向误差随三轴磁强计噪声标准差的变化曲线；

图3(a)是测距误差随磁矩大小的初始估算误差的变化曲线；

图3(b)是测向误差随磁矩大小的初始估算误差的变化曲线。

具体实施方式

下面结合附图对本发明具体实施方式做进一步说明。

本发明的一种基于磁梯度张量分量比值的磁偶极子目标测距测向方法的实现原理具体为：

磁性体磁场是无源无旋场，故其磁梯度张量矩阵是迹为0的对称阵，即磁梯度张量矩阵只有5个独立分量。令

和

其中B_x、B_y和B_z为磁性体磁场的三分量。

三轴磁强计的三个正交测量轴通常不相交于一点，也就是说三轴磁强计的x_s轴与其z_s轴之间存在δ_x的位置偏差，三轴磁强计的y_s轴与其z_s轴之间存在δ_y的位置偏差。为了减少测量误差，采用如图1所示的三轴磁强计测量轴对称配置的十字阵列测量阵列中心处的磁梯度张量的独立分量值，四个三轴磁强计的z_s轴与十字阵列坐标系的z轴一致，且其分别放置在点A、B、C和D处。

点A和点C之间的长度为L_x，点B和点D之间的长度为L_y，则磁梯度张量矩阵的5个独立分量a、b、d、e和f的测量值分别为

式中，B_Ax、B_Ay和B_Az为点A处的磁性体磁场的三分量，B_Bx、B_By和B_Bz为点B处的磁性体磁场的三分量，B_Cx、B_Cy和B_Cz为点C处的磁性体磁场的三分量，B_Dx、B_Dy和B_Dz为点D处的磁性体磁场的三分量。

对于磁偶极子磁性体而言，其磁梯度张量矩阵的5个独立分量a、b、d、e和f的表达式分别为

式中，μ为磁导率，x、y和z为磁梯度张量的测量点相对于磁偶极子的空间坐标，r为磁梯度张量的测量点与磁偶极子之间的距离，m_x、m_y和m_z为磁偶极子磁矩的三分量。

根据式(2)，求得独立分量之间的比值分别为

式中，k_ax＝(3r²-5x²)x，k_ay＝(r²-5x²)y，k_az＝(r²-5x²)z，k_bx＝(r²-5y²)x，k_by＝(3r²-5y²)y，k_bz＝(r²-5y²)z，k_dx＝k_ay，k_dy＝k_bx，k_dz＝-5xyz，k_ex＝k_dz，k_ey＝k_bz，k_ez＝y(r²-5z²)，k_fx＝k_az，k_fy＝k_dz，k_fz＝(r²-5z²)x。

由于m_x、m_y和m_z不可能全部为零，因此由式(3)可得

式(4)可改写为

式中，X＝[x,y,z]^T∈R³为阵列中心点相当于磁偶极子的位置坐标向量。

式(5)可改写为

式(6)可改写为

式(7)可改写为

磁偶极子的磁梯度张量场还存在式(12)所示的关系。

h₁(X)＝xv₃₁+yv₃₂+zv₃₃＝0 (12)

式中，v₃₁、v₃₂和v₃₃是磁梯度张量矩阵的第三个特征值对应的特征向量元素。

根据探测目标的类型和归一化强度u确定探测距离的范围，即有

式中，

λ₁、λ₂和λ₃为磁梯度张量矩阵的三个特征值，m_max和m_min为磁偶极子磁矩大小m的最大值和最小值。

通过求解式(14)所示的约束最优化问题计算磁偶极子目标与三轴磁强计十字阵列中心之间的距离。

式中，F(X)＝|f₁|+|f₂|+|f₃|+|f₄|。

采用序列二次规划(SQP)算法求解式(14)所示的约束优化问题，得到磁偶极子目标相对于十字阵列中心点的距离值

再由式(15)计算出m_x、m_y和m_z。

式中，矩阵[K_ij]的元素分别为K₁₁＝(3r²-5x²)x、K₁₂＝(r²-5x²)y、K₁₃＝(r²-5x²)z、K₂₁＝(r²-5y²)x、K₂₂＝(3r²-5y²)y、K₂₃＝(r²-5y²)z、K₃₁＝K₁₂、K₃₂＝K₂₁、K₃₃＝-5xyz、K₄₁＝K₃₃、K₄₂＝K₂₃、K₄₃＝(r²-5z²)y、K₅₁＝K₁₃、K₅₂＝K₃₃和K₅₃＝(r²-5z²)x。

求解式(16)所示的关于k的多项式，得到k的六个根。

A₆k⁶+A₅k⁵+A₄k⁴+A₃k³+A₂k²+A₁k+A₀＝0 (16)

式中，A₆＝d²(a+2b)-e²(a-b)+2def，A₅＝-2d[(a-b)(a+2b)+(d²+e²+f²)]，A₄＝(a-b)²(a+2b)+d²(4a-7b)+(f²-2e²)(a-b)+6def，A₃＝-4d[(a-b)²+(-d²+e²+f²)]，A₂＝(a-b)²(2a+b)+d²(4b-7a)+(2f²-e²)(a-b)+6def，A₁＝2d[(a-b)(2a+b)-(d²+e²+f²)]，A₀＝d²(2a+b)+f²(a-b)+2def。

舍弃k的复数根，将k的剩余实数根代入式(17)计算q值。

根据十字对称阵列的中心点位置处于偶极子磁性目标的上方或下方以及建立的测量坐标系，判断z的正负，以此确定式(18)中的正负号，并将k、q和t的值代入式(18)计算z值。

再根据式(19)和式(20)分别计算出

x＝kqz (19)

y＝qz (20)

得到位置坐标的实数集合{x_l,y_l,z_l|l＝1,2,…,N_r}，N_r为k的实数根的个数，将这些空间坐标的实数解代入式(21)计算出误差集合{Δ_l}。

从{Δ_l}寻找最小值对应的l，即

则k的值选取为

k＝k_L (23)

由此也确定出q和z的值。

根据式(24)确定x和y的符号，即

式中，sign(·)表示求符号值。磁偶极子目标相当于阵列中心点P的位置坐标x_Q和y_Q的符号分别为sign(x_Q)＝-sign(x)和sign(y_Q)＝-sign(y)。

由表1所描述的计算方法计算出磁偶极子目标的方位角ψ∈[0,2π)。

表1方位角ψ的计算公式

将距离值r代入式(25)和式(26)还可计算出磁偶极子磁场大小t和磁矩大小m。

式中，

和

都为磁梯度张量不变量。

结合图1，本发明具体实施方式包括以下步骤：

步骤1、将四个同规格的三轴磁强计按照测量轴对称配置方式安装在由非磁性材料制作而成的十字基座上，构成如图1所示的三轴磁强计十字阵列。

步骤2、读取四个三轴磁强计的测量输出，根据式(1)计算十字阵列中心处的磁梯度张量的独立分量值a、b、d、e和f。

式中，δ_x和δ_y分别为三轴磁强计的x_s轴和y_s轴偏离z_s轴的位置误差，L_x为十字阵列的点A和点C之间的长度，L_y为点B和点D之间的长度，B_Ax、B_Ay和B_Az为点A处的磁性体磁场的三分量，B_Bx、B_By和B_Bz为点B处的磁性体磁场的三分量，B_Cx、B_Cy和B_Cz为点C处的磁性体磁场的三分量，B_Dx、B_Dy和B_Dz为点D处的磁性体磁场的三分量。

步骤3、设位置坐标向量X＝[x,y,z]^T∈R³，构造关于X的方程式，按式(2)建立优化目标函数F(X)。

F(X)＝|f₁(X)|+|f₂(X)|+|f₃(X)|+|f₄(X)| (2)

式中，

k_ax＝(3r²-5x²)x，k_ay＝(r²-5x²)y，k_az＝(r²-5x²)z，k_bx＝(r²-5y²)x，k_by＝(3r²-5y²)y，k_bz＝(r²-5y²)z，k_dx＝k_ay，k_dy＝k_bx，k_dz＝-5xyz，k_ex＝k_dz，k_ey＝k_bz，k_ez＝y(r²-5z²)，k_fx＝k_az，k_fy＝k_dz，k_fz＝(r²-5z²)x，η_ab＝a/b，η_ad＝a/d，η_ae＝a/e和η_af＝a/f。

步骤4、根据探测目标的类型确定磁矩大小的范围[m_min,m_max]，求解式(3)所示的约束最优化问题计算磁偶极子目标与三轴磁强计十字阵列中心之间的距离

式中，h₁(X)＝xv₃₁+yv₃₂+zv₃₃，v₃₁、v₃₂和v₃₃是磁梯度张量矩阵的第三个特征值对应的特征向量元素。

式中，μ为磁导率，

λ₁、λ₂和λ₃为磁梯度张量矩阵的三个特征值。

步骤5、按式(5)计算出m_x、m_y和m_z。

式中，矩阵[K_ij]的元素分别为K₁₁＝(3r²-5x²)x、K₁₂＝(r²-5x²)y、K₁₃＝(r²-5x²)z、K₂₁＝(r²-5y²)x、K₂₂＝(3r²-5y²)y、K₂₃＝(r²-5y²)z、K₃₁＝K₁₂、K₃₂＝K₂₁、K₃₃＝-5xyz、K₄₁＝K₃₃、K₄₂＝K₂₃、K₄₃＝(r²-5z²)y、K₅₁＝K₁₃、K₅₂＝K₃₃和K₅₃＝(r²-5z²)x。

步骤6、求解式(6)所示的关于k的多项式，得到k的六个根。

A₆k⁶+A₅k⁵+A₄k⁴+A₃k³+A₂k²+A₁k+A₀＝0 (6)

式中，A₆＝d²(a+2b)-e²(a-b)+2def，A₅＝-2d[(a-b)(a+2b)+(d²+e²+f²)]，A₄＝(a-b)²(a+2b)+d²(4a-7b)+(f²-2e²)(a-b)+6def，A₃＝-4d[(a-b)²+(-d²+e²+f²)]，A₂＝(a-b)²(2a+b)+d²(4b-7a)+(2f²-e²)(a-b)+6def，A₁＝2d[(a-b)(2a+b)-(d²+e²+f²)]，A₀＝d²(2a+b)+f²(a-b)+2def。

步骤7、舍弃k的复数根，将k的剩余实数根代入式(7)计算q值。

步骤8、根据十字对称阵列的中心点位置处于偶极子磁性目标的上方或下方以及建立的测量坐标系，判断z的正负，以此确定式(8)中的正负号，并将k、q和t的值代入式(8)计算z值。

再根据式(9)和式(10)分别计算出

x＝kqz (9)

y＝qz (10)

得到位置坐标的实数集合{x_l,y_l,z_l|l＝1,2,…,N_r}，N_r为k的实数根的个数。

步骤9、将{x_l,y_l,z_l|l＝1,2,…,N_r}代入式(11)计算出误差集合{Δ_l}。

步骤10、从{Δ_l}寻找最小值对应的l，即

则k的值选取为

k＝k_L (13)

由此也确定出q和z的值。

步骤11、根据式(14)确定x和y的符号，即

式中，sign(·)表示求符号值。磁偶极子目标相当于阵列中心点P的位置坐标x_Q和y_Q的符号分别为sign(x_Q)＝-sign(x)和sign(y_Q)＝-sign(y)。

由表1所描述的计算方法计算出磁偶极子目标的方位角ψ∈[0,2π)。

表1方位角ψ的计算公式

步骤12、将距离值r代入式(15)和式(16)得到磁偶极子磁场和磁矩的大小。

式中，

和

为检验本发明提出的测距测向方法的可行性，进行如下数值仿真实验。磁偶极子的磁矩三分量分别为m_x＝8×10⁸A·m、m_y＝6×10⁸A·m和m_z＝3×10⁸A·m，其位置坐标真值为x_Q＝180m、y_Q＝220m和z_Q＝-150m，三轴磁强计十字阵列在x和y方向上的基线长分别为L_x＝0.8m和L_y＝0.8m，测量轴的位置偏差为δ_x＝2mm和δ_y＝2mm。三轴磁强计在每个轴的噪声是相互独立的高斯过程，其均值为0，标准差为σ_m。磁偶极子磁矩大小的初始估算误差为ρ。

仿真实例一：设磁偶极子磁矩大小的初始估算误差为10％，即

和

在不同的三轴磁强计噪声标准差的情况下，进行磁偶极子目标测距与测向的50次蒙特卡洛仿真实验。得到三轴磁强计十字阵列对磁偶极子目标的测距与测向误差随三轴磁强计噪声标准差的变化曲线如图2(a)和2(b)所示，其中图2(a)为σ_m从0nT增加到0.8nT时的距离测量误差曲线，图2(b)为σ_m从0nT增加到0.8nT时的方向角测量误差曲线。由图2(a)和2(b)可知，距离测量误差随σ_m的增加而呈增大趋势，但增大的幅度不大；方向角测量误差随σ_m的增加而呈线性增大趋势。

仿真实例二：设三轴磁强计噪声标准差σ_m为0.1nT，在不同的磁偶极子磁矩大小的初始估算误差的情况下，进行磁偶极子目标测距与测向的50次蒙特卡洛仿真实验。得到三轴磁强计十字阵列对磁偶极子目标的测距与测向误差随磁偶极子磁矩大小的初始估算误差的变化曲线如图3(a)和3(b)所示，其中图3(a)为ρ从2％增加到20％时的距离测量误差曲线，图3(b)为ρ从2％增加到20％时的方向角测量误差曲线。由图3(a)和3(b)可知，在ρ较大时距离测量误差随ρ的增加而呈线性增大趋势，在ρ小时距离测量误差随ρ的增加而几乎不变；方向角测量误差随ρ的增加而变化不大。

以测量轴对称配置的三轴磁强计十字阵列为磁梯度张量测量装置，进行了磁偶极子磁性目标的测距测向的蒙特卡洛仿真实验。仿真实验给出了在不同的三轴磁强计噪声标准差的情况下三轴磁强计十字阵列对磁偶极子目标的测距与测向误差随三轴磁强计噪声标准差的变化曲线，结果表明距离测量误差随三轴磁强计测量噪声标准差的增加而呈增大趋势，但增大的幅度不大；方向角测量误差随轴磁强计测量噪声标准差的增加而呈线性增大趋势。仿真实验给出了在不同的磁偶极子磁矩大小的初始估算误差的情况下，三轴磁强计十字阵列对磁偶极子目标的测距与测向误差随磁偶极子磁矩大小的初始估算误差的变化曲线，结果表明在磁矩大小的初始估算误差较大时距离测量误差随初始估算误差的增加而呈线性增大趋势，在磁矩大小的初始估算误差小时距离测量误差随初始估算误差的增加而几乎不变；方向角测量误差随初始估算误差的增加而变化不大。仿真实验证明了本发明方法的可行性。

相比于一般的基于磁梯度张量的磁偶极子参数反演方法，本发明提出的方法无需磁梯度张量测量***测量磁性目标的磁场矢量，也无需磁梯度张量测量***移动以测量空间多个点的磁梯度张量值，拓展了方法的应用对象；本发明所提的磁梯度张量测量方案所需的三轴磁强计相对较少，构造简单，降低了测量***的复杂性和重量，便于轻型平台搭载使用，三轴磁强计测量轴对称配置也降低了磁梯度张量差分式测量的误差，具有很好的实际应用价值；磁梯度张量分量之间的比值运算也消减了介质磁导率的影响。

Claims (1)

1.一种磁偶极子目标测距测向方法，其特征在于，包括以下步骤：

步骤1、将四个同规格的三轴磁强计按照测量轴对称配置方式安装在由非磁性材料制作而成的十字基座上，构成三轴磁强计十字阵列；四个三轴磁强计位置点分别为A、B、C和D点，其中A点和C点关于y轴对称，B点和D点关于x轴对称；

步骤2、读取四个三轴磁强计的测量输出，计算十字阵列中心处的磁梯度张量的独立分量值a、b、d、e和f，具体为：

式中，δ_x和δ_y分别为三轴磁强计的x_s轴和y_s轴偏离z_s轴的位置误差，L_x为十字阵列的点A和点C之间的长度，L_y为点B和点D之间的长度，B_Ax、B_Ay和B_Az为点A处的磁性体磁场的三分量，B_Bx、B_By和B_Bz为点B处的磁性体磁场的三分量，B_Cx、B_Cy和B_Cz为点C处的磁性体磁场的三分量，B_Dx、B_Dy和B_Dz为点D处的磁性体磁场的三分量；

步骤3、设位置坐标向量X＝[x,y,z]^T∈R³，构造关于X的方程式，建立优化目标函数F(X)，具体为：

F(X)＝|f₁(X)|+|f₂(X)|+|f₃(X)|+|f₄(X)|

式中，

k_ax＝(3r²-5x²)x，k_ay＝(r²-5x²)y，k_az＝(r²-5x²)z，k_bx＝(r²-5y²)x，k_by＝(3r²-5y²)y，k_bz＝(r²-5y²)z，k_dx＝k_ay，k_dy＝k_bx，k_dz＝-5xyz，k_ex＝k_dz，k_ey＝k_bz，k_ez＝y(r²-5z²)，k_fx＝k_az，k_fy＝k_dz，k_fz＝(r²-5z²)x，η_ab＝a/b，η_ad＝a/d，η_ae＝a/e和η_af＝a/f，x、y和z为磁梯度张量的测量点相对于磁偶极子的空间坐标，r为磁梯度张量的测量点与磁偶极子之间的距离；

步骤4、根据探测目标的类型确定磁矩大小的范围[m_min,m_max]，求解约束最优化问题计算磁偶极子目标与三轴磁强计十字阵列中心之间的距离

约束最优化问题具体为：

式中，h₁(X)＝xv₃₁+yv₃₂+zv₃₃，v₃₁、v₃₂和v₃₃是磁梯度张量矩阵的第三个特征值对应的特征向量元素；

式中，μ为磁导率，

λ₁、λ₂和λ₃为磁梯度张量矩阵的三个特征值；

步骤5、计算得到磁偶极子磁矩的三分量m_x、m_y和m_z，具体为：

式中，矩阵[K_ij]的元素分别为K₁₁＝(3r²-5x²)x、K₁₂＝(r²-5x²)y、K₁₃＝(r²-5x²)z、K₂₁＝(r²-5y²)x、K₂₂＝(3r²-5y²)y、K₂₃＝(r²-5y²)z、K₃₁＝K₁₂、K₃₂＝K₂₁、K₃₃＝-5xyz、K₄₁＝K₃₃、K₄₂＝K₂₃、K₄₃＝(r²-5z²)y、K₅₁＝K₁₃、K₅₂＝K₃₃和K₅₃＝(r²-5z²)x；

步骤6、求解关于k的多项式，得到k的六个根，具体为：

A₆k⁶+A₅k⁵+A₄k⁴+A₃k³+A₂k²+A₁k+A₀＝0

式中，A₆＝d²(a+2b)-e²(a-b)+2def，A₅＝-2d[(a-b)(a+2b)+(d²+e²+f²)]，A₄＝(a-b)²(a+2b)+d²(4a-7b)+(f²-2e²)(a-b)+6def，A₃＝-4d[(a-b)²+(-d²+e²+f²)]，A₂＝(a-b)²(2a+b)+d²(4b-7a)+(2f²-e²)(a-b)+6def，A₁＝2d[(a-b)(2a+b)-(d²+e²+f²)]，A₀＝d²(2a+b)+f²(a-b)+2def。

步骤7、舍弃k的复数根，根据k的剩余实数根计算q值，具体为：

步骤8、根据十字对称阵列的中心点位置处于偶极子磁性目标的上方或下方以及建立的测量坐标系确定z的正负，并根据k、q和t的值计算z值，具体为：

其中t是磁偶极子磁场大小，当十字对称阵列的中心点位置处于偶极子磁性目标的上方，z为负，当十字对称阵列的中心点位置处于偶极子磁性目标的下方，z为正。

分别计算出x和y，具体为：

x＝kqz

y＝qz

得到位置坐标的实数集合{x_l,y_l,z_l|l＝1,2,…,N_r}，N_r为k的实数根的个数；

步骤9、利用{x_l,y_l,z_l|l＝1,2,…,N_r}计算出误差集合{Δ_l}，具体为：

步骤10、从{Δ_l}寻找最小值对应的l，即

则k的值选取为：

k＝k_L

然后得到q和z的值；

步骤11、确定x和y的符号，具体为：

式中，sign(·)表示求符号值，磁偶极子目标相当于阵列中心点P的位置坐标x_Q和y_Q的符号分别为sign(x_Q)＝-sign(x)和sign(y_Q)＝-sign(y)；

计算磁偶极子目标的方位角ψ∈[0,2π)，具体为：

当sign(x_Q)为“+”且sign(y_Q)为“+”时，ψ＝arctan(1/k)；

当sign(x_Q)为“+”且sign(y_Q)为“-”时，ψ＝2π+arctan(1/k)；

当sign(x_Q)为“-”且sign(y_Q)为“+”或“-”时，ψ＝π+arctan(1/k)；

其中，当|x|-0≤A时且sign(y_Q)为“+”时，ψ＝π/2，；当|x|-0≤A时且sign(y_Q)为“-”时，ψ＝3π/2，A为给定值，A≤10^-4；

步骤12、利用距离值r计算得到磁偶极子磁场和磁矩的大小，具体为：

式中，

Images