CN102184554A - 基于激活区域背景感知的红外目标跟踪方法 - Google Patents

Abstract

本发明公开了一种基于激活区域背景感知的红外目标跟踪方法，主要解决现有跟踪技术实时性较差、模板易受干扰影响和缺少有效的跟踪窗自适应调整策略的问题。其实现步骤是：在粒子滤波框架下首先通过建立激活区域，划分出候选目标选取和背景感知的有效区域，利用协方差算子提取表观特征；其次通过激活区域背景感知，分辨目标当前所处的场景状态，并提取目标位置和尺度观测集；最后根据场景状态制定模板更新策略，避免将干扰引入模板，并引入目标位置观测与表观特征融合对候选目标进行筛选，以提高跟踪精度，同时利用目标尺度观测实现跟踪窗的自适应调整。本发明具有较强的抗干扰能力，实现了强杂波干扰环境下对红外目标的快速自适应精确跟踪。

Description

基于激活区域背景感知的红外目标跟踪方法

技术领域

本发明属于制导技术领域，涉及红外目标跟踪，具体地说是一种基于激活区域背景感知的红外目标跟踪方法，可用于红外成像制导、监控等军用或民用系统。

背景技术

红外目标跟踪，由于目标图像细节特征少、信噪比低等特点，一直以来是视频跟踪研究中的热点和难点。其中，细节特征少主要是指红外图像缺少颜色和纹理信息，物体之间的可区分性较低，所以红外目标跟踪比可见光更易受到干扰，尤其是形状类似的发热源带来的干扰更严重。此外，低信噪比又会使得在特征提取过程中参入较多错误信息，影响目标跟踪效果，因此对红外目标跟踪方法的抗干扰能力和跟踪精度提出了更高的要求。

目前，红外目标跟踪方法主要有：全局特征匹配、均值漂移和粒子滤波等。其中：

全局特征匹配方法，仅依靠特征而未利用任何空间信息，对类似物造成的干扰跟踪稳定性差，容易产生“跳变”失跟，且需对全局进行匹配，十分耗时；

均值漂移方法，虽然计算简单，运行速度快，但是因具有局部收敛特性而产生明显的跟踪滞后效应，存在跟踪系统误差，且抗干扰能力较差的问题；

粒子滤波方法，对目标运动非线性和噪声非高斯性的情况，虽然具有良好的跟踪精度和稳定性，但是若仅依靠单一特征也只适用于相对简单的场景，而在复杂环境下，跟踪稳定性则会明显下降。

近年来，有学者提出采用多特征融合的方法来提高目标特征描述的抗干扰性能，如Porikli等人发表的论文《Covariance tracking using model update based on Lie Algebra》。这种方法是利用协方差算子将多个特征融合到一个协方差矩阵中，不仅比联合特征直方图降低了计算复杂度，而且具有较强的目标描述能力，非常适合于解决红外目标跟踪中的抗干扰问题。但是这种方法依然存在以下不足：一是协方差算子的快速方法是针对全局的，严重影响了实时性；二是传统的模板更新方法是根据前一个连续时间段内的估计状态来预测当前帧的模板，当存在干扰时，这种不加区分的更新方法显然会将含有干扰的状态引入模板，从而导致偏离，使得在经过干扰物之后的跟踪精度明显下降，甚至发散；三是缺少有效的跟踪窗自适应方法，不能完全解决径向运动或非刚性形变导致的目标尺度变化问题。

发明内容

针对上述问题，本发明在粒子滤波框架下，采用协方差算子提取表观特征，提出了一种基于激活区域背景感知的红外目标跟踪方法，以提高抗干扰能力，实现强杂波干扰环境下对目标的快速自适应精确跟踪。

实现本发明的技术关键是：在跟踪过程中，通过建立激活区域，划分出候选目标选取和背景感知的有效区域，从而在不降低精度的前提下避免全局所带来的庞大计算量，利用协方差算子提取表观特征；其次，通过激活区域背景感知，分辨目标当前所处的场景状态有无干扰，并提取目标位置和尺度观测集；最后，根据场景状态制定模板更新策略，避免将干扰引入特征模板，并引入目标位置观测与表观特征融合对候选目标进行筛选，以提高跟踪精度，同时利用目标尺度观测实现跟踪窗的自适应调整。

本发明的具体实现方法包括：

(1)初始化步骤：

(1a)根据目标的初始状态，产生k-1时刻的初始样本粒子集其中，i表示样本粒子的序号，N表示样本粒子总数，k表示时刻，初始时刻为k＝1，

表示k-1时刻第i个粒子的状态估计列向量；

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

(1c)初始化特征模板库J＝1，其中，j表示特征模板的序号，J表示模板总数，且具有上限J＜J_max，模板库在k时刻的更新位置μ_k∈[0，J_max]，μ₁＝2，设场景状态Λ＝0，即假设无干扰；

(2)建立激活区域步骤：

(2a)对目标状态进行预测：即通过粒子传递得到预测粒子集

表示k时刻第i个粒子的状态预测列向量，根据粒子的状态预测列向量和跟踪窗的状态预测列向量B_k-1|k就可以确定出候选目标集

其中

为候选目标，即表示以第i个粒子

为中心、B_k-1|k为长宽所界定出的矩形区域，B_k-1|k＝(r_k-1|k，c_k-1|k)^T，r_k-1|k和c_k-1|k分别表示k时刻跟踪窗长和宽的预测值；

(2b)由候选目标集确定出激活区域AR；

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

(3a)设激活区域AR为W×H维灰度图，提取其对应的W×H×d维特征图F_AR，其中d表示特征维数；

(3b)在F_AR的基础上构造特征向量积分图IP和特征协方差矩阵积分图IQ；

(3c)根据特征向量积分图IP、特征协方差矩阵积分图IQ和候选目标集

提取候选目标的特征集

其中每个候选目标的特征又是由其分块特征构成的集合

t表示候选目标区域分块序号，分块总数为5，

表示第i个候选目标的第t个分块所对应的特征协方差矩阵；

(4)激活区域背景感知步骤：

采用迭代Otsu法对激活区域AR进行分割，通过对分割标记图进行计算，得到目标位置观测集尺度观测集

和场景状态Λ，其中，l表示疑似目标块的序号，L_k为疑似目标块数，场景状态Λ为1时表示有干扰，为0时表示无干扰，下标center和rc无具体的物理含义，仅表示所属变量分别为位置观测和尺度观测；

(5)计算权值步骤：

(5a)计算特征权值：通过计算模板库中所有元素的Log-Euclidean均值，得到当前k时刻的模板预测值：

然后，求取各候选目标特征Cⁱ与模板预测值

的距离，进而计算各候选目标的特征权值下标f表示所属变量为特征权值；

(5b)计算空间权值：在无干扰时，去除目标位置观测集

中的杂波，得到位置观测量Z_k，进而计算各候选目标的空间权值

否则，认为位置观测不可信，空间权值均置1，下标p表示所属变量为空间权值；

(5c)特征权值和空间权值进行融合，得到融合权值

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

根据最大似然准则对模板库进行更新，即用融合权值最大粒子的协方差矩阵C^β更新模板库中的相应元素

再确定下一时刻的更新位置μ_k+1，其中β为融合权值最大粒子的序号，μ_k表示k时刻的模板库更新位置；

(7)跟踪窗自适应调整步骤：

(7a)在无干扰时，去除目标尺度观测集中的杂波，得到尺度观测量Z_rc；

(7b)采用卡尔曼滤波方法对跟踪窗进行自适应调整，估计当前k时刻跟踪窗状态列向量B_k|k，否则，认为尺度观测不可信，跟踪窗保持不变；

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

利用融合权值

对预测粒子集

进行重采样，得到更新粒子集

进而估计当前k时刻的目标状态X_k，再结合跟踪窗状态估计列向量B_k|k，确定出目标的估计范围Target；

(9)输出步骤：

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

本发明具有以下优点：

(1)本发明由于利用协方差算子对目标进行描述，融合多特征，充分挖掘该区域特征的方差及互信息，对干扰、光照变化、以及目标形变有较强的鲁棒性；

(2)本发明由于通过建立激活区域以划分出候选目标选取和背景感知的有效区域，从而在不降低精度的前提下避免全局所带来的庞大计算量；

(3)本发明由于通过背景感知以分辨目标当前所处的场景状态有无干扰，进而避免将干扰引入特征模板，并引入目标位置观测与表观特征融合对候选目标进行筛选，从而有效提高了跟踪精度；

(4)本发明由于提取目标尺度观测，并采用卡尔曼滤波方法实现了跟踪窗的自适应精确调整。

附图说明

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

图2是现有积分图中任意矩形区域示意图；

图3是本发明方法中的激活区域示意图；

图4是本发明方法中的目标区域划分示意图；

图5是本发明方法中的激活区域背景感知示意图；

图6是本发明方法的跟踪窗自适应调整实验结果图；

图7是本发明方法的交叉干扰实验结果图；

图8是本发明方法的同向干扰实验结果图；

图9是本发明方法的复杂背景干扰实验结果图；

图10是本发明方法的海面背景实验结果图；

图11是本发明方法的天空背景实验结果图。

具体实施方式

一、基础理论介绍

1.协方差算子

协方差算子，先是提取一个区域内所有像素的底层视觉特征，例如灰度、梯度等，然后计算协方差矩阵，通过协方差矩阵中的对角元素来反映该区域内特征的自相关，非对角元素反应特征间的互相关信息。这种融合机制可充分挖掘和利用现有特征信息，非常适合于红外目标跟踪。设I为W×H维灰度图，F是从I中提取的W×H×d维特征图，所以F中点(x，y)处对应一个d维的特征列向量：

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

其中，d为特征维数，d＝5，I(x，y)为点(x，y)处的灰度值，和

分别表示点(x，y)处x方向的一阶和二阶梯度，和

分别表示点(x，y)处y方向的一阶和二阶梯度。

对于任意矩形区域

可以用该区域对应的d×d维特征协方差矩阵C_R来描述：

$C_R=\frac{1}{n-1}\underset{m=1}{\overset{n}{\mathrm{\Σ}}}(z_m-\mathrm{\μ}){(z_m-\mathrm{\μ})}^T---2)$

其中， 表示d维实数列向量，m表示区域内像素的序号，z_m为第m个像素所对应的d维特征列向量，n表示区域内像素总数，μ为特征均值列向量，

2.积分图

当计算大量区域的协方差矩阵时，由于这些区域通常具有重叠部分，使得同一位置点在计算其隶属不同区域的协方差矩阵时被多次重复统计，造成计算冗余。对于矩形区域的‘累加和’，积分图可以有效减少计算量，其复杂度正比于区域的大小。在得到原图对应的积分图之后，统计积分图内任意矩形区域的时间相同。对于灰度图I，其积分图定义为：

$\mathrm{II}(x^{\′},y^{\′})=\underset{1\≤x\≤x^{\′},1\≤y\≤y^{\′}}{\mathrm{\Σ}}I(x,y)---3)$

R为积分图II中任意矩形区域，如图2所示。区域R内所有元素的‘累加和’为：

R(x′，y′；x″，y″)＝II(x″，y″)-II(x″，y′)-II(x′，y″)+II(x′，y′)4)

根据式2)，协方差矩阵中的元素C_R(a，b)定义为：

$C_R(a,b)=\frac{1}{n-1}\underset{m=1}{\overset{n}{\mathrm{\Σ}}}(z_m\left(a\right)-\mathrm{\μ}\left(a\right))(z_m\left(b\right)-\mathrm{\μ}\left(b\right))---5)$

将式5)中的μ展开并化简，可得：

$C_R(a,b)=\frac{1}{n-1}[\underset{m=1}{\overset{n}{\mathrm{\Σ}}}{z}_m\left(a\right)z_m\left(b\right)-\frac{1}{n}\underset{m=1}{\overset{n}{\mathrm{\Σ}}}{z}_m\left(a\right)\underset{m=1}{\overset{n}{\mathrm{\Σ}}}{z}_m\left(b\right)]---6)$

可以看出，式6)由

和

计算得来，而这两项可通过积分图来快速获取，所以在特征图F的基础之上构造特征向量积分图IP和特征协方差矩阵积分图IQ：

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

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

对于积分图IP和IQ中的任意一点(x，y)，分别为以下形式：

IP_x，y＝[IP(x，y，1)L  IP(x，y，d)]^T    9)

${\mathrm{IQ}}_{x,y}=\left[\begin{array}{ccc}\mathrm{IQ}(x,y,1,1)& L& \mathrm{IQ}(c,y,1,d)\\ M& O& M\\ \mathrm{IQ}(x,y,d,1)& L& \mathrm{IQ}(x,y,d,d)\end{array}\right]---10)$

因此，任意区域R(x′，y′；x″，y″)的协方差算子可由下式计算得到：

$C_R(x^{\′},y^{\′};x^{\′\′},y^{\′\′})=\frac{1}{n-1}[{\mathrm{IQ}}_{x^{\′\′},y^{\′\′}}-{\mathrm{IQ}}_{x^{\′},y^{\′\′}}-{\mathrm{IQ}}_{x^{\′\′},y^{\′}}+{\mathrm{IQ}}_{x^{\′},y^{\′}}]$ 11)

$-\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]$

其中，n表示区域内像素总数，n＝(x″-x′)·(y″-y′)。

二、基于激活区域背景感知的红外目标跟踪方法

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

步骤1.初始化样本粒子、跟踪窗和特征模板库

1.1)初始化粒子，令初始时刻k＝1，首先根据目标的初始状态X₀，产生N个样本粒子 其中，i表示样本粒子的序号，

表示k-1时刻第i个粒子的状态估计列向量，N(X₀，∑)表示以X₀为均值∑为方差的高斯分布，X₀目标初始位置，∑为过程噪声方差；

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

1.3)初始化特征模板库J＝1，其中，j表示特征模板的序号，J表示模板总数，且具有上限J＜J_max，模板库在k时刻的更新位置μ_k∈[0，J_max]，μ₁＝2，设场景状态Λ＝0，即假设无干扰。

步骤2.目标状态预测及激活区域AR的确定

2.1)通过粒子传递得到k时刻粒子的状态预测列向量：

$p_{k-1|k}^i=p_{k-1|k-1}^i+w^i---12)$

其中，wⁱ为状态噪声，

表示以

为均值∑为方差的高斯分布，根据粒子的状态预测列向量

和跟踪窗的状态预测列向量B_k-1|k就可以确定出候选目标集

其中

表示以第i个粒子

为中心、B_k-1|k为长宽所界定出的矩形区域：

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

其中，

和

分别表示k时刻预测粒子

的横坐标和纵坐标，r_k-1|k和c_k-1|k分别表示k时刻跟踪窗长和宽的预测值；

2.2)建立激活区域AR

在整幅图像I中，一定存在某个矩形区域Γ，满足

即所有候选目标区域的并集为Γ的子集，当Γ取最小值时，如图3中虚线框所围成的矩形区域，即为激活区域AR：

$\mathrm{AR}=\min \left\{\mathrm{\Γ}\right|\underset{i=1}{\overset{N}{U}}{P}_{B_{k-1|k}}^i\⋐\mathrm{\Γ}\}---14)$

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

3.1)提取特征图

激活区域AR为W×H维灰度图，F_AR是从AR中提取的W×H×d维特征图，所以F_AR中点(x，y)处对应一个d维的特征列向量：

$F_{\mathrm{AR}}(x,y)={[\mathrm{AR}(x,y),\frac{\∂\mathrm{AR}(x,y)}{\∂x},\frac{\∂\mathrm{AR}(x,y)}{\∂y},\frac{{\∂}^2\mathrm{AR}(x,y)}{{\∂x}^2},\frac{{\∂}^2\mathrm{AR}(x,y)}{{\∂y}^2}]}^T---15)$

其中，d为特征维数，d＝5，AR(x，y)为点(x，y)处的灰度值，和

分别表示点(x，y)处x方向的一阶和二阶梯度，

和

分别表示点(x，y)处y方向的一阶和二阶梯度；

3.2)提取候选目标的协方差算子

构造特征向量积分图IP和特征协方差矩阵积分图IQ，根据候选目标集

即可提取各候选目标的特征集

作为其表观特征的描述，其中考虑到对目标特征描述的抗局部遮挡性能，将每个候选目标区域分为5个子区域，如图4所示，图4中未被遮挡区域是协方差的统计范围，其中图4(a)为候选目标的整体，图4(b)为候选目标的左半部，图4(c)为候选目标的右半部，图4(d)为候选目标的上半部，图4(e)为候选目标的下半部，根据相应的统计区域得到每个候选目标的特征描述集

具体计算过程与上述“基础理论介绍”相同。

步骤4.激活区域背景感知

跟踪背景感知，主要是针对激活区域进行的，参照图5，本步骤的具体实现如下：

4.1)采用迭代Otsu法对激活区域进行目标分割，得到一幅对应的标记图；

4.2)统计标记图中分割块数L′、各分割块的像素个数M_l和对应坐标集

其中，l表示分割块的序号，l＝1，2，L，L′，m表示分割块内像素的序号，m＝1，2，L，M_l；

4.3)删除M_l大于像素个数上限M_sup或小于像素个数下限M_inf的分割块，及其对应坐标集

得到疑似目标块的个数：L_k＝L′-L_del，其中L_del为删除块数；

4.4)提取各疑似目标块的坐标极值：

和块平均像素个数

其中，l＝1，2，L，L_k；

4.5)提取各块的中心坐标：

构造目标位置观测集为：

其中

表示向下取整，T表示向量转置，下标center表示所属变量为位置观测；

4.6)提取各块的长参数：

和宽参数：构造目标尺度观测集为：

其中res为跟踪窗余量，res＝10，下标rc表示所属变量为尺度观测；

4.7)根据疑似目标块的个数L_k和块平均像素个数

的变化率，修正干扰状态：若L_k＜L_k-1、

且此前Λ＝0，则修正干扰状态为有干扰Λ＝1；若L_k＞L_k-1、

且此前Λ＝1，则修正干扰状态为无干扰Λ＝0，其中α为变化率阈值，α＝0.7。

步骤5.求取粒子权值

5.1)特征权值

通过计算模板库中所有元素的Log-Euclidean均值，得到当前k时刻模板的各个分量的预测值：

${\stackrel{\‾}{T}}_t=\exp \left\{\frac{1}{J}\underset{j=1}{\overset{J}{\mathrm{\Σ}}}\log \left(T_t^j\right)\right\},t=\mathrm{1,2},L,5---16)$

从而构成模板 $\stackrel{\‾}{T}={\left\{{\stackrel{\‾}{T}}_t\right\}}_{t=1}^5;$

计算特征权值：

${\mathrm{\ω}}_f^i=\frac{1}{\sqrt{2\mathrm{\π}{R}_f}}\exp \{-\frac{{\left(D_f^i\right)}^2}{2R_f}\}---17)$

其中，R_f为特征观测噪声方差，为候选目标与模板各对应块距离融合得到的总体距离：

$D_f^i=\underset{\mathrm{\γ}}{\min }[\underset{t=1}{\overset{5}{\mathrm{\Σ}}}\mathrm{\ρ}(C_t^i,{\stackrel{\‾}{T}}_t)-\mathrm{\ρ}(C_{\mathrm{\γ}}^i,{\stackrel{\‾}{T}}_{\mathrm{\γ}})]---18)$

式18)中协方差矩阵间的距离

可计算如下：

$\mathrm{\ρ}(C_t^i,{\stackrel{\‾}{T}}_t)={\left(\underset{a=1}{\overset{d}{\mathrm{\Σ}}}{\ln }^2{\mathrm{\λ}}_a(C_t^i,{\stackrel{\‾}{T}}_t)\right)}^{1/2}---19)$

其中，

为协方差矩阵

和

的广义特征值：

a＝1，2，L，d，x_a为对应的特征相向量，

5.2)空间权值

为了去除位置观测集中的杂波，设k-1时刻目标位置估计值为X_k-1，那么认为在L₂范数意义下与上一时刻估计位置X_k-1距离最小的点即为当前k时刻的目标位置观测量：

$Z_k=Z_{\mathrm{center}}^{\mathrm{\θ}}---20)$

其中，

l表示分割块的序号，l＝1，2，L，L_k，L_k为疑似目标块的个数；

在有干扰时，认为目标位置观测不可信，将空间权值均置1；在无干扰时，将k时刻第i个预测粒子

与目标位置观测量Z_k的似然，作为对应粒子的空间权值 ${\mathrm{\ω}}_p^i:$

${\mathrm{\ω}}_p^i=\left\{\begin{array}{cc}\frac{1}{\sqrt{2\mathrm{\π}{R}_p}}\exp \{-\frac{{\left|\right|p_{k-1|k}^i-Z_k\left|\right|}_2^2}{2R_p}\}& ,\mathrm{\Λ}=0\\ 1& ,\mathrm{\Λ}=1\end{array}\right.---21)$

其中，R_p为空间观测噪声方差；

5.3)特征权值

和空间权值

进行融合，得到融合权值ωⁱ：

${\mathrm{\ω}}^i={\mathrm{\ω}}_f^i\·{\mathrm{\ω}}_p^i---22)$

步骤6.特征模板库的更新策略

6.1)根据最大似然准则对模板库进行更新，即用融合权值最大粒子的协方差矩阵C^β更新模板库中的相应元素

$T^{{\mathrm{\μ}}_k}=C^{\mathrm{\β}}---23)$

其中，

μ_k表示k时刻的模板库更新位置；

6.2)求取k+1时刻模板库的更新位置μ_k+1：

${\mathrm{\μ}}_{k+1}=\left\{\begin{array}{ccc}{\mathrm{\μ}}_k+1& ,& {\mathrm{\μ}}_k\≤J_{\max }-1,\mathrm{\Λ}=0\\ 1& ,& {\mathrm{\μ}}_k=J_{\max },\mathrm{\Λ}=0\\ {\mathrm{\μ}}_k& ,& \mathrm{\Λ}=1\end{array}\right.---24)$

步骤7.跟踪窗自适应调整

7.1)在无干扰时，去除目标尺度观测集

中的杂波，得到尺度观测量Z_rc：

$Z_{\mathrm{rc}}=Z_{\mathrm{rc}}^{\mathrm{\λ}}---25)$

其中，

l表示分割块的序号，l＝1，2，L，L_k，L_k为疑似目标块的个数；

7.2)采用卡尔曼滤波方法对跟踪窗进行自适应调整，估计当前k时刻跟踪窗状态列向量B_k|k：

7.2.1)利用卡尔曼滤波预测方程，得到k时刻跟踪窗状态预测列向量B_k-1|k：

B_k-1|k＝F_rcB_k-1|k-1

其中，B_k-1|k-1表示k-1时刻跟踪窗状态估计列向量，B_k-1|k-1＝(r_k-1|k-1，c_k-1|k-1)^T，r_k-1|k-1和c_k-1|k-1分别表示k-1时刻跟踪窗长和宽的估计值，F_rc表示状态转移矩阵，F_rc＝Id，Id表示单位矩阵，即状态转移矩阵F_rc取作单位矩阵；

7.2.2)利用卡尔曼滤波更新方程，得到k时刻跟踪窗状态估计列向量B_k|k：

$B_{k|k}=\left\{\begin{array}{cc}{B}_{k-1|k}+K_{\mathrm{rc}}(Z_{\mathrm{rc}}-B_{k-1|k})& ,\mathrm{\Λ}=0\\ B_{k-1|k}& ,\mathrm{\Λ}=1\end{array}\right.---26)$

其中，K_rc表示增益矩阵，计算式如下：

$K_{\mathrm{rc}}=P_{\mathrm{rc}}^{k-1}{H}_{\mathrm{rc}}^T{(R_{\mathrm{rc}}+H_{\mathrm{rc}}{P}_{\mathrm{rc}}^{k-1}{H}_{\mathrm{rc}})}^{-1}---27)$

其中，H_rc表示观测矩阵，H_rc＝Id，R_rc表示观测噪声方差矩阵，R_rc＝σId，σ为尺度观测相对误差的放大量，σ＝exp{η||Z_rc-B_k-1k||₂/||B_k-1k||₂}，η表示放大系数，η＝50，

表示k-1时刻的状态协方差矩阵，其初始化为

迭代式为：

$P_{\mathrm{rc}}^k=(\mathrm{Id}-K_{\mathrm{rc}}{H}_{\mathrm{rc}})P_{\mathrm{rc}}^{k-1}+Q_{\mathrm{rc}}---28)$

其中，Q_rc表示状态噪声方差矩阵。

步骤8.目标状态更新

利用融合权值

对粒子进行系统重采样得到更新粒子集

进而估计当前帧目标位置：

$X_k=\frac{1}{N}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{p}_{k|k}^i---29)$

结合当前k时刻目标位置X_k和跟踪窗状态估计列向量B_k|k，确定出目标的估计范围：

$T\arg \mathrm{et}=\left\{(x,y)\right||x-x_k|\≤\frac{r_{k|k}}{2},|y-y_k|\≤\frac{c_{k|k}}{2}\}---30)$

其中，x_k和y_k分别表示k时刻目标位置X_k的横坐标和纵坐标，r_k|k和c_k|k分别表示k时刻跟踪窗长和宽的估计值。

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

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

1.仿真条件

实验环境：Intel Pentium D CPU 2.8Ghz，1GB内存，Matlab7.4仿真实验平台。实验数据：图6～图8中的数据来自Terravic Research Infrared Database中选取的三组红外视频序列；图9～图11为实测数据。实验中粒子数N＝100，过程噪声方差∑＝diag([16，16])，特征观测噪声方差R_f＝0.1，空间观测噪声方差R_p＝0.1，跟踪窗过程噪声方差Q_rc＝diag([10，10])，模板库元素数量上限J_m＝50。

仿真实验中，本发明将与以下两种传统方法进行对比：

方法一：是普通基于协方差算子特征提取并采用传统自适应跟踪窗的粒子滤波跟踪方法，该方法对跟踪窗分别采用95％、100％和105％的3种缩放比例进行动态调整；

方法二：普通基于协方差算子特征提取并采用固定窗的粒子滤波跟踪方法。

2.仿真内容与结果

仿真一：对比本发明和方法一，在目标发生径向运动或非刚性形变时的跟踪性能和跟踪窗自适应调整效果，结果如图6所示，其中图6(a)为用本发明对第1、81、123和150帧的跟踪结果，图6(b)为用现有方法一对第1、81、123和150帧的跟踪结果。

仿真二：对比本发明、方法一和方法二，在目标遭受强杂波干扰情况下的跟踪性能。

交叉干扰情况的跟踪结果如图7所示，其中图7(a)为用本发明对第1、55、73和100帧的跟踪结果，图7(b)为用现有方法一对第1、55、73和100帧的跟踪结果，图7(c)为用现有方法一对第1、55、73和100帧的跟踪结果。

同向干扰情况的跟踪结果如图8所示，其中图8(a)为用本发明对第1、32、115和230帧的跟踪结果，图8(b)为用现有方法一对第1、32、115和230帧的跟踪结果，图8(c)为用现有方法一对第1、32、115和230帧的跟踪结果。

复杂背景干扰情况的跟踪结果如图9所示，其中图9(a)为用本发明对第1、25、40和46帧的跟踪结果，图9(b)为用现有方法一对第1、25、40和46帧的跟踪结果，图9(c)为用现有方法一对第1、25、40和46帧的跟踪结果。

仿真三：对比本发明、方法一和方法二，在海面背景和天空背景情况下的跟踪性能。

海面背景情况的跟踪结果如图10所示，其中图10(a)为用本发明对第1、40、80和100帧的跟踪结果，图10(b)为用现有方法一对第1、40、80和100帧的跟踪结果，图10(c)为用现有方法一对第1、40、80和100帧的跟踪结果。

天空背景情况的跟踪结果如图11所示，其中图11(a)为用本发明对第1、20、34和36帧的跟踪结果，图11(b)为用现有方法一对第1、20、34和36帧的跟踪结果，图11(c)为用现有方法一对第1、20、34和36帧的跟踪结果。

从图6中可以看出，现有方法一对于目标缩小的情况，跟踪框引入了过多的背景；而对于目标变大的情况，跟踪框又不能充分包含目标，跟踪精度明显低于本发明。

从图7、图8和图9中可以看出，现有方法一抗干扰能力较弱，容易误将干扰纳入目标范围，跟踪精度明显低于本发明；现有方法二对尺度基本不变的目标，抗干扰能力较强，但终究不能解决目标尺度变化问题，且跟踪精度仍低于本发明。

从图10和图11中可以看出，三种方法均能实现对目标的跟踪，而本发明的跟踪精度更高。

对图6～图11中的场景分别进行100次蒙特卡罗实验，统计平均跟踪误差Err、平均失跟率TLR及平均每帧运行时间RT，结果如表1所示，表中“——”表示该方法不能进行此项统计。

表1

由表1中的统计数据可以看出：在跟踪误差方面，本发明比方法一至少降低了50％，比方法二至少降低了10％，明显提高了跟踪精度；在失跟率方面，本发明也明显低于方法一和方法二；在运行时间方面，本发明仅为方法一和方法二的10％左右。综上可以得出，本发明在跟踪精度和实时性方面均优于传统方法。

Claims (6)

1.基于激活区域背景感知的红外目标跟踪方法，包括：

(1)初始化步骤：

(1a)根据目标的初始状态，产生k-1时刻的初始样本粒子集

其中，i表示样本粒子的序号，N表示样本粒子总数，k表示时刻，初始时刻为k＝1，

表示k-1时刻第i个粒子的状态估计列向量；

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

(1c)初始化特征模板库

J＝1，其中，j表示特征模板的序号，J表示模板总数，且具有上限J＜J_max，模板库在k时刻的更新位置μ_k∈[0，J_max]，μ₁＝2，设场景状态Λ＝0，即假设无干扰；

(2)建立激活区域步骤：

(2a)对目标状态进行预测：即通过粒子传递得到预测粒子集

表示k时刻第i个粒子的状态预测列向量，根据粒子的状态预测列向量

和跟踪窗的状态预测列向量B_k-1|k就可以确定出候选目标集

其中

为候选目标，即表示以第i个粒子为中心、B_k-1|k为长宽所界定出的矩形区域，B_k-1|k＝(r_k-1|k，c_k-1|k)^T，r_k-1|k和c_k-1|k分别表示k时刻跟踪窗长和宽的预测值；

(2b)由候选目标集确定出激活区域AR；

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

(3a)设激活区域AR为W×H维灰度图，提取其对应的W×H×d维特征图F_AR，其中d表示特征维数；

(3b)在F_AR的基础上构造特征向量积分图IP和特征协方差矩阵积分图IQ；

(3c)根据特征向量积分图IP、特征协方差矩阵积分图IQ和候选目标集

提取候选目标的特征集

其中每个候选目标的特征又是由其分块特征构成的集合

t表示候选目标区域分块序号，分块总数为5，

表示第i个候选目标的第t个分块所对应的特征协方差矩阵；

(4)激活区域背景感知步骤：

采用迭代Otsu法对激活区域AR进行分割，通过对分割标记图进行计算，得到目标位置观测集

尺度观测集

和场景状态Λ，其中，l表示疑似目标块的序号，L_k为疑似目标块数，场景状态Λ为1时表示有干扰，为0时表示无干扰，下标center和rc无具体的物理含义，仅表示所属变量分别为位置观测和尺度观测；

(5)计算权值步骤：

(5a)计算特征权值：通过计算模板库中所有元素的Log-Euclidean均值，得到当前k时刻的模板预测值：

然后，求取各候选目标特征Cⁱ与模板预测值

的距离，进而计算各候选目标的特征权值

下标f表示所属变量为特征权值；

(5b)计算空间权值：在无干扰时，去除目标位置观测集中的杂波，得到位置观测量Z_k，进而计算各候选目标的空间权值

否则，认为位置观测不可信，空间权值均置1，下标p表示所属变量为空间权值；

(5c)特征权值和空间权值进行融合，得到融合权值

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

根据最大似然准则对模板库进行更新，即用融合权值最大粒子的协方差矩阵C^β更新模板库中的相应元素

再确定下一时刻的更新位置μ_k+1，其中β为融合权值最大粒子的序号，μ_k表示k时刻的模板库更新位置；

(7)跟踪窗自适应调整步骤：

(7a)在无干扰时，去除目标尺度观测集

中的杂波，得到尺度观测量Z_rc；

(7b)采用卡尔曼滤波方法对跟踪窗进行自适应调整，估计当前k时刻跟踪窗状态列向量B_k|k，否则，认为尺度观测不可信，跟踪窗保持不变；

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

利用融合权值

对预测粒子集

进行重采样，得到更新粒子集

进而估计当前k时刻的目标状态X_k，再结合跟踪窗状态估计列向量B_k|k，确定出目标的估计范围Target；

(9)输出步骤：

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

2.根据权利要求1所述的目标跟踪方法，其中步骤(2b)所述的由候选目标集确定出激活区域AR，按下式计算得到：

$\mathrm{AR}=\min \left\{\mathrm{\Γ}\right|\underset{i=1}{\overset{N}{U}}{P}_{B_{k-1|k}}^i\⋐\mathrm{\Γ}\}$

其中，Γ表示整幅图像I中的任意矩形区域，

为候选目标，即表示以第i个粒子

为中心、B_k-1|k为跟踪窗长宽所界定出的矩形区域，B_k-1|k＝(r_k-1|k，c_k-1|k)^T，r_k-1|k和c_k-1|k分别表示k时刻跟踪窗长和宽的预测值。

3.根据权利要求1所述的目标跟踪方法，其中步骤(4)所述的计算得到目标位置观测集

尺度观测集

和场景状态Λ，按如下步骤进行：

(4a)统计分割块数L′、各分割块的像素个数M_l和对应坐标集

其中，l表示分割块的序号，l＝1，2，L，L′，m表示分割块内像素的序号，m＝1，2，L，M_l；

(4b)删除M_l大于像素个数上限M_sup或小于像素个数下限M_inf的分割块，及其对应坐标集

得到疑似目标块的个数：L_k＝L′-L_del，其中L_del为删除块数；

(4c)提取各疑似目标块的坐标极值：

和块平均像素个数其中，l＝1，2，L，L_k；

(4d)提取各块的中心坐标：

构造目标位置观测集为：

其中表示向下取整，T表示向量转置，下标center表示所属变量为位置观测；

(4e)提取各块的长参数：

和宽参数：构造目标尺度观测集为：其中res为跟踪窗余量，res＝10，下标rc表示所属变量为尺度观测；

(4f)根据疑似目标块的个数L_k和块平均像素个数

的变化率，修正干扰状态：

若L_k＜L_k-1、

且此前Λ＝0，则修正干扰状态为有干扰Λ＝1；

若L_k＞L_k-1、

且此前Λ＝1，则修正干扰状态为无干扰Λ＝0，其中α为变化率阈值，α＝0.7。

4.根据权利要求1所述的目标跟踪方法，其中步骤(5b)所述的计算空间权值，按如下步骤进行：

(5b1)去除目标位置观测集中的杂波，获得目标位置观测量：

$Z_k=Z_{\mathrm{center}}^{\mathrm{\θ}}$

其中，

l表示分割块的序号，l＝1，2，L，L_k，L_k为疑似目标块的个数，X_k-1为k-1时刻目标位置估计值；

(5b2)在有干扰时，认为目标位置观测不可信，将空间权值均置1；在无干扰时，将k时刻第i个预测粒子与目标位置观测量Z_k的似然，作为对应粒子的空间权值

${\mathrm{\ω}}_p^i=\left\{\begin{array}{cc}\frac{1}{\sqrt{2\mathrm{\π}{R}_p}}\exp \{-\frac{{\left|\right|p_{k-1|k}^i-Z_k\left|\right|}_2^2}{2R_p}\}& ,\mathrm{\Λ}=0\\ 1& ,\mathrm{\Λ}=1\end{array}\right.$

其中，R_p为空间观测噪声方差。

5.根据权利要求1所述的目标跟踪方法，其中步骤(7a)所述的在无干扰时，去除目标尺度观测集

中的杂波，得到尺度观测量Z_rc，按下式计算得到：

$Z_{\mathrm{rc}}=Z_{\mathrm{rc}}^{\mathrm{\λ}}$

其中，

l表示分割块的序号，l＝1，2，L，L_k，L_k为疑似目标块的个数，B_k-1|k＝(r_k-1|k，c_k-1|k)^T，r_k-1|k和c_k-1|k分别表示k时刻跟踪窗长和宽的预测值。

6.根据权利要求1所述的目标跟踪方法，其中步骤(7b)所述的采用卡尔曼滤波方法对跟踪窗进行自适应调整，估计当前k时刻跟踪窗状态列向量B_k|k，按如下步骤进行：

(7b1)利用卡尔曼滤波预测方程，得到k时刻跟踪窗状态预测列向量B_k-1|k：

B_k-1|k＝F_rcB_k-1|k-1

其中，B_k-1|k-1表示k-1时刻跟踪窗状态估计列向量，B_k-1|k-1＝(r_k-1|k-1，c_k-1|k-1)^T，r_k-1|k-1和c_k-1|k-1分别表示k-1时刻跟踪窗长和宽的估计值，F_rc表示状态转移矩阵，F_rc＝Id，Id表示单位矩阵，即状态转移矩阵F_rc取作单位矩阵；

(7b2)利用卡尔曼滤波更新方程，得到k时刻跟踪窗状态估计列向量B_k|k：

$B_{k|k}=\left\{\begin{array}{cc}{B}_{k-1|k}+K_{\mathrm{rc}}(Z_{\mathrm{rc}}-B_{k-1|k})& ,\mathrm{\Λ}=0\\ B_{k-1|k}& ,\mathrm{\Λ}=1\end{array}\right.$

其中，K_rc表示增益矩阵，

H_rc表示观测矩阵，H_rc＝Id，

表示k-1时刻的状态协方差矩阵，其初始化为

迭代式为

Q_rc表示状态噪声方差矩阵，R_rc表示观测噪声方差矩阵，R_rc＝σId，σ为尺度观测相对误差的放大量，σ＝exp{η||Z_rc-B_k-1|k||₂/||B_k-1k||₂}，η表示放大系数，η＝50。

Images