被动声呐通过接收并处理目标辐射源信号获得目标的方位信息^[1]，具有隐蔽性强的优点，是解决水下探测问题的重要设备。而利用不同时刻观测到的方位信息，确定目标位置和航迹的过程，称为纯方位运动目标分析^[2]。由于水下环境复杂，噪声较大，且由被动声呐获得的目标信息较少，在作战时的被跟踪目标常通过机动来提高生存能力，而粒子滤波算法适用于非线性系统的状态估计、对噪声的高斯性依赖不强，因此常用于纯方位水下目标跟踪^[3]。

粒子滤波算法用一批带权值粒子的加权和来估计目标当前时刻的状态。在多次迭代之后，大部分粒子的权值会趋近于零，其余部分粒子的权值增大到趋近于1，造成粒子退化的现象。解决该现象的一种方式是重采样，其核心在于复制权值大的粒子，抛弃权值小的粒子，它能大大提升有效粒子的数量，但会导致粒子贫化的问题。遗传算法在选择、杂交、变异等算子的作用下，遵循适者生存的原理模拟生物界的进化过程，具有全局寻优的能力^[4]，在重采样时加入遗传算法能在提升有效粒子数量的前提下增加粒子的多样性。已有学者对遗传粒子滤波做了一些研究。李秀英^[5]在重采样时引入杂交和突变避免了传统重采样导致的粒子枯竭，使定位精度得到提升。文献[6]在智能城市运动目标跟踪时，将量子遗传算法与粒子滤波算法相结合，在对粒子进行二进制编码时使用了量子比特和量子叠加态，在选择、交叉、变异等操作时使用了量子旋转门。于春娣等^[7]在粒子变异时将高权值父粒子繁殖的子粒子对称地分布在父粒子两侧，从而增加粒子多样性，保证算法精度。戴银娟等^[8]在重采样时将权值方差低于阈值的粒子进行交叉和变异操作，提高了预测结果的精度。文献[9]利用基因突变和重构粒子集的方法，设计了一种优化机制，以改善用于融合定位的粒子滤波算法中的粒子退化问题。

针对纯方位机动目标跟踪问题，若遗传重采样时由大权值父粒子变异产生的新生代粒子分布方向单一、变异幅度过小或过大，会导致目标机动后接近目标真实状态的粒子少、跟踪误差大和收敛时间长的问题，本文提出一种在进化时调整新生代粒子的分布方式及变异的幅度，实现对纯方位机动目标更精确跟踪的改进遗传粒子滤波算法。

1 基于粒子滤波的纯方位目标跟踪 1.1 纯方位运动目标跟踪模型

本文研究的是单个运动观测站，假定它与水下目标处于同一平面，跟踪态势可简化为二维平面^[10]，被动跟踪系统在k时刻的状态方程和观测方程分别如下：

$ {{\boldsymbol{x}}_k} = {\boldsymbol{\phi} _{k/k - 1}}{{\boldsymbol{x}}_{k - 1}} + \boldsymbol{\varGamma} {{\boldsymbol{w}}_{k - 1}} $

(1)

$ {\hat z_k} = a\tan \left( {\frac{{{x_k} - {x_{wk}}}}{{{y_k} - {y_{wk}}}}} \right) + {v_k} $

(2)

式中： $ {\boldsymbol{x}} = {\left\{ {x,y,{v_x},{v_y}} \right\}^{\rm{T}}} $ 和 $ \hat z $ 为目标的状态向量和方位角测量值； $ \left( {{v_x},{v_y}} \right) $ 为目标速度； $ \left( {x,y} \right) $ 和 $ \left( {{x_w},{y_w}} \right) $ 为目标和观测站所在位置；w和v分别为过程演化噪声和观测噪声；a为常数；状态转移阵Φ和噪声转移阵Г分别为

$ \begin{aligned} &\;\; \boldsymbol\varPhi = \left( {\begin{array}{*{20}{c}} 1&0&T&0 \\ 0&1&0&T \\ 0&0&1&0 \\ 0&0&0&1 \end{array}} \right) \\ &\boldsymbol\varGamma = \left( {\begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} {{T^2}/2} \\ 0 \\ T \\ 0 \end{array}}&{\begin{array}{*{20}{c}} 0 \\ {{T^2}/2} \\ 0 \\ T \end{array}} \end{array}} \right) \end{aligned} $

式中T是观测周期。水下纯方位运动目标跟踪就是根据被动声呐接收到的带噪声的方位角序列值来估计目标当前状态的过程。

1.2 跟踪流程

1)初始化。根据先验信息给出目标初始状态以及过程演化噪声 $ {{\boldsymbol{w}}_0} $ 和观测噪声 $ {v_0} $ ，按照高斯分布在初始状态附近采样N个粒子 $ {\boldsymbol{x}}_0^i \sim p\left( {{{\boldsymbol{x}}_0}} \right) $ ，i=1, 2, …, N，粒子的权值 $ w_0^i = 1/N $ ；

2)在k时刻，按照式(1)的状态转换方程得到粒子 $ {\boldsymbol{x}}_k^i $ ；

3)用式(2)的测量方程计算k时刻各个粒子的方位角估计值 $ \hat z_k^i $ ，并用正态分布概率密度函数计算方位角估计值与实际值之间的似然估计 $ \hat p\left( {{z_k}\left| {\hat z_k^i} \right.} \right) $ ：

$ \hat p\left( {{z_k}\left| {\hat z_k^i} \right.} \right) = \frac{1}{{\sqrt {2{\text π}{v_k}} }}{{\rm{e}} ^{ - \frac{{{{\left( {{z_k} - \hat z_k^i} \right)}^2}}}{{2{{\left( {{v_k}} \right)}^2}}}}} $

4)计算粒子的重要性权重 $ w{_k^{i^*}} $ ，归一化后得到权重 $ w_k^i $ ：

$ w{_k^{i^*}} = w_{_{k - 1}}^i \cdot \hat p\left( {{z_k}\left| {\hat z_k^i} \right.} \right) $

$ w_k^i = \frac{{w{{_k^i}^*}}}{{\sum\limits_{i = 1}^N {w{{_k^i}^*}} }} $

5)计算有效采样大小的估计值 $ {\hat N_{\rm{eff}}} $ ，其数值越小表明退化越严重。当 $ {\hat N_{\rm{eff}}} $ < $ {N_{\rm{th}}} $ 时，进行重采样：

$ {\hat N_{\rm{eff}}} = \frac{1}{{\sum\limits_{i = 1}^N {{{\left( {w_k^i} \right)}^2}} }} $

6)计算k时刻目标状态的估计值 $ {\hat {\boldsymbol{x}}_k} $ ：

$ {\hat {\boldsymbol{x}}_k} = \sum\limits_{i = 1}^N {w_k^i} {\boldsymbol{x}}_k^i $

7)令k=k+1，重复步骤2)~7)。

步骤5)中的重采样是防止粒子退化的一个重要的手段，常用的有系统重采样、分层重采样、多项式重采样和留数重采样等^[11]，它们的宗旨是复制权值大的粒子，抛弃权值小的粒子。大权值粒子的过度复制和小权值粒子的简单剔除，是重采样算法中粒子多样性丧失的主要原因^[12]。粒子集中在某一时刻目标的高似然区域，若目标持续机动一段时间，接近目标真实状态的粒子数将很少，跟踪的收敛时间将较长。遗传算法在选择、杂交、变异等算子的作用下，遵循适者生存的原理模拟生物界的进化过程，具有全局寻优的能力，在重采样时加入遗传算法能在提升有效粒子数量的前提下增加粒子的多样性^[13]。

2 改进遗传粒子滤波

本文改进遗传重采样的基本思想，按照粒子的适应度选择出较适合环境的粒子，选择完成后进化，其中，粒子的适应度由其权值决定，在进化时将新生代粒子随机分配到以父粒子为中心均匀划分的各个扇面中，而非单一地划分到父粒子的45°或235°方向，父粒子进化得到的新生代粒子的个数与扇面的个数相差1，从而扩大粒子分布范围，提高新生代粒子的质量；由于变异幅度过小会导致新生代粒子与父粒子的差异较小，无法保证粒子的多样性，变异幅度过大会使遗传算法沦为单纯的随机选择算法，且水下目标位置和速度的变化幅度相差较大，因此本文将扩散因子 $ \gamma $ 拆分成位置扩散因子 $ {\gamma _{\boldsymbol{x}}} $ 和速度扩散因子 $ {\gamma _{\boldsymbol{v}}} $ ，分别对它们进行取值。当前时刻的目标状态由父粒子和新生代粒子共同估计，并将所有粒子用于下一时刻的迭代。

2.1 粒子的选择

选择粒子的原则是适应度越高的粒子能进化出越多的新生代粒子，反之则越少，适应度小于阈值的粒子将直接被抛弃。本文借鉴文献[14]的简单重采样法中求粒子分解数的方法计算进化的新生代粒子数量，即用式(3)比较父粒子的权值与1/N的整数倍，将倍数 $ n_k^i $ 作为k时刻第i个粒子的新生代粒子数。

$ n_k^i = \left[ {w_k^i \cdot N + 0.5} \right] $

(3)

为了保证进化后粒子总数保持不变且权值和为1，需要根据前一个粒子的进化情况对当前粒子的权值进行少量(小于1/2N)的调整。

$ \overline w _k^i = w_k^i + \left( {w_{k - 1}^{i - 1} - n_k^{i - 1}/N} \right) $

2.2 粒子的进化

传统的变异算子如式(4)，其中 $ \rho $ 是步长。经过该变异后，新生代粒子分布在以父粒子为中心、 $ \rho $ 为最大扩散半径的45°和235°方向内，分布方向单一、分布区域较集中导致新生代粒子的涵盖面积较小。

$ \overline x=x+ρ\cdot N(0,1) $

式中N(0,1)是期望为0、方差为1的标准正态分布。

图1为粒子分布方式图。图1(a)显示了一个父粒子按传统变异算子进化生成30个子粒子的结果。由于该分布方式的新生代粒子的分布区域较集中，在目标机动时可能会导致缺少与目标真实状态接近的粒子，需要花费较长时间才能跟踪上目标。

为了使新生代粒子均匀分布在父粒子周围，扩大粒子分布范围，提高新生代粒子的质量，本文对变异算子做了改进，具体方式为：

1)若 $ n_k^i = 0 $ ，则抛弃该粒子。

2)若 $ n_k^i = $ 1，则复制该父粒子1次。

3)若 $ n_k^i \geqslant 2 $ 时，先复制1次该父粒子。为了将其余的 $ n_k^i - 1 $ 个新生代粒子均匀地分配到父粒子周围，按式(5)将以父粒子为中心点的0~360°分为 $ n_k^i - 1 $ 个区间，每个区间都将进化出1个新生代粒子，在各区间内随机选取1个角度值，存入数组 $ \theta \left[ {n_k^i - 1} \right] $ 中，并打乱顺序。同时调整速度的方向，最终得到的 $ n_k^i $ 个子代粒子分布如式(6)：

$ \theta \left[ i \right]\left. { \in \left( {\left( {i - 1} \right) \cdot \frac{{2{\text π}}}{{n_k^i - 1}},i \cdot \frac{{2{\text π}}}{{n_k^i - 1}}} \right.} \right],i = 1,2,\cdots,n_k^i - 1 $

(5)

$ \begin{gathered} {\boldsymbol{x}}_k^i,{\boldsymbol{x}}_k^i + {\boldsymbol{\alpha}} \left( 1 \right)\frac{{\boldsymbol{\gamma}} }{{{M_k}}},{\boldsymbol{x}}_k^i + {\boldsymbol{\alpha}} \left( 2 \right)\frac{{2 \cdot {\boldsymbol{\gamma}} }}{{{M_k}}},\cdots, \hfill \\ {\boldsymbol{x}}_k^i + {\boldsymbol{\alpha}} \left( {n_k^i - 1} \right)\frac{{\left( {n_k^i - 1} \right) \cdot {\boldsymbol{\gamma}} }}{{{M_k}}} \hfill \\ \end{gathered} $

(6)

式中

$ {\boldsymbol{\gamma}} = \left( {\begin{array}{*{20}{c}}{{\gamma _{\boldsymbol{x}}}}&0&0&0 \\ 0&{{\gamma _{\boldsymbol{x}}}}&0&0 \\ 0&0&{{\gamma _{\boldsymbol{v}}}}&0 \\ 0&0&0&{{\gamma _{\boldsymbol{v}}}} \end{array}} \right) $

$ {\boldsymbol{\alpha}} \left( i \right) ={\left( {\sin \left( {\theta \left( i \right)} \right),\cos\left( {\theta \left( i \right)} \right),\sin \left( {\theta \left( i \right)} \right),\cos\left( {\theta \left( i \right)} \right)} \right)^{\rm{T}}} $

$ {\boldsymbol{\gamma}} $ / $ {M_k} $ 代表变异步长，步长决定了新生代粒子的扩散范围，其中 $ {M_k} = \max \left( {n_k^i,i = 1,2,\cdots,N} \right) $ 。因为目标在位置和速度上的变异幅度相差较大，故本文将扩散因子 $ {\boldsymbol{\gamma}} $ 分为2个：位置扩散因子 $ {\gamma _{\boldsymbol{x}}} $ 和速度扩散因子 $ {\gamma _{\boldsymbol{v}}} $ 。本文的扩散方式使权值越大的父粒子得到子粒子越多，且扩散半径越大。图1(b)显示了一个父粒子按本文的改进变异算子进化生成30个子粒子的结果，相较于传统变异算子将子粒子分布在父粒子两侧的方式，本文的分布方式使子粒子覆盖的区域更广，从而增加了粒子的多样性，更有可能得到接近目标真实运动状态的子粒子，经过不断地迭代，跟踪误差将逐渐减小。

3 仿真分析

以正北为0°，采样周期T为1 s，共采样600 s。观测站是运动的，它的初始坐标为(0，0)，前300 s和后300 s运动方向分别为60°和135°，速度保持9 m/s。目标与观测站的初始距离为8 500 m，初始方位角为45°。目标运动情况如表1所示，速度为14 m/s，在61~120 s和181~240 s两段时间内目标机动。

粒子数N取500个，设目标方位角、速度、航向、相对观测站距离误差符合零均值高斯分布，且 $ {\sigma _\beta } $ =0.1°， $ {\sigma _{\boldsymbol{v}}} $ =0.2 m/s， $ {\sigma _{\rm{course}}} $ =0.1°， $ {\sigma _r} $ =0.1 m，将收敛时间定义为目标距离误差缩小至10%的时间^[15]，并用均方根误差(root mean square error，RMSE)衡量各算法的性能，均方根误差越小说明算法的性能越好。为保证算法性能测量的准确性，本文将每个算法仿真30次。分别采用系统重采样粒子滤波(system resampling particle filter，SRPF)、传统遗传粒子滤波(genetic resampling particle filter，GRPF)和本文改进遗传粒子滤波(improved genetic resampling particle filter，IGRPF)得到的纯方位目标跟踪结果以及各算法在迭代过程中的距离误差百分比如图2所示，位置均方根误差RMSEx和速度均方根误差RMSEv如图3所示。

$ {\rm{RMSE}} = \sqrt {\frac{1}{{\rm{MC}}}\sum\limits_{i = 1}^{\rm{MC}} {{{\left( {{X_i} - {{\hat X}_i}} \right)}^2}} } $

(7)

式中：MC为仿真次数，本文取30； $ {X_i} $ 和 $ {\hat X_i} $ 分别表示第i次仿真的参数真值和估计值。

3.1 性能比较

从图2的各算法跟踪结果中可以看出，SRPF的距离误差在600 s内始终未收敛，GRPF和IGRPF算法的收敛时间分别为468 s和399 s，可见本文的IGRPF算法能有效提升纯方位机动目标跟踪的收敛速度；结合图3和表2可知，本文的IGRPF算法在跟踪机动目标时得出的位置均方跟误差相对SRPF和GRPF分别减少了25.2%和13.0%，速度均方根误差相对SRPF减少了4.1%，与GRPF基本一致，由此验证了本文的改进遗传粒子滤波算法扩大了新生代粒子的分布范围、增加了粒子的多样化变异幅度。高质量的新生代粒子有效防止目标机动时状态估计时陷入孤立点，并使得更多的粒子接近目标真实运动状态，从而更准确地跟踪机动目标，也说明本文的改进遗传重采样方法提升了系统的整体估计精度，改善了估计性能。

3.2 仿真中粒子分析

图4~6是3种算法在目标第一次机动不久后(取t=150 s)的粒子分布情况，能够反映粒子在不断预测修正和更新后的状态。由图4~6可见，SRPF和GRPF由于粒子变化范围较小、变异幅度缺乏多样性的问题造成了较大的跟踪误差，分别是945.9 m和839.7 m；而本文的IGRPF算法在遗传重采样时，权重较小的父粒子被抛弃，权重越大的父粒子经过选择和变异算子处理后生成了更多的子粒子，子粒子加权累加后得到的目标估计位置误差更小，为724.4 m。图6(a)中父粒子的半径与权重的大小成正比，举例的父粒子权重较大，进化后在图6(b)中生成了5个子粒子，扩大了粒子的分布范围。经过不断预测修正和更新之后，IGRPF算法的粒子能更快靠近目标的真实位置。

4 结论

本文针对遗传粒子滤波算法在重采样时粒子变异方向和幅度单一，在处理纯方位机动目标跟踪问题上易造成误差较大和收敛时间长的问题，提出以下改进：

1)调整了新生代粒子的分布方式，扩大了其分布的范围。

2)根据纯方位运动目标跟踪模型将扩散因子分为位置扩散因子和速度扩散因子，使粒子变异的幅度更合理。

经仿真实验验证，本文的改进遗传粒子滤波算法能有效提升纯方位机动目标跟踪的收敛速度并减小跟踪误差，提升系统的估计性能，以便更精准地完成水下目标探测。然而，本算法仅在仿真数据上进行验证，今后将逐步开展试验，以验证算法的普适性和有效性。