CN103473765A - 一种多光谱遥感图像的配准与校正方法、系统 - Google Patents

Abstract

本发明属于图像处理技术领域，提供了一种多光谱遥感图像的配准与校正方法、系统。该方法及系统是选取经初步几何校正的多光谱遥感图像中的一个作为基准图像，其余作为待配准图像，首先在基准图像中提取最佳图像子块，之后将最佳图像子块与每一待配准图像进行双尺度空间的亚像素级匹配以得到匹配点，之后根据匹配点计算投影变换矩阵，并将相应的光谱图像投影到基准图像的坐标空间，之后采用必要的图像插值后，得到在基准图像的坐标空间配准和校正后的多光谱图像。相对于现有的多光谱遥感图像配准方法，计算效率高，配准准确度高，特别适合于大幅面的多光谱遥感图像的配准及校正。

Description

一种多光谱遥感图像的配准与校正方法、系统

技术领域

本发明属于图像处理技术领域，尤其涉及一种多光谱遥感图像的配准与校正方法、系统。 

背景技术

在多光谱遥感图像成像中，经过初步几何校正的、具有不同光谱的大幅面遥感图像之间一般都存在有一定的像差，像差类型有：微小的空间平移误差、微小的旋转误差、微小的尺度缩放误差，以及由该三种误差共同作用得到的误差。在具体图像中一般表现为，在图像的不同区域，各光谱图像之间的垂直和水平方向的失配误差一般是不同的。经常存在的情况是，两个不同光谱的图像在局部观察表现为平移和旋转误差，即通过平移和旋转可以在局部基本消除像差，但在整幅图像上存在尺度缩放误差，表现在两个不同光谱图像的左上角区域可能存在1-2个像素的像差，随着向图像右下角区域观察，像差逐渐增加，像差可能会达到5-8个，或10个像素以上，在某些波段之间的像差甚至可能会达到40个像素以上，且在水平和垂直方向上像差的变化一般是不同的。并且，像差的大小和变化情况在各个不同光谱的图像之间一般也各不相同。 

现有技术提出的针对多光谱遥感图像进行配准的方法主要有以下几种： 

一、利用人工手动选取特征点对的方式实现对多光谱遥感图像的配准。该方法十分耗时，计算效率低，不能满足大数据量的应用需求，且受人工主观因素的影响，配准准确度低。 

二、利用基于灰度的自动配准算法实现对多光谱遥感图像的配准。该类方法的原理是：直接利用灰度信息建立各图像之间的相似性度量，然后利用优化搜索策略寻找配准的最佳变换参数。相对于前述方法，该类方法具有较高的配 准精度，但由于图像中的所有像素点均参与运算，计算量大、计算效率低，不适合对大幅面多光谱遥感图像的配准。 

三、利用基于特征的自动配准算法实现对多光谱遥感图像的配准。该类方法的原理是：首先提取图像中的角点、边缘、直线等显著特征，然后确定两幅图像特征集合中各个特征点的对应关系，从而实现配准。相对于前述方法，该类方法计算量小、速度快，得到了广泛应用。举例来说，可有如下两种典型的实现方式： 

1、首先基于特征不变特性变换（Scale Invariant Feature Transform，SIFT）算法提取特征点进行粗匹配，然后用Harris算子或Canny算子提取特征点建立不规则三角网进行精匹配。由于该方式需要计算大量的特征点，因此仍旧存在计算效率低的问题，且无法保证特征点对的匹配准确性和分布均匀性，对大幅面遥感图像的算法计算量大，配准精度得不到保证。 

2、在大幅面遥感图像上先建立二级规则网格，用优化及均匀性准则选取特征网格，然后在选取的特征网格中并行采用SIFT算法提取特征点，利用相关原理和特征点之间的空间距离约束建立匹配的特征点对，并应用最小二乘方法计算图像变换参数实现图像的配准。该方式由于采用了网格技术，改善了特征点对分布的均匀性，但该方法无法保证选取的特征网格整体分布的均匀性，因此依然会存在特征点对分布的整体不均匀性，从而影响其配准精度。另外该方式需要根据图像人为设定SIFT参数的范围，且依然需要计算大量的特征点，对大幅面遥感图像的算法效率依然较低。 

同时，对于前述两种典型的实现方式，由于基于SIFT算法进行图像配准适宜于具有同一模式传感器的成像图像，允许待配准的各图像之间可以存在尺度、平移、缩放的差异。例如在相同模式可见光相机之间成像的图像，或者在相同模式红外相机之间成像的图像，特征的相似性较强。但对不同模式的可见光成像与红外成像的图像之间，特征的相似性较弱。因此在多模图像之间，基于SIFT变换、Harris角点或Canny边缘的匹配方法具有较大的局限性。而在多 光谱遥感图像中，各个不同成像波段所使用的传感器不同，属多模成像模式。因此在不同波段成像的多光谱图像之间使用SIFT变换、Harris角点或Canny边缘进行位置匹配并不很合理。 

综上所述，现有技术提出的针对多光谱遥感图像进行配准的方法存在计算效率低、配准准确度低的问题。 

发明内容

本发明实施例的目的在于提供一种多光谱遥感图像的配准与校正方法，旨在解决现有针对多光谱遥感图像配准方法或者采用人工手动选取特征点对的方式实现、或者利用基于灰度的自动配准算法实现、或者利用基于特征的自动配准算法实现，三种方式均存在计算效率低和/或配准准确度低的问题。 

本发明实施例是这样实现的，一种多光谱遥感图像的配准与校正方法，所述方法包括以下步骤： 

选择多光谱遥感图像中的一个光谱图像作为基准图像，而将所述多光谱遥感图像中的其它光谱图像作为待配准图像，所述多光谱遥感图像包括多个经初步几何校正的光谱图像； 

在所述基准图像中提取具有整体最佳匹配特征的四个最佳图像子块； 

将所述四个最佳图像子块分别与每一待配准图像进行双尺度空间的亚像素级匹配，得到与每一待配准图像对应的四对匹配点； 

根据所述四对匹配点，计算相应的待配准图像到所述基准图像的投影变换矩阵，并将相应的待配准图像投影到所述基准图像的坐标空间，之后采用必要的图像插值后，得到在所述基准图像的坐标空间配准和校正后的多光谱图像。 

本发明实施例的另一目的在于提供一种多光谱遥感图像的配准与校正系统，所述系统包括： 

选择模块，用于选择多光谱遥感图像中的一个光谱图像作为基准图像，而将所述多光谱遥感图像中的其它光谱图像作为待配准图像，所述多光谱遥感图 像包括多个经初步几何校正的光谱图像； 

提取模块，用于在所述基准图像中提取具有整体最佳匹配特征的四个最佳图像子块； 

匹配模块，用于将所述四个最佳图像子块分别与每一待配准图像进行双尺度空间的亚像素级匹配，得到与每一待配准图像对应的四对匹配点； 

配准及校正模块，用于根据所述四对匹配点，计算相应的待配准图像到所述基准图像的投影变换矩阵，并将相应的待配准图像投影到所述基准图像的坐标空间，之后采用必要的图像插值后，得到在所述基准图像的坐标空间配准和校正后的多光谱图像。 

本发明提出的多光谱遥感图像的配准与校正方法是选取经初步几何校正的多光谱遥感图像中的一个作为基准图像，其余作为待配准图像，首先在基准图像中提取最佳图像子块，之后将最佳图像子块与每一待配准图像进行双尺度空间的亚像素级匹配以得到匹配点，之后根据匹配点计算投影变换矩阵，并将相应的光谱图像投影到基准图像的坐标空间，之后采用必要的图像插值后，得到在基准图像的坐标空间配准和校正后的多光谱图像。相对于现有的多光谱遥感图像配准方法，计算效率高，配准准确度高，特别适合于大幅面的多光谱遥感图像的配准及校正。 

附图说明

图1是本发明实施例提供的多光谱遥感图像的配准与校正方法的流程图； 

图2是本发明实施例中，在基准图像中提取四个最佳图像子块的详细流程图； 

图3a是本发明中经初步几何校正的多光谱遥感图像中的一个光谱图像的示例图； 

图3b是本发明中经初步几何校正的多光谱遥感图像中的另一个光谱图像的示例图； 

图3c是图3a和图3b中有效图像区域及其外边界示意图； 

图3d是本发明中有效环形区域、有效环形区域的中心线，以及中心线上形成一四边形的四个顶点的示例图； 

图3e是本发明中无重叠网格化的基准图像、基准图像的有效子图像块、以及中心分别在四个顶点的四个图像分析窗口的示例图； 

图4a是本发明中一灰度图像块的示例图； 

图4b是图4a的灰度直方图； 

图5本发明实施例中，计算滑动过程中的最大整体显著性度量以及最大整体显著性度量所对应的位置的详细流程图； 

图6本发明实施例中，最佳图像子块与待配准图像进行双尺度空间的亚像素级匹配以得到四对匹配点的详细流程图； 

图7是本发明提供的多光谱遥感图像的配准与校正系统的结构图； 

图8是图7中提取模块的结构图； 

图9是图8中第五计算子模块的结构图； 

图10是图7中匹配模块的结构图。 

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。 

为了解决现有技术存在的问题，本发明提出的多光谱遥感图像的配准与校正方法是选取经初步几何校正的多光谱遥感图像中的一个作为基准图像，其余作为待配准图像，首先在基准图像中提取最佳图像子块，之后将最佳图像子块与每一待配准图像进行双尺度空间的亚像素级匹配以得到匹配点，之后根据匹配点计算投影变换矩阵，并将相应的光谱图像投影到基准图像的坐标空间，之后采用必要的图像插值后，得到在基准图像的坐标空间配准和校正后的多光谱 图像。 

图1示出了本发明提供的多光谱遥感图像的配准与校正方法的流程，包括以下步骤： 

S1：选择多光谱遥感图像中的一个光谱图像作为基准图像，而将多光谱遥感图像中的其它光谱图像作为待配准图像。 

其中，多光谱遥感图像包括多个经初步几何校正的光谱图像；基准图像的选择没有特别限制，可选取某一特定波段的光谱图像作为基准图像，也可选取与待配准图像波段相近、且对比度和亮度都较高的光谱图像作为基准图像，以保证多模图像之间匹配的区域一致性较好，图像之间匹配的鲁棒性更强。 

S2：在基准图像中提取具有整体最佳匹配特征的四个最佳图像子块。 

其中，具有整体最佳匹配特性的四个最佳图像子块是指：四个最佳图像子块的中心点中的任意三点不共线，任意三点所构成的三角形的面积都尽可能大，即四个中心点在基准图像中形成的四边形的面积要尽可能大，形状要尽可能接近正方形。这样，后续由其计算出的四对匹配点计算得到的投影变换矩阵在对图像进行变换之后的误差可以控制到尽可能小。另外为了达到好的匹配定位效果，要求每一最佳图像子块作为互信息匹配的局部区域特性要显著。 

进一步地，如图2所示，步骤S2包括以下步骤： 

S21：确定基准图像中的有效图像区域，并在有效图像区域内，沿有效图像区域外边界确定一有效环形区域，并计算有效环形区域的中心线，以及中心线上构成一四边形的四个顶点。 

如图3a和图3b示出了经初步几何校正的多光谱遥感图像中的2个光谱图像的示例。其中间的四边形区域即为有效图像区域，四边形区域外围的黑色区域为无效图像区域，可通过对基准图像进行图像阈值化处理得到有效图像区域，并通过图像形态学边界提取方法得到有效图像区域外边界，如图3c中的阴影部分即为处理后得到的有效图像区域。一般地，有效图像区域是一四边形区域。 

在图3c所示的有效图像区域内，沿有效图像区域外边界确定一个有效环形 区域，当有效图像区域是一个四边形区域时，有效环形区域是一个四边形环形区域，如图3d的阴影部分所示。一般地，可通过确定有效环形区域的环宽度的方法确定有效环形区域，对于大幅面遥感图像，有效环形区域的环宽度一般可确定为512像素或1024像素，也可以根据情况取其它宽度值。 

有效环形区域的中心线是与有效环形区域外边界的距离为环宽度一半的线。在有效环形区域是一个四边形区域时，中心线也是一个四边形，如图3d中的虚线所示，该四边形可由其四个顶点d1，d2，d3，d4决定。 

S22：对基准图像进行网格化规则分块，之后对分块后各子块中的有效子图像块进行标定，获得基准图像的有效子图像块的集合。 

本发明中，对基准图像进行网格化规则分块是指把基准图像切分为大小相同的子块，子块的大小一般可选128×128或64×64，切分后的各子块之间可以有重叠，也可以无重叠。若被切分后的子块中的所有像素点均为步骤S21中有效图像区域中的像素点，则子块即为有效子图像块，否则为无效子图像块。有效子图像块的全体称为有效子图像块的集合。可利用步骤S21中得到的有效图像区域的模板，对有效子图像块进行标定，此处不再赘述。 

S23：分别以中心线上的四个顶点为中心确定基准图像上的四个图像分析窗口D_k(0)，k=1，2，3，4。 

本发明中，分别以步骤S21中得到的有效环形区域的中心线上的四个顶点d1，d2，d3，d4为中心确定四个图像分析窗口D₁(0)、D₂(0)、D₃(0)、D₄(0)。其中，d1是图像分析窗口D₁(0)的中心，d2是图像分析窗口D₂(0)的中心，d3是图像分析窗口D₃(0)的中心，d4是图像分析窗口D₄(0)的中心。图像分析窗口区域D₁(0)、D₂(0)、D₃(0)、D₄(0)的大小可根据有效环形区域的环宽度确定，若环宽度为512，则图像分析窗口的大小可取512×512，若环宽度为1024，则图像分析窗口的大小可取1024×1024。 

如图3e示出了一个无重叠网格化的基准图像、基准图像的有效子图像块、以及中心分别在d1，d2，d3，d4的四个图像分析窗口D₁(0)、D₂(0)、D₃(0)、 D₄(0)的实例。在每一图像分析窗口中，全阴影的子块表示该子块中的每一个像素都在有效图像区域中，因此相应子块为有效图像子块，而部分阴影的子块中有部分像素不在有效图像区域中，因此相应子块是无效图像子块。 

S24：利用基准图像的有效子图像块的集合，计算四个图像分析窗口D_k(0)分别的有效子图像块的集合S_k(0)={S_ki(0)}，计算图像分析窗口D_k(0)的有效子图像块的集合中每一图像块S_ki(0)的显著性度量F_ki(0)，计算每一图像分析窗口D_k(0)的显著性度量F_k(0)，并计算四个图像分析窗口D_k(0)的整体显著性度量F(0)。 

其中，区域的显著性度量是判断一个图像区域是否更适合用于匹配的重要指标。在现有技术中一般是采用图像直方图的信息熵作为度量指标。但在使用互信息方法进行图像匹配时，直方图信息熵度量值选出的图像块往往比较复杂，匹配效果受到影响。为此，本发明中，定义一种新的、用于互信息匹配的图像块Ｓ的显著性度量F的计算公式如下： 

$F={|m4-m3|}^2\·\frac{{|m2-m1|}^2}{1+{\mathrm{\σ}}_1^2+{\mathrm{\σ}}_2^2}---\left(1\right)$

公式（1）中，假设图像块为Ｓ，m0是图像块Ｓ的所有像素灰度值的均值；m1是图像块Ｓ中灰度值小于m0的所有像素灰度值的均值，m2是图像块Ｓ中灰度值大于m0的所有像素灰度值的均值；σ₁是图像块Ｓ中灰度值小于m0的所有像素灰度值的标准差，σ₂是图像块Ｓ中灰度值大于m0的所有像素灰度值的标准差；m3是图像块Ｓ中灰度值小于m1的所有像素灰度值的均值，m4是图像块Ｓ中灰度值大于m2的所有像素灰度值的均值。图4a为一灰度图像块S的实例，图4b是图4a的灰度直方图、以及计算出的m0,m1,m2,m3,m4的值，可见，一般满足：m3<m1<m0<m2<m4。在公式（1）的显著性度量F中，项|m4-m3|²可以表征图像整体的灰度动态范围和对比度的特性；项

是一种类fisher可分性特征，表征了在将图像的像素按照整体均值划分为两类时这两类的可分 性特性。 

定义图像分析窗口D_k(0)，k=1，2，3，4中的有效子图像块集合S_k(0)为S_k(0)={S_ki(0)}，其中，S_ki(0)是在步骤S22中进行网格化规则分块后标定的有效子图像块，并且S_ki(0)的几何中心点处在图像分析窗口D_k(0)中，S_k(0)是满足上述条件的图像集合。又设集合S_k(0)中共有N_k个图像块，对集合S_k(0)中的每个图像块S_ki(0)，使用公式（1）计算其显著性度量值，记为F_ki(0)，则图像分析窗口D_k(0)的显著性度量F_k(0)的计算公式如下所示： 

F_k(0)=max{F_ki(0)，i=1，2，...，N_k}        （2） 

即图像分析窗口D_k(0)的显著性度量值F_k(0)是相应图像分析窗口中的所有有效子图像块的显著性度量的最大值。四个图像分析窗口D_k(0)的整体显著性度量F(0)的计算公式如下所示： 

F(0)=min{F₁(0)，F₂(0)，F₃(0)，F₄(0)}      （3） 

即四个图像分析窗口D_k(0)的整体显著性度量F(0)是该四个图像分析窗口D_k(0)的显著性度量的最小值。 

S25：将四个图像分析窗口D_k(0)的中心分别沿中心线、按顺时针或逆时针方向以等滑动步长逐步滑动到相邻的下一顶点，并计算滑动过程中的最大整体显著性度量F(θ*)，以及最大整体显著性度量所对应的位置θ*。 

进一步地，如图5所示，步骤S25包括以下步骤： 

S251：将中心线被四个顶点分割后的四条中间线分别均分为若干段，以得到相应的滑动步长。 

S252：将四个图像分析窗口的中心分别沿中心线、按顺时针或逆时针方向以相应的滑动步长向相邻的下一顶点滑动。 

S253：图像分析窗口每滑动固定步数，计算相应的整体显著性度量，进而得到相应的图像分析窗口滑动到下一顶点过程中的最大整体显著性度量，以及最大整体显著性度量所对应的位置。 

举例来说，参照图3e所示，若将四个图像分析窗口D_k(0)的中心沿有效环 形区域的中心线，按顺时针方向以等滑动步长逐步移动到相邻的下一顶点，即是指：图像分析窗口D₁(0)的中心为d1，图像分析窗口D₂(0)的中心为d2，图像分析窗口D₃(0)的中心为d3，图像分析窗口D₄(0)的中心为d4；图像分析窗口D₁(0)的中心从顶点d1逐步滑动到顶点d2，图像分析窗口D₂(0)的中心从顶点d2逐步滑动到顶点d4，图像分析窗口D₃(0)的中心从顶点d3逐步滑动到顶点d1，图像分析窗口D₄(0)的中心从顶点d4逐步滑动到顶点d3。等滑动步长是指：将4段中心线d1→d2，d2→d4，d4→d3，d3→d1按照其自身的长度分别均分为N_θ段，N_θ为当前顶点滑动到下一顶点的滑动步数，4段中心线各自分别的均分长度即为四个图像分析窗口D_k(0)每次滑动的步长；令θ=0，1，…，N_θ-1，四个图像分析窗口D_k(0)按照θ的取值，每次按照自身的步长沿中心线滑动相对相同的距离，参数θ表示滑动过程中的步数。 

假设滑动到相对位置θ时的四个图像分析窗口记为D_k(θ)，k=1，2，3，4；计算在θ位置时的四个图像分析窗口记为D_k(θ)的整体显著性度量F(θ)的具体方法如步骤S24所述，此处不再赘述。 

假设D₁(θ)，D₂(θ)，D₃(θ)，D₄(θ)对应从d1→d2，d2→d4，d4→d3，d3→d1滑动的次数为N_θ，每一步滑动由步骤S24中的方法计算得到D_k(θ)的整体显著性度量为F(θ)，则移动过程中的最大整体显著性度量F(θ*)为： 

F(θ*)=max{F(θ)|θ=0，1，…，N_θ-1}       （4） 

该最大整体显著性度量F(θ*)所对应的位置θ*为： 

${\mathrm{\&theta;}}^{*}=\underset{\mathrm{\&theta;}}{\arg \max }\left\{F\left(\mathrm{\&theta;}\right)\right|\mathrm{\&theta;}=\mathrm{0,1},...,N_{\mathrm{\&theta;}}-1\}---\left(5\right)$

S26：获取具有最大整体显著性度量的位置时、图像分析窗口所处的区域D_k(θ*)，k=1，2，3，4，及其有效子图像块的集合S_k(θ*)中、具有最大显著性度量的图像块S_km(θ*)，得到的四个图像块S_b1=S_1m(θ*)，S_b2=S_2m(θ*)，S_b3=S_3m(θ*)，S_b4=S_4m(θ*)即为最终获得的四个最佳图像子块。 

由前述步骤S21至步骤S25可见：1、算法在有效环形区域上进行循环优化，保证了用于后续匹配的四个子图块整体上具有较好匹配显著性特性，可以 有效提高图像块互信息配准的鲁棒性；2、由前述方法得到的四个最佳图像子块是沿图像外边界的分布均匀的四个图像分析窗口中得到的，保证了四个最佳图像子块在大范围内处在一个近似正方形的四边形顶点上，可有效减轻投影变换过程中带来的误差。 

S3：将四个最佳图像子块分别与每一待配准图像进行双尺度空间的亚像素级匹配，得到与每一待配准图像对应的四对匹配点。 

本发明中，由于待配准的图像是经初步几何校正后的多光谱图像，这些图像之间的像差中存在微小的空间平移误差、旋转误差和尺度缩放误差。这些像差的具体特点是：空间平移误差在局部观察比较明显；旋转误差在局部可以观察到，但一般较小；尺度缩放误差在局部基本观察不到，但在大尺度上不同区域对比观察时，尺度缩放误差体现在区域随位置变化而具有不同的平移误差。另外，考虑到多光谱图像的各个光谱图像之间具有多模成像的特点，因此使用参考图像和待配准图像之间的互信息最大的方法实现配准以寻找4对匹配点对。为了描述更清楚，首先简述互信息最大方法进行图像匹配的过程： 

假设参考图像为A，待配准图像为B。参考图像为A作为悬浮图在待配准图像B上滑动，滑动过程中在待配准图像B中与参考图像A在相同窗口中的子图像块为Ba。参考图像A的直方图信息熵为Ｈ(A)，子图像块Ba的直方图信息熵为Ｈ(Ba)，参考图像Ａ和子图像块Ba的联合直方图信息熵为Ｈ(A，Ba)，则参考图像Ａ和子图像块Ba的互信息I(A，Ba)为： 

I(A，Ba)=Ｈ(A)+Ｈ(Ba)-Ｈ(A，Ba)        （6） 

在参考图像Ａ在待配准图像B上滑动过程中取得最大值I(A，Ba)的子图像块Ba所在的位置就是参考图像A在待配准图像B中的匹配位置。 

为了实现亚像素级的配准精度，本发明使用双尺度多分辨互信息匹配的方法，如图6所示，步骤S3进一步包括以下步骤： 

S31：计算最佳图像子块S_bk的中心在基准图像中的位置坐标P_bk=(X_bk，Y_bk)，根据位置坐标构造相应的最佳图像子块S_bk的第一旋转子图像集合S_k ^r，之后对 第一旋转子图像集合S_k ^r中的每一旋转子图像进行插值放大，插值放大后的每一旋转子图像构成第二旋转子图像集合S_k ^zr。 

本发明中，最佳图像子块S_bk的第一旋转子图像集合S_k ^r={S_bk(φ)}是最佳图像子块S_bk的围绕其中心旋转插值后的一系列图像块的集合。在实际中为了克服图像旋转后的边界不完整问题，在基准图像中以最佳图像子块S_bk的中心坐标P_bk为中心选取一块适当大于最佳图像子块S_bk的图像块进行一系列旋转，然后在旋转后的图像的中心截取与最佳图像子块S_bk具有相同大小的图像块作为最佳图像子块S_bk旋转后的子图。第一旋转子图像集合S_k ^r={S_bk(φ)}中的图像子块数量由旋转匹配精度和原图像的最大失配角度决定。例如，若最大失配角度为±5°，旋转匹配精度为0.5°，则第一旋转子图像集合S_k ^r中有25个旋转子图像，分别对应角度φ={-5°，-4.5°，…，-0.5°，0，0.5°，…，4.5°，5°}。 

本发明中，插值放大倍率为z的放大后的第二旋转子图像集合S_k ^zr={S_bk ^z(φ)}中的每一个子图像S_bk ^z(φ)是第一旋转子图像集合S_k ^r中对应的子图像S_bk(φ)插值放大z倍后的图像。插值放大倍数z决定了配准精度，如z分别为4、8、16时，则配准的精度分别为0.25、0.125、0.0625像素。 

S32：在每一待配准图像中，以四个位置坐标P_bk=(X_bk，Y_bk)为中心确定相应的四个第一坐标区域D_ck=[X_bk-dX，X_bk+dX]×[Y_bk-dY，Y_bk+dY]，并令第一旋转子图像集合S_k ^r中的每一旋转子图像的中心在第一坐标区域D_ck中进行滑动互信息匹配计算，得到每一待配准图像的四个粗匹配坐标点P_ck=(X_ck，Y_ck)，k=1，2，3，4。 

本发明中，第一坐标区域D_ck=[X_bk-dX，X_bk+dX]×[Y_bk-dY，Y_bk+dY]中的dX和dY由最大平移失配像差决定，一般可取为最大平移失配误差的1.5倍左右。例如，若最大平移失配像差为20，则可选取dX=dY=30。由于不同光谱图像之间的平移失配像差不同，因此dX和dY也可针对不同波段选取不同的值。 

本发明中，若假设S_n(x，y)是待配准图像中的中心坐标位置在(x，y)的且与基准图像中的最佳图像子块S_bk具有相同大小的图像子块，又假设I(S_bk(φ)，S_n(x，y))是S_n(x，y)和S_bk的互信息，则粗匹配坐标点P_ck=(X_ck，Y_ck)可由下面公式（7）和 公式（8）计算得到： 

$P_{\mathrm{ck}}=(X_{\mathrm{ck}},Y_{\mathrm{ck}})=\underset{(x,y)}{\arg \max }\left\{I(x,y)\right|(x,y)\&Element;D_{\mathrm{ck}}\}---\left(8\right)$

S33：在每一待配准图像中，以四个粗匹配坐标点P_ck=(X_ck，Y_ck)为中心确定相应的四个第一图像区域块C_nk，之后对四个第一图像区域块C_nk分别进行插值放大，得到每一待配准图像的四个第二图像区域块C_nzk，k=1，2，3，4。 

本发明中，以粗匹配坐标点P_ck=(X_ck，Y_ck)为中心确定的第一图像区域块C_nk的尺寸只需要比基准图像中的最佳图像子块S_bk的尺寸大4到6个像素即可。插值放大方法以及放大倍率z要求与步骤S31中构造第二旋转子图像集合S_k ^zr时使用方法及倍率参数一致。 

S34：在每一待配准图像中，以粗匹配坐标点P_ck=(X_ck，Y_ck)为中心，在相应的第一图像区域块C_nk中确定第二坐标区域D_czk=[X_ck-dcX，X_ck+dcX]×[Y_ck-dcY，Y_ck+dcY]，并令第二旋转子图像集合S_k ^zr中的子图像中心在第二坐标区域D_czk中进行互信息滑动匹配，得到每一待配准图像的四个亚像素级的精细匹配坐标点P_nk=(X_nk，Y_nk)。 

本发明中，第二坐标区域D_czk=[X_ck-dcX，X_ck+dcX]×[Y_ck-dcY，Y_ck+dcY]中的dcX和dcY一般可取放大倍率为z的1.5倍。例如若z=16，则dcX=dcY=24。 

本发明中，第二旋转子图像集合S_k ^zr中的子图像在第二坐标区域D_czk中的互信息匹配方法，以及亚像素级精细匹配坐标点P_nk=(X_nk，Y_nk)的计算方法与步骤S32中的方法类似，不再赘述。 

S35：由待配准图像中的四个精细匹配坐标点与基准图像中的四个最佳图像子块的中心位置坐标分别进行一一对应组合，得到与每一待配准图像对应的四对匹配点。 

本发明中，由计算得到的第n个待配准图像中的4个匹配点P_nk=(X_nk，Y_nk)分别与基准图像中的四个最佳图像子块的中心位置坐标点P_bk=(X_bk，Y_bk)进行组合，得到待配准图像与基准图像的4对匹配点对(P_b1，P_n1)、(P_b2，P_n2)、(P_b3，P_n3)、 (P_b4，P_n4)。 

S4：根据四对匹配点，计算相应的待配准图像到基准图像的投影变换矩阵Mnb，并将相应的待配准图像投影到基准图像的坐标空间，之后采用必要的图像插值后，得到在基准图像的坐标空间配准和校正后的多光谱图像。 

本发明中，由4对匹配点对(P_b1，P_n1)、(P_b2，P_n2)、(P_b3，P_n3)、(P_b4，P_n4)计算光谱图像到基准图像的投影变换矩阵Mnb的方法，以及将光谱图像投影到基准图像的坐标空间的方法，以及必要的图像插值方法等，对图像处理和计算机视觉领域的技术人员是习知内容，此处不再赘述。 

图7示出了本发明提供的多光谱遥感图像的配准与校正系统的结构，为了便于说明，仅示出了与本发明相关的部分。 

具体地，本发明提供的多光谱遥感图像的配准与校正系统包括：选择模块1，用于选择多光谱遥感图像中的一个光谱图像作为基准图像，而将多光谱遥感图像中的其它光谱图像作为待配准图像，该多光谱遥感图像包括多个经初步几何校正的光谱图像；提取模块2，用于在基准图像中提取具有整体最佳匹配特征的四个最佳图像子块；匹配模块3，用于将四个最佳图像子块分别与每一待配准图像进行双尺度空间的亚像素级匹配，得到与每一待配准图像对应的四对匹配点；配准及校正模块4，用于根据四对匹配点，计算相应的待配准图像到基准图像的投影变换矩阵，并将相应的待配准图像投影到基准图像的坐标空间，之后采用必要的图像插值后，得到在基准图像的坐标空间配准和校正后的多光谱图像。 

进一步地，如图8所示，提取模块2可包括：第一计算子模块21，用于确定基准图像中的有效图像区域，并在有效图像区域内，沿有效图像区域外边界确定一有效环形区域，并计算有效环形区域的中心线，以及中心线上构成一四边形的四个顶点；标定子模块22，用于对基准图像进行网格化规则分块，之后对分块后各子块中的有效子图像块进行标定，获得基准图像的有效子图像块的集合；第二计算子模块23，用于分别以中心线上的四个顶点为中心确定基准图 像上的四个图像分析窗口；第四计算子模块24，用于利用基准图像的有效子图像块的集合，计算四个图像分析窗口分别的有效子图像块的集合，计算图像分析窗口的有效子图像块的集合中每一图像块的显著性度量，计算每一图像分析窗口的显著性度量，并计算四个图像分析窗口的整体显著性度量；第五计算子模块25，用于将四个图像分析窗口的中心分别沿中心线、按顺时针或逆时针方向以等滑动步长逐步滑动到相邻的下一顶点，并计算滑动过程中的最大整体显著性度量，以及最大整体显著性度量所对应的位置；获取子模块26，用于获取具有最大整体显著性度量的位置时、图像分析窗口所处的区域，及其有效子图像块的集合中、具有最大显著性度量的图像块，得到的四个图像块即为最终获得的四个最佳图像子块。其中，各子模块的执行过程如上所述，在此不赘述。 

更进一步地，如图9所示，第五计算子模块25可包括：划分子模块251，用于将中心线被四个顶点分割后的四条中间线分别均分为若干段，以得到相应的滑动步长；滑动子模块252，用于将四个图像分析窗口的中心分别沿中心线、按顺时针或逆时针方向以相应的滑动步长向相邻的下一顶点滑动；第六计算子模块253，用于当图像分析窗口每滑动固定步数时，计算相应的整体显著性度量，进而得到相应的图像分析窗口滑动到下一顶点过程中的最大整体显著性度量，以及最大整体显著性度量所对应的位置。其中，各子模块的执行过程如上所述，在此不赘述。 

进一步地，如图10所示，匹配模块3可包括：构造子模块31，用于计算最佳图像子块的中心在基准图像中的位置坐标，根据位置坐标构造相应的最佳图像子块的第一旋转子图像集合，之后对第一旋转子图像集合中的每一旋转子图像进行插值放大，插值放大后的每一旋转子图像构成第二旋转子图像集合；粗匹配子模块32，用于在每一待配准图像中，以四个位置坐标为中心确定相应的四个第一坐标区域，并令第一旋转子图像集合S_k ^r中的每一旋转子图像的中心在第一坐标区域中进行滑动互信息匹配计算，得到每一待配准图像的四个粗匹配坐标点；第七计算子模块33，用于在每一待配准图像中，以四个粗匹配坐标 点为中心确定相应的四个第一图像区域块，之后对四个第一图像区域块分别进行插值放大，得到每一待配准图像的四个第二图像区域块；精准匹配子模块34，用于在每一待配准图像中，以粗匹配坐标点为中心，在相应的第一图像区域块中确定第二坐标区域，并令第二旋转子图像集合中的子图像中心在第二坐标区域中进行互信息滑动匹配，得到每一待配准图像的四个亚像素级的精细匹配坐标点；组合子模块35，用于由待配准图像中的四个精细匹配坐标点与基准图像中的四个最佳图像子块的中心位置坐标分别进行一一对应组合，得到与每一待配准图像对应的四对匹配点。其中，各子模块的执行过程如上所述，在此不赘述。 

本发明提出的多光谱遥感图像的配准与校正方法是选取经初步几何校正的多光谱遥感图像中的一个作为基准图像，其余作为待配准图像，首先在基准图像中提取最佳图像子块，之后将最佳图像子块与每一待配准图像进行双尺度空间的亚像素级匹配以得到匹配点，之后根据匹配点计算投影变换矩阵，并将相应的光谱图像投影到基准图像的坐标空间，之后采用必要的图像插值后，得到在基准图像的坐标空间配准和校正后的多光谱图像。相对于现有的多光谱遥感图像配准方法，计算效率高，配准准确度高，特别适合于大幅面的多光谱遥感图像的配准及校正。 

本领域普通技术人员可以理解实现上述实施例方法中的全部或部分步骤是可以通过程序来控制相关的硬件完成，所述的程序可以在存储于一计算机可读取存储介质中，所述的存储介质，如ROM/RAM、磁盘、光盘等。 

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。 

Claims (10)

1.一种多光谱遥感图像的配准与校正方法，其特征在于，所述方法包括以下步骤：

选择多光谱遥感图像中的一个光谱图像作为基准图像，而将所述多光谱遥感图像中的其它光谱图像作为待配准图像，所述多光谱遥感图像包括多个经初步几何校正的光谱图像；

在所述基准图像中提取具有整体最佳匹配特征的四个最佳图像子块；

将所述四个最佳图像子块分别与每一待配准图像进行双尺度空间的亚像素级匹配，得到与每一待配准图像对应的四对匹配点；

根据所述四对匹配点，计算相应的待配准图像到所述基准图像的投影变换矩阵，并将相应的待配准图像投影到所述基准图像的坐标空间，之后采用必要的图像插值后，得到在所述基准图像的坐标空间配准和校正后的多光谱图像。

2.如权利要求1所述的多光谱遥感图像的配准与校正方法，其特征在于，所述在所述基准图像中提取具有整体最佳匹配特征的四个最佳图像子块的步骤进一步包括以下步骤：

确定所述基准图像中的有效图像区域，并在所述有效图像区域内，沿有效图像区域外边界确定一有效环形区域，并计算所述有效环形区域的中心线，以及所述中心线上构成一四边形的四个顶点；

对所述基准图像进行网格化规则分块，之后对分块后各子块中的有效子图像块进行标定，获得所述基准图像的有效子图像块的集合；

分别以所述四个顶点为中心确定所述基准图像上的四个图像分析窗口；

利用所述基准图像的所述有效子图像块的集合，计算所述四个图像分析窗口分别的有效子图像块的集合，计算所述图像分析窗口的所述有效子图像块的集合中每一图像块的显著性度量，计算每一图像分析窗口的显著性度量，并计算所述四个图像分析窗口的整体显著性度量；

将所述四个图像分析窗口的中心分别沿所述中心线、按顺时针或逆时针方向以等滑动步长逐步滑动到相邻的下一顶点，并计算滑动过程中的最大整体显著性度量，以及所述最大整体显著性度量所对应的位置；

获取具有所述最大整体显著性度量的位置时、所述图像分析窗口所处的区域及其有效子图像块的集合中、具有所述最大显著性度量的图像块，得到的四个图像块即为所述四个最佳图像子块。

3.如权利要求2所述的多光谱遥感图像的配准与校正方法，其特征在于，所述计算所述图像分析窗口的所述有效子图像块的集合中每一图像块的显著性度量的步骤表示为：

$F={|m4-m3|}^2\&CenterDot;\frac{{|m2-m1|}^2}{1+{\mathrm{\&sigma;}}_1^2+{\mathrm{\&sigma;}}_2^2}$

其中，F是所述每一图像块的显著性度量，Ｓ是所述每一图像块，m0是所述图像块Ｓ的所有像素灰度值的均值，m1是所述图像块Ｓ中灰度值小于所述m0的所有像素灰度值的均值，m2是所述图像块Ｓ中灰度值大于所述m0的所有像素灰度值的均值，σ₁是所述图像块Ｓ中灰度值小于所述m0的所有像素灰度值的标准差，σ₂是所述图像块Ｓ中灰度值大于所述m0的所有像素灰度值的标准差，m3是所述图像块Ｓ中灰度值小于所述m1的所有像素灰度值的均值，m4是所述图像块Ｓ中灰度值大于所述m2的所有像素灰度值的均值；

所述每一图像分析窗口的显著性度量是相应图像分析窗口中的所有所述有效子图像块的显著性度量的最大值；

所述四个图像分析窗口的整体显著性度量是所述四个图像分析窗口的显著性度量的最小值。

4.如权利要求2所述的多光谱遥感图像的配准与校正方法，其特征在于，所述将所述四个图像分析窗口的中心分别沿所述中心线、按顺时针或逆时针方向以等滑动步长逐步滑动到相邻的下一顶点，并计算滑动过程中的最大整体显著性度量，以及所述最大整体显著性度量所对应的位置的步骤进一步包括以下步骤：

将所述中心线被所述四个顶点分割后的四条中间线分别均分为若干段，以得到相应的滑动步长；

将所述四个图像分析窗口的中心分别沿所述中心线、按顺时针或逆时针方向以相应的滑动步长向相邻的下一顶点滑动；

所述图像分析窗口每滑动固定步数，计算相应的整体显著性度量，进而得到相应的图像分析窗口滑动到下一顶点过程中的最大整体显著性度量，以及所述最大整体显著性度量所对应的位置。

5.如权利要求4所述的多光谱遥感图像的配准与校正方法，其特征在于，所述得到相应的图像分析窗口滑动到下一顶点过程中的最大整体显著性度量，以及所述最大整体显著性度量所对应的位置的步骤表示为：

F(θ*)=max{F(θ)|θ=0，1，…，N_θ-1}

${\mathrm{\&theta;}}^{*}=\underset{\mathrm{\&theta;}}{\arg \max }\left\{F\left(\mathrm{\&theta;}\right)\right|\mathrm{\&theta;}=\mathrm{0,1},...,N_{\mathrm{\&theta;}}-1\}$

其中，F(θ*)为所述最大整体显著性度量，θ*为所述最大整体显著性度量所对应的位置，θ为滑动过程中的步数，F(θ)为在θ位置时的整体显著性度量，N_θ为当前顶点滑动到下一顶点的滑动步数。

6.如权利要求1所述的多光谱遥感图像的配准与校正方法，其特征在于，所述将所述四个最佳图像子块分别与每一待配准图像进行双尺度空间的亚像素级匹配，得到与每一待配准图像对应的四对匹配点的步骤进一步包括以下步骤：

计算所述最佳图像子块的中心在所述基准图像中的位置坐标，根据所述位置坐标构造相应的最佳图像子块的第一旋转子图像集合，之后对所述第一旋转子图像集合中的每一旋转子图像进行插值放大，插值放大后的所述每一旋转子图像构成第二旋转子图像集合；

在每一待配准图像中，以四个所述位置坐标为中心确定相应的四个第一坐标区域，并令所述第一旋转子图像集合中的每一旋转子图像的中心在所述第一坐标区域中进行滑动互信息匹配计算，得到每一待配准图像的四个粗匹配坐标点；

在每一待配准图像中，以所述四个粗匹配坐标点为中心确定相应的四个第一图像区域块，之后对四个所述第一图像区域块分别进行插值放大，得到每一待配准图像的四个第二图像区域块；

在每一待配准图像中，以所述粗匹配坐标点为中心，在相应的所述第一图像区域块中确定第二坐标区域，并令所述第二旋转子图像集合中的子图像中心在所述第二坐标区域中进行互信息滑动匹配，得到每一待配准图像的四个亚像素级的精细匹配坐标点；

由所述待配准图像中的四个所述精细匹配坐标点与所述基准图像中的四个所述最佳图像子块的中心位置坐标分别进行一一对应组合，得到与每一待配准图像对应的四对匹配点。

7.一种多光谱遥感图像的配准与校正系统，其特征在于，所述系统包括：

选择模块，用于选择多光谱遥感图像中的一个光谱图像作为基准图像，而将所述多光谱遥感图像中的其它光谱图像作为待配准图像，所述多光谱遥感图像包括多个经初步几何校正的光谱图像；

提取模块，用于在所述基准图像中提取具有整体最佳匹配特征的四个最佳图像子块；

匹配模块，用于将所述四个最佳图像子块分别与每一待配准图像进行双尺度空间的亚像素级匹配，得到与每一待配准图像对应的四对匹配点；

配准及校正模块，用于根据所述四对匹配点，计算相应的待配准图像到所述基准图像的投影变换矩阵，并将相应的待配准图像投影到所述基准图像的坐标空间，之后采用必要的图像插值后，得到在所述基准图像的坐标空间配准和校正后的多光谱图像。

8.如权利要求7所述的多光谱遥感图像的配准与校正系统，其特征在于，所述提取模块包括：

第一计算子模块，用于确定所述基准图像中的有效图像区域，并在所述有效图像区域内，沿有效图像区域外边界确定一有效环形区域，并计算所述有效环形区域的中心线，以及所述中心线上构成一四边形的四个顶点；

标定子模块，用于对所述基准图像进行网格化规则分块，之后对分块后各子块中的有效子图像块进行标定，获得所述基准图像的有效子图像块的集合；

第二计算子模块，用于分别以所述中心线上的四个顶点为中心确定所述基准图像上的四个图像分析窗口；

第四计算子模块，用于利用所述基准图像的有效子图像块的集合，计算所述四个图像分析窗口分别的有效子图像块的集合，计算所述图像分析窗口的有效子图像块的集合中每一图像块的显著性度量，计算每一图像分析窗口的显著性度量，并计算所述四个图像分析窗口的整体显著性度量；

第五计算子模块，用于将所述四个图像分析窗口的中心分别沿所述中心线、按顺时针或逆时针方向以等滑动步长逐步滑动到相邻的下一顶点，并计算滑动过程中的最大整体显著性度量，以及所述最大整体显著性度量所对应的位置；

获取子模块，用于获取具有所述最大整体显著性度量的位置时、所述图像分析窗口所处的区域及其有效子图像块的集合中、具有最大显著性度量的图像块，得到的四个图像块即为所述四个最佳图像子块。

9.如权利要求8所述的多光谱遥感图像的配准与校正系统，其特征在于，所述第五计算子模块包括：

划分子模块，用于将所述中心线被所述四个顶点分割后的四条中间线分别均分为若干段，以得到相应的滑动步长；

滑动子模块，用于将所述四个图像分析窗口的中心分别沿所述中心线、按顺时针或逆时针方向以相应的滑动步长向相邻的下一顶点滑动；

第六计算子模块，用于当所述图像分析窗口每滑动固定步数时，计算相应的整体显著性度量，进而得到相应的图像分析窗口滑动到下一顶点过程中的最大整体显著性度量，以及所述最大整体显著性度量所对应的位置。

10.如权利要求7所述的多光谱遥感图像的配准与校正系统，其特征在于，所述匹配模块包括：

构造子模块，用于计算所述最佳图像子块的中心在所述基准图像中的位置坐标，根据所述位置坐标构造相应的最佳图像子块的第一旋转子图像集合，之后对所述第一旋转子图像集合中的每一旋转子图像进行插值放大，插值放大后的所述每一旋转子图像构成第二旋转子图像集合；

粗匹配子模块，用于在每一待配准图像中，以所述四个位置坐标为中心确定相应的四个第一坐标区域，并令所述第一旋转子图像集合中的每一旋转子图像的中心在所述第一坐标区域中进行滑动互信息匹配计算，得到每一待配准图像的四个粗匹配坐标点；

第七计算子模块，用于在每一待配准图像中，以所述四个粗匹配坐标点为中心确定相应的四个第一图像区域块，之后对所述四个第一图像区域块分别进行插值放大，得到每一待配准图像的四个第二图像区域块；

精准匹配子模块，用于在每一待配准图像中，以所述粗匹配坐标点为中心，在相应的第一图像区域块中确定第二坐标区域，并令所述第二旋转子图像集合中的子图像中心在所述第二坐标区域中进行互信息滑动匹配，得到每一待配准图像的四个亚像素级的精细匹配坐标点；

组合子模块，用于由待配准图像中的四个所述精细匹配坐标点与四个所述基准图像中的四个所述最佳图像子块的中心位置坐标分别进行一一对应组合，得到与每一待配准图像对应的四对匹配点。

Images