CN105631872B - 基于多特征点的遥感图像配准方法 - Google Patents

Abstract

本发明公开了一种基于多特征点的遥感图像配准方法，主要解决传统多源和多光谱遥感图像配准精度低的缺点。其实现步骤为：1.输入两幅遥感图像；2.构造输入图像的各向异性尺度空间；3.在输入图像的各向异性尺度空间中分别使用Harris和Hessian算子进行特征点检测；4.分别将输入图像的两种特征点检测结果结合并生成特征向量；5.将输入图像的特征向量进行匹配，删除错误匹配的特征点对，输出配准结果。本发明具有配准精度高的优点，可用于多源遥感图像和多光谱遥感图像的配准。

Description

基于多特征点的遥感图像配准方法

技术领域

本发明属于图像处理技术领域，涉及遥感图像配准，可应用于多源、多光谱遥感图像的配准。

背景技术

图像配准是图像处理中非常关键的一个步骤，它是指将不同时刻、不同视角或者是不同传感器获取的同一场景的两幅或者多幅图像进行叠加对准的过程。图像配准技术已经被广泛地用于图像变化检测、图像拼接、医学领域以及模式识别领域。图像配准的方法大致可以分为两类：基于灰度信息的图像配准方法和基于特征信息的图像配准方法。

基于灰度的图像配准是以整幅图像像素间的灰度信息为依据，建立待配准图像和参考图像之间的相似性度量函数，利用某一种搜索算法，寻找出使得相似性度量函数达到最优值的变换模型参数。这种算法容易实现，在配准之前不用提取图像的特征，只需要获取其灰度信息，但是其应用范围较窄，不能直接用于校正图像的非线性形变，且在最优变换的搜索过程中需要巨大的计算量。

基于特征的图像配准方法是目前图像配准最常用的方法之一，其最大的优点在于能够将对整个图像进行的各种分析转化为对图像特征信息，即特征点、特征曲线、边缘、较小的区域的分析，从而大大减小了图像处理过程的运算量，而且对灰度变化、图像变形以及遮挡等都有较好的适应能力，以及能够实现在复杂成像条件下图像的快速、精确配准。经典的特征点提取算子有：Harris算子，Hessian算子等，但是Harris算子主要针对角点特征检测，对于纹理特征较多的图像检测效果不佳导致图像配准不精确，Hessian算子主要针对纹理特征检测，对于角点信息较多的图像配准结果不准确。

发明内容

本发明的目的在于针对已有技术的不足，提出一种基于多特征点的遥感图像配准方法，以提高配准效果，满足多传感器和多光谱图像的配准要求。

为实现上述目的，本发明的技术方案包括如下：

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

(2)构造参考遥感图像I₁和待配准遥感图像I₂的各向异性尺度空间图像：

(2a)计算各向异性尺度空间各层的尺度值；

其中，σ_i表示各向异性尺度空间的第i层图像尺度值，σ₀表示尺度参数的初始基准值，i＝0,1,2,...,L-1，i表示各向异性尺度空间层的序号，L表示各向异性尺度空间层的总数；

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

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

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

(2e)计算参考遥感图像I₁和待配准遥感图像I₂各向异性尺度空间的第i层图像的扩散系数矩阵；

(2f)计算参考遥感图像I₁和待配准遥感图像I₂各向异性尺度空间的第i+1层图像：

Iⁱ⁺¹＝(I-(t_i+1-t_i)·A(Iⁱ))^-1Iⁱ，

其中，Iⁱ⁺¹表示参考遥感图像I₁或者待配准遥感图像I₂各向异性尺度空间的第i+1层图像，A(Iⁱ)表示各向异性尺度空间的第i层图像的扩散矩阵，I表示一个与A(Iⁱ)同等大小的单位矩阵，t_i和t_i+1分别表示各向异性尺度空间的第i层和第i+1层的时间度量值，Iⁱ表示各向异性尺度空间的第i层图像，(·)^-1表示逆矩阵操作；

(2g)判断i≥L-1是否成立，若成立，得到参考遥感图像I₁和待配准遥感图像I₂的各向异性尺度空间图像，否则，令i＝i+1，返回步骤(2e)；

(3)分别在参考遥感图像I₁和待配准遥感图像I₂的各向异性尺度空间图像上使用Hessian算子进行纹理特征检测，得到参考遥感图像I₁的第一特征点集P₁和待配准遥感图像I₂的第一特征点集Q₁，特征点集中保存的是特征点的坐标信息；

(4)分别在参考遥感图像I₁和待配准遥感图像I₂的各向异性尺度空间图像上使用Harris算子进行角点特征检测，得到参考遥感图像I₁的第二特征点集P₂和待配准遥感图像I₂的第二特征点集Q₂；

(5)将参考遥感图像I₁的第一特征点集P₁与第二特征点集P₂进行结合，删除重复的坐标点后得到参考遥感图像I₁的最终特征点集P；将待配准遥感图像的第一特征点集Q₁与第二特征点集Q₂进行结合，删除重复的坐标点后得到待配准遥感图像的最终特征点集Q；

(6)分别生成参考遥感图像I₁的特征点集P的特征向量D₁和待配准遥感图像I₂的特征点集Q的特征向量D₂；

(7)将参考遥感图像I₁的特征向量D₁和待配准遥感图像I₂的特征向量D₂进行匹配，得到初始匹配点对；

(8)使用随机抽样一致性算法RANSAC，提纯初始匹配点对，剔出误匹配点，得到待配准遥感图像I₂到参考遥感图像I₁的仿射变换参数；

(9)根据仿射变换参数值，对待配准遥感图像I₂进行仿射变换，得到仿射变换后的图像F₁；

(10)将仿射变换后的图像F₁与参考遥感图像I₁进行融合，得到融合后的图像。

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

第一，本发明是通过各向异性扩散方程构建尺度空间，并在各向异性尺度空间上检测特征点，各向异性扩散方程具有选择性的扩散平滑、较好的兼顾噪声去除和边缘保护。

第二，本发明使用两种不同的检测方法进行特征点检测，使得检测的特征不仅有纹理信息，而且还有角点信息，使得特征点信息更加全面、丰富和稳定，配准结果更加精确。

附图说明

图1为本发明的实现流程图；

图2为用本发明对多源遥感图像的配准结果图；

图3为用本发明对多光谱遥感图像的配准结果图。

具体实施方式

参考附图1，本发明的实现步骤如下：

步骤1，输入两幅图像

通过传感器拍摄的两幅图像，其中一幅作为参考遥感图像I₁，另一幅作为待配准遥感图像I₂。

步骤2，构建参考遥感图像I₁和待配准遥感图像I₂的各向异性尺度空间图像

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

其中，σ_i表示参考遥感图像I₁或者待配准遥感图像I₂的各向异性尺度空间的第i层图像的尺度值，σ₀表示尺度参数的初始基准值，设置为1.6，i表示参考遥感图像I₁或者待配准遥感图像I₂各向异性尺度空间层的序号，i＝0,1,2,...,L-1，L表示遥感图像I₁或者待配准遥感图像I₂各向异性尺度空间包含的图像总数，设置为8；

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

其中，t_i是进化时间，其表示参考遥感图像I₁或者待配准遥感图像I₂的各向异性尺度空间的第i层图像尺度的时间度量值；

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

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

(2e)计算参考遥感图像I₁和待配准遥感图像I₂各向异性尺度空间的第i层图像的扩散系数矩阵；

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

(2f)计算参考遥感图像I₁和待配准遥感图像I₂各向异性尺度空间的第i+1层图像：

Iⁱ⁺¹＝(I-(t_i+1-t_i)·A(Iⁱ))^-1Iⁱ，

其中，A(Iⁱ)表示参考遥感图像I₁或者待配准遥感图像I₂各向异性尺度空间的第i层图像Ii的扩散矩阵，I是一个与A(Iⁱ)同等大小的单位矩阵，t_i和t_i+1分别是各向异性尺度空间的第i层和第i+1层的时间度量值，(·)^-1表示逆矩阵操作，Iⁱ⁺¹表示参考遥感图像I₁或者待配准遥感图像I₂各向异性尺度空间的第i+1层图像；

(2g)判断i≥L-1是否成立，若成立，得到参考遥感图像I₁或者待配准遥感图像I₂的各向异性尺度空间，否则，令i＝i+1，返回步骤(2e)。

步骤3，分别在参考遥感图像I₁和待配准遥感图像I₂的各向异性尺度空间图像上使用Hessian算子进行纹理特征检测，得到参考遥感图像I₁的第一特征点集P₁和待配准遥感图像I₂的第一特征点集Q₁。

(3a)计算参考遥感图像I₁和待配准遥感图像I₂各向异性尺度空间图像沿x轴方向的一阶微分和二阶微分：

其中，Iⁱ表示参考遥感图像I₁或者待配准遥感图像I₂各向异性尺度空间的第i层图像，表示图像Iⁱ沿x轴正方向的一阶微分，表示图像Iⁱ沿x轴正方向的二阶微分，表示相关操作，i＝0,1,...,L-1；

(3b)计算参考遥感图像I₁和待配准遥感图像I₂各向异性尺度空间图像沿y轴方向的一阶微分和二阶微分：

其中，表示图像Iⁱ沿y轴正方向的一阶微分，表示图像Iⁱ沿y轴正方向的二阶微分；

(3c)计算参考遥感图像I₁和待配准遥感图像I₂各向异性尺度空间图像的混合微分：

其中，表示图像Iⁱ的混合微分；

(3d)在一阶微分，二阶微分和混合微分的基础上，分别计算参考遥感图像I₁和待配准遥感图像I₂的Hessian矩阵：

其中，σⁱ是σ_i的整数值，表示各向异性尺度空间图像每个像素的Hessian矩阵；

(3e)通过计算参考遥感图像I₁和待配准遥感图像I₂的Hessian矩阵局部极大值进行特征点检测：

将每一个像素的Hessian矩阵行列式值与它所有的相邻点的Hessian矩阵行列值进行比较，理论上，比较的范围是当前尺度、上一尺度和下一尺度上的3个大小为σ_i*σ_i的矩形窗口。但是，为了加快搜索速度，窗口大小固定为3*3，则搜索空间是一个边长为3个像素的立方体：中间的检测点和它同尺度的8个相邻点，以及和上、下相邻尺度对应的9*2个像素点—共26个点比较，当每个像素的Hessian矩阵行列式值大于它所有的相邻点的Hessian矩阵行列值时，即为极值点，该极值点为检测出的纹理特征。

步骤4，分别在参考遥感图像I₁和待配准遥感图像I₂的各向异性尺度空间图像上使用Harris算子进行角点特征检测，得到参考遥感图像I₁的第二特征点集P₂和待配准遥感图像I₂的第二特征点集Q₂：

因为现有的许多角点检测技术，都是沿着水平和垂直两个方向进行，或者水平、垂直和对角线4个方向进行，容易将边缘特征点误检测为角点特征，由于Harris角点检测法是沿着泰勒公式展开，它可以沿着各个方向进行，可以有效地避免将边缘特征误检测为角点特征。其主要步骤如下：

(4a)计算参考遥感图像I₁和待配准遥感图像I₂的各向异性尺度空间各层图像的Harris矩阵：

其中，g是标准差为的高斯函数，表示各向异性尺度空间的第i层图像Iⁱ沿x轴正方向的一阶微分，表示图像Iⁱ沿y轴正方向的一阶微分，u(X,σ_i)表示图像Iⁱ中位置为X的像素的Harris矩阵，X表示像素位置坐标，*表示卷积操作，在这个过程中采用高斯函数，是因为距离中心点的位置权重越大，越远离中心点的位置，权重越小，这样可以降低噪声的影响，提高角点检测的精度；

(4b)计算各向异性尺度空间各层图像每个像素点的Harris角点响应：

R(X,σ_i)＝det(u(X,σ_i))-d·tr(u(X,σ_i))²

其中，det(·)表示矩阵的行列式，d表示任意常数，取值为0.04-0.06之间的任意数，tr(·)是矩阵的迹，即矩阵主对角线元素的总和，R(X,σ_i)表示各向异性尺度空间的第i层图像Iⁱ中位置为X的像素的Harris角点响应；

(4c)根据每个像素的Harris角点响应R判断像素点是否是角点：若Harris角点响应R大于阈值d_sh＝0.8，则该像素点为Harris角点，并将其提取出来；否则，视为非角点。

步骤5，分别对参考遥感图像I₁和待配准遥感图像I₂的特征点集进行结合：

(5a)将参考遥感图像I₁的第一特征点集P₁与第二特征点集P₂进行结合，删除重复的坐标点后得到参考遥感图像I₁的最终特征点集P；

(5b)将待配准遥感图像I₂的第一特征点集Q₁与第二特征点集Q₂进行结合，删除重复的坐标点后得到待配准遥感图像的最终特征点集Q。

步骤6，分别生成参考遥感图像I₁的特征点集P的特征向量D₁和待配准遥感图像I₂的特征点集Q的特征向量D₂：

(6a)在以特征点为中心，半径为6σ_i的圆形区域内，计算特征点邻域像素的梯度幅度和梯度方向，得到梯度方向直方图；将梯度方向直方图中大于梯度最大值0.8倍的数值对应的梯度方向，作为特征点的主方向；确定主方向的模值和方向，再将坐标轴旋转到特征点的主方向上，确保特征向量的旋转不变性；

(6b)在以特征点为中心，半径为12σ_i的圆形区域内，将笛卡尔坐标转换为对数极坐标，并将圆形区域沿半径方向分为3个区间，分别为内圆，中间圆和外圆；将内圆作为一个整体，中间圆和外圆沿半径方向等分为8个区间，共形成17个子区域；在每个子区域中计算所有像素在8个方向上的梯度幅度和梯度方向，每个子区域得到一个8维梯度直方图；最后将17个子区域的梯度方向向量组合形成一个136维的特征点的描述子，该136维特征点描述子为特征点的特征向量。

步骤7，将参考遥感图像I₁的特征向量D₁和待配准遥感图像I₂的特征向量D₂进行匹配，得到初始匹配点对：

(7a)计算参考遥感图像I₁中每个特征点的特征向量与待配准遥感图像I₂中每个特征点的特征向量的欧氏距离，并将欧氏距离从小到大进行排序，取最近邻欧氏距离和次近邻欧式距离；

(7b)将最近邻欧氏距离和次近邻欧式距离的比值设置为r，当r小于阈值d_ratio＝0.8时则被选为最优匹配点对，获得初始匹配点集。

步骤8，使用随机抽样一致性算法RANSAC，提纯初始匹配点对，剔出误匹配点，得到待配准遥感图像I₂到参考遥感图像I₁的仿射变换参数：

(8a)初始化循环次数N为0，初始化最优点集Con中包含的匹配点对个数为0；

(8b)令N＝N+1，判断N>500是否成立，若成立，则执行步骤(8h)；否则，执行步骤(8c)；

(8c)从初始匹配点集中随机抽取3个不同的匹配点对，并根据这3个匹配点对，计算变换矩阵M；

(8d)根据初始匹配点集计算满足变换矩阵M的点集C，返回点集C中的元素个数；

(8e)判断点集C中元素个数是否大于阈值T_sh，若成立，执行步骤(8f)；否则，执行步骤(8b)，T_sh设置为初始匹配点集中点对总数的0.05倍；

(8f)判断点集C中元素个数是否大于当前最优点集Con中元素个数，若成立，执行步骤(8g)；否则，执行步骤(8b)；

(8g)用点集C代替当前最优点集Con，使用仿射变换模型，计算满足点集Con的仿射变换参数，返回步骤(8b)；

(8h)输出当前最优点集Con和仿射变换参数。

步骤9，根据仿射变换参数值，对待配准遥感图像I₂进行仿射变换，得到仿射变换后的图像F₁。

步骤10，将仿射变换后的图像F₁与参考遥感图像I₁进行融合，得到融合后的图像。

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

1.仿真条件：

硬件平台为：Intel(R)Core(TM)i5CPU 2.20GHz；

软件平台为：Windows 8.1，Matlab 2010a

2.仿真实验内容：

仿真实验参数设置：取参考遥感图像I₁和待配准遥感图像I₂各向异性尺度空间总的层数L为8；基准尺度值σ₀为1.6，Harris函数阈值d_sh为0.8。

本发明的仿真实验输入的测试遥感图像分为两类：一是多光谱遥感图像对，二是多源遥感图像对。

仿真实验一：本发明仿真实验一是对两幅多源图像进行图像配准，这两幅配准实验图像分别是传感器A拍摄的图像和传感器B拍摄的图像，仿真结果如图2所示，其中：图2(a)是传感器A拍摄的图像，即参考遥感图像I₁，图2(b)是传感器B拍摄的图像，即待配准遥感图像I₂，图2(c)是使用本发明的方法对参考遥感图像和待配准遥感图像的配准结果，图2(d)是使用本发明的方法对参考图像和配准后的图像的融合结果。

仿真实验二：本发明仿真实验二是两幅多光谱图像进行图像配准，这两幅配准实验图像是大小均为512*512的TM影像，其中图3(a)是12波段的参考遥感图像I₁，图3(b)是9波段的待配准遥感图像I₂，图3(c)是使用本发明对参考遥感图像和待配准遥感图像的配准结果，图3(d)是使用本发明对参考图像和配准后的图像的融合结果。

由图2(c)，图3(c)的配准图可以看出，本发明找到的多源待配准遥感图像可以和多源参考遥感图像完全对齐，多光谱待配准遥感图像也可以和输入的多光谱参考遥感图像完全对齐，结果没有移位，完全配准。

由图2(d)，图3(d)的融合图可以看出，本发明找到的多源遥感图像对和多光谱遥感图像对的重叠区域可以很精确地融合，结果几乎没有错位。

Claims (8)

1.一种基于多特征点的遥感图像配准方法，包括如下步骤：

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

(2)构造参考遥感图像I₁和待配准遥感图像I₂的各向异性尺度空间图像：

(2a)计算各向异性尺度空间各层的尺度值；

其中，σ_i表示各向异性尺度空间的第i层图像尺度值，σ₀表示尺度参数的初始基准值，i＝0,1,2,...,L-1，i表示各向异性尺度空间层的序号，L表示各向异性尺度空间层的总数；

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

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

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

(2e)计算参考遥感图像I₁和待配准遥感图像I₂各向异性尺度空间的第i层图像的扩散系数矩阵；

(2f)计算参考遥感图像I₁和待配准遥感图像I₂各向异性尺度空间的第i+1层图像：

其中，表示参考遥感图像I₁各向异性尺度空间的第i+1层图像，表示待配准遥感图像I₂各向异性尺度空间的第i+1层图像，表示参考遥感图像I₁各向异性尺度空间的第i层图像的扩散矩阵，E₁表示一个与同等大小的单位矩阵，表示待配准遥感图像I₂各向异性尺度空间的第i层图像的扩散矩阵，E₂表示一个与同等大小的单位矩阵，t_i和t_i+1分别表示各向异性尺度空间的第i层和第i+1层的时间度量值，表示参考遥感图像I₁各向异性尺度空间的第i层图像，表示待配准遥感图像I₂各向异性尺度空间的第i层图像，(·)^-1表示逆矩阵操作；

(2g)判断i≥L-1是否成立，若成立，得到参考遥感图像I₁和待配准遥感图像I₂的各向异性尺度空间图像，否则，令i＝i+1，返回步骤(2e)；

(3)分别在参考遥感图像I₁和待配准遥感图像I₂的各向异性尺度空间图像上使用Hessian算子进行纹理特征检测，得到参考遥感图像I₁的第一特征点集P₁和待配准遥感图像I₂的第一特征点集Q₁，特征点集中保存的是特征点的坐标信息；

(4)分别在参考遥感图像I₁和待配准遥感图像I₂的各向异性尺度空间图像上使用Harris算子进行角点特征检测，得到参考遥感图像I₁的第二特征点集P₂和待配准遥感图像I₂的第二特征点集Q₂；

(5)将参考遥感图像I₁的第一特征点集P₁与第二特征点集P₂进行结合，删除重复的坐标点后得到参考遥感图像I₁的最终特征点集P；将待配准遥感图像的第一特征点集Q₁与第二特征点集Q₂进行结合，删除重复的坐标点后得到待配准遥感图像的最终特征点集Q；

(6)分别生成参考遥感图像I₁的特征点集P的特征向量D₁和待配准遥感图像I₂的特征点集Q的特征向量D₂；

(7)将参考遥感图像I₁的特征向量D₁和待配准遥感图像I₂的特征向量D₂进行匹配，得到初始匹配点对；

(8)使用随机抽样一致性算法RANSAC，提纯初始匹配点对，剔出误匹配点，得到待配准遥感图像I₂到参考遥感图像I₁的仿射变换参数；

(9)根据仿射变换参数值，对待配准遥感图像I₂进行仿射变换，得到仿射变换后的图像F₁；

(10)将仿射变换后的图像F₁与参考遥感图像I₁进行融合，得到融合后的图像。

2.根据权利要求1所述的方法，其中步骤(2b)中将尺度空间值转换到时间度量值，按照下式进行：

其中，t_i表示各向异性尺度空间的第i层图像尺度的时间度量值。

3.根据权利要求1所述的方法，其中步骤(2e)中计算参考遥感图像I₁和待配准遥感图像I₂各向异性尺度空间的第i层图像的扩散系数矩阵，按如下公式计算：

其中，gⁱ表示各向异性尺度空间中第i层图像的扩散系数矩阵，表示各向异性尺度空间的第i层图像Iⁱ使用标准差为1的高斯滤波后的图像，表示高斯滤波后的图像的梯度幅度，|·|表示取模操作，K表示对比度因子，K的取值为梯度幅度的统计直方图70％的梯度幅度值。

4.根据权利要求1所述的方法，其中步骤(3)中使用Hessian算子在各向异性尺度空间图像上进行纹理特征检测，按如下步骤进行：

(3a)计算各向异性尺度空间图像沿x轴方向的一阶微分和二阶微分：

其中，Iⁱ表示各向异性尺度空间的第i层图像，表示图像Iⁱ沿x轴正方向的一阶微分，表示图像Iⁱ沿x轴正方向的二阶微分，表示相关操作；

(3b)计算各向异性尺度空间图像沿y轴方向的一阶微分和二阶微分：

其中，表示图像Iⁱ沿y轴正方向的一阶微分，表示图像Iⁱ沿y轴正方向的二阶微分；

(3c)计算各向异性尺度空间图像的混合微分：

其中，表示图像Iⁱ的混合微分；

(3d)在一阶微分，二阶微分和混合微分的基础上，计算Hessian矩阵：

其中，σⁱ是σ_i的整数值，表示各向异性尺度空间图像每个像素的Hessian矩阵；

(3e)通过Hessian矩阵局部极大值进行特征点检测：

将每一个像素的Hessian矩阵行列式值与它同尺度的8个相邻点，以及上、下相邻尺度对应的9*2个像素点—共26个点的Hessian矩阵行列值进行比较，当其值大于它的相邻点的Hessian矩阵行列式时，即为极值点,该极值点为检测出的纹理特征点。

5.根据权利要求1所述的方法，其中步骤(4)中使用Harris算子在各向异性尺度空间图像上进行角点特征检测，按如下步骤进行：

(4a)计算各向异性尺度空间各层图像的Harris矩阵：

其中，g是标准差为的高斯函数，表示各向异性尺度空间的第i层图像Iⁱ沿x轴正方向的一阶微分，表示图像Iⁱ沿y轴正方向的一阶微分，u(X,σ_i)表示图像Iⁱ中位置为X的像素的Harris矩阵，X表示像素位置坐标，*表示卷积操作；

(4b)计算各向异性尺度空间各层图像每个像素点的Harris角点响应：

R(X,σ_i)＝det(u(X,σ_i))-d·tr(u(X,σ_i))²

其中，det(·)表示矩阵的行列式，d表示任意常数，取值为0.04-0.06之间的任意数，tr(·)是矩阵的迹，R(X,σ_i)表示各向异性尺度空间的第i层图像Iⁱ中位置为X的像素的Harris角点响应；

(4c)根据每个像素的Harris角点响应判断像素点是否是角点：若Harris角点响应R大于某一阈值d_sh，则该像素点为Harris角点，并将其提取出来；否则，视为非角点，设置阈值d_sh为0.8。

6.根据权利要求1所述的方法，其中步骤(6)中分别生成参考遥感图像I₁的特征点集P的特征向量D₁和待配准遥感图像I₂的特征点集Q的特征向量D₂，按如下步骤进行：

(6a)在以特征点为中心，半径为6σ_i的圆形区域内计算特征点邻域像素的梯度幅度和梯度方向，得到梯度方向直方图；将梯度方向直方图中大于最大值0.8倍的数值对应的梯度方向，作为特征点的主方向；确定主方向的模值和方向，将坐标轴旋转到特征点的主方向上，确保特征向量的旋转不变性；

(6b)在以特征点为中心，半径为12σ_i的圆形区域内，将笛卡尔坐标转换为对数极坐标，并将圆形区域沿半径方向分为3个区间，分别为内圆，中间圆和外圆；将内圆作为一个整体，中间圆和外圆沿半径方向等分为8个区间，共形成17个子区域；在每个子区域中计算所有像素在8个方向上的梯度幅度和梯度方向，每个子区域得到一个8维梯度直方图；最后将17个子区域的梯度方向向量组合形成一个136维的特征点的描述子，该136维特征点描述子为特征点的特征向量。

7.根据权利要求1所述的方法，其中步骤(7)中将参考遥感图像I₁的特征向量D₁和待配准遥感图像I₂的特征向量D₂进行匹配，按如下步骤进行：

(7a)计算参考遥感图像I₁中每个特征点的特征向量与待配准遥感图像I₂中每个特征点的特征向量的欧氏距离，并将欧氏距离从小到大进行排序，取最近邻欧氏距离和次近邻欧式距离；

(7b)将最近邻欧氏距离和次近邻欧式距离的比值设置为r，当r小于阈值d_ratio时则被选为最优匹配点对，获得初始匹配点集，d_ratio设置为0.8。

8.根据权利要求1所述的方法，其中步骤(8)中使用随机抽样一致性算法RANSAC，提纯初始匹配点对，剔出误匹配点，得到待配准遥感图像I₂到参考遥感图像I₁的仿射变换参数，按如下步骤进行：

(8a)初始化循环次数N为0，初始化最优点集Con中包含的匹配点对个数为0；

(8b)令N＝N+1，判断N>500是否成立，若成立，执行步骤(8h)；否则，执行步骤(8c)；

(8c)从初始匹配点集中随机抽取3个不同的匹配点对，并根据这3个匹配点对，计算变换矩阵M；

(8d)根据初始匹配点集计算满足变换矩阵M的点集C，返回点集C中的元素个数；

(8e)判断点集C中元素个数是否大于阈值T_sh，若成立，执行步骤(8f)；否则，执行步骤(8b)，T_sh设置为初始匹配点集中点对总数的0.05倍；

(8f)判断点集C中元素个数是否大于当前最优点集Con中元素个数，若成立，执行步骤(8g)；否则，执行步骤(8b)；

(8g)用点集C代替当前最优点集Con，使用仿射变换模型，计算满足点集Con的仿射变换参数，返回步骤(8b)；

(8h)输出当前最优点集Con和仿射变换参数。