基于优化粒子的粒子滤波目标跟踪方法与流程

本发明属于目标跟踪技术领域，特别涉及一种视频目标跟踪方法，可用于智能视频监控、机器人导航、道路交通管制系统。

背景技术：

在视频跟踪研究中，基于粒子滤波的目标跟踪方法由于其在解决非线性问题方面的优异特性，以及相对开放的跟踪框架对各种特征描述方法的兼容能力，在视频跟踪中得到了广泛的应用。然而，粒子滤波方法的重采样过程会导致粒子多样性差，粒子匮乏问题严重，这将使得粒子利用率低，并且对跟踪算法的跟踪精度与计算效率造成很大的影响。因此，如何解决粒子滤波算法中的粒子匮乏问题，提高粒子利用率，已成为视频跟踪研究中的重点问题。

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

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

普通粒子滤波方法，其跟踪性能对粒子数的依赖性较大，并且重采样过程中舍弃了大量权值较小的粒子，复制了权值较大的粒子，使所有粒子都趋向于同一个粒子，导致粒子的多样性丧失，粒子匮乏问题严重，在目标状态估计时只利用了那些权值较大粒子，而权值较小的粒子都未被利用，这将导致对目标的状态进行估计有很大偏差，甚至致使跟踪失败，同时大量的运算时间都浪费在了未被利用的权值较小的粒子上，降低了计算效率。

技术实现要素：

本发明的目的在于针对上述已有技术中的不足，提出一种基于优化粒子的粒子滤波目标跟踪方法，以解决粒子滤波方法在目标跟踪中出现的粒子匮乏以及粒子利用率低的问题，降低跟踪性能对粒子数的依赖程度，减小目标状态估计的偏差，提高目标跟踪的精度。

实现本发明的技术关键是：利用粒子群优化方法对粒子滤波跟踪方法中的粒子进行优化，即在粒子滤波框架中，通过粒子群优化方法优化权值较小的粒子在空间中的位置，使得所有粒 子的分布都集中于高似然区域，减少重采样过程中对权值较小粒子的舍弃，从而增加粒子的多样性，解决粒子匮乏问题，提高粒子利用率。其实现步骤包括如下：

(1)初始化：

(1a)读入k-1时刻的图像I_k-1，根据目标的初始状态X₀，初始化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)初始化目标的灰度直方图M作为特征模板；

(2)目标状态预测：

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

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

(3)提取候选目标集中候选目标对应的特征集其中V_i表示第i个候选目标 对应的特征；

(4)求取候选目标特征集与特征模板M之间的距离集并根据距离集计算候选目标的权值集其中d_i表示第i个候选目标特征与特征模板M之间的距离，ω_i表示第i个候选目标对应的权值；

(5)粒子优化：

(5a)初始化：

(5a1)初始化个体最优值集群体最优值群体最优值对应的权值为ω_max，其中为k时刻预测粒子集中权值最大的粒子，ω_max为权值最大粒子对应的权值；

(5a2)初始化k时刻预测粒子集对应的速度集k时刻预测粒子集对应的权值集为其中v_i表示k时刻预测粒子集中第i个预测粒子对应的速度，其取值是均值为0，方差为1的正态分布的随机数，表示k时刻预测粒子集中第i个预测粒子对应的权值；

(5a3)初始化迭代次数T＝1，设置最大迭代次数T_max；

(5b)判断是否满足ω_min>＝α·ω_max或T>T_max：

若是，令k时刻优化粒子集该优化粒子集对应的权值集直接转到步骤(6)；否则，执行步骤(5c)开始第T次迭代；

其中α是优化系数α＝0.6，ω_max是k时刻预测粒子集对应的权值集中最大的权值，ω_min是k时刻预测粒子集对应的权值集中最小的权值；

(5c)对k时刻预测粒子集进行优化，得到k时刻优化粒子集和该优化粒子集对应的权值集

(5d)令T＝T+1、k时刻预测粒子集k时刻预测粒子集对应的权值集 返回步骤(5b)；

(6)对权值集进行归一化，并根据归一化后的结果对优化粒子集进行重采样，得到k时刻的粒子集其中表示k时刻第i个粒子；

(7)根据k时刻的粒子集估计k时刻的目标状态X_k，并确定k时刻的目标估计范围Target；

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

本发明具有以下优点：

本发明由于利用粒子群优化方法对粒子滤波过程中权值较小的粒子进行优化，减少了重采样过程中对权值较小粒子的舍弃，从而增加了粒子的多样性，解决了粒子匮乏问题，提高了粒子利用率，减少了跟踪性能对粒子数的依赖程度；同时由于增加了粒子的多样性，提高了粒子利用率，因而可以用较少粒子对目标进行准确描述，从而提高了计算效率。

附图说明

图1是本发明的实现流程图；

图2是用本发明及现有方法对boy视频序列进行跟踪的实验结果图；

图3是用本发明及现有方法对football视频序列进行跟踪的实验结果图；

图4是用本发明及现有方法对women视频序列进行跟踪的实验结果图；

图5是用本发明及现有方法对lemming视频序列进行跟踪的实验结果图；

具体实施方式

以下参照附图，对本发明的技术方案和效果进行进一步说明：

参照图1，本发明的具体实施步骤如下：

步骤1.初始化。

1.1)令初始时刻k＝1，根据目标的初始状态X₀，初始化k-1时刻的粒子集

${\left\{S_{k-1}^i\right\}}_{i=1}^N=\{S_{k-1}^1,S_{k-1}^2,...,S_{k-1}^i,...,S_{k-1}^N\};---<1>$

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

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

1.3)初始化目标的灰度直方图M作为特征模板。

步骤2.目标状态预测。

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

${\left\{q_k^i\right\}}_{i=1}^N=\{q_k^1,q_k^2,...,q_k^i,...,q_k^N\},---<2>$

其中表示k时刻第i个预测粒子，其通过对k-1时刻第i个粒子传递得到，传递公式如下：

$q_k^i=S_{k-1}^i+w^i,---<3>$

其中，wⁱ为状态噪声，服从均值为方差为Ψ的高斯分布；

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

${\left\{P_k^i\right\}}_{i=1}^N=\{P_k^1,P_k^2,...,P_k^i,...,P_k^N\},---<4>$

其中表示k时刻第i个预测粒子对应的候选目标，其确定公式如下：

$P_k^i=\left\{(x,y)\right||x-{\mathrm{\&upsi;}}_k^i|\&le;\frac{r_{k-1}}{2},|y-v_k^i|\&le;\frac{c_{k-1}}{2}\},---<5>$

其中，和分别表示k时刻第i个预测粒子的横坐标和纵坐标，r_k-1和c_k-1分别表示k-1时刻目标跟踪窗B_k-1的长度和宽度值。

步骤3.提取候选目标特征集。

常用的目标特征有灰度特征、颜色特征、纹理特征和形状特征。

现有的提取目标特征的方法中针对颜色特征的提取方法有灰度直方图法、颜色直方图法、颜色矩法、颜色相关图法等；针对纹理特征的提取方法有局部二值特征LBP统计直方图法、梯度法、灰度共生矩阵法、自回归纹理模型法、小波变换法等；针对形状特征的提取方法有边界特征法、傅里叶形状描述符法、几何参数法等。

本实例使用但不限于现有方法中的灰度直方图法作为目标特征，并组成候选目标特征集

${\left\{V_i\right\}}_{i=1}^N=\{V_1,V_2,...,V_i,...V_N\},---<6>$

其中V_i表示第i个候选目标对应的灰度直方图，其确定公式如下：

$V_i=C\&CenterDot;\underset{i=1}{\overset{n}{\&Sigma;}}k\left(\right|\left|\frac{y-x_i}{h}\right|{|}^2)\mathrm{\&delta;}\&lsqb;b\left(x_i\right)-u\&rsqb;,---<7>$

其中C为归一化因子u为灰度值u＝0,1,2,…,255，n为候选目标包含的像素数，y为候选目标中心位置，x_i为候选目标范围内任意像素点的位置，b(x_i)表示像素点x_i的的灰度值，δ(·)是克罗内克函数，k(·)是带宽为h的核函数的轮廓函数。

步骤4.计算权值集。

4.1)计算候选目标特征集与特征模板M之间的距离集

${\left\{d_i\right\}}_{i=1}^N=\{d_1,d_2,...,d_i,...d_N\},---<8>$

其中d_i表示第i个候选目标特征与特征模板之间的距离，ρ_i是第i个候选目标特征与特征模板的之间的巴氏系数，

4.2)根据距离集计算候选目标权值集

${\left\{{\mathrm{\&omega;}}_i\right\}}_{i=1}^N=\{{\mathrm{\&omega;}}_1,{\mathrm{\&omega;}}_2,...,{\mathrm{\&omega;}}_i,...{\mathrm{\&omega;}}_N\},---<9>$

其中ω_i表示第i个候选目标的权值，其计算公式为：

${\mathrm{\&omega;}}_i=\frac{1}{\sqrt{2\mathrm{\&pi;}R}}\exp \{-\frac{{\left(d_i\right)}^2}{2R}\},---<10>$

其中，d_i表示第i个候选目标特征与特征模板之间的距离，R为特征观测噪声方差。

步骤5.粒子优化。

5.1.)初始化：

5.1.1)初始化个体最优值集群体最优值群体最优值对应的权值为ω_max，其中为k时刻预测粒子集中权值最大的粒子，ω_max为权值最大粒子对应的权值；

5.1.2)初始化k时刻预测粒子集对应的速度集初始化k时刻预测粒子集 对应的权值集为其中v_i表示k时刻预测粒子集中第i个预测粒子 对应的速度，其取值是均值为0，方差为1的正态分布的随机数，表示k时刻预测粒子集 中第i个预测粒子对应的权值；

5.1.3)初始化迭代次数T＝1，设置最大迭代次数T_max；

5.2)判断是否满足ω_min>＝α·ω_max或T>T_max：

若是，则令k时刻优化粒子集该优化粒子集对应的权值集 直接转到步骤(6)；否则，执行步骤5.3)开始第T次迭代；

其中α是优化系数α＝0.6，ω_max是k时刻预测粒子集对应的权值集中最大的权值，ω_min是k时刻预测粒子集对应的权值集中最小的权值；

5.3)对k时刻预测粒子集进行优化，得到k时刻优化粒子集和该优化粒子集对应的权值集

5.3.1)令粒子序号i＝1；

5.3.2)判断是否满足i<＝N，若是，则执行步骤5.3.3)，否则，执行步骤5.4)；

5.3.3)判断k时刻预测粒子集中第i个粒子对应的权值ω_i是否满足ω_i<α·ω_max：若是，则令k时刻优化粒子集中第i个优化粒子该优化粒子对应权值并执行步骤5.3.4)，否则，直接转到步骤5.3.9)；

5.3.4)更新k时刻预测粒子集中第i个粒子对应的速度v_i：

$v_i=w\&CenterDot;v_i+c_1\&CenterDot;r_1\&CenterDot;(p_i-q_k^i)+c_2\&CenterDot;r_2\&CenterDot;(p_{gd}-q_k^i),---<11>$

其中，w为惯性权值w＝0.9，c₁为与个体最优值有关的学习因子c₁＝2，c₂为与全局最优值有关的学习因子c₂＝2，r₁和r₂为(0,1)之间两个不同的随机数，p_id为个体最优值集中的第i个个体最优值p_id；

5.3.5)用步骤5.3.4)得到的速度v_i，对k时刻预测粒子集中第i个粒子进行优化，得到k时刻优化粒子集中第i个优化粒子

$Q_k^i=q_k^i+v_i;---<12>$

5.3.6)计算k时刻优化粒子集中第i个优化粒子对应的权值ω_i'：

${\mathrm{\&omega;}}_i^{\&prime;}=\frac{1}{\sqrt{2\mathrm{\&pi;}R}}\exp \{-\frac{{\left(d_{new}\right)}^2}{2R}\},---<13>$

其中d_new表示k时刻优化粒子集中第i个优化粒子对应的候选目标特征与特征模板之间的距离；

5.3.7)将第i个优化粒子对应的权值ω_i'与第i个预测粒子对应的权值进行比较，更新个体最优值集中的第i个个体最优值p_i：

如果则更新个体最优值否则p_i值不变；

5.3.8)将第i个优化粒子对应的权值ω_i'与群体最优值对应的权值ω_max进行比较，更新群体最优值p_gd：

如果ω_i'>ω_max，则更新群体最优值其对应的权值ω_max＝ω_i'，否则p_gd值不变，其对应权值ω_max不变；

5.3.9)令i＝i+1，返回步骤5.3.2)；

5.4)令T＝T+1、k时刻预测粒子集k时刻预测粒子集对应的权值集 返回步骤5.2)。

步骤6.重采样。

现有的重采样方法有多项式重采样、系统重采样、残差重采样等，本实例使用但不限于现有方法中的系统重采样方法对优化粒子集进行重采样，其步骤如下：

6.1)对候选目标的权值集进行归一化，得到归一化后的权值集

${\left\{{\mathrm{\&omega;}}_i^{\&prime;\&prime;}\right\}}_{i=1}^N=\{{\mathrm{\&omega;}}_1^{\&prime;\&prime;},{\mathrm{\&omega;}}_2^{\&prime;\&prime;},...,{\mathrm{\&omega;}}_3^{\&prime;\&prime;},...{\mathrm{\&omega;}}_N^{\&prime;\&prime;}\}---<14>$

其中，ω_i”为ω_i'对应的归一化结果，

6.2)根据步骤6.1)得到的归一化后的权值集对k时刻优化粒子集进行重采样，得到k时刻的粒子集其中表示k时刻第i个粒子的状态区间。

步骤7.目标状态更新。

7.1)根据k时刻的粒子集估计k时刻的目标状态X_k：

$X_k=\frac{1}{N}\underset{i=1}{\overset{N}{\&Sigma;}}{S}_k^i,---<15>$

其中，N表示粒子总数，表示k时刻第i个粒子；

7.2)根据k时刻的目标状态X_k和目标跟踪窗B_k，确定k时刻的目标估计范围Target：

$T\arg et=\left\{(x,y)\right||x-x_k|\&le;\frac{r_k}{2},|y-y_k|\&le;\frac{c_k}{2}\},---<16>$

其中，x_k和y_k分别表示k时刻目标状态X_k的横坐标和纵坐标，r_k和c_k分别表示k时刻目标跟踪窗的长度和宽度值。

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

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

1.实验条件。

实验环境：计算机采用Intel(R)Core(TM)i5-2004CPU@3.1Ghz，4GB内存，软件采用Matlab2015b仿真实验平台。

实验参数：过程噪声方差Ψ＝diag([15,15])，特征观测噪声方差R＝0.01.

2.实验方法：本发明方法、现有普通粒子滤波方法；

3.实验内容与结果。

实验1：用所述两种方法，在粒子总数N＝50的情况下，对boy视频序列进行跟踪，结果如图2所示，其中：

图2(a)为用本发明对boy视频序列第1、92、105、125、132、183、192、200帧进行跟踪的结果图；

图2(b)为用普通粒子滤波方法对boy视频序列第1、92、105、125、132、183、192、200帧进行跟踪的结果图；

图2(c)为用本发明对boy视频序列第1、92、105、125、132、183、192、200帧进行跟踪结果的局部放大图；

图2(d)为用普通粒子滤波方法对boy视频序列第1、92、105、125、132、183、192、200帧进行跟踪结果的局部放大图；

从图2中可以看出，在目标快速运动过程中产生旋转和运动模糊时，本发明的跟踪结果明显优于普通粒子滤波方法。

实验2：用所述两种方法，在粒子总数N＝50的情况下，对football视频序列进行跟踪，结果如图3所示，其中：

图3(a)为用本发明对football视频序列第1、10、20、30、40、50、60、70帧进行跟踪的结果图；

图3(b)为用普通粒子滤波方法对football视频序列第1、10、20、30、40、50、60、70帧进行跟踪的结果图；

图3(c)为用本发明对football视频序列第1、10、20、30、40、50、60、70帧进行跟踪结果的局部放大图；

图3(d)为用普通粒子滤波方法对football视频序列第1、10、20、30、40、50、60、70帧进行跟踪结果的局部放大图；

从图3中可以看出，在目标快速运动且遇到相似背景时，本发明的跟踪结果明显优于普通粒子滤波方法。

实验3：用所述两种方法，在粒子总数N＝50的情况下，对women视频序列进行跟踪，结果如图4所示，其中：

图4(a)为用本发明对women视频序列第1、110、130、150、180、200、220、240帧进行跟踪的结果图；

图4(b)为用普通粒子滤波方法对women视频序列第1、20、40、60、80、100、120、125帧进行跟踪的结果图；

图4(c)为用本发明对women视频序列第1、110、130、150、180、200、220、240帧进行跟踪结果的局部放大图；

图4(d)为用普通粒子滤波方法对women视频序列第1、20、40、60、80、100、120、125帧进行跟踪结果的局部放大图；

从图4中可以看出，在目标遇到遮挡情况时，本发明可以实现对目标的准确跟踪，而普通粒子滤波方法在第120帧以后出现失跟，无法继续对目标进行跟踪。

实验4：用所述两种方法，在粒子总数N＝50的情况下，对lemming视频序列进行跟踪，结果如图5所示，其中：

图5(a)为用本发明对lemming视频序列第1、30、50、70、90、110、130、150帧进行跟踪的结果图；

图5(b)为用普通粒子滤波方法对lemming视频序列第1、5、15、20、25、30、32帧进行跟踪的结果图；

图5(c)为用本发明对lemming视频序列第1、30、50、70、90、110、130、150帧进行跟踪结果的局部放大图；

图5(d)为用普通粒子滤波方法对lemming视频序列第1、5、15、20、25、30、32帧进行跟踪结果的局部放大图；

从图5中可以看出，在目标进入复杂背景时，本发明可以实现对目标的准确跟踪，而普通粒子滤波方法在30帧以后出现失跟，无法继续对目标进行跟踪。

对图2的场景分别在粒子总数N＝10，N＝20，N＝50，N＝100，N＝500的情况下进行50次蒙特卡罗实验，统计平均跟踪误差Err和平均每帧运行时间RT，结果如表1所示：

表1图2场景跟踪误差与运行时间统计表

对图3的场景分别在粒子总数N＝10，N＝20，N＝50，N＝100，N＝500的情况下进行50次蒙特卡罗实验，统计平均跟踪误差Err和平均每帧运行时间RT，结果如表2所示：

表2图3场景跟踪误差与运行时间统计表

由表1、表2中的统计数据可以看出：

1)粒子数N＝10时本发明对两组视频的跟踪误差比普通粒子滤波方法分别降低了65.8％、53.4％；粒子数N＝20时本发明的跟踪误差比普通粒子滤波方法分别降低了55.1％、53.9％；粒子数N＝50时本发明的跟踪误差比普通粒子滤波方法分别降低了47.9％、51.9％；粒子数N＝100时本发明的跟踪误差比普通粒子滤波方法分别降低了44.7％、51.8％；粒子数N＝500时本发明的跟踪误差比普通粒子滤波方法分别降低了44.3％、50.5％；

2)在粒子数相同时虽然本发明比普通粒子滤波方法运行时间长，但是跟踪精度比普通粒子滤波方法高很多；

3)在用时相当的情况下，本发明跟踪精度远优于普通粒子滤波方法；

4)本发明跟踪性能对粒子数的依赖性非常小，但是普通粒子滤波方法对粒子数有很大依赖性；

5)随着粒子数的增加，本发明与普通粒子滤波方法的跟踪误差趋于平稳，此时，本发明的跟踪误差远小于普通粒子滤波方法。

统计图2场景的第105帧、图3场景的第20帧、图4场景的第120帧、图5场景的第30帧分别在粒子总数N＝10，N＝20，N＝50，N＝100，N＝500的情况下用所述两种方法进行跟踪时的粒子利用率，结果如表3所示：

表3粒子利用率统计表

由表3的统计数据可以看出：在粒子数N＝10的情况下，本发明对四种场景跟踪时粒子利用率比普通粒子滤波方法分别提高了167％、50％、80％、28.6％，在粒子数N＝20的情况下，本发明对四种场景跟踪时粒子利用率比普通粒子滤波方法分别提高了200％、72.7％、58.3％、24.4％，在粒子数N＝50的情况下，本发明对四种场景跟踪时粒子利用率比普通粒子滤波方法 分别提高了320％、53.3％、20.6％、42.9％，在粒子数N＝100的情况下，本发明对四种场景跟踪时粒子利用率比普通粒子滤波方法分别提高了190％、22％、27.9％、47.4％，在粒子数N＝500的情况下，本发明对四种场景跟踪时粒子利用率比普通粒子滤波方法分别提高了367％、39.1％、28.7％、38.5％。

综上可以得出，本发明粒子利用率高于普通粒子滤波方法，而且跟踪性能对粒子数的依赖性较小，在用时相当的情况下本发明的跟踪精度远优于普通粒子滤波方法。