CN103985141B - 基于hsv颜色协方差特征的目标跟踪方法 - Google Patents

Abstract

本发明公开了一种基于HSV颜色协方差特征的目标跟踪方法，主要解决现有跟踪技术针对彩色图像目标提取的特征存在信息冗余大和特征间独立性差的问题。其实现步骤是：在粒子滤波框架下首先通过目标状态预测得到候选目标，提取候选目标图像的色相、饱和度、明度、拉普拉斯算子响应作为表观特征，融合构建协方差算子；其次通过计算候选目标特征与模板的相似度权值，调整跟踪窗，并更新特征模版；最后根据权值融合候选目标，更新目标状态，实现对目标的有效跟踪。本发明能够有效减少目标特征的信息冗余，增强特征之间的独立性，提高目标特征描述的精确性和鲁棒性，可实现对彩色图像目标的实时精确跟踪。

Description

基于HSV颜色协方差特征的目标跟踪方法

技术领域

本发明属于跟踪监测技术领域，特别涉及一种彩色图像的目标跟踪方法，可用于视频跟踪、目标监控等系统。

背景技术

彩色图像，由于其图像细节丰富，感光器件成本低，在实际中应用范围广，因此针对彩色图像的目标跟踪有着广阔的应用前景和实际意义。丰富的图像细节，既给目标特征描述提供了大量素材，同时也带来了特征甄别选取的难度。如何充分挖掘目标的图像信息，去除其中的冗余，有效地组织和建立多特征描述模型，成为彩色图像目标跟踪研究中的重点问题。

目前，针对彩色目标跟踪的特征描述模型主要有：联合特征直方图，协方差算子等。其中：

联合特征直方图，存在大量的冗余信息，计算复杂度高，对光照变化敏感；

协方差算子，是一个高效精炼的多特征融合描述模型，具有较强的鲁棒性，但其传统定义中所选取的特征主要针对目标在各个方向上纹理的描述，同方向上不同阶次的梯度在某些情况下会存在一定程度的信息冗余，而且缺乏对彩色图像中丰富的颜色信息的挖掘。

传统的彩色图像目标跟踪方法，包括基于联合特征直方图的目标跟踪方法和基于协方差矩阵的目标跟踪方法。其中，基于联合直方图的跟踪方法选取红色r、绿色g和蓝色b构成联合颜色直方图；基于协方差矩阵的跟踪方法中对目标所提取的特征包括灰度值、x方向的一阶梯度、x方向的二阶梯度、y方向的一阶梯度、y方向的二阶梯度。这些传统方法在实际情况下，由于所提取的特征间具有较高的相似性，特征间的相互独立性较差，存在信息冗余，导致目标特征描述的准确性和鲁棒性较低，使得跟踪系统的稳定性较差。

发明内容

本发明的目的在于针对上述已有技术的不足，在粒子滤波框架下，利用协方差算子特征描述模型，融合了色相H、饱和度S、明度V和拉普拉斯算子响应，提出了一种基于HSV颜色协方差特征的目标跟踪方法，以提高目标特征描述的准确性和鲁棒性，实现对彩色图像中运动目标的实时精确跟踪。

实现本发明的技术关键是：在跟踪过程中，通过提取目标图像的色相、饱和度、明度、拉普拉斯算子响应，获取目标图像的表观特征，通过表观特征对候选目标建立协方差矩阵，以减少信息冗余，增强特征独立性，从而提升目标特征描述的精确性和鲁棒性，提高跟踪精度。其实现步骤包括如下：

(1)初始化步骤：

(1a)读入k-1时刻的彩色图I_k-1，根据目标的初始状态，产生k-1时刻的初始粒子集其中，表示k-1时刻第i个粒子的状态估计值，N表示粒子总数，i表示粒子的序号，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时刻的彩色图Ik，通过对k-1时刻彩色图中粒子集的传递，得到k时刻彩色图中的预测粒子集其中为k时刻第i个粒子的状态预测值；

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

(3)候选目标特征提取步骤：

(3a)对于k时刻彩色图I_k，提取其对应的W×H×d维特征图F，其中W表示彩色图I_k的宽，H表示彩色图I_k的高，d表示特征维数；

(3b)在特征图F的基础上得到特征向量积分图IP和特征向量乘积积分图IQ；

(3c)根据特征向量积分图IP、特征向量乘积积分图IQ和候选目标集提取候选目标的特征集其中Cⁱ表示第i个候选目标的特征协方差矩阵；

(4)计算特征权值步骤：

(4a)求取各候选目标特征集与特征模板M之间的距离集其中ρⁱ表示第i个候选目标的特征协方差矩阵与特征模板之间的距离；

(4b)根据距离集计算候选目标的特征权值集其中ωⁱ表示第i个候选目标的特征权值；

(5)目标跟踪窗调整步骤：

(5a)对特征权值最大的候选目标的目标跟踪窗进行缩放，分别产生缩小目标跟踪窗的候选目标和放大目标跟踪窗的候选目标其中β为特征权值最大粒子的序号，下标small和big无具体的物理含义，仅表示所属变量分别为缩小目标跟踪窗的候选目标和放大目标跟踪窗的候选目标；

(5b)分别提取缩放目标跟踪窗的候选目标和的特征协方差矩阵和

(5c)根据和计算和对应的特征权值和

(5d)根据与缩放目标跟踪窗前候选目标的特征权值ω^β的对比结果，更新k时刻的目标跟踪窗B_k，并更新特征权值最大的候选目标的特征协方差矩阵C^β及最大特征权值ω^β；

(6)特征模板更新步骤：

通过计算特征模板M和特征权值最大的粒子对应的特征协方差矩阵C^β的加权对数-欧几里得均值，来对特征模板进行更新；

(7)目标状态更新步骤：

(7a)利用重采样算法，根据特征权值集对k时刻预测粒子集进行重采样，得到k时刻的更新粒子集其中表示k时刻第i个粒子的状态估计值；

(7b)根据k时刻的更新粒子集估计k时刻的目标状态Xk；

(7c)根据k时刻的目标状态X_k和目标跟踪窗B_k，确定出k时刻目标的估计范围Target；

(8)输出步骤：

输出k时刻的目标的估计范围Target，若下一时刻观测信息到达，令k＝k+1，转到步骤(2)进行迭代，否则，目标跟踪过程结束。

本发明具有以下优点：

(1)本发明由于利用协方差算子对目标进行描述，融合多特征，从而对干扰、光照变化、以及目标形变具有较强的鲁棒性；

(2)本发明由于提取色相H、饱和度S、明度V和拉普拉斯算子响应，充分挖掘目标彩色图像的丰富信息，并减少信息冗余，增强特征间的相互独立性，从而实现对目标特征的精确描述，增强了目标跟踪的稳定性。

附图说明

图1是本发明的整体流程图；

图2是本发明方法的特征提取信息比较实验结果示意图；

图3是本发明方法的目标非刚性形变实验结果图；

图4是本发明方法的光照变化实验结果图；

图5是本发明方法的遮挡干扰实验结果图。

具体实施方式

参照图1，本发明的具体实施过程包括以下步骤：

步骤1.初始化粒子、目标跟踪窗和特征模板。

1.1)初始化粒子：令初始时刻k＝1，根据目标的初始状态X₀，产生N个粒子，组成k-1时刻的粒子集其中表示k-1时刻第i个粒子的状态估计值，服从均值为X₀方差为Ψ的高斯分布，X₀为目标的初始状态，Ψ为过程噪声方差，i表示粒子的序号，N表示粒子总数；

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

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

步骤2.目标状态预测。

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

2.1.1)分别得到k时刻每一个粒子的状态预测值：

$q_k^i=p_{k-1}^i+w^i,---1)$

其中，为k时刻第i个粒子的状态预测值，为k-1时刻第i个粒子的状态估计值，wⁱ为状态噪声，服从均值为方差为Ψ的高斯分布；

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

${\left\{q_k^i\right\}}_{i=1}^N=\{q_k^1,q_k^2,\·\·\·,q_k^i,\·\·\·,q_k^N\};---2)$

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

2.2.1)分别计算每一个候选目标区域：

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

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

2.2.2)用步骤2.2.1)中所得的N个候选目标，组成候选目标集

${\left\{P_k^i\right\}}_{i=1}^N=\{P_k^1,P_k^2,\·\·\·,P_k^i,\·\·\·{,P}_k^N\}.---4)$

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

3.1)对于k时刻的W×H维彩色图I_k，提取其W×H×d维特征图F：

F＝{F(x,y)|F(x,y)＝[h(x,y),s(x,y),v(x,y),l(x,y)]^T,1≤x≤W,1≤y≤H}，5)

其中，d为特征维数，d＝4，h(x,y)为彩色图I_k中像素点(x,y)处的色相值，s(x,y)为彩色图I_k中像素点(x,y)处的饱和度值，v(x,y)为彩色图I_k中像素点(x,y)处的明度值，l(x,y)为彩色图I_k中像素点(x,y)处的拉普拉斯算子响应值，T表示向量转置；

3.2)根据特征图F计算得到特征向量积分图IP和特征向量乘积积分图IQ：

3.2.1)根据特征图F计算特征向量积分图IP和特征向量乘积积分图IQ中的点：

$\mathrm{IP}(x^{\text{'}},y^{\text{'}},a)=\underset{1\≤x\≤x^{\text{'}},1\≤y\≤y^{\text{'}}}{\mathrm{\Σ}}F(x,y,a),a=1,\·\·\·,d,---6)$

$\mathrm{IQ}(x^{\text{'}},y^{\text{'}},a,b)=\underset{1\≤x\≤x^{\text{'}},1\≤y\≤y^{\text{'}}}{\mathrm{\Σ}}F(x,y,a)F(x,y,b),a,b=1,\·\·\·,d,---7)$

其中，a为第一组特征序号，b为第二组特征序号，d为特征维数，F(x,y,a)为特征图F中的点，F(x,y,b)为特征图F中的点，IP(x′,y′,a)为特征向量积分图IP中的点，IQ(x′,y′,a,b)为特征向量乘积积分图IQ中的点；

将积分图IP和IQ中的任意一点IP_x,y和IQ_x,y分别表示如下：

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

3.2.2)将步骤3.2.1)所得的点组成集合，得到特征向量积分图IP和特征向量乘积积分图IQ：

IP＝{IP_x,y|1≤x≤W,1≤y≤H}， 10)

IQ＝{IQ_x,y|1≤x≤W,1≤y≤H}， 11)

其中，IP_x,y为特征向量积分图IP中的点，IQ_x,y为特征向量乘积积分图IQ中的点，W表示IP和IQ的宽，H表示积分图IP和IQ的高；

3.3)根据特征向量积分图IP、特征向量乘积积分图IQ和候选目标集计算候选目标的特征集

3.3.1)分别计算每一个候选目标的特征协方差矩阵：

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

其中，Cⁱ为第i个候选目标的特征协方差矩阵，(x′,y′)为第i个候选目标对应矩形区域左上角的顶点坐标，(x′′,y′′)为第i个候选目标对应矩形区域右下角的顶点坐标，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.3.2)用步骤3.3.1)中所得的N个特征协方差矩阵，组成候选目标特征集

${\left\{C^i\right\}}_{i=1}^N=\{C^1,C^2,\·\·\·,C^i,\·\·\·,C^N\}.---13)$

步骤4.求取粒子特征权值。

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

4.1.1)分别计算每一个候选目标特征与特征模板M的之间距离：

${\mathrm{\ρ}}^i={\left(\underset{a=1}{\overset{d}{\mathrm{\Σ}}}{\ln }^2{\mathrm{\λ}}_a\right)}^{1/2},---14)$

其中，ρⁱ表示第i个候选目标特征与特征模板之间的距离，λ_a为特征协方差矩阵Ci和特征模板M的广义特征值：λ_a＝Mx_a(Cⁱx_a)^-1，Cⁱ为第i个候选目标的特征协方差矩阵，x_a为广义特征向量，a表示特征序号，a＝1,2,…,d；

4.1.2)用步骤4.1.1)中所得的N个距离，组成距离集

${\left\{{\mathrm{\ρ}}^i\right\}}_{i=1}^N=\{{\mathrm{\ρ}}^1,{\mathrm{\ρ}}^2,\·\·\·,{\mathrm{\ρ}}^i,\·\·\·{\mathrm{\ρ}}^N\};---15)$

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

4.2.1)分别计算每一个候选目标的特征权值：

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

其中，ωⁱ表示第i个候选目标的特征权值，R为特征观测噪声方差；

4.2.2)用步骤4.2.1)中所求的N个候选目标的特征权值，组成候选目标特征权值集

${\left\{{\mathrm{\ω}}^i\right\}}_{i=1}^N=\{{\mathrm{\ω}}^1,{\mathrm{\ω}}^2,\·\·\·,{\mathrm{\ω}}^i,\·\·\·{\mathrm{\ω}}^N\}.---17)$

步骤5.目标跟踪窗调整。

5.1)对特征权值最大的候选目标的目标跟踪窗进行缩放，产生缩小目标跟踪窗的候选目标区域和放大目标跟踪窗的候选目标区域

$P_{k,\mathrm{small}}^{\mathrm{\β}}=\left\{(x,y)\right||x-{\mathrm{\υ}}_k^{\mathrm{\β}}|\≤\frac{0.95r_{k-1}}{2},|y-v_k^{\mathrm{\β}}|\≤\frac{0.95c_{k-1}}{2}\},---18)$

$P_{k,\mathrm{big}}^{\mathrm{\β}}=\left\{(x,y)\right||x-{\mathrm{\υ}}_k^{\mathrm{\β}}|\≤\frac{1.05r_{k-1}}{2},|y-v_k^{\mathrm{\β}}|\≤\frac{1.05c_{k-1}}{2}\},---19)$

其中，β为特征权值最大粒子的序号，下标small和big无具体的物理含义，仅表示所属变量分别为缩小目标跟踪窗的候选目标和放大目标跟踪窗的候选目标，和分别表示k时刻第β个预测粒子的横坐标和纵坐标，r_k-1和c_k-1分别表示k-1时刻目标跟踪窗的长度和宽度值；

5.2)分别计算缩放目标跟踪窗的候选目标和的特征协方差矩阵和

5.2.1)计算缩小目标跟踪窗的候选目标对应的特征协方差矩阵

$\begin{array}{c}{C}_{\mathrm{small}}^{\mathrm{\β}}=\frac{1}{n_1-1}[{\mathrm{IQ}}_{x_2,y_2}-{\mathrm{IQ}}_{x_1,y_2}-{\mathrm{IQ}}_{x_2,y_1}+{\mathrm{IQ}}_{x_1,y_1}\\ -\frac{1}{n_1}({\mathrm{IP}}_{x_2,y_2}-{\mathrm{IP}}_{x_1,y_2}-{\mathrm{IP}}_{x_2,y_1}+{\mathrm{IP}}_{x_1,y_1}){({\mathrm{IP}}_{x_2,y_2}-{\mathrm{IP}}_{x_1,y_2}-{\mathrm{IP}}_{x_2,y_1}+{\mathrm{IP}}_{x_1,y_1})}^T]\end{array},---20)$

其中，(x₁,y₁)为缩小目标跟踪窗的候选目标对应矩形区域左上角的顶点坐标，(x₂,y₂)为缩小目标跟踪窗的候选目标对应矩形区域右下角的顶点坐标，n₁为区域内像素总数n₁＝(x₂-x₁)·(y₂-y₁)，为特征向量乘积积分图IQ中的不同点，为特征向量积分图IP中的不同点；

5.2.2)计算放大目标跟踪窗的候选目标对应的特征协方差矩

$\begin{array}{c}{C}_{\mathrm{big}}^{\mathrm{\β}}=\frac{1}{n_2-1}[{\mathrm{IQ}}_{x_4,y_4}-{\mathrm{IQ}}_{x_4,y_4}-{\mathrm{IQ}}_{x_2,y_1}+{\mathrm{IQ}}_{x_3,y_3}\\ -\frac{1}{n_2}({\mathrm{IP}}_{x_4,y_4}-{\mathrm{IP}}_{x_3,y_4}+{\mathrm{IP}}_{x_4,y_3}+{\mathrm{IP}}_{x_3,y_3}){({\mathrm{IP}}_{x_4,y_4}-{\mathrm{IP}}_{x_3,y_4}-{\mathrm{IP}}_{x_4,y_4}+{\mathrm{IP}}_{x_3,y_3})}^T]\end{array},---21)$

其中，(x₃,y₃)为缩小目标跟踪窗的候选目标对应矩形区域左上角的顶点坐标，(x₄,y₄)为缩小目标跟踪窗的候选目标对应矩形区域右下角的顶点坐标，n₂为区域内像素总数n₂＝(x₄-x₃)·(y₄-y₃)，为特征向量乘积积分图IQ中的不同点，为特征向量积分图IP中的不同点；

5.3)根据特征协方差矩阵和计算对应的缩放目标跟踪窗的候选目标和的特征权值和

5.3.1)计算特征协方差矩阵和与特征模板M之间的距离ρ_small和ρ_big；

${\mathrm{\ρ}}_{\mathrm{small}}={\left(\underset{a=1}{\overset{d}{\mathrm{\Σ}}}{\ln }^2{\mathrm{\λ}}_{a,\mathrm{small}}\right)}^{1/2},---22)$

${\mathrm{\ρ}}_{\mathrm{big}}={\left(\underset{a=1}{\overset{d}{\mathrm{\Σ}}}{\ln }^2{\mathrm{\λ}}_{a,\mathrm{big}}\right)}^{1/2},---23)$

其中，λ_a,small为特征协方差矩阵和特征模板M的广义特征值：λ_a,big为特征协方差矩阵和特征模板M的广义特征值：x_a,small为λ_a,small对应的广义特征向量，x_a,big为λ_a,big对应的广义特征向量，a表示特征序号，a＝1,2,…,d，下标small和big无具体的物理含义；

5.3.2)根据距离ρ_small和ρ_big计算对应的缩放目标跟踪窗的候选目标和的特征权值和

${\mathrm{\ω}}_{\mathrm{small}}^{\mathrm{\β}}=\frac{1}{\sqrt{2\mathrm{\πR}}}\exp \{-\frac{{\left({\mathrm{\ρ}}_{\mathrm{small}}\right)}^2}{2R}\},---24)$

${\mathrm{\ω}}_{\mathrm{big}}^{\mathrm{\β}}=\frac{1}{\sqrt{2\mathrm{\πR}}}\exp \{-\frac{{\left({\mathrm{\ρ}}_{\mathrm{small}}\right)}^2}{2R}\}---\left(25\right)$

其中，R为特征观测噪声方差，ρ_small为特征协方差矩阵与特征模板M之间的距离，ρ_big为特征协方差矩阵与特征模板M之间的距离；

5.4)根据和与缩放目标跟踪窗前候选目标的特征权值ω^β比较的结果，更新k时刻的目标跟踪窗B_k、特征协方差矩阵C^β及最大特征权值ω^β：

$B_k=\left\{\begin{array}{cc}0.95B_{k-1}& ,{\mathrm{\ω}}_{\mathrm{small}}^{\mathrm{\β}}={\mathrm{\ω}}_{\max }\\ B_{k-1}& ,{\mathrm{\ω}}^{\mathrm{\β}}={\mathrm{\ω}}_{\max }\\ 1.05B_{k-1}& ,{\mathrm{\ω}}_{\mathrm{big}}^{\mathrm{\β}}={\mathrm{\ω}}_{\max }\end{array}\right.,---26)$

$C^{\mathrm{\β}}=\left\{\begin{array}{cc}{C}_{\mathrm{small}}^{\mathrm{\β}}& ,{\mathrm{\ω}}_{\mathrm{small}}^{\mathrm{\β}}={\mathrm{\ω}}_{\max }\\ C^{\mathrm{\β}}& ,{\mathrm{\ω}}^{\mathrm{\β}}={\mathrm{\ω}}_{\max }\\ C_{\mathrm{big}}^{\mathrm{\β}}& ,{\mathrm{\ω}}_{\mathrm{big}}^{\mathrm{\β}}={\mathrm{\ω}}_{\max }\end{array}\right.,---27)$

${\mathrm{\ω}}^{\mathrm{\β}}=\left\{\begin{array}{cc}{\mathrm{\ω}}_{\mathrm{small}}^{\mathrm{\β}}& ,{\mathrm{\ω}}_{\mathrm{small}}^{\mathrm{\β}}={\mathrm{\ω}}_{\max }\\ {\mathrm{\ω}}^{\mathrm{\β}}& ,{\mathrm{\ω}}^{\mathrm{\β}}={\mathrm{\ω}}_{\max }\\ {\mathrm{\ω}}_{\mathrm{big}}^{\mathrm{\β}}& ,{\mathrm{\ω}}_{\mathrm{big}}^{\mathrm{\β}}={\mathrm{\ω}}_{\max }\end{array}\right.,---28)$

其中， ${\mathrm{\ω}}_{\max }=\max \{{\mathrm{\ω}}_{\mathrm{small}}^{\mathrm{\β}},{\mathrm{\ω}}^{\mathrm{\β}},{\mathrm{\ω}}_{\mathrm{big}}^{\mathrm{\β}}\}.$

步骤6.特征模板更新。

6.1)计算特征模板M和特征权值最大粒子的协方差矩阵C^β的加权对数-欧几里得均值M′：

M′＝exp[γlog(M)+(1-γ)log(C^β)]， 29)

其中，β为特征权值最大的粒子序号，γ为更新速率因子,0＜γ＜1。

6.2)用步骤6.1)所得均值M′对特征模板M进行更新：

M＝M′。 30)

步骤7.目标状态更新。

7.1)利用重采样算法，根据特征权值集对k时刻预测粒子集重采样，得到k时刻的更新粒子集

7.2)根据k时刻的更新粒子集估计k时刻的目标状态Xk：

$X_k=\frac{1}{N}\underset{i=1}{\overset{1}{\mathrm{\Σ}}}{p}_k^i---31)$

其中，N表示粒子总数，表示k时刻第i个更新粒子的状态估计值；

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

$T\arg \mathrm{et}=\left\{(x,y)\right||x-x_k|\≤\frac{r_k}{2},|y-y_k|\≤\frac{c_k}{2}\},---32)$

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

步骤8输出和迭代判断。

8.1)输出步骤7.3)所得的目标估计范围Target；

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

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

1.仿真条件。

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

仿真参数：粒子数N＝100，过程噪声方差Ψ＝diag([16,16])，特征观测噪声方差R＝0.1，更新速率因子γ＝0.8。

2.仿真方法。

方法1：本发明方法；

方法2：现有基于联合特征直方图特征提取的粒子滤波跟踪方法，该方法选取r(红色)、g(绿色)和b(蓝色)构成联合颜色直方图；

方法3：现有基于传统协方差算子特征提取的粒子滤波跟踪方法，该方法选取灰度值、x方向的一阶和二阶梯度、y方向的一阶和二阶梯度等5维特征融合构成协方差算子。

3.仿真内容与结果。

仿真1：用所述三种方法，对目标特征进行提取，结果如图2所示，其中：

图2(a)为目标的原始彩色图像；

图2(b)为用方法1提取的一组特征：色相、饱和度、明度和拉普拉斯算子响应；

图2(c)为用方法2提取的一组特征：红色、绿色和蓝色；

图2(d)为用方法3提取的一组特征：灰度、x方向一阶梯度、x方向二阶梯度、y方向一阶梯度和y方向二阶梯度。

从图2可以看出，方法2提取的特征中，红色、绿色和蓝色之间存在较强的相似性，信息冗余大，特征间独立性差；方法3提取的特征中，x方向一阶梯度和二阶梯度，y方向一阶梯度和二阶梯度等特征之间也存在相似性，且特征能量较弱，易受干扰；而本发明提取的特征中，色相、饱和度、明度和拉普拉斯算子响应等特征之间有显著差异，独立性强，信息冗余小，且特征能量较强，抗干扰能力强。

仿真2：用所述三种方法，在目标非刚性形变情况下进行跟踪，结果如图3所示，其中：

图3(a)为用方法1对第1、26、34和67帧的跟踪结果；

图3(b)为用方法2对第1、26、34和67帧的跟踪结果；

图3(c)为用方法3对第1、26、34和67帧的跟踪结果。

从图3可以看出，方法2和方法3对目标的非刚性形变的适应性较差，跟踪结果更容易发生偏离；而本发明的跟踪精度更高。

仿真3：用所述三种方法，在光照变化情况下进行跟踪，结果如图4所示，其中：

图4(a)为用方法1对第1、21、61和81帧的跟踪结果；

图4(b)为用方法2对第1、21、61和81帧的跟踪结果；

图4(c)为用方法3对第1、21、61和81帧的跟踪结果。

从图4可以看出，方法2对光照变化造成目标特征整体发生偏移的适应性较差；现有方法3在光照变暗时的跟踪结果较好，但在光照变亮时跟踪结果发生了较大的偏离；而本发明对光照变化具有较强的鲁棒性。

仿真4：用所述三种方法，在遮挡干扰情况下进行跟踪，结果如图5所示，其中：

图5(a)为用方法1对第1、27、48和63帧的跟踪结果；

图5(b)为用方法2对第1、27、48和63帧的跟踪结果；

图5(c)为用方法3对第1、27、48和63帧的跟踪结果。

从图5可以看出，方法2和方法3在遮挡干扰情况下，均会发生较大程度的偏离；而本发明对遮挡干扰具有较强的鲁棒性。

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

表1

由表1中的统计数据可以看出：在跟踪误差方面，本发明在目标非刚性形变情况下，比方法2和方法3降低了30％；在光照变化情况下，本发明比方法2和方法3降低了60％～80％；在遮挡干扰情况下，本发明比方法2和方法3降低了15％～50％；在运行时间方面，本发明与方法2、方法3相当。综上可以得出，在用时相当的情况下，本发明在跟踪精度和鲁棒性方面均优于方法2和方法3。

Claims (7)

1.一种基于HSV颜色协方差特征的目标跟踪方法，包括：

(1)初始化步骤：

(1a)读入k-1时刻的彩色图I_k-1，根据目标的初始状态，产生k-1时刻的初始粒子集其中，表示k-1时刻第i个粒子的状态估计值，N表示粒子总数，i表示粒子的序号，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时刻预测粒子集和目标跟踪窗B_k-1，确定k时刻候选目标集其中为k时刻第i个候选目标，表示以为中心、B_k-1为长宽所界定出的矩形区域；

(3)候选目标特征提取步骤：

(3a)对于k时刻彩色图I_k，提取其对应的W×H×d维特征图F：

F＝{F(x,y)|F(x,y)＝[h(x,y),s(x,y),v(x,y),l(x,y)]^T,1≤x≤W,1≤y≤H}，

其中W表示彩色图I_k的宽，H表示彩色图I_k的高，d表示特征维数,d＝4，h(x,y)为彩色图I_k中像素点(x,y)处的色相值，s(x,y)为彩色图I_k中像素点(x,y)处的饱和度值，v(x,y)为彩色图I_k中像素点(x,y)处的明度值，l(x,y)为彩色图I_k中像素点(x,y)处的拉普拉斯算子响应值，T表示向量转置；

(3b)在特征图F的基础上得到特征向量积分图IP和特征向量乘积积分图IQ：

(3b1)根据特征图F计算特征向量积分图IP和特征向量乘积积分图IQ中的点：

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

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

其中，a为第一组特征序号，b为第二组特征序号，d为特征维数，F(x,y,a)为特征图F中的点，F(x,y,b)为特征图F中的点，IP(x′,y′,a)为特征向量积分图IP中的点，IQ(x′,y′,a,b)为特征向量乘积积分图IQ中的点；

将积分图IP和IQ中的任意一点IP_x,y和IQ_x,y分别表示如下：

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

(3b2)将步骤(3b1)所得的点组成集合，得到特征向量积分图IP和特征向量乘积积分图IQ：

IP＝{IP_x,y|1≤x≤W,1≤y≤H}，

IQ＝{IQ_x,y|1≤x≤W,1≤y≤H}，

其中，IP_x,y为特征向量积分图IP中的点，IQ_x,y为特征向量乘积积分图IQ中的点，W表示IP和IQ的宽，H表示积分图IP和IQ的高；

(3c)根据特征向量积分图IP、特征向量乘积积分图IQ和候选目标集提取候选目标的特征集

(3c1)分别计算每一个候选目标的特征协方差矩阵：

$\begin{array}{c}{C}^i=\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},$

其中，Cⁱ为第i个候选目标的特征协方差矩阵，(x′,y′)为第i个候选目标对应矩形区域左上角的顶点坐标，(x″,y″)分别为第i个候选目标对应矩形区域右下角的顶点坐标，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中的不同点，T表示向量转置；

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

${\left\{C^i\right\}}_{i=1}^N=\{C^1,C^2,...,C^i,...,C^N\};$

(4)计算特征权值步骤：

(4a)求取各候选目标特征集与特征模板M之间的距离集其中ρⁱ表示第i个候选目标的特征协方差矩阵与特征模板之间的距离；

(4b)根据距离集计算候选目标的特征权值集其中ωⁱ表示第i个候选目标的特征权值；

(5)目标跟踪窗调整步骤：

(5a)对特征权值最大的候选目标的目标跟踪窗进行缩放，分别产生缩小目标跟踪窗的候选目标和放大目标跟踪窗的候选目标其中β为特征权值最大粒子的序号，下标small和big无具体的物理含义，仅表示所属变量分别为缩小目标跟踪窗的候选目标和放大目标跟踪窗的候选目标；

(5b)分别提取缩放目标跟踪窗的候选目标和的特征协方差矩阵和

(5c)根据和计算和对应的特征权值和

(5d)根据与缩放目标跟踪窗前候选目标的特征权值ω^β的对比结果，更新k时刻的目标跟踪窗B_k，并更新特征权值最大的候选目标的特征协方差矩阵C^β及最大特征权值ω^β；

(6)特征模板更新步骤：

通过计算特征模板M和特征权值最大的粒子对应的特征协方差矩阵C^β的加权对数-欧几里得均值，来对特征模板进行更新；

(7)目标状态更新步骤：

(7a)利用重采样算法，根据特征权值集对k时刻预测粒子集进行重采样，得到k时刻的更新粒子集其中表示k时刻第i个粒子的状态估计值；

(7b)根据k时刻的更新粒子集估计k时刻的目标状态X_k；

(7c)根据k时刻的目标状态X_k和目标跟踪窗B_k，确定出k时刻目标的估计范围Target；

(8)输出步骤：

输出k时刻的目标的估计范围Target，若下一时刻观测信息到达，令k＝k+1，转到步骤(2)进行迭代，否则，目标跟踪过程结束。

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

(2b1)分别计算每一个候选目标区域：

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

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

(2b2)用步骤(2b1)中所得的N个候选目标，组成候选目标集

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

3.根据权利要求1所述的方法，其中所述步骤(4a)中计算候选目标特征集与特征模板M之间的距离集通过如下步骤计算：

(4a1)分别计算每一个候选目标特征与特征模板M的之间距离：

${\mathrm{\ρ}}^i={\left(\underset{a=1}{\overset{d}{\Σ}}{\ln }^2{\mathrm{\λ}}_a\right)}^{1/2},$

其中，ρⁱ表示第i个候选目标特征与特征模板之间的距离，λ_a为特征协方差矩阵Cⁱ和特征模板M的广义特征值：λ_a＝Mx_a(Cⁱx_a)^-1，Cⁱ为第i个候选目标的特征协方差矩阵，x_a为广义特征向量，a表示特征序号，a＝1,2,…,d；

(4a2)用步骤(4a1)中所得的N个距离，组成距离集

${\left\{{\mathrm{\ρ}}^i\right\}}_{i=1}^N=\{{\mathrm{\ρ}}^1,{\mathrm{\ρ}}^2,...,{\mathrm{\ρ}}^i,...{\mathrm{\ρ}}^N\}.$

4.根据权利要求1所述的方法，其中所述步骤(4b)中根据距离集计算候选目标特征权值集通过如下步骤进行：

(4b1)分别计算每一个候选目标的特征权值：

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

其中，ωⁱ表示第i个候选目标的特征权值，R为特征观测噪声方差，ρⁱ为第i个候选目标特征与特征模板之间的距离；

(4b2)用步骤(4b1)中所求的N个候选目标的特征权值，组成候选目标特征权值集

${\left\{{\mathrm{\ω}}^i\right\}}_{i=1}^N=\{{\mathrm{\ω}}^1,{\mathrm{\ω}}^2,...,{\mathrm{\ω}}^i,...{\mathrm{\ω}}^N\}.$

5.根据权利要求1所述的方法，其中所述步骤(5d)中更新k时刻的目标跟踪窗B_k、特征协方差矩阵C^β及最大特征权值ω^β，通过如下公式更新：

$B_k=\left\{\begin{array}{cc}0.95B_{k-1},& {\mathrm{\ω}}_{small}^{\mathrm{\β}}={\mathrm{\ω}}_{max}\\ B_{k-1},& {\mathrm{\ω}}^{\mathrm{\β}}={\mathrm{\ω}}_{max}\\ 1.05B_{k-1},& {\mathrm{\ω}}_{big}^{\mathrm{\β}}={\mathrm{\ω}}_{\max }\end{array}\right.,$

$C^{\mathrm{\β}}=\left\{\begin{array}{cc}{C}_{small}^{\mathrm{\β}},& {\mathrm{\ω}}_{small}^{\mathrm{\β}}={\mathrm{\ω}}_{max}\\ C^{\mathrm{\β}},& {\mathrm{\ω}}^{\mathrm{\β}}={\mathrm{\ω}}_{max}\\ C_{big}^{\mathrm{\β}},& {\mathrm{\ω}}_{big}^{\mathrm{\β}}={\mathrm{\ω}}_{\max }\end{array}\right.,$

${\mathrm{\ω}}^{\mathrm{\β}}=\left\{\begin{array}{cc}{\mathrm{\ω}}_{small}^{\mathrm{\β}},& {\mathrm{\ω}}_{small}^{\mathrm{\β}}={\mathrm{\ω}}_{max}\\ {\mathrm{\ω}}^{\mathrm{\β}},& {\mathrm{\ω}}^{\mathrm{\β}}={\mathrm{\ω}}_{max}\\ {\mathrm{\ω}}_{big}^{\mathrm{\β}},& {\mathrm{\ω}}_{big}^{\mathrm{\β}}={\mathrm{\ω}}_{\max }\end{array}\right.,$

其中，

6.根据权利要求1所述的方法，其中所述步骤(7b)中根据k时刻的更新粒子集估计k时刻的目标状态X_k，通过如下公式计算：

$X_k=\frac{1}{N}\underset{i=1}{\overset{N}{\Σ}}{p}_k^i,$

其中N表示样本粒子总数，表示k时刻的第i个更新粒子的状态估计值。

7.根据权利要求1所述的方法，其中所述步骤(7c)中根据k时刻的目标状态X_k和目标跟踪窗B_k，确定k时刻目标的估计范围Target，其通过如下公式确定：

$T\arg et=\left\{(x,y)\right||x-x_k|\≤\frac{r_k}{2},|y-y_k|\≤\frac{c_k}{2}\},$

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