船舶静电场来源于船舶的腐蚀防腐电流，主要包括静态电场信号(static electric field，SE)和静态电位信号(underwater electric potential，UEP)，船舶静电场具有明显的分布特征^[1-2]，包含着船舶的位置信息，且受环境的影响相对较小，可以用来对船舶进行跟踪定位，作为声学跟踪的一个有效补充。但是对UEP信号的测量需要参考电极，参考电极要求离测量电极有一定的距离，这为实际中的电位传感器布设增加了难度。SE信号的测量需要三轴电场传感器，对制作工艺要求较高，特别是三轴传感器的正交性及与姿态传感器的对准。利用电位差(difference of UEP，DUEP)DUEP对舰船进行跟踪，即采用多组电位传感器组，每一组传感器以其中一个电位传感器的UEP信号为基准，该组中其余电位传感器的UEP信号与其做差，将差值作为观测信号。DUEP信号直接将传感器组中的一个电极作为基准电极，相比UEP观测信号，省去了参考电极的负担；相比SE观测信号，DUEP信号对传感器要求相对简单，电场信号的定义为E=dU/dl，可见DUEP信号与SE信号是很接近的，但是利用DUEP作为观测信号，电位传感器阵列可以有更自由的组合方式。

另外，针对实际电场跟踪问题，有2个问题需要考虑：1)对于船舶电场跟踪的应用场合，水下武器系统而言，大部分系统平时是处于睡眠状态，当值更引信检测到目标时，才会引发各系统动作。电场跟踪系统也是一样，平时处于睡眠状态，当确认目标出现以后，跟踪系统开始工作，对目标进行跟踪。因此，船舶电场跟踪的第一步工作是对船舶目标的检测，文献[3]提出了一种基于功率谱的船舶静电场实时检测算法，本质上是通过设定的浮动阈值进行分段功率谱进行滑动检测，本文将引入此方法将其与跟踪检测算法相融合。2)所有船舶电场可利用的观测信号都面临一个先验信息缺失的问题，为此本文引入滤波器组的方法，以卡尔曼滤波为代表的滤波估值算法可以根据传感器量测信息对目标的相关参数进行实时的估计，在目标跟踪领域获得越来越广泛的研究和应用^[4-11]。并以渐进扩展卡尔曼滤波(progressive extended kalman filter，PEKF)为基本滤波单元，以解决基于DUEP信号的船舶跟踪中先验信息缺失的问题。

1 船舶跟踪问题 1.1 状态空间模型

船舶状态空间模型的形式为^[12]：

$ \left\{ {\begin{array}{*{20}{l}} {{\mathit{\boldsymbol{x}}_k} = \mathit{\boldsymbol{a}}({\mathit{\boldsymbol{x}}_{k - 1}}) + {\mathit{\boldsymbol{w}}_{k - 1}}}\\ {{\mathit{\boldsymbol{y}}_k} = \mathit{\boldsymbol{h}}({\mathit{\boldsymbol{x}}_k}) + {\mathit{\boldsymbol{v}}_k}} \end{array}} \right. $

(1)

式中：x_k∈Rⁿ为k时刻的n维状态向量；y_k∈R^m为k时刻m维观测向量；状态转移函数a:Rⁿ→Rⁿ；观测函数h:Rⁿ→R^m；w_k-1∈Rⁿ为状态噪声向量，v_k∈R^m为观测噪声向量，且w_k-1~N(0, Q_k-1)，v_k~N(0, R_k)，Q_k-1为过程噪声协方差矩阵，R_k为观测噪声协方差矩阵，{w_k-1}和{v_k}不相关。

1.1.1 观测方程

在船舶静电场的各种建模方法中，点电流法建模精度高，更加贴近船舶实际^[13]。

依据点电流法，船舶可以等效为N个点电流等间距构成的均匀线阵模型(简称点电流阵列)来描述。阵列中每个点电流的电流密度为I_pi，间隔为l_d。海水中某场点P_f处的电位信号U可等效为N个点电流产生的电位叠加。

$ U = \frac{1}{{4\pi \sigma }}\sum\limits_{i = 1}^N {{I_{{\rm{p}}i}}} K({I_{{\rm{p}}i}},{P_{\rm{f}}}) $

(2)

式中：K(I_pi, P_f)为点电流I_pi和场点P_f之间的距离函数，在空气—海水—海床3层均匀介质(如图 1所示)条件下为：

$ \begin{array}{l} K({I_{{\rm{p}}i}},P) = \frac{1}{{\sqrt {{r^2} + {{(z - {z_i})}^2}} }} + \frac{1}{{\sqrt {{r^2} + {{(z + {z_i} - 2h)}^2}} }} + \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \sum\limits_{m = 1}^\infty {{k^m}} \left[ {\frac{1}{{\sqrt {{r^2} + {{(z + {z_i} - 2h + 2mH)}^2}} }} + } \right.\\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \frac{1}{{\sqrt {{r^2} + {{(z + {z_i} - 2h - 2mH)}^2}} }} + \\ \left. {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \frac{1}{{\sqrt {{r^2} + {{(z - {z_i} + 2mH)}^2}} }} + \frac{1}{{\sqrt {{r^2} + {{(z - {z_i} - 2mH)}^2}} }}} \right] \end{array} $

(3)

式中：H为海水深度；h为传感器深度；k=(σ-σ₁)/(σ+σ₁)为海底反射系数；σ为海水电导率；σ₁为海床电导率；m为反射层数，实际计算中其上限值可取10~20，r²=(x-x_i)²+(y-y_i)²。

对一个电位传感器组，以其中一个电位传感器为基准，其电位为U₀，其余电位传感器的信号为U_j。则在第j(j=1, 2, …, J)个电位传感器量测到的目标信号可建模如下：

$ {y^{(j)}} = {U_j} - {U_0} + v_k^{(j)} - v_k^{(0)} $

(4)

式中$v_k^{(j)}$为电场传感器的背景噪声。进一步可得该电位传感器组的观测模型为：

$ {\mathit{\boldsymbol{y}}_k} = \mathit{\boldsymbol{h}}({\mathit{\boldsymbol{x}}_k}) + {\mathit{\boldsymbol{v}}_k} = \left[ {\begin{array}{*{20}{c}} {{h^{(1)}}({\mathit{\boldsymbol{x}}_k})}\\ {{h^{(2)}}({\mathit{\boldsymbol{x}}_k})}\\ \vdots \\ {{h^{(J)}}({\mathit{\boldsymbol{x}}_k})} \end{array}} \right] + \left[ {\begin{array}{*{20}{c}} {v_k^{(1)}}\\ {v_k^{(2)}}\\ \vdots \\ {v_k^{(J)}} \end{array}} \right] $

(5)

1.1.2 状态方程

船舶目标的状态量x_k包含目标位置、速度、点电流阵列：

$ {\mathit{\boldsymbol{x}}_k} = {[{\mathit{\boldsymbol{r}}_k},{\mathit{\boldsymbol{V}}_k},{\mathit{\boldsymbol{I}}_p}]^{\rm{T}}} $

(6)

式中：r_k=[x, y, z]^T；V_k=[V_x, V_y]^T(忽略z方向上运动)；I_p=[I_p1, I_p2, …, I_pN]^T。

船舶目标具有低动态特性，其运动状态$\mathit{\boldsymbol{x}}_k^{(r,V)}$=[r_k, V_k]，可简单的建模为离散白噪声加速度(discrete white noise acceleration, DWNA)模型^[12]，其运动状态方程为：

$ \mathit{\boldsymbol{x}}_k^{(r,V)} = {\mathit{\boldsymbol{\varPhi }}_m}\mathit{\boldsymbol{x}}_{k - 1}^{(r,V)} + {\mathit{\boldsymbol{\varGamma }}_m}\mathit{\boldsymbol{w}}_{k - 1}^{(r,V)} $

(7)

式中：Φ_m为运动状态转移矩阵；Γ_m为噪声增益矩阵。

$ {\mathit{\boldsymbol{\varPhi }}_m} = \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{I}}_{3 \times 3}}}&\begin{array}{l} {T_s}{\mathit{\boldsymbol{I}}_{2 \times 2}}\\ {{\bf{0}}_{1 \times 2}} \end{array}\\ {{{\bf{0}}_{2 \times 3}}}&{{\mathit{\boldsymbol{I}}_{2 \times 2}}} \end{array}} \right];{\mathit{\boldsymbol{\varGamma }}_m} = \left[ {\begin{array}{*{20}{c}} {\frac{{T_s^2}}{2}{\mathit{\boldsymbol{I}}_{2 \times 2}}}\\ {{{\bf{0}}_{1 \times 2}}}\\ {{T_s}{\mathit{\boldsymbol{I}}_{2 \times 2}}} \end{array}} \right] $

DWNA模型假设加速度为高斯白噪声，即：

$ \mathit{\boldsymbol{ \boldsymbol{\varSigma} }} = {\rm{E}}[{\mathit{\boldsymbol{w}}_k}\mathit{\boldsymbol{w}}_k^{\rm{T}}] = {\rm{diag}} (\sigma _x^2,\sigma _\gamma ^2) $

(8)

式中σ_x、σ_y分别为x、y方向上的加速度噪声强度。

点电流阵列强度量值为常量，加一个比较小的噪声即可。则全状态方程为：

$ {\mathit{\boldsymbol{x}}_k} = \mathit{\boldsymbol{F}}{\mathit{\boldsymbol{x}}_{k - 1}} + {\mathit{\boldsymbol{w}}_{k - 1}} $

(9)

式中：$\mathit{\boldsymbol{F}} = \left[ \begin{array}{l}{\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_m}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\bf{0}_{5 \times N}}\\{\bf{0}_{N \times 5}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\mathit{\boldsymbol{I}}_{N \times N}}{\kern 1pt} \end{array} \right],\mathit{\boldsymbol{Q}} = \left[ \begin{array}{c}{\mathit{\boldsymbol{ \boldsymbol{\varGamma} }}_m}\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}{\mathit{\boldsymbol{ \boldsymbol{\varGamma} }}_m}\\\alpha {\mathit{\boldsymbol{I}}_{N \times N}}\end{array} \right]$为噪声协方差矩阵；α为一个很小的常数。

2 基于滤波器组的滤波跟踪方法 2.1 基于最大似然选择法的滤波器组基本原理

船舶电场跟踪中，滤波估值方法对方位信息最为敏感，也是决定跟踪效果的关键因素。滤波器组的基本原理是在不知目标真实方位的情况下，假定目标从多个不同方位而来，同时起始多个滤波器，即设定多个滤波器组的初值：

$ \left\{ {\begin{array}{*{20}{l}} {\mathit{\boldsymbol{\hat x}}_0^{(j)} = {{[{{(\mathit{\boldsymbol{\hat r}}_0^{(j)})}^{\rm{T}}},{{(\mathit{\boldsymbol{v}}_0^{(j)})}^{\rm{T}}},\mathit{\boldsymbol{I}}_p^{(j)}]}^{\rm{T}}}}\\ {\mathit{\boldsymbol{\hat r}}_0^{(j)} = {{[{{\mathit{\boldsymbol{\hat r}}}_{xy,0}}{\rm{cos}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\psi _j},{{\mathit{\boldsymbol{\hat r}}}_{xy,0}}{\rm{sin}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\psi _j},z]}^{\rm{T}}}}\\ {\mathit{\boldsymbol{v}}_0^{(i)} = {{[ v{\rm{cos}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\psi _j},v{\rm{sin}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\psi _j}]}^{\rm{T}}}} \end{array}} \right. $

(10)

式中ψ_j为假定的目标方位角。

假定多个初值后，利用最大似然选择法，选择出最接近真实目标轨迹的跟踪轨迹。

由上述分析可知：

$ \left\{ {\begin{array}{*{20}{l}} {{x_0} \in \{ x_0^{(j)}\} _{j = 1}^N}\\ {x_0^{(j)}\backsim N(\hat x_0^{(j)},P_0^{(j)})}\\ {P(x_0^{(j)}) = \mu _0^{(j)} = \frac{1}{N}} \end{array}} \right. $

(11)

式中：x₀为真实目标初值；$x_0^{(j)}$为第j个初值假设；$\hat x_0^{(j)}$与$P_0^{(j)}$分别为第j个假设初值的均值和协方差矩阵；$x_0^{(j)}$为真实初值的先验概率。

根据贝叶斯定理，在给定当前观测值序列y_1:k条件下，候选初值$x_0^{(j)}$为真实初值的后验概率为：

$ \begin{array}{l} \mu _k^{(j)} = p(x_0^{(j)}|{y_{1:k}}) = \frac{{p({y_k}|x_0^{(j)},{y_{1:k - 1}})\mu _{k - 1}^{(j)}}}{{\sum\limits_{i = 1}^N p ({y_k}|x_0^{(i)},{y_{1:k - 1}})\mu _{k - 1}^{(i)}}} = \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \frac{{L_k^{(j)}\mu _{k - 1}^{(j)}}}{{\sum\limits_{i = 1}^N {L_k^{(i)}} \mu _{k - 1}^{(i)}}} = \frac{{\prod\limits_{i = 1}^k {L_i^{(j)}} }}{{\sum\limits_{i = 1}^N {\prod\limits_{m = 1}^k {L_m^{(i)}} } }},j = 1,2, \cdots ,N \end{array} $

(12)

式中：p(y_k|$x_0^{(j)}$, y_1:k-1)为初值假设$x_0^{(j)}$的似然函数，记L_j(k)=p(y_k|$x_0^{(j)}$, y_1:k-1)，根据高斯假设L_j(k)~N($\mathit{\boldsymbol{z}}_k^{(j)}$; 0，$\mathit{\boldsymbol{S}}_k^{(j)}$)，其中$\mathit{\boldsymbol{z}}_k^{(j)}$与$\mathit{\boldsymbol{S}}_k^{(j)}$分别为k时刻假设初值$x_0^{(j)}$对应滤波器的新息及其协方差矩阵。

由于k时刻的归一化操作不会改变分子项$\prod\nolimits_{i = 1}^k {L_j^{(i)}} $，j=1, 2，…, N的相对大小，因此最大似然选择算法可直接利用似然函数进行判别，故对似然函数$L_k^{(j)}$取对数：

$ \begin{array}{*{20}{c}} {{\rm{ln}}L_k^{(j)} = {\rm{ln}}N(\mathit{\boldsymbol{z}}_k^{(j)};0,\mathit{\boldsymbol{S}}_k^{(j)}) = }\\ { - \frac{1}{2}(m{\rm{ln}}2\pi + {\rm{ln}}|S_k^{(j)}| + {{(\mathit{\boldsymbol{z}}_k^{(j)})}^{\rm{T}}}{{(S_k^{(j)})}^{ - 1}}\mathit{\boldsymbol{z}}_k^{(j)})} \end{array} $

(13)

进而可得对数似然函数$\tilde \mu _k^{(j)}$为：

$ \tilde \mu _k^{(j)} = {\rm{ln}}\prod\limits_{i = 1}^k {L_i^{(j)}} = \sum\limits_{i = 1}^k {{\rm{ln}}} L_i^{(j)} $

(14)

利用各模型中的最大似然值进行归一化，即：

$ \tilde \mu _k^{(j)} \Leftarrow \tilde \mu _k^{(j)} - \mathop {{\rm{max}}}\limits_j (\tilde \mu _k^{(j)}) $

(15)

2.2 渐进扩展卡尔曼滤波

PEKF相比于传统卡尔曼滤波有更好的收敛性和稳定性^[14-15]，其基本步骤为：

1) 时间更新。

$ \left\{ {\begin{array}{*{20}{l}} {{\mathit{\boldsymbol{P}}_{k|k - 1}} = \mathit{\boldsymbol{F}}{\mathit{\boldsymbol{P}}_{k - 1|k - 1}}{\mathit{\boldsymbol{F}}^{\rm{T}}} + {\mathit{\boldsymbol{Q}}_{k - 1}}}\\ {{{\mathit{\boldsymbol{\hat x}}}_{k|k - 1}} = \mathit{\boldsymbol{F}}{{\mathit{\boldsymbol{\hat x}}}_{k - 1|k - 1}}} \end{array}} \right. $

式中：${{\mathit{\boldsymbol{\hat x}}}_{k|k - 1}}$为线性化展开点；H_k为h(x_k)在${{\mathit{\boldsymbol{\hat x}}}_{k|k - 1}}$的Jacobian矩阵。

2) 观测更新。

① 初始化${\left. {{{\hat x}_0}} \right|_{i = 1}} = {{\hat x}_{k|k - 1}},{\left. {{\mathit{\boldsymbol{P}}_0}} \right|_{i = 1}} = {\mathit{\boldsymbol{P}}_{k|k - 1}}$

② For i=1, 2, …, N_pu

$ {{\mathit{\boldsymbol{S}}_{i - 1}} = {\mathit{\boldsymbol{H}}_{i - 1}}{\mathit{\boldsymbol{P}}_{i - 1}}\mathit{\boldsymbol{H}}_{i - 1}^{\rm{T}} + {\mathit{\boldsymbol{R}}_k}/\delta \lambda } $

$ {{\mathit{\boldsymbol{K}}_i} = {\mathit{\boldsymbol{P}}_{i - 1}}\mathit{\boldsymbol{H}}_{i - 1}^{\rm{T}}\mathit{\boldsymbol{S}}_{i - 1}^{ - 1}} $

$ {{{\mathit{\boldsymbol{\hat x}}}_i} = {{\mathit{\boldsymbol{\hat x}}}_{i - 1}} + {\mathit{\boldsymbol{K}}_i}({\mathit{\boldsymbol{y}}_k} - \mathit{\boldsymbol{h}}({{\hat x}_{i - 1}}))} $

$ {{\mathit{\boldsymbol{P}}_i} = (\mathit{\boldsymbol{I}} - {\mathit{\boldsymbol{K}}_i}{\mathit{\boldsymbol{H}}_{i - 1}}){\mathit{\boldsymbol{P}}_{i - 1}}} $

③ ${{\mathit{\boldsymbol{\hat x}}}_{k|k}} = {\left. {{{\mathit{\boldsymbol{\hat x}}}_1}} \right|_{i = n}},{\mathit{\boldsymbol{P}}_{k|k}} = {\left. {{\mathit{\boldsymbol{P}}_1}} \right|_{i = n}}$

当n=1，PEKF退化为EKF。

综上所述，可给出FBS-PEKF跟踪算法步骤如下:

1) 确定假设目标距离数K及每个距离上起始的滤波器数量N_F，确定滤波器的初始条件$\mathit{\boldsymbol{\hat x}}_0^{(i,j)}$, $\mathit{\boldsymbol{P}}_0^{(i,j)}$，i=1, 2, …K，j=1, 2, …, N_F;

2) 从k=1时刻开始，对每个初值模型进行PUEKF滤波，得到各个滤波器的输出, 即{$\mathit{\boldsymbol{\hat x}}_k^{(i,j)},\mathit{\boldsymbol{P}}_k^{(i,j)},z_{i,j}^{(k)},S_{i,j}^{(k)}$}；

3) 利用式(14)更新k-1时刻权值，得到当前时刻各初值模型的似然值$\tilde \mu _{i,j}^{(k)}$。

4) 输出k时刻最大似然选择的滤波器结果$\mathit{\boldsymbol{\hat x}}_0^{({i_{\rm \max }},{j_{\rm \max }})},\mathit{\boldsymbol{P}}_0^{({i_{\rm \max }},{j_{\rm \max }})}$，且

$ \{ {i_{{\rm{max}}}},{j_{{\rm{max}}}}\} = \mathop {{\rm{max}}}\limits_{i,j} (\tilde \mu _k^{(i,j)}) $

(16)

5) 对似然值进行归一化。

6) 回到步骤2)或终止滤波。

3 基于功率谱的静电场检测技术

舰船静电场能量主要集中在DC~0.1 Hz的低频段。海洋环境电场主要是由于海浪运动切割地磁场产生的感应电场，其能量主要集中在0.1 Hz以下，在短时间(15~30 min)内，海浪的变化可视为平稳随机过程。即海洋环境电场的能量在短时间内不会发生明显变化，而当舰船出现时，目标周围电场将引起低频频段内的能量发生改变，这为以低频频段内的能量和作为特征检测舰船电场提供了可能。静电场检测算法的基本步骤为：

1) 计算特征值。假设检测开始时刻为t，对检测时刻t以及之前的N-1个数据{x(n), x(n-1), …, x(n-N+1)}进行功率谱估计，取特征量为谱估计后所得序列中与特征频段对应的M点频谱值的平均：

$ P = \frac{1}{M}\sum\limits_{i = 1}^M {{P_k}} $

(17)

2) 确定阈值。取检测时刻t之前某段时间特征量的平均值为阈值，计算检测时刻t之前w个特征量的平均值P_Ave=(P_n-1+P_n-2+…+P_n-w)/w，取检测阈值为T₀=βP_Ave，其中β为阈值因子，β的取值可根据要求的虚警概率和检测概率确定。

3) 判决条件。当t时刻的特征量P_t>T₀时认为满足判别条件，即目标出现一次，若不满足条件则回到步骤1)，对t+1时刻继续进行滑动检测。

4 仿真结果及分析

设计如图 2所示的仿真实验，采用9个点电源模拟舰船，每个电场传感器组包含4个电位传感器。模拟舰船的点电源阵列为I_p=[4.07, -0.68, …, -1.39]^T(由某型船舶Beasy数据利用点电源建模法反演得到)，l_d=15 m；以传感器阵列为基础建立右手坐标系，传感器所在线为Y轴，传感器中心点为坐标原点，X轴过原点垂直于Y轴，Z轴向上；电位传感器组的位置为o₁=[0 48 0;0 50 0;0.8 49 0;-0.8 49 0]，o₂=[0-48 0;0-50 0;0.8-49 0;-0.8 -49 0]；舰船起始位置r₀=[1 040，-30，-4.57]^T，初始速度为V₀=[-10, -2]^T；过程加速度噪声强度σ_x=σ_y=0.001，α=0.001；观测噪声协方差矩阵R_k=σ²I_2×2，σ=1×10^-9；传感器深度H_sensor=50 m；海水深度H_sea=51 m；海水电导率σ₁=4 S/m；海水电导率σ₁=0.04 S/m。

滤波算法中假设目标点电源数量为3，点电源初值为[-10, 30, -20]，l_d=33 m。对功率谱滑动检测器进行配置：信号采样率100 Hz，滑动窗长度为500，检测门限窗数为4，滑动间隔为1 s，连续3 s检测到目标则认为目标存在。对滤波器组进行配置：信号采样率1 Hz；检测到目标后假定目标水面初始距离为900 m；在该距离上起始N_F=4个PEKF滤波器，各个滤波器的假设方位分别为π/4、3π/4、5π/4、7π/4。

选择位置分量在k时刻的均方根误差(root mean square errors at k moment of r position component，${\rm RMSE}_k^r$)作为衡量算法优劣的标准:

$ {\rm{ RMSE}} _k^r = \sqrt {\frac{1}{{MC}}\sum\limits_{i = 1}^{MC} {{\rm{ tr }}} [({r_k} - {{\hat r}_k}){{({r_k} - {{\hat r}_k})}^{\rm{T}}}} $

(18)

如图 3~图 6所示滤波器组的跟踪结果。其中图 3是阈值——特征值曲线；图 4是滤波器组中4个滤波器的${\rm RMSE}_k^r$；图 5是滤波器组的对数似然函数值曲线；图 6是基于最大似然法选择出的最优的跟踪轨迹。

从跟踪结果中可以得到以下几点结论：

1) 基于DUEP对船舶目标进行跟踪是可行的，且具有较高的精度；

2) 滑动检测算法有效，能够及时的检测到目标，并起始滤波(图 3)；

3) 最大似然选择法能够选择正确的跟踪轨迹，对数似然函数值最大的曲线对应的为轨迹1(图 5、图 6)；

4) 从跟踪结果来看，只有轨迹1是收敛的，表明基于DUEP的船舶电场跟踪也是对方位信息最为敏感(图 4、图 6)；

5) 由于船舶在纵向和横向的尺度不同，因此对其在纵向和横向的跟踪精度要求是不同的，从图 6中可以看出，在船舶横向上，跟踪精度是很高的，${\rm RMSE}_k^r$误差只要来源于纵向误差。这个是符合实际需求的。

5 结论

1) 利用船舶DUEP信号可有效实现船舶的信号检测和跟踪，相比于UEP和SE，DUEP量级更大抗干扰性更强，对方位信息最为敏感，且对传感器的布设要求低。

2) 基于最大似然选择法的滤波器组方法可有效解决静态电位差跟踪中初始位置、速度等先验信息缺失问题。基于静态电位差的船舶电场跟踪方法是可行的，且利用工程实现，对利用船舶电场信号对船舶进行跟踪具有重要的意义。