CN102722697B - 一种无人飞行器视觉自主导引着陆的目标跟踪方法 - Google Patents

Abstract

本发明公开了一种无人飞行器视觉自主导引着陆的目标跟踪方法，首先将目标在连续两帧图像间的运动幅值按照“coarse-to-fine”的顺序划分为多个分辨率等级，通过离线训练的方式模拟不同等级下的先验运动并计算其对应的先验误差Jacobian矩阵，同时，由于每一层Jacobian矩阵的求取结合了训练时的先验知识，可以保证算法在迭代搜索目标的过程中能够有效地跳出局部极值，避免跟踪失败；采用模板图像目标区域的稀疏特征即FAST角点处的灰度值描述目标，相较于传统Lucas-Kanade算法在计算中常常用目标区域的所有像素点稠密的表达目标，运算复杂度大大降低。

Description

一种无人飞行器视觉自主导引着陆的目标跟踪方法

技术领域

本发明涉及一种基于多分辨率运动先验的逆向复合目标跟踪方法，特别适用于无人飞行器视觉自主导引着陆过程中的目标稳定跟踪，属于数字图像处理领域。

背景技术

无人飞行器视觉自主着陆是无人飞行器控制研究领域里的热点问题。它利用数字图像处理技术获取定位参数，具有设备简单，代价小，获取的信息量大等诸多优点。相较于GPS和惯性导航，是完全自主和无源的。预定着陆点模板图像与机载实时图像的快速稳定匹配和跟踪是实现精确着陆控制的前提。无人机着陆过程中，目标匹配图像往往存在旋转、尺度及视角变化、局部遮挡和运动模糊等，因此对图像匹配和跟踪算法的实时性和鲁棒性均提出了极高的要求。

自Lucas和Kanade首次提出基于仿射变换模型的图像配准算法(AffineImage Alignment)[1](参见Lucas B，Kanade T.An iterative image registrationtechnique with an application to stereo vision[C].Proceedings ofthe7th InternationalJoint Conference on Artificial Intelligence，1981，2：674-679.)以来，Lucas-Kanade算法在目标跟踪，光流场估计，图像拼接，运动估计，图像配准，人脸识别等领域获得广泛的应用。

Lucas-Kanade算法的基本原理是：求取几何变换参数使得模板图像T(x)的像素点x=(x，y)经过变换W后映射在输入图像I中的变形图像I与原始模板图像间的差方和(SSD)最小，可用数学表达式描述为：

$\hat{p}=\underset{p}{\arg \min }\underset{x}{\mathrm{\Σ}}{[I\left(W(x;p)\right)-T\left(x\right)]}^2---\left(1\right)$

上式也称为Lucas-Kanade算法的最小二乘目标函数。式中W(x；p)为几何变换函数，T(x)为模板在坐标x处的像素点灰度值，I(W(x；p))表示输入图像在坐标x经变换后的新坐标W(x；p)处的灰度值。

最小化目标函数实际上是一个非线性最优化问题。为了求解这一最优化问题，假设已知变换参数初始值然后反复迭代累加Δp，并确保每次迭代关于Δp的近似目标表达式：

$O\left(\mathrm{\Δp}\right)=\underset{x}{\mathrm{\Σ}}{[I\left(W(x;p+\mathrm{\Δp})\right)-T\left(x\right)]}^2---\left(2\right)$

最小。每次迭代结束后，对变换参数p进行一次更新：

p←p+Δp  (3)

如此反复执行最小化和参数更新操作，直到迭代收敛为止。

在算法迭代过程中，为了求取变换参数增量Δp，首先将(2)式在Δp=0处附近进行一阶Tayor展开，可得：

$O\left(\mathrm{\Δp}\right)=\underset{x}{\mathrm{\Σ}}{[I\left(W(x;p)\right)+\▿I\frac{\∂W}{\∂p}\mathrm{\Δp}-T\left(x\right)]}^2---\left(4\right)$

其中，表示输入图像I在新坐标W(x；p)处的梯度向量，为变换函数W(x；p)关于其参数向量p的Jacobian矩阵，定义W(x；p)=(W_x，W_y)^T，则有：

$\frac{\∂W}{\∂p}=\left(\begin{array}{cccc}\frac{{\∂W}_x}{{\∂p}_1}& \frac{{\∂W}_x}{{\∂p}_2}& ...& \frac{{\∂W}_x}{{\∂p}_n}\\ \frac{{\∂W}_y}{{\∂p}_1}& \frac{{\∂W}_y}{{\∂p}_2}& ...& \frac{{\∂W}_y}{{\∂p}_n}\end{array}\right)---\left(5\right)$

由式(4)易知，Δp可通过求取该方程的最小二乘解得到。因而，求(4)式关于Δp的偏导，可得：

$2\underset{x}{\mathrm{\Σ}}{[\▿I\frac{\∂W}{\∂p}]}^T[I\left(W(x;p)\right)+\▿I\frac{\∂W}{\∂p}\mathrm{\Δp}-T\left(x\right)]=0---\left(6\right)$

解方程(6)可直接得到Δp闭合形式的解

$\mathrm{\Δp}=H^{-1}\underset{x}{\mathrm{\Σ}}{[\▿I\frac{\∂W}{\∂p}]}^T[T\left(x\right)-I\left(W(x;p)\right)]---\left(7\right)$

其中，H为n×n维Hessian矩阵，如下式所示：

$H=\underset{x}{\mathrm{\Σ}}{[\▿I\frac{\∂W}{\∂p}]}^T[\▿I\frac{\∂W}{\∂p}]---\left(8\right)$

由Lucas-Kanade算法的具体步骤可知，和的计算均与几何变换参数向量p相关，而p的计算是一个迭代更新的过程，W(s；p)在迭代过程中不断变化导致每次迭代都需要重新计算Hessian矩阵，算法效率较低。

为了提高算法效率，后续改进算法主要从如何降低每次迭代过程中Hessian矩阵的计算量入手，主要有前向复合算法[2](参见Shum HY，Szeliski R.Construction of panoramic image mosaics with global and local alignment[J].International Journal of Computer Vision，2000，16(1)：63-84.)和逆向复合算法[3](参见Baker S，Matthews I.Equivalence and efficiency of image alignmentalgorithms[C].In Proceedings of the IEEE Conference on Computer Vision andPattem Recognition.Kauai：IEEE，2001，1090-1097.)两种。

前向复合算法不需要在每一步迭代过程中重新计算但仍然需要在每步迭代过程中重新计算Hessian矩阵H，其参数更新的计算复杂度相较基本Lucas-Kanade算法优化有限，计算效率相当。逆向复合算法中，Hessian矩阵是个常量，只需在迭代开始前预先计算一次即可，因而逆向复合算法效率大大高于前面两种算法，可以实现很高的实时性能。

无人飞行器自主视觉导引着陆阶段，由于机载摄像机相对于目标的高速运动，目标在实时视频相邻帧图像中的位置和形状都将发生较大幅度的变化，而原始的逆向复合跟踪算法由于原理上的限制，只能对在相邻帧图像之间存在较小运动的目标实现稳定跟踪，无法在这种高动态飞行状态下实现对地面目标的可靠跟踪。

发明内容

有鉴于此，本发明提供了一种无人飞行器视觉自主导引着陆的目标跟踪方法，能够在不增加算法在线计算复杂度的前提下实现高动态飞行环境下对地目标的稳定跟踪，适应目标存在尺度缩放、旋转、光照变化等极端环境。

本发明的一种无人飞行器视觉自主导引着陆的目标跟踪方法，包括如下步骤：

步骤1、机载摄像机采集着陆目标点模板图像并对模板图像进行仿射光照归一化处理，得到归一化之后的模板图像I_norm(x)，其中x表示模板图像中像素点的坐标；

步骤2、对归一化后的模板图像I_norm(x)进行特征提取，得到具有N个特征点的特征模板图像；

步骤3、离线方式训练模板图像在m组不同运动范围下对应的先验误差Jacobian矩阵，具体方法如下：

设定无人飞行器的运动幅度范围S，将运动幅度范围S从小到大均分成m段运动范围，针对特征模板图像在每一段运动范围，采用透视变换方法对特征模板图像随机变形：即将特征模板图像的4个顶点在该段运动范围内分别随机移动一次作为一次随机变形，得到一幅变形后的特征模板图像，将特征模板图像按照上述方式移动N_p次，且满足N_p>>N，得到特征模板图像与随机变形后的N_p幅图像之间的透视变换关系矩阵 $\left[\begin{array}{cccc}\mathrm{\δ}{\mathrm{\μ}}_1^1& \mathrm{\δ}{\mathrm{\μ}}_1^2& ...& \mathrm{\δ}{\mathrm{\μ}}_1^{N_p}\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ \mathrm{\δ}{\mathrm{\μ}}_j^1& \mathrm{\δ}{\mathrm{\μ}}_j^2& ...& \mathrm{\δ}{\mathrm{\μ}}_j^{N_p}\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ \mathrm{\δ}{\mathrm{\μ}}_n^1& \mathrm{\δ}{\mathrm{\μ}}_n^2& ...& \mathrm{\δ}{\mathrm{\μ}}_n^{N_p}\end{array}\right],$ 计算随机变形后的N_p幅特征模板图像相对于变形前的特征模板图像的对应特征点之间的灰度误差向量矩阵 $\left[\begin{array}{cccc}{e}_1^1& e_1^2& ...& e_1^{N_p}\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ e_j^1& e_j^2& ...& e_j^{N_p}\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ e_N^1& e_N^2& ...& e_N^{N_p}\end{array}\right];$ 其中，透视变换关系矩阵中的第k列 ${\mathrm{\δ\μ}}^k={[{\mathrm{\δ\μ}}_1^k,{\mathrm{\δ\μ}}_2^k,...,{\mathrm{\δ\μ}}_n^k]}^T$ 表示第k次变换得到的图像与特征模板图像之间的透视变换关系对应的参数向量，灰度误差向量矩阵中的第k列表示第k次变换得到的图像与特征模板图像对应特征点之间的灰度误向量，k∈[1，N_p]，n=8；

然后再根据公式 $\left[\begin{array}{ccccc}{a}_{11}& ...& a_{1j}& ...& a_{1N}\\ .& .& .& .& .\\ .& .& .& .& .\\ .& .& .& .& .\\ a_{i1}& ...& a_{\mathrm{ij}}& ...& a_{\mathrm{iN}}\\ .& .& .& .& .\\ .& .& .& .& .\\ .& .& .& .& .\\ a_{n1}& ...& a_{\mathrm{nj}}& ...& a_{\mathrm{nN}}\end{array}\right]\left[\begin{array}{cccc}{e}_1^1& e_1^2& ...& e_1^{N_p}\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ e_j^1& e_j^2& ...& e_j^{N_p}\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ e_N^1& e_N^2& ...& e_N^{N_p}\end{array}\right]=\left[\begin{array}{cccc}{\mathrm{\δ\μ}}_1^1& {\mathrm{\δ\μ}}_1^2& ...& {\mathrm{\δ\μ}}_1^{N_p}\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ {\mathrm{\δ\μ}}_j^1& {\mathrm{\δ\μ}}_j^2& ...& {\mathrm{\δ\μ}}_j^{N_p}\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ {\mathrm{\δ\μ}}_n^1& {\mathrm{\δ\μ}}_n^2& ...& {\mathrm{\δ\μ}}_n^{N_p}\end{array}\right]$ 求得每一段运动范围对应的先验误差Jacobian矩阵A_h[q]，q∈[1，m]，其中

$A_h\left[q\right]=\left[\begin{array}{ccccc}{a^q}_{11}& ...& {a^q}_{1j}& ...& {a^q}_{1N}\\ .& .& .& .& .\\ .& .& .& .& .\\ .& .& .& .& .\\ {a^q}_{i1}& ...& {a^q}_{\mathrm{ij}}& ...& {a^q}_{\mathrm{iN}}\\ .& .& .& .& .\\ .& .& .& .& .\\ .& .& .& .& .\\ {a^q}_{n1}& ...& {a^q}_{\mathrm{nj}}& ...& {a^q}_{\mathrm{nN}}\end{array}\right],i\∈[1,n];$

步骤4、实际跟踪过程之初，对机载摄像机在t₀时刻获得的图像采用步骤1的方法进行仿射光照归一化处理，然后确定步骤2得到的特征模板图像在t₀时刻图像帧中的位置，将该位置在t₀时刻图像帧中的区域图像定义为目标跟踪区域，采用步骤2的特征提取方法对目标跟踪区域图像进行特征提取，得到目标跟踪模板图像I(x,t₀)；

步骤5、对无人飞行器获得的t时刻的实时输入图像采用步骤1的方法进行仿射光照归一化处理，得到当前输入图像；然后再进行在线跟踪，具体步骤如下：

5-1、通过最优化目标函数O(μ(t))=||I(f(x；μ(t))，t)-I(x,t₀)||²计算目标跟踪模板图像在t时刻的变换参数μ(t)，使该变换参数μ(t)对目标跟踪模板图像I(x,t₀)进行变换后得到的图像I(f(x；μ(t))，t)与目标跟踪模板图像I(x，t₀)的差异最小；

5-2、采用步骤5-1计算得到的变换参数μ(t)对目标跟踪模板图像I(xt₀)进行变换得到图像I(f(x；μ(t))，t₀)，以t时刻的变换参数μ(t)代替t+τ时刻的变换参数μ(t+τ)，计算图像i(f(x；μ(t))，t₀)在t+τ时刻的输入图像I(f(x；μ(t))，t+τ)，其中τ是无人飞行器采集图像的最小时间间隔；

5-3、计算输入图像I(f(x；μ(t))，t+τ)与目标跟踪模板图像I(x,t₀)之间的对应特征点之间的灰度误差图像e(t)=I(x,t₀)-I(f(x；μ(t))，t+τ)；令q=m，r=1；

5-4、计算Δμ(t+τ)=A_h[q]e(t)；

5-5、首先计算当前输入图像和目标跟踪模板图像I(x，t₀)之间的透视变换矩阵F(μ(t))，再计算t+τ时刻的输入图像相对于目标跟踪模板图像的透视变换矩阵：F(μ(t+τ))=F(μ(t))F^-1(Δμ(t+τ))；

5-6、判断Δμ(t+τ)与ε的关系：

当Δμ(t+τ)≥ε时，计算e(t)=I(x,t₀)-I(F(μ(t+τ))x，t+τ)，r=r+1，其中I(F(μ(t+τ))x，t+τ)表示目标跟踪模板图像I(x,t₀)经过透视变换矩阵F(μ(t+τ))变换的图像；

再判断r与nt_iter的关系：

当r≤n_iter时，执行步骤5-4；

当r>n_iter时，执行步骤5-7；

当Δμ(t+τ)<ε，q=q-1；

当q>1时，执行步骤5-4；

当q=1时，执行步骤5-7；

其中n_iter为迭代次数，取5≤n_iter≤15；ε为跟踪精度阈值矢量，其维数与t+τ时刻的Δμ(t+τ)相同，ε=[ε₁，ε₂，...，ε_n]^T，其中ε₁=ε₂=...=ε_n；

5-7、将F(μ(t+τ))发送给无人飞行器控制系统，无人飞行器根据F(μ(t+τ))进行目标跟踪。

本发明的一种无人飞行器视觉自主导引着陆的目标跟踪方法具有如下有益效果：

(1)针对无人飞行器视觉自主着陆飞行阶段，目标区域在实时图像中的尺度及位置剧烈变化这一高动态环境下，逆向复合算法难以实现稳定跟踪的问题，首先将目标在连续两帧图像间的运动幅值按照“coarse-to-fine”的顺序划分为多个分辨率等级，通过离线训练的方式模拟不同等级下的先验运动并计算其对应的先验误差Jacobian矩阵。基于多分辨率分层思想的目标搜索策略一方面可以大大增加算法的收敛范围，使得高动态环境下的目标跟踪成为现实；同时，由于每一层Jacobian矩阵的求取结合了训练时的先验知识，可以保证算法在迭代搜索目标的过程中能够有效地跳出局部极值，避免跟踪失败。

(2)采用模板图像目标区域的稀疏特征即FAST角点处的灰度值描述目标，相较于传统Lucas-Kanade算法在计算中常常用目标区域的所有像素点稠密的表达目标，运算复杂度大大降低。

(3)通过对图像进行光照归一化处理以提高算法在不同光照环境下跟踪结果的稳定性。

附图说明

图1为本发明的生成透视变换的先验运动示意图。

图2为本发明的步骤5的流程图。

具体实施方式

下面结合附图并举实施例，对本发明进行详细描述。

步骤1、机载摄像机采集着陆目标点模板图像并对模板图像进行仿射光照归一化处理，得到归一化之后的模板图像在x点处的像素点灰度值I_norm(x)，其中x表示模板图像中像素点的坐标；

对于视频序列，可将目标运动看作是二维空间与时间维度上的三维运动，为此，将t时刻实时输入图像在x=(x，y)处像素点的灰度值用I(x，t)表示。将某t_x时刻的实时视频图像选作参考图像，对于在参考图像中选定的目标区域，可以用含有N个元素的点集表示，并将该点集对应的图像灰度值列向量I(x,t_x)=[I(x₁,t_x)，I(x₂，t_x)，...，I(x_N,t_x)]^T定义为模板图像。

分别计算出输入模板图像I(x,t₀)的像素点灰度值的均值和标准差μ(i)和σ(I)，归一化之后的图像在x点处的像素灰度值为：

$I_{\mathrm{norm}}\left(x\right)=\frac{I\left(x\right)-\mathrm{\&mu;}\left(I\right)}{\mathrm{\&sigma;}\left(I\right)}$

其中I_norm(x)表示归一化之后的图像在x点处的像素点灰度值，I(x)表示目标模板图像在x点处的像素点灰度值。

步骤2、对归一化后的模板图像I_norm(x)进行特征提取，得到具有N个特征点的特征模板图像；

将模板图像均分成s×s个局部子区域，针对每一个局部子区域，分别以该局部子区域中的每一个像素点为圆心选取一个半径为3个像素点的Bresenham圆，当圆周上有连续9个以上像素点的灰度值与圆心点处像素点灰度值的差值同时大于或同时小于阈值T时，该圆心的像素点则确定为候选角点，用所有符合上述条件的像素点灰度值来描述的图像即为特征模板图像，s取值的大小决定了特征模板图像上FAST-9特征点的数目，一般情况下，640*480大小图像，取s=160-200；T取值的大小决定了FAST-9特征点的灵敏度，一般情况下，取50<T<100。T取值越大，特征点越少；T取值越小，特征点越多。

步骤3：离线方式训练模板图像在m组不同运动范围下对应的先验误差Jacobian矩阵A_h[q]，q∈[1，m]：

对于模板图像，人为地给出一个先验运动变换“扰动量”Δμ，用以模拟无人飞行器在接近目标的过程中由于存在复杂的滚动、俯仰和偏航等姿态变化引起的目标模板在后续视频图像中的变形。

在一定距离和角度范围内，可将地面目标近似成平面物体，由视觉成像几何原理，这种情况下的几何变形采用透视变换来实现图像的变形，对于透视变换来说，包含8个独立参数变量，因此，随机变换通过在一定运动范围内随机移动模板的4个顶点实现。

如图1所示，本实施例中，设计特征模板图像4个顶点随机移动的范围为0-100个像素点，即设无人飞行器的运动幅度范围S为100个像素点，将100个像素点均分成5段，即运动范围分为0-20像素点，20-40像素点、40-60像素点、60-80像素点和80-100像素点共5段，在每一段运动范围中，将4个顶点在运动范围中分别移动一次(假如在0-20个像素点的运动范围内，即是4个顶点的移动范围均在0-20个像素点范围内随机移动)，得到一个变形后的图像，如此随机地移动N_p(N_p>>N)次，得到N_p个变形后的图像，则得到特征模板图像与随机变形后的N_p幅图像之间的透视变换关系矩阵可表示为：

$\left[\begin{array}{cccc}\mathrm{\&delta;}{\mathrm{\&mu;}}_1^1& \mathrm{\&delta;}{\mathrm{\&mu;}}_1^2& ...& \mathrm{\&delta;}{\mathrm{\&mu;}}_1^{N_p}\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ \mathrm{\&delta;}{\mathrm{\&mu;}}_j^1& \mathrm{\&delta;}{\mathrm{\&mu;}}_j^2& ...& \mathrm{\&delta;}{\mathrm{\&mu;}}_j^{N_p}\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ \mathrm{\&delta;}{\mathrm{\&mu;}}_n^1& \mathrm{\&delta;}{\mathrm{\&mu;}}_n^2& ...& \mathrm{\&delta;}{\mathrm{\&mu;}}_n^{N_p}\end{array}\right],$ 并可计算随机变形后的特征模板图像相对于变形前的特征模板图像的对应特征点之间的灰度误差向量矩阵 $\left[\begin{array}{cccc}{e}_1^1& e_1^2& ...& e_1^{N_p}\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ e_j^1& e_j^2& ...& e_j^{N_p}\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ e_N^1& e_N^2& ...& e_N^{N_p}\end{array}\right];$ 其中，透视变换关系矩阵中的第k列表示第k次变换得到的图像与特征模板图像之间的透视变换关系对应的参数向量，灰度误差向量矩阵中的第k列表示第k次变换得到的图像与特征模板图像对应特征点之间的灰度误向量，k∈[1，N_p]，由于特征模板图像移动的4个顶点具有8个坐标，因此n=8；

再根据公式 $\left[\begin{array}{ccccc}{a}_{11}& ...& a_{1j}& ...& a_{1N}\\ .& .& .& .& .\\ .& .& .& .& .\\ .& .& .& .& .\\ a_{i1}& ...& a_{\mathrm{ij}}& ...& a_{\mathrm{iN}}\\ .& .& .& .& .\\ .& .& .& .& .\\ .& .& .& .& .\\ a_{n1}& ...& a_{\mathrm{nj}}& ...& a_{\mathrm{nN}}\end{array}\right]\left[\begin{array}{cccc}{e}_1^1& e_1^2& ...& e_1^{N_p}\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ e_j^1& e_j^2& ...& e_j^{N_p}\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ e_N^1& e_N^2& ...& e_N^{N_p}\end{array}\right]=\left[\begin{array}{cccc}{\mathrm{\&delta;\&mu;}}_1^1& {\mathrm{\&delta;\&mu;}}_1^2& ...& {\mathrm{\&delta;\&mu;}}_1^{N_p}\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ {\mathrm{\&delta;\&mu;}}_j^1& {\mathrm{\&delta;\&mu;}}_j^2& ...& {\mathrm{\&delta;\&mu;}}_j^{N_p}\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ {\mathrm{\&delta;\&mu;}}_n^1& {\mathrm{\&delta;\&mu;}}_n^2& ...& {\mathrm{\&delta;\&mu;}}_n^{N_p}\end{array}\right]$ 求得每一段运动范围对应的先验误差Jacobian矩阵A_h[q]，q∈[1，m]，其中

$A_h\left[q\right]=\left[\begin{array}{ccccc}{a^q}_{11}& ...& {a^q}_{1j}& ...& {a^q}_{1N}\\ .& .& .& .& .\\ .& .& .& .& .\\ .& .& .& .& .\\ {a^q}_{i1}& ...& {a^q}_{\mathrm{ij}}& ...& {a^q}_{\mathrm{iN}}\\ .& .& .& .& .\\ .& .& .& .& .\\ .& .& .& .& .\\ {a^q}_{n1}& ...& {a^q}_{\mathrm{nj}}& ...& {a^q}_{\mathrm{nN}}\end{array}\right],i\&Element;[1,n];$

可求取矩阵A_h在如下最小二乘意义下的解：

$O\left(A_h\right)={\mathrm{\&Sigma;}}_{k=1}^{N_p}(\mathrm{\&Delta;}{\mathrm{\&mu;}}^k-A_h{e}^k)$

其中，k∈[1，N_P]， ${\mathrm{\&Delta;\&mu;}}^k={[{\mathrm{\&delta;\&mu;}}_1^k,{\mathrm{\&delta;\&mu;}}_2^k,...,{\mathrm{\&delta;\&mu;}}_n^k]}^T,e^k={[e_1^k,e_2^k,...e_N^k]}^T.$

令 $U={[{\mathrm{\&Delta;\&mu;}}^1,{\mathrm{\&Delta;\&mu;}}^2,...,{\mathrm{\&Delta;\&mu;}}^{N_p}]}_{n\&times;N_p},E={[e^1,e^2,...,e^{N_p}]}_{N\&times;N_p},$ 得如下方程组：

A_hE=U

并最终可求得矩阵A_h的最小二乘解：

A_h=UE^T(EE^T)^-1

步骤4、在实际跟踪过程中，对机载摄像机在t₀时刻获得的图像采用步骤1的方法进行仿射光照归一化处理，然后确定步骤2得到的特征模板图像在t₀时刻图像帧中的位置，将该位置在t₀时刻图像帧中的区域图像定义为目标跟踪区域，采用步骤2的特征提取方法对目标跟踪区域图像进行特征提取，得到目标跟踪模板图像，具体方法如下：

(1)提取模板图像的SURF特征点和SURF特征描述子；

(2)提取实时图像的SURF特征点和SURF特征描述子；

(3)比较特征描述子之间的欧式距离，确定正确匹配对；

(4)利用RANSAC算法选定模板图像在t₀时刻图像帧中的位置，并确定目标跟踪区域；

(5)对于在t₀时刻图像帧中选定的目标跟踪区域，用含有N个元素的稀疏FAST-9特征点集x={x₁，x₂，...，x_N}表示，并将该点集对应的图像灰度值列向量I(x,t₀)=[I(x₁，t₀)，I(x₂，t₀)，...，(x，t₀)]^T定义为目标跟踪模板图像。

步骤5、对无人飞行器获得的t时刻的实时输入图像采用步骤1的方法进行仿射光照归一化处理，得到当前输入图像；然后采用“coarse-to-fine”多分辨率分层训练的思想，按照运动范围由大到小的顺序，分别对当前输入图像进行搜索，最后计算用于控制无人飞行器跟踪目标的透视变换矩阵，具体步骤如下：

5-1、通过最优化目标函数O(μ(t))=||If(x；μ(t))，t)-I(x，t₀)||²计算目标跟踪模板图像在t时刻的变换参数μ(t)，使该变换参数μ(t)对目标跟踪模板图像I(x，t₀)进行变换后得到的图像I(f(x；μ(t))，t)与目标跟踪模板图像I(x,t₀)的差异最小；

在无人飞行器自主视觉导引着陆飞行过程中，由于机载摄像机与目标区域之间存在的空间相对运动(位置和姿态)，目标区域对应的模板图像将在实时图像中产生相对于参考图像的平移运动以及几何变形等，我们利用几何变换(仿射或者透视变换)x→f(x；μ(t))来描述这种相对运动。其中μ(t)=[μ₁(t)，μ₂(t)，...，μ_n(t)]^T为描述这种几何变换的n维参数列向量(N>n)，且变换函数f对μ、x均可微。对于目标模板，有f(x；μ(t₀))=x。

根据环境光照不变性假设可知，物理场景中的空间点在跟踪过程中亮度保持不变，这样，视频序列中的目标跟踪问题便转化为下式所示的的几何变换参数μ(t)的求解问题：

I(f(x；μ(t))，t)=I(x，t₀)  (9)

式中I(f(x；μ(t))，t)=[I(f(x₁；μ(t))，t)，I(f(x₂；μ(t))，t，...，I(f(x_N；μ(t)),t)]^T。

相应的，t时刻的变换参数μ(t)可通过最优化如下目标函数得到：

O(μ(t))=||I(f(x；μ(t))，t)-I(x，t₀)||²  (10)

5-2、采用步骤5-1计算得到的变换参数μ(t)对目标跟踪模板图像I(x,t₀)进行变换得到图像I(f(x；μ(t))，t₀)，以t时刻的变换参数μ(t)代替t+τ时刻的变换参数μ(t+τ)，计算图像I(f(x；μ(t))，t₀)在t+τ时刻的输入图像I(f(x；μ(t))，t+τ)，其中τ是无人飞行器采集图像的最小时间间隔；

对于视频序列，目标跟踪问题可以进一步描述为：假定t时刻的几何变换参数μ(t)已知，即I(f(x；μ(t))，t)=I(x,t₀)，求取t+τ时刻的变换参数μ(t+τ)，使得该时刻的实时输入图像I(f(x；μ(t+τ))，t+τ)=I(x,t₀)。

可将式(10)所示的目标函数改写为如下关于时间t和变换参数向量μ(t)的增量形式：

O(Δμ(t+τ))=||I(f(x；μ(t+τ))，t+τ)-I(x，t₀)||²  (11)=||I(f(x；μ(t)+Δμ(t+τ))，t+τ)-I(x,t₀)||²

求解上述方程，可得t时刻变换参数增量Δμ(t)的最小二乘解：

Δμ(t+τ)=A(t)[I(x,t₀)-I(f(x；μ(t))，t+τ)]  (12)

式中A(t)=(R^T(μ(t)，t)R(μ(t)，t))^-1R^T(μ(t)，t)。

上式中， $R(\mathrm{\&mu;}\left(t\right),t)={[I_{{\mathrm{\&mu;}}_1}(f(x;\mathrm{\&mu;}\left(t\right)),t)I_{{\mathrm{\&mu;}}_2}(f(x;\mathrm{\&mu;}\left(t\right)),t),...,I_{{\mathrm{\&mu;}}_n}(f(x;\mathrm{\&mu;}\left(t\right)),t)]}_{N\&times;n}$ 为向量I(f(x；μ(t))，t)关于向量μ(t)的N×n维Jacobian矩阵。

根据逆向复合模板跟踪原理，式(11)的目标函数改写为：

O(Δμ(t+τ))=||I(f(x；μ(t))，t+τ)-I(f(x；Δμ(t+τ))，t₀)||²  (13)

同样，可由上式求得t时刻变换参数增量Δμ(t)的最小二乘解：

Δμ(t+τ)=A(t₀)[I(x，t₀)-I(f(x；μ(t))，t+τ)]  (14)

式中，A(t₀)=(R^T(x,t₀)R(x,t₀))^-1R^T(x,t₀)，通常将其称为跟踪误差Jacobian矩阵。

式中， $R(x,t_0)={[I_{{\mathrm{\&mu;}}_1}(x,t_0),I_{{\mathrm{\&mu;}}_2}(x,t_0),...,I_{{\mathrm{\&mu;}}_n}(x,t_0)]}_{N\&times;n}$ 为向量I(f(x；μ(t))，t)关于向量μ(t)的在t₀时刻的N×n维Jacobian矩阵。

5-3计算输入图像I(f(x；μ(t))，t+τ)与目标跟踪模板图像I(x，t₀)之间的对应特征点之间的灰度误差图像e(t)=I(x，t₀)-I(f(x；μ(t))，t+τ)；令q=m，r=1；

5-4、计算Δμ(t+τ)=A_h[q]e(t)；

5-5、首先计算当前输入图像和目标跟踪模板图像I(x，t₀)之间的透视变换矩阵F(μ(t))，再计算t+τ时刻的输入图像相对于目标跟踪模板图像的透视变换矩阵：F(μ(t+τ))=F(μ(t))F^-1(Δμ(t+τ))；

其中，F^-1(Δμ(t+τ))的计算方法如下：由于透视变换矩阵 $F\left(\mathrm{\&mu;}\left(t\right)\right)=\left[\begin{array}{ccc}{p}_{1\left(t\right)}& p_{2\left(t\right)}& p_{3\left(t\right)}\\ p_{4\left(t\right)}& p_{5\left(t\right)}& p_{6\left(t\right)}\\ p_{7\left(t\right)}& p_{8\left(t\right)}& 1\end{array}\right]$ 为包含8个独立变量的3x3可逆单应矩阵，则对应的t时刻透视变换对应的参数向量μ(t)=(p_1(t)，p_2(t)，p_3(t)，p_4(t)，p_5(t)，p_6(t)，p_7(t)，p_8(t))^T，因此，可将μ(t)中向量还原到F(μ(t))中，求得F(μ(t))；同理，在步骤5-4求得的Δμ(t+τ)的基础上，根据上述原理计算得到F(Δμ(t+τ))，再对其求逆矩阵，可最终求得F^-1(Δμ(t+τ))的值。

5-6、判断Δμ(t+τ)与ε的关系：

当Δμ(t+τ)≥ε时，计算e(t)=I(x，t₀)-I(F(μ(t+τ))x,t+τ)，r=r+1，其中I(F(μ(t+τ))x，t+τ)表示目标跟踪模板图像I(xt₀)经过透视变换矩阵F(μ(t+τ))变换的图像；

再判断r与n_iter的关系：

当r≤n_iter时，执行步骤5-4；

当r>n_iter时，执行步骤5-7；

当Δμ(t+τ)<ε，计算F(μ(t+τ))，q=q-1；

当q>1时，执行步骤5-4；

当q=1时，执行步骤5-7；

其中n_iter为迭代次数，取5≤n_iter≤15；ε为跟踪精度阈值矢量，其维数与t+τ时刻的Δμ(t+τ)相同，ε=[ε₁，ε₂，...，ε_n]^T，其中ε₁=ε₂=...=ε_n；取0.5<ε_n<0.8。

5-7、将F(μ(t+τ))发送给无人飞行器控制系统，无人飞行器根据F(μ(t+τ))进行目标跟踪。

在上述的步骤5-1至5-7中，首先利用最大运动范围对应的Jacobian矩阵求取t+τ时刻的变换参数的增量Δμ(t+τ)，相当于在一个较大的区域范围内搜索目标。再将Δμ(t+τ)与设定的跟踪精度阈值矢量ε进行比较，如果Δμ(t+τ)大于等于ε，说明该当前输入图像属于此时Jacobian矩阵对应的运动范围，利用该Jacobian矩阵使算法在有限的迭代次数内收敛，使Δμ(t+τ)逐渐减小，直至Δμ(t+τ)小于ε。求取此时透视变换矩阵F(μ(t+τ))，控制无人飞行器进行目标跟踪。但此时的跟踪精度不高，只相当于圈定了一个包含目标的较小区域。继续在此较小区域内利用低一层的Jacobian矩阵重复迭代，进一步逼近目标的精确位置。依此类推，直至算法在最底层Jacobian矩阵对应的区域范围内收敛，即可认为完成了对目标的精确跟踪。这种基于多分辨率分层思想的目标搜索策略一方面可以大大增加算法的收敛范围，使的高动态环境下的目标跟踪成为现实；另一方面，由于离线已经训练出不同运动范围对应的Jacobian矩阵，在本步骤中只需多次迭代寻找当前输入图像对应的最小运动范围即可，避免跟踪失败。

综上所述，以上仅为本发明的较佳实施例而已，并非用于限定本发明的保护范围。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (3)

1.一种无人飞行器视觉自主导引着陆的目标跟踪方法，其特征在于，包括如下步骤：

步骤1、机载摄像机采集着陆目标点模板图像并对模板图像进行仿射光照归一化处理，得到归一化之后的模板图像I_norm(x)，其中x表示模板图像中像素点的坐标；

步骤2、对归一化后的模板图像I_norm(x)进行特征提取，得到具有N个特征点的特征模板图像；

步骤3、离线方式训练模板图像在m组不同运动范围下对应的先验误差Jacobian矩阵,具体方法如下：

设定无人飞行器的运动幅度范围S，将运动幅度范围S从小到大均分成m段运动范围，针对特征模板图像在每一段运动范围，采用透视变换方法对特征模板图像随机变形：即将特征模板图像的4个顶点在该段运动范围内分别随机移动一次作为一次随机变形，得到一幅变形后的特征模板图像，将特征模板图像按照上述方式移动N_p次，且满足N_p>>N，得到特征模板图像与随机变形后的N_p幅图像之间的透视变换关系矩阵 $\left[\begin{array}{cccc}\mathrm{\&delta;}{\mathrm{\&mu;}}_1^1& \mathrm{\&delta;}{\mathrm{\&mu;}}_1^2& ...& \mathrm{\&delta;}{\mathrm{\&mu;}}_1^{N_p}\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ \mathrm{\&delta;}{\mathrm{\&mu;}}_j^1& \mathrm{\&delta;}{\mathrm{\&mu;}}_j^2& ...& \mathrm{\&delta;}{\mathrm{\&mu;}}_j^{N_p}\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ \mathrm{\&delta;}{\mathrm{\&mu;}}_n^1& \mathrm{\&delta;}{\mathrm{\&mu;}}_n^2& ...& \mathrm{\&delta;}{\mathrm{\&mu;}}_n^{N_p}\end{array}\right],$ 计算随机变形后的N_p幅特征模板图像相对于变形前的特征模板图像的对应特征点之间的灰度误差向量矩阵 $\left[\begin{array}{cccc}{e}_1^1& e_1^2& ...& e_1^{N_p}\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ e_j^1& e_j^2& ...& e_j^{N_p}\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ e_N^1& e_N^2& ...& e_N^{N_p}\end{array}\right],$ 其中，透视变换关系矩阵中的第k列表示第k次变换得到的图像与特征模板图像之间的透视变换关系对应的参数向量，灰度误差向量矩阵中的第k列表示第k次变换得到的图像与特征模板图像对应特征点之间的灰度误差向量，k∈[1,N_p]，n＝8；

然后再根据公式 $\left[\begin{array}{ccccc}{a}_{11}& ...& a_{1j}& ...& a_{1N}\\ .& .& .& .& .\\ .& .& .& .& .\\ .& .& .& .& .\\ a_{i1}& ...& a_{\mathrm{ij}}& ...& a_{\mathrm{iN}}\\ .& .& .& .& .\\ .& .& .& .& .\\ .& .& .& .& .\\ a_{n1}& ...& a_{\mathrm{nj}}& ...& a_{\mathrm{nN}}\end{array}\right]\left[\begin{array}{cccc}{e}_1^1& e_1^2& ...& e_1^{N_p}\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ e_j^1& e_j^2& ...& e_j^{N_p}\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ e_N^1& e_N^2& ...& e_N^{N_p}\end{array}\right]=\left[\begin{array}{cccc}\mathrm{\&delta;}{\mathrm{\&mu;}}_1^1& \mathrm{\&delta;}{\mathrm{\&mu;}}_1^2& ...& \mathrm{\&delta;}{\mathrm{\&mu;}}_1^{N_p}\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ \mathrm{\&delta;}{\mathrm{\&mu;}}_j^1& \mathrm{\&delta;}{\mathrm{\&mu;}}_j^2& ...& \mathrm{\&delta;}{\mathrm{\&mu;}}_j^{N_p}\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ \mathrm{\&delta;}{\mathrm{\&mu;}}_n^1& \mathrm{\&delta;}{\mathrm{\&mu;n}}_n^2& ...& \mathrm{\&delta;}{\mathrm{\&mu;}}_n^{N_p}\end{array}\right]$ 求得每一段运动范围对应的先验误差Jacobian矩阵A_h[q],q∈[1,m],其中

$A_h=\left[\begin{array}{ccccc}{a}_{11}& ...& a_{1j}& ...& a_{1N}\\ .& .& .& .& .\\ .& .& .& .& .\\ .& .& .& .& .\\ a_{i1}& ...& a_{\mathrm{ij}}& ...& a_{\mathrm{iN}}\\ .& .& .& .& .\\ .& .& .& .& .\\ .& .& .& .& .\\ a_{n1}& ...& a_{\mathrm{nj}}& ...& a_{\mathrm{nN}}\end{array}\right],i\&Element;[1,n];$

步骤4、实际跟踪过程之初，对机载摄像机在t₀时刻获得的图像采用步骤1的方法进行仿射光照归一化处理，然后确定步骤2得到的特征模板图像在t₀时刻图像帧中的位置，将该位置在t₀时刻图像帧中的区域图像定义为目标跟踪区域，采用步骤2的特征提取方法对目标跟踪区域图像进行特征提取，得到目标跟踪模板图像I(x,t₀)；

步骤5、对无人飞行器获得的t时刻的实时输入图像采用步骤1的方法进行仿射光照归一化处理，得到当前输入图像；然后再进行在线跟踪，具体步骤如下：

5-1、通过最优化目标函数O(μ(t))＝||I(f(x；μ(t)),t)-I(x,t₀)||²计算目标跟踪模板图像在t时刻的变换参数μ(t)，使该变换参数μ(t)对目标跟踪模板图像I(x,t₀)进行变换后得到的图像I(f(x；μ(t)),t)与目标跟踪模板图像I(x,t₀)的差异最小；

5-2、采用步骤5-1计算得到的变换参数μ(t)对目标跟踪模板图像I(x,t₀)进行变换得到图像I(f(x；μ(t)),t₀)，以t时刻的变换参数μ(t)代替t+τ时刻的变换参数μ(t+τ)，计算图像I(f(x；μ(t)),t₀)在t+τ时刻的输入图像I(f(x；μ(t)),t+τ)，其中τ是无人飞行器采集图像的最小时间间隔；

5-3、计算输入图像I(f(x；μ(t))，t+τ)与目标跟踪模板图像I(x,t₀)之间的对应特征点之间的灰度误差图像e(t)＝I(x，t₀)-I(f(x；μ(t))，t+τ)；令q＝m，r＝1；

5-4、计算Δμ(t+τ)＝A_h[q]e(t)；

5-5、首先计算当前输入图像和目标跟踪模板图像I(x,t₀)之间的透视变换矩阵F(μ(t))，再计算t+τ时刻的输入图像相对于目标跟踪模板图像的透视变换矩阵：F(μ(t+τ))＝F(μ(t))F^-1(Δμ(t+τ))；

5-6、判断Δμ(t+τ)与ε的关系：

当Δμ(t+τ)≥ε时，计算e(t)＝I(x,t₀)-I(F(μ(t+τ))x,t+τ)，r＝r+1，其中I(F(μ(t+τ))x,t+τ)表示目标跟踪模板图像I(x,t₀)经过透视变换矩阵F(μ(t+τ))变换的图像；

再判断r与n_iter的关系：

当r≤n_iter时，执行步骤5-4；

当r＞n_iter时，执行步骤5-7；

当Δμ(t+τ)＜ε，q＝q-1；

当q＞1时，执行步骤5-4；

当q＝1时，执行步骤5-7；

其中n_iter为迭代次数，取5≤n_iter≤15；ε为跟踪精度阈值矢量，其维数与t+τ时刻的Δμ(t+τ)相同，ε＝[ε₁,ε₂,...,ε_n]^T，其中ε₁＝ε₂＝...＝ε_n；

5-7、将F(μ(t+τ))发送给无人飞行器控制系统，无人飞行器根据F(μ(t+τ))进行目标跟踪。

2.如权利要求1所述的一种无人飞行器视觉自主导引着陆的目标跟踪方法，其特征在于，所述步骤2中采用FAST-9角点检测方法对模板图像进行特征提取，具体为：

将模板图像均分成s×s个局部子区域，针对每一个局部子区域，分别以该局部子区域中的每一个像素点为圆心选取一个半径为3个像素点的Bresenham圆，当圆周上有连续9个以上像素点的灰度值与圆心点处像素点灰度值的差值同时大于或同时小于阈值T时，该圆心的像素点则确定为候选角点，用所有符合上述条件的像素点灰度值来描述的图像即为特征模板图像，其中，s取值的大小与FAST-9特征点的数目成正比；T取值的大小与FAST-9特征点的灵敏度成反比。

3.如权利要求1所述的一种无人飞行器视觉自主导引着陆的目标跟踪方法，其特征在于，所述步骤4中的目标跟踪区域的确定方法为：

步骤301、提取特征模板图像的SURF特征点和SURF特征描述子；

步骤302、提取实时图像的SURF特征点和SURF特征描述子；

步骤303、比较特征模板图像的特征描述子和实时图像的特征描述子之间的欧式距离，确定正确匹配对；

步骤304、利用RANSAC算法选定特征模板图像在实时图像中的t₀时刻图像帧中的位置，并确定目标跟踪区域。