基于各向异性梯度尺度空间的遥感图像配准方法与流程

本发明属于图像处理技术领域，更进一步涉及对遥感图像进行配准处理技术领域中的一种基于各向异性梯度尺度空间的遥感图像配准方法。本发明可应用于对多光谱、多源不同时刻获取的遥感图像的配准。

背景技术：

图像配准是指在不同时段，对同一场景从不同视角使用相同或者不同传感器拍摄的有重叠区域图像进行几何矫正的过程。图像配准技术是近年来发展迅速的图像处理技术之一。图像配准技术已经广泛用于各个领域，如航空航天技术、图像镶嵌、地理信息系统、图像融合、三维重建、目标识别和变化检测等等。目前随着遥感技术的发展，由于传感器的不同物理特性所产生的遥感图像不断增多，因此综合利用各种图像进行数据提取和分析已成为遥感领域的一个重要手段。同时由于各种传感器之间物理特性和成像方式的不同，在数据应用和数据融合时不同几何特性和不同分辨率的图像间必须进行严格的配准。

目前遥感图像配准主要分为两类：基于区域灰度的配准方法和基于特征的配准方法。其中，基于区域灰度的常用图像配准方法有：互相关(CC)、不变矩、基于FFT的相位相关法和互信息(MI)等。但是基于灰度的图像配准方法存在着以下缺点：①对图像灰度变化比较敏感，尤其是非线性的光照变化，这大大降低了算法的性能；②计算复杂度高；③对目标的旋转、形变和遮挡比较敏感。为了克服其存在的缺点，人们提出了基于特征的图像配准方法。基于图像特征的方法首先从图像中提取边缘、角点、轮廓和区域中心等特征，然后使用特征之间的相关决定图像的最优对齐。这些显著特征可以大大压缩图像的信息量，使得计算量减小，速度更快，而且对图像的灰度变化具有鲁棒性。目前来说点特征是最为常用且效率较高的一种方法。二维图像中的角点、拐点、交叉点等是图像的明显特征。这些点具有平移、旋转和缩放不变性，几乎不受光照条件的影响。它只需用图像中大约0.05％的像素点就可以表示整幅图像的数据信息，其信息含量高、计算速度快，使实时处理成为可能。因此，基于特征的图像配准方法在遥感图像配准领域也得到了广泛的应用。

经典的SIFT算法具有对可见光图像旋转、缩放、部分的光照和视角变换保持不变性的优点。当SIFT方法检测特征点数目增加时，生成特征点主方向和生成特征点描述子步骤花费时间迅速增加，限制了SIFT算法在实际中的应用。由于SIFT算法的优点，国内外学者提出了许多基于SIFT的提升算法。然而当使用基于SIFT的方法匹配遥感图像时出现了许多不正确的匹配点，正确匹配率(CMT)迅速下降。原因在于拍摄时间、光谱和获取传感器的不同导致遥感图像对同一区域的像素亮度显著不同，同时像素对之间的亮度映射可能是线性、非线性，甚至会出现亮度反转情况。

Sedaghat在其发表的论文“Uniform Robust Scale-Invariant Feature Matching for Optical Remote Sensing Images”(《IEEE Transactions on Geoscience and Remote Sensing》,2011,pages 4516-4527)中提出了一种UR-SIFT算法。该算法较好的控制了尺度空间提取的特征点的数量、质量和分布，提高了具有局部变换的部分遥感图像的配准精度，但是该算法仍然存在的不足之处是，不能精确的配准图像亮度有较大非线性变化的遥感图像对。

Gong在其发表的论文“A Novel Coarse-to-Fine Scheme for Automatic Image Registration Based on SIFT and Mutual Information”(《IEEE Transactions on Geoscience and Remote Sensing》,2014,pages 4328-4338)中公开了一种由粗到精的配准方法。该方法先使用SIFT算法获得图像对之间的初始变换关系，然后结合初始变换关系和互信息方法获得更加精确的图像配准。该方法仍然存在的不足之处是，由于该方法是基于互信息的方法，所以计算复杂度高，而且配准精度依懒于sift算法获得的初始解，当图像亮度存在较大非线性变化的情况下，sift算法并不能获得很好的初始解，依赖sift算法获得的初始解的互信息方法也不能实现精确的配准。

西安电子科技大学在其申请的专利“基于互信息选块和sift特征的遥感图像配准方法”(专利申请号：CN201410379927，公开号：CN104200461A)中提出了一种基于sift和互信息的遥感图像配准方法。该方法的实现过程为：随机从参考遥感图像和待配准遥感图像中选取图像对，计算每一对图像的互信息；对互信息进行降序排列；选取前n个互信息较大的子图像对；对选取的前n个互信息较大的图像对提取sift特征，并进行粗匹配；删除错误匹配点，进行细匹配；计算配准参数和互信息值；选取互信息最大的配准参数作为最终的配准参数。该方法虽然能够加快图像配准的速度，但是，该方法仍然存在的不足之处是，在特征点生成、特征点主方向生成和特征点的描述符生成阶段仍采用sift方法，在遥感图像对的亮度非线性变化较大的情况下正确匹配率迅速下降，同时该方法采用传统的最近邻和次近邻距离比匹配准则，当遥感图像中存在较多重复特征情况下正确配对的特征点数目也迅速下降。

技术实现要素：

本发明的目的是克服上述已有技术的不足，提出一种基于各向异性梯度尺度空间的遥感图像配准方法，以提高特征匹配的准确度，实现对图像亮度有较大非线性变化的遥感图像对的配准。

本发明实现上述目的的思路是：首先根据非线性扩散原理构建参考遥感图像和待配准遥感图像的各向异性尺度空间，再对参考遥感图像和待配准遥感图像的各向异性尺度空间的梯度幅度图像上使用Harris算子进行角点检测，接着生成点特征的对图像亮度的非线性变化有较强鲁棒性的描述子，最后使用改进的特征点匹配准则进行特征点匹配，得到最后的遥感图像对配准结果图。

本发明的步骤包括如下：

(1)输入参考遥感图像和待配准遥感图像；

(2)构造各向异性扩散的尺度空间：

(2a)计算参考遥感图像和待配准遥感图像的各向异性尺度空间各层的尺度值；

(2b)将尺度值转换到时间度量值；

(2c)对输入的参考遥感图像和待配准遥感图像，采用标准差为σ₀的高斯滤波，得到参考遥感图像和待配准遥感图像各向异性尺度空间第0层图像；

(2d)将各向异性尺度空间层的序号i初始化为零；

(2e)按照下式，分别计算参考遥感图像和待配准遥感图像各向异性尺度空间的第i层图像的扩散系数矩阵：

其中，cⁱ表示参考遥感图像或者待配准遥感图像各向异性尺度空间的第i层图像的扩散系数矩阵，表示参考遥感图像或者待配准遥感图像各向异性尺度空间的第i层图像Iⁱ应用标准差为1的高斯滤波后的图像，表示高斯滤波后的图像的梯度幅度，|·|表示取模操作，K表示对比度因子，K的取值为梯度幅度的统计直方图70％百分位的梯度幅度值；

(2f)按照下式，对参考遥感图像和待配准遥感图像各向异性尺度空间的第i层图像做行扩散：

其中，表示沿参考遥感图像或者待配准遥感图像各向异性尺度空间的第i层图像Iⁱ的行扩散后的图像，I₁表示行数和列数均与参考遥感图像或者待配准遥感图像的列数相同的单位方阵，t_i和t_i+1分别表示参考遥感图像或者待配准遥感图像各向异性尺度空间的第i层和第i+1层的时间度量值，A₁(Iⁱ)表示编码参考遥感图像或者待配准遥感图像各向异性尺度空间的第i层图像扩散系数cⁱ的矩阵，Iⁱ表示参考遥感图像或者待配准遥感图像各向异性尺度空间的第i层图像，(·)^-1表示逆矩阵操作；

(2g)按照下式，对参考遥感图像和待配准遥感图像各向异性尺度空间的第i层图像做列扩散：

其中，表示沿参考遥感图像或者待配准遥感图像各向异性尺度空间的第i层图像Iⁱ的列扩散后的图像，I₂表示行数和列数均与参考遥感图像或者待配准遥感图像的行数相同的单位方阵；

(2h)按照下式，计算参考遥感图像和待配准遥感图像各向异性尺度空间的第i+1层图像：

其中，Iⁱ⁺¹表示参考遥感图像或者待配准遥感图像各向异性尺度空间的第i+1层图像；

(2i)判断i≥N-2是否成立，若成立，得到参考遥感图像和待配准遥感图像的各向异性尺度空间，否则，令i＝i+1，执行步骤(2e)；其中，N表示参考遥感图像和待配准遥感图像各向异性尺度空间层的总数；

(3)计算梯度幅度图像：

使用索贝尔算子Sobel，计算参考遥感图像和待配准遥感图像各向异性尺度空间的梯度幅度图像；

(4)计算梯度幅度图像的差分：

(4a)按照下式，使用索贝尔算子Sobel，计算参考遥感图像和待配准遥感图像各向异性尺度空间的梯度幅度图像的x轴正方向的差分：

其中，表示参考遥感图像或者待配准遥感图像各向异性尺度空间的梯度幅度图像GIⁿ沿水平向右的x轴正方向的差分，GIⁿ表示参考遥感图像或者待配准遥感图像各向异性尺度空间的梯度幅度图像的第n层图像，表示相关操作，n＝0,1,...,N-1，N表示参考遥感图像和待配准遥感图像各向异性尺度空间层的总数；

(4b)按照下式，使用索贝尔算子Sobel，计算参考遥感图像和待配准遥感图像各向异性尺度空间的梯度幅度图像的y轴正方向的差分：

其中，表示参考遥感图像或者待配准遥感图像各向异性尺度空间的梯度幅度图像GIⁿ沿竖直向下的y轴正方向的差分；

(5)按照下式，计算梯度幅度图像的梯度幅度：

其中，GGIⁿ表示参考遥感图像或者待配准遥感图像各向异性尺度空间的梯度幅度图像的第n层图像GIⁿ的梯度幅度，表示开根号操作；

(6)按照下式，计算梯度幅度图像的梯度角度：

其中，AGGIⁿ表示参考遥感图像或者待配准遥感图像各向异性尺度空间的梯度幅度图像的第n层图像GIⁿ的梯度角度，arctan(·)表示四象限反正切操作；

(7)检测特征点：

使用Harris角点检测算子，在参考遥感图像和待配准遥感图像的各向异性尺度空间的梯度幅度图像上检测特征点，得到特征点集合：

其中，C_IR表示在参考遥感图像的各向异性尺度空间的梯度幅度图像上检测的特征点集合，C_IS表示在待配准遥感图像的各向异性尺度空间的梯度幅度图像上检测的特征点集合，R表示在参考遥感图像的各向异性尺度空间的梯度幅度图像上检测的特征点的总数，S表示在待配准遥感图像的各向异性尺度空间的梯度幅度图像上检测的特征点的总数；

(8)生成特征点主方向：

(8a)将水平方向角度在[0,2π]内36等分；

(8b)由特征点集合的半径为6σ_n的圆形邻域区域内像素的梯度方向AGGIⁿ(X)，确定圆形区域内像素所在的等分角度范围，其中，σ_n表示参考遥感图像或者待配准遥感图像各向异性尺度空间的第n层的尺度值，AGGIⁿ(X)表示特征点周围圆形邻域内位置坐标为X的像素的梯度角度；

(8c)分别累加位于每个等分角度范围内的所有像素的梯度幅度值GGIⁿ(X)，形成梯度方向直方图，其中，GGIⁿ(X)表示特征点周围圆形邻域内位置坐标为X的像素的梯度幅度；

(8d)将梯度方向直方图中大于最大值0.8倍的数值对应的梯度方向，作为特征点的主方向；

(9)生成特征点描述子：

(9a)将半径为ρ的特征点的圆形邻域逆时针旋转θ_m度，圆形邻域沿半径方向将ρ分为3个区间，内圆邻域半径为0.25ρ，中间圆邻域半径为0.73ρ，外圆邻域半径为ρ，圆形邻域沿角度方向将[0,2π]等分为8个区间，内圆作为一个整体，特征点周围圆形邻域被划分为了17个不同位置的面积相等的子区域，其中，ρ的取值为12σ_n，θ_m表示该特征点的主方向；

(9b)将特征点圆形邻域内像素的笛卡尔坐标转换为对数极坐标，对数极坐标角度方向水平向右，在[0,2π]的范围内等分为8个区间，对数极坐标对数长度方向竖直向下，在范围内非等分为3个区间，其中，ρ表示特征点周围圆形邻域的半径；

(9c)对数极坐标网格中每个子区域内的所有像素根据其梯度幅度GGIⁿ(X)和梯度方向AGGIⁿ(X)计算梯度方向直方图，每一个子区域形成了一个8维的梯度方向向量，依次拼接17子区域的梯度方向向量就形成了一个136维的特征点的描述子，其中，梯度方向直方图的角度在[0,2π]的范围内等分为8个区间；

(9d)采用步骤(9a)、步骤(9b)、步骤(9c)的相同方法，生成特征点集合C_IR的描述子集合D_IR和特征点集合C_IS的描述子集合D_IS，其中，D_IR表示参考遥感图像的各向异性尺度空间的梯度幅度图像上检测到的特征点的描述子集合，D_IS表示待配准遥感图像的各向异性尺度空间的梯度幅度图像上检测到的特征点的描述子集合，L₁表示参考遥感图像的各向异性尺度空间的梯度幅度图像上检测到的特征点的描述子总数，L₂表示待配准遥感图像的各向异性尺度空间的梯度幅度图像上检测到的特征点的描述子总数；

(10)特征点匹配：

(10a)将待配准遥感图像特征点的描述子序号k初始化为1；

(10b)按照下式，计算待配准遥感图像特征点的描述子与参考遥感图像每一个特征点的描述子的欧式距离组成的向量：

其中，d_k表示待配准遥感图像特征点的描述子与参考遥感图像的每一个特征点的描述子的欧式距离组成的向量，d_k,j表示配准遥感图像特征点的描述子与参考遥感图像特征点的描述子的欧式距离，表示待配准遥感图像特征点的第k个描述子，表示参考遥感图像特征点的第j个描述子，j＝1,2,...,L₁，L₁表示参考遥感图像特征点描述子的个数，||·||表示取范数操作；

(10c)对向量d_k中的L₁个元素从小到大排序得到向量s_k；

(10d)如果则描述子对应的特征点和与描述子的欧式距离为s_k,0的参考遥感图像的描述子对应的特征点匹配；如果且则描述子对应的特征点和与描述子的欧式距离为d_k的前E个元素的参考遥感图像的描述子对应的特征点匹配；如果则描述子对应的特征点不和参考遥感图像的任何特征点匹配；

其中，s_k,0表示向量s_k的第一个元素，s_k,1表示向量s_k的第二个元素，T_L和T_H表示最近邻和次近邻距离比阈值，E的取值范围为[2,3]中的整数；

(10e)判断k≥L₂是否成立，若成立，得到待配准遥感图像和参考遥感图像的特征点匹配关系C_IR,IS，执行步骤(11)；否则，令k＝k+1，执行步骤(10b)；

其中，C_IR,IS表示参考遥感图像的各向异性尺度空间的梯度幅度图像上检测到的特征点和在待配准遥感图像的各向异性尺度空间的梯度幅度图像上检测到的特征点的配对关系集合，L₃表示特征点匹配阶段匹配的点对总数；

(11)删除错误匹配的特征点对：

使用随机抽样一致RANSAC算法，删除错误的匹配点对，得到待配准遥感图像到参考遥感图像的几何变换关系；

(12)配准参考图像和待配准图像：

根据几何变换关系，使用双线性插值方法配准参考遥感图像和待配准遥感图像。

本发明与现有技术相比具有如下优点：

第一，由于本发明在各向异性尺度空间的梯度幅度图像上进行特征点检测，使用各向异性尺度空间的梯度幅度图像的梯度幅度和梯度角度计算特征点的主方向和特征点的描述符，克服了现有技术不能应对遥感图像对之间的亮度非线性变化较大的问题，使得本发明提高了遥感图像对的正确匹配率。

第二，由于本发明在生成特征点的主方向和生成特征点描述子阶段并未使用高斯权重窗口，克服了现有技术在生成特征点主方向和生成特征点描述子阶段计算复杂度高的问题，使得本发明提高了算法的运行效率。

第三，由于本发明在特征点匹配阶段对最近邻和次近邻距离比设置了两个阈值，克服了现有技术使用一个阈值进行特征点匹配时正确配对的特征点数目少的问题，使得本发明提高了正确配对的特征点数目和配准精度。

附图说明

图1为本发明的流程图；

图2为现有技术的特征点周围圆形邻域区域划分和对数极坐标网格示意图；

图3为本发明步骤11的随机采样一致RANSAC算法流程图；

图4为本发明仿真实验的第一组多光谱遥感图像对和配准结果图；

图5为本发明仿真实验的第二组多光谱遥感图像对和配准结果图；

图6为本发明仿真实验的多源遥感图像对和配准结果图。

具体实施方式

下面结合附图对本发明做进一步的详细描述。

参照附图1，本发明的步骤如下。

步骤1，输入遥感图像对。

输入参考遥感图像和待配准遥感图像。

步骤2，构造各向异性扩散的尺度空间。

(2a)计算参考遥感图像和待配准遥感图像的各向异性尺度空间各层的尺度值：

σ_e(m,u)＝σ₀2^m+u/L

其中，σ_e表示参考遥感图像或者待配准遥感图像各向异性尺度空间的第e层的尺度值，σ₀表示参考遥感图像或者和待配准遥感图像的基准尺度值，由不同σ_e组成的各向异性尺度空间共G组，每组L层，m表示各向异性尺度空间的组索引，m＝0,1,...,G-1，u表示每组内层的索引，u＝0,1,...,L-1，e＝0,1,2,...,N-1。

(2b)将尺度值转换到时间度量值：

其中，t_e表示参考遥感图像或者待配准遥感图像各向异性尺度空间的第e层的时间度量值，σ_e表示参考遥感图像或者待配准遥感图像各向异性尺度空间的第e层的尺度，e＝0,1,2,...,N-1，N表示参考遥感图像和待配准遥感图像各向异性尺度空间层的总数，m表示各向异性尺度空间的组索引，m＝0,1,...,G-1，G表示参考遥感图像和待配准遥感图像各向异性尺度空间组的总数，u表示每组内层的索引，u＝0,1,...,L-1，L表示参考遥感图像和待配准遥感图像各向异性尺度空间每组层的总数。

(2c)对输入的参考遥感图像和待配准遥感图像，采用标准差为σ₀的高斯滤波，得到参考遥感图像和待配准遥感图像各向异性尺度空间第0层图像。

(2d)将各向异性尺度空间层的序号i初始化为零。

(2e)按照下式，分别计算参考遥感图像和待配准遥感图像各向异性尺度空间的第i层图像的扩散系数矩阵：

其中，cⁱ表示参考遥感图像或者待配准遥感图像各向异性尺度空间的第i层图像的扩散系数矩阵，表示参考遥感图像或者待配准遥感图像各向异性尺度空间的第i层图像Iⁱ应用标准差为1的高斯滤波后的图像，表示高斯滤波后的图像的梯度幅度，|·|表示取模操作，K表示对比度因子，K的取值为梯度幅度的统计直方图70％百分位的梯度幅度值。

(2f)对参考遥感图像和待配准遥感图像各向异性尺度空间的第i层图像Iⁱ的行做一维扩散：

第一步：计算编码各向异性尺度空间的第i层图像Iⁱ的扩散系数矩阵cⁱ第h行的矩阵A_1,h(Iⁱ)。

其中，A_1,h(Iⁱ)表示编码参考遥感图像或者待配准遥感图像各向异性尺度空间的第i层图像Iⁱ的扩散系数矩阵cⁱ第h行的矩阵，cⁱ表示参考遥感图像或者待配准遥感图像各向异性尺度空间的第i层图像Iⁱ的扩散系数矩阵，表示扩散系数矩阵cⁱ的第h行，第J列的元素，J＝0,1,...,Q-1，表示矩阵cⁱ第h行，第w列的元素的左右2邻域的列坐标的集合，w＝0,1,...,Q-1，Q表示参考遥感图像或者待配准遥感图像的列数。

第二步：按照下式，采用Thomas算法解各向异性行扩散方程。

其中，表示沿参考遥感图像或者待配准遥感图像的各向异性尺度空间的第i层图像Iⁱ的第h行扩散后的结果，I是大小和输入参考遥感图像或者待配准遥感图像的列数相同的单位矩阵，t_i和t_i+1分别是参考遥感图像或者待配准遥感图像各向异性尺度空间的第i层和第i+1层的时间度量值，A_1,h(Iⁱ)表示编码参考遥感图像或者待配准遥感图像各向异性尺度空间的第i层图像Iⁱ的扩散系数矩阵cⁱ第h行的矩阵，表示参考遥感图像或者待配准遥感图像的各向异性尺度空间的第i层图像Iⁱ的第h行，(·)^-1表示逆矩阵操作。

第三步：对于h＝0,1,...,P-1，重复执行本步骤中的第一步和第二步，得到其中，表示沿参考遥感图像或者待配准遥感图像各向异性尺度空间的第i层图像Iⁱ的行做一维扩散后的结果，P表示参考遥感图像或者待配准遥感图像的行数。

(2g)对参考遥感图像和待配准遥感图像各向异性尺度空间的第i层图像Iⁱ的列做一维扩散：

第一步：计算编码各向异性尺度空间的第i层图像Iⁱ的扩散系数矩阵cⁱ第v列的矩阵A_2,v(Iⁱ)。

其中，A_2,v(Iⁱ)表示编码参考遥感图像或者待配准遥感图像各向异性尺度空间的第i层图像Iⁱ的扩散系数矩阵cⁱ第v列的矩阵，cⁱ表示参考遥感图像或者待配准遥感图像各向异性尺度空间的第i层图像Iⁱ的扩散系数矩阵，表示矩阵cⁱ的第M行，第v列的元素，M＝0,1,...,P-1，表示矩阵cⁱ第v列，第f行的元素的上下2邻域的行坐标的集合，f＝0,1,...,P-1，P表示参考遥感图像或者待配准遥感图像的行数。

第二步：按照下式，采用Thomas算法解各向异性列扩散方程。

其中，表示沿参考遥感图像或者待配准遥感图像的各向异性尺度空间的第i层图像Iⁱ的第v列扩散后的结果，I是大小和输入参考遥感图像或者待配准遥感图像行数相同的单位矩阵，t_i和t_i+1分别是各向异性尺度空间的第i层和第i+1层的时间度量值，A_2,v(Iⁱ)表示编码参考遥感图像或者待配准遥感图像各向异性尺度空间的第i层图像Iⁱ的扩散系数矩阵cⁱ第v列的矩阵，表示参考遥感图像或者待配准遥感图像的各向异性尺度空间的第i层图像Iⁱ的第v列，(·)^-1表示逆矩阵操作。

第三步：对于v＝0,1,...,Q-1，重复执行本步骤中的第一步和第二步，得到表示沿参考遥感图像或者待配准遥感图像各向异性尺度空间的第i层图像Iⁱ的列做一维扩散后的结果，Q表示参考遥感图像或者待配准遥感图像的列数。

(2h)按照下式，计算参考遥感图像和待配准遥感图像各向异性尺度空间的第i+1层图像：

其中，Iⁱ⁺¹表示参考遥感图像或者待配准遥感图像各向异性尺度空间的第i+1层图像，表示参考遥感图像或者待配准遥感图像各向异性尺度空间的第i层图像经过行扩散后的图像，表示参考遥感图像或者待配准遥感图像各向异性尺度空间的第i层图像经过列扩散后的图像。

(2i)判断i≥N-2是否成立，若成立，结束迭代，得到参考遥感图像和待配准遥感图像各向异性尺度空间，否则，令i＝i+1，执行步骤(2d)；其中，i表示各向异性尺度空间层的序号，N表示参考遥感图像和待配准遥感图像各向异性尺度空间层的总数。

步骤3，计算梯度幅度图像。

使用索贝尔算子Sobel，计算参考遥感图像和待配准遥感图像各向异性尺度空间的梯度幅度图。

按照下式，计算参考遥感图像和待配准遥感图像各向异性尺度空间的各层图像沿x轴正方向的差分：

其中，Iⁿ表示参考遥感图像或者待配准遥感图像各向异性尺度空间的第n层图像，表示坐标原点在左上角的图像Iⁿ沿水平向右x轴正方向的差分，表示相关操作，n＝0,1,...,N-1，N表示参考遥感图像和待配准遥感图像各向异性尺度空间层的总数。

按照下式，计算参考遥感图像和待配准遥感图像各向异性尺度空间的各层图像沿y轴正方向的差分：

其中，表示坐标原点在左上角的图像Iⁿ沿竖直向下y轴正方向的差分。

按照下式，计算参考遥感图像和待配准遥感图像各向异性尺度空间各层图像的梯度幅度：

其中，GIⁿ表示参考遥感图像或者待配准遥感图像各向异性尺度空间的第n层图像Iⁿ的梯度幅度，表示参考遥感图像或者待配准遥感图像各向异性尺度空间的第n层图像Iⁿ沿水平向右x轴正方向的差分，表示参考遥感图像或者待配准遥感图像各向异性尺度空间的第n层图像Iⁿ沿竖直向下y轴正方向的差分，n＝0,1,...,N-1，N表示参考遥感图像和待配准遥感图像各向异性尺度空间层的总数。

步骤4，计算梯度幅度图像的差分。

按照下式，使用索贝尔算子Sobel，计算参考遥感图像和待配准遥感图像各向异性尺度空间的梯度幅度图像的x轴正方向的差分：

其中，表示参考遥感图像或者待配准遥感图像各向异性尺度空间的梯度幅度图像GIⁿ沿水平向右的x轴正方向的差分，GIⁿ表示参考遥感图像或者待配准遥感图像各向异性尺度空间的梯度幅度图像的第n层图像，表示相关操作，n＝0,1,...,N-1，N表示参考遥感图像和待配准遥感图像各向异性尺度空间层的总数。

按照下式，使用索贝尔算子Sobel，计算参考遥感图像和待配准遥感图像各向异性尺度空间的梯度幅度图像的y轴正方向的差分：

其中，表示参考遥感图像或者待配准遥感图像各向异性尺度空间的梯度幅度图像GIⁿ沿竖直向下的y轴正方向的差分，GIⁿ表示输入参考遥感图像或者待配准遥感图像各向异性尺度空间的梯度幅度图像的第n层图像，表示相关操作，n＝0,1,...,N-1，N表示参考遥感图像和待配准遥感图像各向异性尺度空间层的总数。

步骤5，计算梯度幅度图像的梯度幅度。

按照下式，计算梯度幅度图像的梯度幅度：

其中，GGIⁿ表示参考遥感图像或者待配准遥感图像各向异性尺度空间的梯度幅度图像的第n层图像GIⁿ的梯度幅度，表示参考遥感图像或者待配准遥感图像各向异性尺度空间的梯度幅度图像的第n层图像GIⁿ沿水平向右x轴正方向的差分，表示参考遥感图像或者待配准遥感图像各向异性尺度空间的梯度幅度图像的第n层图像GIⁿ沿竖直向下y轴正方向的差分，表示开根号操作，n＝0,1,...,N-1，N表示参考遥感图像和待配准遥感图像各向异性尺度空间层的总数。

步骤6，计算梯度幅度图像的梯度角度。

按照下式，计算梯度幅度图像的梯度角度：

其中，AGGIⁿ表示参考遥感图像或者待配准遥感图像各向异性尺度空间的梯度幅度图像的第n层图像GIⁿ的梯度角度，表示参考遥感图像或者待配准遥感图像各向异性尺度空间的梯度幅度图像的第n层图像GIⁿ沿水平向右x轴正方向的差分，表示参考遥感图像或者待配准遥感图像各向异性尺度空间的梯度幅度图像的第n层图像GIⁿ沿竖直向下y轴正方向的差分，arctan(·)表示四象限反正切操作，n＝0,1,...,N-1，N表示参考遥感图像和待配准遥感图像各向异性尺度空间层的总数。

步骤7，检测特征点。

使用Harris角点检测算子，在参考遥感图像或者待配准遥感图像的各向异性尺度空间的梯度幅度图像上检测特征点。

按照下式，计算参考遥感图像或者待配准遥感图像各向异性尺度空间的梯度幅度图像各层图像的每个像素的Harris矩阵：

其中，u(X,σ_n)表示参考遥感图像或者待配准遥感图像各向异性尺度空间的梯度幅度图像的第n层图像GIⁿ位置为X的像素的Harris矩阵，X表示像素位置坐标，σ_n表示参考遥感图像或者待配准遥感图像各向异性尺度空间的第n层的尺度值，是x轴方向和y轴方向标准差都是的高斯函数，*表示卷积操作。

计算参考遥感图像或者待配准遥感图像各向异性尺度空间的梯度幅度图像的各层图像的每个像素的Harris函数：

R(X,σ_n)＝det(u(X,σ_n))-D·tr(u(X,σ_n))²

其中，R(X,σ_n)表示参考遥感图像或者待配准遥感图像各向异性尺度空间的梯度幅度图像的第n层图像GIⁿ位置为X的像素的Harris函数，det(·)表示方阵的行列式，D表示任意常数，tr(·)表示矩阵主对角线元素的总和。

对于R(X,σ_n)＞T_th且R(X,σ_n)大于其8邻域内任意一点的Harris函数的点作为特征点，得到参考遥感图像和待配准遥感图像的各向异性尺度空间的梯度幅度图像上检测的特征点的集合：

其中，T_th表示Harris函数阈值。

步骤8，生成特征点主方向。

将水平方向角度在[0,2π]内36等分。

由特征点集合的半径为6σ_n的圆形邻域区域内像素的梯度方向AGGIⁿ(X)，确定圆形区域内像素所在的等分角度范围，其中，σ_n表示参考遥感图像或者待配准遥感图像各向异性尺度空间的第n层的尺度值，AGGIⁿ(X)表示特征点周围圆形邻域内位置坐标为X的像素的梯度角度，n＝0,1,...,N-1，N表示参考遥感图像和待配准遥感图像各向异性尺度空间层的总数。

分别累加位于每个等分角度范围内的所有像素的梯度幅度值GGIⁿ(X)，形成梯度方向直方图，其中，GGIⁿ(X)表示特征点周围圆形邻域内位置坐标为X的像素的梯度幅度，n＝0,1,...,N-1，N表示参考遥感图像和待配准遥感图像各向异性尺度空间层的总数。

将梯度方向直方图中大于最大值0.8倍的数值对应的梯度方向，作为特征点的主方向。

步骤9，生成特征点描述子。

(9a)将半径为ρ的特征点的圆形邻域逆时针旋转θ_m度，得到附图2(a)所示的相同半径的圆形邻域，圆形邻域沿半径方向将ρ分为3个区间，内圆邻域半径为0.25ρ，中间圆邻域半径为0.73ρ，外圆邻域半径为ρ，圆形邻域沿角度方向将[0,2π]等分为8个区间，内圆作为一个整体，特征点周围圆形邻域被划分为了17个不同位置的面积相等的子区域，其中，ρ的取值为12σ_n，σ_n表示特征点所在的参考遥感图像或者待配准遥感图像各向异性尺度空间的第n层的尺度值，n＝0,1,...,N-1，N表示参考遥感图像和待配准遥感图像各向异性尺度空间层的总数，θ_m表示该特征点的主方向。

(9b)将特征点圆形邻域内像素的笛卡尔坐标转换为对数极坐标，对数极坐标角度方向水平向右，附图2(a)中标注的阿拉伯数字区域对应附图2(b)中标注的阿拉伯数字区域，在[0,2π]的范围内等分为8个区间，对数极坐标对数长度方向竖直向下，在范围内非等分为3个区间，其中，ρ表示特征点周围圆形邻域的半径。

(9c)对数极坐标网格中每个子区域内的所有像素根据其梯度幅度GGIⁿ(X)和梯度方向AGGIⁿ(X)计算梯度方向直方图，每一个子区域形成了一个8维的梯度方向向量，依次拼接17子区域的梯度方向向量就形成了一个136维的特征点的描述子，其中，梯度方向直方图的角度在[0,2π]的范围内等分为8个区间，GGIⁿ(X)表示特征点周围圆形邻域内位置坐标为X的像素的梯度幅度，AGGIⁿ(X)表示特征点周围圆形邻域内位置坐标为X的像素的梯度角度。

(9d)采用步骤(9a)、步骤(9b)、步骤(9c)的相同方法，生成特征点集合C_IR的描述子集合D_IR和特征点集合C_IS的描述子集合D_IS，其中，D_IR表示参考遥感图像的各向异性尺度空间的梯度幅度图像上检测到的特征点的描述子集合，D_IS表示待配准遥感图像的各向异性尺度空间的梯度幅度图像上检测到的特征点的描述子集合，L₁表示参考遥感图像的各向异性尺度空间的梯度幅度图像上检测到的特征点的描述子总数，L₂表示待配准遥感图像的各向异性尺度空间的梯度幅度图像上检测到的特征点的描述子总数。

步骤10，特征点匹配。

匹配在参考遥感图像的各向异性尺度空间梯度幅度图像上检测到的特征点和在待配准遥感图像的各向异性尺度空间梯度幅度图像上检测到的特征点：

(10a)将待配准遥感图像特征点的描述子序号k初始化为1。

(10b)按照下式，计算待配准遥感图像特征点的描述子与参考遥感图像每一个特征点的描述子的欧式距离组成的向量：

其中，d_k表示待配准遥感图像特征点的描述子与参考遥感图像的每一个描述子的欧式距离组成的向量，d_k,j表示配准遥感图像特征点的描述子与参考遥感图像特征点的描述子的欧式距离，表示待配准遥感图像特征点的第k个描述子，表示参考遥感图像特征点的第j个描述子，j＝1,2,...,L₁，L₁表示参考遥感图像特征点描述子的个数，||·||表示取范数操作。

(10c)对向量d_k中的L₁个元素从小到大排序得到向量s_k。

(10d)如果则描述子对应的特征点和与描述子的欧式距离为s_k,0的参考遥感图像的描述子对应的特征点匹配；如果且则描述子对应的特征点和与描述子的欧式距离为d_k的前E个元素的参考遥感图像的描述子对应的特征点匹配；如果则描述子对应的特征点不和参考遥感图像的任何特征点匹配；

其中，s_k,0表示向量s_k的第一个元素，s_k,1表示向量s_k的第二个元素，T_L和T_H表示最近邻和次近邻距离比阈值，E的取值范围为[2,3]中的整数。

(10e)判断k≥L₂是否成立，若成立，得到待配准遥感图像和参考遥感图像的特征点匹配关系C_IR,IS，执行步骤(11)；否则，令k＝k+1，执行本步骤的步骤(10b)；

其中，k表示待配准遥感图像特征点的描述子序号，C_IR,IS表示参考遥感图像的各向异性尺度空间的梯度幅度图像上检测到的特征点和在待配准遥感图像的各向异性尺度空间的梯度幅度图像上检测到的特征点的配对关系集合，L₃表示特征点匹配阶段匹配的点对总数。

步骤11，删除错误匹配的特征点对。

参照附图3，使用随机抽样一致RANSAC算法，删除错误的匹配点对，得到待配准遥感图像到参考遥感图像的几何变换关系，具体步骤如下：

(11a)初始化循环次数Num为0，初始化最优一致点集C中包含的匹配点对个数为0。

(11b)Num＝Num+1，判断Num>1000是否成立，若成立，执行步骤(11i)；否则，执行步骤(11c)。

(11c)从特征点匹配关系C_IR,IS中随机选取3个不同的匹配特征点对；

其中，C_IR,IS表示参考遥感图像的各向异性尺度空间的梯度幅度图像上检测到的特征点和在待配准遥感图像的各向异性尺度空间的梯度幅度图像上检测到的特征点的配对关系集合。

(11d)使用仿射变换模型，计算满足3个不同的匹配特征点对的变换矩阵T1。

(11e)计算特征点匹配关系C_IR,IS中满足变换矩阵T1的一致点集Con；

其中，C_IR,IS表示参考遥感图像的各向异性尺度空间的梯度幅度图像上检测到的特征点和在待配准遥感图像的各向异性尺度空间的梯度幅度图像上检测到的特征点的配对关系集合。

(11f)判断一致点集Con中包含的点对个数大于一致点对阈值TH是否成立，若成立，执行步骤(11g)；否则，执行步骤(11b)；

其中，TH表示一致点对阈值，TH取值为L₃的0.1倍，L₃表示特征点匹配阶段匹配的点对总数。

(11g)判断一致点集Con中包含的点对个数大于当前最优一致点对个数是否成立，若成立，执行步骤(11h)；否则，执行步骤(11b)。

(11h)一致点集Con代替当前最优一致点集C，使用仿射变换模型，计算满足一致点集C的变换矩阵T2，执行步骤(11b)。

(11i)输出当前最优一致点集C和变换矩阵T2。

输出的最优一致点集C中包含参考遥感图像和待配准遥感图像正确匹配的特征点对；变换矩阵T2即为获得的待配准遥感图像到参考遥感图像的几何变换关系T。

步骤12，配准参考图像和待配准图像。

根据待配准遥感图像到参考遥感图像的几何变换关系T，使用双线性插值方法配准参考遥感图像和待配准遥感图像。

下面结合仿真图对本发明的效果做进一步的说明。

1.仿真条件：

本发明仿真实验在Intel(R)Core(TM)i3 CPU M 380 2.53GHz Windows 7系统下，Matlab 2010a运行平台上，完成本发明的仿真实验。

2.仿真实验内容：

本发明的仿真实验中取参考遥感和待配准遥感图像尺度空间的组数G为3，组内层数L为3，因此，各向异性尺度空间的层的总数N为9；基准尺度值σ₀为1.6，仿真实验中输入参考遥感图像和待配准遥感图像数据转换成了0-1之间的浮点数据类型，Harris函数阈值T_th为0.005，特征点匹配中阈值T_L和T_H分别设置为0.6和0.97，图像间几何变换模型是仿射变换。本发明的仿真实验输入的测试遥感图像分为两类：(1)多光谱遥感图像对，(2)多源遥感图像对。

附图4(a)和附图4(b)是第一组多光谱遥感图像对P-A，其中，附图4(a)是参考图像，大小为761×748像素，传感器类型是Landsat-7ETM+，5波段，获取时间是2000/7/24；附图4(b)是待配准图像，大小为761×748像素，传感器类型是Landsat 4–5TM，3波段，获取时间是1999/7/6。

附图5(a)和附图5(b)是第二组多光谱遥感图像对P-B，其中，附图5(a)是参考图像，大小为761×748像素，传感器类型是Landsat-7ETM+，3波段，获取时间是2003/4/12；附图5(b)是待配准图像，大小为761×748像素，传感器类型是Landsat 4–5TM，5波段，获取时间是2006/06/15。

附图6(a)和附图6(b)是多源遥感图像对P-C，其中，附图6(a)是参考图像，大小为800×800像素，传感器类型是ALOS-PALSAR，获取时间是2010/6/5；附图6(b)是待配准图像，大小为800×800像素，传感器类型是Landsat ETM+，获取时间是1999/6/26。

本发明的仿真实验中首先对比经典的SIFT算法，还对比了文中提出的双阈值最近邻匹配准则和经典的单阈值最近邻匹配准则，因此，对于每个测试图像对应用两种不同的方法和文中提出的方法进行对比，这两种方法是SIFT算法和Harris各向异性梯度尺度空间+最近邻距离比匹配+RANSAC算法(Harris-NDGSS-NNDR-RANSAC)，分别简称为SIFT和HNNR，文中提到的完整算法是Harris各向异性梯度尺度空间+双阈值最近邻距离比匹配+RANSAC算法(Harris-NDGSS-DNNDR-RANSAC)，简称为HNDR。本发明仿真实验中使3种不同的方法对同一个图像对尽可能产生同样数量的特征点，仿真实验结果下表所示。

3.匹配性能评估方法：

使用正确匹配的特征点对个数N_c和均方根误差RMSE作为评估配准性能的标准。N_c是使用随机抽样一致RANSAC算法后得到的正确匹配的特征点对个数。

均方根误差RMSE是根据下式计算得到的：

其中，x_r和y_r分别表示参考遥感图像上手动选取的第r个点的列坐标和行坐标，和分别表示待配准遥感图像上手动选取的第r个点经过待配准遥感图像到参考遥感图像的几何变换关系T变换后的列坐标和行坐标，r＝1,2,...,N_GT，N_GT表示从参考遥感图像和待配准遥感图像上手动选取的真实匹配点对总数。

4.仿真实验结果及分析：

从附录图4(a)和图4(b)可以看出，参考遥感图像和待配准遥感图像的大部分相同区域亮度反转；如下表所示，SIFT算法不能实现正确配准，但是基于各向异性梯度尺度空间的HNNR和HNDR算法能检测较多的正确匹配点对，验证了文章提出的算法能配准图像亮度反转的遥感图像对。附图5(a)和附图5(b)也是多光谱遥感图像对，两幅图个别区域图像亮度反转或者非线性变化；如下表所示，基于SIFT的方法也能精确实现配准，但是基于各向异性梯度尺度空间的HNNR和HNDR算法能检测更多的正确匹配点对，也验证了文章提出的算法能配准图像亮度非线性变化或者亮度反转的遥感图像对。附图6(a)和附图6(b)是多源遥感图像对，相同区域图像亮度主要呈现线性变化模式，如下表所示，3种方法都能实现正确配准，且SIFT算法也能实现较多的正确匹配对，这是因为SIFT算法主要检测图像纹理结构，在纹理边缘保存较好的图像上能提取更多有价值的点特征，也证明了在一定范围的遥感图像配准任务中可以使用SIFT算法。其中，附图4(c)、附图5(c)和附图6(c)是待配准遥感图像向参考遥感图像坐标系变换后的配准结果，附图4(d)、附图5(d)和附图6(d)是参考遥感图像和变换后的待配准遥感图像的融合。

采用HNNR算法和HNDR算法，对遥感图像对进行的配准结果如下表所示。由下表可以看出，HNDR算法正确匹配的特征点个数都多于HNNR算法正确匹配的特征点个数，HNDR算法和HNNR算法的不同之处在于HNDR算法采用的是文中提出的双阈值最近邻距离比匹配准则进行特征点匹配的，而HNNR算法采用的是经典的最近邻距离比匹配准则，由此验证了双阈值最近邻距离比匹配方法能更加有效的利用特征点描述子信息进行特征点匹配。通过比较均方根误差RMSE，可以看出提出的方法能够实现子像素精度的配准。总之，HNDR算法不仅能较好的完成多源遥感图像的配准任务，同时克服了现有算法配准图像亮度存在较大非线性变化甚至图像亮度反转的多光谱遥感图像对正确率低的问题。

其中，Pair表示测试的遥感图像对，Method表示匹配遥感图像对的不同方法，N_IR和N_IS分别表示参考遥感图像和待配准遥感图像检测到的特征点个数，N_c表示正确匹配的特征点对个数，RMSE表示均方根误差。