CN104637070A - 基于概率假设密度的目标数变化的视频跟踪算法 - Google Patents

Abstract

本发明公开了一种基于概率假设密度的目标数变化的视频跟踪算法，包括以下步骤：采用运动检测方法，实现运动目标的快速有效分割产生目标量测；利用运动检测结果设计新生目标强度更新方法以产生新生目标强度函数；对PHD进行预测并更新；对更新后高斯元进行修剪；状态抽取，输出估计的目标数和状态；采用auction的航迹识别方法进行“航迹-状态估计”关联以识别目标航迹，输出带有身份标识的航迹。本发明能够有效跟踪数目变化的多目标，尤其能够跟踪位置未知的新生目标，并能够产生目标航迹。该算法改进了GM-PHD滤波的性能，能够在目标数发生变化和目标发生交叉情况下，取得较为鲁棒的估计状态，并能够识别目标航迹。

Description

基于概率假设密度的目标数变化的视频跟踪算法

技术领域

本发明涉及的是一种目标跟踪领域的算法，具体是基于概率假设密度的目标数变化的视频跟踪算法。

背景技术

视频目标跟踪广泛应用于军事和民用方面，如在军事方面的各种防卫系统：战场监视系统，机载火力控制系统等，民用方面包括汽车、个人GPS导航系统等。视频目标跟踪的目的是根据目标的特征如位置、速度、色彩、形状、纹理等以及目标描述模型在连续图像帧间建立对应关系，进而获得目标位置和形状等目标状态信息。一般视频目标跟踪方法可以分成两类：(1)数据驱动算法，直接建立目标模板和候选区域的相似度函数，利用匹配运算或优化算法将跟踪问题转化为局部的迭代寻优问题，如Mean-shift跟踪；(2)基于模型的概率跟踪算法，该类方法将目标跟踪问题转换为目标状态估计问题，目标状态的所有信息均包含于状态的后验密度中，目标状态后验密度可根据贝叶斯滤波来求解，如粒子滤波跟踪。

在单目标跟踪场景下，贝叶斯滤波首先预测目标状态的概率密度，然后利用当前帧的新量测对预测状态概率密度进行更新获得状态后验概率密度。在线性高斯条件下，贝叶斯滤波的闭合形式的递推方法等同于卡尔曼滤波。在非线性非高斯条件下，贝叶斯滤波的有效的实现方法是粒子滤波算法。

在多目标跟踪场景下，1994年Mahler首次系统地针对多传感器多目标跟踪问题提出了有限集统计理论，并以此为数学理论基础提出随机有限集框架下的多目标贝叶斯滤波。2000年Mahler结合点过程理论给出了Possion条件下近似的随机有限集的一阶矩递推滤波公式，从而得到了多目标贝叶斯滤波的近似即PHD滤波器。2005年，Ba-Ngu Vo等给出PHD滤波器在线性高斯条件下的解析解，GM-PHD方法即高斯混合PHD(GM-PHD)，已用于各种传感器的目标跟踪问题。GM-PHD借鉴了高斯混合滤波算法的思想，具有计算量小，强大的去除杂波能力，状态抽取可靠，易于处理多目标跟踪中目标的消亡和遮挡等优点。但该方法有两个缺点：(1)由于GM-PHD的输出为状态估计的随机集，无法给出目标身份(航迹)特征信息；(2)无法跟踪位置未知的新生目标。

发明内容

本发明针对现有技术中存在的上述不足，提出一种基于概率假设密度的目标数变化的视频跟踪算法。该算法能够有效跟踪数目变化的多目标，尤其能够跟踪位置未知的新生目标， 并能够产生目标航迹。该算法改进了GM-PHD滤波的性能，能够在目标数发生变化和目标发生交叉情况下，取得较为鲁棒的估计状态，并能够识别目标航迹。

本发明是通过以下技术方案实现的，本发明包括以下步骤：

第一步，采用运动检测方法，实现运动目标的快速有效分割，提取前景目标的质心，尺寸作为目标量测，并建立目标模型。

所述的目标模型包括目标状态模型和观测模型。

第二步，利用运动检测结果设计新生目标强度更新方法，用以产生新生目标强度函数。

第三步，使用步骤2得到的更新后的新生目标强度进行PHD预测。

所述的PHD预测，具体是：对新产生目标预测，对孵化目标预测以及对继续存在目标进行预测。

第四步，对检测到的目标进行更新，即使用步骤1得到的量测随机集对预测PHD进行更新。

所述的PHD的更新，具体是：对未检测到的目标预测高斯元的更新和由量测对检测到的目标预测高斯元的更新。

第五步，对以上输出的更新后高斯元进行修剪，主要任务是对于距离很近的高斯元进行合并，对极小权值的高斯元进行删除。

第六步，状态抽取，对于高斯元的权值大于阈值(一般取0.5)，其所对应的期望值作为目标状态，输出估计的目标数和状态。

第七步，采用auction的航迹识别方法进行“航迹-状态估计”关联以识别目标航迹，输出带有身份标识的航迹。

与现有的技术相比，本发明的有益效果是：该算法能够有效跟踪数目变化的多目标，尤其能够跟踪位置未知的新生目标，并能够产生目标航迹。该算法改进了GM-PHD滤波的性能，能够在目标数发生变化和目标发生交叉情况下，取得较为鲁棒的估计状态，并能够识别目标航迹。

附图说明

图1是目标数变化的目标跟踪系统框图。

图2是运动检测结果(第290,312,340,382,735,930帧)。

图3是原GM-PHD跟踪结果(第290,312,340,382,735,930帧)。

图4是本文算法跟踪结果(第290,312,340,382,735,930帧)。

图5是OSPA距离比较曲线。

图6是实施例中参数的设置。

具体实施方式

下面结合附图对本发明的实施例作详细说明：本实施例在以本发明技术方案为前提下进行实施，给出了详细的实施方式和具体的操作过程，但本发明的保护范围不限于下述的实施例。

实施例

本实施例采用行人检测视频数据库“OneShopOneWait2front”来实施。图6给出了实施例中主要参数的设置值。本实施例采用一个w×h矩形表示目标区域，则目标状态向量可以表示为，x_k＝(po_x,k,po_y,k,v_x,k,v_y,k,w,h)^T，其中po_k＝(po_x,k,po_y,k)为矩形的中心，v_k＝(v_x,k,v_y,k)为目标中心在图像的水平和竖直方向的速度，目标的一个观测向量为z_k＝(po_k,w,h)^T。假设k时刻描述后验PHD的高斯元参数为设k+1时刻量测随机集为 其中，J_k，分别为k时刻后验PHD的高斯元数目，第i个高斯元的均值，权值及相应的协方差。假设k时刻新生目标集为其高斯元参数为 ${\{{\mathrm{\ω}}_{\mathrm{\γ},k}^{\left(i\right)},m_{\mathrm{\γ},k}^{\left(i\right)},p_{\mathrm{\γ},k}^{\left(i\right)}\}}_{i=1}^{J_{\mathrm{\γ},k}}.$

第一步，采用运动检测方法计算目标位置，大小等参数，实现运动目标的快速有效分割，提取前景目标的质心坐标，长和宽作为目标量测向量，并建立目标状态模型和观测模型。

所述的目标状态模型采用线性高斯常速模型来建模，具体是：

x_k＝Fx_k-1+u_k             (1) 

其中，F为状态转移矩阵,u_k为零均值协方差为Q_v的高斯白噪声序列，

$F=\left(\begin{array}{ccc}{I}_2& A& 0_2\\ 0_2& I_2& 0_2\\ 0_2& 0_2& I_2\end{array}\right)---\left(2\right)$

$Q_v={\mathrm{\σ}}_v^2=\left(\begin{array}{ccc}\frac{{\mathrm{\Δ}}^4}{4}{I}_2& \frac{{\mathrm{\Δ}}^3}{2}{I}_2& 0_2\\ \frac{{\mathrm{\Δ}}^3}{2}{I}_2& {{\mathrm{\Δ}}^{2}I}_2& 0_2\\ 0_2& 0_2& \mathrm{\Δ}\end{array}\right)---\left(3\right)$

式中I_n和0_n分别表示n×n单位矩阵和零矩阵，Δ为采样周期，σ_v为系统噪声标准差，

$A=\left(\begin{array}{cc}\mathrm{\Δ}& 0\\ 0& \mathrm{\Δ}\end{array}\right).$

所述的观测模型采用线性高斯观测模型，具体是：

z_k＝Hx_k+w_k       (4) 

其中，H为观测矩阵，w_k为零均值协方差矩阵为Q_w的高斯白噪声。H和Q_w定义如下，其中 ${\mathrm{\σ}}_w^2=[{\mathrm{\σ}}_x^2,{\mathrm{\σ}}_y^2,{\mathrm{\σ}}_w^2,{\mathrm{\σ}}_h^2],$ σ_w为观测噪声标准差。

$H=\left(\begin{array}{ccc}{I}_2& 0_2& 0_2\\ 0_2& 0_2& I_2\end{array}\right)---\left(5\right)$

$Q_w={\mathrm{\σ}}_w^2{I}_4---\left(6\right)$

第二步，利用运动检测结果设计新生目标强度更新方法，用以产生新生目标强度函数。

第三步，使用步骤2得到的更新后的新生目标强度进行PHD预测。

所述的PHD预测，具体是：对新产生目标预测，对孵化目标预测以及对继续存在目标进行预测。

第四步，对检测到的目标进行更新，即使用步骤1得到的量测随机集对预测PHD进行更新。

所述的PHD的更新，具体是：对未检测到的目标预测高斯元的更新和由量测对检测到的目标预测高斯元的更新。

第五步，对以上输出的更新后高斯元进行修剪，主要任务是对于距离很近的高斯元进行合并，对极小权值的高斯元进行删除。

第六步，状态抽取，对于高斯元的权值大于阈值(一般取0.5)，其所对应的期望值作为目标状态，输出估计的目标数和状态。

第七步，采用auction的航迹识别方法进行“航迹-状态估计”关联以识别目标航迹，输出带有身份标识的航迹T_k+1.注意，对于状态数与航迹数不等(N_k+1>L_k)时，X_k+1可能包含新生目标，则需要计算X_k+1中未关联的状态与步骤2得到的准新生目标的巴氏距离，保留巴氏距离小的状态为新目标航迹增加至T_k+1中。

所述的新目标强度更新方法，具体是：

(1)遍历计算与上一帧量测集的巴氏距离，在连续5帧图像中重复上述操作，根据延时决策，保留下的量测记为

(2)遍历计算与上一帧航迹的巴氏距离，距离小的量测可能由上一帧幸存的目标产生，因此，剔除中由幸存目标产生的量测，保留准新目标量测，记为以准新生目标量测值构造新生目标高斯元，并更新新生目标强度函数为  ${\{{\mathrm{\ω}}_{\mathrm{\γ},k+1}^{\left(i\right)},m_{\mathrm{\γ},k+1}^{\left(i\right)},p_{\mathrm{\γ},k+1}^{\left(i\right)}\}}_{i=1}^{J_{\mathrm{\γ},k+1}}.$

所述的对新生目标进行预测，具体是：

对j＝1,…,J_g,k+1，i＝1,…,J_g,k+1有 ${\mathrm{\ω}}_{k+1|k}^{\left(i\right)}={\mathrm{\ω}}_{\mathrm{\γ},k+1}^{\left(j\right)},m_{k+1|k}^{\left(i\right)}=m_{\mathrm{\γ},k+1}^{\left(j\right)},P_{k+1|k}^{\left(i\right)}=P_{\mathrm{\γ},k+1}^{\left(i\right)},$ 即使用步骤2得到的更新后的新生目标强度预测PHD参数；

所述的对孵化目标进行预测，具体是：

对j＝1,…,J_b,k+1，l＝1,…,J_k，i＝J_g,k+1,…,(J_g,k+1+J_kJ_b,k+1)按照下式预测其权值，均值及协方差，

$m_{k+1|k}^{\left(i\right)}=d_{\mathrm{\β},k}^{\left(j\right)}{+F}_{\mathrm{\β},k}^{\left(j\right)}{m}_k^{\left(l\right)}---\left(7\right)$

$P_{k+1|k}^{\left(i\right)}=Q_{\mathrm{\β},k}^{\left(j\right)}+F_{\mathrm{\β},k}^j{P}_k^{\left(l\right)}{\left(F_{\mathrm{\β},k}^j\right)}^T---\left(8\right)$

${\mathrm{\ω}}_{k+1|k}^{\left(i\right)}={\mathrm{\ω}}_k^{\left(l\right)}{\mathrm{\ω}}_{\mathrm{\β},k+1}^{\left(j\right)}---\left(9\right)$

其中，分别为(1)式中的状态转移矩阵及孵化目标状态噪声协方差。

所述的对继续存在目标进行预测，具体是：

设其生存概率为p_S，则对j＝1,…,J_k，

i＝(J_g,k+1+J_kJ_b,k+1),…,(J_g,k+1+J_kJ_b,k+1+J_k)按照以下公式预测其权值，均值及协方差，

${\mathrm{\ω}}_{k+1|k}^{\left(i\right)}=p_S{\mathrm{\ω}}_k^{\left(j\right)}---\left(10\right)$

$m_{k+1|k}^{\left(i\right)}={F_{k}m}_k^{\left(j\right)}---\left(11\right)$

$P_{k+1|k}^{\left(i\right)}=Q_k^{\left(j\right)}+F_k{P}_k^{\left(j\right)}{\left(F_k\right)}^T---\left(12\right)$

所述的PHD更新，具体是：

利用量测Z_k+1，及(4)式的观测矩阵H和量测噪声R更新。设p_D为检测概率，对未检测到的目标用下式更新：对j＝1,…,J_k+1|k(其中J_k+1|k＝J_g,k+1+J_kJ_b,k+1+J_k)有,

${\mathrm{\ω}}_{k+1}^{\left(j\right)}={(1-p_D)\mathrm{\ω}}_{k+1|k}^{\left(j\right)}---\left(13\right)$

$m_{k+1}^{\left(j\right)}=m_{k+1|k}^{\left(j\right)}---\left(14\right)$

$P_{k+1}^{\left(j\right)}=P_{k+1|k}^{\left(j\right)}---\left(15\right)$

对检测到的目标更新，即使用第一步得到的量测随机集对预测PHD进行更新，对每一个z∈Z_k+1且j＝1,…,J_k+1|k，更新预测PHD得后验PHD高斯元的权值，均值及协方差为：

${\mathrm{\ω}}_{k+1}^{\left(j\right)}\left(z\right)=\frac{p_D{\mathrm{\ω}}_{k+1|k}^{\left(j\right)}N(z;H{m}_{k|k+1}^{\left(j\right)},\mathrm{nR}+H{P}_{k+1|k}^j{H}^T)}{{\mathrm{\κ}}_k\left(z\right)+p_D{\mathrm{\Σ}}_{l=1}^{J_{k+1|k}}{\mathrm{\ω}}_{k+1|k}^{l}N(z;H{m}_{k|k+1}^{\left(l\right)},\mathrm{nR}+H{P}_{k+1|k}^l{H}^T)}---\left(16\right)$

$m_{k+1}^{\left(j\right)}=m_{k+1|k}^{\left(j\right)}+P_{k+1|k}^{\left(j\right)}{H}^T{[R+H{P}_{k+1|k}^{\left(j\right)}{H}^T]}^{-1}(z-H{m}_{k+1|k}^{\left(j\right)})---\left(17\right)$

$P_{k+1}^{\left(j\right)}=P_{k+1|k+1}^{\left(j\right)}---\left(18\right)$

其中，κ_k(z)为泊松分布的杂波RFS的强度函数。 

本实施例中，基于auction的航迹识别算法具体操作过程为：

所述的基于auction的航迹识别算法是为了克服PHD滤波的结果是不带有目标身份或航迹特征的目标状态随机集的缺点。本实施例中，将目标航迹识别问题可以看作二维分配问题，即将当前k+1时刻的目标状态估计X_k+1分配给上一个时刻k已确认的目标航迹T_k，而最优二维分配的求解可转化为求取关联代价函数和的极小值。假设X_k+1包含N_k+1个状态估计T_k包含L_k个航迹则“航迹-状态估计”关联的目标函数可以定义为，

$C=\min {\mathrm{\Σ}}_{m=1}^{N_{k+1}}{\mathrm{\Σ}}_{n=1}^{L_k}{c}_{m,n}{a}_{m,n}---\left(19\right)$

满足约束条件： ${\mathrm{\Σ}}_{m=1}^{N_{k+1}}{a}_{m,n}=1,n=\mathrm{1,2},\·\·\·,L_k---\left(20\right)$

${\mathrm{\Σ}}_{n=1}^{L_k}{a}_{m,n}=1,m=\mathrm{1,2},\·\·\·,N_{k+1}---\left(21\right)$

其中，a_m,n为二值变量，状态估计与航迹关联时取1，否则为0，c_m,n为二者相关联的代价。

本实施例设计的代价函数的原理为采用巴氏距离衡量目标的位置和大小的相似性，相似性越大则两个目标区域的距离越小，关联代价越小，反之越大。假设k+1时刻的一个目标 状态为 $x_{k+1}^i={({\mathrm{po}}_{x,k+1}^i,p{o}_{y,k+1}^i,v_{x,k+1}^i,v_{y,k+1}^i,w_{k+1}^i,h_{k+1}^i)}^T,$ k时刻的一条航迹为  $t_k^j={({\mathrm{po}}_{x,k}^j,v_{x,k}^j,p{o}_{y,k}^j,v_{y,k}^j,w_k^j,h_k^j)}^T,$ 和的相似度可以通过即巴氏距离来度量，

$d(x_{k+1}^i,t_k^j)=\sqrt{1-\mathrm{\ρ}[x_{k+1}^i,t_k^j]}---\left(22\right)$

其中为Bhattacharyya系数，A_s的取值分别为向量中的位置和尺寸值 $\{p{o}_{x,k+1}^i,p{o}_{y,k+1}^i,w_{k+1}^i,h_{k+1}^i\},$ B_s取值为 $\{p{o}_{x,k}^j,p{o}_{y,k}^j,w_k^j,h_k^j\}.$

则定义和相关联的代价c_i,j为，

$c_{i,j}=\exp (-d(x_{k+1}^i,t_k^j))---\left(23\right)$

当确立了代价函数，关联目标函数可通过auction优化算法求解。

1)计算关联矩阵，根据式(23)遍历计算所有状态估计与确认航迹的关联代价，生成N_k+1×L_k的二维关联矩阵。

2)初始化X_k+1中所有未分配成功的状态估计并将航迹价格(price)设定为0。

3)根据公式(24)对任一未分配成功的状态找到与其对应的“最佳”航迹如果所有状态分配成功，则结束计算转到步骤7。

$c_{\mathrm{mn}}-P_n=\underset{j=1,\·\·\·,L\left(k\right)}{\max }(c_{\mathrm{mj}}-P_j)---\left(24\right)$

4)解除上一个循环计算后分配给航迹的状态而将分配给

5)按照下式更新航迹的价格

P_n＝P_n+d_n+ε       (25) 

其中，d_n为状态估计的“最佳”和“次最佳”分配代价的差值，而ε为一个设定常数或函数用以防止死循环。

6)返回至步骤3。

7)输出航迹T_k+1。

实施效果

图2给出了本实施例运动检测算法得到的检测后的二值化图像，目标量测(质心位置和尺寸)由检测图中光斑计算得到。图3给出了GM-PHD直接对运动检测得到的量测滤波后 的结果，其中“*”指示了目标质心位置，矩形框为目标尺度大小。在使用量测对预测PHD的高斯元进行更新后，更新后的高斯元只有当其权值大于状态抽取的阈值时方能保留，而小于该阈值的高斯元均会被删除。当新生目标强度函数未知时，由新生目标量测更新后的高斯元的权值会很小而被当作杂波删除，进而发生丢跟问题。因此原GM-PHD在新生目标的图像(见图3第312，340,382，930帧)中出现了丢跟问题。此外，原GM-PHD输出的是目标状态随机集，无法给目标航迹，图3中没有目标航迹标识，仅指示了目标位置和尺度估计。图4给出了本实例的算法的跟踪结果，可见本实施例的算法既能跟踪位置未知的新生目标，又能够给出目标航迹。该图中“*”，“o”，“Δ”和“”分别给出了四个目标的航迹，确认了目标在视频序列前后帧的身份关系。

OSPA距离可以实现目标状态估计与真实状态在最优配对情况下两个集合之间的距离，常用来衡量多目标跟踪误差。图5给出了本实施例的算法和原GM-PHD的OSPA跟踪误差曲线。OSPA对目标数估计错误有严重判罚，即当目标数估计与真值有偏差时，OSPA距离会表现出显著增大态势。从图5中可以看出本实施例算法产生的OSPA距离比原GM-PHD显著降低，算法的误差较小。可见，本实施系统改进了GM-PHD滤波的性能，能够在目标数发生变化和目标发生交叉情况下，取得较为鲁棒的估计状态，并能够识别目标航迹。

Claims (8)

1.一种基于概率假设密度的目标数变化的视频跟踪算法，其特征在于，包括以下步骤：

第一步，采用运动检测方法，实现运动目标的快速有效分割，提取前景目标的质心，尺寸作为目标量测，并建立目标模型；

第二步，利用运动检测结果设计新生目标强度更新方法，用以产生新生目标强度函数；

第三步，使用步骤2得到的更新后的新生目标强度进行PHD预测；

第四步，对检测到的目标进行更新，即使用步骤1得到的量测随机集对预测PHD进行更新；

第五步，对以上输出的更新后高斯元进行修剪，主要任务是对于距离很近的高斯元进行合并，对极小权值的高斯元进行删除；

第六步，状态抽取，对于高斯元的权值大于阈值(一般取0.5)，其所对应的期望值作为目标状态，输出估计的目标数和状态；

第七步，采用auction的航迹识别方法进行“航迹-状态估计”关联以识别目标航迹，输出带有身份标识的航迹。

2.根据权利要求1所述的基于概率假设密度的目标数变化的视频跟踪算法，其特征是，所述的目标模型包括目标状态模型和观测模型，具体是，

x_k＝Fx_k-1+u_k，

z_k＝Hx_k+w_k，

其中，F为状态转移矩阵,u_k为零均值协方差为Q_v的高斯白噪声序列，H为观测矩阵，w_k为零均值协方差矩阵为Q_w的高斯白噪声。

3.根据权利要求1所述的基于概率假设密度的目标数变化的视频跟踪算法，其特征是，所述的更新的新生目标强度更新方法，具体是：计算当前帧量测和上一帧量测的巴氏距离，保留连续5帧图像中巴氏距离小的量测计算与上一帧航迹的巴氏距离，选择巴氏距离大的量测为准新生目标量测并由准新生目标量测构造新生目标高斯元，与上一帧的新生目标高斯元合并以得到更新的新生目标强度。

4.根据权利要求1所述的基于概率假设密度的目标数变化的视频跟踪算法，其特征是，所述的PHD预测，具体是：对更新后的新生目标预测，对孵化目标预测以及对继续存在目标进行预测。

5.根据权利要求1所述的基于概率假设密度的目标数变化的视频跟踪算法，其特征是，所述的PHD的更新，具体是：对未检测到的目标预测高斯元的更新和由量测对检测到的目标预测高斯元的更新。

6.根据权利要求1所述的基于概率假设密度的目标数变化的视频跟踪算法，其特征是，所述的基于auction的航迹识别方法，具体是：将目标航迹识别问题可以看作二维分配问题，即将当前k+1时刻的目标状态估计X_k+1分配给上一个时刻k已确认的目标航迹T_k，而最优二维分配的求解可转化为求取关联代价函数和的极小值。

7.根据权利要求6所述的基于auction的航迹识别方法，其特征是，所述的关联代价函数采用巴氏距离衡量目标的位置和大小的相似性来建立，具体是，

$c_{i,j}=\exp (-d(x_{k+1}^i,t_k^j))$

其中，c_i,j为航迹和状态相关联的代价，为巴氏距离，为Bhattacharyya系数，A_s的取值分别为向量中的位置和尺寸值B_s取值为

8.根据权利要求6所述的基于auction的航迹识别方法，其特征是，所述的最优二维分配的求解采用auction优化算法，具体是，遍历计算所有状态估计与确认航迹的关联代价，生成N_k+1×L_k的二维关联矩阵，并利用auction优化算法实现当前帧状态和确认航迹的最优配对，从而完成当前帧多目标状态的航迹识别。