CN106600624A - 基于粒子群的粒子滤波视频目标跟踪方法 - Google Patents

Abstract

本发明公开了一种基于粒子群的粒子滤波视频目标跟踪方法，主要解决点粒子描述精确性与量测模糊性之间不匹配的问题。其实现步骤是：1.在粒子滤波框架下将普通点粒子拓展为粒子群，通过目标状态预测得到粒子群中心；2.确定每一个粒子群中的粒子，并提取对应的候选目标；3.通过计算候选目标特征与模板之间的距离，对粒子群中心进行更新，并得到粒子群对应权值；4.对粒子群中心进行重采样，进而融合得到目标的状态估计，实现对目标的可靠跟踪。本发明更为高效的实现了对目标状态空间的覆盖，提高了目标跟踪的性能，提高了对大动态范围小目标的搜索和捕捉性能，可用于视频监控、无人机航拍跟踪、机器人导航等系统。

Description

基于粒子群的粒子滤波视频目标跟踪方法

技术领域

本发明属于跟踪监测技术领域，特别涉及一种视频目标跟踪方法，可用于视频监控、无人机航拍跟踪、机器人导航等系统。

背景技术

在视频跟踪研究中，基于粒子滤波的目标跟踪方法由于其在解决非线性问题方面的优异特性，以及相对开放的跟踪框架对各种特征描述方法的兼容能力，在视频跟踪中得到了广泛的应用。然而，由于采样粒子的空间覆盖能力与实际量测存在的模糊性不相匹配，致使普通粒子滤波跟踪方法对于运动动态范围大的小目标的捕获性能变差。因此，如何提高粒子采样效率，利用较少的粒子实现对目标可能存在区域的有效覆盖，同时与量测模糊性相匹配，已成为视频跟踪研究中的难点问题。

目前，典型的目标跟踪方法有：均值漂移，普通粒子滤波等。其中：

均值漂移方法，虽然简单易实现，计算量小，但因具有局部收敛特性容易陷入局部最优解而产生明显的跟踪滞后效应，对小目标和快速移动的大动态范围目标易产生失跟，同时由于需要解析形式推导均值漂移向量，对特征描述方法有一定要求，使得许多特征方法无法应用到其中；

普通粒子滤波方法的跟踪性能较大程度依赖于所采用的粒子数，且将所用粒子抽象为空间体积为零的点。这使得大量粒子映射到量测空间中，只占据若干个与量测尺度相当的区域，存在点粒子描述精确性与量测模糊性之间不匹配的问题，导致计算量的浪费。

发明内容

本发明的目的在于针对上述已有技术中的不足，提出一种基于粒子群的粒子滤波视频目标跟踪方法，已解决点粒子描述精确性与量测模糊性之间不匹配的问题，从而提升粒子采样效率及对大动态范围小目标的搜索和捕捉性能。

实现本发明的技术方案旨在在跟踪过程中，通过将点粒子拓展为粒子群，以保证估计精度与箱粒子覆盖能力之间的平衡，从而更为高效的实现对目标状态空间的覆盖，提高对目标的搜索和捕捉性能，以及跟踪的实时性。

本发明基于粒子群的粒子滤波视频目标跟踪方法，其步骤包括如下：

(1)初始化粒子群：

(1a)读入k-1时刻的图像I_k-1，初始化k-1时刻的粒子群中心集为其中，表示k-1时刻第i个粒子群的中心，i为粒子群序号，取值为1,2,···,N，N表示粒子群总数，k表示时刻，初始时刻为k＝1；

(1b)初始化目标跟踪窗：B_k-1＝(r_k-1,c_k-1)^T，其中r_k-1和c_k-1分别表示k-1时刻目标跟踪窗的长度和宽度值，T表示向量转置运算；

(1c)根据目标初始状态与目标跟踪窗B_k-1，初始化目标的特征协方差矩阵M作为特征模板；

(2)粒子群状态预测：

(2a)读入k时刻的图像I_k，通过对k-1时刻图像中粒子群中心集的传递，得到k时刻图像中的预测粒子群中心集为其中表示k时刻的第i个预测粒子群中心；

(2b)根据k时刻的预测粒子群中心集确定k时刻的粒子群集其中表示k时刻第i个粒子群中的第j个粒子，j为粒子序号，取值为1,2,···,M，M表示每个粒子群中的粒子总数；

(3)粒子群状态更新

(3a)根据k时刻粒子群集和目标跟踪窗B_k-1，确定k时刻的候选目标集其中，表示k时刻第i个粒子群中的第j个粒子所对应的候选目标，它表示以为中心、B_k-1为长宽所界定出的矩形区域；

(3b)提取k时刻的候选目标集对应的候选目标特征集其中表示k时刻第i个粒子群中的第j个粒子所对应候选目标的特征协方差矩阵；

(3c)根据候选目标特征集确定权值集其中表示k时刻第i个粒子群中的第j个粒子所对应的权值；

(3d)根据权值集确定k时刻的更新粒子群中心集并确定粒子群权值集其中表示k时刻第i个更新粒子群中心，表示k时刻第i个粒子群所对应的权值；

(4)重采样：

利用重采样算法，根据粒子群权值集对k时刻的更新粒子群中心集进行重采样，得到k时刻的粒子群中心集其中表示k时刻第i个粒子群的中心；

(5)目标状态估计：

(5a)根据k时刻的粒子群中心集估计k时刻的目标状态X_k；

(5b)根据k时刻的目标状态X_k和目标跟踪窗B_k，确定出k时刻目标T_k，并输出；

(6)判断是否迭代：

检查下一时刻的信息是否到达，若是，令k＝k+1，返回步骤(2)进行迭代，否则，目标跟踪过程结束。

其中步骤(1a)中初始化k-1时刻的粒子群中心集为通过如下步骤确定：

(1a1)根据目标的初始状态X₀，产生k-1时刻第i个粒子群中心其中，服从均值为X₀方差为Ψ₁的高斯分布，X₀为目标的初始状态，Ψ₁为过程噪声方差；

(1a2)用步骤(1a1)所得的N个粒子群中心，组成粒子群中心集

其中步骤(2a)中通过对k-1时刻图像中粒子群中心集的传递，得到k时刻图像中的预测粒子群中心集为通过如下步骤确定：

(2a1)对k-1时刻第i个粒子群中心进行传递，得到k时刻第i个预测粒子群中心

其中，v_k-1为状态噪声，其服从均值为状态噪声协方差为Ψ₁的高斯分布；

(2a2)用步骤(2a1)中所得的N个粒子群中心，组成粒子群中心集

其中步骤(2b)中根据k时刻的预测粒子群中心集确定k时刻的粒子群集通过如下步骤确定：

(2b1)根据k时刻第i个预测粒子群中心确定k时刻第i个粒子群其中，服从均值为方差为Ψ₂的高斯分布，Ψ₂为粒子群过程噪声方差；

(2b2)用步骤(2b1)得到的N个粒子群，组成k时刻的粒子群集

其中步骤(3a)中根据k时刻粒子群集和目标跟踪窗B_k-1，确定k时刻的候选目标集通过如下步骤确定：

(3a1)确定k时刻第i个粒子群中的第j个粒子所对应的候选目标

其中，和分别表示k时刻i个粒子群中的第j个粒子的横坐标和纵坐标,i为粒子群序号，取值为1,2,···,N，j为粒子序号，取值为1,2,···,M，；

(3a2)用步骤(3a1)中所得的N×M候选目标，组成k时刻的候选目标集

其中步骤(3b)中提取k时刻的候选目标集对应的候选目标特征集通过如下步骤确定：

(3b1)计算图像I_k中(x,y)处的像素点所对应的特征向量F(x,y)：

其中，I_k(x,y)为图像I_k像素点(x,y)处的灰度值，和分别表示图像I_k像素点(x,y)处x方向的一阶和二阶梯度，和分别表示图像I_k像素点(x,y)处y方向的一阶和二阶梯度，T表示向量转置，x的取值为1,2···W，y的取值为1,2···H，W表示图像I_k的宽，H表示图像I_k的高；

(3b2)用步骤(3b1)中得到的W×H个特征向量，组成特征图F：

F＝{F(x,y)|x＝1,2,···,W,y＝1,2,···,H}；

(3b3)分别计算特征向量积分图IP中的点IP(x′,y′,a)和特征向量乘积积分图IQ中的点IQ(x′,y′,a,b)：

其中，F(x,y,a)为特征图F中的点F(x,y)处的第a维特征，F(x,y,b)为特征图F中F(x,y)处的第b维特征，a为第一组特征序号，其取值为1,2…d，b为第二组特征序号，其取值为1,2…d，d为特征维数；

(3b4)将步骤(3b3)所得的点IP(x′,y′,a)和IQ(x′,y′,a,b)分别表示成矩阵形式IP_x,y和IQ_x,y：

IP_x,y＝[IP(x,y,1)…IP(x,y,d)]^T，

(3b5)将步骤(3b4)所得的点IP_x,y和IQ_x,y分别组成集合，得到特征向量积分图IP和特征向量乘积积分图IQ：

IP＝{IP_x,y|x＝1,2···W₁,y＝1,2···H₁}，

IQ＝{IQ_x,y|x＝1,2···W₁,y＝1,2···H₁}，

其中，W₁表示IP和IQ的宽，H₁表示积分图IP和IQ的高；

(3b6)计算候选目标集中每一个候选目标所对应的特征协方差矩阵

其中，，(x′,y′)为候选目标对应矩形区域左上角的顶点坐标，(x″,y″)为候选目标对应矩形区域右下角的顶点坐标，n为区域内像素总数n＝(x″-x′)·(y″-y′)，IQ_x″,y″、IQ_x′,y″、IQ_x″,y′、IQ_x′,y′为特征向量乘积积分图IQ中的四个不同点，IP_x″,y″、IP_x′,y″、IP_x″,y′、IP_x′,y′为特征向量积分图IP中的四个不同点；

(3b7)用步骤(3b6)中所得的N×M个特征协方差矩阵，组成候选目标的特征集

其中步骤(3c)中根据候选目标特征集确定权值集通过如下步骤确定：

(3c1)计算k时刻第i个粒子群中的第j个粒子所对应的候选目标特征与特征模板M的之间距离

(3c2)根据距离计算k时刻第i个粒子群中的第j个粒子所对应的权值

其中，R为特征观测噪声方差；

(3c3)用步骤(3c2)中所得的N×M个距离，组成权值集

其中步骤(3d)中根据权值集确定k时刻的更新粒子群中心集并确定粒子群权值集通过如下步骤确定：

(3d1)确定k时刻第i个粒子群中权值最大的粒子的序号β⁽ⁱ⁾：

(3d2)将k时刻第i个粒子群的中心更新为β⁽ⁱ⁾对应的粒子，并将第i个粒子群的权值更新为β⁽ⁱ⁾对应粒子的权值；

(3d3)用步骤(3d2)得到的N个粒子群中心，组成粒子群中心集

(3d4)用步骤(3d2)得到的N个粒子群权值，组成粒子群权值集

本发明具有以下优点：

1)本发明通过将普通的点粒子拓展为粒子群，解决了采样粒子的覆盖能力与量测模糊性不匹配的问题，从而提高了对目标状态空间的覆盖能力；

2)本发明依靠箱粒子较强的空间覆盖能力，提高了采样效率，从而提升了对大动态范围小目标的搜索和捕捉性能。

3)本发明与普通粒子滤波跟踪方法相比较，跟踪精度更高。

附图说明

图1，是本发明的整体流程框图。

图2，为粒子数为50时，普通粒子滤波方法和本发明方法对足球视频序列第2、10、20、30、40、50、60和70帧进行跟踪的结果图。

图3，为粒子数为100时，普通粒子滤波方法和本发明方法对足球视频序列第2、10、20、30、40、50、60和70帧进行跟踪的结果图。

图4，为粒子数为200时，普通粒子滤波方法和本发明方法对足球视频序列第2、10、20、30、40、50、60和70帧进行跟踪的结果图。

图5，为粒子数为50时，普通粒子滤波方法和本发明方法对学生视频序列第2、55、104、159、199和252帧进行跟踪的结果图。

图6，为粒子数为100时，普通粒子滤波方法和本发明方法对学生视频序列第2、55、104、159、199和252帧进行跟踪的结果图。

图7，为粒子数为200时，普通粒子滤波方法和本发明方法对学生视频序列第2、55、104、159、199和252帧进行跟踪的结果图。

图8，为粒子数为50时，普通粒子滤波方法和本发明方法对行人视频序列第2、80、100、130、150和240帧进行跟踪的结果图。

图9，为粒子数为100时，普通粒子滤波方法和本发明方法对行人视频序列第2、80、100、130、150和240帧进行跟踪的结果图。

图10，为粒子数为200时，普通粒子滤波方法和本发明方法对行人视频序列第2、80、100、130、150和240帧进行跟踪的结果图。

图11，为本发明对无人机航拍视频1第1、50、100、150、200和250帧的跟踪结果图。

图12，为本发明对无人机航拍视频1第1、50、100、150、200和250帧的跟踪结果的局部放大图。

图13，为本发明对无人机航拍视频2第1、50、100、150、200和250帧的跟踪结果图。

图14，为本发明对无人机航拍视频2第1、50、100、150、200和250帧的跟踪结果的局部放大图。

具体实施方式

针对上述技术方案，现举一较佳实施例并结合图示进行具体说明。参照图1，本发明的具体实施过程包括以下步骤：

步骤1.初始化粒子群。

1.1)令初始时刻k＝1，读入k-1时刻的图像I_k-1，将普通点粒子拓展为粒子群，并初始化k-1时刻的粒子群中心集为

1.1.1)根据目标的初始状态X₀，产生k-1时刻第i个粒子群中心其中，服从均值为X₀方差为Ψ₁的高斯分布，X₀为目标的初始状态，Ψ₁为过程噪声方差，i为粒子群序号，取值为1,2,···,N，N表示粒子群总数；

1.1.2)用步骤1.1.1)所得的N个粒子群中心，组成粒子群中心集

1.2)初始化目标跟踪窗：B_k-1＝(r_k-1,c_k-1)^T，其中r_k-1和c_k-1分别表示k-1时刻目标跟踪窗的长度和宽度值，T表示向量转置；

1.3)根据目标初始状态X₀与目标跟踪窗B_k-1，初始化目标的特征协方差矩阵M作为特征模板。

步骤2.粒子群状态预测。确定每个粒子群包含的粒子。

2.1)读入k时刻的图像I_k，通过对k-1时刻图像中粒子群中心集的传递，得到k时刻图像中的预测粒子群中心集为

2.1.1)对k-1时刻第i个粒子群中心进行传递，得到k时刻第i个预测粒子群中心

其中，v_k-1为状态噪声，其服从均值为状态噪声协方差为Ψ₁的高斯分布；

2.1.2)用步骤2.1.1)中所得的N个粒子群中心，组成粒子群中心集

2.2)根据k时刻的预测粒子群中心集确定k时刻的粒子群集

2.2.1)根据k时刻第i个预测粒子群中心确定k时刻第i个粒子群其中表示k时刻第i个粒子群中的第j个粒子，j为粒子序号，取值为1,2,···,M，M表示每个粒子群中的粒子总数，服从均值为方差为Ψ₂的高斯分布，Ψ₂为粒子群过程噪声方差；

2.2.2)用步骤2.2.1)得到的N个粒子群，组成k时刻的粒子群集

步骤3.粒子群状态更新。根据每个粒子群的粒子的候选目标、特征提取、权值计算对粒子群进行位置和权值更新。

3.1)根据k时刻粒子群集和目标跟踪窗B_k-1，确定k时刻的候选目标集

3.1.1)确定k时刻第i个粒子群中的第j个粒子所对应的候选目标

其中，和分别表示k时刻i个粒子群中的第j个粒子的横坐标和纵坐标,i为粒子群序号，取值为1,2,···,N，j为粒子序号，取值为1,2,···,M；

3.1.2)用步骤3.1.1)中所得的N×M候选目标，组成k时刻的候选目标集

3.2)提取k时刻的候选目标集对应的候选目标特征集

3.2,1)计算图像I_k中(x,y)处的像素点所对应的特征向量F(x,y)：

其中，I_k(x,y)为图像I_k像素点(x,y)处的灰度值，和分别表示图像I_k像素点(x,y)处x方向的一阶和二阶梯度，和分别表示图像I_k像素点(x,y)处y方向的一阶和二阶梯度，T表示向量转置，x的取值为1,2···W，y的取值为1,2···H，W表示图像I_k的宽，H表示图像I_k的高；

3.2.2)用步骤3.2.1)中得到的W×H个特征向量，组成特征图F：

F＝{F(x,y)|x＝1,2,···,W,y＝1,2,···,H}； 8)

3.2.3)分别计算特征向量积分图IP中的点IP(x′,y′,a)和特征向量乘积积分图IQ中的点IQ(x′,y′,a,b)：

其中，F(x,y,a)为特征图F中的点F(x,y)处的第a维特征，F(x,y,b)为特征图F中F(x,y)处的第b维特征，a为第一组特征序号，其取值为1,2…d，b为第二组特征序号，其取值为1,2…d，d为特征维数；

3.2.4)将步骤3.2.3)所得的点IP(x′,y′,a)和IQ(x′,y′,a,b)分别表示成矩阵形式IP_x,y和IQ_x,y：

IP_x,y＝[IP(x,y,1)…IP(x,y,d)]^T， 11)

3.2.5)将步骤3.2.4)所得的点IP_x,y和IQ_x,y分别组成集合，得到特征向量积分图IP和特征向量乘积积分图IQ：

IP＝{IP_x,y|x＝1,2···W₁,y＝1,2···H₁}， 13)

IQ＝{IQ_x,y|x＝1,2···W₁,y＝1,2···H₁}， 14)

其中，W₁表示IP和IQ的宽，H₁表示积分图IP和IQ的高。

3.2.6)计算候选目标集中每一个候选目标所对应的特征协方差矩阵

其中，，(x′,y′)为候选目标对应矩形区域左上角的顶点坐标，(x″,y″)为候选目标对应矩形区域右下角的顶点坐标，n为区域内像素总数n＝(x″-x′)·(y″-y′)，IQ_x″,y″、IQ_x′,y″、IQ_x″,y′、IQ_x′,y′为特征向量乘积积分图IQ中的四个不同点，IP_x″,y″、IP_x′,y″、IP_x″,y′、IP_x′,y′为特征向量积分图IP中的四个不同点；

3.2.7)用步骤3.2.6)中所得的N×M个特征协方差矩阵，组成候选目标的特征集

3.3)据候选目标特征集确定权值集

3.3.1)计算k时刻第i个粒子群中的第j个粒子所对应的候选目标特征与特征模板M的之间距离

3.3.2)根据距离计算k时刻第i个粒子群中的第j个粒子所对应的权值

其中，R为特征观测噪声方差；

3.3.3)用步骤3.3.2)中所得的N×M个距离，组成权值集

3.4)根据权值集确定k时刻的更新粒子群中心集并确定粒子群权值集

3.4.1)确定k时刻第i个粒子群中权值最大的粒子的序号β⁽ⁱ⁾：

3.4.2)将k时刻第i个粒子群的中心更新为β⁽ⁱ⁾对应的粒子，并将第i个粒子群的权值更新为β⁽ⁱ⁾对应粒子的权值；

3.4.3)用步骤3.4.2)得到的N个粒子群中心，组成粒子群中心集

3.4.3)用步骤3.4.2)得到的N个粒子群权值，组成粒子群权值集

步骤4.粒子群重采样。

现有的重采样方法有多项式重采样、系统重采样、残差重采样等，本实例使用但不限于现有方法中系统重采样方法对箱粒子集进行重采样，再对重采样后的箱粒子集进行放大，得到最后的箱粒子集，其过程如下：

利用重采样算法，根据粒子群权值集对k时刻的更新粒子群中心集进行重采样，得到k时刻的粒子群中心集其中表示k时刻第i个粒子群的中心。

步骤5.估计目标状态。

5.1)根据k时刻的粒子群中心集估计k时刻的目标状态X_k：

其中，N表示粒子总数；

5.2)根据k时刻的目标中心位置和目标跟踪窗，确定k时刻的目标T_k：

其中，(x,y)为目标范围内的任意像素点，x和y分别为该像素点的横坐标和纵坐标，x_k和y_k分别表示k时刻目标状态X_k的横坐标和纵坐标，r_k和c_k分别表示k时刻目标跟踪窗的长度和宽度值；

5.3)输出步骤5.2)所得的目标T_k。

步骤6判断是否迭代。

检查下一时刻的信息是否到达，若是，令k＝k+1，返回步骤2进行迭代，否则，目标跟踪过程结束。

本发明的效果可通过以下实验进一步说明：

1.实验条件。

实验环境：计算机采用Intel Core i5-2400CPU 3.1Ghz，4GB内存，软件采用Matlab R2014a仿真实验平台。

2.实验内容与结果。

实验1：用普通粒子滤波方法和本发明方法，对football视频序列进行跟踪，结果如图2至图4所示，图中白色的框表示本发明的跟踪结果，灰色的框表示普通粒子滤波方法的跟踪结果，其中：

图2，为粒子数为50时，普通粒子滤波方法和本发明方法对足球视频序列第2、10、20、30、40、50、60和70帧进行跟踪的结果图；

图3，为粒子数为100时，普通粒子滤波方法和本发明方法对足球视频序列第2、10、20、30、40、50、60和70帧进行跟踪的结果图；

图4，为粒子数为200时，普通粒子滤波方法和本发明方法对足球视频序列第2、10、20、30、40、50、60和70帧进行跟踪的结果图；

从图2至图4可以看出，本发明对足球视频序列的跟踪结果与实际的真实位置几乎重合，普通粒子滤波方法的跟踪却与真实位置有一定的跟踪差距，因此，本发明的跟踪精度要比普通粒子滤波方法的跟踪精度高。

实验2：用普通粒子滤波方法和本发明方法，对学生视频序列进行跟踪，结果如图5至图7所示，白色的框表示本发明的跟踪结果，灰色的框表示普通粒子滤波方法的跟踪结果，其中：

图5，为粒子数为50时，普通粒子滤波方法和本发明方法对学生视频序列第2、55、104、159、199和252帧进行跟踪的结果图；

图6，为粒子数为100时，普通粒子滤波方法和本发明方法对学生视频序列第2、55、104、159、199和252帧进行跟踪的结果图；

图7，为粒子数为200时，普通粒子滤波方法和本发明方法对学生视频序列第2、55、104、159、199和252帧进行跟踪的结果图；

从图5至图7可以看出，本发明对学生视频序列的跟踪结果与实际的真实位置几乎重合，普通粒子滤波方法的跟踪却与真实位置有一定的跟踪差距，因此，本发明的跟踪精度要比普通粒子滤波方法的跟踪精度高。

实验3：用普通粒子滤波方法和本发明方法，对行人视频序列进行跟踪，结果如图8至图10所示，白色的框表示本发明的跟踪结果，灰色的框表示普通粒子滤波方法的跟踪结果，其中：

图8，为粒子数为50时，普通粒子滤波方法和本发明方法对行人视频序列第2、80、100、130、150和240帧进行跟踪的结果图；

图9，为粒子数为100时，普通粒子滤波方法和本发明方法对行人视频序列第2、80、100、130、150和240帧进行跟踪的结果图；

图10，为粒子数为200时，普通粒子滤波方法和本发明方法对行人视频序列第2、80、100、130、150和240帧进行跟踪的结果图；

从图8至图10可以看出，在目标被车辆遮挡时，本发明可以对目标进行准确跟踪，但是普通粒子滤波方法无法在遮挡时对目标进行跟踪，出现了失跟现象。

实验4：用本发明方法，在粒子数为100时，对无人机航拍视频序列1进行跟踪，结果如图11至图12所示，其中：

图11，为本发明对无人机航拍视频1第1、50、100、150、200和250帧的跟踪结果图；

图12，为本发明对无人机航拍视频1第1、50、100、150、200和250帧的跟踪结果的局部放大图；

从图11至图12可以看出，本发明方法可以对无人机航拍视频1中的大动态范围小目标进行准确跟踪。

实验5：用本发明方法，在粒子数为100时，对无人机航拍视频序列2进行跟踪，结果如图13至图14所示，其中：

图13，为本发明对无人机航拍视频2第1、50、100、150、200和250帧的跟踪结果图；

图14，为本发明对无人机航拍视频2第1、50、100、150、200和250帧的跟踪结果的局部放大图；

从图13至图14可以看出，本发明方法可以对无人机航拍视频2中的大动态范围小目标进行准确跟踪。

实验六

对图2至图4中的场景分别进行100次蒙特卡罗实验，统计平均跟踪误差Err和平均每帧运行时间RT，结果如表1所示。

表1跟踪误差和运行时间统计表

对图5至图7中的场景分别进行100次蒙特卡罗实验，统计平均跟踪误差Err和平均每帧运行时间RT，结果如表2所示。

表2跟踪误差和运行时间统计表

由表1、表2中的统计数据可以看出：在运行时间相当的情况下，图2至图4中的场景，本发明比普通粒子滤波在粒子数50、100、200时跟踪误差分别降低了30.6％、34.4％、35.2％；图5至图7中的场景，本发明比普通粒子滤波在粒子数50、100、200时跟踪误差分别降低了31.9％、32.7％、33.5％。

综上可以得出，在用时相当的情况下，本发明比普通粒子滤波算法的跟踪精度要高，并且可以实现对大动态范围小目标的准确跟踪。

Claims (8)

1.一种基于粒子群的粒子滤波视频目标跟踪方法，包括：

(1)初始化粒子群：

(1a)读入k-1时刻的图像I_k-1，初始化k-1时刻的粒子群中心集为其中，表示k-1时刻第i个粒子群的中心，i为粒子群序号，取值为1,2,…,N，N表示粒子群总数，k表示时刻，初始时刻为k＝1；

(1b)初始化目标跟踪窗：B_k-1＝(r_k-1,c_k-1)^T，其中r_k-1和c_k-1分别表示k-1时刻目标跟踪窗的长度和宽度值，T表示向量转置运算；

(1c)根据目标初始状态与目标跟踪窗B_k-1，初始化目标的特征协方差矩阵M作为特征模板；

(2)粒子群状态预测：

(2a)读入k时刻的图像I_k，通过对k-1时刻图像中粒子群中心集的传递，得到k时刻图像中的预测粒子群中心集为其中表示k时刻的第i个预测粒子群中心；

(2b)根据k时刻的预测粒子群中心集确定k时刻的粒子群集其中表示k时刻第i个粒子群中的第j个粒子，j为粒子序号，取值为1,2,…,M，M表示每个粒子群中的粒子总数；

(3)粒子群状态更新

(3a)根据k时刻粒子群集和目标跟踪窗B_k-1，确定k时刻的候选目标集其中，表示k时刻第i个粒子群中的第j个粒子所对应的候选目标，它表示以为中心、B_k-1为长宽所界定出的矩形区域；

(3b)提取k时刻的候选目标集对应的候选目标特征集其中表示k时刻第i个粒子群中的第j个粒子所对应候选目标的特征协方差矩阵；

(3c)根据候选目标特征集确定权值集其中表示k时刻第i个粒子群中的第j个粒子所对应的权值；

(3d)根据权值集确定k时刻的更新粒子群中心集并确定粒子群权值集其中表示k时刻第i个更新粒子群中心，表示k时刻第i个粒子群所对应的权值；

(4)重采样：

利用重采样算法，根据粒子群权值集对k时刻的更新粒子群中心集进行重采样，得到k时刻的粒子群中心集其中表示k时刻第i个粒子群的中心；

(5)目标状态估计：

(5a)根据k时刻的粒子群中心集估计k时刻的目标状态X_k；

(5b)根据k时刻的目标状态X_k和目标跟踪窗B_k，确定出k时刻目标T_k，并输出；

(6)判断是否迭代：

检查下一时刻的信息是否到达，若是，令k＝k+1，返回步骤(2)进行迭代，否则，目标跟踪过程结束。

2.根据权利要求1所述的方法，其中步骤(1a)中初始化k-1时刻的粒子群中心集为通过如下步骤确定：

(1a1)根据目标的初始状态X₀，产生k-1时刻第i个粒子群中心其中，服从均值为X₀方差为Ψ₁的高斯分布，X₀为目标的初始状态，Ψ₁为过程噪声方差；

(1a2)用步骤(1a1)所得的N个粒子群中心，组成粒子群中心集

${\left\{x_{k-1}^i\right\}}_{i=1}^N=\{x_{k-1}^1,x_{k-1}^2,...,x_{k-1}^i,...,x_{k-1}^N\}.$

3.根据权利要求1所述的方法，其中步骤(2a)中通过对k-1时刻图像中粒子群中心集的传递，得到k时刻图像中的预测粒子群中心集为通过如下步骤确定：

(2a1)对k-1时刻第i个粒子群中心进行传递，得到k时刻第i个预测粒子群中心

$y_k^{\left(i\right)}=x_{k-1}^{\left(i\right)}+v_{k-1},$

其中，v_k-1为状态噪声，其服从均值为状态噪声协方差为Ψ₁的高斯分布；

(2a2)用步骤(2a1)中所得的N个粒子群中心，组成粒子群中心集

${\left\{y_k^{\left(i\right)}\right\}}_{i=1}^N=\{y_k^{\left(1\right)},y_k^{\left(2\right)},...,y_k^{\left(i\right)},...,y_k^{\left(N\right)}\}.$

4.根据权利要求1所述的方法，其中步骤(2b)中根据k时刻的预测粒子群中心集确定k时刻的粒子群集通过如下步骤确定：

(2b1)根据k时刻第i个预测粒子群中心确定k时刻第i个粒子群其中，服从均值为方差为Ψ₂的高斯分布，Ψ₂为粒子群过程噪声方差；

(2b2)用步骤(2b1)得到的N个粒子群，组成k时刻的粒子群集

${\left\{{\left\{s_k^{(i,j)}\right\}}_{j=1}^M\right\}}_{i=1}^N=\{{\left\{s_k^{(1,j)}\right\}}_{j=1}^M,{\left\{s_k^{(2,j)}\right\}}_{j=1}^M,...,{\[s_k^{(i,j)}\]}_{j=1}^M,...,{\left\{s_k^{(N,j)}\right\}}_{j=1}^M\}.$

5.根据权利要求1所述的方法，其中步骤(3a)中根据k时刻粒子群集和目标跟踪窗B_k-1，确定k时刻的候选目标集通过如下步骤确定：

(3a1)确定k时刻第i个粒子群中的第j个粒子所对应的候选目标

$P_k^{(i,v)}=\left\{(x,y)\right||x-{\mathrm{\υ}}_k^{(i,j)}|\≤\frac{r_{k-1}}{2},|y-o_k^{(i,j)}|\≤\frac{c_{k-1}}{2}\},$

其中，和分别表示k时刻i个粒子群中的第j个粒子的横坐标和纵坐标,i为粒子群序号，取值为1,2,…,N，j为粒子序号，取值为1,2,…,M，；

(3a2)用步骤(3a1)中所得的N×M候选目标，组成k时刻的候选目标集

${\left\{{\left\{p_k^{(i,j)}\right\}}_{j=1}^M\right\}}_{i=1}^N=\{{\left\{p_k^{(1,j)}\right\}}_{j=1}^M,{\left\{p_k^{(2,j)}\right\}}_{j=1}^M,...,{\left\{p_k^{(i,j)}\right\}}_{j=1}^M,...,{\left\{p_k^{(N,j)}\right\}}_{j=1}^M\}.$

6.根据权利要求1所述的方法，其中步骤(3b)中提取k时刻的候选目标集对应的候选目标特征集通过如下步骤确定：

(3b1)计算图像I_k中(x,y)处的像素点所对应的特征向量F(x,y)：

$F(x,y)={\[I_k(x,y),\frac{\∂I_k(x,y)}{\∂x},\frac{\∂I_k(x,y)}{\∂y},\frac{{\∂}^2{I}_k(x,y)}{\∂x^2},\frac{{\∂}^2{I}_k(x,y)}{\∂y^2}\]}^T,$

其中，I_k(x,y)为图像I_k像素点(x,y)处的灰度值，和分别表示图像I_k像素点(x,y)处x方向的一阶和二阶梯度，和分别表示图像I_k像素点(x,y)处y方向的一阶和二阶梯度，T表示向量转置，x的取值为1,2…W，y的取值为1,2…H，W表示图像I_k的宽，H表示图像I_k的高；

(3b2)用步骤(3b1)中得到的W×H个特征向量，组成特征图F：

F＝{F(x,y)|x＝1,2,…,W,y＝1,2,…，H}；

(3b3)分别计算特征向量积分图IP中的点IP(x′,y′,a)和特征向量乘积积分图IQ中的点IQ(x′,y′,a,b)：

$IP(x^{\′},y^{\′},a)=\underset{1\≤x\≤x^{\′},1\≤y\≤y^{\′}}{\Σ}F(x,y,a),$

$IQ(x^{\′},y^{\′},a,b)=\underset{1\≤x\≤x^{\′},1\≤y\≤y^{\′}}{\Σ}F(x,y,a)F(x,y,b),$

其中，F(x,y,a)为特征图F中的点F(x,y)处的第a维特征，F(x,y,b)为特征图F中F(x,y)处的第b维特征，a为第一组特征序号，其取值为1,2…d，b为第二组特征序号，其取值为1,2…d，d为特征维数；

(3b4)将步骤(3b3)所得的点IP(x′,y′,a)和IQ(x′,y′,a,b)分别表示成矩阵形式IP_x,y和IQ_x,y：

IP_x,y＝[IP(x,y,1)…IP(x,y,d)]^T，

(3b5)将步骤(3b4)所得的点IP_x,y和IQ_x,y分别组成集合，得到特征向量积分图IP和特征向量乘积积分图IQ：

IP＝{IP_x,y|x＝1,2…W₁,y＝1,2…H₁}，

IQ＝{IQ_x,y|x＝1,2…W₁,y＝1,2…H₁}，

其中，W₁表示IP和IQ的宽，H₁表示积分图IP和IQ的高；

(3b6)计算候选目标集中每一个候选目标所对应的特征协方差矩阵

$\begin{array}{c}{C}_k^{(i,j)}=\frac{1}{n-1}\[{\mathrm{IQ}}_{x^{\′\′},y^{\′\′}}-{\mathrm{IQ}}_{x^{\′},y^{\′\′}}-{\mathrm{IQ}}_{x^{\′\′},y^{\′}}+{\mathrm{IQ}}_{x^{\′},y^{\′}}\\ -\frac{1}{n}({\mathrm{IP}}_{x^{\′\′},y^{\′\′}}-{\mathrm{IP}}_{x^{\′},y^{\′\′}}-{\mathrm{IP}}_{x^{\′\′},y^{\′}}+{\mathrm{IP}}_{x^{\′},y^{\′}}){({\mathrm{IP}}_{x^{\′\′},y^{\′\′}}-{\mathrm{IP}}_{x^{\′},y^{\′\′}}-{\mathrm{IP}}_{x^{\′\′},y^{\′}}+{\mathrm{IP}}_{x^{\′},y^{\′}})}^T\]\end{array},$

其中，，(x′,y′)为候选目标对应矩形区域左上角的顶点坐标，(x″,y″)为候选目标对应矩形区域右下角的顶点坐标，n为区域内像素总数n＝(x″-x′)·(y″-y′)，IQ_x″,y″、IQ_x′,y″、IQ_x″,y′、IQ_x′,y′为特征向量乘积积分图IQ中的四个不同点，IP_x″,y″、IP_x′,y″、IP_x″,y′、IP_x′,y′为特征向量积分图IP中的四个不同点；

(3b7)用步骤(3b6)中所得的N×M个特征协方差矩阵，组成候选目标的特征集

${\left\{{\left\{C_k^{(i,j)}\right\}}_{j=1}^M\right\}}_{i=1}^N=\{{\left\{C_k^{(1,j)}\right\}}_{j=1}^M,{\left\{C_k^{(2,j)}\right\}}_{j=1}^M,...,{\left\{C_k^{(i,j)}\right\}}_{j=1}^M,...,{\left\{C_k^{(N,j)}\right\}}_{j=1}^M\}.$

7.根据权利要求1所述的方法，其中步骤(3c)中根据候选目标特征集确定权值集通过如下步骤确定：

(3c1)计算k时刻第i个粒子群中的第j个粒子所对应的候选目标特征与特征模板M的之间距离

${\mathrm{\ρ}}_k^{(i,j)}=\left|\right|log\left(C_k^{(i,j)}\right)-log\left(M\right)\left|\right|,i=1,2,...,N,j=1,2,...,M;$

(3c2)根据距离计算k时刻第i个粒子群中的第j个粒子所对应的权值

$w_k^{(i,j)}=\frac{1}{\sqrt{2\mathrm{\π}R}}\exp \left\{-\frac{{\left({\mathrm{\ρ}}_k^{(i,j)}\right)}^2}{2R}\right\},$

其中，R为特征观测噪声方差；

(3c3)用步骤(3c2)中所得的N×M个距离，组成权值集

${\left\{{\left\{w_k^{(i,j)}\right\}}_{j=1}^M\right\}}_{i=1}^N=\{{\left\{w_k^{(1,j)}\right\}}_{j=1}^M,{\left\{w_k^{(2,j)}\right\}}_{j=1}^M,...,{\left\{w_k^{(3,j)}\right\}}_{j=1}^M,...,{\left\{w_k^{(N,j)}\right\}}_{j=1}^M\}.$

8.根据权利要求1所述的方法，其中步骤(3d)中根据权值集确定k时刻的更新粒子群中心集并确定粒子群权值集通过如下步骤确定：

(3d1)确定k时刻第i个粒子群中权值最大的粒子的序号β⁽ⁱ⁾：

${\mathrm{\β}}^{\left(i\right)}=\underset{j}{\mathrm{argmax}}\left(w_k^{(i,j)}\right),$

(3d2)将k时刻第i个粒子群的中心更新为β⁽ⁱ⁾对应的粒子，并将第i个粒子群的权值更新为β⁽ⁱ⁾对应粒子的权值；

(3d3)用步骤(3d2)得到的N个粒子群中心，组成粒子群中心集

${\left\{g_k^{\left(i\right)}\right\}}_{i=1}^N=\{g_k^{\left(1\right)},g_k^{\left(2\right)},...,g_k^{\left(i\right)},...,g_k^{\left(N\right)}\},$

(3d4)用步骤(3d2)得到的N个粒子群权值，组成粒子群权值集