1. 引言

近年来，随着海洋资源的不断开发，我国海洋安全防御维护面临较为严峻的挑战 [1] ，如岛屿主权、资源争端等 [2] ，且在近海水域常常有敌方水下潜航器、舰艇等的蓄意入侵与破坏 [3] ，因此，对水下运动目标体进行定位与识别对构建海洋防卫安全体系有着重大的意义 [4] 。

当前在水下目标定位探测方面，主要采用的方法是声学探测技术，但海洋环境不断复杂化以及水声效应的影响，再加上水下目标的消声隐身效果随着科技的发展也不断优化 [5] ，仅利用声学技术对水下运动目标进行定位受到了较大的限制，往往不能取得较好的效果。因此，需要发展新型的水下目标体探测技术。

水下目标体（如舰艇等）一般是由不同的合金材料制成的 [6] ，其化学特性有一定的差异，也导致了不同材料实际在海水中的电位是不同的，互相之间会发生电偶腐蚀，产生腐蚀电流。而目标体上为防止腐蚀也会安装保护系统，从而产生了保护电流，于是在水下运动目标周围会形成一个闭合导电回路 [7] ，进而产生了一个极低频电磁信号 [8] ，这也是利用极低频电磁场对水下目标进行定位的研究基础。

在舰艇极低频电磁场信号探测方面，梁成浩等采用船体缩比模型研究信号特征，为极低频电场的检测提供了新途径 [9] 。渠晓东研究了水下目标极低频电磁信号的产生机理，推导了海水中的水平电偶极子源的电磁场计算公式，研究了极低频电磁波在海水–空气分界面上的反射、折射现象、多径效应、传播规律和参数影响 [10] 。张朔宁等以时谐垂直电偶极子作为水下目标产生极低频电磁场模型，模拟了等效电性源空中磁场分量、海水中电场分量的空间分布特征，分析了轴频电磁场的可探测范围以及目标方位对于电磁场的影响 [11] [12] 。庞鑫等对舰船轴频磁场进行了建模，并将南海试验得到实测数据与模型分析数据进行比较分析，发现两者特征基本吻合，验证了水平时谐电偶极子模型模拟轴频磁场的可行性 [13] 。陈帅等提出了一种将目标电场与磁场信息融合的探测方法，能够提高信号的探测效率 [14] 。徐震寰等采用了自适应傅里叶分解的方法，提高了提取舰船极低频电磁场的速度与精度 [15] 。岳瑞永等仿真计算了水下目标产生的磁场，研究了其特性以及传播规律，为极低频磁场的识别提供了一定的理论支撑 [16] 。而在基于极低频电磁场进行目标定位方面，李涛等采用遗传算法对目标位置进行求解，但在积分时较为复杂 [17] ，张雪峰等将粒子群优化算法与遗传算法相结合 [18] ，对位置参数进行求解，取得了较好的效果 [19] 。孙宝全等设计了渐进更新扩展卡尔曼滤波方法来通过静态电场来跟踪舰船目标 [20] 。王向军等采用了扩展卡尔曼滤波器，融合了极低频电场的幅值以及相位，为目标定位提供了新的思路 [21] 。

本文提出了一种基于水下目标极低频电磁场的目标定位方法，用水平电偶极子模型模拟水下真实目标体，引入了反演方法，通过已知的电磁场各个分量数据反演计算目标的位置信息。同时，改变了多个参数，设计了多组试验验证了方法的可行性与适用范围，取得了较好的定位效果。

2. 方法原理

2.1. 原理概述

本文主要研究的是舰船、潜艇等目标体在水下航行时，产生的以螺旋桨旋转频率为基频的极低频电磁场，其场源可以用水平电偶极子来模拟 [22] 。于是，可以假设电偶极子在层状介质中运动，来模拟海洋中真实目标体的运动情况。由于缺少实际测量的电磁场数据，因此使用Key文章中的方法来计算层状介质中水平电偶极子产生的电磁场各个分量数据 [23] 。假设一层状模型如图1所示。

用矢量位A来表示电磁场。

$\{\begin{array}{l}E=i\omega A+\nabla \left(\frac{\nabla \cdot A}{\mu \sigma }\right)\hfill \\ B=\nabla \times A\hfill \end{array}$(1)

于是可以得到各个电磁场分量的表达式 [23]

$\{\begin{array}{l}{E}_x=-\frac{\sin 2\phi }{4\pi \mu {\sigma }_i}{\displaystyle {\int }_0^{\infty }\left({\text{<mover accent="true"> A <mo> ^ </mo> </mover>}}_y\left(\lambda ,z\right)+\frac{\partial {\hat{\Lambda }}_z\left(\lambda ,z\right)}{\partial z}\right)}\left({\lambda }^3{J}_0\left(\lambda r\right)-\frac{2{\lambda }^2}{r}{J}_1\left(\lambda r\right)\right)\text{d}\lambda \\ E_y=\frac{i\omega }{2\pi }{\displaystyle {\int }_0^{\infty }{\text{<mover accent="true"> A <mo> ^ </mo> </mover>}}_y}\left(\lambda ,z\right)\lambda J_0\left(\lambda r\right)\text{d}\lambda -\frac{{\cos }^2\phi }{2\pi \mu {\sigma }_i}{\displaystyle {\int }_0^{\infty }\left({\text{<mover accent="true"> A <mo> ^ </mo> </mover>}}_y\left(\lambda ,z\right)+\frac{\partial {\hat{\Lambda }}_z\left(\lambda ,z\right)}{\partial z}\right){\lambda }^3{J}_0\left(\lambda r\right)\text{d}\lambda }\\ +\frac{\cos 2\phi }{2\pi \mu {\sigma }_i}{{\displaystyle \int }}_0^{\infty }\left({\text{<mover accent="true"> A <mo> ^ </mo> </mover>}}_y\left(\lambda ,z\right)+\frac{\partial {\text{<mover accent="true"> Λ <mo> ^ </mo> </mover>}}_z\left(\lambda ,z\right)}{\partial z}\right)\frac{{\lambda }^2}{r}{J}_1\left(\lambda r\right)\text{d}\lambda \\ E_z=-\frac{\cos \phi }{2\pi \mu {\sigma }_i}{\displaystyle {\int }_0^{\infty }\left({\lambda }^2{\hat{\Lambda }}_z\left(\lambda ,z\right)+\frac{\partial {\text{<mover accent="true"> A <mo> ^ </mo> </mover>}}_y\left(\lambda ,z\right)}{\partial z}\right){\lambda }^2{J}_1\left(\lambda r\right)\text{d}\lambda }\end{array}$(2)

$\{\begin{array}{l}{H}_x=-\frac{\cos 2\phi }{4\pi \mu }\left[{{\displaystyle \int }}_0^{\infty }{\text{<mover accent="true"> Λ <mo> ^ </mo> </mover>}}_z\left(\lambda ,z\right)\left({\lambda }^3{J}_0\left(\lambda r\right)-\frac{2{\lambda }^2}{r}{J}_1\left(\lambda r\right)\right)\text{d}\lambda \right]\\ -\frac{1}{2\pi \mu }\left[{{\displaystyle \int }}_0^{\infty }\left(\frac{1}{2}{\lambda }^2{\text{<mover accent="true"> Λ <mo> ^ </mo> </mover>}}_z\left(\lambda ,z\right)+\frac{\partial {\text{<mover accent="true"> A <mo> ^ </mo> </mover>}}_y\left(\lambda ,z\right)}{\partial z}\right)\lambda J_0\left(\lambda r\right)\text{d}\lambda \right]\\ H_y=\frac{\sin 2\phi }{4\pi \mu }{{\displaystyle \int }}_0^{\infty }\left({\text{<mover accent="true"> Λ <mo> ^ </mo> </mover>}}_z\left(\lambda ,z\right)\right)\left({\lambda }^3{J}_0\left(\lambda r\right)-\frac{2{\lambda }^2}{r}{J}_1\left(\lambda r\right)\right)\text{d}\lambda \\ H_z=-\frac{\sin \phi }{2\pi \mu }{{\displaystyle \int }}_0^{\infty }\left({\text{<mover accent="true"> A <mo> ^ </mo> </mover>}}_y\left(\lambda ,z\right)\right){\lambda }^2{J}_1\left(\lambda r\right)\text{d}\lambda \end{array}$(3)

所有电磁场分量均表示成了Hankel变换的形式。以上是场源位于源点，方向沿着y轴的场分量结果，对于任意水平方向场源，无需再次推导，可通过坐标旋转获得相应的电磁场值。

引入反演方法，在单个目标定位时，定义模型参数 $m\left[x,y,z\right]$，数据向量 $d\left[\left|E_x\right|,\left|E_y\right|,\left|E_z\right|,\left|H_x\right|,\left|H_y\right|,\left|H_z\right|\right]$，反演迭代公式 [24] 为

$m_{k+1}={\left[\mu {\partial }^{\text{T}}\partial +{\left(W{J}_k\right)}^{\text{T}}W{J}_k\right]}^{-1}{\left(W{J}_k\right)}^{\text{T}}W{\hat{d}}_k$(4)

其中J为雅可比矩阵。当拟合差达到阈值1时，终止迭代。

运动目标体的位置随着时间变化，假设它的运动轨迹为时间的函数，每一个时刻对应一个目标的位置以及观测基站采集到的观测数据，即极低频电磁场的各个分量，利用迭代公式通过电磁场的各个分量数据迭代得到目标的位置，即x、y、z坐标，对于每一个时刻重复这一步骤，最终可以得到水下目标的估计运动轨迹。

2.2. 单个运动目标体定位算例

为了验证文中方法的可行性与准确性，设计了如图2所示的理论模型。第一层为空气层，观测基站位于其与海水的分界面上，第二层为海水层，假设海水的电导率 ${\sigma }_1$为3.3 s/m，层厚为300 m，第三层为海床，电导率 ${\sigma }_2$为0.04 s/m。

利用水平电偶极子模拟的水下运动目标位于海水层中，初始位置坐标为(100 m, 100 m, 100 m)，在海水中作曲线运动，运动轨迹方程设置为式(1)。

$\{\begin{array}{l}x=50+10\times t\\ y=50+0.5\times {\left(t-5\right)}^2\\ z=100\sin \left(0.25t\right)+200\end{array}$(5)

假设其频率为1 Hz，电偶极矩的大小为100 A·m。观测基站的坐标为(1500 m, 0 m, 0 m)，采集电场与磁场共六个分量的数据，采样间隔为1 s，为了更加接近真实环境的情况，假设仪器采集到的电磁场各个分量数据受到噪声的影响，噪声水平为1%，即在计算电磁场各个分量模值时加入1%的高斯噪声。通过采集到的电磁场各个分量数据模值计算目标在各个时刻的位置坐标值，得到的最终结果见图3，图4。

图3为利用前文中介绍的方法得到的目标最终定位结果，从x-y平面与x-z平面两个角度展示了目标估计运动轨迹与真实运动轨迹的对比，共记录了31组目标的电磁场各分量数据，测量的时间间隔为1 s，对应图中31个目标位置。如图3所示，连续的实线表示预设的目标在第二层内的运动轨迹，红圈表示由各个时刻采集到的电磁场数据计算的目标体坐标，所有的坐标点组成了目标的估计轨迹。对比两条轨迹可以发现，不管是图3(a)中的x-y平面，还是图3(b)中的x-z平面，估计得到的目标位置基本与其真实轨迹相拟合，虽然由于噪声影响，个别点出现偏离真实轨迹的情况，但也都在允许的误差范围内。

为了更加直观的对比目标在各个时刻的真实坐标值与估计坐标值，图4分别展示了目标定位的绝对误差与相对误差结果，图中的星号、圆圈以及三角符号分别代表了目标体在x、y、与z方向上的误差结果，从图4(a)中可以看出，目标估计位置在x与y方向的偏差全都小于5 m，其中在x方向上的最大绝对偏差为4.01 m，在y方向上的最大绝对偏差3.31 m，且存在很多时刻的估计位置误差接近0，同时在z方向上，大多数位置的估计误差在5 m以内，仅有少部分点的误差接近10 m，最大的绝对偏差仅有9.11 m，说明文中方法是能够准确的通过目标电磁场各个分量数据计算得到目标的位置信息的。

图4(b)展示了目标估计位置的相对误差结果，可以看出目标在各个时刻，各个方向上的估计位置坐标的相对误差均小于5%，在x方向上的平均相对误差要稍大，但也仅有1.19%，这是由于目标在x方向上的位置坐标值本身就要小于其他两个方向上的坐标值，在计算平均相对误差时，分母变小，使得最终计算得到的相对误差要略大，同时，在y与z方向上的平均相对误差分别为0.69%与2.17%，其中在y方向上的平均相对误差值是最小的，也与目标在y方向上的真实位置坐标值相对较大有一定的关系，而在z方向上的平均相对误差是最大的，这是由于在z方向上的绝对误差相对较大。但从总体上来看，各个时刻估计得到的坐标位置误差都保持在一个很小的水平，同样也说明了本文方法是能够精确计算单个水平电偶极子在层状介质中的位置的。

为进一步验证方法的可行性与应用范围，改变部分参数，设计了下述多个试验。

3. 定位算例分析

3.1. 不同电偶极矩的偶极子定位结果

前文中的算例偶极子的电偶极矩设置为100 A·m，为了研究电偶极矩大小对于定位结果的影响，将目标电偶极矩设置为1000 A·m，其它参数保持不变，计算目标在层状介质中的位置坐标，得到的最终结果如图5所示。

从图5中可以看出，在偶极子的电偶极矩设置为1000 A·m后，计算得到的目标估计轨迹依然与其真实轨迹基本一致，能够很好地反映目标在层状介质中的真实运动情况，和前文中电偶极矩的大小为100 A·m的算例相对比，两者基本有一致的曲线趋势，直角坐标系下目标估计位置点基本能够落在目标实际运动轨迹上。同时，目标位置的最大偏差在x、y与z方向上分别为5.27 m、3.95 m与14.1 m，除了个别点的误差由于所加噪声影响，绝对误差相对较大以外，其他点都能保持一个较小的误差水平，其平均相对误差在各个方向上分别为1.33%、0.75%、2.47%，都小于3%，说明在电偶极矩增大之后，本文方法依然能够取得很好的效果。

为了进一步验证这个结论，在多个电偶极矩大小上都采用本文中方法对偶极子进行位置的确定，将所得的最终结果放入表格中进行对比，如表1所示。不同电偶极矩的定位结果最大绝对误差存在一定的差距，这是由于加了噪声的原因，导致有个别时刻采集到偶极子的电磁场数据与其真实值有一定的差距，进而体现在误差上可能有相对较大的值，但也仍在可以接受的误差范围内。而不同电偶极矩定位结果的平均相对误差在各个方向上的值实际上互相是接近的，x方向上的平均相对误差都在1.5%左右，y方向上的值在0.8%左右，而z方向上的误差值在2.3%左右，这也说明了不管目标的电偶极矩大小是多少，文中方法都可以计算出目标准确的位置。

3.2. 不同频率的偶极子定位结果

下面研究频率大小对于定位结果的影响，将目标频率设置为0.1 Hz，其它参数保持不变，依然使用文中方法进行计算，得到的最终目标位置结果如图6所示。

从图6中可以看出，在改变目标的频率之后，也能够计算出目标准确的位置坐标，与前文中电偶极矩的大小为1 Hz的算例一样，在直角坐标系下目标估计坐标点组成的轨迹能够很好地拟合预设定的运动轨迹。且结果的最大绝对误差在x、y与z方向上分别为10.48 m、3.94 m与15.72 m，其平均相对误差在各个方向上分别为1.79%、0.94%、2.9%，都与前文算例中的结果相差不大。说明在频率减小之后，本文方法依然能够得到目标体准确的位置信息。

在多个频率上计算的误差结果如表2所示。和前文中不同电偶极矩的算例类似，在最大绝对误差存在一定的差距，但也都在一个较小的误差范围内。而在各个方向上的平均相对误差结果中我们可以看到，虽然每组数据的目标频率不同，但其误差百分比值都较为接近，绝大多数都在3%以内，这也说明了本文方法在各个频率上都能有很好的应用，可以得到较好的结果，且适用范围较大。

3.3. 观测基站位置对定位结果的影响

前文中观测基站的位置为(1500 m, 0 m, 0 m)，最终的定位结果取得了较好的效果。为研究观测基站相对目标位置的远近对定位结果的影响，选取了多个观测基站位置，其他的参数不变，其频率为1 Hz，电偶极矩的大小为100 A·m，得到的多组结果如表3所示。

从表3中的结果可以明显看出，随着观测基站的x坐标值变大，在x方向上目标的位置坐标值误差有了明显的增大，最大绝对误差从4.01 m增大到了43.59 m，平均相对误差也从1.19%增大到了16.29%。从图7中在x方向上目标位置坐标随时间的变化可以更加明显的看出误差的变化。如图7(a)所示，当观测基站坐标为(3000 m, 0 m, 0 m)时，大部分时刻估计的目标位置坐标还是能够与目标真实的运动轨迹相拟合，当观测基站坐标为(5000 m, 0 m, 0 m)时，如图7(b)所示，有部分点明显偏离了其真实位置，而当观测基站坐标变为(8000 m, 0 m, 0 m)时，大部分点都与其真实运动轨迹有一定的偏离，且误差较大。

为了探究误差较大的原因，将观测基站位置分别设置在(10,000 m, 0 m, 0 m)与(0 m, 10,000 m, 0 m)进行了两次计算，得到的结果如图8所示。从图8中的目标估计轨迹结果图可以看出，在图8(a)中观测基站位于(10,000 m, 0 m, 0 m)处时，计算得到的目标估计x坐标值的误差较大，而在y方向与z方向上的定位结果并没有受到观测基站在x方向上的距离较远产生的影响；同理，在图8(b)中，由于观测基站位于(0 m, 10,000 m, 0 m)，在y方向上的距离目标较远，因此计算得到的y坐标值与其真实的位置有一定的偏差，可以在图中明显的看出在y方向上的误差结果要高于其他两个方向上的误差。

而产生这种误差的原因是在反演过程中求解雅可比矩阵时，各个方向是分别求取的，在某个方向上观测基站与目标距离较远时，其电磁场线要更加稀疏，当改变该方向的值时，电磁场的变化不明显，从而导致该方向的估计坐标值不够准确。而如果要解决误差较大的问题，可以采用组合观测基站的方法，在两个位置上同时采集目标产生的电磁场各个分量数据，利用采集到的两组数据来计算目标的位置坐标。同样使用上文中位置坐标为(10,000 m, 0 m, 0 m)与(0 m, 10,000 m, 0 m)的观测基站，但是同时使用这两个观测基站采集到的数据对目标位置进行计算，得到的结果如图9所示。

从图9中的结果可以看出，虽然定位使用的两个观测基站分别在x方向上与y方向上距离目标较远，但由于使用了两组观测基站采集到的数据，最终也能得到目标的准确位置，得到的目标估计运动轨迹与其真实运动情况完全吻合。其中，在各个方向上的绝对误差分别为2.65 m、0.99 m与12.32 m，而各个方向上的平均相对误差相比之前算例中仅使用一个观测基站得到的误差结果要更小，在x、y、z方向上的值分别为0.34%、0.54%、1.45%，这是由于反演过程中增加了模型参数的数量，使得最后的反演结果要更加准确。说明使用两个观测基站对目标进行定位可以避免单个观测基站出现的定位不够准确的问题，且假设存在多个观测基站的情况，只要有一个观测基站在该方向的距离没有超过最大探测距离，就能够准确地计算出在该方向上的位置坐标。

4. 结论

本文提出了一种基于水下运动目标体产生极低频电磁场的单个目标定位方法，且假设具体的模型进行了仿真试验，最后得到的计算结果证明文中方法能够较为准确的估计目标体的位置坐标，且其估计运动轨迹与其实际运动情况相吻合。为了研究各个参数对于定位结果产生的影响，设计了多组试验，分别改变了目标偶极子的电偶极矩、频率及目标与观测基站的位置，得到了不同参数的多组定位结果，从结果可以得出，目标偶极子的电偶极矩与频率的大小不会对定位结果产生影响，观测基站与目标的在某个方向上的距离越远，定位结果越不准确，但仅影响该方向上的定位误差，且误差可以通过增加观测基站的数量来进行消除。总体来看，本文方法能够较为准确地确定单个电偶极子在层状介质中的位置坐标，且电偶极子本身的特性，如电偶极矩、频率对定位结果无影响，适用范围较大，可以为水下运动目标体的定位提供一定的参考价值。

参考文献