CN101980044B - 未知测量噪声分布下的多目标跟踪方法 - Google Patents

Abstract

本发明涉及目标跟踪技术领域，特别是未知测量噪声分布下的多目标跟踪方法，本发明采用风险值衡量粒子的优劣，通过风险评估函数评估粒子后再进行权值更新，整个过程不依赖于测量噪声分布，可以在未知的复杂测量噪声分布下稳定跟踪多个目标，比一般的粒子概率假设密度滤波技术鲁棒性和稳定性更高。

Description

未知测量噪声分布下的多目标跟踪方法

技术领域

本发明涉及目标跟踪技术领域，特别是未知测量噪声分布下的多目标跟踪方法。

背景技术

目标跟踪的目的是利用跟踪滤波算法实时地估计运动目标的位置、速度和加速度等运动状态。和单目标跟踪相比，多目标跟踪由于面临着目标数时变、观测值与目标对应关系的不确定性和观测杂波等情况，变得更加复杂。

现有的多目标跟踪技术主要有两大类，一类是传统的多目标跟踪技术，另一类是基于随机有限集(Random Finite Set)的多目标跟踪技术。传统多目标跟踪技术需要通过数据关联(Data Association)确定观测值与各个目标之间的对应关系，因此也称为基于关联的多目标跟踪技术，这类技术存在着一些不足之处：首先，跟踪的精度很大程度上取决于数据关联的准确程度，同时，数据关联需要耗费大量的时间，影响跟踪的速度。其次，难以处理目标的出现和消失，即目标数目的变化。

基于随机有限集的多目标跟踪技术是在随机有限集基础上通过有限集统计理论实现的。与传统方法不同的是，该技术将多目标状态和观测分别定义成一个随机有限集变量，通过贝叶斯最优递推的方式估计多目标后验概率分布，并得到多目标状态。该技术可以有效克服上述传统多目标跟踪技术的不足，避免数据关联，且很容易处理目标个数的变化。在实现时，通常采用贝叶斯递推传递多目标后验分布的一阶矩——概率假设密度(Probability Hypothesis Density)来近似最优贝叶斯递推，这称为概率假设密度滤波。这种滤波方法使得基于随机有限集的多目标方法能够在计算上实现，例如：Ba-Ngu Vo,W.K.Ma.The Gaussian Mixture Probability Hypothesis Density Filter.IEEE Transaction on signalprocessing,2006,54(11):4091-4104.和Ba-Ngu Vo,S.Singh,A.Doucet.Sequential Monte Carlomethods for Multi-target Filtering with Random Finite Sets.IEEE Transaction on Aerospaceand electronic systems,2005,41(4):1224-1245.这两篇文献公开的技术分别在线性高斯和更广泛意义的非线性非高斯情况下实现了概率假设密度滤波，其中后者的应用领域更为广泛，称为粒子概率假设密度滤波(Particle Probability Hypothesis Density Filtering)。粒子概率假设密度滤波技术利用带权粒子集近似多目标的概率假设密度，通过贝叶斯递推传递和更新每个粒子，达到近似各个时刻概率假设密度的目的。具体来说，可分为四个步骤实现：(1)把上一时刻的粒子集通过重要性函数传递到当前时刻，并计算每个粒子对应的权值。(2).在接收到最新的观测值后，更新每个粒子的权值。(3).从粒子中提取目标的个数和目标状态值。(4).根据粒子的权值对粒子进行重采样操作，得到当前时刻的粒子集。

但是，按以上四个步骤实现的粒子概率假设密度滤波技术进行跟踪滤波时，需要依赖观测噪声的概率分布来计算似然函数才能更新粒子权值，而实际跟踪过程中，测量噪声往往比较复杂，难以准确得到测量噪声的概率分布，如果采用与实际测量噪声不符、不准确的噪声分布更新粒子权值，将使得算法的跟踪性能降低甚至发散。另一方面，滤波更新需要在多维观测空间中计算似然函数，这将耗费大量的时间，降低跟踪速度。

发明内容

本发明的目的是提供一种未知测量噪声分布下的多目标跟踪方法，以便通过消除滤波算法对测量噪声的依赖性克服粒子概率假设密度滤波技术的不足，提高跟踪的鲁棒性、稳定性和速度。

本发明的技术方案是，未知测量噪声分布下的多目标跟踪方法，其方法特征是：假设步骤101已经跟踪得到了k-1时刻的状态粒子集为

则进行下一步目标跟踪需要102的步骤：

对前一时刻(k-1)的粒子集，按公式（1）从建议分布

和p_k(·|Z_k)中随机抽取粒子：

${\tilde{x}}_k^{\left(i\right)}~\left\{\begin{array}{c}{q}_k(\·|x_{k-1}^{\left(i\right)},Z_k),i=1,...L_{k-1}\\ p_k(\·|Z_k),i=L_{k-1}+1,...L_{k-1}+J_k\end{array}\right.---\left(1\right)$

其中，～表示从概率分布中随机采样，

为k-1时刻第i个状态粒子，

为对应的权值，L_k-1为k-1时刻的粒子数，J_k为抽取的新目标状态粒子数，Z_k为观测值，

为抽取得到的预测粒子；

然后通过公式(2)计算所抽取的粒子对应的权值，最终得到103步中k时刻的预测粒子集： ${\{{\tilde{x}}_k^{\left(i\right)},{\tilde{w}}_{k|k-1}^{\left(i\right)}\}}_{i=1}^{L_{k-1}+J_k}:$

${\tilde{w}}_{k|k-1}^{\left(i\right)}=\left\{\begin{array}{c}\frac{{\mathrm{\β}}_{k|k-1}\left({\tilde{x}}_k^{\left(i\right)}\right|x_{k-1}^{\left(i\right)})+P_{s,k}\left({\tilde{x}}_k^{\left(i\right)}\right)f_{k|k-1}\left({\tilde{x}}_k^{\left(i\right)}\right|x_{k-1}^{\left(i\right)})w_{k-1}^{\left(i\right)}}{q_k\left({\tilde{x}}_k^{\left(i\right)}\right|x_{k-1}^{\left(i\right)},Z_k)},i=1,...,L_{k-1}\\ \frac{{\mathrm{\γ}}_k\left({\tilde{x}}_k^{\left(i\right)}\right)}{J_k{P}_k\left({\tilde{x}}_k^{\left(i\right)}\right|Z_k)},i=L_{k-1}+1,...,L_{k-1}+J_k\end{array}\right.---\left(2\right)$

其中，γ_k(·)表示k时刻新生目标概率假设密度，β_k|k-1(·|·)为k时刻衍生目标的概率假设密度，p_s,k(·)为目标的存活概率，f_k|k-1(·|·)为单目标转移概率密度；

随后进行步骤104，利用风险函数

计算每个粒子的风险值，然后通过风险评估函数

评估粒子的风险；其中，

表示第i个粒子的风险值，z为观测值，h(·)为测量方程，‖·‖表示向量2-范数，评估函数的参数λ的取值范围为：0.9<λ<1.1,本发明实施过程中，取λ＝1；

在经过步骤105接收到当前k时刻的最新观测值后，步骤106利用评估结果按照公式(3)更新粒子的权值：

${\tilde{w}}_k^{\left(i\right)}=[1-p_{D,k}\left({\tilde{x}}_k^{\left(i\right)}\right)+\underset{z\&Element;Z_k}{\mathrm{\&Sigma;}}\frac{p_{D,k}\left({\tilde{x}}_k^{\left(i\right)}\right)\mathrm{\&mu;}\left(R_k^{\left(i\right)}\right)}{{\mathrm{\&kappa;}}_k\left(z\right)+\underset{j=1}{\overset{L_{k-1}+J_k}{\mathrm{\&Sigma;}}}{p}_{D,k}\left({\tilde{x}}_k^{\left(j\right)}\right)\mathrm{\&mu;}\left(R_k^{\left(i\right)}\right){\tilde{w}}_{k|k-1}^{\left(j\right)}}]{\tilde{w}}_{k|k-1}^{\left(i\right)}---\left(3\right)$

步骤107从公式(3)中得到更新后的粒子集其中，取κ_k(z)＝rU(z)，为目标的检测率，Z_k为当前时刻接收到的测量值集合，z∈Z_k为测量集合中的一个测量值，κ_k(z)为杂波的概率假设密度，r为杂波个数，服从泊松分布，U(z)为每个杂波的分布，取为均匀分布，

为粒子风险评估后的结果；

步骤108中，将所有的粒子权值累加得到目标个数的估计值：

利用k-means聚类算法对粒子进行聚类，从所得到的聚类中选择峰值最大的前

个类，提取其峰值对应的状态点作为多目标状态估计值；

步骤109步：重采样粒子集

得到新的粒子集调整重采样后的粒子权值得到110步中的当前时刻的最终粒子集合

然后返回101步骤，将最终粒子集合

做为101步骤的(k-1)的粒子集，重复下一时刻的跟踪，其中，L_k为重采样后的粒子个数。

所述建议分布

和p_k(·|Z_k)的分别取单目标转移概率密度f_k|k-1(·|·)和新生目标概率假设密度γ_k(·)，即

和p_k(·|Z_k)＝γ_k(·)。

本发明具有的积极效果是：在实际跟踪环境中，测量噪声的分布复杂且难以估计，使得依赖于测量噪声分布的传统概率假设密度滤波技术适应性和稳定性差。本发明采用风险值衡量粒子的优劣，通过风险评估函数评估粒子后再进行权值更新，整个过程不依赖于测量噪声分布，可以在未知的复杂测量噪声分布下稳定跟踪多个目标，比一般的粒子概率假设密度滤波技术鲁棒性和稳定性更高。同时，粒子更新过程在一维空间中进行，能节省大量的计算时间，提高跟踪的速度。

附图说明

下面结合实施例附图对本发明作进一步说明：

图1是本发明实施例方法步骤过程图；

图2为第一个仿真实验中，在已知测量噪声分布下对四个运动目标进行跟踪得到的X方向运动轨迹，直线表示四个目标的X方向真实运动轨迹，圆圈‘o’表示一般的粒子概率假设密度滤波跟踪技术(图中表示为SMC-PHDF)跟踪目标得到的X方向运动轨迹，叉号‘×’表示本发明的跟踪技术(图中表示为RE-PHDF)跟踪目标得到的X方向运动轨迹。

图3为第一个仿真实验中，在已知测量噪声分布下对四个运动目标进行跟踪得到的Y方向运动轨迹，直线表示四个目标的Y方向真实运动轨迹，圆圈‘o’表示一般的粒子概率假设密度滤波跟踪技术(图中表示为SMC-PHDF)跟踪目标得到的Y方向运动轨迹，叉号‘×’表示本发明的跟踪技术(图中表示为RE-PHDF)跟踪目标得到的Y方向运动轨迹。

图4为第一个仿真实验中，在已知测量噪声分布下对四个运动目标进行跟踪得到的目标个数估计，直线表示真实的目标个数，圆圈‘o’表示一般的粒子概率假设密度滤波跟踪技术(图中表示为SMC-PHDF)估计的目标个数，叉号‘×’表示本发明的跟踪技术(图中表示为RE-PHDF)估计的目标个数。

图5为第一个仿真实验中，在已知测量噪声分布下对四个运动目标进行跟踪得到的跟踪误差，误差采用Wasserstein距离定义，圆圈‘o’表示一般的粒子概率假设密度滤波跟踪技术(图中表示为SMC-PHDF)的跟踪误差，叉号‘×’表示本发明的跟踪技术(图中表示为RE-PHDF)的跟踪误差。

图6为第二个仿真实验中，在未知测量噪声分布下对四个运动目标进行跟踪得到的X方向运动轨迹，标注符号意义与图2相同。

图7为第二个仿真实验中，在未知测量噪声分布下对四个运动目标进行跟踪得到的Y方向运动轨迹，标注符号意义与图3相同。

图8为第二个仿真实验中，在已知测量噪声分布下对四个运动目标进行跟踪得到的目标个数估计，标注符号意义与图4相同。

具体实施方式

如图1所示，假设101步已经跟踪得到了k-1时刻的状态粒子集为

则进行一步目标跟踪需要102的步骤：

对前一时刻(k-1)的粒子集，按公式（1）从建议分布

和p_k(·|Z_k)中随机抽取粒子：

${\tilde{x}}_k^{\left(i\right)}~\left\{\begin{array}{c}{q}_k(\&CenterDot;|x_{k-1}^{\left(i\right)},Z_k),i=1,...L_{k-1}\\ p_k(\&CenterDot;|Z_k),i=L_{k-1}+1,...L_{k-1}+J_k\end{array}\right.---\left(1\right)$

其中，～表示从概率分布中随机采样，

为k-1时刻第i个状态粒子，为对应的权值，L_k-1为k-1时刻的粒子数，J_k为抽取的新目标状态粒子数，Z_k为观测值，

为抽取得到的预测粒子；所述建议分布

和p_k(·|Z_k)的选取可以有多种不同的形式，本发明按照普遍的做法，分别取单目标转移概率密度f_k|k-1(·|·)和新生目标概率假设密度γ_k(·)，即

和p_k(·|Z_k)=γ_k(·)。

然后通过公式(2)计算所抽取的粒子对应的权值，最终得到103步中k时刻的预测粒子集： ${\{{\tilde{x}}_k^{\left(i\right)},{\tilde{w}}_{k|k-1}^{\left(i\right)}\}}_{i=1}^{L_{k-1}+J_k}:$

${\tilde{w}}_{k|k-1}^{\left(i\right)}=\left\{\begin{array}{c}\frac{{\mathrm{\&beta;}}_{k|k-1}\left({\tilde{x}}_k^{\left(i\right)}\right|x_{k-1}^{\left(i\right)})+P_{s,k}\left({\tilde{x}}_k^{\left(i\right)}\right)f_{k|k-1}\left({\tilde{x}}_k^{\left(i\right)}\right|x_{k-1}^{\left(i\right)})w_{k-1}^{\left(i\right)}}{q_k\left({\tilde{x}}_k^{\left(i\right)}\right|x_{k-1}^{\left(i\right)},Z_k)},i=1,...,L_{k-1}\\ \frac{{\mathrm{\&gamma;}}_k\left({\tilde{x}}_k^{\left(i\right)}\right)}{J_k{P}_k\left({\tilde{x}}_k^{\left(i\right)}\right|Z_k)},i=L_{k-1}+1,...,L_{k-1}+J_k\end{array}\right.---\left(2\right)$

其中，γ_k(·)表示k时刻新生目标概率假设密度，β_k|k-1(·|·)为k时刻衍生目标的概率假设密度，p_s,k(·)为目标的存活概率，f_k|k-1(·|·)为单目标转移概率密度；本发明中，取建议分布

和p_k(·|Z_k)分别为单目标转移概率密度f_k|k-1(·|·)和新生目标概率假设密度γ_k(·)，即

和p_k(·|Z_k)＝γ_k(·)。

随后进行步骤104，利用风险函数

计算每个粒子的风险值，然后通过风险评估函数

评估粒子的风险；

其中，

表示第i个粒子的风险值，z为观测值，h(·)为测量方程，‖·‖表示向量2-范数，评估函数的参数λ的取值范围为：0.9<λ<1.1,本发明实施过程中，取λ=1；

在经过步骤105接收到当前k时刻的最新观测值后，步骤106利用评估结果按照公式(3)更新粒子的权值：

${\tilde{w}}_k^{\left(i\right)}=[1-p_{D,k}\left({\tilde{x}}_k^{\left(i\right)}\right)+\underset{z\&Element;Z_k}{\mathrm{\&Sigma;}}\frac{p_{D,k}\left({\tilde{x}}_k^{\left(i\right)}\right)\mathrm{\&mu;}\left(R_k^{\left(i\right)}\right)}{{\mathrm{\&kappa;}}_k\left(z\right)+\underset{j=1}{\overset{L_{k-1}+J_k}{\mathrm{\&Sigma;}}}{p}_{D,k}\left({\tilde{x}}_k^{\left(j\right)}\right)\mathrm{\&mu;}\left(R_k^{\left(i\right)}\right){\tilde{w}}_{k|k-1}^{\left(j\right)}}]{\tilde{w}}_{k|k-1}^{\left(i\right)}---\left(3\right)$

步骤107从公式(3)中得到更新后的粒子集

其中，取κ_k(z)＝rU(z)，为目标的检测率，Z_k为当前时刻接收到的测量值集合，z∈Z_k为测量集合中的一个测量值，κ_k(z)为杂波的概率假设密度，r为杂波个数，服从泊松分布，U(z)为每个杂波的分布，取为均匀分布，

为粒子风险评估后的结果。

步骤108中，将所有的粒子权值累加得到目标个数的估计值：

利用k-means聚类算法对粒子进行聚类，从所得到的聚类中选择峰值最大的前

个类，提取其峰值对应的状态点作为多目标状态估计值；

步骤109步：重采样粒子集得到新的粒子集调整重采样后的粒子权值得到110步中的当前时刻的最终粒子集合

然后返回101步骤，将最终粒子集合

做为101步骤的(k-1)的粒子集，重复下一时刻的跟踪，其中，L_k为重采样后的粒子个数。

下面给出仿真实例说明。根据跟踪过程中对测量噪声先验知识的掌握程度，给出两个仿真实验来验证本发明的可行性和相比一般技术的优点。实验均在Intel Pentium 41.80GHz处理器,512M内存的PC机上，使用MATLAB 7.1平台实现，整个仿真的时间步长为40步，每步持续一个时间单位(一个时间单位为2秒)。在XY二维平面区域[-300300]×[-300300]m²中，先后出现四个做匀加速(constant acceleration,CA)的运动目标，每个目标的运动方程(CA模型)和测量方程为：

$\left\{\begin{array}{c}{X}_k={\mathrm{FX}}_{k-1}+W_{k-1}\\ Z_k={\mathrm{HX}}_k+V_k\end{array}\right.$

其中：

$F=\left(\begin{array}{cccccc}1& T& \frac{T^2}{2}& 0& 0& 0\\ 0& 1& T& 0& 0& 0\\ 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 1& T& \frac{T^2}{2}\\ 0& 0& 0& 0& 1& T\\ 0& 0& 0& 0& 0& 1\end{array}\right),$ $H=\left(\begin{array}{cccccc}1& 0& 0& 0& 0& 0\\ 0& 0& 0& 1& 0& 0\end{array}\right)$ 采样周期T=1(时间单位)过程噪声为零均值的高斯白噪声：W_k-1～N(0,Q),Q＝diag(0.5,0.1,0.01,0.5,0.1,0.01)，其中：diag(·)表示对角矩阵。目标状态向量为：X_K＝[x_k,xv_k,xa_k,y_k,yv_k,ya_k]^T，测量向量为：Z_k＝[x_k,y_k]^T，其中：x_k、xv_k、xa_k分别是目标在X方向的位置、速度和加速度，y_k、yv_k、ya_k分别是目标在Y方向的位置、速度和加速度。跟踪滤波时每个目标的粒子数取为300，目标运动区域内的杂波数服从均值r=4的泊松分布，且每个杂波在目标运动区域内服从均匀分布，即杂波的概率假设密度为k＝r/400²。目标的检测概率Pd=0.98，目标的存活概率Ps=0.99，不考虑目标的衍生情况。新生目标的概率假设密度为

其中：

Q_r＝diag(1,4,0.1,1,4,0.1)，diag(·)表示对角矩阵。为方便叙述，在仿真实验中,SMC-PHDF表示一般的粒子概率假设滤波技术，RE-PHDF表示本发明的滤波技术。

仿真实验一

首先，考虑一种比较理想的噪声情况：测量噪声是高斯分布，且能准确获得其先验知识。假设实际测量噪声的概率分布如下：

V_k～N(0,R)， $R=\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right)$

且进行滤波估计时完全掌握了测量噪声的这一概率分布，即掌握的测量噪声概率分布先验知识为：N(0,R_prj)，其中： $R_{\mathrm{prj}}=R=\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right).$ 实验中，SMC-PHDF采用这个获取的测量噪声分布计算似然函数，并更新粒子权值，实验仿真中取λ＝1。在该噪声环境下，对场景中的四个目标进行跟踪，得到的跟踪轨迹如图2、3所示，可以看出两种跟踪滤波技术都可以较准确的跟踪场景中的所有目标，且从图4的目标个数的估计来看，两种技术对于实际目标个数的估计也比较准确，且两者的估计效果非常接近。为了更准确的比较两种算法的跟踪性能，采用Wasserstein距离来度量跟踪技术的跟踪误差。Wasserstein距离定义为：任意两个非空集合

和X，它们的Wasserstein距离：

其中：

表示对所有满足条件的传输矩阵C取最小，一个传输矩阵C满足条件：C_i,j≥0，

|X|表示集合X的元素个数。当

和X任意一个为空时，Wasserstein距离没有定义。图5是两种跟踪技术的多目标跟踪误差(Wasserstein距离)，从误差曲线可以更直观的看出，RE-PHDF的跟踪精度非常接近SMC-PHDF，甚至在某些点上的精度高于SMC-PHDF。表1给出了经过20次蒙特卡罗仿真得到的算法平均运行时间比较。

表1.两种算法的平均运行时间比较

算法	SMC-PHDF	RE-PHDF
			时间(s)	78.2	38.7

从表中的数据看出：在精度相似的情况下，RE-PHDF算法的运行时间要明显比SMC-PHDF少，这在对实时性要求较高的实际应用中具有重要意义。

仿真实验二

在实际跟踪环境中，通常无法准确掌握测量噪声分布的先验知识，实验二在未知测量噪声分布情况下，比较本发明的跟踪技术和一般技术的性能。假设实际测量噪声的概率分布为混合高斯形式：V_k～0.2N(0,R₁)+0.6N(0,R₂)+0.2N(0,R₃)，其中 $R_1=\left(\begin{array}{cc}1& 0\\ 0& 2\end{array}\right),$ $R_2=\left(\begin{array}{cc}0.1& 0\\ 0& 0.2\end{array}\right),$ $R_3=\left(\begin{array}{cc}4& 0\\ 0& 6\end{array}\right).$ 但是滤波估计时没有准确的获取到测量噪声的这一先验概率分布，假设所获取的测量噪声的先验概率分布为：N(0,R_prj)，其中 $R_{\mathrm{prj}}=R_2=\left(\begin{array}{cc}0.1& 0\\ 0& 0.2\end{array}\right).$ 实验中，SMC-PHDF采用这个获取的测量噪声概率分布计算似然函数，并更新粒子权值，实验仿真中取λ＝1。采用SMC-PHDF和本发明的跟踪技术进行跟踪实验，得到的跟踪结果如图6、7、8所示。其中图6、7是两种跟踪技术对场景中目标进行跟踪得到的轨迹，图8是两种跟踪技术估计的目标个数。从图8可以看出，SMC-PHDF技术除了前几个时间步准确估计到了实际目标个数以外，其它时间步上估计出的都是错误的目标个数。而本发明的跟踪技术较准确的估计到了实际目标个数，在目标个数估计的性能上要优于SMC-PHDF技术。同时，从图6、7的跟踪轨迹看，SMC-PHDF技术虽然可以在目标出现以后的起始几个时间步上跟踪目标轨迹，但是由于测量噪声先验知识的不足，造成误差的积累，导致算法很快就丢失目标，甚至所有目标——即跟踪算法发散(如果跟踪算法从某时刻开始直到仿真结束持续丢失所有目标，则认为跟踪算法发散)。而本发明的技术能以较高的精度持续跟踪场景中出现的所有目标，直到目标消失，并且对于目标的多次交叉同样可以取得较好的跟踪效果。通过20次蒙特卡罗仿真对两种技术进行稳定性比较，采用多次试验中的算法发散次数来衡量两种跟踪技术的稳定性，比较结果如表2所示：

表2.两种算法的发散次数比较

算法	SMC-PHDF	RE-PHDF
			发散次数	16	0

从中可以看出：在掌握不充分的测量噪声先验知识情形下，SMC-PHDF极易发散，性能难以保证，而本发明的技术具有更高的稳定性和鲁棒性。

Claims (1)

1.未知测量噪声分布下的多目标跟踪方法，其特征在于：

步骤101：假设已经跟踪得到了k-1时刻的状态粒子集为

则进行下一步目标跟踪需要步骤102；

步骤102：对k-1时刻的状态粒子集，按公式（1）从建议分布

和p_k(·|Z_k)中随机抽取粒子：

${\tilde{x}}_k^{\left(i\right)}~\left\{\begin{array}{c}{q}_k(\&CenterDot;|x_{k-1}^{\left(i\right)},Z_k),i=1,...L_{k-1}\\ p_k(\&CenterDot;|Z_k),i=L_{k-1}+1,...L_{k-1}+J_k\end{array}\right.---\left(1\right)$

其中，～表示从概率分布中随机采样，

为k-1时刻第i个状态粒子，

为对应的权值，L_k-1为k-1时刻的粒子数，J_k为抽取的新目标状态粒子数，Z_k为k时刻测量值集合，

为抽取得到的预测粒子，所述建议分布

和p_k(·|Z_k)分别取单目标转移概率密度f_k|k-1(·|·)和新生目标概率假设密度γ_k(·)，即 $q_k(\&CenterDot;|x_{k-1}^{\left(i\right)},Z_k)=f_{k|k-1}(\&CenterDot;|x_{k-1}^{\left(i\right)})$ 和p_k(·|Z_k)=γ_k(·)；

然后通过公式(2)计算所抽取的粒子对应的权值；

步骤103：得到k时刻的预测粒子集：

${\tilde{w}}_{k|k-1}^{\left(i\right)}=\left\{\begin{array}{c}\frac{{\mathrm{\&beta;}}_{k|k-1}\left({\tilde{x}}_k^{\left(i\right)}\right|x_{k-1}^{\left(i\right)})+P_{s,k}\left({\tilde{x}}_k^{\left(i\right)}\right)f_{k|k-1}\left({\tilde{x}}_k^{\left(i\right)}\right|x_{k-1}^{\left(i\right)})w_{k-1}^{\left(i\right)}}{q_k\left({\tilde{x}}_k^{\left(i\right)}\right|x_{k-1}^{\left(i\right)},Z_k)},i=1,...,L_{k-1}\\ \frac{{\mathrm{\&gamma;}}_k\left({\tilde{x}}_k^{\left(i\right)}\right)}{J_k{P}_k\left({\tilde{x}}_k^{\left(i\right)}\right|Z_k)},i=L_{k-1}+1,...,L_{k-1}+J_k\end{array}\right.---\left(2\right)$

其中，γ_k(·)表示k时刻新生目标概率假设密度，β_k|k-1(·|·)为k时刻衍生目标的概率假设密度，p_s,k(·)为目标的存活概率，f_k|k-1(·|·)为单目标转移概率密度；

步骤104：利用风险函数

计算每个粒子的风险值，然后通过风险评估函数

评估粒子的风险；其中，

表示第i个粒子的风险值，z为测量值，h(·)为测量方程，‖·‖表示向量2-范数，评估函数的参数λ的取值范围为：0.9<λ<1.1；

步骤105：接收到当前k时刻的最新观测值；

步骤106：利用评估结果按照公式(3)更新粒子的权值：

${\tilde{w}}_k^{\left(i\right)}=[1-p_{D,k}\left({\tilde{x}}_k^{\left(i\right)}\right)+\underset{z\&Element;Z_k}{\mathrm{\&Sigma;}}\frac{p_{D,k}\left({\tilde{x}}_k^{\left(i\right)}\right)\mathrm{\&mu;}\left(R_k^{\left(i\right)}\right)}{{\mathrm{\&kappa;}}_k\left(z\right)+\underset{j=1}{\overset{L_{k-1}+J_k}{\mathrm{\&Sigma;}}}{p}_{D,k}\left({\tilde{x}}_k^{\left(j\right)}\right)\mathrm{\&mu;}\left(R_k^{\left(i\right)}\right){\tilde{w}}_{k|k-1}^{\left(j\right)}}]{\tilde{w}}_{k|k-1}^{\left(i\right)}---\left(3\right)$

步骤107：从公式(3)中得到更新后的粒子集

其中，取κ_k(z)=rU(z)，为目标的检测率，Z_k为当前k时刻接收到的测量值集合，z∈Z_k为k时刻测量值集合中的一个测量值，κ_k(z)为杂波的概率假设密度，r为杂波个数，服从泊松分布，U(z)为每个杂波的分布，取为均匀分布，

为粒子风险评估后的结果；

步骤108：将所有的粒子权值累加得到目标个数的估计值：

利用k-means聚类算法对粒子进行聚类，从所得到的聚类中选择峰值最大的前

个类，提取其峰值对应的状态点作为多目标状态估计值；

步骤109：重采样粒子集

得到新的粒子集

步骤110：调整重采样后的粒子权值,得到当前时刻的最终粒子集合

然后返回步骤101，将最终粒子集合

作为步骤101的k-1时刻的状态粒子集，重复下一时刻的跟踪，其中，L_k为重采样后的粒子个数。

Images