在遥感影像的很多具体应用中，要求的配准精度很高[1]Dai X and Khorram S.The effects of image misregistration on the accuracy of remotely sensed change detection.IEEE Trans.Geosci.Remote Sensing，1998，vol.36， pp.1566-1577；[2]Townshend J R G，Justice CO，Gurney C and McManus J.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 and Flusser J.Image registration methods：a survey.Imaging and Vision Computing，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 decomposition of 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 Image Registration for Applications in Remote Sensing.IEEE Transactions on Geoscience and Remote Sensing，2005，43(9)：2127-2137)中首先在参考图 像中提取Harris[8](Harris C and Stephens M J.A combined comer and edge detector.In Alvey Vision Conference，1988，pp.147-152)角点，然后 采用相似和模糊组合不变矩[9](Flusser J and Suk T.Degraded image analysis：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 performance automatic image registration for remote sensing[Ph.D.Dissertation]. George Mason Univertsity，USA，1999)中首先采用某些度量(梯度、熵和 可配准性)来选取纹理丰富的区域，然后只对这些区域采用和文献[10] 类似的操作。上述这些方法强烈依赖于相似性度量标准，可能会面临两 方面的问题。第一，我们知道，一般的相似性度量是建立在统计最优化 基础上的，统计数据(对应于子图像大小)越多，可靠性相对越高。但 在遥感图像中，考虑到局部干扰和变形，子区域不宜选得太大，而区域 选择太小，可靠性将会降低；第二，相似性度量阈值的选取也很困难。 度量标准定得太低，可靠性和精确性都会降低，而度量标准定得太高， 可靠性和精度相对提高了，但检测到的特征数量却大大减少了。
还有一类方法是基于整体特征匹配的。前面两类方法都是对每个特 征逐个进行匹配，而这类方法是对变换模型空间直接进行搜索，来确定 所有特征的整体最优匹配。如[12](Netanyahu N S，Moigne JL and Masek J G.Georegistration of Landsat Data via Robust Matching of Multiresolution Features.IEEE Trans.Geoscience and Remote Sensing，2004，42(7)： 1586-1600)在参考图像和输入图像中分别提取小波极值点，采取分支有 界(Branch-and-bound)搜索策略，采用Hausdorff距离[13](Huttenlocher D P.Klanderman G A and Rucklidge W J.Comparing images using the Hausdorff distance.IEEE Transactions on Pattern Analysis and Machine Intellinence，1993，15：850-863)作为点集的相似性度量，来估计最优的刚 体变换模型，迭代松弛法[14](Ranade S and Rosenfeld A.Point pattern matching by relaxation.Pattern Recognition，1980，12：269-275)也是这类 方法的代表。采用这类方法的前提是，要求参考图像和输入图像的特征 集大部分是匹配的，而且参考图像和输入图像的特征集满足某个整体映 射模型。对不同时相的遥感图像来说，通常场景变化比较大，干扰比较 严重，因此从参考图像和输入图像中提取的特征集误差很大，因此很难 应用这类方法。而且这类方法的另一个致命缺点是计算代价非常大。
造成控制点定位误差的原因是多种多样的，除了计算过程造成的各 种误差外，对多时相影像配准，还存在场景变化问题。这样，通过某些 相似性度量方法(甚至手工方法)确定的“高精度”对应点，由于场景的 自然变化，实际上其中某些对应点并不能反映场景的真实几何变形，而 是由场景的实际变化造成的。毫无疑问，这类点不应该看作是高精度对 应点。但这种情况目前还没有好的办法来区分。
通过以上分析可以看出，尽管文献中有很多种控制点自动选取方法， 但是，这些方法由于精度和鲁棒性等原因，很难真正应用到实际问题中， 目前，遥感影像的控制点选择仍然以手工方式为主。

为了实现大面幅遥感影像的高精度配准，通常需要大量的高精度控 制点，然而，传统的控制点选择方法，其定位精度很难达到子像素级， 为了解决传统技术的问题，本发明的目的是使控制点的定位精度达到子 像素级，为此，本发明提出了一种拟合重采样的遥感影像高精度控制点 自动选择方法。
本发明遥感影像高精度控制点自动选择方法，包括以下两个主要步 骤：
子图像高精度拟合配准；
从已配准子图像中选择控制点。
根据本发明的实施例，所述的子图像高精度配准，包括以下步骤： 对子图像粗配准；利用子图像粗配准信息，对子图像特征点检测；利用 检测得到的子图像特征点进行匹配；利用匹配结果进行统计拟合变换模 型。
根据本发明的实施例，所述的子图像粗配准，采用基于傅立叶Fourier 变换的方法估计影像的近似旋转变换参数和近似平移变换参数。
根据本发明的实施例，所述的子图像特征点检测，是首先检测图像 边缘点，然后从图像边缘点中根据Harris角点度量函数选择特征点。
根据本发明的实施例，所述的子图像特征点检测，采用高斯差分 DOG(Difference of Gaussian，高斯差分)算子的极值点作为特征点。
根据本发明的实施例，所述的子图像特征点匹配，是采用基于特征 相关的匹配方法。
根据本发明的实施例，所述的子图像特征点匹配，是采用基于 SIFT(Scale Invariant Feature Transform)描述子的匹配方法。
根据本发明的实施例，所述基于特征相关的匹配方法，采用灰度相 关系数作为相似性度量标准，在粗配准参数所确定的范围内搜索对应特 征点，并且要求灰度相关系数大于某一预先设定的阈值。
根据本发明的实施例，所述采用基于SIFT描述子的匹配方法，将描 述子欧氏距离的最近邻确定为对应特征点，而且要求描述子最近邻对应 的特征点位于粗配准参数所确定的一定范围内。
根据本发明的实施例，所述统计拟合变换模型，是对已经确定对应 特征点，采用最小二乘方法来拟合子图像变换模型。
根据本发明的实施例，所述从已配准子图像中选择控制点，是采用 变换模型，从子图像中根据实际应用需要选点作为控制点。
根据本发明的实施例，所述基于特征相关的匹配方法步骤是：
(1)在参考图像f1中采用Canny算子检测图像边缘；
(2)对每个边缘点计算Harris角点强度S，如果S大于某个给定的阈 值HT，将这个点作为一个特征点；
(3)对所有特征点，采用局部非极大值抑制算法，保留半径为r的范 围内角点强度最大的点作为后续匹配的侯选特征点。
(4)采用粗配准方法，计算出输入图像f2相对于参考图像f1的近似旋 转角度θ0和平移量(x0，y0)，并将输入图像反方向旋转θ0后得到归整后的 输入图像f2′，在参考图像和归整后的输入图像中仅存在近似的平移偏差 (x0，y0)；
(5)采用模板匹配方法搜索f1中每个侯选特征点在f2′中的匹配点；对 f1中的某个侯选特征点(xi，yi)，只需要在f2′中(xi+x0，yi+y0)周围很小的邻 域R进行；在x和y方向仅取R＝±3像素作为搜索范围。然后，通过对相 关函数的表面拟合，将特征点的匹配位置插值到子像素精度，在f1和f2′中 得到K对对应特征点(p1i，p2i′)(i＝1，2，…，K)；
(6)根据第(4)步计算出来的旋转角度θ0，将归整后输入图像中特征点 p2i的坐标转换为输入图像坐标p2i，这样组成了f1和f2中特征对 (p1i，p2i)(i＝1，2，…，K)；
(7)应用所有特征对(p1i，p2i)(i＝1，2，…，K)，采用最小二乘近似方法，拟 合出最小二乘意义下的子图像配准变换模型参数。
根据本发明的实施例，所述基于SIFT描述子的匹配方法步骤是：
(1)分别在参考图像和输入图像中检测DOG极值点作为特征点，并 计算每个特征点的旋转和尺度不变描述子SIFT；
(2)采用基于Fourier变换的粗配准方法，计算出输入图像f2相对于 参考图像f1的近似旋转角度θ0和平移量(x0，y0)，对参考图像中的任意一 点p1(x，y)，在输入图像中近似对应于p2(u，v)，其中：
u＝xcosθ0+ysinθ0+x0
v＝-xsinθ0+ycosθ0+y0
(3)对参考图像中的每个特征点p1(xi，yi)，计算它的SIFT描述子与输 入图像中所有特征描述子的欧氏距离，假设最近邻距离对应于输入图像 中的特征点p2(uj，vj)，而根据上式计算出来的p1(xi，yi)在输入图像中的近 似对应位置为p2(uj′，vj′)，如果p2(uj，vj)与p2(uj′，vj′)之间的欧氏距离小于某个 设定的阈值RT，那么，将p1(xi，yi)与p2(uj，vj)作为一对匹配特征点；否则， 认为p1(xi,yi)在输入图像中没有对应特征点；
(4)采用非极小值抑制算法，保留半径为r的范围内描述子距离最短 的点作为用于后续模型拟合的特征点；
(5)应用所有特征对(p1i，p2i)(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 phase correlation image alignment method.In Proceedings of the IEEE 1975 International 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变换的平移和旋转图像配准方法。 这种方法分两个步骤来完成，第一步先计算旋转角度，然后再确定平移 量。
考虑被配准的两幅图像f1(x，y)和f2(x，y)，其中f2(x，y)是f1(x，y)经过 平移(x0，y0)和旋转θ0后的图像，即
f2(x，y)＝f1(xcosθ0+ysinθ0-x0，-xsinθ0+cosθ0-y0)
根据Fourier变换的平移和旋转特性，f1(x，y)和f2(x，y)的Fourier变换之 间的关系为：
F 2 ( ξ , η ) = e - j 2 π ( ξx 0 + ηy 0 ) F 1 ( ξ cos θ 0 + η sin θ 0 , - ξ sin θ 0 + η cos θ 0 )
若M1和M2代表F1和F2的幅度谱，那么，从上式可以得到：
M2(ξ，η)＝M1(ξcosθ0+ηsinθ0，-ξsinθ0+ηcosθ0)
如果从F1和F2的幅度谱考虑，从上式很容易看出，两个频谱幅度是 相同的，只是M2相对于M1旋转了θ0。将幅度谱转化到极坐标空间后， 旋转位移可以简化为类似于平移的表达形式：
M1(ρ，θ)＝M2(ρ，θ-θ0)
这样就可以采用经典的相位相关法计算出角度θ0。
根据计算出来的角度θ0，将f2(x，y)旋转-θ0得到归整后的图像 f2′(x，y)，这时，在f1(x，y)和f2′(x，y)之间仅有平移误差，可在图像空间再 一次应用相位相关法求得平移量。这样，经过两次相位相关运算，就求 出了旋转角度θ0和平移量(x0, y0)。
大量实验表明，基于Fourier变换的配准方法，其角度计算误差与分 度间隔有关，一般不会超过1～2个分度间隔，平移误差一般不超过1～ 2个像素范围。
从后面的讨论我们将会看到，首先计算粗配准参数，对后续精确匹 配过程具有重要的意义。下面介绍实验中用到的两种配准方法。
2基于特征相关的匹配方法
所述的子图像特征点匹配，是采用基于特征相关的匹配方法。
基于特征相关的图像配准方法，其基本思想是：只在参考图像中提 取特征，而在输入图像中一定范围内采用类似于模板匹配的方法，根据 某种相似性度量标准，搜索参考图像中的特征在输入图像中的对应特征。 由于我们的实验对象是不同时相的同波段影像，因此我们采用了基于灰 度相关的相似性度量。
所述基于特征相关的匹配方法，采用灰度相关系数作为相似性度量 标准，在粗配准参数所确定的范围内搜索对应特征点，并且要求灰度相 关系数大于某一预先设定的阈值。
在基于灰度相关的相似性度量标准中，相关系数被证明是比较鲁棒 的，而且相关系数对灰度的线性变化和加性高斯噪声不太敏感。
灰度相关系数定义为：
C coef = Σ i = 1 M Σ j = 1 N ( f 1 ( x i , y j ) - μ 1 ) ( f 2 ( x i , y j ) - μ 2 ) Σ i = 1 M Σ j = 1 N ( f 1 ( x i , y j ) - μ 1 ) 2 Σ i = 1 M Σ j = 1 N ( f 2 ( x i , y j ) - μ 2 ) 2
首先在参考图像中提取特征点集，对每个特征点Pi，以Pi为中心， 裁剪大小为M×N像素的一个窗口，称为模板图像f1，窗口中的像素的灰 度均值为μ1，然后在输入图像中(一定搜索范围R内)搜索模板图像的 最佳匹配位置。在每个搜索位置，裁剪和模板图像一样大小的浮动图像 f2，灰度均值为μ2，计算f1和f2的灰度相关系数Ccoef。当所有位置搜索 完毕，如果最大相关系数大于某一设定的阈值Tc，那么将最大相关系数 对应的浮动图像中心Pj作为Pi的匹配位置。
为了提高匹配的可靠性，一般希望选择的模板图像具有比较丰富的 纹理，而不是灰度变化比较平坦的区域。基于这样的考虑，我们采用了 如下特征点选择算法：
所述的子图像特征点检测，是首先检测图像边缘点，然后从图像边 缘点中根据Harris角点度量函数选择特征点。
一般认为，图像边缘包含了丰富的信息，对应于梯度变化较大的区 域。因此，我们的选择策略是从图像边缘点中提取特征点。这里采用了 Canny边缘检测算法[17](Canny J.A computational approach to edge detection，IEEE Transactions on Pattern Analysis and Machine Intelligence， 1986(8)：679-698)。
当边缘检测完成后，取边缘点中具有较大“强度”的那些点作为特征 点。通常，点的强度是根据像素灰度值沿两个方向的变化来测量的。 Schmid[18](Schmid C，Mohr R and Bauckhage C.Evaluation of interest point detectors.International Journal of Computer Vision，2000，37(2)： 151-172)曾对不同的特征点检测算法进行了评价和比较，得出的结论是， Harris算法[8]性能最优，而且信息熵最高(意味着纹理信息比较丰富， 后续的匹配可靠性较高)，因此，本发明采用了Harris算法提取特征点。
概括起来，基于特征相关的子图像配准方法可归纳为以下步骤。
(1)在参考图像f1中采用Canny算子检测图像边缘。
(2)对每个边缘点计算Harris角点强度S，如果S大于某个给定的阈 值HT，那么将这个点作为一个特征点。
(3)对所有特征点，采用局部非极大值抑制算法，仅保留半径为r的 范围内角点强度最大的点作为后续匹配的侯选特征点。
(4)首先采用粗配准方法，计算出输入图像f2相对于参考图像f1的近 似旋转角度θ0和平移量(x0，y0)，并将输入图像反方向旋转θ0后得到归整 后的输入图像f2′，这时，在参考图像和归整后的输入图像中仅存在近似 的平移偏差(x0，y0)。
(5)采用模板匹配方法搜索f1中每个侯选特征点在f2′中的匹配点。对 f1中的某个侯选特征点(xi，yi)，只需要在f2′中(xi+x0，yi+y0)周围很小的邻 域R进行。我们在x和y方向仅取R＝±3像素作为搜索范围。然后，通过 对相关函数的表面拟合，将特征点的匹配位置插值到子像素精度。最后， 在f1和f2′中得到K对对应特征点(p1i，p2i)(i＝1，2，…，K)。
(6)根据第(4)步计算出来的旋转角度θ0，将归整后输入图像中特征点 p2i′的坐标转换为输入图像坐标p2i，这样组成了f1和f2中特征对 (p1i，p2i)(i＝1，2，…，K)。
(7)所述统计拟合变换模型，是对已经确定对应特征点，采用最小二 乘方法来拟合子图像变换模型。
应用所有特征对(p1i，p2i)(i＝1，2，…，K)，采用最小二乘近似方法，拟合 出最小二乘意义下的子图像配准变换模型参数(本发明对应仿射变换模 型)。
3基于SIFT描述子的特征匹配方法
所述的子图像特征点匹配，是采用基于描述子SIFT的匹配方法。
这种方法沿用了基于特征的图像配准一般框架：首先在参考图像和 输入图像中分别独立提取特征点，然后采用基于不变量描述子的相似性 度量方法进行特征点匹配。
SIFT描述子[19](Lowe D G.Distinctive image features from scale-invariant keypoints.International Journal of Computer Vision，2004， 60(2)，91-110)，是图像匹配和目标识别等应用中性能最为杰出的描述 子之一[20](Mikolajczyk K and Schmid C.A performance evaluation of local descriptors.IEEE Trans.Pattern Analysis and Machine Intelligence， 2005，27(10)：1615-1630)，SIFT描述子在干扰和变形比较严重的遥感图 像配准中同样具有很强的鲁棒性和很强的特征辨析性能。
为了和前面讲的基于特征相关的子图像配准方法进行比较，我们对 基于SIFT特征的配准方法也进行了实验。
基于SIFT描述子的子图像配准方法可以归纳为以下步骤：
(1)所述的子图像特征点检测，采用高斯差分DOG算子的极值点作 为特征点。
采用文献[19]中的方法，分别在参考图像和输入图像中检测DOG (Difference of Gaussian)极值点作为特征点，并计算每个特征点的旋转和 尺度不变SIFT描述子。
(2)采用基于Fourier变换的粗配准方法，计算出输入图像f2相对于 参考图像f1的近似旋转角度θ0和平移量(x0，y0)，那么，对参考图像中的 任意一点p1(x，y)，在输入图像中近似对应于p2(u，v)，其中：
u＝xcosθ0+ysinθ0+x0
v＝-xsinθ0+ycosθ0+y0
(3)所述采用基于描述子SIFT的匹配方法，将描述子欧氏距离的最 近邻确定为对应特征点，而且要求描述子最近邻对应的特征点位于粗配 准参数所确定的一定范围内。
对参考图像中的每个特征点p1(xi，yi)，计算它的SIFT描述子与输入 图像中所有特征描述子的欧氏距离，假设最近邻距离对应于输入图像中 的特征点p2(uj，vj)，而根据上式计算出来的p1(xi，yi)在输入图像中的近似 对应位置为p2(uj′，vj′)，如果p2(uj，vj)与p2(uj′，vj′)之间的欧氏距离小于某个设 定的阈值RT，那么，将p1(xi，yi)与p2(uj，vj)作为一对匹配特征点。否则，认 为p1(xi，yi)在输入图像中没有对应特征点。
需要注意的是，这里我们没有对描述子的最近邻距离设定阈值，一 方面是因为SIFT描述子是一个128维的高维向量，最近邻距离会相差 很大，很难设定某个固定阈值；另外一个原因是，鉴于SIFT描述子的 高辨析能力，如果在众多(通常在400×400大小的图像中可以检测到数 百到数千DOG极值点)待比较的特征点中，其最近邻距离对应的特征 点位于预先设定的比较小的范围内的话，那么这也间接地说明，正确匹 配的可能性是很大的。
(4)采用类似于非极大值(这里是非极小值)抑制算法，仅保留半径 为r的范围内描述子距离最短的点作为用于后续模型拟合的特征点。
(5)所述统计拟合变换模型，是对已经确定对应特征点，采用最小二 乘方法来拟合子图像变换模型。
应用所有特征对(p1i，p2i)(i＝1，2，…，K)，采用最小二乘近似方法，拟合 出最小二乘意义下的子图像配准变换模型参数。
从上面介绍的两种子图像配准方法来看，子图像的粗配准具有重要 的意义，这体现在：第一，对基于特征相关的方法来说，减小了相关搜 索范围，提高了计算效率；第二，对基于描述子的方法来说，提高了匹 配可靠性；第三，也是最重要的，无论相关方法还是描述子方法，这样 确定的特征对应点的误差被控制在一定的范围内，即对应点最大误差不 会超过预先设定的搜索范围，这在一定程度上保证了后续模型拟合的精 确性和可靠性。
所述从已配准子图像中选择控制点，是采用变换模型，从子图像中 根据实际应用需要选点作为控制点。
当采用基于特征相关的方法或者基于SIFT描述子的匹配方法拟合 出子图像变换模型后，可根据应用需要，利用变换模型，从子图像中选 择适当的对应点作为大面幅遥感影像的控制点。
4实验结果及讨论
在介绍实验内容之前，我们先讨论一下图像配准的精度评价问题。
4.1图像配准的精度评价概述
对于图像配准问题来说，提供一个最终的配准误差估计是很有必要 的，否则在实际应用中很难被用户接受。但通常情况下由于很难得到真 实数据，因此，这个问题并不是一个很简单的问题。这里我们介绍几个 与我们后续实验有关的误差评价方法，尽管这些方法无法可靠地反映最 终的配准误差，至少也可以提供一些分析的依据。
配准误差常采用以下几种方法来评价。
●最简单的度量方法是计算控制点(CPs)的均方误差(CPE)。设参 考图像中的控制点集为P1(xi，yi)，对应的输入图像中的控制点集为 P2(ui，vi)，i＝1，2，…，K，K表示控制点的数量。假设通过所有控制点 计算出来的映射模型为u＝T1(x，y)，v＝T2(x，y)，那么，CPE定义为：
CPE = Σ i = 1 K ( 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更有意 义。然而，确切地说，控制点的误差与图像的配准误差是不同的 两个概念，很多情况下不具有可比性。
●三角验证：对同一场景至少有三幅图像的情况，可以采用一致性 检查方法。对三幅图像的情况，首先根据控制点计算两两图像之 间的映射函数Tij，假设f2＝T12(f1)，f3＝T23(f2)，那么，根据模型传递 关系，f3＝T23(T12(f1))，另外，f3′＝T13(f1)。从理论上来说，f3＝f3′， 这样，通过度量f3与f3′的误差来验证映射模型的正确性。f3与f3′的 误差可通过选择一些测试点进行度量。例如，在我们的实验中， 在图像f1中沿两个方向每隔10个像素选择一个测试点，最后来测 试这些点在f3与f3′的误差，通过计算均方误差来度量这些测试点 的整体误差。
●相互验证：采用多种方法进行一致性检查。用不同特点的方法分 别配准，如果这两种方法的差别很小，那么可认为该方法具有较好的配 准精度(尽管不能严格保证)。我们对前面提到的两种配准方法进行了 相互验证。但是，如果两种方法差别比较大，那么无法确定到底哪种方 法不准确或两种方法都不准确。
4.2参数设置及实验结果
(1)参数设置
①Canny边缘检测算法采用Matlab中图像处理工具箱中标准函数的 缺省设置。
②Harris角点检测算法中，高斯函数标准差为1，角点度量的阈值HT 设为1000；
③特征点局部非极大值抑制半径设为5象素；
④模板图像大小设置为50×50像素；
⑤模板匹配的搜索范围R设置为±3像素；
⑥灰度相关系数阈值Tc设为0.7；
⑦基于SIFT描述子方法中，描述子最近邻点的位置误差阈值RT设 为5像素。
(2)模型参数
实验中，子图像映射模型采用仿射变换。第一幅图像为参考图像， 像素点坐标表示为(x，y)，第二幅图像为输入图像，像素点坐标表示为(u，v)， 模型变换关系表示为：
u＝a11x+a12y+tx
v＝a21x+a22y+ty
后面统计表中，“模型参数”中每一项代表上述对应变换的模型系数。
(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”表示测 试点的最大均方根误差，Ratio1表示均方根误差大于1个像素的测试点 占所有测试点的比例。
“TPE”是采用随机选取的一半点计算变换模型，而用另一半点进行 测试的均方根误差，“TPE_Max”表示测试点的最大均方根误差，Ratio2 表示均方根误差大于1个像素的测试点占所有测试点的比例。
表中“两种方法相互验证”是指，在输入图像(第一幅图像)中沿x和 y方向每隔10个像素选择一个测试点(对400×400的图像，共选取 39×39＝1521个测试点)，构成一个测试点集P，假设两种方法求出的变 换模型分别为T1和T2，那么表中“RMS2”对应于点集T1(P)和T2(P)的均方根 误差。
表中“三角验证”是指，从第一幅图像中采取和以上相同的方式均匀 选取1521个点作为测试点集P，假设第一幅图像到第二幅图像、第二幅 图像到第三幅图像，第一幅图像到第三幅图像的变换模型分别为T12、T23 和T13，那么表中“RMS3”对应于点集T23(T12(P))和点集T13(P)的均方根误差。
表4中，“和理论数据比较”是指，已知理论变换模型为T12，实际计 算出来的模型为T12′，采用以上相同的测试点集P，“RMS4”对应于点集 T12(P)和点集T12′(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-1 Vs.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)， 图中：
最上面一行为同一场景不同时间的三幅Landsart3波段影像，分别 于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中，每一对图像已由专业人员采用手工选取控制点的方法进行 了配准，也即参考图像和输入图像的坐标理论上是完全对应的(不考虑 配准误差时)。已经知道了两幅图像之间的“真实”变换模型， (a11，a12，tx，a21，a22，ty)＝(1，0，0，0，1，0)，我们试图来验证该方法的正确性。
如图4所示：是利用本发明对Landsat影像进行配准的结果(场景 4)，图中：
上面一行为同一场景不同时间的两幅Landsart 3波段影像，分别于 2001年5月和2001年7月获得，并且已经经过手工配准，大小均为 400×400像素。下面左右图分别为基于特征相关和基于SIFT的配准方法 的对应特征点显示。括号中数字，是采用非极大值抑制后的特征点数量， 这里的抑制半径为10个像素，和后面表中数据不同，表中的抑制半径为 5像素。
如图5所示：是利用本发明对Landsat影像进行配准的结果(场景5)，
图中：
上面一行为同一场景不同时间的两幅Landsart3波段影像，分别于 2001年5月和2001年7月获得，并且已经经过手工配准，大小均为 400×400像素。下面左右图分别为基于特征相关和基于SIFT的配准方法 的对应特征点显示。括号中数字，是采用非极大值抑制后的特征点数量， 这里的抑制半径为10个像素，和后面表中数据不同，表中的抑制半径为 5像素。
表4中，“和理论模型相比较”是指，将表中方法计算出来的模型与 理论模型采用和前面相同的测试点集来计算误差。从“两种方法相互比 较”和“和理论模型相比较”可看出，最终配准确实达到了子像素级精度。 但同时也发现，在“和理论模型相比较”中，误差有些偏大，这个结果似 乎稍微有点奇怪。和前面三个场景相比较，这些图像之间相隔时间较短 (只有大约两个月)，噪声也相对较小，另外，特征点数量很多，分布也 很好，从理论上说，计算误差应该比场景1，2，3小一些，但实际情况 不是如此，我们推测，是“真实”模型不精确所致。或许，我们可以认为， 我们提出的自动方法比手工方法更精确。
综合以上分析并结合图表数据可以看到，特征点大都具有一定的定 位误差，甚至误差还比较大(几个象素)，而且大于1个象素误差的特征 点还占有相当的比例。我们不知道造成这些误差的具体原因，或许是计 算过程的各种因素造成的，或许是场景真实变化造成的，但值得庆幸的 是，最终的配准精度都达到了亚像素级精度。另外，我们提出的先粗配 准，然后采用小范围搜索的策略来确定特征点对应，有以下几个主要优 点：第一，有效地增加了特征点数量。比如，当我们限定搜索范围后， 对于基于特征相关的配准方法，可以适当降低相关阈值，这样既可以增 加特征点数量，而且特征点的定位误差可以限定在一定的范围内，而不 影响后续模型的估计；第二，对相关运算来说，减小搜索范围，意味着 运算量的减少；第三，鲁棒性强，可靠性高。从场景图像看，各种噪声 还是比较严重的，但无论是粗配准还是后续的精配准，都具有很强的适 应性。这些特点很适合于复杂的遥感图像。
(4)当然，当各种噪声大到一定程度后，这些方法得出的配准结果将 无法满足子像素精度要求，场景6就是一个例子。
如图6所示：是利用本发明对Landsat影像进行配准的结果(场景 6)，图中：
最上面一行为同一场景不同时间的三幅Landsart3波段影像，分别 于1996年，1997年和2002年获得，大小均为400×400像素。中间和最 下面一行分别为基于特征相关和基于SIFT的配准方法的对应特征点显 示。图中“6-1Vs.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)

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