CN110910425A - 一种针对抵近飞行过程的目标跟踪方法 - Google Patents

Abstract

本发明公开一种针对抵近飞行过程的目标跟踪方法，包含：S1、根据初始目标位置与尺度及跟踪器预装参数进行跟踪其初始化；S2、求解全局运动参数，对上一帧目标中心位置进行全局运动补偿，作为搜索窗口的中心位置；S3、提取位置样本和尺度样本；压缩位置样本和尺度样本；分别计算位置滤波频域响应图和尺度滤波频域响应图；计算位置滤波空域响应图和尺度滤波空域响应图，进行位置估计和尺度估计；自适应更新尺度增量因子；S4、提取位置样本和尺度样本；更新位置模型目标模板和尺度模型目标模板；构建位置特征投影矩阵和尺度特征投影矩阵；更新位置模型和尺度模型。本发明能确保特征点在背景上均匀分布，提高抵近飞行过程中目标跟踪的成功率。

Description

一种针对抵近飞行过程的目标跟踪方法

技术领域

本发明涉及计算机视觉、数字图像处理领域，特别涉及一种针对抵近飞行过程的目标跟踪方法。

背景技术

抵近飞行过程是指飞行器飞行并接近目标的过程，如精确制导武器的末制导过程。在抵近飞行过程中，目标跟踪的主要难点在于快速运动、运动模糊和尺度变化，具体说明如下：

(1)快速运动是指目标在相邻帧图像上的相对位移超出阈值，相对位移定义为目标在相邻帧图像上的位移与目标尺寸的比值。目标跟踪算法往往根据目标尺寸确定搜索范围，具体是将上一帧的目标区域扩大若干倍作为搜索区域。因此，不同于OTB(OnlineTracking Benchmark)属性标准采用的绝对位移，本发明采用目标在相邻帧图像上的相对位移作为快速运动的衡量指标。快速运动往往是由探测平台的抖动引起的。对此，现有目标跟踪算法一般通过在前端加入稳像处理模块，或者在跟踪过程中扩大搜索区域，来消除快速运动对算法跟踪性能的影响。然而，这两种方法均显著加大了算法的计算成本，限制了算法的应用范围。

(2)运动模糊是指由目标或飞行器的运动引起的目标成像模糊。成像模糊使得目标表观发生变化，在复杂场景中，将进一步降低目标与背景的区分度。传统的模板匹配算法，计算简单、易于实现，但已经无法应对实际的应用需求。因此，需要通过特征提取来增大目标与背景的区分度。近年来，梯度方向直方图(HOG，H,Histograms of OrientedGradients)特征及其改进特征在行人检测、人脸识别等领域获得了广泛的应用，其特征表达性能获得了业界的一直认可。

(3)尺度变化是指当前帧目标尺寸与第一帧目标尺寸的比值超出阈值。在抵近飞行的全过程中，目标尺度的变化分为两个阶段。第一阶段，飞行器距目标较远，目标尺度没有明显变化；第二阶段，飞行器距目标较近，目标尺度剧烈变化。现有的跟踪算法在对目标尺度的预测中，采用固定的尺度等级。这种方式，在抵近飞行的第一阶段，目标尺度基本不变，不仅浪费运算资源，还可能因干扰而错误地估计目标尺度；在第二阶段，目标尺度剧烈变化，极有可能超出尺度等级范围，无法准确估计目标尺度。目标跟踪算法若不能正确估计目标尺度，将无法准确维护目标模型，将会出现跟踪框漂移，甚至跟丢目标。

发明内容

本发明的目的在于提供一种针对抵近飞行过程的目标跟踪方法，能够确保特征点在背景上均匀分布，在图像尺寸已知的情况下，几乎没有时间消耗；本发明还提高了算法效率，并提高尺度估计的准确率，从而提高抵近飞行过程中目标跟踪的成功率。

为了达到上述目的，本发明通过以下技术方案实现：

一种针对抵近飞行过程的目标跟踪方法，包含以下步骤：

S1、跟踪器初始化：根据初始目标位置与尺度及跟踪器预装参数进行跟踪器初始化；

S2、确定搜索窗口中心位置：将上一帧提取均匀特征点作为局部特征，在当前帧提取匹配特征点，求解全局运动参数，对上一帧目标中心位置进行全局运动补偿，作为搜索窗口的中心位置；

步骤S3、位置滤波和尺度滤波：根据搜索窗口中心位置和上一帧目标尺度，提取位置样本和尺度样本；压缩位置样本和尺度样本；分别计算位置滤波频域响应图和尺度滤波频域响应图；计算位置滤波空域响应图和尺度滤波空域响应图，并进行位置估计和尺度估计；自适应更新尺度增量因子；

步骤S4、更新位置模型和尺度模型：根据当前帧目标位置和尺度，采用与滤波模块一样的特征提取方法，提取位置样本和尺度样本；更新位置模型目标模板和尺度模型目标模板；分别构建位置特征投影矩阵和尺度特征投影矩阵；分别更新位置模型和尺度模型。

优选地，所述步骤S1中，进一步包含：

步骤S1.1、输入初始目标位置p₀、初始目标尺寸sz₀、初始目标尺寸缩放因子κ_resize、目标区域扩展因子κ_padding、位置滤波器预期响应标准差因子κ_σ,trans、尺度等级n_S、插值后尺度等级n_SI、尺度滤波器预期响应标准差因子κ_σ,scale与尺度增量因子α；

步骤S1.2、初始化当前目标位置p_t和当前目标尺度因子s_t：

步骤S1.3、初始化目标基尺寸sz_base和跟踪窗口基尺寸sz_winbase：

步骤S1.4、初始化位置滤波预期响应：根据式(3)计算位置滤波预期响应标准σ_trans；生成位移矩阵T＝{(t_x,t_y)}；根据式(4)计算空域位置滤波预期响应g_trans(t_x,t_y)；根据式(5)并通过二维傅里叶变换获得频域位置滤波器预期响应G_trans；

σ_trans＝κ_σ,trans×prod(sz_base)/r_f (3)

其中，t_x表示x表示方向的位移，t_y表示y方向的位移；

步骤S1.5、初始化位置滤波汉宁窗：

其中，h和w分别表示位置滤波器的高度和宽度；

步骤S1.6、初始化尺度滤波预期响应：根据式(7)计算尺度滤波预期响应标准差σ_scale；根据式(8)生成尺度因子序列s和s_I；根据式(9)计算空域尺度滤波预期响应g_scale(s_i)；根据式(10)并通过二维傅里叶变换获得频域尺度滤波器预期响应G_scale；

σ_scale＝κ_σ,scale×n_SI (7)

步骤S1.7、初始化尺度滤波汉宁窗：

步骤S1.8、输出当前目标位置p_t、当前目标尺度因子s_t、目标基尺寸sz_base、跟踪窗口基尺寸sz_winbase、位置滤波预期响应G_trans、位置滤波汉宁窗h_trans、尺度因子序列s、插值后的尺度因子序列s_I、尺度滤波预期响应G_scale和尺度滤波汉宁窗h_scale。

优选地，所述步骤S2中，进一步包含以下过程：

输入τ-1时刻的图像I(τ-1)，尺寸为[W,H]，采样步长为s，对输入图像进行均匀特征点提取的结果为：

{(x_i,y_i)|i＝0,1,…,(mn-1)}； (12)

其中，m＝floor((W+s-1)/s)；n＝floor((H+s-1)/s)；(x_i,y_i)为第i个特征点的坐标，如下：

其中，floor(·)为向下取整函数；mod(a,b)为a被b整除的余数。

优选地，所述步骤S2中，采用PyrLK光流法在当前帧提取匹配特征点，进一步包含以下过程：

(1)输入τ-1时刻图像I(τ-1)、τ时刻的图像I(τ)并输入图像I(τ-1)上的待跟踪点集{(x_i,y_i)|i＝0,1,…,(mn-1)}；

(2)以递归方式分别对输入图像I(τ-1)和图像I(τ)构建图像金字塔，如下：

(2-a)输入图像I，尺寸为W×H；所述输入图像I包含图像I(τ-1)和图像I(τ)；

(2-b)令层数l＝0，金字塔的第l层图像I_(l)＝I，尺寸W_(l)×H_(l)＝W×H；

(2-c)令l＝l+1，若l＜L，转到步骤(2-d)，否则，转到步骤(2-e)；L表示金字塔的层数；

(2-d)根据式(14)，由金字塔的第l-1层图像I_(l-1)，计算第l层图像I_(l)，然后转到步骤(2-c)，循环进行；

(2-e)输出图像集{I_(l)},l＝0,1,…,L-1；

(3)对于待跟踪点集中的每一个点，逐层求解，计算金字塔光流，如下：

(3-a)输入图像I(τ-1)、I(τ)的对应的图像金子塔Pyr(τ-1)＝{I_(l)(τ-1)}和Pyr(τ)＝{I_(l)(τ)}，l＝0,1,…,L-1；

(3-b)对任一点p∈{(x_i,y_i)|i＝0,1,…,(mn-1)}，设其坐标为(x,y)，按照以下步骤，计算其光流矢量；

(3-c)点p在金字塔Pyr(τ-1)的第(l-1)层图像上的坐标为：

(3-d)由上至下逐层求解，以上一层的光流估计结果作为下一层光流估计的起点，获得每一层的光流矢量v_(l)；

(3-e)最终根据式(16)求得点p的光流矢量：

优选地，所述步骤S2中，基于最小中值平方算法求解全局运动参数，进一步包含以下过程：

采用八参数透视模型对全局运动建模，即：

其中，(x_i,y_i)和(x′_i,y′_i)分别表示图像上的点在全局运动补偿前、后的坐标位置；s_i为尺度因子；a₀,…,a₇分别为全局运动参数；P称作单应性矩阵；

所述最小中值平方算法的计算步骤如下：

(a)输入匹配特征点对G₁＝{(x_i,y_i),(x′_i,y′_i)},i＝1,…,N，N表示匹配特征点对的数量；令k＝0；

(b)k←k+1；从已获得的N对匹配特征点对中随机选4对，求解线性方程组，获得单应性矩阵

(c)利用单应性矩阵

将图像I(τ)上的所有特征点(x_i,y_i)映射到图像I(τ+1)中，得到特征点位置的估计值

(d)求特征点位置估计残差的中值

(e)若k＜K，则重复上述步骤(b)、(c)、(d)，由此获得K组数据：

其中，K表示参数估计的次数，为预先设定的常数；

(f)求使得残差中值最小的矩阵估计：

优选地，参数K的确定过程包含：

设真实数据模型的错误率为ε，则一次估计中所取的四个数据点均为内点的概率为(1-ε)⁴；

要保证K次估计中至少有一次估计所取的数据点均为内点的概率p较大，则K需要满足：

确定p和ε值后，并有式(23)计算得出参数K的值。

优选地，所述步骤S2中，所述全局运动补偿进一步包含：

利用获得的全局运动参数，根据式(24)对(τ-1)时刻的目标位置p_τ-1＝(x_τ-1,y_τ-1)进行全局运动补偿，获得补偿后的位置p′_τ-1＝(x′_τ-1,y′_τ-1)，作为搜索窗口的中心位置：

优选地，所述步骤S3中，提取位置样本和尺度样本的方法包含：

步骤S31.1、提取位置滤波窗口和尺度滤波窗口：

设经全局运动补偿后，上一帧的目标位置为p′_t-1、目标尺度为s_t-1；

跟踪其初始化模块中确定的目标基尺寸sz_base、跟踪窗口基尺寸sz_winbase、尺度等级n_S、尺度因子序列s；

预装参数位置模型尺寸sz_trans和尺度模型尺寸sz_scale；

以上一帧的目标位置p′_t-1为中心，提取大小为sz_winbase×s_t-1的图像块，即为位置滤波窗口；

以上一帧的目标位置p′_t-1为中心，提取大小为sz_base×s_t-1×s_i,s_i∈s的图像块集，即为尺度滤波窗口，共n_S个窗口；

步骤S31.2、窗口归整化：采用双线性插值法，将位置滤波窗口缩放到位置模型尺寸sz_trans；将尺度滤波器窗口缩放到尺度模型尺寸sz_scale；

步骤S31.3、特征提取：对缩放后的位置滤波窗口进行FHOG特征提取，记作z_t,trans；分别对缩放后的n_S个尺度滤波窗口进行FHOG提取，将获得的n_S个FHOG特征压平，并在尺度方向进行特征拼接，获得多尺度FHOG特征，记作z_t,scale；

和/或，所述步骤S3中，压缩位置样本和尺度样本的方法包含：

步骤S32.1、设上一帧模型更新阶段获得的位置特征投影矩阵为P_t-1,trans，尺度特征投影矩阵为

步骤S32.2、根据式(25)压缩位置样本，根据式(26)压缩尺度样本：

和/或，所述步骤S3中，计算位置滤波频域响应图和尺度滤波频域响应图的方法包括：

步骤S33.1、根据式(27)～(28)分别将压缩后的位置样本和尺度样本变换到频域；其中，

和

分别表示二维傅里叶变换和一维傅里叶变换：

步骤S33.2、设上一帧模型更新阶段获得的位置滤波器的分子和分母分别为

和

尺度滤波器的分子和分母分别为

和

预装参数归整化权重为λ；

表示位置样本压缩后的维度，压缩前维度为d；根据式(29)～(30)分别计算位置滤波频域响应图

和尺度滤波频域响应图

和/或，所述步骤S3中，计算位置滤波空域响应图和尺度滤波空域响应图的方法如下：

步骤S34.1、频率响应图升采样：根据频域特性将位置滤波频域响应图

升采样到位置模型尺寸，记作Y_t,trans；将尺度滤波频域响应图

升采样到插值后尺度等级n_SI尺寸，记作Y_t,scale；

步骤S34.2、根据式(31)～(32)分别将插值后的位置滤波频域响应图Y_t,trans和尺度滤波频域响应图Y_t,scale变换到空域；其中，

和

分别表示二维傅里叶反变换和一维傅里叶反变换：

和/或，所述步骤S3中，位置估计和尺度估计的方法包含：

步骤S35.1、位置估计：根据式(33)，通过最大化位置滤波空域响应y_t,trans，估计位移矢量

其中，T＝{(t_x,t_y)}为位移矩阵；根据式(34)，获得当前帧目标位置估计p_t：

步骤S35.2、尺度估计：根据式(35)，通过最大化尺度滤波空域响应y_t,scale，获得当前帧目标尺度因子s；根据式(36)，获得当前帧目标尺度估计s_t：

s_t＝s_t-1×s (36)

步骤S35.3、输出当前帧目标位置p_t和当前帧尺寸sz_t＝floor(sz_base×s_t)；其中，floor(·)表示向下取整；

和/或，所述步骤S3中，自适应更新尺度增量因子的方法包含：

若连续N帧满足式(37)：

lg(s_t)/lg(α)≥(n_SI-3)， (37)

令α←α+0.03，p₀＝p_t，s₀＝s_t，并更新尺度因子序列s、插值后的尺度因子序列s_I、尺度滤波预期响应G_scale和尺度滤波汉宁窗h_scale。

优选地，所述步骤S4中，更新位置模型目标模板和尺度模型目标目标模板的方法包含：

步骤S41.1、设提取的位置样本和尺度样本分别为f_t,trans和f_t,scale；

步骤S41.2、根据式(38)～(39)，分别更新位置模型目标模板u_t,trans和尺度模型目标模板u_t,scale：

其中，η表示模型学习率；

和/或，所述步骤S4中，构建位置特征投影矩阵和尺度特征投影矩阵的方法包含：

步骤S42.1、构建位置特征投影矩阵：对位置模型目标模板u_t,trans的自相关矩阵进行奇异值分解，取前d_trans个最大的特征值对应的特征向量，构成位置特征投影矩阵P_t,trans；其中，d_trans为算法预装参数，表示压缩后位置特征的维度；

步骤S42.2、构建尺度特征投影矩阵：分别对尺度模型目标模板u_t,scale和尺度样本f_t,scale进行奇异值分解，对应地获得尺度特征投影矩阵

和

和/或，所述步骤S4中，更新位置模型和尺度模型的方法包含：

步骤S43.1、更新位置模型：根据式(40)分别压缩目标模板和位置样本；根据式(41)将其变换到频域；根据式(42)更新位置模型；其中，

为位置模型的分子，

为位置模型的分子，η为模型学习率：

其中，

表示位置滤波预期响应G_trans的共轭，

表示

的共轭；

步骤S43.2、更新尺度模型：根据式(43)分别压缩目标模板和尺度样本；根据式(44)将其变换到频域；根据式(45)更新尺度模型，其中，

为位置模型的分子，

为位置模型的分子，η为模型学习率：

其中，

表示位置滤波预期响应G_scale的共轭，

表示

的共轭。

与现有技术相比，本发明的有益效果在于：

(1)本发明利用金字塔光流法(PyrLK)和最小中值平方法进行相邻帧间的全局运动估计，从而对前一帧位置进行全局运动补偿；本发明提取均匀特征点作为金字塔光流法的输入：一方面，能够确保特征点在背景上均匀分布；另一方面，在图像尺寸已知的情况下，几乎没有时间消耗；

(2)本发明通过对前一帧位置进行全局运动补偿，确定搜索区域的中心位置。一方面，不需要在前端增加稳像模块，就可以缩小位置滤波的搜索区域，显著减小了算法的运算量；另一方面，减小了搜索区域中心位置与当前帧目标位置的距离，使得位置滤波与尺度滤波的并行计算成为可能，显著提高了算法效率；

(3)本发明在尺度滤波中，采用自适应尺度增量因子，提高了尺度估计的准确率，从而增大了目标跟踪的成功率。

附图说明

图1为本发明的针对抵近飞行过程的目标跟踪算法流程图；

图2为本发明的金字塔跟踪图解。

具体实施方式

为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

如图1-图2所示，本发明提供了一种针对抵近飞行过程的目标跟踪方法，包括以下步骤：

步骤S1、跟踪器初始化：根据初始目标位置和尺度及跟踪器预装参数进行跟踪器初始化。

示例地，所述步骤S1中，进一步包含：如图1所示，当时间t＝1时，开始跟踪，记外部输入的目标初始位置和尺寸分别为p₀和sz₀，此时要根据初始目标位置和尺度及跟踪器预装参数进行跟踪其初始化；具体如下：

步骤S1.1、输入初始目标位置p₀、初始目标尺寸sz₀、初始目标尺寸缩放因子κ_resize、目标区域扩展因子κ_padding、位置滤波器预期响应标准差因子κ_σ,trans、尺度等级n_S、插值后尺度等级n_SI、尺度滤波器预期响应标准差因子κ_σ,scale、尺度增量因子α。

步骤S1.2、根据下式(1)初始化当前目标位置p_t和当前目标尺度因子s_t：

步骤S1.3、根据下式(2)初始化目标基尺寸sz_base和跟踪窗口基尺寸sz_winbase：

步骤S1.4、初始化位置滤波预期响应：首先，根据式(3)计算位置滤波预期响应标准σ_trans；然后，生成位移矩阵T＝{(t_x,t_y)}；接着，根据式(4)计算空域位置滤波预期响应g_trans(t_x,t_y)；最后，根据式(5)，通过二维傅里叶变换获得频域位置滤波器预期响应G_trans；如下：

σ_trans＝κ_σ,trans×prod(sz_base)/r_f (3)

其中，t_x表示x表示方向的位移，t_y表示y方向的位移。

步骤S1.5、根据下式(6)初始化位置滤波汉宁窗：

其中，h和w分别表示位置滤波器的高度和宽度，即模型尺寸。

步骤S1.6、初始化尺度滤波预期响应：首先，根据式(7)计算尺度滤波预期响应标准差σ_scale；然后，根据式(8)生成尺度因子序列s和s_I；接着，根据式(9)计算空域尺度滤波预期响应g_scale(s_i)；最后，根据式(10)，通过二维傅里叶变换获得频域尺度滤波器预期响应G_scale；如下：

σ_scale＝κ_σ,scale×n_SI (7)

步骤S1.7、根据下式(11)初始化尺度滤波汉宁窗：

步骤S1.8、输出当前目标位置p_t、当前目标尺度因子s_t、目标基尺寸sz_base、跟踪窗口基尺寸sz_winbase、位置滤波预期响应G_trans、位置滤波汉宁窗h_trans、尺度因子序列s、插值后的尺度因子序列s_I、尺度滤波预期响应G_scale和尺度滤波汉宁窗h_scale。

步骤S2、确定搜索窗口中心位置：首先，在上一帧提取均匀特征点作为局部特征；然后，采用PyrLK光流法在当前帧提取匹配特征点；接着，基于LMedS算法求解全局运动参数；最后，对上一帧目标中心位置进行全局运动补偿，作为搜索窗口的中心位置。

示例地，所述步骤S2中进一步包含：

步骤S2.1、提取均匀特征点：输入τ-1时刻的图像I(τ-1)，尺寸为[W,H]，采样步长为s，则对输入图像进行均匀特征点提取的结果为：

{(x_i,y_i)|i＝0,1,…,(mn-1)}； (12)

其中，m＝floor((W+s-1)/s)，n＝floor((H+s-1)/s)，W和H分别表示图像的宽度和高度；(x_i,y_i)为第i个特征点的坐标，如下：

其中，floor(·)为向下取整函数；mod(a,b)为a被b整除的余数。

步骤S2.2、基于PyrLK光流法匹配特征点；

如图2所示，输入前一帧图像、当前帧图像以及前一帧图像上的待跟踪点集，金字塔跟踪过程如下：首先，以递归方式自底向上对前一帧图像和当前帧图像构建图像金字塔；然后，对于待跟踪点集中的每一个点，自顶向下逐层求解，获得每一层的光流矢量；最后根据金字塔各层图像之间的比例关系，合成这些光流矢量，获得光流矢量，从而获得待跟踪点的匹配特征点。

所述步骤S2.2进一步包含：

(1)输入τ-1时刻的图像I(τ-1)、τ时刻(当前时刻)的图像I(τ)和图像I(τ-1)上的待跟踪点集{(x_i,y_i)|i＝0,1,…,(mn-1)}；

(2)首先，以递归方式分别对输入图像I(τ-1)和图像I(τ)构建图像金字塔，如下：

(2-a)输入图像I(例如I(τ-1)或I(τ))，尺寸为W×H；

(2-b)令层数l＝0，金字塔的第l层图像I_(l)＝I，尺寸W_(l)×H_(l)＝W×H；

(2-c)令l＝l+1，若l＜L，转到步骤(2-d)，否则，转到步骤(2-e)；其中，L表示金字塔的层数，层数增加将会增加内存和时间成本，一般取3-5层；

(2-d)根据式(14)，由金字塔的第l-1层图像I_(l-1)，计算第l层图像I_(l)，然后转到步骤(2-c)，循环进行；

(2-e)输出图像集{I_(l)},l＝0,1,…,L-1。

(3)然后，对于待跟踪点集中的每一个点，逐层求解，计算金字塔光流，如下：

(3-a)输入图像I(τ-1)、I(τ)的图像金子塔Pyr(τ-1)＝{I_(l)(τ-1)}和Pyr(τ)＝{I_(l)(τ)}，l＝0,1,…,L-1；

(3-b)对任一点p∈{(x_i,y_i)|i＝0,1,…,(mn-1)}，设其坐标为(x,y)，按照以下步骤，计算其光流矢量；

(3-c)点p在金字塔Pyr(τ-1)的第(l-1)层图像上的坐标为：

(3-d)由上至下(即自顶向下)逐层求解，以上一层的光流估计结果作为下一层光流估计的起点，获得每一层的光流矢量v_(l)，如图2所示；其中，图像金字塔的构是由高分辨率到低分辨率以递归的方式构建的，图像尺寸由大到小，即由金字塔的底部向金字塔的顶部方向依次构建。而光流的求解过程是由金字塔的顶部向金字塔的底部依次求解，最终得到在原始图像尺寸上的光流矢量。

(3-e)最终，根据式(16)求得点p的光流矢量：

步骤S2.3、基于LMedS算法求解全局运动参数：采用八参数透视模型对全局运动建模，即：

其中，(x_i,y_i)和(x′_i,y′_i)分别表示图像上的点在全局运动补偿前、后的坐标位置；s_i为尺度因子；a₀,…,a₇分别为全局运动参数；P也称作单应性矩阵；本发明采用抗噪性能较好的最小中值平方(Least median of squares,LMedS)算法进行全局运动参数估计，该算法的步骤如下：

(a)输入匹配特征点对G₁＝{(x_i,y_i),(x′_i,y′_i)},i＝1,…,N，N表示匹配特征点对的数量；k为用于计数的参数，令k＝0；

(b)k←k+1；从已获得的N对匹配特征点对中随机选4对(任意三点不共线)，求解线性方程组，获得单应性矩阵

(c)利用矩阵

将图像I(τ)上的所有特征点(x_i,y_i)映射到图像I(τ+1)中，得到特征点位置的估计值

即：

(d)求特征点位置估计残差的中值

即：

(e)若k＜K，则重复步骤(b)、(c)、(d)，由此获得K组数据：

(f)求使得残差中值最小的矩阵估计：

(g)采用LM(Levenberg-Marquardt)算法，利用所有内点精确优化模型参数。其中，K表示参数估计的次数，为预先设定的常数，大于步骤(g)中算出的值。下面为确定参数K的过程：设真实数据模型的错误率(即错匹配点占特征点集的比例)为ε，则一次估计中所取的四个数据点均为内点的概率为(1-ε)⁴；要保证“K次估计中至少有一次估计所取的数据点均为内点”的概率p较大，则K需要满足：

令p＝0.995,ε＝0.45，由式(23)计算得K≈56；因此，即使匹配点对中有45％的错匹配，也仅需进行56次估计就可以获得全局最优的参数估计。

步骤S2.4、全局运动补偿：利用前面获得的全局运动参数，根据式(24)对(τ-1)时刻的目标位置p_τ-1＝(x_τ-1,y_τ-1)进行全局运动补偿，获得补偿后的位置p′_τ-1＝(x′_τ-1,y′_τ-1)，作为搜索窗口的中心位置：

步骤S3、位置滤波和尺度滤波：首先，根据搜索窗口中心位置和上一帧目标尺度，提取位置样本和尺度样本；然后，压缩位置样本和尺度样本；接着，分别计算位置滤波频域响应图和尺度滤波频域响应图；再计算位置滤波空域响应图和尺度滤波空域响应图，并进行位置估计和尺度估计；最后，自适应更新尺度增量因子。

所述步骤S3中，所述提取位置样本和尺度样本的具体方法包括：

(31.1)提取位置滤波窗口和尺度滤波窗口：

设经全局运动补偿后，上一帧的目标位置为p′_t-1、目标尺度为s_t-1；跟踪其初始化模块中确定的目标基尺寸sz_base、跟踪窗口基尺寸sz_winbase、尺度等级n_S、尺度因子序列s；预装参数位置模型尺寸sz_trans和尺度模型尺寸sz_scale；

以p′_t-1为中心，提取大小为sz_winbase×s_t-1的图像块，即为位置滤波窗口；

以p′_t-1为中心，提取大小为sz_base×s_t-1×s_i,s_i∈s的图像块集，即为尺度滤波窗口，共n_S个窗口。

(31.2)窗口归整化：采用双线性插值法，将位置滤波窗口缩放到位置模型尺寸sz_trans；将尺度滤波器窗口缩放到尺度模型尺寸sz_scale。

(31.3)特征提取：对缩放后的位置滤波窗口进行FHOG特征提取，记作z_t,trans；分别对缩放后的n_S个尺度滤波窗口进行FHOG提取，将获得的n_S个FHOG特征压平，并在尺度方向进行特征拼接，获得多尺度FHOG特征，记作z_t,scale。

所述步骤S3中，所述压缩位置样本和尺度样本的具体方法包括：

(32.1)设上一帧模型更新阶段获得的位置特征投影矩阵为P_t-1,trans，尺度特征投影矩阵为

(32.2)根据式(25)压缩位置样本，根据式(26)压缩尺度样本：

所述步骤S3中，计算位置滤波频域响应图和尺度滤波频域响应图的具体方法包括：

(33.1)根据式(27)～(28)，分别将压缩后的位置样本和尺度样本变换到频域；其中，

和

分别表示二维傅里叶变换和一维傅里叶变换：

(33.2)设上一帧模型更新阶段获得的位置滤波器的分子和分母分别为

和

尺度滤波器的分子和分母分别为

和

预装参数归整化权重为λ；

表示位置样本压缩后的维度，压缩前维度为d。根据式(29)～(30)，分别计算位置滤波频域响应图

和尺度滤波频域响应图

所述步骤S3中，所述计算位置滤波空域响应图和尺度滤波空域响应图的具体方法如下：

(34.1)频率响应图升采样：根据频域特性，将位置滤波频域响应图

升采样到位置模型尺寸，记作Y_t,trans；将尺度滤波频域响应图

升采样到插值后尺度等级n_SI尺寸，记作Y_t,scale。

(34.2)根据式(31)～(32)，分别将插值后的位置滤波频域响应图Y_t,trans和尺度滤波频域响应图Y_t,scale变换到空域；其中，

和

分别表示二维傅里叶反变换和一维傅里叶反变换：

所述步骤S3中，所述位置估计和尺度估计的具体方法如下：

(35.1)位置估计：首先，根据式(33)，通过最大化位置滤波空域响应y_t,trans，估计位移矢量

其中，T＝{(t_x,t_y)}为位移矩阵；然后，根据式(34)，获得当前帧目标位置估计p_t：

(35.2)尺度估计：根据式(35)，通过最大化尺度滤波空域响应y_t,scale，获得当前帧目标尺度因子s；然后，根据式(36)，获得当前帧目标尺度估计s_t：

s_t＝s_t-1×s (36)

(35.3)输出当前帧目标位置和尺寸：当前帧目标位置为p_t，当前帧目标尺寸为sz_t＝floor(sz_base×s_t)；其中，floor(·)表示向下取整。

所述步骤S3中，所述自适应更新尺度增量因子的具体方法如下：

(36.1)若连续N帧满足式(37)：

lg(s_t)/lg(α)≥(n_SI-3)， (37)

则令α←α+0.03，p₀＝p_t，s₀＝s_t，并更新尺度因子序列s，插值后的尺度因子序列s_I，尺度滤波预期响应G_scale，尺度滤波汉宁窗h_scale。

步骤S4、更新位置模型和尺度模型：首先，根据当前帧目标位置和尺度，采用与滤波模块一样的特征提取方法，提取位置样本和尺度样本；然后，更新位置模型目标模板和尺度模型目标模板；接着，分别构建位置特征投影矩阵和尺度特征投影矩阵；最后，分别更新位置模型和尺度模型。

所述步骤S4中，所述更新位置模型目标模板和尺度模型目标目标模板的具体方法如下：

(41.1)设提取的位置样本和尺度样本分别为f_t,trans和f_t,scale。

(41.2)根据式(38)～(39)，分别更新位置模型目标模板u_t,trans和尺度模型目标模板u_t,scale：

其中，η表示模型学习率。

所述步骤S4中，所述构建位置特征投影矩阵和尺度特征投影矩阵的具体方法如下：

(42.1)构建位置特征投影矩阵：对位置模型目标模板u_t,trans的自相关矩阵进行奇异值分解，取前d_trans个最大的特征值对应的特征向量，构成位置特征投影矩阵P_t,trans；其中，d_trans为算法预装参数，表示压缩后位置特征的维度。

(42.2)构建尺度特征投影矩阵：分别对尺度模型目标模板u_t,scale和尺度样本f_t,scale进行奇异值分解，获得尺度特征投影矩阵

和

所述步骤S4中，所述更新位置模型和尺度模型的具体方法如下：

(43.1)更新位置模型：首先，根据式(40)分别压缩目标模板和位置样本；然后，根据式(41)将其变换到频域；最后，根据式(42)更新位置模型，其中，

为位置模型的分子，

为位置模型的分子，η为模型学习率：

其中，

表示位置滤波预期响应G_trans的共轭，

表示

的共轭。

(43.2)更新尺度模型：首先，根据式(43)分别压缩目标模板和尺度样本；然后，根据式(44)将其变换到频域；最后，根据式(45)更新尺度模型，其中，

为位置模型的分子，

为位置模型的分子，η为模型学习率：

其中，

表示位置滤波预期响应G_scale的共轭，

表示

的共轭。

综上所述，本发明利用金字塔光流法(PyrLK)和最小中值平方法进行相邻帧间的全局运动估计，从而对前一帧位置进行全局运动补偿；本发明提取均匀特征点作为金字塔光流法的输入：一方面，能够确保特征点在背景上均匀分布；另一方面，在图像尺寸已知的情况下，几乎没有时间消耗；本发明通过对前一帧位置进行全局运动补偿，确定搜索区域的中心位置。一方面，不需要在前端增加稳像模块，就可以缩小位置滤波的搜索区域，显著减小了算法的运算量；另一方面，减小了搜索区域中心位置与当前帧目标位置的距离，使得位置滤波与尺度滤波的并行计算成为可能，显著提高了算法效率；本发明在尺度滤波中，采用自适应尺度增量因子，提高了尺度估计的准确率，从而增大了目标跟踪的成功率。

尽管本发明的内容已经通过上述优选实施例作了详细介绍，但应当认识到上述的描述不应被认为是对本发明的限制。在本领域技术人员阅读了上述内容后，对于本发明的多种修改和替代都将是显而易见的。因此，本发明的保护范围应由所附的权利要求来限定。

Claims (9)

1.一种针对抵近飞行过程的目标跟踪方法，其特征在于，包含以下步骤：

S1、跟踪器初始化：根据初始目标位置与尺度及跟踪器预装参数进行跟踪器初始化；

S2、确定搜索窗口中心位置：将上一帧提取均匀特征点作为局部特征，在当前帧提取匹配特征点，求解全局运动参数，对上一帧目标中心位置进行全局运动补偿，作为搜索窗口的中心位置；

步骤S3、位置滤波和尺度滤波：根据搜索窗口中心位置和上一帧目标尺度，提取位置样本和尺度样本；压缩位置样本和尺度样本；分别计算位置滤波频域响应图和尺度滤波频域响应图；计算位置滤波空域响应图和尺度滤波空域响应图，并进行位置估计和尺度估计；自适应更新尺度增量因子；

步骤S4、更新位置模型和尺度模型：根据当前帧目标位置和尺度，采用与滤波模块一样的特征提取方法，提取位置样本和尺度样本；更新位置模型目标模板和尺度模型目标模板；分别构建位置特征投影矩阵和尺度特征投影矩阵；分别更新位置模型和尺度模型。

2.如权利要求1所述的针对抵近飞行过程的目标跟踪方法，其特征在于，所述步骤S1中，进一步包含：

步骤S1.1、输入初始目标位置p₀、初始目标尺寸sz₀、初始目标尺寸缩放因子κ_resize、目标区域扩展因子κ_padding、位置滤波器预期响应标准差因子κ_σ,trans、尺度等级n_S、插值后尺度等级n_SI、尺度滤波器预期响应标准差因子κ_σ,scale与尺度增量因子α；

步骤S1.2、初始化当前目标位置p_t和当前目标尺度因子s_t：

步骤S1.3、初始化目标基尺寸sz_base和跟踪窗口基尺寸sz_winbase：

步骤S1.4、初始化位置滤波预期响应：根据式(3)计算位置滤波预期响应标准σ_trans；生成位移矩阵T＝{(t_x,t_y)}；根据式(4)计算空域位置滤波预期响应g_trans(t_x,t_y)；根据式(5)并通过二维傅里叶变换获得频域位置滤波器预期响应G_trans；

σ_trans＝κ_σ,trans×prod(sz_base)/r_f (3)

其中，t_x表示x表示方向的位移，t_y表示y方向的位移；

步骤S1.5、初始化位置滤波汉宁窗：

其中，h和w分别表示位置滤波器的高度和宽度；

步骤S1.6、初始化尺度滤波预期响应：根据式(7)计算尺度滤波预期响应标准差σ_scale；根据式(8)生成尺度因子序列s和s_I；根据式(9)计算空域尺度滤波预期响应g_scale(s_i)；根据式(10)并通过二维傅里叶变换获得频域尺度滤波器预期响应G_scale；

σ_scale＝κ_σ,scale×n_SI (7)

步骤S1.7、初始化尺度滤波汉宁窗：

步骤S1.8、输出当前目标位置p_t、当前目标尺度因子s_t、目标基尺寸sz_base、跟踪窗口基尺寸sz_winbase、位置滤波预期响应G_trans、位置滤波汉宁窗h_trans、尺度因子序列s、插值后的尺度因子序列s_I、尺度滤波预期响应G_scale和尺度滤波汉宁窗h_scale。

3.如权利要求2所述的针对抵近飞行过程的目标跟踪方法，其特征在于，所述步骤S2中，进一步包含以下过程：

输入τ-1时刻的图像I(τ-1)，尺寸为[W,H]，采样步长为s，对输入图像进行均匀特征点提取的结果为：

{(x_i,y_i)|i＝0,1,…,(mn-1)}； (12)

其中，m＝floor((W+s-1)/s)；n＝floor((H+s-1)/s)；(x_i,y_i)为第i个特征点的坐标，如下：

其中，floor(·)为向下取整函数；mod(a,b)为a被b整除的余数。

4.如权利要求3所述的针对抵近飞行过程的目标跟踪方法，其特征在于，所述步骤S2中，采用PyrLK光流法在当前帧提取匹配特征点，进一步包含以下过程：

(1)输入τ-1时刻图像I(τ-1)、τ时刻的图像I(τ)并输入图像I(τ-1)上的待跟踪点集{(x_i,y_i)|i＝0,1,…,(mn-1)}；

(2)以递归方式分别对输入图像I(τ-1)和图像I(τ)构建图像金字塔，如下：

(2-a)输入图像I，尺寸为W×H；所述输入图像I包含图像I(τ-1)和图像I(τ)；

(2-b)令层数l＝0，金字塔的第l层图像I_(l)＝I，尺寸W_(l)×H_(l)＝W×H；

(2-c)令l＝l+1，若l＜L，转到步骤(2-d)，否则，转到步骤(2-e)；L表示金字塔的层数；

(2-d)根据式(14)，由金字塔的第l-1层图像I_(l-1)，计算第l层图像I_(l)，然后转到步骤(2-c)，循环进行；

(2-e)输出图像集{I_(l)},l＝0,1,…,L-1；

(3)对于待跟踪点集中的每一个点，逐层求解，计算金字塔光流，如下：

(3-a)输入图像I(τ-1)、I(τ)的对应的图像金子塔Pyr(τ-1)＝{I_(l)(τ-1)}和Pyr(τ)＝{I_(l)(τ)}，l＝0,1,…,L-1；

(3-b)对任一点p∈{(x_i,y_i)|i＝0,1,…,(mn-1)}，设其坐标为(x,y)，按照以下步骤，计算其光流矢量；

(3-c)点p在金字塔Pyr(τ-1)的第(l-1)层图像上的坐标为：

(3-d)由上至下逐层求解，以上一层的光流估计结果作为下一层光流估计的起点，获得每一层的光流矢量v_(l)；

(3-e)最终根据式(16)求得点p的光流矢量：

5.如权利要求4所述的针对抵近飞行过程的目标跟踪方法，其特征在于，所述步骤S2中，基于最小中值平方算法求解全局运动参数，进一步包含以下过程：

采用八参数透视模型对全局运动建模，即：

其中，(x_i,y_i)和(x′_i,y′_i)分别表示图像上的点在全局运动补偿前、后的坐标位置；s_i为尺度因子；a₀,…,a₇分别为全局运动参数；P称作单应性矩阵；

所述最小中值平方算法的计算步骤如下：

(a)输入匹配特征点对G₁＝{(x_i,y_i),(x′_i,y′_i)},i＝1,…,N，N表示匹配特征点对的数量；令k＝0；

(b)k←k+1；从已获得的N对匹配特征点对中随机选4对，求解线性方程组，获得单应性矩阵

(c)利用单应性矩阵

将图像I(τ)上的所有特征点(x_i,y_i)映射到图像I(τ+1)中，得到特征点位置的估计值

(d)求特征点位置估计残差的中值

(e)若k＜K，则重复上述步骤(b)、(c)、(d)，由此获得K组数据：

其中，K表示参数估计的次数，为预先设定的常数；

(f)求使得残差中值最小的矩阵估计：

6.如权利要求5所述的针对抵近飞行过程的目标跟踪方法，其特征在于，参数K的确定过程包含：

设真实数据模型的错误率为ε，则一次估计中所取的四个数据点均为内点的概率为(1-ε)⁴；

要保证K次估计中至少有一次估计所取的数据点均为内点的概率p较大，则K需要满足：

确定p和ε值后，并有式(23)计算得出参数K的值。

7.如权利要求6所述的针对抵近飞行过程的目标跟踪方法，其特征在于，所述步骤S2中，所述全局运动补偿进一步包含：

利用获得的全局运动参数，根据式(24)对(τ-1)时刻的目标位置p_τ-1＝(x_τ-1,y_τ-1)进行全局运动补偿，获得补偿后的位置p′_τ-1＝(x′_τ-1,y′_τ-1)，作为搜索窗口的中心位置：

8.如权利要求7所述的针对抵近飞行过程的目标跟踪方法，其特征在于，所述步骤S3中，提取位置样本和尺度样本的方法包含：

步骤S31.1、提取位置滤波窗口和尺度滤波窗口：

设经全局运动补偿后，上一帧的目标位置为p′_t-1、目标尺度为s_t-1；

跟踪其初始化模块中确定的目标基尺寸sz_base、跟踪窗口基尺寸sz_winbase、尺度等级n_S、尺度因子序列s；

预装参数位置模型尺寸sz_trans和尺度模型尺寸sz_scale；

以上一帧的目标位置p′_t-1为中心，提取大小为sz_winbase×s_t-1的图像块，即为位置滤波窗口；

以上一帧的目标位置p′_t-1为中心，提取大小为sz_base×s_t-1×s_i,s_i∈s的图像块集，即为尺度滤波窗口，共n_S个窗口；

步骤S31.2、窗口归整化：采用双线性插值法，将位置滤波窗口缩放到位置模型尺寸sz_trans；将尺度滤波器窗口缩放到尺度模型尺寸sz_scale；

步骤S31.3、特征提取：对缩放后的位置滤波窗口进行FHOG特征提取，记作z_t,trans；分别对缩放后的n_S个尺度滤波窗口进行FHOG提取，将获得的n_S个FHOG特征压平，并在尺度方向进行特征拼接，获得多尺度FHOG特征，记作z_t,scale；

和/或，所述步骤S3中，压缩位置样本和尺度样本的方法包含：

步骤S32.1、设上一帧模型更新阶段获得的位置特征投影矩阵为P_t-1,trans，尺度特征投影矩阵为

步骤S32.2、根据式(25)压缩位置样本，根据式(26)压缩尺度样本：

和/或，所述步骤S3中，计算位置滤波频域响应图和尺度滤波频域响应图的方法包括：

步骤S33.1、根据式(27)～(28)分别将压缩后的位置样本和尺度样本变换到频域；其中，

和

分别表示二维傅里叶变换和一维傅里叶变换：

步骤S33.2、设上一帧模型更新阶段获得的位置滤波器的分子和分母分别为

和

尺度滤波器的分子和分母分别为

和

预装参数归整化权重为λ；

表示位置样本压缩后的维度，压缩前维度为d；根据式(29)～(30)分别计算位置滤波频域响应图

和尺度滤波频域响应图

和/或，所述步骤S3中，计算位置滤波空域响应图和尺度滤波空域响应图的方法如下：

步骤S34.1、频率响应图升采样：根据频域特性将位置滤波频域响应图

升采样到位置模型尺寸，记作Y_t,trans；将尺度滤波频域响应图

升采样到插值后尺度等级n_SI尺寸，记作Y_t,scale；

步骤S34.2、根据式(31)～(32)分别将插值后的位置滤波频域响应图Y_t,trans和尺度滤波频域响应图Y_t,scale变换到空域；其中，

和

分别表示二维傅里叶反变换和一维傅里叶反变换：

和/或，所述步骤S3中，位置估计和尺度估计的方法包含：

步骤S35.1、位置估计：根据式(33)，通过最大化位置滤波空域响应y_t,trans，估计位移矢量

其中，T＝{(t_x,t_y)}为位移矩阵；根据式(34)，获得当前帧目标位置估计p_t：

步骤S35.2、尺度估计：根据式(35)，通过最大化尺度滤波空域响应y_t,scale，获得当前帧目标尺度因子s；根据式(36)，获得当前帧目标尺度估计s_t：

s_t＝s_t-1×s (36)

步骤S35.3、输出当前帧目标位置p_t和当前帧尺寸sz_t＝floor(sz_base×s_t)；其中，floor(·)表示向下取整；

和/或，所述步骤S3中，自适应更新尺度增量因子的方法包含：

若连续N帧满足式(37)：

lg(s_t)/lg(α)≥(n_SI-3)， (37)

令α←α+0.03，p₀＝p_t，s₀＝s_t，并更新尺度因子序列s、插值后的尺度因子序列s_I、尺度滤波预期响应G_scale和尺度滤波汉宁窗h_scale。

9.如权利要求8所述的针对抵近飞行过程的目标跟踪方法，其特征在于，所述步骤S4中，更新位置模型目标模板和尺度模型目标目标模板的方法包含：

步骤S41.1、设提取的位置样本和尺度样本分别为f_t,trans和f_t,scale；

步骤S41.2、根据式(38)～(39)，分别更新位置模型目标模板u_t,trans和尺度模型目标模板u_t,scale：

其中，η表示模型学习率；

和/或，所述步骤S4中，构建位置特征投影矩阵和尺度特征投影矩阵的方法包含：

步骤S42.1、构建位置特征投影矩阵：对位置模型目标模板u_t,trans的自相关矩阵进行奇异值分解，取前d_trans个最大的特征值对应的特征向量，构成位置特征投影矩阵P_t,trans；其中，d_trans为算法预装参数，表示压缩后位置特征的维度；

步骤S42.2、构建尺度特征投影矩阵：分别对尺度模型目标模板u_t,scale和尺度样本f_t,scale进行奇异值分解，对应地获得尺度特征投影矩阵

和

和/或，所述步骤S4中，更新位置模型和尺度模型的方法包含：

步骤S43.1、更新位置模型：根据式(40)分别压缩目标模板和位置样本；根据式(41)将其变换到频域；根据式(42)更新位置模型；其中，

为位置模型的分子，

为位置模型的分子，η为模型学习率：

其中，

表示位置滤波预期响应G_trans的共轭，

表示

的共轭；步骤S43.2、更新尺度模型：根据式(43)分别压缩目标模板和尺度样本；根据式(44)将其变换到频域；根据式(45)更新尺度模型，其中，

为位置模型的分子，

为位置模型的分子，η为模型学习率：

其中，

表示位置滤波预期响应G_scale的共轭，

表示

的共轭。

Images