CN111028292B - 一种亚像素级图像匹配导航定位方法 - Google Patents

Abstract

本发明提供了一种亚像素级图像匹配导航定位方法，首先用哈尔小波对图像进行压缩，并利用仿射不变性特征来表征图像，提取边缘特征，然后结合高斯亚像素拟合进行特征点的亚像素级匹配定位。本发明具有很强的鲁棒性，能够减少因图像视角偏差较大所引起的匹配性能降低问题，所提出的高精度亚像素图像匹配方法能够高精度的重构出真实环境的纹理等细节信息，提高无人机定位和环境构建的准确性；可以大大降低计算量，提高匹配效率。

Description

一种亚像素级图像匹配导航定位方法

技术领域

本发明涉及视觉导航领域，尤其是一种图像匹配导航定位方法。通过具有仿射不变性的特征点来表征图像，并结合高斯亚像素拟合原理及特征描述符简化方法达到快速亚像素级定位，解决微小型无人机视觉导航对图像匹配算法高精度的需求。

背景技术

无人机视觉导航是指无人机借助视觉传感器感知环境，利用计算机对图像进行处理、分析和识别，通过与先验图像的匹配，来对无人机自身位姿进行估计，以实现测距、悬停、自主着陆和规划路径避障等任务。图像匹配是无人机视觉导航过程的核心技术之一，其匹配的速度和精度对导航系统的准确度和效率有着直接影响。目前大多数图像匹配算法是在某个离散化空间上估计一个视差集合，如整数像素视差，这对于只需要少量特征点的应用场合是足够的。而无人机视觉导航中往往需要处理超高分辨率的航拍图像，图像间存在光照条件、视角和尺度差异，并且要求亚像素级定位精度。因此，鲁棒性好且具有亚像素精度的图像匹配是提高无人机视觉导航系统整体性能的另一个关键问题。

在亚像素级图像匹配领域，国内外已经有一系列的研究成果。大致可以分为4类：1)扩展相位相关法：该类方法是对相位相关法的改进，不仅对灰度差和窄带噪声不敏感，而且还保留了相位相关法的优点。扩展相位相关法可分为时域求解法和频域求解法。在频域求解时匹配精度高，计算量小。2)基于插值的方法：该类方法主要通过灰度内插实现亚像素级匹配精度。分为相似性函数内插法和图像灰度直插法。该类方法的优点是精度高。3)最优化方法：该类方法选用适当的数学模型，通过迭代对参数进行优化，通过逐步优化得到亚像素精度的匹配。该类方法能够支持各种变换模型，匹配精度也比较高。4)基于特征点的方法。该类方法主要基于图像中的显著性特征(如点特征、线特征和区域特征)，通常结合内插法或相位相关法来使配准精度达到亚像素级。该类方法在保持鲁棒性的同时具有较高的匹配精度，因此被广泛应用于亚像素图像配准。

但无论哪种方法都各有优势和不足：1)扩展的相位相关法只能处理固定形式的图像变化，对于更复杂的非线性变换就无法达到匹配效果。2)插值法的缺点是计算量大，且匹配精度依赖于内插算法的质量。3)最优化法中寻找全局最优化的概率需要提高。 4)基于特征点的方法不适用于图像视角变化较大的情况。

综合以上研究与分析，各种方法在匹配精度、运行时间、和图像之间差异性造成的鲁棒性方面，存在较大的差异。因此，将图像匹配精度提高到亚像素级，减少算法复杂度，降低误匹配率同时增强算法对图像差异的鲁棒性，仍然是一个大的挑战。

发明内容

为了克服现有技术的不足，本发明提供一种亚像素级图像匹配导航定位方法。本发明解决的问题是能够为无人机提供一种亚像素级图像匹配定位方法，首先用哈尔小波对图像进行压缩，并利用仿射不变性特征来表征图像，提取边缘特征，然后结合高斯亚像素拟合进行特征点的亚像素级匹配定位。解决微小型无人机视觉导航对图像匹配算法高精度的需求。

目前大多数应用在无人机视觉导航中的匹配算法都只满足像素级精度，无法高精度的重构出真实环境的纹理等细节信息。此外，当无人机在视觉导航中常常会出现图像之间存在较大视角差异，这种情况下图像匹配性能会急剧下降，对无人机定位和环境构建的准确性进一步产生严重影响。仿射不变特征由于抗视角变换能力强，被越来越多应用在图像差异较大时的图像匹配。因此，这里采用具有仿射不变性的特征点来表征图像，并结合高斯亚像素拟合原理达到特征点的亚像素级定位，最终实现亚像素级精度的图像匹配定位。

本发明解决其技术问题所采用的技术方案的具体步骤如下：

步骤1：采用哈尔小波变换对输入图像进行压缩，首先执行步骤1.1；

1.1)将输入匹配图像进行小波变换，执行完以后转到步骤1.2；

1.2)运动符号编码器对图像进行编码，生成压缩图像，执行完以后转到步骤1.3)；

1.3)对编码完成的图像进行解码处理，执行完以后转到步骤1.4；

1.4)进行小波系数的反量化，执行完以后转到步骤1.5；

1.5)生成最终的重构图像，执行完以后转到步骤2.1；

假设原始图像的大小为2M×2N，记为I(x,y)，则经过哈尔小波变换可将原始图像I(x,y)分为四幅图，A为变换后的近似图像、B为水平细节图像、C为垂直细节图像、 D为对角线细节图像，通过在近似图像A中进行图像特征提取与匹配能够有效降低计算量，提高快速性；

步骤2：对压缩处理后的图像进行仿射变换模拟；

仿射模型含有六个自由度，其中两个平移变量tx、ty，一个相机旋转变量ψ，一个均匀缩放系数s，另外两个变量用于定义相机主光轴方向，称为经度φ和纬度θ；

2.1)对每一幅图像进行仿射变换,模拟由于摄像机位置变化产生的仿射形变；图像的仿射形变取决于经纬度的变化，因此在对原始图像做仿射变换时只对摄像机的经纬度进行模拟，在对经纬度采样时按照如下规则进行：

①纬度角θ和倾斜程度t满足一一对应关系，即θ＝arccos(1/t)，t的采样值选为 t＝1,a,a²,...,aⁿ，指数n为仿射变换的模拟次数，该值越大则模拟出的仿射变换的情况越多；

②经度角φ通过一组采样数据来实现模拟：0,b/t,...,kb/t(kb/t＜180°)，即一个t对应多个经度角。选择b为角度数，k是一个满足

的最大整数。执行完以后转到步骤3.1；

步骤3)：对每一幅模拟图像提取边缘特征点；

3.1)采用中值滤波法对图像进行噪声去除，以除去孤立的点，执行完以后转到步骤3.2；

3.2)采用SUSAN算法提取图像边缘轮廓,执行完以后转到步骤3.3；

3.3)将边缘像素连接起来，执行完以后转到步骤3.4；

3.4)在图像边缘上提取图像的特征点，得到粗定位边缘点，执行完以后转到步骤4.1；

步骤4)：亚像素级特征提取

4.1)沿着粗定位边缘点的梯度方向寻找距离最近的邻近点，当梯度方向与水平方向夹角θ小于45°时，Li＝1/cos(θ)为距离长度，执行完以后转到步骤4.2；

4.2)对邻近点进行灰度线性插值处理；最后通过邻近点的灰度值计算粗定位特征点的灰度差值；根据高斯亚像素拟合原理，图像边缘点灰度值变化呈高斯分布，高斯曲线的表达式如下：

式中，μ是中心点坐标，σ是高斯函数的标准差，u是粗定位点的灰度差值；执行完以后转到步骤4.3；

4.3)采用最小二乘法拟合并计算出μ和σ，假设粗定位边缘坐标为p(x₀,y₀)，则根据σ和求取得边缘点梯度方向夹角θ求出亚像素级边缘点的坐标为 p′(x₀+σcos(θ),y₀+σsin(θ))；执行完以后转到步骤5.1；

步骤5)：亚像素特征匹配

包括构建特征向量和特征向量的匹配，采用以下步骤来实现：

5.1)以亚像素边缘点为特征点(x,y)，按照下式计算其梯度和模值：

式中，L是每个特征点所在的模拟图像；在以特征点为中心的邻域窗口内进行采样，采用直方图统计领域像素的梯度方向，并用高斯函数对直方图进行平滑；梯度直方图描述的范围是0°～360°，每10°一个柱，共36个柱，以直方图的峰值作为该特征点的主方向，执行完以后转到步骤5.2；

5.2)在特征点所在的尺度空间中以关键点为中心产生一个圆形区域，计算像素梯度，将圆形区域旋转至主方向，对梯度做高斯加权；执行完以后转到步骤5.3；

5.3)将圆形区域分为2×2个扇形子区域，在每个扇形区域内统计8个方向的梯度直方图，执行完以后转到步骤5.4；

5.4)将4个扇形的梯度直方图合并得到32维的描述符然后进行归一化处理，执行完以后转到步骤5.5；

5.5)计算两幅图像中特征向量的夹角，若图像1中特征向量与图像2中特征向量的最小夹角为θ_min，次小夹角为θ_sec，且存在(θ_min/θ_sec)＜t，则对应的特征点为一对匹配点；执行完以后转到步骤5.6.1；反之则特征点不是一对匹配点，继续寻找下一对匹配点，执行完以后转到步骤5.1；

5.6)采用改进的RANSAC算法去除误匹配；

5.6.1)用比值θ_min/θ_sec对初始匹配点对进行排序并构成一个集合，比值最小的特征点对是集合的第一个元素，根据比值大小升序排序，执行完以后转到步骤5.6.2；

5.6.2)在排好序匹配对中选取前S个为一个数据集，在前S-1个数据集中选择3个匹配对，然后选取第S个数据为第4个匹配对，利用四个匹配对计算单应矩阵H；执行完以后转到步骤5.6.3；

5.6.3)将步骤5.6.2中的四组点剔除后，其他点根据单应矩阵H，计算相应投影点；执行完以后转到步骤5.6.4；

5.6.4)计算其他点和投影点之间的误差e，并和误差极限δ比较，e<δ为内点，反之为外点；执行完以后转到步骤5.6.5；

5.6.5)计算内点数目t并和设定内点数目阈值T比较，如t＞T，则内点数目更新为t，执行完以后转到步骤5.6.6；若t≤T，则令S等于S-1，并转到步骤5.6.2；

5.6.6)更新内点数目t后，再次计算单应矩阵H，获取新的内点；执行完以后转到步骤5.6.7；

5.6.7)若迭代次数小于给定的最大迭代次数，返回单应矩阵H，反之模型不符合，剔除匹配点对，结束流程；

即用比值θ_min/θ_sec对初始匹配点对进行排序并构成一个集合，比值最小的特征点对是集合的第一个元素；每次按顺序从集合中选取4对特征点进行变换矩阵的估计。

本发明的有益效果在于：

1)利用仿射不变性特征来表征图像，通过改变相机主光轴方向的经度和纬度这两个角度参数来模拟不同的视场，然后对获取的所有模拟图像对利用某一尺度不变方法进行特征提取与匹配来获得的，将生成许多不同经纬度参数下的模拟图像，足以模拟不同的视点下的仿射形变图像，保证了尺度一定下的仿射不变性。抗视角变换能力强和其他特征相比较具有很强的鲁棒性。当无人机在视觉导航中常常会出现图像之间存在较大视角差异，这种情况下图像匹配性能会急剧下降，对无人机定位和环境构建的准确性进一步产生严重影响。本方法具有很强的鲁棒性，能够减少因图像视角偏差较大所引起的匹配性能降低问题。

2)目前大多数应用在无人机视觉导航中的匹配算法都只满足像素级精度，此方法提出的高精度亚像素图像匹配方法能够高精度的重构出真实环境的纹理等细节信息，提高无人机定位和环境构建的准确性。

3)在亚像素特征匹配时采用简化的描述符进行匹配，提高图像匹配速度。步骤5.4) 只需要在4个扇形区域上统计直方图，统计次数降低，另外算法的描述符维度大幅度降低，在进行描述符相似度的计算时，使用该描述符可以大大降低计算量，提高匹配效率。

附图说明

图1是本发明步骤流程示意图。

图2是本发明图像的小波变换分解示意图。

图3是本发明描述符扇形子区域图。

图4是本发明待匹配基准图。

具体实施方式

下面结合附图和实施例对本发明进一步说明。

实施例如图1所示，具体实施步骤如下

步骤1)：采用哈尔小波变换对输入图像进行压缩，首先执行步骤1.1；

1.1)将输入匹配图像进行小波变换，执行完以后转到步骤1.2；

1.2)运动符号编码器对图像进行编码，生成压缩图像，执行完以后转到步骤1.3；

1.3)对编码完成的图像进行解码处理，执行完以后转到步骤1.4；

1.4)进行小波系数的反量化，执行完以后转到步骤1.5；

1.5)生成最终的重构图像，执行完以后转到步骤2.1；

假设原始图像的大小为2M×2N，记为I(x,y)，则经过哈尔小波变换可将原始图像I(x,y)分为四幅图，如图2所示。图2中变换后的四幅图像保留了原始图像的主要结构。 A为变换后的近似图像、B为水平细节图像、C为垂直细节图像、D为对角线细节图像。通过在近似图像A中进行图像特征提取与匹配能够有效降低计算量，提高快速性。

步骤2)：对压缩处理后的图像进行仿射变换模拟。

仿射模型含有六个自由度，其中两个平移变量tx、ty，一个相机旋转变量ψ，一个均匀缩放系数s，另外两个变量用于定义相机主光轴方向，称为经度φ和纬度θ。通过改变相机主光轴方向的经度和纬度这两个角度参数来模拟不同的视场，然后对获取的所有模拟图像对利用某一尺度不变方法进行特征提取与匹配来获得的，将生成许多不同经纬度参数下的模拟图像，足以模拟不同的视点下的仿射形变图像，保证了尺度一定下的仿射不变性。

对每一幅图像进行仿射变换模拟由于摄像机位置变化产生的仿射形变。图像的仿射形变主要取决于经纬度的变化，因此在对原始图像做仿射变换时只对摄像机的经纬度进行模拟，为了模拟出尽可能多的形变图像，在对经纬度采样时按照如下规则进行：

仿射模型含有六个自由度，其中两个平移变量tx、ty，一个相机旋转变量ψ，一个均匀缩放系数s，另外两个变量用于定义相机主光轴方向，称为经度φ和纬度θ。

2.1)对每一幅图像进行仿射变换模拟由于摄像机位置变化产生的仿射形变。图像的仿射形变主要取决于经纬度的变化，因此在对原始图像做仿射变换时只对摄像机的经纬度进行模拟，为了模拟出尽可能多的形变图像，在对经纬度采样时按照如下规则进行：

①纬度角θ和倾斜程度t满足一一对应关系，即θ＝arccos(1/t)，t的采样值选为 t＝1,a,a²,...,aⁿ。指数n为仿射变换的模拟次数，该值越大则模拟出的仿射变换的情况越多。

②经度角φ通过一组采样数据来实现模拟：0,b/t,...,kb/t(kb/t＜180°)，即一个t对应多个经度角。选择b为角度数，k是一个满足

的最大整数。执行完以后转到步骤3.1；

步骤3)：对每一幅模拟图像提取边缘特征点。

3.1)采用中值滤波法对图像进行噪声去除，以除去孤立的点。执行完以后转到步骤3.2.1；

3.2)对图像进行边缘提取和边缘连接。

3.2.1)采用SUSAN算法提取图像边缘轮廓,执行完以后转到步骤3.2.2)；

边缘轮廓是图像最基本的特征，它蕴含了丰富的内在信息。SUSAN算子是一个简单有效的边缘检测算法，具有较高的算法效率，而且抗噪声能力强。采用SUSAN算法提取图像边缘轮廓。事实上，有噪声影响或遮挡的情况下，得到的边缘像素常常存在孤立或划分成小段连续的情况。

3.2.2)有将断点处的边缘像素连接起来。执行完以后转到步骤3.3.1；

连接准则为：灰度的相似性距离的连贯性。包含下面所述的①和②：

①点的灰度值相近。可用梯度的大小和方向度量像素间的灰度的相似性。边缘点(x1，y1)和(x2，y2)的梯度幅度G₁、G₂和方向θ₁、θ₂分别满足表达式(1)和式(2)。

|||G₁|-|G₂|||≤G₀ (1)

|||θ₁|-|θ₂|||≤θ₀ (2)

其中，G₀和θ₀分别为梯度幅度阈值和方向角度阈值；

②距离的连贯性。像素点(x2，y2)和像素点(x1，y1)连接另一必要条件是点 (x2，y2)在点(x1，y1)的邻域内；

3.3)在图像上提取图像的特征点：

3.3.1)判断边缘轮廓线是否封闭，若封闭，根据式(3)计算几何重心，并对边缘上的边缘点作标记，以免重复计算，行完以后转到步骤3.3.3)；若非封闭，转步骤3.3.2；

对于一个封闭的的不规则平面图形，设边缘点的坐标为(x_i，y_i)，1≤i≤n，n为边缘点的个数，其几何中心坐标(x₀，y₀)计算如下：

3.3.2)若不是封闭曲线，首先进行分段，然后求每一段的几何重心，并对边缘点作标记。执行完以后转到步骤3.3.3；

3.3.3)以几何重心为极点，对边缘点进行极坐标化。执行完以后转到步骤3.3.4；

为了方便计算边缘点上各点距离几何重心的长度及位置，需要将边缘点极坐标化。以几何重心为极点，转换公式见式(4)，极坐标化后的边缘点形成了极径-幅角的曲线：

3.3.4)在极径-极角的曲线上取极大(小)值点。执行完以后转到步骤3.3.5；

极坐标化使问题得到简化，方便在局部寻找极值点，然后再进一步确定特征点。下面分别介绍极值点和特征点的抽取过程。

①极值点的确定若幅角区间[θ₁，θ₂)内的极大值点用P_max表示，则P_max可以用下式描述：

②极值点的确定若幅角区间[θ₁，θ₂)内的极大值点用P_min表示，则P_min可以用下式描述

3.3.5)根据局部非极大(小)值抑制原则消除非极大(小)值点，剩余的点就是特征点。执行完以后转到步骤3.4)；

最终的特征点是按局部区域内非极大(小)值抑制的原则获取的。因此，每一区间的极大(小)值点和相邻区间内的极大(小)值点进行比较，若该极大(小)值点的极径大(小)于相邻两个区间极大(小)值点的极径，则认为是特征点。

3.4)对图像中每一个边缘轮廓线重复步骤3.3.1～3.3.5；执行完以后转到步骤4.1；

在模拟后的图像中进行特征提取，就可以得到具有抗视角变换能力的像素级边缘特征。相较于传统的特征提取算法，该算法在图像存在光照变化、尺度变化和较大视角差异时，仍然能够对图像特征进行准确、充分地提取。

步骤4)：亚像素级特征提取

由于边缘提取算法只能获取像素级精度的特征点，意味着整数像素点之间的信息丢失了。因此需要在此基础上进行亚像素级特征提取，从而获得更高的精度。该算法在本该环节中，采用以下步骤实现：

4.1)沿着粗定位边缘点(即上一步骤获取的像素级边缘特征)的梯度方向寻找距离最近的邻近点(并非都是整像素点)，当梯度方向与水平方向夹角θ小于45°时， Li＝1/cos(θ)为距离长度，执行完以后转到步骤4.2；

4.2)对邻近点进行灰度线性插值处理；最后通过邻近点的灰度值计算粗定位特征点的灰度差值。根据高斯亚像素拟合原理，图像边缘点灰度值变化呈高斯分布，高斯曲线的表达式如下：

式中，μ是中心点坐标，σ是高斯函数的标准差，u是粗定位点的灰度差值；执行完以后转到步骤4.3；

4.3)采用最小二乘法拟合并计算出μ和σ，假设粗定位边缘坐标为p(x₀,y₀)，则根据σ和求取得边缘点梯度方向夹角θ可以求出亚像素级边缘点的坐标为 p′(x₀+σcos(θ),y₀+σsin(θ))。执行完以后转到步骤5.1；

步骤5)：亚像素特征匹配

该步骤包括构建特征向量和特征向量的匹配。拟采用以下步骤来实现：

5.1)以亚像素边缘点为特征点(x,y)，按照式(6)计算其梯度和模值。

式中，L是每个特征点所在的模拟图像。在以特征点为中心的邻域窗口内进行采样，采用直方图来统计领域像素的梯度方向，并用高斯函数对直方图进行平滑。梯度直方图描述的范围是0°～360°，每10°一个柱，共36个柱，以直方图的峰值作为该特征点的主方向。如图3所示，执行完以后转到步骤5.2；

5.2)在特征点所在的尺度空间中以关键点为中心产生一个圆形区域，计算像素梯度，将圆形区域旋转至主方向，对梯度做高斯加权。执行完以后转到步骤5.3；

5.3)将圆形区域分为2×2个扇形子区域，在每个扇形区域内统计8个方向的梯度直方图如图3所示。执行完以后转到步骤5.4；

5.4)将4个扇形的梯度直方图合并得到32维的描述符然后进行归一化处理，以降低光照变化对描述符的影响。执行完以后转到步骤5.5；

5.5)计算两幅图像中特征向量的夹角，若图像1中特征向量与图像2中特征向量的最小夹角为θ_min，次小夹角为θ_sec，且存在(θ_min/θ_sec)＜t，则对应的特征点为一对匹配点。执行完以后转到步骤5.6.1；反之则特征点不是一对匹配点，继续寻找下一对匹配点，执行完以后转到步骤5.1。

5.6)采用改进的RANSAC算法去除误匹配。

5.6.1)用比值θ_min/θ_sec对初始匹配点对进行排序并构成一个集合，比值最小的特征点对是集合的第一个元素，根据比值大小升序排序，执行完以后转到步骤5.6.2；

5.6.2)在排好序匹配对中选取前S个为一个数据集，在前S-1个数据集中选择3个匹配对，然后选取第S个数据为第4个匹配对，利用四个匹配对计算单应矩阵H；执行完以后转到步骤5.6.3；

5.6.3)将上述四组点剔除后，其他点根据单应矩阵H，计算相应投影点；执行完以后转到步骤5.6.4；

5.6.4)计算其他点和投影点之间的误差e，并和误差极限δ比较，e<δ为内点，反之为外点；执行完以后转到步骤5.6.5；

5.6.5)计算内点数目t并和设定内点数目阈值T比较，t＞T,内点数目更新为t，执行完以后转到步骤5.6.6；若t≤T，则令S等于S-1，并转到步骤5.6.2；

5.6.6)更新内点数目t后，再次计算单应矩阵H，以及获取新的内点；执行完以后转到步骤5.6.7；

5.6.7)若迭代次数小于给定的最大迭代次数，返回单应矩阵H，反之模型不符合，剔除匹配点对。执行完以后结束流程。

即用比值θ_min/θ_sec对初始匹配点对进行排序并构成一个集合，比值最小的特征点对是集合的第一个元素。每次按顺序从集合中选取4对特征点进行变换矩阵的估计。相较于传统RANSAC算法中随机选取特征点对，该算法能有效避免陷入局部最优，保证了最终图像匹配的准确度。

(二)实施例

为验证实际应用的精确性与高效性，进行实验。实验采用某地区上空的航拍遥感影像，图像分辨率为1280×1920像素，位移图像通过对原始图像平移得到。通过对存在整像素平移关系的两幅图像做10倍下采样的方法，便可得到存在亚像素平移关系的原始图像和待配准图像。

为方便计算，本文只对原始图像如图4所示，进行水平方向的平移操作。依次对原始图像进行步长为20像素的水平平移操作，一共生成10帧位移图像。相对于原始图像，第n帧位移图像的水平位移dx＝1+15×(n-1)像素；再对原始图像和20帧平移图像做10倍下采样，便得到大小均为128×192像素，平移关系依次为向右平移0.1， 1.6，3.1，4.6，6.1，7.6，9.1，10.6，12.1,13.6,15.1,16.6,18.1,19.6,21.1,22.6,24.1,25.6,27.1,28.6 像素的20帧位移图像。

采用数字图像相关法，可以轻易求得平移关系依次为0.1，1.6，3.1，4.6，6.1，7.6，9.1，10.6，12.1，13.6，15.1,16.6,18.1,19.6,21.1,22.6,24.1,25.6,27.1,28.6像素的20帧位移图像与原始图像的整像素配准位置依次为0，2，3，5，6，8，9，11，12，14， 15，17，18，20，21，23，24，26，27，29像素。在事先获得整像素配准位置的前提下，我们分别采用曲面拟合法和梯度法对原始图像和20帧位移图像进行亚像素配准，配准结果见表1。显然，在使用梯度法求解亚像素配准位置时，位移图像正是以整像素配准位置作为其一阶泰勒级数的展开位置。

配准误差如下表：

平移关系	本算法	绝对误差
			0.1	0.08	-0.02
1.6	1.61	0.01
			3.1	3.12	0.02
4.6	4.60	0.00
			6.1	6.10	0.00
7.6	7.59	-0.01
			9.1	9.11	0.01
10.6	10.58	-0.02
			12.1	12.09	0.01
13.6	13.60	0.00
			15.1	15.08	-0.02
16.6	16.58	-0.02
			18.1	18.12	0.02
19.6	19.59	-0.01
			21.1	21.11	0.01
22.6	22.62	0.02
			24.1	24.10	0.00
25.6	25.58	-0.02
			27.1	27.09	-0.01
28.6	28.61	0.01

匹配定位误差控制在0.02亚像素级以内。

将此方法应用在无人机景象匹配中，当无人机飞到预定区域时，机载图像传感器实时获取当地景象作为实测图像送到机载计算机中和基准图像进行匹配，由于基准图的地理坐标位置是已知的，因此，根据与实测图的匹配位置，即可以确定当前飞行器的准确位置。由于选用基于亚像素的匹配，在实时定位中将景象匹配精度提高到亚像素级。

Claims (1)

1.一种亚像素级图像匹配导航定位方法，其特征在于包括下述步骤：

步骤1：采用哈尔小波变换对输入图像进行压缩，首先执行步骤1.1；

1.1)将输入匹配图像进行小波变换，执行完以后转到步骤1.2；

1.2)运动符号编码器对图像进行编码，生成压缩图像，执行完以后转到步骤1.3)；

1.3)对编码完成的图像进行解码处理，执行完以后转到步骤1.4；

1.4)进行小波系数的反量化，执行完以后转到步骤1.5；

1.5)生成最终的重构图像，执行完以后转到步骤2.1；

假设原始图像的大小为2M×2N，记为I(x,y)，则经过哈尔小波变换可将原始图像I(x,y)分为四幅图，A为变换后的近似图像、B为水平细节图像、C为垂直细节图像、D为对角线细节图像，通过在近似图像A中进行图像特征提取与匹配能够有效降低计算量，提高快速性；

步骤2：对压缩处理后的图像进行仿射变换模拟；

仿射模型含有六个自由度，其中两个平移变量tx、ty，一个相机旋转变量ψ，一个均匀缩放系数s，另外两个变量用于定义相机主光轴方向，称为经度φ和纬度θ；

2.1)对每一幅图像进行仿射变换,模拟由于摄像机位置变化产生的仿射形变；图像的仿射形变取决于经纬度的变化，因此在对原始图像做仿射变换时只对摄像机的经纬度进行模拟，在对经纬度采样时按照如下规则进行：

①纬度角θ和倾斜程度t满足一一对应关系，即θ＝arccos(1/t)，t的采样值选为t＝1,a,a²,...,aⁿ，指数n为仿射变换的模拟次数，该值越大则模拟出的仿射变换的情况越多；

②经度角φ通过一组采样数据来实现模拟：0,b/t,...,kb/t(kb/t＜180°)，即一个t对应多个经度角，选择b为角度数，k是一个满足

的最大整数，执行完以后转到步骤3.1；

步骤3)：对每一幅模拟图像提取边缘特征点；

3.1)采用中值滤波法对图像进行噪声去除，以除去孤立的点，执行完以后转到步骤3.2；

3.2)采用SUSAN算法提取图像边缘轮廓,执行完以后转到步骤3.3；

3.3)将边缘像素连接起来，执行完以后转到步骤3.4；

3.4)在图像边缘上提取图像的特征点，得到粗定位边缘点，执行完以后转到步骤4.1；

步骤4)：亚像素级特征提取

4.1)沿着粗定位边缘点的梯度方向寻找距离最近的邻近点，当梯度方向与水平方向夹角θ小于45°时，Li＝1/cos(θ)为距离长度，执行完以后转到步骤4.2；

4.2)对邻近点进行灰度线性插值处理；最后通过邻近点的灰度值计算粗定位特征点的灰度差值；根据高斯亚像素拟合原理，图像边缘点灰度值变化呈高斯分布，高斯曲线的表达式如下：

式中，μ是中心点坐标，σ是高斯函数的标准差，u是粗定位点的灰度差值；执行完以后转到步骤4.3；

4.3)采用最小二乘法拟合并计算出μ和σ，假设粗定位边缘坐标为p(x₀,y₀)，则根据σ和求取得边缘点梯度方向夹角θ求出亚像素级边缘点的坐标为p′(x₀+σcos(θ),y₀+σsin(θ))；执行完以后转到步骤5.1；

步骤5)：亚像素特征匹配

包括构建特征向量和特征向量的匹配，采用以下步骤来实现：

5.1)以亚像素边缘点为特征点(x,y)，按照下式计算其梯度和模值：

式中，L是每个特征点所在的模拟图像；在以特征点为中心的邻域窗口内进行采样，采用直方图统计领域像素的梯度方向，并用高斯函数对直方图进行平滑；梯度直方图描述的范围是0°～360°，每10°一个柱，共36个柱，以直方图的峰值作为该特征点的主方向，执行完以后转到步骤5.2；

5.2)在特征点所在的尺度空间中以关键点为中心产生一个圆形区域，计算像素梯度，将圆形区域旋转至主方向，对梯度做高斯加权；执行完以后转到步骤5.3；

5.3)将圆形区域分为2×2个扇形子区域，在每个扇形区域内统计8个方向的梯度直方图，执行完以后转到步骤5.4；

5.4)将4个扇形的梯度直方图合并得到32维的描述符然后进行归一化处理，执行完以后转到步骤5.5；

5.5)计算两幅图像中特征向量的夹角，若图像1中特征向量与图像2中特征向量的最小夹角为θ_min，次小夹角为θ_sec，且存在(θ_min/θ_sec)＜t，则对应的特征点为一对匹配点；执行完以后转到步骤5.6.1；反之则特征点不是一对匹配点，继续寻找下一对匹配点，执行完以后转到步骤5.1；

5.6)采用改进的RANSAC算法去除误匹配；

5.6.1)用比值θ_min/θ_sec对初始匹配点对进行排序并构成一个集合，比值最小的特征点对是集合的第一个元素，根据比值大小升序排序，执行完以后转到步骤5.6.2；

5.6.2)在排好序匹配对中选取前S个为一个数据集，在前S-1个数据集中选择3个匹配对，然后选取第S个数据为第4个匹配对，利用四个匹配对计算单应矩阵H；执行完以后转到步骤5.6.3；

5.6.3)将步骤5.6.2中的四组点剔除后，其他点根据单应矩阵H，计算相应投影点；执行完以后转到步骤5.6.4；

5.6.4)计算其他点和投影点之间的误差e，并和误差极限δ比较，e<δ为内点，反之为外点；执行完以后转到步骤5.6.5；

5.6.5)计算内点数目t并和设定内点数目阈值T比较，如t＞T，则内点数目更新为t，执行完以后转到步骤5.6.6；若t≤T，则令S等于S-1，并转到步骤5.6.2；

5.6.6)更新内点数目t后，再次计算单应矩阵H，获取新的内点；执行完以后转到步骤5.6.7；

5.6.7)若迭代次数小于给定的最大迭代次数，返回单应矩阵H，反之模型不符合，剔除匹配点对，结束流程；

即用比值θ_min/θ_sec对初始匹配点对进行排序并构成一个集合，比值最小的特征点对是集合的第一个元素；每次按顺序从集合中选取4对特征点进行变换矩阵的估计。

Images