0 引言

水下铁磁性目标在地磁场中会引起磁异常，通过磁传感器阵列接收水下目标的磁异常信号并进行处理，可解算出目标的位置和磁矩信息^[1]，是水下目标探测的重要手段，在军事和民用领域都具有广泛应用。

国内外学者对铁磁性目标的位置、姿态及磁性参数求解的算法研究较多。张坚等^[2]利用正交基函数（OBF）算法结合 BP 神经网络在船舶磁场信号检测中取得了较好的效果，能够很好地对噪声进行过滤并突出磁异常信号的特征；杨勇等^[3]利用 OBF 算法大幅地提高了目标磁信号的信噪比，能够有效地检测出磁异常，从而易于检测磁性目标；WANG 等^[4]研究了三轴磁传感器接收到的磁异常信号的特征，提出一种正交能量比（OER）算法，通过单个磁传感器实现磁矩方向估计。ALIMI 等^[5]提出了一种基于 2 阶段 LM 算法的铁磁性目标检测与跟踪算法，适用于对 2~4 个三轴磁强计阵列附近移动的铁磁性物体的定位和磁矩估计；WANG 等^[6]针对便捷式瞬变电磁系统提出了一种利用磁梯度张量和高斯–牛顿法的地下目标快速定位的新算法。在这些方法中，高斯–牛顿法是一种常用的磁性目标反演定位数值迭代算法，但其只适用于小残量问题，在低信噪比的情况下容易失效且严重依赖于初始值的选择。LM 算法是对高斯–牛顿法的改进，在低信噪比条件下也能达到一定的精度，但若初始值选取不当，也会导致目标的定位精度降低甚至算法失效。针对高斯–牛顿法和 LM 算法严重依赖于初始值的问题，GE等^[7]利用粒子群算法（Particle Swarm Optimization，PSO）对初始值不敏感的特性获得粗略解，然后利用高斯–牛顿法得到更精确的解，开发了一种基于粒子群优化和高斯–牛顿方法相结合的混合算法（PSO-GN）；孙文^[8]利用粒子群算法的全局收敛性寻找初始值，并将初始值代入 LM 算法进行精确求解，提出一种基于 PSO 算法与 LM 算法的混合优化算法（PSO-LM），用于解决磁性目标定位问题。但粒子群算法容易陷入局部最优，导致收敛精度低或不易收敛，而且参数选择不当可能会增加运行时间、降低算法效率。

磁传感器接收的磁感应强度是目标的位置和磁矩的函数。由于未知参数多，在测量中存在的噪声干扰因素导致传感器测量值与真实值有误差，铁磁性目标的位置和磁矩事先未知，选取的初始值也可能远离目标的真实位置，这些因素使得传统高斯–牛顿法、LM 算法和智能优化算法都不能得到令人满意的定位精度和效率。所以，寻找一种对初始值敏感度较弱且能满足实时性要求的算法，仍然是一项富有挑战性的工作。本文将铁磁性目标产生的磁异常等效为磁偶极子，通过引入信赖域搜索技术对 LM 算法进行改进，并基于改进 LM 算法提出了一种改进 LM-GN 算法估计目标参数。该算法克服了高斯–牛顿法和 LM 算法严重依赖于初始值的缺点，在初始值离真实值较远的情况下，也能实现较准确的定位效果，降低了解对初始值的依赖性，提高了收敛速度，能够满足实时性需求，具有一定的实际工程应用价值。

1 磁定位模型分析

当磁传感器与磁性目标的距离大于 3 倍以上磁性目标几何尺寸时，可以将该磁性目标看作磁偶极子^[9]，见图1。

如图1 所示，设磁性目标的中心坐标为$Q\left(x_0，y_0，z_0\right)$，其磁矩矢量为 m，磁传感器的坐标为$P（x，y，z）$，磁性目标的中心点到传感器的矢径为 r，其距离为 r。由磁偶极子模型可知， P 处磁感应强度 B 的表达式^[10]为

式中，${\mu }_0$为真空中的磁导率，大小为$4\pi \times {10}^{-7}\mathrm{H}/\mathrm{m}$。

将式（1）在空间直角坐标系下展开得到式（2），式中$\left(B_x，B_y，B_z\right)$和 $\left(m_x，m_y，m_z\right)$分别代表了磁感应强度矢量和磁矩矢量在各坐标轴上的投影分量。

综上所述，对于一个三轴磁传感器，只能建立 3 个方程组，而模型中未知数的个数为 6 个。因此，至少需要 2 个传感器进行测量，构成 6 个方程才能实现方程组的求解，从而获得目标的位置和磁矩参数。假设传感器数量为 $N（N\ge 2）$，那么式（2）将变为由 3N 个方程构成的高度非线性、多极值的超定方程组。不妨设 $B_{x1}，B_{y1}，B_{z1}，\cdots ，B_{xN}，B_{yN}，B_{zN}$分别为测量阵列测得的实际磁感应强度值，而经过磁偶极子模型的理论计算值为${\stackrel{-}{B}}_{x1}，{\stackrel{-}{B}}_{y1}，{\stackrel{-}{B}}_{z1}，\cdots ，{\stackrel{-}{B}}_{xN}，{\stackrel{-}{B}}_{yN}，{\stackrel{-}{B}}_{zN}$。基于最小二乘准则，最小化实际测量值和理论计算值的平方误差，得到如下磁性目标定位目标函数：

使得式（3）成立的目标位置和磁矩值，即为需要的磁性目标参数的最优估计值。

2 非线性最小二乘问题的经典优化求解方法

迭代法是解决非线性最小二乘问题的一种有效的方法，其主要思想是先给定所求方程的一个初始解，然后依据迭代公式设定的步长不断对初始值进行修正，直到满足迭代要求为止^[11]。经典的优化方法包括高斯–牛顿法、LM 算法等。

2.1 高斯–牛顿法

牛顿法迭代公式为

其中

$S\left(x\right)=\sum _{i=1}^m{F}_i\left(x\right){\nabla }^2{F}_i\left(x\right)$

由于 $S（x）$中的残量函数 $F_i（x）$的 Hessian 矩阵 ${\nabla }^2{F}_i（x）$ 计算量较大，如果忽略这一项，便得到求解非线性最小二乘问题的高斯–牛顿迭代公式：

高斯–牛顿法仅需要计算一阶导数，以目标函数的 Jacobian 矩阵来构造 Hessian 矩阵的近似矩阵，减少了计算量，提高了运行效率。但是，若在迭代过程中$J_k^{\mathrm{T}}{J}_k$不可逆，算法将发散。另外，由于忽略了 Hessian 矩阵的二阶偏导项，因此当残差较大或方程组的非线性强度高时，高斯–牛顿法会带来较大的误差。

2.2 LM 算法

高斯–牛顿法在迭代过程中要求矩阵 $J_k^{\mathrm{T}}{J}_k$列满秩，而这一条件限制了它的应用。LM 算法^[12-13] 通过求解下述优化模型来获取搜索方向：

式中 ${\lambda }_k>0$。由最优性条件知 d_k 满足

求得

从而得到求解非线性最小二乘问题的 LM 算法迭代公式：

LM 算法的关键是在迭代过程中如何调整阻尼因子${\lambda }_k$。${\lambda }_k$过大，算法迭代速度变慢，效率低；${\lambda }_k$过小，矩阵$J_k^{\mathrm{T}}{J}_k+{\lambda }_{k}I$接近奇异，算法稳定性差，易发散。FAN 和 PAN 证明了选取合理的阻尼因子策略，可使得 LM 算法在满足局部误差界条件下，具有局部二阶收敛性，为提高磁性目标参数估计精度奠定了理论基础^[14]。为了满足阻尼因子随迭代过程变化的要求，本文选择${\lambda }_k=\mu {∥F\left(x_k\right)∥}^{\sigma }$，式中$\mu >0，\sigma \in \left[\mathrm{1，2}\right]$。由于磁性目标产生的磁异常值会随着距离的三次方衰减，当目标较远时，磁传感器测得的磁感应强度值通常比较小，此时$∥F\left(x_k\right)∥$远小于 1，因此，${\lambda }_k$关于$\sigma$ 是单调递减的。当$\sigma$=1时，则 ${\lambda }_k$较大，又由式（8），从而产生较小的步长。为了避免过小的步长导致收敛速度慢，本文选择 $\sigma$=2，使其在保证收敛的同时，进一步提高算法的迭代效率。

2.3 改进 LM 算法

LM 算法改进了$J_k^{\mathrm{T}}{J}_k$的奇异问题，但只具有局部收敛性。为了降低 LM 算法对初值的敏感性，引入信赖域技术^[15]，其基本思想是在当前迭代点 x_k 的附近用一个简单函数近似目标函数 f，并用该近似函数在 x_k 的某个领域内的极小值点作为下一个迭代点。

为了验证当前迭代步的有效性，引入取舍指标：

式中，r_k 定义了目标函数的实际减小量与预测减小量之比，反映了目标函数与近似函数的近似程度。

当取舍指标r_k大于设定的（非负）阈值时，说明近似函数与目标函数的近似程度好，接受当前迭代步 d_k；否则，不接受该迭代步。利用r_k来调节参数μ_k，若r_k小于设定的（非负）阈值，则增大 μ_k； 若r_k大于设定的（非负）阈值，则减小μ_k；否则，μ_k不变。

3 改进 LM-GN 算法

改进 LM 算法由于使用了信赖域搜索技术使得它比传统高斯–牛顿法和 LM 算法的收敛范围更广^[14]，降低了算法对初始值的依赖程度，但却减慢了收敛速度。因此，本文基于改进 LM 算法提出了一种混合算法，称为改进 LM-GN 算法，该混合算法利用改进 LM 算法不严重依赖于初始值的性能和高斯–牛顿法收敛速度快的优点，提高磁性目标参数估计的精度和效率。

算法流程图如图2 所示。

改进 LM-GN 算法通过设置一个非负阈值κ，当目标函数值大于这一阈值时，说明迭代点离真值较远，此时使用对初始值依赖程度较低的改进 LM 算法进行迭代；当迭代到一定程度，即目标函数值小于或等于这一阈值时，说明该迭代点已经在空间最优解的附近，此时使用收敛速度较快的高斯–牛顿法继续进行迭代计算，最终得到最佳的磁性目标定位结果。算法 3.1（改进 LM-GN 算法）具体步骤如下。

Step 1 给定变量初值 x₀，迭代步数 k=1，常数 m，$0<p_0<p_1<p_2<1，c_1，c_2>1$以及收敛精度 ε、阈值κ，设置${\mu }_0$满足${\mu }_0>$m；

Step 2 计算 $g_k=\nabla f\left(x_k\right)$，若$∥g_k∥\le \epsilon$，停算，输出近似解 x_k，否则，转到 Step 3；

Step 3 计算$f\left(x_k\right)$，若 $f\left(x_k\right)>\kappa$，则转 Step 4，否则，转 Step 8；

Step 4 令${\lambda }_k={\mu }_k{∥F\left(x_k\right)∥}^2$，根据式（8）计算 d_k；

Step 5 按式（10）计算 r_k，并根据 r_k 选择是否接受 d_k：

$x_{k+1}=\left\{\begin{array}{c}{x}_k+d_k{r}_k>p_0\\ x_k\text{o}\text{t}\text{h}\text{e}\text{r}\text{w}\text{i}\text{s}\text{e}\end{array}\right.$

Step 6 调整参数 ${\mu }_k$：

${\mu }_{k+1}=\left\{\begin{array}{c}{c}_1{\mu }_k{r}_k<p_1\\ {\mu }_k{p}_1<r_k<p_2\\ max\left\{\frac{{\mu }_k}{c_2},m\right\}{r}_k>p_2\end{array}\right.$

Step 7 令 $k=k+1$，转到 Step 2；

Step 8 根据 $d_k=-{\left(J_k^{\mathrm{T}}{J}_k\right)}^{-1}{g}_k$计算 d_k；

Step 9 令$x_{k+1}=x_k+d_k，k=k+1$，转 Step 2。

式中：p₀ 为迭代步取舍指标；p₁ 为增大阻尼因子时 r_k 取值的上限阈值；p₂ 为减小阻尼因子时 r_k 取值的下限阈值。为了防止序列接近解时步长过大，令 m 为参数${\mu }_k$的下限阈值。

4 仿真实验分析

本节通过仿真实验比较改进 LM-GN 算法、改进 LM 算法、LM 算法和高斯–牛顿法的磁性目标参数估计性能。实验中磁传感器测量阵列和目标的空间分布如图3 所示。

假设磁性目标的实际参数为 $\left(x_0，y_0，z_0，m_{x_0}\right.\text{，}\left.m_{y_0}，m_{z_0}\right)$，通过优化算法估计得到的结果为$（x，y，z，\left.m_x，m_y，m_z\right)$，定义定位误差为

定义磁矩误差为

算法参数设置如下：收敛精度阈值$\epsilon ={10}^{-21}$， LM 算法中的阻尼因子更新方式为${\lambda }_k=\mu {∥F\left(x_k\right)∥}^2$，其中$\mu ={10}^{-4}$，改进 LM 算法中的 $c_1=c_2=5，p_0=0.25，p_1=0.5，p_2=0.75，m={10}^{-6}，{\mu }_0={10}^{-4}$，判断阈值κ为10^-15。

目标参数的真值设为[10 10 10 800 700 800]^T，在无噪声、标准差分别为0.01 nT 和0.1 nT 的 0 均值高斯白噪声的情况下，为了分析高斯–牛顿法和 LM 算法对初始值的依赖程度，根据参数初值离真值的距离选取3组不同的参数初值分别为[8 8 8 500 500 500]^T、[1 1 1 500 500 500]^T和[0 0 1 800 700 800]^T进行仿真实验。高斯–牛顿法、LM 算法、改进 LM 算法和改进 LM-GN 算法的参数估计结果分别如表1、表2 和表3 所示。

由表1、表2 和表3 结果可以看出，在同一迭代初始值条件下，随着噪声标准差的不断增大，4 种算法的目标参数估计结果与真值的差值逐渐增大，表明噪声对磁性目标定位有较大的影响。另外，与 LM 算法相比，高斯–牛顿法更容易受到噪声的干扰，这是由于随着噪声标准差的不断增大，该优化问题由小残差问题转变为大残差问题，此时 Hessian 矩阵中的 S（x）的作用不容忽视，高斯–牛顿迭代步中舍去这一项将导致该算法发散，而 LM 算法中通过${\lambda }_{k}I$修正了 Hessian 矩阵的近似项，在保证了算法收敛性的同时也避免了$J_k^{\mathrm{T}}{J}_k$接近奇异时，搜索方向的模 $∥d_k∥$过大的问题。在同一噪声标准差不同迭代初值条件下，这 4 种方法所得解的正确性依赖于选取的迭代初值与真值的靠近程度，初值越靠近真实值，算法的估计结果越准确。对于高斯–牛顿法和 LM 算法，若其中存在几个参数远离真实值，则不能保证其迭代的收敛性和解的正确性，由表1 可知即使在无噪声的情况下这 2 种算法也会迭代失败，可见高斯–牛顿法和 LM 算法的计算结果严重依赖于初始值的选取。另外，高斯–牛顿法对初始值的依赖程度强于 LM 算法，这是因为方程组的非线性程度越高，S（x）的值越大，而高斯 –牛顿法的迭代公式舍去 S（x），导致矩阵接近奇异，在远离真实值的初始值开始得到的搜索方向不是最有效的下降方向，进而陷入初始值附近的局部最优解，导致迭代失败。在相同的噪声标准差条件下，只有在选取的初始值每个参数都接近真实值的情况下，这 2 种算法才会收敛到全局最优解。然而在磁性目标定位问题中，目标的位置和磁矩真值事先未知，所以初始值的选取具有随机性，并不能保证选取的初始值靠近真实值。而改进 LM 算法和改进 LM-GN 算法在上述不同的初始值的情况下都能迭代成功，说明引入信赖域技术之后，算法在一定程度上降低了解对初始值的依赖程度。

下面仍考虑无噪声、标准差分别为 0.01nT 和 0.1nT 的 0 均值高斯白噪声的情况，为评估不同算法对初始值的依赖程度，在解空间内随机给定初始值的情况下，分别进行 1 000 次的 Monte Carlo 实验，并将定位误差是否超过 1 m 作为定位失败或成功的评判指标。得到 4 种算法的定位正确率如表4 所示。本文选取解空间为$x\in \left[\mathrm{0，20}\right]，y\in \left[\mathrm{0，20}\right]，z\in \left[\mathrm{0，20}\right]，m_x\in \left[\mathrm{0，1000}\right]，m_y\in \left[\mathrm{0，1000}\right]，m_y\in \left[\mathrm{0，1000}\right]$。

由表4 结果可以看出：在随机给定初始值的情况下，高斯–牛顿法和 LM 算法的定位正确率都很低，特别是高斯–牛顿法，即使在无噪声的情况下，定位正确率仅有 25.4%，LM 算法也只达到了 63.1%，说明这 2 种算法严重依赖于初始值的选取。而改进 LM-GN 算法和改进 LM 算法的定位正确率相当，且高于高斯–牛顿法和 LM 算法。即使在噪声较大时，定位正确率也几乎能达到 100%。这说明引入信赖域搜索技术使得算法具有全局收敛性，降低了算法对初始值的依赖程度。随着噪声的不断增大， 4 种算法的定位正确率都在逐渐减小，这是因为求解非线性最小二乘问题时代入的 $B_{x1}，B_{y1}，B_{z1}，\cdots ，B_{xN}，B_{yN}，B_{zN}$ 是带有噪声的实际测量的磁感应强度三分量值，因此求解式（3）所得结果仍不可避免地要受到噪声的影响。

设目标参数的真值仍为[10 10 10 800 700 800]^T，算法参数选取同上，由于高斯–牛顿法与 LM 算法对初始值敏感，为了能够比较这 4 种算法的收敛速度，需要一个合理的初始值，不妨假设初始值选为[3 3 3 200 200 200]^T，在噪声标准差为 0.1 nT 的条件下，高斯–牛顿法、LM 算法、改进 LM 算法和改进 LM-GN 算法同时运行一次得到的各参数值的迭代情况如图4 所示，得到目标函数值达到同样终止条件所需迭代次数如表5 所示。

由图4 可知，改进 LM-GN 算法在迭代初期与改进 LM 算法迭代相似，但在迭代一定步数之后与改进 LM 算法有明显的不同，改进 LM-GN 算法各参数收敛到目标值的速度更快，运行时达到相同的终止条件所需要的迭代步数明显更少。这是因为改进 LM-GN 算法中设置的判断阈值起了作用，在目标函数值较大时，即迭代点离真值较远，使用的是改进 LM 算法，在一定程度上减弱了解对初始值的依赖性，直至目标函数值达到判断阈值，此时迭代点在收敛解的附近，使用高斯–牛顿法加快了算法的收敛速度。改进 LM 算法和 LM 算法迭代步中求逆部分增加了 ${\lambda }_{k}I$项，使得计算得到的步长更小，因此，在避免矩阵奇异的同时却降低了收敛速度。综上所述，由改进 LM-GN 算法得到的搜索方向更接近于目标下降的方向，迭代步数更少，收敛速度更快。

如表5 所示，在噪声为 0.1nT 时，改进 LM-GN 算法和高斯–牛顿法只需要 10 步，而改进 LM 算法需 16 步，LM 算法需要 48 步才能达到相同精度。综上所述，改进 LM-GN 算法在一定程度上比改进 LM 算法具有更快的收敛速度，达到相同的终止条件所需的迭代步数更少，算法降低解对初始值依赖程度的同时，又加快了收敛速度。

接下来，仿真实验仍在上述参数设置下，在无噪声条件下，这 4 种算法迭代得到的解与真值相同。考虑不同噪声标准差条件下 4 种算法达到充分迭代的终止条件下运行 100 次，定位误差和磁矩误差计算结果如表6 和表7 所示，运行 1 次所需时间如表8 所示。

由表6 和表7 可知，在初始值选取合适的情况下，4 种算法的迭代解达到了相似的精度，说明这 4 种优化算法在良好初始值的条件下都能够解决磁性目标的定位问题。随着噪声标准差的不断增大，算法得到的解精度逐渐降低，可见噪声对磁性目标参数估计结果影响较大。

在初始值选取合适的情况下，由表8 可知，高斯–牛顿法收敛速度最快，改进 LM 算法和 LM 算法由于步长的减小导致收敛速度变慢，由图4 可知改进 LM 算法比 LM 算法迭代步数少，但表8 却显示其收敛速度慢于 LM 算法，这是因为改进 LM 算法中需要计算 r_k 和更新 μ_k，从而导致计算时间比 LM 算法多了约 0.14%。改进 LM-GN 算法所需时间多于高斯–牛顿法，但相对于 LM 算法和改进 LM 算法，其收敛速度明显加快，表明改进 LM-GN 算法在降低解对初始值依赖程度的同时也具有较快的收敛速度。

为了分析改进 LM-GN 算法对运动目标的定位性能，假设目标 z 轴坐标不变，设为 10 m，磁矩分量不变，设为$（\mathrm{800，700，800}）$，目标由初始位置$（-8，-8）$分别在 x、y 方向上均以 1 m/s 的速度做匀速直线运动，直到终点位置$（\mathrm{15，15}）$处停止。考虑标准差为 0.1nT 的 0 均值高斯白噪声，在解空间内随机生成初始值的情况下，进行 1 000 次 Monte-Carlo 实验，计算磁性目标运动到每个点处的定位正确率，结果如图5 所示。

由图5 可知，改进 LM-GN 算法在$x，y\in [-\mathrm{3，13}]$的空间范围内都能达到较高的正确率，说明该算法降低了解对初始值的依赖程度，而在其他空间范围内的迭代成功率不够高，主要是由目标距离磁传感器较远使得信噪比降低导致的。

为了全面评价所提算法的定位性能，进一步将其与 PSO-LM^[8]算法和 PSO-GN^[7]算法进行性能比较，其中 PSO 算法的种群规模设置为 30，迭代次数设置为 10。目标参数的真值仍为[10 10 10 800 700 800]^T，表9 为在噪声标准差为 0.1nT 的条件下，3 种算法在定位误差、磁矩误差和运行时间方面的比较结果。

从表9 可以看出，改进 LM-GN 算法在定位误差和磁矩误差方面能达到与 PSO-LM算法和 PSO-GN 算法相似的精度，但在运行时间上明显少于这 2 种算法，表明了本文方法在运行效率方面的优势。

5 结束语

本文基于运用信赖域技术修正的 LM 算法提出了一种改进 LM-GN 算法解决磁性目标定位问题。针对高斯–牛顿法和 LM 算法在求解非线性最小二乘问题时严重依赖初始值的缺点，通过引入信赖域搜索技术，得到了改进 LM 算法，并根据改进 LM 算法收敛速度下降这一不足，将改进 LM 算法与高斯–牛顿法结合，提出了一种改进 LM-GN 算法。仿真结果表明：在初始值随机选取的情况下，若初始值接近真实值，则高斯–牛顿法和 LM 算法都能够收敛到全局最优解；但若选取的初始值远离真实值，则 2 种方法均迭代失败。而本文提出的改进 LM-GN 算法不仅克服了传统方法严重依赖初始值的缺点，在一定的解空间中能够有效地估计出目标的位置和磁矩，降低了对初始值的依赖程度，且计算速度快，迭代步数少，能够满足实时性需求。下一步工作将开展外场试验，并用实测数据验证本文方法的有效性。

参考文献