一种联合概率数据关联的视频多目标快速跟踪方法

摘要

本发明针对监控范围较大、目标外观特征少的视频多目标数据关联及跟踪问题，仅利用目标运动特征，提出了一种基于改进的联合概率数据关联JPDA的复杂情况下视频多目标快速跟踪方法，即一种联合概率数据关联的视频多目标快速跟踪方法。本方法采用简化的穆尔蒂算法求JPDA的最优K个联合事件，大大降低了计算复杂度；根据JPDA的关联概率讨论目标的运动情况，分析在多目标新出现、遮挡、消失、分离等复杂情况下当前帧量测与跟踪目标的数据关联问题，获取复杂运动的多目标跟踪轨迹。本发明提出的方法能实现较大监控范围下视频多目标的快速跟踪，并能大大提高跟踪性能。

说明书

技术领域

本发明属于视频多目标跟踪技术领域，涉及一种联合概率数据关联的视频多目标快速跟踪方法。

背景技术

在多运动目标的视频监控中，需要根据检测结果进行数据关联来匹配、识别出连续帧间的不同运动目标，从而实现多目标的跟踪。视频监控中往往需要对较大范围的场景进行监控，常出现运动目标外观特征相似或目标区域较小的情况，此时数据关联仅能依靠目标的运动特征完成，而传统的数据关联方法应用于视频多目标跟踪还存在很多问题。

目前解决多运动目标的数据关联问题，主要的方法有：1)根据最近邻方法，计算落入跟踪门限内的量测，判定目标的运动情况，并根据有效量测直接估计、更新当前运动状态，这种方法计算量小，但在目标数目较多、运动情况复杂时，抗干扰能力差，容易产生错误关联；2)杂波环境下多目标数据关联方法如联合概率数据关联(JointProbability Data Assocaition，JPDA)、多假设跟踪(Multiple Hypothesis Tracking，MHT)，目前多用于红外、雷达小目标或机动目标跟踪系统，虽然和视频监控跟踪有相似性，但须满足一对一关联的约束，而在视频监控跟踪系统中，多运动目标常发生新出现、消失、遮挡、分离等复杂运动情况，即出现一对多或多对一的关联情况，因此需进一步研究这些经典的多目标数据关联方法如何应用于视频多目标跟踪；3)采用优化算法分析当前帧检测区域与跟踪目标的最优关联，如采用的图优化、赋权二分图等方法，但此类方法通常需要获取遮挡或分离区域外观特征进行匹配优化计算，在目标区域小、外观信息少的情况下，难以得到目标外观特征，从而使得关联准确性低。

发明内容

本发明所要解决的技术问题是提供一种联合概率数据关联的视频多目标快速跟踪方法，能提高跟踪的有效性、可靠性和实时性。

本发明是通过以下技术方案实现的：

一种联合概率数据关联的视频多目标快速跟踪方法，其特征在于，包括以下步骤：

步骤1：根据当前帧图像的检测结果产生确认矩阵；

步骤2：采用简化的murty算法获得确认矩阵对应的联合事件及参数；

步骤3：计算k时刻所有量测的联合事件的条件概率；k的取值从跟踪开始时刻计算；

步骤4：计算量测与目标t的关联概率β_k^j，t，β_k^0，t；用于评估量测与目标关联的可能性大小；β_k^j，t为k时刻有效量测j和目标t的关联概率，β_k^0，t为k时刻虚量测和目标t的关联概率；

步骤5：通过Kalman滤波器，得到目标的状态估计及协方差阵：根据步骤4计算的关联概率，分析目标的运动状态及与上一帧目标的关联情况，包括五种情况：正常、新出现、消失、遮挡、分离，并根据上述五种情况得到目标的状态估计值，即获得目标在当前帧x、y方向的位置和速度，完成视频多目标跟踪；协方差阵是目标的估计状态和真实状态的偏差，表征跟踪的精确性；

步骤5完成后返回到步骤1对下一帧图像进行跟踪；

所述步骤1的确认矩阵为Ω；表示当前帧有效量测与目标跟踪门间的关系，定义为：

$\Omega ={\left\{{\omega }_j^t\right\}}_{j=\mathrm{1,2},...,m}^{t=\mathrm{0,1,2},...N};$

式中，t表示目标，j表示量测，ω_j^t表示目标t与量测j的关系；N表示跟踪目标的个数，m表示量测的个数，确认矩阵Ω的行表示当前帧量测，列表示跟踪目标，其中${\omega }_j^t=1$表示当前帧量测j落入跟踪目标t的跟踪门内；${\omega }_j^t=0$表示量测j没有落入目标t的跟踪门内，并令t＝0表示虚警；对应Ω第一列的所有元素为1，用于表示当前帧任一量测都有可能源于当前检测错误产生的虚警；

所述步骤2的具体步骤为：

包括如下步骤：

步骤a：根据确认矩阵，构建目标、量测的赋权二分图，得到待匹配矩阵

步骤b：应用匈牙利算法求待匹配矩阵对应的赋权二分图中目标、量测的最优匹配，并得到匹配结果矩阵；

步骤c：如未在步骤a的待匹配矩阵中增加虚目标或虚量测，则转步骤d；如增加了，则在匹配结果矩阵中删去，得到确认矩阵的一个最优匹配关系，即获得确认矩阵的一个可行矩阵；进入下步骤d；

步骤d：将当前步骤a中求得的待匹配矩阵减去当前的匹配结果矩阵，得到新的需优化的待匹配矩阵；

重复执行步骤b到步骤d K次，得到对应于确认矩阵Ω的K个最优可行矩阵，即得到确认矩阵对应的联合事件；

再得到对应的联合事件θ_k，i中的参数：假量测数Φ(θ_k，i)、量测关联指示器τ_j(θ_k，i)及目标检测指示器δ_t(θ_k，i)；

量测关联指示器：${\tau }_j\left({\theta }_{k,i}\right)=\left(\begin{array}{c}1,t_j>0\\ 0,t_j=0\end{array}\right),j=\mathrm{1,2},...,m_k$，其中t_j是联合事件θ_k，i对应的可行矩阵中与量测j关联的目标的取值，即τ_j(θ_k，i)表明量测j是否和一个真实目标关联；

目标检测指示器：，其中N表示目标数；则δ_t(θ_k，i)表示在联合事件θ_k，i对应的可行矩阵中目标t是否被检测到；

假量测数Φ(θ_k，i)表示联合事件θ_k，i对应的可行矩阵中假量测的数目，根据量测关联指示器得到：$\Phi \left({\theta }_{k,i}\right)=\underset{j=1}{\overset{m_k}{\Sigma }}[1-{\tau }_j\left({\theta }_{k,i}\right)];$

所述步骤3中，k时刻所有量测的联合事件的条件概率为：

$P\left({\theta }_{k,i}\right|Z^k)=\frac{\Phi \left({\theta }_{k,i}\right)!}{{\mathrm{cV}}^{\Phi \left({\theta }_{k,i}\right)}}\underset{j=1}{\overset{m_k}{\Pi }}{\left({\Lambda }_{k,j}\right)}^{{\tau }_j\left({\theta }_{k,i}\right)}\underset{t=1}{\overset{n}{\Pi }}{\left(P_D^t\right)}^{{\delta }_t\left({\theta }_{k,i}\right)}{(1-P_D^t)}^{1-{\delta }_t\left({\theta }_{k,i}\right)};$

其中，c为归一化常数，V表示跟踪门体积，Λ_k，j表示滤波残差似然函数，P_D^t表示目标t的检测概率；

检测概率P_D^t为常数，$P_D^t=0.98;$

滤波残差似然函数为：

${\Lambda }_{k,j}=\frac{1}{{\left(2\pi \right)}^{w/2}{\left|S_k^{t_j}\right|}^{1/2}}\exp \left[{(z_{k,j}-{\hat{z}}_{k|k-1}^{t_j})}^T{S}_k^{{t_j}^{-1}}(z_{k,j}-{\hat{z}}_{k|k-1}^{t_j})\right];$

其中代表目标t_j的预测位置：${\hat{z}}_{k|k-1}^{t_j}=C{\hat{x}}_{k|k-1},$$C=\left(\begin{array}{cccc}1& 0& 0& 0\\ 0& 0& 1& 0\end{array}\right);$是Kalman滤波器中一步预测状态向量：${\hat{x}}_{k|k-1}=A{x}_{k-1|k-1}$$A=\left(\begin{array}{cccc}1& \mathrm{\Delta t}& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& \mathrm{\Delta t}\\ 0& 0& 0& 1\end{array}\right),$是状态向量x在k-1时刻的估计值，与k时刻的状态向量估计值的计算相同，代表相应于目标t_j的残差协方差矩阵，z_k，j表示k时刻第j个量测，从当前帧(k时刻)的检测结果中读取；残差协方差矩阵S_k定义为：$S_k=C_k{P}_{k|k-1}{C}_k^T$；各目标的残差协方差矩阵计算相同，即的计算与S_k相同；

P_k|k-1表示Kalman滤波器中根据k-1时刻的协方差阵预测得到k时刻的预测协方差阵：

P_k|k-1＝AP_k-1|k-1A^T+Q_k-1；初始值$P_{0|0}=\left(\begin{array}{cccc}0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right);$

其中矩阵Q_k-1表示k-1时刻的值，与k时刻的值计算相同，k时刻的值为：$Q_k=G_k{Q}^{\prime }{G}_k^T$，G_k也为k时刻的值，定义为：$G_k=\left(\begin{array}{cc}\frac{k^2}{2}& 0\\ k& 0\\ 0& \frac{k^2}{2}\\ 0& k\end{array}\right),$Q′为常数，取：

$Q^{\prime }=\left(\begin{array}{cc}4& 0\\ 0& 4\end{array}\right);$

P_k-1|k-1是k-1时刻的状态估计协方差阵；

跟踪门体积V：其中β为常数，β＝9.5；

所述步骤4的β_k^j，t，β_k^0，t计算公式如下：

${\beta }_k^{j,t}=P\left({\theta }_k^{j,t}\right|Z^k)=\underset{i=1}{\overset{n_k}{\Sigma }}P\left({\theta }_{k,i}\right|Z^k){\hat{\omega }}_j^t\left({\theta }_{k,i}\right),j\in [1,m_k];$

${\beta }_k^{0,t}=1-\underset{j=1}{\overset{m_k}{\Sigma }}{\beta }_k^{j,t};$

n_k表示k时刻联合事件的个数，即可行矩阵的个数，m_k表示有效量测的个数；

所述步骤5的具体步骤为：

由步骤4中计算的关联概率得到k时刻目标t状态估计：

${\hat{x}}_{k|k}^t=\underset{j=0}{\overset{m_k}{\Sigma }}{\beta }_k^{j,t}{\hat{x}}_{k|k,j}^t=\underset{j=0}{\overset{m_k}{\Sigma }}{\beta }_k^{j,t}(A{\hat{x}}_{k-1|k-1}^t+K_k^t(z_k-\mathrm{CA}{\hat{x}}_{k-1|k-1}^t));$

K_k^t为目标t在k时刻的Kalman滤波增益矩阵，$K_k^t=P_{k|k-1}^t{C}^T{[{\mathrm{CP}}_{k|k-1}^t{C}^T+R]}^{-1},$P_k|k-1^t表示目标t的一步预测协方差矩阵，各目标的预测协方差矩阵计算相同，即P_k|k-1^t与前文中的P_k|k-1的计算相同；R是常数，取为：R＝[0.150 0.150]；表示目标t在上一帧的状态估计；初始值${\hat{x}}_{0|0}^t=\left(\begin{array}{cccc}{p}_{x_0}& 0& p_{y_0}& 0\end{array}\right),$其中分别表示目标初始位置的x方向、y方向的值；

以下分别针对五种情况说明：

1)正常：无须修正；

2)新出现：新目标的初始状态估计为${\hat{x}}_{k|k}^t={\left(\begin{array}{cccc}{p}_x& 0& p_y& 0\end{array}\right)}^T;$p_x、p_y分别表示当前帧检测结果中目标位置在x方向、y方向的值，所述的当前帧即k时刻对应帧；

3)消失：终止该目标的跟踪；

4)遮挡：遮挡目标t₁、t₂状态估计值和修正为：${\hat{x}}_{k|k}^{t_1}={\left(\begin{array}{cccc}{p}_x& 0& p_y& 0\end{array}\right)}^T,$${\hat{x}}_{k|k}^{t_2}={\left(\begin{array}{cccc}{p}_x& 0& p_y& 0\end{array}\right)}^T;$

5)分离：表示一个目标分成2个或多个，无须修正状态估计值；根据目标t的状态估计，计算状态估计协方差阵P_k|k^t，以进行下一帧跟踪：

$P_{k|k}^t=P_{k|k-1}^t-(1-{\beta }_k^{0,t})K_k^t{S}_k^t{\left(K_k^t\right)}^T+\underset{j=0}{\overset{m_k}{\Sigma }}{\beta }_k^{j,t}[{\hat{x}}_{k|k,j}^t{\left({\hat{x}}_{k|k,j}^t\right)}^T-{\hat{x}}_{k|k}^t{\left({\hat{x}}_{k|k}^t\right)}^T];$

P_k|k-1^t表示目标t的预测协方差阵，计算方法同前述的P_k|k-1的计算方法；式中K_k^t、S_k^t分别为目标t在k时刻的Kalman滤波增益矩阵、残差协方差矩阵。

当没有任何量测源于目标t时，计算对于没有任何量测源于目标t的状态估计协方差阵为：

$P_{k|k,0}^t=P_{k|k-1}^t;$

有益效果：

与现有技术相比，本发明的优越性体现在：

1、仅利用目标的运动特征，通过改进的联合概率数据关联(JPDA)算法实现监控视频中多运动目标的快速跟踪；

2、为避免传统JPDA算法在目标、量测较多时的关联匹配呈指数增加，提出一种简化的murty算法【中文：穆尔蒂】，快速计算JPDA算法中确认矩阵对应的最优K个联合事件，能大大降低计算复杂度，显著提高跟踪算法的实时性；

3、传统JPDA算法产生可行事件需满足量测、目标“一对一”的关联约束，并由此计算关联概率、得到目标的状态估计。而在视频监控中，多目标存在新出现、消失、遮挡、分离等运动情况，量测和目标的关联情况复杂，因此，必须对传统JPDA算法中计算目标状态估计值进行改进、修正。本发明提出根据JPDA算法中的关联概率判定目标的运动情况，分析在目标正常、新出现、消失、遮挡、分离(前景检测不准确造成目标碎片)等复杂情况下当前帧量测与跟踪目标的数据关联问题，并据此推导目标状态估计，实现复杂运动情况下视频多目标的跟踪，能大大提高跟踪的有效性、可靠性。

附图说明：

图1为本发明的总体流程图；

图2为第k帧量测与目标跟踪门的关系示意图；

图3为图1对应的赋权二分图；

图4为图3增加一虚量测的示意图；

图5为目标的运动情况示意图，a)新出现，b)遮挡，c)分离

图6为不同方法在监控视频中的跟踪误差曲线。【a)监控场景；b)传统JPDA算法的跟踪误差曲线；d)本发明的方法的跟踪误差曲线】说明：由于本发明主要针对传统JPDA算法进行了改进，因此，改为误差与传统JPDA对比即可。

具体实施方式

实施例1：

本发明提出的视频多目标跟踪方法的基本思想是：首先，采用文献【中文题目：采用自适应核密度估计的基于运动的背景减除法，作者：Mittal.A，Paragios.N.英文题目：Motion-Based Background Subtraction using Adaptive Kernel Density Estimation，发表刊物：In：Proceeding of IEEE Conference on Computer Vision and Pattern Recognition，Washington DC，USA：IEEE Pres，2004，302～309】提出的背景检测方法【背景减除法是一种常用的已知算法】检测当前帧的目标，并将检测出的目标用其外接矩形框表示。由于在监控范围较大、目标外观特征少的情况下，仅能利用目标运动特征进行跟踪。本发明中，将视频中所有目标的运动特性认为是服从某种运动模型，在建立运动模型后，只须在连续帧中通过运动模型匹配(也称数据关联)计算出各个目标的运动向量值，即可完成视频多目标的跟踪。

选取目标外接矩形框中心点的运动特征建立运动模型，用于代表目标的运动特征，具体方式如下：

视频监控中，可认为目标运动过程和观测过程都是线性的，因此我们只针对线性情况分析目标运动特性，建立目标运动模型，线性运动的状态方程和测量方程为：

x_k＝Ax_k-1+w_k              (1)

z_k＝Cx_k+v_k                (2)

公式(1)和(2)为针对所有的目标的运动模型。

式(1)中，w为系统噪声反映线性系统模型的精确程度，具有均值为零的高斯分布，其协方差矩阵为Q，k时刻的值取为：$Q_k=G_k{Q}^{\prime }{G}_k^T$，G_k也为k时刻的值，定义为：$G_k=\left(\begin{array}{cc}\frac{k^2}{2}& 0\\ k& 0\\ 0& \frac{k^2}{2}\\ 0& k\end{array}\right),$Q′为常数，取：$Q^{\prime }=\left(\begin{array}{cc}4& 0\\ 0& 4\end{array}\right);$

v为观测噪声，是均值为零的白噪声序列，和w互不相关，其协方差矩阵为R，是常数，根据某像素点在视频中的方差确定，实验视频中取为：R＝[0.150 0.150]。

同时，定义状态向量x_k＝[p_x，v_x，p_y，v_y]^T和观测向量z_k＝[p_x，p_y]^T，其中p_x、v_x、p_y、v_y分别表示第k帧目标外接矩形框中心点x方向坐标值、速度值及y方向坐标值、速度值。

由于在视频中，通常相邻帧间时间间隔Δt很小，目标运动可近似认为是匀速运动，则根据匀速运动动力学方程得到式(1)中的状态转移矩阵A及式(2)中的观测矩阵C：

$A=\left(\begin{array}{cccc}1& \mathrm{\Delta t}& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& \mathrm{\Delta t}\\ 0& 0& 0& 1\end{array}\right)$$C=\left(\begin{array}{cccc}1& 0& 0& 0\\ 0& 0& 1& 0\end{array}\right)---\left(3\right)$

A、C为常量。建立如上所述的运动模型后，再通过匹配各个目标在上一帧和当前帧的运动模型，计算出当前帧中各目标的状态向量x_k，则可根据状态向量中的p_x、p_y分量，得到该目标外接矩形框中心点在当前帧的位置，从而通过连续帧运动状态的关联性，即匹配程度完成多个目标在视频中的跟踪，并获得各个目标的运动轨迹。

因此，在建立目标运动模型后，需要匹配帧间不同目标的运动特征来完成跟踪，这是数据关联问题。本发明提出基于改进的联合概率数据关联算法跟踪视频多目标，包括5个步骤，具体如下：

1、产生确认矩阵：一本步骤中需确认当前帧的检测结果中哪些是对应目标的，即确定哪些检测信息是有效的。

由于当前帧检测到的可能包括目标、干扰(如光线变化、树叶晃动等背景中的动态变化)，则本步骤中需确认当前帧的检测结果中哪些是对应目标的，即确定哪些检测信息是有效的。

当前帧的检测信息，在数据关联问题中，即为量测，通过目标检测，则可得到目标的量测，用观测向量表示：z_k＝[p_x，p_y]^T；同时上一帧目标的跟踪结果，提供了目标的运动状态向量：x_k＝[p_x，v_x，p_y，v_y]^T。本步骤中根据上一帧目标的状态估计得到跟踪门，用于判定当前帧的有效量测，并用确认矩阵表示当前帧的有效量测和多目标跟踪门之间的关系。

设跟踪门限为t_g，如目标跟踪门d_k²(滤波残差向量)的范数小于跟踪门限，则判定当前帧的某量测有效：

$d_k^2\le t_g---\left(4\right)$

$d_k^2={\tilde{z}}_{k|k-1}^T{S}_k^{-1}{\tilde{z}}_{k|k-1}---\left(5\right)$

式中S_k分别是Kalman滤波器【即卡尔曼滤波器】得到的滤波残差向量、残差协方差矩阵。S_k是Kalman滤波器常用的量。

具体如下：滤波残差向量定义为当前帧对应的量测z_k与预测观测量之差：

${\tilde{z}}_{k|k-1}=z_k-{\hat{z}}_{k|k-1}---\left(6\right)$

$=z_k-C{\hat{x}}_{k|k-1}$

C为式(3)定义的观测矩阵，是Kalman滤波器中一步预测状态向量：${\hat{x}}_{k|k-1}={\mathrm{Ax}}_{k-1|k-1},$是状态向量x在k-1时刻的估计值，与k时刻的状态向量估计值的计算相同，见第5点关于状态估计的计算。

残差协方差矩阵S_k定义为：$S_k={\mathrm{CP}}_{k|k-1}{C}^T---\left(7\right)$P_k|k-1为Kalman滤波器中的一步预测协方差矩阵(即根据前一时刻(k-1时刻)的协方差矩阵预测得到当前时刻(k时刻)的预测协方差矩阵)：

P_k|k-1＝AP_k-1|k-1A^T+Q_k-1               (8)

初始值$P_{0|0}=\left(\begin{array}{cccc}0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right)$

P_k-1|k-1是k-1时刻的状态估计协方差阵，与k时刻的状态估计协方差阵P_k|k的计算相同，见式(23)、(24)。

则确认矩阵Ω表示当前帧有效量测与目标跟踪门间的关系，定义为：

$\Omega ={\left\{{\omega }_j^t\right\}}_{j=\mathrm{1,2},...,m}^{t=\mathrm{0,1,2},...N}---\left(9\right)$

式中，t表示目标，j表示量测，ω_j^t表示目标t与量测j的关系。确认矩阵Ω的行表示当前帧量测，列表示跟踪目标，其中${\omega }_j^t=1$表示当前帧量测j落入跟踪目标t的跟踪门内；${\omega }_j^t=0$表示量测j没有落入目标t的跟踪门内，并令t＝0表示虚警(量测对应的是干扰而不是目标，即误把干扰当目标)，对应Ω第一列的所有元素为1，用于表示当前帧任一量测都有可能源于当前检测错误产生的虚警，

如图2，当前帧两个有效量测z₁、z₂与两个目标t₁、t₂的跟踪门及虚警t₀形成确认矩阵Ω：

下文中的量测均指有效量测。

2、产生所有联合(可行)事件及其参数：分析得到所有量测与各目标的所有可能的匹配情况，建立一一对应关系。这些参数用于后续步骤计算关联概率。

多目标的数据关联的难点在于：量测可能在不同目标的跟踪门内，即该量测可能源自多个目标。为了分析量测与目标匹配的各种可能，可通过步骤1产生的确认矩阵，分析得到所有量测与各目标的所有可能的匹配情况，用联合事件表示所有量测与各目标匹配的一种可能，并可根据匹配情况得到联合事件的参数。由于当目标、量测数目较多时，计算确认矩阵对应的所有联合事件会产生“组合爆炸”问题，因此本步骤中采用简化的murty算法计算确认矩阵对应的K个最优联合事件。具体内容如下：

由k时刻的确认矩阵Ω分析得到所有联合事件：${\theta }_k={\left\{{\theta }_{k,i}\right\}}_{i=1}^{n_k},$n_k表示集合θ_k中元素的个数，其中第i个联合事件为${\theta }_{k,i}=\underset{j=1}{\overset{m_k}{\cap }}{\theta }_{k,i}^{j,t_j},$表示当前帧m_k个量测与各个目标匹配的一种可能。而表示第j个量测与目标t关联的事件为θ_k^j，t，称为关联事件。【${\theta }_{k,i}=\underset{j=1}{\overset{m_k}{\cap }}{\theta }_{k,i}^{j,t_j}$中∩符号表示“所有”，本公式代表一个集合】

依据产生联合事件的两个基本假设：每个量测有唯一的源；对应一个给定目标，最多有一个量测以其为源，对确认矩阵Ω拆分，得到与联合事件对应的可行矩阵：

$\hat{\Omega }\left({\theta }_{k,i}\right)=\left[{\hat{\omega }}_j^t\left({\theta }_{k,i}\right)\right]---\left(11\right)$

其中描述在第i个联合事件中，量测j是否源于目标t。可见，计算联合事件θ_k，i，只需得到其对应的可行矩阵

本发明提出一种简化的murty算法【穆尔蒂算法】，快速计算确认矩阵产生最优K个的可行矩阵，(murty算法的基本思想是将给定的一个问题及解，划分为几个将解空间对应划分的子问题，通过找出子问题的解来获取原问题的解。因此，对应于获取确认矩阵Ω的最佳K个匹配关系，即是先构造一个包含所有可能匹配的集合，每次找到一个最优匹配，并在匹配集合中删除这个最优匹配，再从剩下的匹配集合中求一个最优匹配，如此循环K次，即得到确认矩阵Ω中最佳K个匹配关系。可见，如何求取每个匹配集合中的最优匹配是问题的关键。)包括如下步骤：

(1)根据确认矩阵，构建目标、量测的赋权二分图，得到待匹配矩阵

赋权二分图的左边节点表示跟踪目标，右边节点表示当前帧量测，t₀表示虚警，连接目标和量测的线段表示一个可能的分配，线段的权值为线段两端目标、量测在确认矩阵中所在行、列的对应元素；

再将确认矩阵转化为代匹配矩阵：如目标、量测数一致，则直接将确认矩阵作为代匹配矩阵；如目标、量测数不一致，则对确认矩阵增加虚目标或虚量测，构成待匹配矩阵。

例如：

如图2对应的确认矩阵(式(10))可用如下的赋权二分图描述：

图中，实线对应权值为1，虚线对应权值为0。采用解决线性分配问题的匈牙利算法可得到赋权二分图中t、z的最优匹配。

由于匈牙利算法要求赋权二分图的左、右节点数相同，因此如跟踪目标数(包括虚警)、量测数不一致时，采用增加虚目标或虚量测的方式，使左、右节点一致，如在图3中，量测少于跟踪目标数，则增加一个虚量测，虚量测与各目标间的权值为0，得到图4所示的赋权二分图：

得到与图4对应的待匹配的矩阵

(2)匈牙利算法求待匹配矩阵对应的赋权二分图中目标、量测的最优匹配，并得到匹配结果矩阵

①基本概念：

二分图表示为G＝(V，E)，其中V表示顶点集，E表示边集。若二分图G的顶点集V可以分为两个不相交的子集X、Y，且每个子集内部顶点间不存在边，而边只存在于两个子集的顶点之间，则称该图为二分图。

a)匹配：设M是二分图_G的边集合中的一个子集，如果M中任意两条边在G中均不邻接，则称M是G的一个匹配。M中的一条边的两个端点叫做在M是配对的。

b)饱和与非饱和：若匹配M的某条边与顶点v关联，则称M饱和顶点v，并且称v是M-饱和的，否则称v是M-不饱和的。

c)交互道：若M是二分图G＝(V，E)的一个匹配。设从图G中的一个顶点到另一个顶点存在一条道路，这条道路是由属于M的边和不属于M的边交替出现组成的，则称这条道路为交互道。

d)可增广道路：若一交互道的两端点为关于M非饱和顶点时，则称这条交互道是可增广道路。若一条边的两端点非饱和，则这条边也是可增广道路。

②匈牙利算法的基本步骤：

根据Hall定理得出的匈牙利算法是求二分图最大匹配的一种算法，主要步骤如下：

a)任给初始匹配。

b)若X已饱和则结束，否则转c)。

c)在X中找一非饱和点x₀，作：V₁←{x₀}(点集V₁初始只取x₀点)；(V₂为空集，表示空集)。

d)若Г(V₁)＝V₂(和V₁邻接的点集为Г(V₁))则因无法匹配而停止；否则任选一点y∈Г(V₁)\V₂【这个公式表示从点集Г(V₁)中去掉V₂中的点】。e)若y已饱和转f)，否则求一条从x₀→y的可增广道路P，M←M⊕E(P)，转b)。(⊕表示对称差：设A、B是两个集合，则定义：$A\oplus B\underset{‾}{\underset{‾}{\Delta }}(A\cup B)\backslash (A\cap B),$运算规则：0⊕0＝1⊕1＝0，1⊕0＝0⊕1＝1；E(P)表示可增广道路P的边集合)f)由于y已饱和，故M中有一条边(y，z)作：V₁←V₁∪{z}，V₂←V₂∪{y}【该公式表示点集运算，点集V₁新增加边(y，z)中的点z；点集V₂新增加边(y，z)中的点y】，转d)。

再将匈牙利算法求得的赋权二分图最优匹配，转化为匹配结果矩阵：赋权二分图中匹配的目标、量测，在匹配结果矩阵中对应元素为1，否则为0。

例如：

如对图4所示的赋权二分图，按以上匈牙利算法的基本步骤求得最优匹配结果：t₁与z₁匹配、t₂与z₂匹配，转化为匹配结果矩阵：

(3)如未在步骤(1)的待匹配矩阵中增加虚目标或虚量测，则转步骤(4)；如增加了，则在匹配结果矩阵中删去，得到确认矩阵的一个最优匹配关系，即获得确认矩阵的一个可行矩阵。

例如：对式(13)删去增加的虚量测所在的行：

$\hat{\Omega }\left({\theta }_{k,1}\right)=\left(\begin{array}{ccc}0& 1& 0\\ 0& 0& 1\end{array}\right)---\left(14\right)$

即得到确认矩阵(式(10))的一个可行矩阵。

(4)将当前步骤(1)中求得的待匹配矩阵减去当前匹配结果矩阵，得到新的需优化的待匹配矩阵

例如，步骤(1)中求得的待匹配矩阵(式(12))减去当前匹配结果矩阵(式(13))：

$\stackrel{·}{\Omega }-{\stackrel{·}{\Omega }}_1=\left(\begin{array}{ccc}1& 1& 0\\ 1& 0& 1\\ 0& 0& 0\end{array}\right)-\left(\begin{array}{ccc}0& 1& 0\\ 0& 0& 1\\ 0& 0& 0\end{array}\right)=\left(\begin{array}{ccc}1& 0& 0\\ 1& 0& 0\\ 0& 0& 0\end{array}\right)---\left(15\right)$

即得到新的待匹配矩阵为：

$\stackrel{··}{\Omega }=\left(\begin{array}{ccc}1& 0& 0\\ 1& 0& 0\\ 0& 0& 0\end{array}\right)---\left(16\right)$

(5)重复执行步骤(2)——(4)K次，即可找到对应于确认矩阵Ω的K个最优可行矩阵，也就得到确认矩阵对应的联合事件。

联合事件的数目K取值越大，计算结果越接近真实值，但是计算复杂度也会增加，实际应用中，K值可根据系统计算能力和需达到的实时性高低共同确定。本发明中，K取动态值，设定为当前帧跟踪目标数和量测数之和。

最后根据求得的可行矩阵，得到对应的联合事件θ_k，i中的参数：假量测数Φ(θ_k，i)、量测关联指示器τ_j(θ_k，i)及目标检测指示器δ_t(θ_k，i)。

量测关联指示器：${\tau }_j\left({\theta }_{k,i}\right)=\left(\begin{array}{c}1,t_j>0\\ 0,t_j=0\end{array}\right),j=\mathrm{1,2},...,m_k,$其中t_j是联合事件θ_k，i对应的可行矩阵中与量测j关联的目标的取值，即τ_j(θ_k，i)表明量测j是否和一个真实目标关联；

目标检测指示器：其中N表示目标数；则δ_t(θ_k，i)表示在联合事件θ_k，i对应的可行矩阵中目标t是否被检测到；

假量测数Φ(θ_k，i)表示联合事件θ_k，i对应的可行矩阵中假量测的数目，根据量测关联指示器得到：$\Phi \left({\theta }_{k,i}\right)=\underset{j=1}{\overset{m_k}{\Sigma }}[1-{\tau }_j\left({\theta }_{k,i}\right)];$

如式(14)对应联合事件θ_k，1的可行矩阵：

可知量测关联的目标取值t_j为0或1，对量测z₁来说，所在行的第2列取值为1，则量测z₁与真实目标t₁关联，即量测关联指示器τ_j(θ_k，1)＝1(量测j＝z₁)；

对目标t₁来说，其所在列对应两个量测z₁、z₂中存在量测z₁，使得目标t₁对应取值1，即表明目标t₁被检测到，则目标检测指示器${\delta }_{t_1}\left({\theta }_{k,1}\right)=\underset{j=1}{\overset{m_k}{\Sigma }}{\hat{\omega }}_j^{t_1}\left({\theta }_{k,1}\right)=1$(量测j＝z₁、z₂，量测数m_k＝2)

假量测数Φ(θ_k，1)表示可行矩阵中不与任何真实目标关联的量测数目，由于量测z₁与真实目标t₁关联、量测z₂与真实目标t₂关联，因此假量测数Φ(θ_k，1)＝0。

3、应用Bayes法则，计算k时刻所有量测的联合事件的条件概率：

$P\left({\theta }_{k,i}\right|Z^k)=\frac{\Phi \left({\theta }_{k,i}\right)!}{c{V}^{\Phi \left({\theta }_{k,i}\right)}}\underset{j=1}{\overset{m_k}{\Pi }}{\left({\Lambda }_{k,j}\right)}^{{\tau }_j\left({\theta }_{k,i}\right)}\underset{t=1}{\overset{n}{\Pi }}{\left(P_D^t\right)}^{{\delta }_t\left({\theta }_{k,i}\right)}{(1-P_D^t)}^{1-{\delta }_t\left({\theta }_{k,i}\right)}---\left(18\right)$

其中，c为归一化常数，V表示跟踪门体积，Λ_k，j表示滤波残差似然函数，P_D^t表示目标t的检测概率。

具体说明如下：

检测概率P_D^t为常数，实验视频中设为：$P_D^t=0.98;$

滤波残差似然函数：

${\Lambda }_{k,j}=\frac{1}{{\left(2\pi \right)}^{w/2}{\left|S_k^{t_j}\right|}^{1/2}}\exp \left[{(z_{k,j}-{\hat{z}}_{k|k-1}^{t_j})}^T{S}_k^{{t_j}^{-1}}(z_{k,j}-{\hat{z}}_{k|k-1}^{t_j})\right]---\left(19\right)$

其中代表目标t_j的预测位置：${\hat{z}}_{k|k-1}^{t_j}=C{\hat{x}}_{k|k-1},$代表相应于目标t_j的残差协方差矩阵。

跟踪门体积V：其中β为常数，实验视频中取β＝9.5。

4、计算量测与目标t的关联概率β_k^j，t，β_k^0，t【β_k^j，t有效量测j和所有目标的关联概率，β_k^0，t为虚量测和所有目标的关联概率】，用于评估量测与目标关联的可能性大小：

${\beta }_k^{j,t}=P\left({\theta }_k^{j,t}\right|Z^k)=\underset{i=1}{\overset{n_k}{\Sigma }}P\left({\theta }_{k,i}\right|Z^k){\hat{\omega }}_j^t\left({\theta }_{k,i}\right),j\in [1,m_k]---\left(20\right)$

${\beta }_k^{0,t}=1-\underset{j=1}{\overset{m_k}{\Sigma }}{\beta }_k^{j,t};$

5、通过Kalman滤波器，得到目标状态估计及协方差阵：根据状态估计获得目标在当前帧xy方向的位置和速度。协方差阵是估计和真实目标状态的偏差，表征跟踪的精确性。计算得到目标状态估计向量获得目标中心点的位置、速度，即完成当前帧目标跟踪。由步骤4中计算的关联概率得到k时刻目标t状态估计：

${\hat{x}}_{k|k}^t=\underset{j=0}{\overset{m_k}{\Sigma }}{\beta }_k^{j,t}{\hat{x}}_{k|k,j}^t=\underset{j=0}{\overset{m_k}{\Sigma }}{\beta }_k^{j,t}(A{\hat{x}}_{k-1|k-1}^t+K_k^t(z_k-\mathrm{CA}{\hat{x}}_{k-1|k-1}^t))---\left(21\right)$

K_k^t为目标t在k时刻的Kalman滤波增益矩阵，$K_k^t=P_{k|k-1}^t{C}^T{[{\mathrm{CP}}_{k|k-1}^t{C}^T+R]}^{-1},$P_k|k-1^t表示目标t的一步预测协方差矩阵，各目标的预测协方差矩阵计算相同，即P_k|k-1^t与前文中的P_k|k-1的计算相同；表示目标t在上一帧的状态估计，初始值${\hat{x}}_{0|0}^t=\left(\begin{array}{cccc}{p}_{x_0}& 0& p_{y_0}& 0\end{array}\right),$其中分别表示目标初始位置的x方向、y方向的值。

在视频监控中，多目标存在新出现、消失、遮挡、分离等运动情况，量测和目标的关联情况复杂，因此，必须对式(21)得到的状态估计进行修正，才能得到准确的目标状态估计值，从而实现复杂运动情况下的视频多目标跟踪。本发明提出根据步骤4式(20)计算的关联概率，分析目标的运动状态及与上一帧目标的关联情况，包括五种情况：正常、新出现、消失、遮挡、分离，并根据不同情况修正式(21)的目标状态估计值，完成视频多目标跟踪。具体内容如下：

从式(20)可见，由关联概率构成的关联矩阵β_k(式(22))中前m_k行的关联概率β_k^j，t(1≤j≤m_k)反映第j个有效量测与目标t的关联情况，而从式(20)中可知，关联矩阵β_k中的最后一行关联概率β_k^0，t是由前m_k行的关联概率β_k^j，t计算得到，表示没有任何量测源于目标t的概率。因此，选取关联矩阵β_k中前m_k行的关联概率构成新的关联矩阵β_k′(式(22))，β_k′的行表示当前帧有效量测、列表示上一帧的跟踪目标。通过关联矩阵β_k′分析目标的五种运动情况：目标正常、新出现、消失、遮挡、分离，判定在各运动情况下能否通过式(21)得到准确有效的状态估计，如存在较大误差，则需修正状态估计。

${\beta }_k=\left(\begin{array}{cc}{\beta }_k^{1,t_1}& ...\\ ...& ...\\ {\beta }_k^{m_k,t_1}& ...\\ {\beta }_k^{0,t_1}& ...\end{array}\right)$${\beta }_k^{\prime }=\left(\begin{array}{cc}{\beta }_k^{1,t_1}& ...\\ ...& ...\\ {\beta }_k^{m_k,t_1}& ...\end{array}\right)---\left(22\right)$

1)正常：如图2所示，当目标t跟踪门内有且只有一个有效量测z_j，关联概率矩阵量测z_j所在行有唯一大于0的元素t时，则当前量测z_j与元素t所在列对应的上一帧跟踪目标t关联，认为当前帧目标t运动“正常”，目标帧间状态稳定，近似为线性运动，因此可采用式(21)计算目标t的状态估计，无须修正。

2)新出现：

由于是当前帧出现的新目标，没有上一帧的状态估计值则无法通过式(21)计算当前帧的状态估计值。因此需判定哪些量测是新目标，并确定初始跟踪状态，判定新目标包括两个步骤：

①初步判定：

a)如图5a)所示，当前帧量测z_j在所有目标的跟踪门外，认为z_j可能代表新目标；

b)如关联概率矩阵中量测z_j所在行的元素全为0，则当前帧量测z_j与所有上一帧目标无关联，认为z_j可能代表新目标；

②确认为“新目标”：

由于新目标可能是目标消失一段时间重新出现，则该目标不能认为是新目标，而应恢复对该消失目标的重新跟踪。因此，将初步判定为“新出现”的目标与前k_s帧(实验中设k_s＝5)内所有状态为“消失”的目标进行步骤1——5关联，如能满足“正常”运动状态的条件，则该新目标实际为“消失”目标重新出现，则更新此目标状态为“正常”；如该目标不能与任一“消失”目标关联，则确认此目标为“新出现”。

通过以上分析，判定出哪些量测是新目标，新增加为跟踪对象，并用其当前量测z_j(z_j＝[p_x  p_y]^T)作为新出现目标的初始状态中的位置向量，速度向量设为0，即得到新目标的初始状态估计${\hat{x}}_{k|k}^t={\left(\begin{array}{cccc}{p}_x& 0& p_y& 0\end{array}\right)}^T$。【修正后的结果】

3)消失：

当目标在当前帧消失时，即使跟踪门内无有效量测时，但根据式(20)可得${\beta }_k^{0,t}=1,$同时目标存在上一帧的状态估计值因此仍可根据式(21)得到消失目标在当前帧的状态估计，显然与目标已消失的实际运动情况不符。因此须根据量测与目标的关联情况，判定目标是否消失，以终止消失目标的跟踪。

如关联概率矩阵β_k′中上帧跟踪目标t所在列的元素全为0，则表明当前帧无量测与上一帧目标t关联，目标t在当前帧可能消失，如连续三帧判定为消失，则确认该目标已消失，终止该目标的跟踪；

4)遮挡：目标在当前帧发生遮挡，只能检测为一个区域，如图5b)所示，第k-1帧的两个运动目标t₁、t₂在当前帧相互遮挡，检测为一个区域z_j，无法分别得到各个目标的量测值，即仅得到遮挡区域中心点的量测值(在图5b)中用“+”号表示)，该量测值在遮挡目标跟踪门内(如图5b)中目标跟踪门与量测的关系)，由于遮挡区域的中心点量测值与各遮挡目标实际中心点位置存在差异，仅采用该量测与各跟踪目标的关联概率进行状态估计必然存在较大误差。如当量测与某一目标关联概率较大，而与其它目标关联概率非常小时，式(21)只能得到关联概率较大的目标状态估计，其它遮挡目标状态估计近似为0。如表1所示，遮挡目标t₂与量测z_j关联概率远远大于遮挡目标t₁与量测z_j的关联概率，此时目标t₁状态估计近似为0(即其位置、速度近似为0)，与目标t₁的真实位置、速度误差很大，因此需判定目标是否遮挡，并重新修正每个遮挡目标的状态估计。

如关联概率矩阵β_k′中量测z_j所在行有多个大于0的元素(如表1所示)，则量测z_j与这些元素所在列表示的上一帧跟踪目标关联，即多个上一帧目标在当前帧发生遮挡。由于量测z_j表示遮挡区域在当前帧的位置，遮挡目标t₁、t₂状态估计中的位置向量修正为量测z_j的值、速度向量修正为0，即：${\hat{x}}_{k|k}^{t_1}={\left(\begin{array}{cccc}{p}_x& 0& p_y& 0\end{array}\right)}^T,$${\hat{x}}_{k|k}^{t_2}={\left(\begin{array}{cccc}{p}_x& 0& p_y& 0\end{array}\right)}^T.$

表1遮挡情况下量测z_j所在行

  t₁  t₂  z_j  2.7939e-8  0.6231

5)分离【一个目标分成2个或多个】：由于目标检测结果无法做到完全准确，因此可能造成检测到的目标不完整(即目标碎片)，认为目标在当前帧发生分离，如图5c)，第k-1帧目标t在当前帧分离为两个检测区域z₁、z₂。由于采用式(21)估计目标状态的本质正是根据跟踪目标与当前多个量测间的关联情况，得到目标状态估计，而目标分离产生的多个目标碎片恰好提供了该跟踪目标的多个量测(如图5c)中目标跟踪门与量测的关系)，因此根据式(21)能准确有效地估计跟踪目标在当前帧的运动状态，无须修正状态估计值。

如关联概率矩阵β_k′中上帧跟踪目标t所在列有多个大于0的元素，则该列表示的上一帧跟踪目标t与这些元素所在行对应的多个量测关联，表示上一帧跟踪目标t在当前帧分离，如当前帧量测z₁、z₂在目标t的跟踪门内，并且关联矩阵中目标t所在的列为表2：

表2分离目标t所在列

  t  z₁  0.0005  z₂  0.0005

从表2中可见，跟踪目标t在当前帧分离为量测z₁、z₂，并且关联概率给出了量测z₁、z₂与目标t的关联性，据此通过式(21)得到目标t的状态估计。

最后，根据目标t的状态估计，计算状态估计协方差阵P_k|k^t，以进行下一帧跟踪：

$P_{k|k}^t=P_{k|k-1}^t-(1-{\beta }_k^{0,t})K_k^t{S}_k^t{\left(K_k^t\right)}^T+\underset{j=0}{\overset{m_k}{\Sigma }}{\beta }_k^{j,t}[{\hat{x}}_{k|k,j}^t{\left({\hat{x}}_{k|k,j}^t\right)}^T-{\hat{x}}_{k|k}^t{\left({\hat{x}}_{k|k}^t\right)}^T]---\left(23\right)$

式中K_k^t、S_k^t分别为目标t在k时刻的Kalman滤波增益矩阵、残差协方差矩阵。P_k|k-1^t表示目标t的预测协方差阵，和上文的P_k|k-1计算方法一样，初值P_0|0^t和上文的P_0|0取值一样。

当没有任何量测源于目标t时，计算协方差阵为：

$P_{k|k,0}^t=P_{k|k-1}^t---\left(24\right)$

本发明提出的方法在配置为奔4处理器(主频3.00Ghz)、1GB内存的台式计算机、MATLAB 7.1的编程环境下对典型监控视频进行实验分析。监控视频使用的是PETS-ECCV 2004数据库的视频，选取第270帧到第503帧共233帧进行实验分析，图像分辨率为384×288，监控场景中有两人不规则运动，并且人体区域较小、外观特征少。

图6a)所示为在监控视频的场景，图6b)-d)给出了不同方法在监控视频中的跟踪误差，跟踪误差是由目标中心真实坐标值与目标状态估计的位置值的欧式距离得到，其中图6b)、6c)分别是采用传统JPDA算法、本文算法跟踪目标1、2的误差结果，可见本文方法的跟踪误差远远小于传统JPDA算法，最大误差值仅不到4个像素，同时表3给出了跟踪误差的统计量，可见本文算法大大提高了跟踪精度，并且跟踪时间达到了约11fps，能实现实时快速跟踪。说明：由于本发明主要针对传统JPDA算法进行了改进，因此，改为误差与传统JPDA对比即可。

表3监控视频1跟踪误差的统计量及跟踪时间