随着水下无人航行器(underwater unmanned vehicle，UUV)的逐渐兴起，无人探测系统的研究引起了学者极大的重视，其智能程度往往直接影响航行器的性能以及任务执行能力。无人探测系统与传统声呐系统一个极大的差别在于没有了人为判别和甄选目标的过程，只通过信号处理系统自动完成。这一过程在水下复杂的环境中进行，需要处理可能出现的目标中断、虚警野值等情况，涉及了阈值的判定、多目标的分辨与关联等多个难点问题，因此也成为提高无人探测系统性能的一个瓶颈。

目前成熟的目标方位跟踪算法包括最近邻关联方法；联合概率数据关联(joint probability data association, JPDA)算法及其快速算法^[1]；递推卡尔曼滤波及由其改进而来的扩展卡尔曼滤波(extended kalman filter, EKF)^[2]、无迹卡尔曼滤波(unscented kalman filter, UKF)^[3]等。其中，最近邻数据关联以及JPDA算法运算量小，易于系统实现，但是都精度不高，在目标交汇、目标数较多的情况下有较大可能无法给出正确的跟踪结果。扩展卡尔曼滤波算法部分解决了这一问题，但是在目标非线性运动情况下可能无法收敛，造成跟踪错误。进一步改进的无迹卡尔曼滤波算法降低了目标跟踪过程无法收敛的可能性，其主要思想是在状态估计值附近进行确定性采样，称为Sigma采样，使样本的均值和协方差接近于目标状态的概率分布，再将这些点代入状态更新和测量更新中，以获得非线性变换后的统计特性，所选取的Sigma采样策略直接影响了目标跟踪的性能。在水下目标跟踪这一时变空变过程中，经过选取的采样策略也会带来一定的误差。

粒子滤波是一种处理非高斯、非线性、非平稳信号参数估计和滤波的有效方法^[4]，在目标跟踪定位领域已经得到广泛应用^[5-6]。这一算法可以利用目前时间之前所有观测信息。在波束形成提供了目标方位和强度信息后，粒子滤波采用一系列具有不同权重的粒子进行建模，对每个目标的每种假设进行验证，从中选取贝叶斯意义上的最大后验概率结果输出。这一方法具有较强的稳定性，已经在实际项目中取得了成功的应用，其中较为著名的有用于语音目标跟踪以及源分离的ManyEars系统^[7]等。哈尔滨工程大学马珊等^[8]、声学所的许枫等^[9]都将粒子滤波成功应用在了水下目标跟踪。已有文献[10]指出在水下环境应用中EKF、UKF都有失效的情况，粒子滤波算法性能具有明显优势。

本文利用波束形成结果包含能量信息的特性，采用一种基于能量值的粒子滤波方法，并结合波束内插以及自动阈值检测的方法，利用方位信息实现了目标自动跟踪。通过数值仿真以及湖上实验与快速JPDA预置跟踪算法比较。

1 粒子滤波基本模型 1.1 粒子滤波跟踪原理

很多实际问题，包括目标跟踪问题，都可以用具有一个初始分布的隐马尔科夫过程模型来描述。

首先假定目标运动是一个马尔科夫过程。令k时刻目标的状态向量为{x^k; k∈N, x^k∈Rⁿ}，其中n为状态向量维数，N为非负整数集合。初始情况下目标的概率分布为p(x₀)，k时刻目标的观测向量为{y^k; k∈N, y^k∈R^m}，其中m为观测向量维数。建立目标状态与观测状态的关系，可以得到隐马尔科夫模型：

$ \begin{array}{l} p\left( {{\mathit{\boldsymbol{x}}^k}\left| {{\mathit{\boldsymbol{x}}^{k - 1}}} \right.} \right)\\ p\left( {{\mathit{\boldsymbol{y}}^k}\left| {{\mathit{\boldsymbol{x}}^{k - 1}}} \right.} \right) \end{array} $

(1)

把时刻k以及之前所有目标j的状态向量和观测向量分别看作一个集合，定义为${\mathit{\boldsymbol{x}}^{0:k}} \buildrel \Delta \over = \left\{ {{\mathit{\boldsymbol{x}}^0},{\mathit{\boldsymbol{x}}^1}, \cdots ,{\mathit{\boldsymbol{x}}^k}} \right\}$和${\mathit{\boldsymbol{y}}^{1:k}} \buildrel \Delta \over = \left\{ {{\mathit{\boldsymbol{y}}^1},{\mathit{\boldsymbol{y}}^2}, \cdots ,{\mathit{\boldsymbol{y}}^k}} \right\}$。在这一模型中，跟踪问题的求解即通过状态向量求目标j的参数f(x)的期望值：

$ E\left( {f_j^k} \right) = \int {f_j^k\left( {{\mathit{\boldsymbol{x}}^{0:k}}} \right)p\left( {{\mathit{\boldsymbol{x}}^{0:k}}\left| {{\mathit{\boldsymbol{y}}^{1:k}}} \right.} \right){\rm{d}}{\mathit{\boldsymbol{x}}^{0:k}}} $

(2)

这一积分无法通过解析的方法求解，粒子滤波采用重要性采样的蒙特卡洛积分方法给出了一个近似解。蒙特卡罗积分用一组在确定概率分布下采样的，独立同分布的样本来估计积分的解值，也就是将问题中的概率p(x^0:k|y^1:k)近似为^[11]

$ p\left( {{\mathit{\boldsymbol{x}}^{0:k}}\left| {{\mathit{\boldsymbol{y}}^{1:k}}} \right.} \right) \approx \frac{1}{N}\sum\limits_{i = 1}^N {\delta \left( {\mathit{\boldsymbol{x}} - {\mathit{\boldsymbol{x}}_i}} \right)} $

(3)

式中：N为采样的样本数目，x_i为对分布p(x^0:k|y^1:k)采样得到的样本值。每个样本可以称为一个“粒子”，所求参数的期望值即这些粒子给出的参数均值。

1.2 序贯重要性采样

给出式(2) 数值解法的原理之后，概率p(x^0:k|y^1:k)的值仍然没有显式的表示，直接运算得到积分值仍然不可行。在跟踪问题中，所关心的往往不是整个状态序列的估计，而是当前目标状态的估计。假设p(x^0:k-1|y^1:k-1)已知，利用贝叶斯定理以及目标状态的马尔科夫性假设，文献[4]给出了递推表达式：

$ \begin{array}{l} p\left( {{\mathit{\boldsymbol{x}}^{0:k}}\left| {{\mathit{\boldsymbol{y}}^{1:k}}} \right.} \right) = \\ \frac{{p\left( {{\mathit{\boldsymbol{y}}^k}\left| {{\mathit{\boldsymbol{x}}^{0:k}},{\mathit{\boldsymbol{y}}^{1:k - 1}}} \right.} \right)p\left( {{\mathit{\boldsymbol{x}}^{0:k}}\left| {{\mathit{\boldsymbol{y}}^{1:k - 1}}} \right.} \right)}}{{p\left( {{\mathit{\boldsymbol{y}}^k}\left| {{\mathit{\boldsymbol{y}}^{1:k - 1}}} \right.} \right)}} = \\ \frac{{p\left( {{\mathit{\boldsymbol{y}}^k}\left| {{\mathit{\boldsymbol{x}}^k}} \right.} \right)p\left( {{\mathit{\boldsymbol{x}}^k}\left| {{\mathit{\boldsymbol{x}}^{k - 1}}} \right.} \right)p\left( {{\mathit{\boldsymbol{x}}^{0:k - 1}}\left| {{\mathit{\boldsymbol{y}}^{1:k - 1}}} \right.} \right)}}{{p\left( {{\mathit{\boldsymbol{y}}^k}\left| {{\mathit{\boldsymbol{y}}^{1:k - 1}}} \right.} \right)}} \propto \\ p\left( {{\mathit{\boldsymbol{y}}^k}\left| {{\mathit{\boldsymbol{x}}^k}} \right.} \right)p\left( {{\mathit{\boldsymbol{x}}^k}\left| {{\mathit{\boldsymbol{x}}^{k - 1}}} \right.} \right)p\left( {{\mathit{\boldsymbol{x}}^{0:k - 1}}\left| {{\mathit{\boldsymbol{y}}^{1:k - 1}}} \right.} \right) \end{array} $

(4)

假设q(x)是一个已知的易于采样的函数，称为建议分布，代入式(2)、(4) 可以得到：

$ \begin{array}{l} E\left( {f_j^k} \right) = \int {f_j^k\left( {{\mathit{\boldsymbol{x}}^{0:k}}} \right)\frac{{p\left( {{\mathit{\boldsymbol{x}}^{0:k}}\left| {{\mathit{\boldsymbol{y}}^{1:k}}} \right.} \right)}}{{q\left( x \right)}}q\left( x \right){\rm{d}}{\mathit{\boldsymbol{x}}^{0:k}}} \approx \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\frac{1}{N}\sum\limits_{i = 1}^N {f_j^k\left( {\mathit{\boldsymbol{x}}_i^{0:k}} \right)\frac{{p\left( {\mathit{\boldsymbol{x}}_i^{0:k}\left| {\mathit{\boldsymbol{y}}_i^{1:k}} \right.} \right)}}{{q\left( {{\mathit{\boldsymbol{x}}_i}} \right)}} \buildrel \Delta \over = } \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\frac{1}{N}\sum\limits_{i = 1}^N {f_j^k\left( {\mathit{\boldsymbol{x}}_i^{0:k}} \right)w_{ji}^k} \end{array} $

(5)

式中：w_ji^k称作目标j的粒子i在k时刻的归一化“权值”，代表了粒子在粒子群模拟目标j概率密度函数中的“重要性”。代入式(4)，得到权值的递推公式：

$ w_{ji}^k \propto \frac{{p\left( {{\mathit{\boldsymbol{y}}^k}\left| {\mathit{\boldsymbol{x}}_i^k} \right.} \right)p\left( {\mathit{\boldsymbol{x}}_i^k\left| {\mathit{\boldsymbol{x}}_i^{k - 1}} \right.} \right)}}{{q\left( {{\mathit{\boldsymbol{x}}_i}\left| {\mathit{\boldsymbol{x}}_i^{k - 1}{\mathit{\boldsymbol{y}}^k}} \right.} \right)}}w_{ji}^{k - 1} \buildrel \Delta \over = W_{ji}^k $

(6)

计算权值w_ji^k可以近似为计算与其成正比的非归一化权值W_ji^k，再归一化：

$ w_{ji}^k = \frac{{W_{ji}^k}}{{\sum\limits_{i = 1}^N {W_{ji}^k} }} $

(7)

1.3 系统重采样

在选取合适的建议分布后，粒子滤波即可根据递推公式不断迭代，逐渐接近目标状态的概率分布。但是由于对目标状态实际分布基本没有先验信息，所选取的建议分布往往与实际情况相差较远。这样，不同粒子权值之间差异越来越大，产生粒子退化现象。这种情况下大量粒子权重很小，对结果几乎不产生什么影响。一方面产生了大量的无效运算，另一方面剩下的粒子如果数量过少，则对目标状态概率分布的描述过于模糊甚至无效化。粒子退化现象将导致算法的不稳定性。

重采样方法解决了粒子退化的问题。其原理是将经过经验积累的离散分布作为更加接近目标状态实际概率分布的建议分布，重新采样得到新的式(3)，重新开始递推过程。重采样后的粒子具有相同的初始权值。也就是将粒子的重要性权值变成了新分布中的粒子数目。

$ {{\mathit{\boldsymbol{\hat p}}}_N}\left( \mathit{\boldsymbol{x}} \right) = \sum\limits_{i = 1}^N {\frac{{{n_i}}}{N}\delta \left( {\mathit{\boldsymbol{x}} - {\mathit{\boldsymbol{x}}_i}} \right)} $

(8)

式中：${{\hat p}_N}\left( \mathit{\boldsymbol{x}} \right)$表示重采样后的粒子分布期望，n_i表示重采样为x_i的粒子个数，δ为狄拉克函数。

常用的重采样算法包括多项式采样、分层采样、系统重采样和残差采样等。文献[12-13]对常用重采样算法性能进行了比较，以上几种重采样方法都是无偏的，且性能差距并不大。考虑到实时性要求带来的运算量问题，本文采用了系统重采样方法。首先对粒子权重进行排序，排序后的粒子为{x₀, x₁, …, x_N}。令：

$ {U^{\left( i \right)}} = \frac{{i - 1}}{N}I + U $

(9)

式中U是均匀分布U((0, I/N])上一次采样的结果。

则满足式(9) 的对x_i采样为

$ {U^{\left( i \right)}} \in \left( {\sum\nolimits_{j = 1}^{i - 1} {{w_j}} ,\sum\nolimits_{j = 1}^i {{w_j}} } \right] $

(10)

重采样操作虽然解决了粒子退化问题，但是也将导致随机的粒子丢失，降低粒子的多样性。频繁的重采样将加重这一问题，甚至可能使粒子滤波跟踪算法收敛到固定值。因此，粒子分布有效性N_eff^[14]定义为

$ {N_{{\rm{eff}}}} \approx {\left( {\sum\limits_{i = 1}^N {w_{ji}^{\left( 2 \right)}} } \right)^{ - 1}} < 0.7N $

(11)

为与时标区别，用“(2)”代表平方运算标志。仅当式(11) 成立时，采取重采样操作，减少了不必要的重采样过程。

2 跟踪算法实现

除了上述的粒子权重更新以及重采样过程，要实现水下目标方位能量值与方位信息结合的跟踪，还需要在粒子滤波基本算法中加入其他部分。

2.1 目标方位预测

假设模型中在t时刻有M个目标，每个目标用N个粒子进行建模，每个粒子分别有位置参数x_ji^t，权重参数w_ji^t，其中目标编号j∈{0, 1, …, M-1}，粒子编号i∈{0, 1, …, N-1}。定义目标状态参数s_ji^t=$\left[ {\mathit{\boldsymbol{x}}_{ji}^t,\mathit{\boldsymbol{\dot x}}_{ji}^t} \right]$。由于目标仅仅具有角度信息，状态参数需要经过距离归一化，即认为任一时刻目标都位于在其方位角与单位圆的交点上，满足条件|x_ji^t|=1。

为将不同时刻之间目标状态建立联系，提高目标状态的连续性，利用目标状态的马尔科夫性假设，每个时刻开始之前，通过每个目标上一时刻状态s_ji^t-1预测目前时刻状态s_ji^t。使用常用的激励-衰减模型^[15]，这一模型用2个变量很好地模拟了不同动态信号源的运动，表示为：

$ \begin{array}{l} \mathit{\boldsymbol{\dot x}}_{ji}^t = a\mathit{\boldsymbol{\dot x}}_{ji}^{t - 1} + b{F_{\rm{x}}}\\ \mathit{\boldsymbol{x}}_{ji}^t = \mathit{\boldsymbol{x}}_{ji}^{t - 1} + \Delta T\mathit{\boldsymbol{\dot x}}_{ji}^t \end{array} $

(12)

式中：a=e^-αΔT控制激励部分，$b = \beta \sqrt {1 - {\alpha ^2}} $控制衰减部分，F_x是正态分布的随机变量且方差为1，ΔT是两次信息更新之间的时间间隔。参数α、β的选取分为3种不同情况：静止目标(α=2, β=0.04)；匀速目标(α=0.05, β=0.2)；加速目标(α=0.5, β=0.2)。本文方法对这3种情况分别赋予了一定出现概率，以保证粒子的多样性。经过实验，3种情况概率分别设置为0.05、0.9、0.05时，可以得到较好的跟踪效果。预测完成后对结果进行距离归一化处理。

2.2 目标方位检测

粒子滤波方位跟踪的原始数据是不同方位波束能量值。UUV平台基阵孔径较小，带来的阵增益也较小，原始信号的信噪比较低。在这种条件下，绝大部分高精度测向算法都很有可能在信道情况较差时失效。因此选用鲁棒性更强的波束形成算法计算波束能量值，其方位分辨率较低，结合阵列和算法性能，在作者工程实践中方位扫描间隔为1°。如果直接将波束能量结果作为粒子滤波的目标参数输入，则可能出现目标跳跃的情况。比如某缓慢行驶的目标在相对UUV方位角60°的方位，可能出现连续40次探测结果均为60°，第41次目标方位角却突然跳变到61°的情况。

这一现象造成了粒子滤波目标跟踪多方面的不利因素。首先，方位结果跳变容易造成目标方位预测中模型跳变，粒子权重不能有效分配。其次，方位结果跳变会造成序贯重要性采样和重采样中目标概率的不准确，所模拟的目标状态概率分布也不符合实际情况。最后，即使跟踪能够有效持续的进行，目标形成的运动轨迹也不平滑，对平台自主决策产生不利影响。

为了解决方位跳变问题，本文在波束能量结果的基础上，通过二次拟合、插值的方法估计波束能量最高值的精确方位。文献[16]指出指向性函数在主瓣附近可以用抛物线函数较好的进行估计。

令波束能量结果的最高值角度为k，使用k及与其相邻2个结果点方位的能量估计值E_k-1、E_k、E_k+1，进行二次曲线拟合。设拟合的二次曲线为y=ax²+bx+c，拟合结果能量最高的角度为k+Δθ，可以得到联立的等式：

$ \left\{ \begin{array}{l} {E_{k - 1}} = a{\left( {k - 1} \right)^2} + b\left( {k - 1} \right) + c\\ {E_k} = a{k^2} + bk + c\\ {E_{k + 1}} = a{\left( {k + 1} \right)^2} + b\left( {k + 1} \right) + c\\ 0 = a{\left( {k + \Delta \theta } \right)^2} + b\left( {k + \Delta \theta } \right) + c \end{array} \right. $

(13)

解得：

$ \Delta \theta = \frac{1}{2}\frac{{{E_{k + 1}} - {E_{k - 1}}}}{{2{E_k} - {E_{k + 1}} - {E_{k - 1}}}} $

(14)

拟合得到的潜在目标波束能量：

$ E = a{\left( {k + \Delta \theta } \right)^2} + b\left( {k + \Delta \theta } \right) + c $

(15)

2.3 基于能量值的潜在信号源真实概率

由式(15) 给出的目标能量信息可以推断其为真实信号源的概率。JPDA、卡尔曼滤波跟踪等方法的目标判定普遍基于硬门限判定，所有强度高于门限的目标赋予存在概率1，所有强度低于门限的目标赋予存在概率0，这一过程损失了波束能量值包含的目标存在概率信息。本文所采用的能量值和方位信息结合的粒子滤波方法则将能量值定量转换为潜在信号源为真实信号源的概率，充分利用了波束能量信息。

假设波束能量检测的输出结果由Q个潜在信号源组成，其中潜在信号源q的波达方向为y_q。令P_q为潜在信号源q是真实信号源(即非虚警)的概率，P_q的值可以看作对波束能量检测输出结果的信心，波束能量越大，潜在信号源真实的概率也就越大。定义P_q为^[17]

$ {P_q} = \left\{ \begin{array}{l} {\upsilon ^2}/2,\;\;\;\;\;\;\;\;\;\;\;\;q = 0,\upsilon \le 1\\ 1 - {\upsilon ^{ - 2}}/2,\;\;\;\;\;\;\;q = 0,\upsilon > 1\\ 0.3,\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;q = 1\\ 0.16,\;\;\;\;\;\;\;\;\;\;\;\;\;\;q = 2\\ 0.03,\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;q = 3 \end{array} \right. $

(16)

式(16) 是一个半经验公式，在出现多个潜在信号源的情况下，为了排除潜在目标之间的相互影响，设定能量较小的潜在信号源存在概率为一个确定值。υ为目标相对强度，其常用的定义是υ=E₀/E_T，E₀为潜在信号源q的波束能量，E_T为常数阈值。但是在被动声纳接收到的信号中，信号源形式以及信道传播过程非常复杂，背景也具有时变空变特性，信号较大的不确定性使常数阈值方法只在部分相对理想的场景下有效。因此，引入了波束均值E和均方差${\hat E}$两个统计量，重新定义目标相对强度：

$ \upsilon = \frac{{\left( {{E_0} - \bar E} \right)}}{{\hat E}} $

(17)

使用相对强度代替固定阈值，在背景噪声整体抬高或降低的时，潜在信号源存在概率随之变化的幅度大大减小，增强了检测跟踪算法的稳定性。

2.4 结合能量值信息的多目标方位关联

对任意一个潜在信号源q，有3种可能的假设：虚警(H₀)；对应一个跟踪中的存在目标(H₁)；对应一个尚未跟踪的新目标(H₂)。在H₁假设下，还需要确定潜在信号源与跟踪中目标的对应关系。

根据目标跟踪的实际物理情况，一个潜在信号源最多与一个跟踪中目标相关联，一个跟踪中目标最多与一个潜在信号源相关联。设跟踪中目标的数目为M，定义一个关联f:{0, 1, …, Q-1}→{-2, -1, 0, 1, …, M-1}为将潜在信号源集合映射到信号源状态集合的函数，信号源状态集合中“-2”表示虚警，“-1”表示尚未跟踪的新目标，其余值表示相应编号的跟踪中目标。

t时刻的潜在信号源q对应虚警和潜在新目标的概率分别为所有可能关联下相应概率之和，分别表示为

$ P_q^t\left( {{H_0}} \right) = \sum\limits_f {P\left( {f\left| {{\mathit{\boldsymbol{y}}^t}} \right.} \right)\delta \left( {f\left( q \right) + 2} \right)} $

(18)

$ P_q^t\left( {{H_2}} \right) = \sum\limits_f {P\left( {f\left| {{\mathit{\boldsymbol{y}}^t}} \right.} \right)\delta \left( {f\left( q \right) + 1} \right)} $

(19)

式中：y^t为t时刻的观测量。对所有可能的关联计算其条件概率P(f|y^t)，可以得到潜在信号源q关联到跟踪中目标j的概率P_qj^t：

$ P_{qj}^t = \sum\limits_f {P\left( {f\left| {{\mathit{\boldsymbol{y}}^t}} \right.} \right)\delta \left( {f\left( q \right) + j} \right)} $

(20)

式中：δ(f(q)-j)为狄拉克函数。为书写方便，去掉式中的t。

对任一关联的条件概率P(f|y)应用贝叶斯定理可得

$ P\left( {f\left| \mathit{\boldsymbol{y}} \right.} \right) = \frac{{p\left( {\mathit{\boldsymbol{y}}\left| f \right.} \right)P\left( f \right)}}{{p\left( \mathit{\boldsymbol{y}} \right)}} $

(21)

由于只有一个正确的关联，有

$ \sum\limits_f {P\left( {f\left| \mathit{\boldsymbol{y}} \right.} \right) = 1} $

(22)

所以通过归一化可以省去的p(y)计算，主要需要关注p(y|f)和P(f)的计算。假设在给定关联下每个观测结果为条件独立的，则

$ P\left( {\mathit{\boldsymbol{y}}\left| f \right.} \right) = \prod\limits_q {p\left( {{\mathit{\boldsymbol{y}}_q}\left| {f\left( q \right)} \right.} \right)} $

(23)

假设虚警和新目标的概率分布为均匀分布，潜在信号源q与跟踪中目标j对应的概率则用波束能量曲线和相应粒子群分布近似估计：

$ p\left( {{\mathit{\boldsymbol{y}}_q}\left| {f\left( q \right)} \right.} \right) = \left\{ \begin{array}{l} 1/4{\rm{ \mathsf{ π} ,}}\;\;\;f\left( q \right) = - 1, - 2\\ \sum\limits_i {{w_{f\left( q \right)}},{}_ip\left( {{\mathit{\boldsymbol{y}}_q}\left| {{\mathit{\boldsymbol{x}}_{ji}}} \right.} \right)} ,\;\;\;f\left( q \right) \ge 0 \end{array} \right. $

(24)

其中p(y_q|x_ji)为以该粒子参数为均值的正态分布。

假设每个潜在信号源的关联之间都是条件独立的，则f是正确关联的概率为

$ P\left( f \right) = \prod\limits_q {P\left( {f\left( q \right)} \right)} $

(25)

其中任意一个潜在信号源q关联正确的概率为

$ P\left( {f\left( q \right)} \right) = \left\{ \begin{array}{l} \left( {1 - {P_q}} \right){P_{{\rm{false}}}},\;\;\;\;\;\;f\left( q \right) = - 2\\ {P_q}{P_{{\rm{new}}}},\;\;\;\;\;\;\;\;\;\;\;\;\;\;f\left( q \right) = - 1\\ {P_q}P\left( {E_j^t\left| {{\mathit{\boldsymbol{y}}^{t - 1}}} \right.} \right),\;\;\;\;f\left( q \right) \ge 0 \end{array} \right. $

(26)

式中P_false和P_new分别是关联到虚警和新目标的先验概率，分别为0.05。P(E_j^t|y^t-1)是目标j真实存在的概率，表达为一个递推式：

$ P\left( {{E_j}\left| {{\mathit{\boldsymbol{y}}^{t - 1}}} \right.} \right) = P_j^{t - 1} + \left( {1 - P_j^{t - 1}} \right)\frac{{{P_o}P\left( {{E_j}\left| {{\mathit{\boldsymbol{y}}^{t - 2}}} \right.} \right)}}{{1 - \left( {1 - {P_o}} \right)P\left( {{E_j}\left| {{\mathit{\boldsymbol{y}}^{t - 2}}} \right.} \right)}} $

(27)

式中：P_o是目标存在但未观测到的先验概率，P_j^t=∑_qP_qj^t是目标j被观测到的概率，即与每一个潜在信号源关联的概率之和。

2.5 目标增加和删除

实际探测过程中，目标随时可能出现和消失。因此，何时加入新的目标以及何时删除已经消失的跟踪中目标是跟踪过程中的重要步骤。

本文算法中将新目标出现分为2个阶段，分别是判断潜在新目标出现并加入跟踪中目标序列；以及确认跟踪中目标存在并加入跟踪结果输出。判断潜在新目标的出现主要依靠潜在目标出现概率P_q^t(H₂)，当P_q^t(H₂)大于设定门限(本文中采用0.3) 时，认为出现了潜在新目标，并为此潜在新目标建立粒子群，加入跟踪中目标序列。确认跟踪中目标存在则主要依靠目标存在概率P(E_j|y^t)。确认目标存在前，其存在概率P(E_j|y^t)不断迭代更新，当目标存在概率大于设定值(本文中采用0.98) 时，认为目标确实存在，将潜在目标加入目标序列并输出，同时设定P(E_j|y^t)为1并不再继续更新。

目标消失情况仅做一个阶段判断。当目标一段时间内没有观察到，也就是P_j^t＜G_dis条件在持续的一段时间中始终满足，G_dis为依据经验得到的距离函数，则认为目标消失，从目标序列中删除相应目标及其粒子。

3 数据验证结果

为说明粒子滤波跟踪方法的性能，选取快速JPDA多目标预置跟踪方法与其进行跟踪结果的比较。

3.1 仿真结果

首先建立仿真信号模型。设定一个目标的方位轨迹之后，认为其波束能量在每一帧中不断起伏变化。把仿真数据中一个目标所有波束能量的值看作一个集合，假定这一集合取值符合正态分布。这样一来，每次的波束能量都可以用正态分布的抽样得到。最后加入高斯白噪声，形成仿真信号。在不同信噪比之下，已知高斯白噪声的能量，可以求出相应信号强度。以求出的信号强度作为能量正态分布的均值，抽样得到仿真信号能量，纯净信号与白噪声相加形成不同信噪比之下的仿真信号。

在高斯白噪声背景下，分别在50°、100°方位加入了相对强度分别为3、2 dB的2个目标源。假设目标强度符合正态分布。得到以灰度图显示的目标方位历程如图 1、2所示。

由图 1、2可知，粒子滤波给出了非常准确的跟踪结果，而预置跟踪则无法得到目标信息。由以上结果可以看出，粒子滤波跟踪比预置跟踪具有更加稳定的性能，针对弱目标探测跟踪性能更强。

为了得到更加明确的仿真结论，定量说明2种算法性能的差异，使用蒙特卡洛方法^[18]对2种算法的性能进行仿真。仿真了多种信噪比情况下的信号，并考虑了较低信噪比下的情况，得到2种算法的性能对比如图 3所示。对低信噪比下目标的探测，粒子滤波给出了明显较好的结果，这一结论在之后的湖试实验中也进一步得到了印证。

3.2 单目标实验结果

选取某次UUV湖试结果。UUV平台两侧各配备了一个32阵元的、有后挡板的均匀阵元间距舷侧阵。方位角度角度定义为以UUV艉部方向为0°角，顺时针递增。

首先验证在只存在单个目标的情况下，粒子滤波方位跟踪的有效性。所选取的时段为晚间，周围较大水域都没有路过行船，只有实验拖船拖曳合作声源航行。合作声源选用了2~5kHz宽带模拟白噪声信号，与UUV都在水深约15 m左右位置，近似认为在同一水平面。UUV本艇以及合作声源都记录了GPS信息，以此解算得到目标相对本艇方位基本呈线性变化，得到以灰度图显示的方位历程图如图 4所示。

信噪比较高时，单目标情况下2种方法都能够成功跟踪到目标。由于湖试条件限制，没有找到低信噪比条件下的单目标数据进一步验证单目标情况下算法性能。

3.3 多目标实验结果

为验证存在多个目标的情况下，粒子滤波方位跟踪的有效性，选择某次湖试一段时间进行测试。

由于湖试过程中没有多个合作声源发射，多目标探测中目标均为非合作声源。UUV在实验基地船旁边固定吊放，测量周围行船的航迹。根据实验记录还原的目标航迹如图 5所示。在起始时刻，目标1在UUV右后方，正在向右前方行驶，航迹如箭头1所示；目标2在UUV右前方，正在向右后方快速行驶，航迹如箭头2所示。在开始记录之后120 s左右，UUV右侧出现目标3，从UUV的右前方向右后方行驶，航迹如箭头3所示。

两种方法的跟踪结果分别如图 6所示，图中曲线表示目标方位角轨迹。

以灰度图显示的方位历程图如图 7所示，0°~179°及180°~359°分别为左侧、右侧所配备的舷侧阵计算得到的波束结果。两阵均有后障板，由于目标位于右阵，加上环境噪声、本艇噪声等因素，造成背景强度略有不同。这一现象并不影响目标跟踪过程。

2种方法的跟踪结果分别如图 8所示，图中不同标号线表示不同目标的方位角轨迹。

图 8中每条线表示一个跟踪结果目标的方位角轨迹，按照目标出现时间分别用一个标号表示。图 8(a)中粒子滤波跟踪得到的1、2、5号目标方位历程曲线分别对应图 5中的目标1、2、3，方位跟踪结果完全正确；3、4号目标所对应的是能量泄漏到左舷的假目标。图 8(b)中预置跟踪得到的1号目标方位历程曲线对应图 5中目标1；3号目标前半段与2号目标后半段方位历程曲线对应图 5中目标2的部分轨迹，在目标轨迹交叉时产生了错误判断；4号目标方位历程曲线对应图 5中目标3的轨迹。

对照图 7中的方位历程及图 8中的跟踪结果，在目标情况复杂情况下，预置跟踪的逻辑判断显得相对混乱，尤其在目标交叉情况下难以判断各个目标分别的走向。能量值和方位信息结合的粒子滤波准确得到了两个目标分别的走向。同时，在200 s之后的低强度目标检测跟踪中，粒子滤波也得到了正确的结果。

4 结论

1) 在数值仿真和蒙特卡洛仿真中，在信噪比较低的情况下，本文跟踪方法相较JPDA方法具有明显更高的检测概率，其稳健性以及对弱目标跟踪能力强。

2) 在湖试实验中，这一方法在单目标情况下以及多目标交叉的复杂情形下均给出了各目标的正确航迹，表现出了很强的实用性和稳定性，多目标交叉的跟踪效果相较JPDA方法得到很大提高。

本文算法对弱目标较强的跟踪能力虽然得到了仿真结果的证实，但是在实验中没有相应的参考结果。因此，还需进一步设计实验验证本文方法的弱目标检测性能。