CN106019253A - 基于箱粒子cphd的多扩展目标跟踪方法 - Google Patents

Abstract

本发明提出了一种基于箱粒子CPHD的多扩展目标跟踪方法，用于解决现有箱粒子PHD多扩展目标跟踪方法中存在的目标数目估计不精确的技术问题，实现步骤包括：对多扩展目标量测集合进行量测划分；初始化箱粒子集和势分布；对初始箱粒子集进行预测，对初始势分布进行预测；对预测箱粒子集进行量测更新，对预测势分布进行量测更新；搜索对预测箱粒子强度贡献最大的量测箱；若搜索到对预测箱粒子强度贡献最大的量测箱，对更新箱粒子集进行收缩；估计目标数目；采用随机子划分方法对收缩箱粒子集进行重采样；提取目标状态。本发明具有目标数目估计精确的特点，可用于车辆编队和人群等***中的跟踪。

Description

基于箱粒子CPHD的多扩展目标跟踪方法

技术领域

本发明属于信息处理技术领域，涉及一种多扩展目标跟踪方法，具体涉及一种基于箱粒子CPHD的多扩展目标跟踪方法，可用于车辆编队和人群等***中的跟踪。

背景技术

随着现代先进传感器如相控阵雷达和逆合成孔径雷达等探测设备分辨率的不断提高，每个运动目标的回波信号可能会同时占据多个距离分辨单元，即每个目标可能产生多个量测，称这样的目标为扩展目标。扩展目标跟踪不仅可以提供目标的运动状态，包括位置、速度、加速度，还可以提供目标的扩展形态，包括目标大小、形状等多方面信息，因此引起了国内外学者的广泛关注。在现实跟踪场景中，存在杂波和虚警，往往无法预知扩展目标的个数，对数目未知且变化的多扩展目标跟踪已成为目标跟踪领域具有挑战的热门课题。

多扩展目标跟踪方法可分为基于数据关联的多扩展目标跟踪方法和基于随机集的多扩展目标跟踪方法两大类。在基于数据关联的多扩展目标跟踪方法中，由于多扩展目标会产生大量量测，每个量测与每个扩展目标进行关联假设，会带来巨大的数据计算量，导致跟踪实时性较差。基于随机集的多扩展目标跟踪方法，将目标状态集合和量测集合均建模成随机集，跟踪过程中避免了复杂的数据关联，具有良好的实时性。Mahler结合扩展目标空间分布模型和随机有限集理论提出的扩展目标PHD滤波算法有效地避免了复杂的数据关联，对多扩展目标进行有效跟踪；随后，Orguner等又提出了多扩展目标的CPHD滤波算法。箱粒子(Box-Particle)滤波是由Abdallah等提出的一种基于区间分析的蒙特卡罗滤波算法，能够有效处理量测的三种不确定性。在实际应用中，箱粒子滤波只需要用几十个的箱粒子代替即可达到粒子滤波中几百几千个粒子的滤波效果，大大降低了运算复杂度。Schikora等给出PHD滤波器的箱粒子实现形式；Petrov等提出箱粒子多扩展目标滤波算法及其应用。

2015年宋骊平等在《控制与决策》第10期的文章《基于箱粒子的多扩展目标PHD滤波》中，提出了箱粒子PHD多扩展目标跟踪技术，文中采用区间分析来处理扩展目标问题，不需要区分真实量测和混入目标扩展范围的杂波量测，只需通过确定目标扩展范围来进行计算，具有良好的抗杂波性，同时能够大大降低计算复杂度。然而，由于此跟踪技术采用对k-1时刻的初始箱粒子强度进行预测和对k时刻的预测箱粒子强度进行量测更新，得到k时刻的更新箱粒子强度，再通过k时刻的更新箱粒子强度对目标数目进行估计，在此过程中，只传播多扩展目标强度，忽略了目标数目的高阶信息，因其目标数目估计的精确度不高，导致多扩展目标跟踪精确度降低。

发明内容

本发明的目的在于克服上述现有技术存在的缺陷，提出了一种基于箱粒子CPHD的多扩展目标跟踪方法，通过在预测和量测更新过程中，同时传播扩展目标强度和目标数目势分布信息，用于解决现有箱粒子PHD多扩展目标跟踪方法中存在的目标数目估计不精确的技术问题。

为实现上述目的，本发明采用的技术方案，包括如下步骤：

(1)将k时刻的多扩展目标量测集合Z_k划分成多个划分单元并利用包含函数[f]分别对多个划分单元进行包含，得到k时刻的量测箱粒子集其中k表示滤波时刻，m_k表示k时刻的划分单元个数，W_j表示第j个划分单元内量测分布的空间范围，表示k时刻的量测箱粒子状态。

(2)假设k-1时刻的存活箱粒子集为初始势分布为p_k-1(n)，将k-1时刻的新生箱粒子集并入存活箱粒子集得到k-1时刻的初始箱粒子集 其中和分别表示k-1时刻的存活箱粒子状态、新生箱粒子状态和初始箱粒子状态，和分别表示k-1时刻的存活箱粒子强度、新生箱粒子强度和初始箱粒子强度，N_k-1,s、N_k-1,b和N_k-1分别表示k-1时刻的存活箱粒子数目、新生箱粒子数目和初始箱粒子数目，n表示目标个数。

(3)对得到的k-1时刻的初始箱粒子集进行预测，得到k时刻的预测箱粒子集同时对假设的k-1时刻的初始势分布p_k-1(n)进行势分布预测，得到k时刻的预测势分布p_k|k-1(n)，其中，和分别表示k时刻的预测箱粒子状态和预测箱粒子强度。

(4)对k时刻的预测箱粒子集和预测势分布p_k|k-1(n)进行量测更新，实现步骤如下：

(4a)根据得到的k时刻的量测箱粒子集和k时刻的预测箱粒子集求得似然函数

(4b)利用求得的似然函数对得到的k时刻的预测箱粒子集进行量测更新，得到k时刻的更新箱粒子集

(4c)利用求得的似然函数对得到的k时刻的预测势分布p_k|k-1(n)进行量测更新，得到k时刻的势分布p_k(n)。

(5)对得到的k时刻的量测箱粒子集进行搜索，并判断是否搜索到对预测箱粒子强度贡献最大的量测箱[z_k]，若是，执行步骤(6)，否则，将得到的k时刻的更新箱粒子集输出。

(6)利用搜索到的量测箱[z_k]，对得到的k时刻的更新箱粒子集进行收缩，得到k时刻的收缩箱粒子集其中表示收缩箱粒子状态。

(7)对得到的k时刻的势分布p_k(n)进行估计，得到目标数目估计

(8)利用随机子划分方法，对得到的k时刻的收缩箱粒子集进行重采样，得到k时刻的重采样箱粒子集为其中和分别表示k时刻的重采样箱粒子状态和重采样箱粒子强度，N表示重采样次数。

(9)对得到的k时刻的重采样箱粒子状态集进行点化，得到k时刻的箱中心粒子状态集对该箱中心粒子状态集进行k-means聚类，得到目标的状态估计其中表示k时刻的箱中心粒子状态。

(10)判断k时刻是否为跟踪最后时刻，若是，则跟踪结束；否则保存k时刻的重采样箱粒子集并作为下个时刻的存活箱粒子集，令k＝k+1，执行步骤(2)。

本发明与现有技术相比，具有以下优点：

本发明由于在对目标数目进行估计时，通过对假设的k-1时刻的初始势分布进行势分布预测和对得到的k时刻的预测势分布进行量测更新，得到k时刻的势分布，再通过k时刻的势分布对目标数目进行估计，保留了目标高阶势信息，与现有箱粒子PHD多扩展目标跟踪方法中采用的对k-1时刻的初始箱粒子强度进行预测和对k时刻的预测箱粒子强度进行量测更新，得到k时刻的更新箱粒子强度，再通过k时刻的更新箱粒子强度对目标数目进行估计的方法对比，有效地提高了目标数目估计精度和稳定性。

附图说明

图1是本发明的流程图；

图2是本发明与现有的箱粒子PHD多扩展目标跟踪方法单次蒙特卡罗仿真跟踪效果对比图；

图3是本发明与现有的箱粒子PHD多扩展目标跟踪方法100次蒙特卡罗仿真的目标数目估计对比图和OSPA距离仿真对比图；

图4是本发明与现有的箱粒子PHD多扩展目标跟踪方法，分别在检测概率为0.95和0.92条件下进行100次蒙特卡罗仿真的目标数目估计对比图。

具体实施方式

以下结合附图和实施例，对本发明作进一步详细说明。

参照图1，本发明包括如下步骤：

(1)对多扩展目标量测集合进行量测划分：

本发明中量测划分采用最大距离划分，将量测集合划分成多个划分单元，由于箱粒子滤波不区分目标扩展范围内目标量测和杂波量测，因此可将得到的划分单元分别用包含函数[f]包含成具有一定大小的矩形箱，矩形箱的大小由划分单元内量测分布的空间范围决定：

${\{\[{\mathrm{zp}}_k^{\left(j\right)}\]\}}_{j=1}^{m_k}=\{\[f\]W_j\}$

其中，m_k为划分单元的个数，W_j表示第j个划分单元内量测分布的空间范围，表示k时刻的量测箱粒子状态。

当目标距离较近，且划分单元内包含的量测数目大于门限值χ时采取子划分。本文中的子划分首先求出量测的聚类中心，再根据划分单元内量测分布进行子划分，最后将划分结果扩展为箱，箱的大小保持之前的划分大小。

(2)初始化箱粒子集和势分布：

假设k-1时刻的存活箱粒子集为初始势分布为p_k-1(n)，在实际场景中，目标和杂波均产生量测，由于杂波是随机的，杂波产生的量测一般会在下一时刻消失，而新生目标则不同，其量测往往能够在下一时刻的传播中存留下来，因此，上一时刻存在量测的地方很有可能会出现新生目标。根据这一假设，在上一时刻的量测箱粒子周围采样产生k-1时刻的新生箱粒子集其状态服从高斯分布权值设置为则k-1时刻的初始箱粒子集为

其中和分别表示k-1时刻的存活箱粒子状态、新生箱粒子状态和初始箱粒子状态，和分别表示k-1时刻的存活箱粒子强度、新生箱粒子强度和初始箱粒子强度，N_k-1,s、N_k-1,b和N_k-1分别表示k-1时刻的存活箱粒子数、新生箱粒子数和初始箱粒子数，n表示目标数目，P_b表示箱粒子新生概率，Σ表示新生箱粒子状态分布协方差。

(3)对初始箱粒子集进行预测和对初始势分布进行预测：

首先，对初始箱粒子集进行预测，箱粒子***预测模型为：

$\[x_{k|k-1}^{\left(i\right)}\]=\[f_{k|k-1}\](\[x_{k-1}^{\left(i\right)}\])$

初始箱粒子强度预测方程为：

$w_{k|k-1}^{\left(i\right)}=P_S\·w_{k-1}^{\left(i\right)}$

其中，f_k|k-1表示状态转移函数，P_S表示箱粒子存活概率。

其次，对初始势分布进行预测，马尔科夫转移矩阵为：

$M(n,n^{\′})=\underset{i=0}{\overset{min\{n,n^{\′}\}}{\Σ}}{p}_{birth}(n-i)\left(\begin{array}{c}{n}^{\′}\\ i\end{array}\right)d^{n^{\′}-i}{(1-d)}^i$

初始势分布预测方程为：

$p_{k|k-1}\left(n\right)=\underset{n^{\′}=0}{\overset{\mathrm{\∞}}{\Σ}}{p}_{k-1}\left(n^{\′}\right)M(n,n^{\′})$

其中，p_birth(n)表示从k-1时刻到k时刻有n个目标新生的概率，d表示目标消失概率，且d＝1-P_S。

(4)对预测箱粒子集进行量测更新和对预测势分布进行量测更新：

首先，根据量测箱粒子集和预测箱粒子集，由以下公式求得似然函数：

$L_k(\[{\mathrm{zp}}_k^{\left(j\right)}\]|\[x_{k|k-1}^{\left(i\right)}\])=\frac{|\[h_{cp}\](\[x_{k|k-1}^{\left(i\right)}\],\[{\mathrm{zp}}_k^{\left(j\right)}\])|}{|\[x_{k|k-1}^{\left(i\right)}\]|}$

$\[h_{cp}\](\[x_{k|k-1}^{\left(i\right)}\],\[{\mathrm{zp}}_k^{\left(j\right)}\])=\[x_{k|k-1}^{\left(i\right)}\]\∩\[{\mathrm{zp}}_k^{\left(j\right)}\]$

其中，h_cp为约束传播算法，在这里用的约束传播采取预测箱粒子状态与量测箱粒子状态重叠交叉的部分作为约束结果输出。

其次，利用求得的似然函数对得到预测箱粒子集进行量测更新，预测箱粒子状态可由下式更新得到：

$\[x_k^{\left(i\right)}\]=\[h_{cp}\](\[x_{k|k-1}^{\left(i\right)}\],\[{\mathrm{zp}}_k^{\left(j\right)}\])$

预测箱粒子强度量测更新可由下式得到：

其中，D和分别表示目标被检测到和未被检测到情况的简记符号，等式右边括号中前半部分为传感器没有检测到目标时的权值更新，后半部分为检测到目标时的权值更新。P_D表示目标检测概率。

最后，利用求得的似然函数对得到预测势分布进行量测更新，更新方程如下式实现：

$p_k\left(n\right)=\frac{L(\[Z_k\]|n)}{L(\[Z_k\])}{p}_{k|k-1}\left(n\right)$

其中，

$L(\[Z_k\]|D)=\underset{j=1}{\overset{m_k}{\Σ}}l(\[{\mathrm{zp}}_k^{\left(s\right)}\]|\[x\])L(\[Z_k\]|a_k^{\left(j\right)}=s)$

$L(\[Z_k\]|a_k^{\left(j\right)}=s)=\frac{1}{n_{k|k-1}}\[\underset{j=1}{\overset{m_k}{\Σ}}{\mathrm{\β}}_k^{\left(j\right)}{\mathrm{\α}}_k^{\left(j\right)}{\mathrm{\σ}}_{j-1}\left(\right\{L_k^{\left(1\right)},...,L_k^{\left(m_k\right)}\}\backslash \{L_k^{\left(s\right)}\left\}\right)\]$

$L(\[Z_k\])=\underset{j=1}{\overset{m_k}{\Σ}}{\mathrm{\α}}^{\left(j\right)}{\mathrm{\β}}^{\left(j\right)}{\mathrm{\σ}}_j\left(\right\{L_k^{\left(1\right)},...,L_k^{\left(m_k\right)}\left\}\right)$

$L(\[Z_k\]|n)=\underset{j=0}{\overset{\min (m,n)}{\Σ}}{\mathrm{\β}}_k^{\left(j\right)}\frac{n!{(1-P_D)}^{n-j}{\mathrm{\σ}}_j\left(\right\{L_k^{\left(1\right)},...,L_k^{\left(m_k\right)}\left\}\right)}{(n-j)!}$

$L_k^{\left(s\right)}=\frac{1}{n_{k|k-1}}\∫P_D(\[x\])w_{k|k-1}(\[x\])l(\[{\mathrm{zp}}_k^{\left(s\right)}\]|\[x\])d\[x\]$

${\mathrm{\α}}^{\left(j\right)}\≡\underset{n=j}{\overset{\mathrm{\∞}}{\Σ}}\frac{n!}{(n-j)!}{p}_{k|k-1}\left(n\right){(1-P_D)}^{n-j}$

${\mathrm{\β}}_k^{\left(j\right)}\≡p_c(m-j)\frac{(m_k-j)!}{m_k!}{\mathrm{\λ}}^{-j}$

p_c(m)表示出现m个虚警的概率，σ_j(·)为初等对称函数，可由下式计算得出：

${\mathrm{\&sigma;}}_j\left(\right\{y_1,...,y_m\left\}\right)=\underset{1\&le;i_1<...<i_j\&le;m}{\&Sigma;}{y}_{i_1}{y}_{i_2}...y_{i_j}$

其中σ₀≡1。

(5)搜索对预测强度贡献最大的量测箱[z_k]：

在扩展目标中，箱粒子是一个包含划分预测范围或者量测划分单元的一定大小的矩形箱，本身就包含很多冗余信息，为了剔除扩展范围内的冗余信息，算法在量测更新之后加入收缩步骤，主要根据相应的量测箱粒子来收缩预测得到的箱粒子。首先需要对得到的量测箱粒子集进行搜索，并判断是否搜索到对预测强度贡献最大的量测箱[z_k]，若是，执行步骤(6)，否则，将得到的k时刻的更新箱粒子集输出，搜索由以下公式进行：

$\&lsqb;z_k\&rsqb;=\underset{w_{j,i}}{\arg \max }\{\&lsqb;{\mathrm{zp}}_k^{\left(j\right)}\&rsqb;,w_{j,i}>0\}$

$w_{j,i}=P_D\frac{L(\&lsqb;Z_k\&rsqb;|D)}{L(\&lsqb;Z_k\&rsqb;)}$

(6)对更新箱粒子集进行收缩：

利用搜索到的对预测强度贡献最大的量测箱[z_k]对k时刻的更新箱粒子集进行收缩，收缩由以下公式进行：

$\&lsqb;{\tilde{x}}_k^{\left(i\right)}\&rsqb;=\&lsqb;h_{cp}\&rsqb;(\&lsqb;x_k^{\left(i\right)}\&rsqb;,\&lsqb;z_k\&rsqb;)$

(7)对量测更新后的势分布进行估计得到目标数目：

对得到的k时刻的势分布p_k(n)进行估计，采用最大后验估计方法，得到目标数目估计为：

${\hat{N}}_k=\underset{n}{\arg \max }{p}_k\left(n\right)$

(8)采用随机子划分方法对收缩箱粒子进行重采样：

区别于传统粒子滤波直接根据粒子权值大小进行复制的重要性重采样策略，本发明算法采用随机子划分重采样方法。首先根据箱粒子权值进行重采样，得到第m个箱粒子所需要采样的次数N^(m),m＝1,…,n，总的采样箱粒子数为然后在位置和速度四个维度中选择数值最大的一维进行子划分，将第m个箱粒子均匀划分为N^(m)个大小合适且在此维度上有均匀差异的箱粒子。这样得到的箱粒子既能满足下一时刻滤波的需要，增加箱粒子多样性，又可以去除箱粒子多余部分。采用这一重采样策略必须注意的是，在每一个时刻都需要根据上一时刻的量测划分结果补入一定量的新生箱粒子，否则算法会迅速发散。由重采样得到的新箱粒子集为：

${\left\{(w_k^{res,\left(i\right)}=1/N,\&lsqb;x_k^{res,\left(i\right)}\&rsqb;)\right\}}_{i=1}^N$

其中，和分别表示k时刻的重采样箱粒子状态和重采样箱粒子强度，N表示重采样次数；

(9)提取目标状态：

首先，将重采样得到的箱粒子集进行点化，求得每个箱中心粒子状态为：

$c_k^{res,\left(i\right)}=mid(\&lsqb;x_k^{res,\left(i\right)}\&rsqb;)$

其中mid([·])为取箱粒子[·]的中心粒子状态。

其次，对得到的点化后的箱中心粒子状态进行K-means聚类，将聚类得出的个中心作为目标的状态估计结果。其中，表示k时刻的箱中心粒子状态。

(10)判断k时刻是否为最后时刻，若是，则跟踪结束；否则保存k时刻的重采样箱粒子集并作为下个时刻的存活箱粒子集，令k＝k+1，执行步骤(2)。

以下结合仿真实验，对本发明的技术效果作进一步描述：

1.仿真条件：

本发明的仿真是在主频3.1GHZ的Intel(R)Core(TM)i3-2100CPU、内存6.00GB的硬件环境和MATLAB R2010a的软件环境下进行的。

模拟一个大小为A＝[-800,1000]×[－1000,1000](m²)且分布有随机噪声的两维观测区域，区域中有4个直线运动的扩展目标，目标数目随时间变化，有目标新生也有目标消失。目标状态方程与量测方程分别为：

[x_k]＝F[x_k－1]+Gω_k－1,

z_k＝Hx_k+v_k,

$\&lsqb;x_k\&rsqb;={\&lsqb;\&lsqb;x_k\&rsqb;,\&lsqb;v_k^x\&rsqb;,\&lsqb;y_k\&rsqb;,\&lsqb;v_k^y\&rsqb;\&rsqb;}^T,$

z_k＝[x_k,y_k]^T,

[zp_k]＝[[xp_k],[yp_k]]^T.

目标状态转移矩阵和观测矩阵分别为：

$F=\left[\begin{array}{cccc}1& T& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& T\\ 0& 0& 0& 1\end{array}\right],G=\left[\begin{array}{cc}{T}^2/2& 0\\ T& 0\\ 0& T^2/2\\ 0& T\end{array}\right],H=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 0& 1& 0\end{array}\right]$

其中，ω_k－1，v_k分别为过程噪声和量测噪声。为目标状态区间，[[x_k],[y_k]]^T为目标位置区间，是目标速度区间。扩展目标的量测是以真实位置为中心，服从高斯分布，而每个目标产生的量测数目则是参数为10的泊松分布，量测扩展范围为12m。在这里量测均为区间量测，是量测划分单元经包含函数形成的量测箱。状态噪声标准差σ_x＝σ_y＝0.1m，量测噪声标准差为对于箱粒子CPHD滤波器，采用99％区间度量和来建模区间噪声。目标存活概率P_S＝0.99，新生概率P_b＝0.01，平均杂波个数r＝3，采样周期T＝1。箱粒子大小为[25,23]，速度区间为[1,1]。最优子模型分配(OSPA)距离参数设置为c＝70和p＝2。在本发明和箱粒子PHD多扩展目标跟踪方法中，每一时刻每个存活的箱粒子数为35个，根据上一时刻的量测划分结果每个划分单元补入的新生箱粒子个数为15。

4个目标相继出现在监测区域，然后相继消失，目标1出现在1时刻，终止于31时刻，初始状态为(-500m,20m/s,-80m,-25m/s)；目标2出现在7时刻，终止于60时刻，初始状态为(-600m,18m/s,80m,10m/s)；目标3出现在25时刻，终止于52时刻，初始状态为(-550m,19m/s,-550m,10m/s)；目标4出现在16时刻，终止于40时刻，初始状态为(200m,20m/s,-40m,-10m/s)。

2.仿真内容：

(1)在检测概率P_D＝0.98条件下，对本发明和箱粒子PHD多扩展目标跟踪方法分别进行跟踪效果、目标数目估计和OSPA距离对比仿真实验，其结果分别如图2，图3(a)和图3(b)所示；

(2)在检测概率P_D＝0.95和P_D＝0.92条件下，对本发明和箱粒子PHD多扩展目标跟踪方法进行目标数目估计对比仿真，其效果如图4(a)和图4(b)所示。

3.仿真结果与分析：

参照图2，图中横坐标表示X方向位置，纵坐标为Y方向位置。线条表示目标在两维空间运动轨迹，符号“o”代表本发明跟踪效果，符号“*”代表箱粒子PHD多扩展目标跟踪方法跟踪效果。从图中可以看出，符号“o”在大部分运动轨迹比符号“*”更吻合线条，即本发明在单次蒙特卡洛仿真的对比中效果更好。

参照图3(a)，横坐标为跟踪时间，纵坐标为目标数目。实线表示真实目标数目，带实黑点的线表示本发明目标数目估计，带空心黑点的线表示箱粒子PHD多扩展目标跟踪方法目标数目估计。图中可以看到，在17s～25s，真实目标数目为3，本发明数目估计大部分时刻在3附近，而箱粒子PHD多扩展目标跟踪方法数目估计大部分时刻在2.7附近；在26s～31s，真实目标数目为4，本发明数目估计大部分时刻在3.9附近，而箱粒子PHD多扩展目标跟踪方法数目估计大部分时刻在3.7附近；通过对比明显看出，本发明目标数目估计更接近真实目标数目，有效提高估计精度。

参照图3(b)，横坐标为跟踪时间，纵坐标为OSPA距离。带实黑点的线表示本发明OSPA距离，带空心黑点的线表示箱粒子PHD多扩展目标跟踪方法OSPA距离。图中可以看到，在17s～40s中，本发明OSPA距离数值大多数时刻处在15到80之间，峰值出现在16s、31s和40s，最大峰值出现在40s值为80左右；箱粒子PHD多扩展目标跟踪方法中OSPA距离数值大多处在40到160之间，峰值出现在29s值为160。大部分时刻本发明OSPA距离要小于箱粒子PHD多扩展目标跟踪方法OSPA距离，即对比可得本发明具有更精确的目标数目估计。

综合图3(a)和图3(b)两幅图本发明和箱粒子PHD多扩展目标跟踪技术的目标数目估计和OSPA距离的对比仿真，可以看出本发明有效地提高了目标数目估计精度。

参照图4，横坐标为跟踪时间，纵坐标为目标数目。实线表示真实目标数目，带实黑点的线表示本发明目标数目估计，带空心黑点的线表示箱粒子PHD多扩展目标跟踪方法目标数目估计。图4(a)中可以看到，在17s～25s，真实目标数目为3，本发明数目估计大部分时刻在2.9附近，而箱粒子PHD多扩展目标跟踪方法数目估计大部分时刻在2.6附近；在26s～31s，真实目标数目为4，本发明数目估计大部分时刻在3.7附近，而箱粒子PHD多扩展目标跟踪方法数目估计大部分时刻在3.4附近。图4(b)中可以看到，在17s～25s，真实目标数目为3，本发明数目估计大部分时刻在2.8附近，而箱粒子PHD多扩展目标跟踪方法数目估计大部分时刻在2.4附近；在26s～31s，真实目标数目为4，本发明数目估计大部分时刻在3.6附近，而箱粒子PHD多扩展目标跟踪方法数目估计大部分时刻在3.0附近。

综合图4(a)和图4(b)，本发明和箱粒子PHD多扩展目标跟踪方法目标数目估计对比，在17s～25s，本发明目标数目估计幅度降低0.1，而箱粒子PHD多扩展目标跟踪方法目标数目估计幅度降低0.2；在26s～31s，本发明目标数目估计幅度降低0.1，而箱粒子PHD多扩展目标跟踪方法目标数目估计幅度降低0.4，可以得出，随着检测概率的降低，本发明目标数目估计幅度降低更小，即本发明对目标数目估计具有良好的稳定性。

Claims (10)

1.一种基于箱粒子CPHD的多扩展目标跟踪方法，包括如下步骤：

(1)将k时刻的多扩展目标量测集合Z_k划分成多个划分单元并利用包含函数[f]分别对多个划分单元进行包含，得到k时刻的量测箱粒子集其中k表示滤波时刻，m_k表示k时刻的划分单元个数，W_j表示第j个划分单元内量测分布的空间范围，表示k时刻的量测箱粒子状态；

(2)假设k-1时刻的存活箱粒子集为初始势分布为p_k-1(n)，将k-1时刻的新生箱粒子集并入存活箱粒子集得到k-1时刻的初始箱粒子集 其中和分别表示k-1时刻的存活箱粒子状态、新生箱粒子状态和初始箱粒子状态，和分别表示k-1时刻的存活箱粒子强度、新生箱粒子强度和初始箱粒子强度，N_k-1,s、N_k-1,b和N_k-1分别表示k-1时刻的存活箱粒子数目、新生箱粒子数目和初始箱粒子数目，n表示目标个数；

(3)对得到的k-1时刻的初始箱粒子集进行预测，得到k时刻的预测箱粒子集同时对假设的k-1时刻的初始势分布p_k-1(n)进行势分布预测，得到k时刻的预测势分布p_k|k-1(n)，其中，和分别表示k时刻的预测箱粒子状态和预测箱粒子强度；

(4)对k时刻的预测箱粒子集和预测势分布p_k|k-1(n)进行量测更新，实现步骤如下：

(4a)根据得到的k时刻的量测箱粒子集和k时刻的预测箱粒子集求得似然函数

(4b)利用求得的似然函数对得到的k时刻的预测箱粒子集进行量测更新，得到k时刻的更新箱粒子集

(4c)利用求得的似然函数对得到的k时刻的预测势分布p_k|k-1(n)进行量测更新，得到k时刻的势分布p_k(n)；

(5)对得到的k时刻的量测箱粒子集进行搜索，并判断是否搜索到对预测箱粒子强度贡献最大的量测箱[z_k]，若是，执行步骤(6)，否则，将得到的k时刻的更新箱粒子集输出；

(6)利用搜索到的量测箱[z_k]，对得到的k时刻的更新箱粒子集进行收缩，得到k时刻的收缩箱粒子集其中表示收缩箱粒子状态；

(7)对得到的k时刻的势分布p_k(n)进行估计，得到目标数目估计

(8)利用随机子划分方法，对得到的k时刻的收缩箱粒子集进行重采样，得到k时刻的重采样箱粒子集为其中和分别表示k时刻的重采样箱粒子状态和重采样箱粒子强度，N表示重采样次数；

(9)对得到的k时刻的重采样箱粒子状态集进行点化，得到k时刻的箱中心粒子状态集对该箱中心粒子状态集进行k-means聚类，得到目标的状态估计其中表示k时刻的箱中心粒子状态；

(10)判断k时刻是否为跟踪最后时刻，若是，则跟踪结束；否则保存k时刻的重采样箱粒子集并作为下个时刻的存活箱粒子集，令k＝k+1，执行步骤(2)。

2.根据权利要求1中所述的基于箱粒子CPHD的多扩展目标跟踪方法，其特征在于：步骤(2)中所述的k-1时刻新生箱粒子集是通过对上一时刻的量测箱粒子周围采样产生的，其状态服从高斯分布权值设置为其中，P_b表示箱粒子新生概率，Σ表示新生箱粒子状态分布协方差。

3.根据权利要求1中所述的基于箱粒子CPHD的多扩展目标跟踪方法，其特征在于：步骤(3)中所述的对得到的k-1时刻的初始箱粒子集进行预测，是通过***预测模型和强度预测方程实现的，其中，f_k|k-1表示状态转移函数，P_S表示箱粒子存活概率。

4.根据权利要求1中所述的基于箱粒子CPHD的多扩展目标跟踪方法，其特征在于：步骤(3)中所述的对假设的k-1时刻的初始势分布p_k-1(n)进行势分布预测，是通过马尔科夫转移矩阵和势分布预测方程实现的，其中，p_birth(n)表示从k-1时刻到k时刻有n个目标新生的概率，d表示目标消失概率。

5.根据权利要求1中所述的基于箱粒子CPHD的多扩展目标跟踪方法，其特征在于：步骤(4a)中所述的似然函数采用以下公式求得：

$L_k\left(\&lsqb;{\mathrm{zp}}_k^{\left(j\right)}\&rsqb;|\&lsqb;x_{k|k-1}^{\left(i\right)}\&rsqb;\right)=\frac{\left|\&lsqb;h_{cp}\&rsqb;\left(\&lsqb;x_{k|k-1}^{\left(i\right)}\&rsqb;,\&lsqb;{\mathrm{zp}}_k^{\left(j\right)}\&rsqb;\right)\right|}{\left|\&lsqb;x_{k|k-1}^{\left(i\right)}\&rsqb;\right|}$

其中，h_cp表示约束传播算法。

6.根据权利要求1中所述的基于箱粒子CPHD的多扩展目标跟踪方法，其特征在于：步骤(4b)中所述的利用求得的似然函数对得到的k时刻的预测箱粒子集进行量测更新，采用以下公式实现：

$\&lsqb;x_k^{\left(i\right)}\&rsqb;=\&lsqb;h_{cp}\&rsqb;(\&lsqb;x_{k|k-1}^{\left(i\right)}\&rsqb;,\&lsqb;{\mathrm{zp}}_k^{\left(j\right)}\&rsqb;)$

其中，D和分别表示目标被检测到和未被检测到情况的简记符号，P_D表示目标检测概率，L([Z_k])、L([Z_k]|D)和分别表示k时刻量测箱粒子集似然函数、目标被检测到的量测箱粒子集似然函数和目标未被检测到的量测箱粒子集似然函数。

7.根据权利要求1中所述的基于箱粒子CPHD的多扩展目标跟踪方法，其特征在于：步骤(4c)中所述的利用求得的似然函数对得到的k时刻的预测势分布p_k|k-1(n)进行量测更新，采用以下公式实现：

$p_k\left(n\right)=\frac{L(\&lsqb;Z_k\&rsqb;|n)}{L(\&lsqb;Z_k\&rsqb;)}{p}_{k|k-1}\left(n\right)$

其中，L([Z_k]|n)表示n个目标存在的量测箱粒子集似然函数。

8.根据权利要求1中所述的基于箱粒子CPHD的多扩展目标跟踪方法，其特征在于：步骤(5)中所述的搜索，采用以下公式实现：

$\&lsqb;z_k\&rsqb;=\arg \underset{w_{j,i}}{max}\{\&lsqb;{\mathrm{zp}}_k^{\left(j\right)}\&rsqb;,w_{j,i}>0\}$

其中，j＝1,…,m_k，

9.根据权利要求1中所述的基于箱粒子CPHD的多扩展目标跟踪方法，其特征在于：步骤(6)中所述的收缩，采用实现。

10.根据权利要求1中所述的基于箱粒子CPHD的多扩展目标跟踪方法，其特征在于：步骤(7)中所述的对得到的k时刻的势分布p_k(n)进行估计，采用最大后验估计方法实现。