CN100587518C - 遥感影像高精度控制点自动选择方法 - Google Patents

Abstract

本发明公开一种大面幅遥感影像高精度控制点自动选择方法，提出的拟合重采样控制点自动选择方法是：对一定大小的子图像进行高精度拟合配准，然后从配准的子图像中重新采样作为大面幅影像的控制点。其最大优点是，最终的控制点不是采用传统方法来直接检测，而是采用间接方法来实现的，避免了传统方法直接提取高精度控制点的困难，同时可以方便地调整最终控制点的分布。本发明对大量实验结果，采用多种不同的精度评价方法进行分析和比较，结果表明，所选择的控制点可以达到子像素级精度，是一种切实可行的并具有很强鲁棒性的控制点选择方法。

Description

遥感影像高精度控制点自动选择方法

技术领域

本发明属于计算机视觉和遥感领域，具体地说，涉及到遥感影像配准中控制点的自动选择。

背景技术

在遥感影像的很多具体应用中，要求的配准精度很高[1]Dai X andKhorram S.The effects of image misregistration on the accuracy of remotelysensed change detection.IEEE Trans.Geosci.Remote Sensing，1998，vol.36，pp.1566-1577；[2]Townshend J R G，Justice C O，Gurney C and McManusJ.The impact of misregistration on the detection of changes in land cover.IEEE Trans.Geosci.Remote Sensing，1992，vol.30，pp.1054-1060，通常要达到亚像素级，如遥感信息融合、资源变化测量、图像镶嵌等。高精度配准意味着映射模型的合理选择和模型参数的精确估计，通常，模型类型根据先验知识来确定，而模型参数的估计则完全依赖于控制点的性能。

遥感影像配准的几何映射模型通常分为两种类型[6](ZitováB andFlusser J.Image registration methods：a survey.Imaging and VisionComputing，2003，21：977～1000)。第一种是全局映射模型，即整幅图像的变换模型只有一个，如常用的刚体变换、多项式变换等。由于遥感图像尺寸大，通常每个方向达到数千到上万像素，场景覆盖范围广，可达到几万到几十万平方公里，成像环境复杂，场景内容变化多样，因此，一般情况下整幅影像采用一个全局映射模型是不恰当的。第二种映射方式是局部映射模型，这种模型可以处理局部变形的情况，即映射函数与像素或图像区域的位置有关，每个位置或区域采用不同的映射函数，如分段线性映射[3]、分段多项式映射[4]，TPS(Thin Plate Spline)映射[5]等。([3]Goshtasby A.Piecewise linear mapping functions for image registration.Pattern Recognition，1986，19：459-466

[4]Goshtasby A.Piecewise cubic mapping functions for image registration.Pattern Recognition，1987，20：525-533

[5]Bookstein F L.Principal warps：Thin-plate splines and the decompositionof deformations.IEEE Trans.Pattern Anal.Mach.Intell.，1989，vol.11，no.6，pp.567-585)

映射模型参数的估计是通过控制点来实现的。对整体变换模型来说，由于控制点的数量通常要远远多于确定模型所需要的最低要求，因此常采用最小二乘拟合的方法来计算模型参数。最小二乘拟合方法使得控制点位置的平方误差和最小，这也意味着，这种方法可以容忍控制点具有一定的误差。经验表明，当控制点数量较多且分布比较均匀时，这种方法是非常有效和非常精确的[6]。而局部映射模型大都采用插值的方法，这类方法要求控制点对应是精确的，而控制点周围像素的映射函数强烈依赖于与邻近控制点之间的关系。确定大量的精确控制点是这类方法的难点所在，目前，在遥感领域，控制点的选择仍然以人工方式为主。

下面，我们首先回顾一下文献中控制点自动匹配方法。

文献中控制点匹配方法很多[6]，但大多数方法都遵循这样一个框架：首先在参考图像和输入图像中分别提取诸如封闭区域、边缘、线段交点、角点等作为图像特征，然后借助于特征的不变量描述或者空间关系确定参考图像和输入图像中的对应特征，将匹配后的特征用代表性的点来表示(如区域重心等)，并将这些点作为对应控制点。但是，由于成像环境存在不可避免的干扰因素，再加上特征点提取算法固有的特点，自动算法提取的控制点位置和实际位置总会存在一定的误差，通常干扰越大，误差越大。在复杂的遥感图像配准中，期望通过这类方法选择的控制点的精度达到亚像素级是不现实的，甚至对一般的视频图像来说，也很难达到这个精度要求。

另外一类控制点选择方法是基于特征相关的。前面提到的方法是分别在参考图像和输入图像中独立提取特征，然后确定这些特征之间的对应关系。而基于特征相关的方法，是只在参考图像中提取特征，而在输入图像中一定范围内采用类似于相关的方法，根据一定的相似性度量标准，搜索参考图像中的特征在输入图像中的最佳对应特征。文献[7](Bentoutou Y，Taleb N，Kpalma K and Ronsin J.An Automatic ImageRegistration for Applications in Remote Sensing.IEEE Transactions onGeoscience and Remote Sensing，2005，43(9)：2127-2137)中首先在参考图像中提取Harris[8](Harris C and Stephens M J.A combined corner andedge detector.In Alvey Vision Conference，1988，pp.147-152)角点，然后采用相似和模糊组合不变矩[9](Flusser J and Suk T.Degraded imageanalysis：An invariant approach.IEEE Trans.Pattern Anal.Mach.Intell.，1998，vol.20，no.6，pp.590-603)作为相似性度量，通过最近邻距离规则确定角点的最佳对应关系。文献[10](Kennedy R E and Cohen W B.Automated designation of tie-points for image-to-image coregistration.International Journal of Remote Sensing，2003，24(17)：3467-3490)首先人工校正图像之间的旋转误差，使得图像之间只存在平移变形。接着，在参考图像中按照一定的步长均匀提取一定大小的网格图像，采用灰度相关方法搜索每个网格在输入图像中的匹配位置，然后，根据相关阈值和空间约束剔除错误匹配对，最后，将网格中心作为控制点来估计模型参数。由于均匀选取的网格图像可能位于平坦区域，这样匹配的结果有时会很不可靠。针对这个问题，文献[11](Chalermwat P.High performanceautomatic image registration for remote sensing[Ph.D.Dissertation].George Mason Univertsity，USA，1999)中首先采用某些度量(梯度、熵和可配准性)来选取纹理丰富的区域，然后只对这些区域采用和文献[10]类似的操作。上述这些方法强烈依赖于相似性度量标准，可能会面临两方面的问题。第一，我们知道，一般的相似性度量是建立在统计最优化基础上的，统计数据(对应于子图像大小)越多，可靠性相对越高。但在遥感图像中，考虑到局部干扰和变形，子区域不宜选得太大，而区域选择太小，可靠性将会降低；第二，相似性度量阈值的选取也很困难。度量标准定得太低，可靠性和精确性都会降低，而度量标准定得太高，可靠性和精度相对提高了，但检测到的特征数量却大大减少了。

还有一类方法是基于整体特征匹配的。前面两类方法都是对每个特征逐个进行匹配，而这类方法是对变换模型空间直接进行搜索，来确定所有特征的整体最优匹配。如[12](Netanyahu N S，Moigne J L and MasekJ G.Georegistration of Landsat Data via Robust Matching of MultiresolutionFeatures.IEEE Trans.Geoscience and Remote Sensing，2004，42(7)：1586-1600)在参考图像和输入图像中分别提取小波极值点，采取分支有界(Branch-and-bound)搜索策略，采用Hausdorff距离[13](Huttenlocher DP，Klanderman G A and Rucklidge W J.Comparing images using theHausdorff distance.IEEE Transactions on Pattern Analysis and MachineIntellinence，1993，15：850-863)作为点集的相似性度量，来估计最优的刚体变换模型，迭代松弛法[14](Ranade S and Rosenfeld A.Point patternmatching by relaxation.Pattern Recognition，1980，12：269-275)也是这类方法的代表。采用这类方法的前提是，要求参考图像和输入图像的特征集大部分是匹配的，而且参考图像和输入图像的特征集满足某个整体映射模型。对不同时相的遥感图像来说，通常场景变化比较大，干扰比较严重，因此从参考图像和输入图像中提取的特征集误差很大，因此很难应用这类方法。而且这类方法的另一个致命缺点是计算代价非常大。

造成控制点定位误差的原因是多种多样的，除了计算过程造成的各种误差外，对多时相影像配准，还存在场景变化问题。这样，通过某些相似性度量方法(甚至手工方法)确定的“高精度”对应点，由于场景的自然变化，实际上其中某些对应点并不能反映场景的真实几何变形，而是由场景的实际变化造成的。毫无疑问，这类点不应该看作是高精度对应点。但这种情况目前还没有好的办法来区分。

通过以上分析可以看出，尽管文献中有很多种控制点自动选取方法，但是，这些方法由于精度和鲁棒性等原因，很难真正应用到实际问题中，目前，遥感影像的控制点选择仍然以手工方式为主。

发明内容

为了实现大面幅遥感影像的高精度配准，通常需要大量的高精度控制点，然而，传统的控制点选择方法，其定位精度很难达到子像素级，为了解决传统技术的问题，本发明的目的是使控制点的定位精度达到子像素级，为此，本发明提出了一种拟合重采样的遥感影像高精度控制点自动选择方法。

本发明遥感影像高精度控制点自动选择方法，包括以下两个主要步骤：

子图像高精度拟合配准；

从已配准子图像中选择控制点。

根据本发明的实施例，所述的子图像高精度配准，包括以下步骤：对子图像粗配准；利用子图像粗配准信息，对子图像特征点检测；利用检测得到的子图像特征点进行匹配；利用匹配结果进行统计拟合变换模型。

根据本发明的实施例，所述的子图像粗配准，采用基于傅立叶Fourier变换的方法估计影像的近似旋转变换参数和近似平移变换参数。

根据本发明的实施例，所述的子图像特征点检测，是首先检测图像边缘点，然后从图像边缘点中根据Harris角点度量函数选择特征点。

根据本发明的实施例，所述的子图像特征点检测，采用高斯差分DOG(Difference of Gaussian，高斯差分)算子的极值点作为特征点。

根据本发明的实施例，所述的子图像特征点匹配，是采用基于特征相关的匹配方法。

根据本发明的实施例，所述的子图像特征点匹配，是采用基于SIFT(Scale Invariant Feature Transform)描述子的匹配方法。

根据本发明的实施例，所述基于特征相关的匹配方法，采用灰度相关系数作为相似性度量标准，在粗配准参数所确定的范围内搜索对应特征点，并且要求灰度相关系数大于某一预先设定的阈值。

根据本发明的实施例，所述采用基于SIFT描述子的匹配方法，将描述子欧氏距离的最近邻确定为对应特征点，而且要求描述子最近邻对应的特征点位于粗配准参数所确定的一定范围内。

根据本发明的实施例，所述统计拟合变换模型，是对已经确定对应特征点，采用最小二乘方法来拟合子图像变换模型。

根据本发明的实施例，所述从已配准子图像中选择控制点，是采用变换模型，从子图像中根据实际应用需要选点作为控制点。

根据本发明的实施例，所述基于特征相关的匹配方法步骤是：

(1)在参考图像f₁中采用Canny算子检测图像边缘；

(2)对每个边缘点计算Harris角点强度S，如果S大于某个给定的阈值H_T，将这个点作为一个特征点；

(3)对所有特征点，采用局部非极大值抑制算法，保留半径为r的范围内角点强度最大的点作为后续匹配的侯选特征点。

(4)采用粗配准方法，计算出输入图像f₂相对于参考图像f₁的近似旋转角度θ₀和平移量(x₀，y₀)，并将输入图像反方向旋转θ₀后得到归整后的输入图像f′₂，在参考图像和归整后的输入图像中仅存在近似的平移偏差(x₀，y₀)；

(5)采用模板匹配方法搜索f₁中每个侯选特征点在f′₂中的匹配点；对f₁中的某个侯选特征点(x_i，y_i)，只需要在f′₂中(x_i+x₀，y_i+y₀)周围很小的邻域R进行；在x和y方向仅取R＝±3像素作为搜索范围。然后，通过对相关函数的表面拟合，将特征点的匹配位置插值到子像素精度，在f₁和f′₂中得到K对对应特征点(p_1i，p′_2i)(i＝1，2，…，K)；

(6)根据第(4)步计算出来的旋转角度θ₀，将归整后输入图像中特征点p′_2i的坐标转换为输入图像坐标p_2i，这样组成了f₁和f₂中特征对(p_1i，p_2i)(i＝1，2，…，K)；

(7)应用所有特征对(p_1i，p_2i)(i＝1，2，…，K)，采用最小二乘近似方法，拟合出最小二乘意义下的子图像配准变换模型参数。

根据本发明的实施例，所述基于SIFT描述子的匹配方法步骤是：

(1)分别在参考图像和输入图像中检测DOG极值点作为特征点，并计算每个特征点的旋转和尺度不变描述子SIFT；

(2)采用基于Fourier变换的粗配准方法，计算出输入图像f₂相对于参考图像f₁的近似旋转角度θ₀和平移量(x₀，y₀)，对参考图像中的任意一点p₁(x，y)，在输入图像中近似对应于p₂(u，v)，其中：

u＝xcosθ₀+ysinθ₀+x₀

v＝-xsinθ₀+ycosθ₀+y₀

(3)对参考图像中的每个特征点p₁(x_i，y_i)，计算它的SIFT描述子与输入图像中所有特征描述子的欧氏距离，假设最近邻距离对应于输入图像中的特征点p₂(u_j，v_j)，而根据上式计算出来的p₁(x_i，y_i)在输入图像中的近似对应位置为p₂(u′_j，v′_j)，如果p₂(u_j，v_j)与p₂(u′_j，v′_j)之间的欧氏距离小于某个设定的阈值R_T，那么，将p₁(x_i，y_i)与p₂(u_j，v_j)作为一对匹配特征点；否则，认为p₁(x_i，y_i)在输入图像中没有对应特征点；

(4)采用非极小值抑制算法，保留半径为r的范围内描述子距离最短的点作为用于后续模型拟合的特征点；

(5)应用所有特征对(p_1i，p_2i)(i＝1，2，…，K)，采用最小二乘近似方法，拟合出最小二乘意义下的子图像配准变换模型参数。

选择遥感影像中一定数量和一定大小的子图像，对每个子图像进行高精度拟合配准，然后从配准的子图像中重新采样任意点，这些点相应地被认为是高精度的，作为整幅影像的控制点，这也即是本发明所指的拟合重采样。从多个子图像中选择出高精度控制点后，后续步骤采用文献中常用的局部映射模型[3-5]来配准完整的大面幅遥感影像。在整个配准过程中，选择控制点是关键步骤。本发明以Landsat影像为例，主要内容是高精度控制点的选择，也可以说，主要是子图像的高精度配准。

具体讲，对子图像的配准，我们采用全局映射模型拟合配准的方法。首先在子图像中提取一定数量且分布较好的特征点(这些特征点允许有一定的误差)，然后用所有这些特征点来拟合全局变换模型，以此来实现子图像的高精度配准。这种图像配准方法依赖于以下假设：第一，模型选择的合理性。对大面幅卫星遥感影像来说，采用全局映射模型通常是不恰当的，但很多情况下，对一定大小的局部子图像来说，采用全局映射模型却是合理的。比如在Landsat影像配准问题中，对局部地形起伏不是很大的地区，采用仿射变换模型或多项式变换模型来描述图像之间的变形被证明是非常精确和有效的；第二，特征点数量应该足够多，而且分布比较均匀。在本发明的具体实现过程中，已经考虑到了这个问题。

当模型选择合理，特征点充足时，这种配准方法可以容忍特征点存在一定的误差，但最终的配准精度能达到亚像素级。相应地，子图像内每个像素点也都达到了亚像素级配准精度。这样我们就可以从已经配准子图像中任意选点作为高精度控制点。

这种方法的特点在于，最终的控制点不是采用传统的方法来直接检测，而是采用间接的方法来实现的，其最大的优点是避免了传统方法提取高精度控制点的困难，同时可以方便地调整最终控制点的分布，比如在特征稀少的区域控制点采样可以密集些，而在特征密集的区域可以采样稀疏些。

在这个整体构架下，本发明中采取先粗配后精配的策略，对两种具体的配准方法进行了实验比较，取得了比较满意的实验结果。从最终实验结果看，本发明提出的拟合重采样思想和先粗配后精配的具体实现策略可以有效和可靠地实现子图像的高精度配准。

附图说明

图1是利用本发明对Landsat影像进行配准的结果(场景1)

图2是利用本发明对Landsat影像进行配准的结果(场景2)

图3是利用本发明对Landsat影像进行配准的结果(场景3)

图4是利用本发明对Landsat影像进行配准的结果(场景4)

图5是利用本发明对Landsat影像进行配准的结果(场景5)

图6是利用本发明对Landsat影像进行配准的结果(场景6)

具体实施方式

下面将结合附图对本发明加以详细说明，应指出的是，所描述的实施例仅旨在便于对本发明的理解，而对其不起任何限定作用。

本发明中的研究对象为经过系统校正(Level 1G)后的Landsat影像。根据卫星导航数据进行系统校正后的Landsat影像定位误差一般在几个像素范围。例如，对Landsat-5数据，导航误差在3～30像素，对Landsat-7数据，为1～3像素。另外，还可能有微小的旋转误差。

不同时相的Landsat成像具有较好的几何一致性，而且本发明的实验场景中，地形起伏属于中等，因此局部地区可采用仿射变换模型或刚体变换(旋转+平移)模型来描述不同影像之间的变形。

下面先介绍粗配准方法。

1基于Fourier变换的粗配准

所述的子图像粗配准，采用基于傅立叶变换的方法估计影像的近似旋转变换参数和近似平移变换参数。

作为相位相关法[15](Kuglin C D and Hines D C.The phasecorrelation image alignment method.In Proceedings of the IEEE 1975International Conference on Cybernetics and Society，New York，1975，pp.163-165)的扩展，De Castro和Morandi[16](Castro E D and Morandi C.Registration of translated and rotated images using finite Fourier transform.IEEE Transactions on Pattern Analysis and Machine Intelligence，1987(9)：700-703)曾提出了一种基于Fourier变换的平移和旋转图像配准方法。这种方法分两个步骤来完成，第一步先计算旋转角度，然后再确定平移量。

考虑被配准的两幅图像f₁(x，y)和f₂(x，y)，其中f₂(x，y)是f₁(x，y)经过平移(x₀，y₀)和旋转θ₀后的图像，即

f₂(x，y)＝f₁(xcosθ₀+ysinθ₀-x₀，-xsinθ₀+cosθ₀-y₀)

根据Fourier变换的平移和旋转特性，f₁(x，y)和f₂(x，y)的Fourier变换之间的关系为：

$F_2(\mathrm{\ξ},\mathrm{\η})=e^{-j2\mathrm{\π}({\mathrm{\ξx}}_0+\mathrm{\η}{y}_0)}{F}_1(\mathrm{\ξ}\cos {\mathrm{\θ}}_0+\mathrm{\η}\sin {\mathrm{\θ}}_0,-\mathrm{\ξ}\sin {\mathrm{\θ}}_0+\mathrm{\η}\cos {\mathrm{\θ}}_0)$

若M₁和M₂代表F₂和F₂的幅度谱，那么，从上式可以得到：

M₂(ξ，η)＝M₁(ξcosθ₀+ηsinθ₀，-ξsinθ₀+ηcosθ₀)

如果从F₁和F₂的幅度谱考虑，从上式很容易看出，两个频谱幅度是相同的，只是M₂相对于M₁旋转了θ₀。将幅度谱转化到极坐标空间后，旋转位移可以简化为类似于平移的表达形式：

M₁(ρ，θ)＝M₂(ρ，θ-θ₀)

这样就可以采用经典的相位相关法计算出角度θ₀。

根据计算出来的角度θ₀，将f₂(x，y)旋转-θ₀得到归整后的图像f′₂(x，y)，这时，在f₁(x，y)和f′₂(x，y)之间仅有平移误差，可在图像空间再一次应用相位相关法求得平移量。这样，经过两次相位相关运算，就求出了旋转角度θ₀和平移量(x₀，y₀)。

大量实验表明，基于Fourier变换的配准方法，其角度计算误差与分度间隔有关，一般不会超过1～2个分度间隔，平移误差一般不超过1～2个像素范围。

从后面的讨论我们将会看到，首先计算粗配准参数，对后续精确匹配过程具有重要的意义。下面介绍实验中用到的两种配准方法。

2基于特征相关的匹配方法

所述的子图像特征点匹配，是采用基于特征相关的匹配方法。

基于特征相关的图像配准方法，其基本思想是：只在参考图像中提取特征，而在输入图像中一定范围内采用类似于模板匹配的方法，根据某种相似性度量标准，搜索参考图像中的特征在输入图像中的对应特征。由于我们的实验对象是不同时相的同波段影像，因此我们采用了基于灰度相关的相似性度量。

所述基于特征相关的匹配方法，采用灰度相关系数作为相似性度量标准，在粗配准参数所确定的范围内搜索对应特征点，并且要求灰度相关系数大于某一预先设定的阈值。

在基于灰度相关的相似性度量标准中，相关系数被证明是比较鲁棒的，而且相关系数对灰度的线性变化和加性高斯噪声不太敏感。

灰度相关系数定义为：

$C_{\mathrm{coef}}=\frac{\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}(f_1(x_i,y_j)-{\mathrm{\μ}}_1)(f_2(x_i,y_j)-{\mathrm{\μ}}_2)}{\sqrt{\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{(f_1(x_i,y_j)-{\mathrm{\μ}}_1)}^2}\sqrt{\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{(f_2(x_i,y_j)-{\mathrm{\μ}}_2)}^2}}$

首先在参考图像中提取特征点集，对每个特征点P_i，以P_i为中心，裁剪大小为M×N像素的一个窗口，称为模板图像f₁，窗口中的像素的灰度均值为μ₁，然后在输入图像中(一定搜索范围R内)搜索模板图像的最佳匹配位置。在每个搜索位置，裁剪和模板图像一样大小的浮动图像f₂，灰度均值为μ₂，计算f₁和f₂的灰度相关系数C_coef。当所有位置搜索完毕，如果最大相关系数大于某一设定的阈值T_c，那么将最大相关系数对应的浮动图像中心P_j作为P_i的匹配位置。

为了提高匹配的可靠性，一般希望选择的模板图像具有比较丰富的纹理，而不是灰度变化比较平坦的区域。基于这样的考虑，我们采用了如下特征点选择算法：

所述的子图像特征点检测，是首先检测图像边缘点，然后从图像边缘点中根据Harris角点度量函数选择特征点。

一般认为，图像边缘包含了丰富的信息，对应于梯度变化较大的区域。因此，我们的选择策略是从图像边缘点中提取特征点。这里采用了Canny边缘检测算法[17](Canny J.A computational approach to edgedetection，IEEE Transactions on Pattern Analysis and Machine Intelligence，1986(8)：679-698)。

当边缘检测完成后，取边缘点中具有较大“强度”的那些点作为特征点。通常，点的强度是根据像素灰度值沿两个方向的变化来测量的。Schmid[18](Schmid C，Mohr R and Bauckhage C.Evaluation of interestpoint detectors.International Journal of Computer Vision，2000，37(2)：151-172)曾对不同的特征点检测算法进行了评价和比较，得出的结论是，Harris算法[8]性能最优，而且信息熵最高(意味着纹理信息比较丰富，后续的匹配可靠性较高)，因此，本发明采用了Harris算法提取特征点。

概括起来，基于特征相关的子图像配准方法可归纳为以下步骤。

(1)在参考图像f₁中采用Canny算子检测图像边缘。

(2)对每个边缘点计算Harris角点强度S，如果S大于某个给定的阈值H_T，那么将这个点作为一个特征点。

(3)对所有特征点，采用局部非极大值抑制算法，仅保留半径为r的范围内角点强度最大的点作为后续匹配的侯选特征点。

(4)首先采用粗配准方法，计算出输入图像f₂相对于参考图像f₁的近似旋转角度θ₀和平移量(x₀，y₀)，并将输入图像反方向旋转θ₀后得到归整后的输入图像f′₂，这时，在参考图像和归整后的输入图像中仅存在近似的平移偏差(x₀，y₀)。

(5)采用模板匹配方法搜索f₁中每个侯选特征点在f′₂中的匹配点。对f₁中的某个侯选特征点(x_i，y_i)，只需要在f′₂中(x_i+x₀，y_i+y₀)周围很小的邻域R进行。我们在x和y方向仅取R＝±3像素作为搜索范围。然后，通过对相关函数的表面拟合，将特征点的匹配位置插值到子像素精度。最后，在f₁和f′₂中得到K对对应特征点(p_1i，p′_2i)(i＝1，2，…，K)。

(6)根据第(4)步计算出来的旋转角度θ₀，将归整后输入图像中特征点p′_2i的坐标转换为输入图像坐标p_2i，这样组成了f₁和f₂中特征对(p_1i，p_2i)(i＝1，2，…，K)。

(7)所述统计拟合变换模型，是对已经确定对应特征点，采用最小二乘方法来拟合子图像变换模型。

应用所有特征对(p_1i，p_2i)(i＝1，2，…，K)，采用最小二乘近似方法，拟合出最小二乘意义下的子图像配准变换模型参数(本发明对应仿射变换模型)。

3基于SIFT描述子的特征匹配方法

所述的子图像特征点匹配，是采用基于描述子SIFT的匹配方法。

这种方法沿用了基于特征的图像配准一般框架：首先在参考图像和输入图像中分别独立提取特征点，然后采用基于不变量描述子的相似性度量方法进行特征点匹配。

SIFT描述子[19](Lowe D G.Distinctive image features fromscale-invariant keypoints.International Journal of Computer Vision，2004，60(2)，91-110)，是图像匹配和目标识别等应用中性能最为杰出的描述子之一[20](Mikolajczyk K and Schmid C.A performance evaluation oflocal descriptors.IEEE Trans.Pattern Analysis and Machine Intelligence，2005，27(10)：1615-1630)，SIFT描述子在干扰和变形比较严重的遥感图像配准中同样具有很强的鲁棒性和很强的特征辨析性能。

为了和前面讲的基于特征相关的子图像配准方法进行比较，我们对基于SIFT特征的配准方法也进行了实验。

基于SIFT描述子的子图像配准方法可以归纳为以下步骤：

(1)所述的子图像特征点检测，采用高斯差分DOG算子的极值点作为特征点。

采用文献[19]中的方法，分别在参考图像和输入图像中检测DOG(Difference ofGaussian)极值点作为特征点，并计算每个特征点的旋转和尺度不变SIFT描述子。

(2)采用基于Fourier变换的粗配准方法，计算出输入图像f₂相对于参考图像f₁的近似旋转角度θ₀和平移量(x₀，y₀)，那么，对参考图像中的任意一点p₁(x，y)，在输入图像中近似对应于p₂(u，v)，其中：

u＝xcosθ₀+ysinθ₀+x₀

v＝-xsinθ₀+ycosθ₀+y₀

(3)所述采用基于描述子SIFT的匹配方法，将描述子欧氏距离的最近邻确定为对应特征点，而且要求描述子最近邻对应的特征点位于粗配准参数所确定的一定范围内。

对参考图像中的每个特征点p₁(x_i，y_i)，计算它的SIFT描述子与输入图像中所有特征描述子的欧氏距离，假设最近邻距离对应于输入图像中的特征点p₂(u_j，v_j)，而根据上式计算出来的p₁(x_i，y_i)在输入图像中的近似对应位置为p₂(u′_j，v′_j)，如果p₂(u_j，v_j)与p₂(u′_j，v′_j)之间的欧氏距离小于某个设定的阈值R_T，那么，将p₁(x_i，y_i)与p₂(u_j，v_j)作为一对匹配特征点。否则，认为p₁(x_i，y_i)在输入图像中没有对应特征点。

需要注意的是，这里我们没有对描述子的最近邻距离设定阈值，一方面是因为SIFT描述子是一个128维的高维向量，最近邻距离会相差很大，很难设定某个固定阈值；另外一个原因是，鉴于SIFT描述子的高辨析能力，如果在众多(通常在400×400大小的图像中可以检测到数百到数千DOG极值点)待比较的特征点中，其最近邻距离对应的特征点位于预先设定的比较小的范围内的话，那么这也间接地说明，正确匹配的可能性是很大的。

(4)采用类似于非极大值(这里是非极小值)抑制算法，仅保留半径为r的范围内描述子距离最短的点作为用于后续模型拟合的特征点。

(5)所述统计拟合变换模型，是对已经确定对应特征点，采用最小二乘方法来拟合子图像变换模型。

应用所有特征对(p_1i，p_2i)(i＝1，2，…，K)，采用最小二乘近似方法，拟合出最小二乘意义下的子图像配准变换模型参数。

从上面介绍的两种子图像配准方法来看，子图像的粗配准具有重要的意义，这体现在：第一，对基于特征相关的方法来说，减小了相关搜索范围，提高了计算效率；第二，对基于描述子的方法来说，提高了匹配可靠性；第三，也是最重要的，无论相关方法还是描述子方法，这样确定的特征对应点的误差被控制在一定的范围内，即对应点最大误差不会超过预先设定的搜索范围，这在一定程度上保证了后续模型拟合的精确性和可靠性。

所述从已配准子图像中选择控制点，是采用变换模型，从子图像中根据实际应用需要选点作为控制点。

当采用基于特征相关的方法或者基于SIFT描述子的匹配方法拟合出子图像变换模型后，可根据应用需要，利用变换模型，从子图像中选择适当的对应点作为大面幅遥感影像的控制点。

4实验结果及讨论

在介绍实验内容之前，我们先讨论一下图像配准的精度评价问题。

4.1图像配准的精度评价概述

对于图像配准问题来说，提供一个最终的配准误差估计是很有必要的，否则在实际应用中很难被用户接受。但通常情况下由于很难得到真实数据，因此，这个问题并不是一个很简单的问题。这里我们介绍几个与我们后续实验有关的误差评价方法，尽管这些方法无法可靠地反映最终的配准误差，至少也可以提供一些分析的依据。

配准误差常采用以下几种方法来评价。

●最简单的度量方法是计算控制点(CPs)的均方误差(CPE)。设参考图像中的控制点集为P₁(x_i，y_i)，对应的输入图像中的控制点集为P₂(u_i，v_i)，i＝1，2，…，K，K表示控制点的数量。假设通过所有控制点计算出来的映射模型为u＝T₁(x，y)，v＝T₂(x，y)，那么，CPE定义为：

$\mathrm{CPE}=\sqrt{\frac{\underset{i=1}{\overset{K}{\mathrm{\Σ}}}{(u_i-T_1(x_i,y_i))}^2+{(v_i-T_2(x_i,y_i))}^2}{K}}$

尽管CPE是常用的误差估计方法，但并不是一个很好的误差估计度量。实际上，CPE仅仅能度量CPs与映射模型的拟合程度，而映射模型本身就是通过CPs计算出来的，这样的度量不能反映图像的整体情况。另外CPE较大，可能是由于CP定位误差引起的，未必真正反映了配准精度很差。

●与CPE类似的另一种方法常称作测试点误差(TPE)，这种方法是只用控制点集中的一部分点来计算映射模型，而用点集中其余控制点计算均方误差来度量统计误差。这种方法似乎比CPE更有意义。然而，确切地说，控制点的误差与图像的配准误差是不同的两个概念，很多情况下不具有可比性。

●三角验证：对同一场景至少有三幅图像的情况，可以采用一致性检查方法。对三幅图像的情况，首先根据控制点计算两两图像之间的映射函数T_ij，假设f₂＝T₁₂(f₁)，f₃＝f₂₃(f₂)，那么，根据模型传递关系，f₃＝T₂₃(T₁₂(f₁))，另外，f′₃＝T₁₃(f₁)。从理论上来说，f₃＝f′₃，这样，通过度量f₃与f′₃的误差来验证映射模型的正确性。f₃与f′₃的误差可通过选择一些测试点进行度量。例如，在我们的实验中，在图像f₁中沿两个方向每隔10个像素选择一个测试点，最后来测试这些点在f₃与f′₃的误差，通过计算均方误差来度量这些测试点的整体误差。

●相互验证：采用多种方法进行一致性检查。用不同特点的方法分别配准，如果这两种方法的差别很小，那么可认为该方法具有较好的配准精度(尽管不能严格保证)。我们对前面提到的两种配准方法进行了相互验证。但是，如果两种方法差别比较大，那么无法确定到底哪种方法不准确或两种方法都不准确。

4.2参数设置及实验结果

(1)参数设置

①Canny边缘检测算法采用Matlab中图像处理工具箱中标准函数的缺省设置。

②Harris角点检测算法中，高斯函数标准差为1，角点度量的阈值H_T设为1000；

③特征点局部非极大值抑制半径设为5象素；

④模板图像大小设置为50×50像素；

⑤模板匹配的搜索范围R设置为±3像素；

⑥灰度相关系数阈值T_c设为0.7；

⑦基于SIFT描述子方法中，描述子最近邻点的位置误差阈值R_T设为5像素。

(2)模型参数

实验中，子图像映射模型采用仿射变换。第一幅图像为参考图像，像素点坐标表示为(x，y)，第二幅图像为输入图像，像素点坐标表示为(u，v)，模型变换关系表示为：

u＝a₁₁x+a₁₂y+t_x

v＝a₂₁x+a₂₂y+t_y

后面统计表中，“模型参数”中每一项代表上述对应变换的模型系数。

(3)实验数据说明

实验内容共安排了6个场景，场景图像大小均为400×400像素。场景1，2，3和6中，每个场景包括三幅不同时间获取的图像。这些子图像是从三幅完整Landsat影像中裁剪的，这三幅图像对应于河北承德某地区，是分别于1996年、1997年和2002年获取的未经配准的原始影像。场景4和5分别包括两幅图像，这些图像是从另外两幅完整Landsat影像中裁剪的，这两幅影像是分别于2001年5月和2001年7月获得的，并且已经由专业人员采用手工选择控制点的方法进行了配准。

对具有三幅图像的场景，计算两两图像之间的配准数据。

图表中有关说明如下：“场景2-1，2-2，2-3”指场景2的第一幅(1996年)、第二幅(1997年)和第三幅(2002年)图像；对两幅图像的场景，前一幅对应于2001年5月，后一幅对应于2001年7月；“Corr”和“SIFT”分别指前面提到的两种配准方法；表中“CPE”是采用所有特征点计算变换模型，用所有特征点作为测试点计算均方根误差，“CPE_Max”表示测试点的最大均方根误差，Ratiol表示均方根误差大于1个像素的测试点占所有测试点的比例。

“TPE”是采用随机选取的一半点计算变换模型，而用另一半点进行测试的均方根误差，“TPE_Max”表示测试点的最大均方根误差，Ratio2表示均方根误差大于1个像素的测试点占所有测试点的比例。

表中“两种方法相互验证”是指，在输入图像(第一幅图像)中沿x和y方向每隔10个像素选择一个测试点(对400×400的图像，共选取39×39＝1521个测试点)，构成一个测试点集P，假设两种方法求出的变换模型分别为T1和T2，那么表中“RMS2”对应于点集T₁(P)和T₂(P)的均方根误差。

表中“三角验证”是指，从第一幅图像中采取和以上相同的方式均匀选取1521个点作为测试点集P，假设第一幅图像到第二幅图像、第二幅图像到第三幅图像，第一幅图像到第三幅图像的变换模型分别为T₁₂、T₂₃和T₁₃，那么表中“RMS3”对应于点集T₂₃(T₁₂(P))和点集T₁₃(P)的均方根误差。

表4中，“和理论数据比较”是指，已知理论变换模型为T₁₂，实际计算出来的模型为T′₁₂，采用以上相同的测试点集P，“RMS4”对应于点集T₁₂(P)和点集T′₁₂(P)的均方根误差。

图表中其它内容可根据以上说明类推。

4.3实验结果讨论

(1)从所有的表格中首先可以看到，采用Fourier变换计算出来的旋转角度近似为0°，说明Landsat影像经过系统校正后的旋转误差是很小的，其次，将x和y方向的近似平移量和最终计算出来的变换模型的平移分量进行比较发现，通过Fourier变换计算出来的粗配准参数误差基本都在1个像素范围内，说明基于Fourier变换的方法是鲁棒的，并且具有较高的精度。

(2)前三个场景的实验结论基本相同，参看表1-表3。从“两种方法相互验证”和“三角验证”的统计结果看，无论是均方误差还是最大均方误差均小于1个像素，而且，从图1～图3的特征点分布看，数量多，分布比较好，因此我们完全有理由认为子图像配准达到了亚像素精度。

如图1所示：是利用本发明对Landsat影像进行配准的结果(场景1)，图中：

最上面一行为同一场景不同时间的三幅Landsart 3波段影像，分别于1996年，1997年和2002年获得，大小均为400×400像素。中间和最下面一行分别为基于特征相关和基于SIFT的配准方法的对应特征点显示。图中“1-1Vs.1-2”是指场景1-1和场景1-2的匹配对应点在场景1-1中的显示，类似，“1-2Vs.1-3”指场景1-2和场景1-3的匹配对应点在场景1-2中的显示。括号中数字，是采用非极大值抑制后的特征点数量，这里的抑制半径为10个像素，和后面表中数据不同，表中的抑制半径为5像素。

如图2所示：是利用本发明对Landsat影像进行配准的结果(场景2)，图中：

最上面一行为同一场景不同时间的三幅Landsart 3波段影像，分别于1996年，1997年和2002年获得，大小均为400×400像素。中间和最下面一行分别为基于特征相关和基于SIFT的配准方法的对应特征点显示。图中“2-1Vs.2-2”是指场景2-1和场景2-2的匹配对应点在场景2-1中的显示，类似，“2-2Vs.2-3”指场景2-2和场景2-3的匹配对应点在场景2-2中的显示。括号中数字，是采用非极大值抑制后的特征点数量，这里的抑制半径为10个像素，和后面表中数据不同，表中的抑制半径为5像素。

如图3所示：是利用本发明对Landsat影像进行配准的结果(场景3)，图中：

最上面一行为同一场景不同时间的三幅Landsart 3波段影像，分别于1996年，1997年和2002年获得，大小均为400×400像素。中间和最下面一行分别为基于特征相关和基于SIFT的配准方法的对应特征点显示。图中“3-1Vs.3-2”是指场景3-1和场景3-2的匹配对应点在场景3-1中的显示，类似，“3-2Vs.3-3”指场景3-2和场景3-3的匹配对应点在场景3-2中的显示。括号中数字，是采用非极大值抑制后的特征点数量，这里的抑制半径为10个像素，和后面表中数据不同，表中的抑制半径为5像素。

但是，如果从常用的CPE和TPE的统计结果看，效果却很差，尤其是基于SIFT的方法结果更差。下面我们进一步分析一下这些结果。

第一，从特征点的定位精度看，大量的特征点定位误差大于1个像素。从CPE看，定位误差大于1个像素的比例最高达到47％，单个特征点最大误差达到差不多5.6像素，而且测试点的均方根误差很多也是大于1个像素的。然而，最终计算出的模型却是高精度的。这也说明，单从特征点的定位精度来判断最终模型的精确程度并不总是合理的；第二，CPE和TPE的统计结果似乎差别不大，这也一定程度上反映了结果的正确性。通常情况下，当模型计算不正确时，一般TPE比CPE要大很多，在这里，两者相差不多，还有好几处TPE比CPE还要小，这也说明，通过部分对应点计算出来的模型和通过所有对应点计算出来的模型相差无几；第三，基于SIFT的特征点定位误差比基于相关的方法要大，这通常也是基于特征的方法的一个普遍特点。但是，表中反映的数据误差也不能完全反映真实的差别。回想前面的参数设置，我们对相关搜索范围设置为±3像素，而将SIFT方法的最近邻匹配位置误差设置为5像素，这也是SIFT方法误差看上去比相关方法误差要大很多的原因之一；第四，由于搜索范围设置得比较小(当然要归功于可靠的粗匹配)，保证了最终的对应点误差不会超过一定的范围，有利于后续匹配的可靠性；第五，从特征点分布看，在这三个场景中，特征点的数量都比较多，而且分布也比较好。

(3)场景4和场景5将实验结果与“真实”结果进行了比较。在场景4和场景5中，每一对图像已由专业人员采用手工选取控制点的方法进行了配准，也即参考图像和输入图像的坐标理论上是完全对应的(不考虑配准误差时)。已经知道了两幅图像之间的“真实”变换模型，(a₁₁，a₁₂，t_x，a₂₁，a₂₂，t_y)＝(1，0，0，0，1，0)，我们试图来验证该方法的正确性。

如图4所示：是利用本发明对Landsat影像进行配准的结果(场景4)，图中：

上面一行为同一场景不同时间的两幅Landsart 3波段影像，分别于2001年5月和2001年7月获得，并且已经经过手工配准，大小均为400×400像素。下面左右图分别为基于特征相关和基于SIFT的配准方法的对应特征点显示。括号中数字，是采用非极大值抑制后的特征点数量，这里的抑制半径为10个像素，和后面表中数据不同，表中的抑制半径为5像素。

如图5所示：是利用本发明对Landsat影像进行配准的结果(场景5)，图中：

上面一行为同一场景不同时间的两幅Landsart 3波段影像，分别于2001年5月和2001年7月获得，并且已经经过手工配准，大小均为400×400像素。下面左右图分别为基于特征相关和基于SIFT的配准方法的对应特征点显示。括号中数字，是采用非极大值抑制后的特征点数量，这里的抑制半径为10个像素，和后面表中数据不同，表中的抑制半径为5像素。

表4中，“和理论模型相比较”是指，将表中方法计算出来的模型与理论模型采用和前面相同的测试点集来计算误差。从“两种方法相互比较”和“和理论模型相比较”可看出，最终配准确实达到了子像素级精度。但同时也发现，在“和理论模型相比较”中，误差有些偏大，这个结果似乎稍微有点奇怪。和前面三个场景相比较，这些图像之间相隔时间较短(只有大约两个月)，噪声也相对较小，另外，特征点数量很多，分布也很好，从理论上说，计算误差应该比场景1，2，3小一些，但实际情况不是如此，我们推测，是“真实”模型不精确所致。或许，我们可以认为，我们提出的自动方法比手工方法更精确。

综合以上分析并结合图表数据可以看到，特征点大都具有一定的定位误差，甚至误差还比较大(几个象素)，而且大于1个象素误差的特征点还占有相当的比例。我们不知道造成这些误差的具体原因，或许是计算过程的各种因素造成的，或许是场景真实变化造成的，但值得庆幸的是，最终的配准精度都达到了亚像素级精度。另外，我们提出的先粗配准，然后采用小范围搜索的策略来确定特征点对应，有以下几个主要优点：第一，有效地增加了特征点数量。比如，当我们限定搜索范围后，对于基于特征相关的配准方法，可以适当降低相关阈值，这样既可以增加特征点数量，而且特征点的定位误差可以限定在一定的范围内，而不影响后续模型的估计；第二，对相关运算来说，减小搜索范围，意味着运算量的减少；第三，鲁棒性强，可靠性高。从场景图像看，各种噪声还是比较严重的，但无论是粗配准还是后续的精配准，都具有很强的适应性。这些特点很适合于复杂的遥感图像。

(4)当然，当各种噪声大到一定程度后，这些方法得出的配准结果将无法满足子像素精度要求，场景6就是一个例子。

如图6所示：是利用本发明对Landsat影像进行配准的结果(场景6)，图中：

最上面一行为同一场景不同时间的三幅Landsart 3波段影像，分别于1996年，1997年和2002年获得，大小均为400×400像素。中间和最下面一行分别为基于特征相关和基于SIFT的配准方法的对应特征点显示。图中“6-1 Vs.6-2”是指场景6-1和场景6-2的匹配对应点在场景6-1中的显示，类似，“6-2Vs.6-3”指场景6-2和场景6-3的匹配对应点在场景6-2中的显示。括号中数字，是采用非极大值抑制后的特征点数量，这里的抑制半径为10个像素，和后面表中数据不同，表中的抑制半径为5像素。

首先我们从特征点的分布图图6看，特征点比前面5个场景明显少，而且特征点的分布也不好，再从表5中的统计数据看，“两种方法相互验证”中误差明显增大，根据这些信息，凭经验我们知道，这两种方法的配准结果可能都很难达到子像素精度。但一个有趣的现象是，从表6“三角验证”中可以看出，采用特征相关方法无论是均方根误差还是最大误差都小于1个像素，似乎说明这种方法达到了子像素级精度，即使观察SIFT方法，也会发现，大于1个像素误差的测试点也不是很多(约9％)，怎么解释这个现象？再次观察特征点分布图，可以发现，两种方法都有一个共同的现象，尽管分布不均匀，但两两图像之间的特征点分布有一定的相似性，这正是原因所在。特征点分布的相似性造成了模型计算的系统误差，导致的结果是，尽管绝对误差较大，但它们之间的相对误差却较小，这样就给出了虚假的“好”结果。这也说明，三角验证方法也并不总是可靠的。

5小结

在卫星遥感影像配准中，确定高精度控制点是一个很重要的问题，传统的控制点直接匹配方法很难达到子像素级精度。本发明中我们提出了一种和传统方法完全不同的思想，即文中所说的“拟合重采样”思想。在这个框架下，采用传统方法确定的图像对应点允许有一定的误差，但最终选择的控制点却可以达到子像素精度。这对于退化原因复杂，干扰比较严重的遥感图像配准来说不失为一个新的解决方案。

在具体实现过程中，我们提出的先粗后精的配准策略较好地保证了最终结果的可靠性，既尽可能多地保证了对应点的数量，同时也减少了运算量。

实验结果验证了该方法的可行性。另外，文中给出的实验采用的都是仿射变换模型，我们也曾采用二次多项式模型进行了实验，得出的结论和仿射变换模型基本相同。这也说明，在这个框架下，根据实际问题，可以采用不同的全局变换模型，也可以采用其它的具体实现方法。

其次，从实验还可以看到，子图像的精确配准强烈依赖于图像对应点的数量和分布。尽管在我们的具体实现中已经采取了缩小搜索范围、降低阈值的策略试图找到尽可能多的对应点，但是并非任何情况下都能得到保证。例如，对于灰度变化比较平缓的区域，特征点数量很少；再比如，当图像干扰比较严重时，很难找到大量的图像对应点(如实验6)。对于这些情况，我们认为：从完整的大型遥感图像配准来说，如果某些子图像中的特征点无论数量还是分布不能满足一定要求的话，可以放弃从这些子图像中选择控制点，而从与这些子图像相邻的已经精确配准的子图像中采取比较密集的控制点采样方法来补偿局部控制点不足的问题；从另外一个角度讲，如果一个子图像中部分区域对应点密集而另外一些区域特征点稀少，这样拟合出来的变换模型可能不能保证子图像中所有点都能达到子像素精度，但通常对应点密集的区域拟合精度会很高，而特征密集的区域通常对应于重要的区域，这样处理也是一种折衷的方法，即特征密集区域控制点定位精度高，而特征稀少的区域控制点定位精度稍低，这在某些具体应用中或许也能满足应用需求。

最后需要说明的是，在自动实现过程中，对每个子图像的配准结果进行可靠的评价是很重要的，我们研究了多种误差估计方法，但是可以看出，没有哪种方法能完全可靠地对结果做出评价，一种可能的方法是采用多个线索来综合判断，但这样会增大计算量。但就目前的情况看，除此之外，从一般意义上讲，似乎还没有更好的方法。

表1：Landsat影像模型参数估计及误差统计(场景1)

表2：Landsat影像模型参数估计及误差统计(场景2)

表3：Landsat影像模型参数估计及误差统计(场景3)

表4：Landsat影像模型参数估计及误差统计(场景4和场景5)

表5：Landsat影像模型参数估计及误差统计(场景6)

上面描述是用于实现本发明及其实施例，因此，本发明的范围不应由该描述来限定。本领域的技术人员应该理解，在不脱离本发明的范围的任何修改或局部替换，均属于本发明权利要求来限定的范围。

Claims (12)

1、一种遥感影像高精度控制点自动选择方法，其特征在于：包括以下两个主要步骤：

选择遥感影像中一定数量和一定大小的子图像；对子图像粗配准；利用子图像粗配准信息，对子图像特征点检测；利用检测得到的子图像特征点进行匹配；利用匹配结果进行统计拟合变换模型；对每个子图像高精度拟合配准；

从已配准子图像中重新采样任意点，这些任意点作为整幅影像的选择控制点，采用局部映射模型来配准完整的大面幅遥感影像。

2、按权利要求1所述的遥感影像高精度控制点自动选择方法，其特征在于：所述的子图像粗配准，采用基于傅立叶变换的方法估计影像的近似旋转变换参数和近似平移变换参数。

3、按权利要求1所述的遥感影像高精度控制点自动选择方法，其特征在于：所述的子图像特征点检测，是首先检测图像边缘点，然后从图像边缘点中根据Harris角点度量函数选择特征点。

4、按权利要求1所述的遥感影像高精度控制点自动选择方法，其特征在于：所述的子图像特征点检测，采用高斯差分DOG算子的极值点作为特征点。

5、按权利要求1所述的遥感影像高精度控制点自动选择方法，其特征在于：所述的子图像特征点匹配，是采用基于特征相关的匹配方法。

6、按权利要求1所述的遥感影像高精度控制点自动选择方法，其特征在于：所述的子图像特征点匹配，是采用基于描述子SIFT的匹配方法。

7、按权利要求5所述的遥感影像高精度控制点自动选择方法，其特征在于：所述基于特征相关的匹配方法，采用灰度相关系数作为相似性度量标准，在粗配准参数所确定的范围内搜索对应特征点，并且要求灰度相关系数大于某一预先设定的阈值。

8、按权利要求6所述的遥感影像高精度控制点自动选择方法，其特征在于：所述采用基于描述子SIFT的匹配方法，将描述子欧氏距离的最近邻确定为对应特征点，而且要求描述子最近邻对应的特征点位于粗配准参数所确定的一定范围内。

9、按权利要求1所述的遥感影像高精度控制点自动选择方法，其特征在于：所述统计拟合变换模型，是对已经确定对应特征点，采用最小二乘方法来拟合子图像变换模型。

10、按权利要求1所述的遥感影像高精度控制点自动选择方法，其特征在于：所述控制点，是采用变换模型，从子图像中根据实际应用需要选点作为控制点。

11、按权利要求5或7所述的遥感影像高精度控制点自动选择方法，其特征在于：所述基于特征相关的子图像配准步骤是：

(1)在参考图像f₁中采用算子Canny检测图像边缘；

(2)对每个边缘点计算Harris角点强度S，如果S大于某个给定的阈值H_T，将这个点作为一个特征点；

(3)对所有特征点，采用局部非极大值抑制算法，保留半径为r的范围内角点强度最大的点作为后续匹配的侯选特征点；

(4)采用粗配准方法，计算出输入图像f₂相对于参考图像f₁的近似旋转角度θ₀和平移量(x₀，y₀)，并将输入图像反方向旋转θ₀后得到归整后的输入图像f′₂，在参考图像和归整后的输入图像中仅存在近似的平移偏差(x₀，y₀)；

(5)采用模板匹配方法搜索f₁中每个侯选特征点在f′₂中的匹配点；对f₁中的某个侯选特征点(x_i，y_i)，在f′₂中(x_i+x₀，y_i+y₀)周围很小的邻域R进行；在x和y方向仅取R＝±3像素作为搜索范围；然后，通过对相关函数的表面拟合，将特征点的匹配位置插值到子像素精度；在f₁和f′₂中得到K对对应特征点(p_1i，p′_2i)(i＝1，2，…，K)；

(6)根据第(4)步计算出来的旋转角度θ₀，将归整后输入图像中特征点p′_2i的坐标转换为输入图像坐标p_2i，组成f₁和f₂中特征对(p_1i，p_2i)(i＝1，2，…，K)；

(7)应用所有特征对(p_1i，p_2i)(i＝1，2，…，K)，采用最小二乘近似方法，拟合出最小二乘意义下的子图像配准变换模型参数。

12、按权利要求6或8所述的遥感影像高精度控制点自动选择方法，其特征在于：所述基于SIFT特征的子图像配准方法步骤是：

(1)分别在参考图像和输入图像中检测高斯差分DOG极值点作为特征点，并计算每个特征点的旋转和尺度不变描述子SIFT；

(2)采用基于傅立叶Fourier变换的粗配准方法，计算出输入图像f₂相对于参考图像f₁的近似旋转角度θ₀和平移量(x₀，y₀)，对参考图像中的任意一点p₁(x，y)，在输入图像中近似对应于p₂(u，v)，其中：

u＝xcosθ₀+ysinθ₀+x₀

v＝-xsinθ₀+ycosθ₀+y_o

(3)对参考图像中的每个特征点p₁(x_i，y_i)，计算它的描述子SIFT与输入图像中所有特征描述子的欧氏距离，假设最近邻距离对应于输入图像中的特征点p₂(u_j，v_j)，计算出的p₁(x_i，y_i)在输入图像中的近似对应位置为p₂(u′_j，v′_j)，如果p₂(u_j，v_j)与p₂(u′_j，v′_j)之间的欧氏距离小于某个设定的阈值R_T，将p₁(x_i，y_i)与p₂(u_j，v_j)作为一对匹配特征点；否则，认为p₁(x_i，y_i)在输入图像中没有对应特征点；

(4)采用非极小值抑制算法，保留半径为r的范围内描述子距离最短的点作为用于后续模型拟合的特征点；

(5)应用所有特征对(p_1i，p_2i)(i＝1，2，…，K)，采用最小二乘近似方法，拟合出最小二乘意义下的子图像配准变换模型参数。

Images