CN104732532A - 一种遥感卫星多光谱图像配准方法 - Google Patents

Abstract

本发明一种遥感卫星多光谱图像配准方法，该方法包括以下步骤：①根据多光谱探测器的距离得到粗配准图像；②基于基准图像对待配准图像进行直方图映射；③分块划分上述两幅图像，针对每对图像块：提取每对图像块的角点、近红外图像的配准还需要提取各块的主方向相位一致特征、基于相位相关和曲面拟合得到各块的正反向亚像素匹配点集；④三角网格划分基准图像和待配准图像，计算变换系数并对待配准图像进行重采样得到配准图像。与现有技术相比，本发明方法对不同反射率地物在不同谱段上表现出的灰度差异大而导致的配准难问题，同时有效解决了大幅面遥感图像不同区域存在不同程度畸变的问题。

Description

一种遥感卫星多光谱图像配准方法

技术领域

本发明涉及一种遥感卫星多光谱图像配准方法，属于遥感图像处理技术领域，特别适用于遥感卫星多光谱图像的配准。

背景技术

随着用户对遥感数据需求的增长以及遥感技术的发展，多光谱相机在遥感卫星上得到广泛应用。多光谱相机具有较宽谱段的探测能力，通常从蓝谱段到近红外谱段，一般包括蓝、绿、红、近红外等谱段，例如，CBERS-1的四个谱段分别为0.45-0.52、0.52-0.59、0.63-0.69、0.77-0.89(单位：μm)，QuickBird的四个谱段分别为0.45-0.52、0.52-0.60、0.63-0.69、0.76-0.90(单位：μm)，SPOT5的四个谱段分别为0.50-0.59、0.61-0.68、0.79-0.89、1.58-1.70(单位：μm)。

通过多光谱多个谱段图像的配准能够实现多光谱图像的融合，合成较高质量的彩色图像，进一步通过对彩色图像的分析，可以为资源普查、灾害监测等提供定量应用的数据。然而，不同谱段图像间的微小位置差异就能导致合成的彩色图像出现“重叠影”，严重影响图像质量，为提高多光谱图像的融合质量，在处理过程中需要各谱段间图像的严格配准。常用的图像配准方法主要包括基于图像灰度与基于图像特征的配准方法，其中，基于图像灰度的方法直接利用图像的灰度信息，建立两幅图像之间的相似性度量，但对于灰度差异大的两幅图像配准，该方法失去价值；基于图像特征的方法，通过提取基准图像和待配准图像的一些共同特征，建立特征之间的对应关系，但大多数的角点、边缘特征提取主要依赖于图像灰度信息，导致基准图像和待配准图像提取的特征很难一致。

由于地面物质对不同波段光的反射率不同，不同波段的图像灰度差别也较大，特别是蓝、绿、红谱段与近红外谱段间的灰度差异更大。而现有的基于灰度与基于特征的配准方法对灰度信息比较敏感，导致配准难度加大，同时由于大气、相机、卫星系统等的影响导致不同谱段间存在一定的畸变。

发明内容

本发明解决的技术问题是：克服现有技术的不足，提供一种遥感卫星多光谱图像配准方法，用于解决遥感卫星多光谱谱段间的大幅面图像配准，特别解决了不同反射率地物在近红外谱段与蓝、绿、红谱段图像中表现出灰度差异大而导致的配准难问题。

本发明的技术方案是：一种遥感卫星多光谱图像配准方法，步骤如下：

1)根据多光谱相机蓝、绿、红、近红外四个谱段探测器间的距离，对获取的四个谱段图像分别进行裁剪，得到粗配准的多光谱图像，具体步骤如下：

多光谱相机焦平面上的蓝谱段探测器到绿、红、近红外谱段探测器的距离分别为d₁、d₂、d₃，则从绿、红、近红外谱段图像的顶部分别裁剪δ₁、δ₂、δ₃个像素，从蓝、绿、红谱段图像的底部分别裁剪δ₃、δ₂、δ₁个像素，其中， τ为多光谱相机的单个探测器尺寸，表示下取整，裁剪后的图像即为粗配准的多光谱图像；

2)将粗配准的红谱段图像作为基准图像f(x,y)，粗配准的蓝、绿、近红外图像中的任意一幅图像作为待配准图像，建立基准图像到待配准图像的直方图映射关系，以建立的直方图映射关系对待配准图像进行映射，得到映射图像g(x,y)，所述f(x,y)与g(x,y)的大小为M×N个像素；M、N为正整数；

3)以大小为m×m个像素的图像块分别划分f(x,y)与g(x,y)，得到的每对图像块表示为f_i与g_i，共得到对图像块，其中，m为正整数，表示上取整；针对每一对图像块f_i与g_i进行如下步骤：

3.1)分别提取图像块f_i与g_i的Harris角点，图像块f_i的所有角点构成f_i的正向角点集；

3.2)当待配准图像为近红外谱段图像时，计算图像块f_i的梯度主方向Φ_i；构造方向为Φ_i的二维Log-Gabor滤波器对f_i与g_i进行相位一致特征检测，分别得到f_i与g_i的主方向相位一致特征PC(f_i)与PC(g_i)；

3.3)将图像块f_i的正向角点集中的每个角点作为待匹配点，在图像块g_i中基于相位相关法和二元二次曲面拟合计算对应的亚像素匹配点，则得到在图像块g_i中的正向亚像素匹配角点集；正向角点集和正向亚像素匹配角点集记为正向匹配角点对集Γ_i；

3.4)计算从g_i提取的每个角点与正向亚像素匹配角点集中所有角点的最小距离，如果最小距离小于设定阈值则剔除g_i中的该角点，否则保留g_i中的该角点，所有符合上述条件的角点构成g_i中的反向角点集；

3.5)将g_i中的反向角点集中的每个角点作为待匹配点，在图像块f_i中基于相位相关法和二元二次曲面拟合计算对应的亚像素匹配点，则得到在f_i中的反向亚像素匹配角点集；反向角点集和反向亚像素匹配角点集构成反向匹配角点对集Γ'_i；正反向匹配角点对集表示为Γ_i+Γ'_i；

3.6)将图像块f_i划分若干个大小为K×K个像素的块区域，K为正整数；遍历所有块区域，判断块区域内是否有匹配的角点对，即是否有集合Γ_i+Γ'_i中的元素，若没有则将该块区域的中心点作为构造点，所有的构造点构成构造点集；

3.7)将构造点集中的每个构造点作为待匹配点，在图像块g_i中基于相位相关法和二元二次曲面拟合计算对应的亚像素匹配点，则得到在g_i中的亚像素匹配点集；构造点集及其亚像素匹配点集构成匹配构造点对集Γ”_i；完整的匹配点集合记为Γ_i+Γ'_i+Γ”_i；

3.8)基于RANSAC法剔除误匹配点对，误匹配点集记为Γ”'_i，则图像块f_i与g_i的匹配点对集 ${\mathrm{\Γ}}_{\mathrm{total}}^i{=\mathrm{\Γ}}_i+{{\mathrm{\Γ}}^{\′}}_i+{{\mathrm{\Γ}}^{\′\′}}_i-{{\mathrm{\Γ}}^{\′\′\′}}_i;$

3.9)重复执行步骤3.1)～步骤3.8)，直到f(x,y)和g(x,y)的对图像块都处理完毕，得到所有图像块的匹配点对集，即整幅图像的匹配点对集Γ_total；

4)在整幅图像的匹配点对集Γ_total的基础上，基于Delaunay三角网格划分基准图像和待配准图像，分别计算基准图像和待配准图像中对应三角网格的变换系数，将变换系数分别作用于待配准图像对应的三角网格内部的每个点，得到配准图像；

5)重复步骤2)～步骤4)，直到蓝、绿、近红外谱段图像都完成配准，分别得到它们的配准图像。

所述步骤3.2)中计算图像块f_i的梯度主方向Φ_i，通过以下步骤实现：

3.2.1)将[0,π)区间划分为n个方向，记为集合n为正整数；

3.2.2)计算图像块f_i中每个像素的梯度幅值mag(x,y)和梯度方向θ(x,y)；

$\mathrm{mag}(x,y)=\sqrt{{\mathrm{\Δx}}^2+{\mathrm{\Δy}}^2},\mathrm{\θ}(x,y)={\mathrm{tg}}^{-1}|\mathrm{\Δy}/\mathrm{\Δx}|;$

其中，Δx＝f_i(x+1,y)-f_i(x-1,y)，Δy＝f_i(x,y+1)-f_i(x,y-1)；当Δx＝0时，θ(x,y)＝π/2，θ(x,y)的范围为[0,π)；

3.2.3)根据划分的n个方向将步骤3.2.2)中得到的该像素的梯度幅值进行分解，具体步骤如下：

计算d＝θ(x,y)×n/π，则该像素的梯度方向介于第个方向和第个方向之间，mod表示求余运算；将该像素的梯度幅值分解到第个方向和第个方向上，得到在第个方向的梯度幅值为在第个方向的梯度幅值为

3.2.4)对图像块f_i中的所有像素重复步骤3.2.2)、步骤3.2.3)，将所有像素的梯度幅值分解到n个方向上；

3.2.5)分别计算n个方向上的幅值和，最大的幅值和记为MAG_max，对应的方向记为属于 $\{0,\frac{\mathrm{\&pi;}}{n},\frac{2\mathrm{\&pi;}}{n},\&CenterDot;\&CenterDot;\&CenterDot;,\frac{(n-1)\mathrm{\&pi;}}{n}\},$ 与相邻的两个方向表示为 ${\mathrm{\&Phi;}}_L=\left\{\begin{array}{cc}{\mathrm{\&Phi;}}_{{\mathrm{MAG}}_{\max }}-\frac{\mathrm{\&pi;}}{n}& {\mathrm{\&Phi;}}_{{\mathrm{MAG}}_{\max }}\&GreaterEqual;\frac{\mathrm{\&pi;}}{n}\\ \frac{(n-1)\mathrm{\&pi;}}{n}& {\mathrm{\&Phi;}}_{{\mathrm{MAG}}_{\max }}=0\end{array}\right.,{\mathrm{\&Phi;}}_R=\left\{\begin{array}{cc}{\mathrm{\&Phi;}}_{{\mathrm{MAG}}_{\max }}+\frac{\mathrm{\&pi;}}{n}& {\mathrm{\&Phi;}}_{{\mathrm{MAG}}_{\max }}<\frac{(n-1)\mathrm{\&pi;}}{n}\\ 0& {\mathrm{\&Phi;}}_{{\mathrm{MAG}}_{\max }}=\frac{(n-1)\mathrm{\&pi;}}{n}\end{array}\right.,$ 对应的幅值和分别为MAG_L、MAG_R，利用MAG_max、MAG_L、MAG_R及其对应的方向Φ_L、Φ_R拟合一元二次曲线，则曲线最大峰值所对应的方向即为图像块f_i的梯度主方向Φ_i。

所述步骤3.3)、步骤3.5)、步骤3.7)中任一项基于相位相关法和二元二次曲面拟合计算对应的亚像素匹配点，通过以下步骤实现：

(a)当待配准图像是蓝谱段或绿谱段图像时，以待匹配点为中心，在图像块f_i、g_i中分别截取大小为m₁×m₁个像素的块区域；当待配准图像是近红外谱段图像时，以待匹配点为中心，在主方向相位一致特征PC(f_i)与PC(g_i)中分别截取大小为m₁×m₁个像素的块区域，m₁取128或256；

(b)基于相位相关计算两个块区域的相位相关系数，最大相位相关系数对应的平移量即为两个块区域的平移量，得到像素级匹配点；

(c)以像素级匹配点为中心的局部区域内的相位相关系数拟合二元二次曲面，所述局部区域的大小为3×3个像素；计算曲面最大值所对应的位置即为亚像素匹配点。

本发明与现有技术相比有益效果为：

(1)本发明提出一整套提取较均匀分布的角点、计算基准图像与待配准图像间的正反向匹配角点对集合、在无角点区域构造点并计算匹配点对集合、划分三角网格并计算变换系数等步骤的自动配准技术方案。

(2)本发明提出基于相位一致特征相位相关的配准方法，首先构造主方向的滤波器提取相位一致特征，然后通过相位相关和局部区域曲面拟合获得鲁棒的亚像素匹配精度，解决了近红外与蓝、绿、红谱段图像灰度差异大带来的配准难问题，同时仅提取主方向的相位一致特征减少了计算量。

(3)本发明提出从基准图像到待配准图像的正向匹配角点集、与从待配准图像到基准图像的反向匹配角点集的双向策略、以及在无明显特征区域构造点并计算匹配点对，可控制整幅图像划分的三角网格密度，有效地解决了由于地形、相机等导致的大幅面图像的不同区域存在不同程度畸变的问题。

附图说明

图1是本发明遥感卫星多光谱图像配准方法流程图；

图2a是本发明多光谱相机焦平面上的蓝谱段探测器到绿、红、近红外谱段探测器的距离，图2b是本发明对四个谱段图像进行裁剪获得粗配准图像的示意图；

图3是本发明有着显著灰度差异的红外与近红外谱段图像的角点匹配结果，其中图3a是红外谱段图像的结果，图3b是近红外谱段图像的结果；

图4是基于本发明方法的红外与近红外谱段图像的匹配结果及划分的三角网格，其中图4a是红外谱段图像的结果，图4b是近红外谱段图像的结果。

具体实施方式

一种遥感卫星多光谱图像配准方法，具体步骤如图1所示，该方法由以下步骤实现：

1、根据多光谱相机蓝、绿、红、近红外四个谱段探测器间的距离，对获取的四个谱段图像分别进行裁剪，得到粗配准的多光谱图像。

对多光谱相机采用分光棱镜获取多光谱图像的方式，四个谱段间没有偏移量，不需要对输入的四个谱段图像进行裁剪，获取的四个谱段图像即为粗配准的多光谱图像；对多光谱相机采用多色CCD探测器集成在一个器件上获取多光谱图像的方式，如图2a所示，焦平面上的蓝谱段探测器到绿、红、近红外谱段探测器的距离分别为d₁、d₂、d₃，由于四个谱段探测器先后依次对地面同一区域成像，导致了谱段间图像的偏移量，如图2b所示，则从绿、红、近红外谱段图像的顶部分别裁剪δ₁、δ₂、δ₃个像素，从蓝、绿、红谱段图像的底部分别裁剪δ₃、δ₂、δ₁个像素，其中，τ为多光谱相机的单个探测器尺寸，表示下取整，裁剪后的图像即为粗配准的多光谱图像。

2、将粗配准的红谱段图像作为基准图像f(x,y)，粗配准的蓝、绿、近红外图像中的任意一幅图像作为待配准图像，建立基准图像到待配准图像的直方图映射关系，以建立的直方图映射关系对待配准图像进行映射，得到映射图像g(x,y)，所述f(x,y)与g(x,y)的大小为M×N个像素；M、N为正整数。

所述建立基准图像到待配准图像的直方图映射关系，其主要目的是减小两幅图像间的灰度差异，具体步骤包括：

(1)计算基准图像的直方图概率密度函数，作为期望直方图，计算期望直方图各级灰度值对应的累积概率密度。

(2)计算待配准图像的直方图概率密度函数，作为原始直方图，计算原始直方图各级灰度值对应的累积概率密度。

(3)基于直方图匹配建立原始直方图与期望直方图的直方图映射关系，将直方图映射关系作用到待配准图像，得到映射图像。

3、以大小为m×m个像素的图像块分别划分f(x,y)和g(x,y)，得到的每对图像块表示为f_i和g_i，共得到对图像块，其中，表示上取整，m一般取值500～2000的整数。针对每一对图像块f_i与g_i进行如下步骤：

3.1、分别提取图像块f_i与g_i的Harris角点，图像块f_i的所有角点构成f_i的正向角点集。

Harris算子(Harris C.“A combined corner and edge detector”,1988)只用到图像灰度的一阶差分，具有计算简单、提取的角点特征均匀合理、可以定量提取特征点，并且算法稳定性好，对图像噪声等具有鲁棒性，特征点的数量可通过修改阈值来控制，因此比较适合大幅面的卫星图像的配准。

3.2、当待配准图像为近红外谱段图像时，计算图像块f_i的梯度主方向Φ_i；构造方向为Φ_i的二维Log-Gabor滤波器，基于该滤波器对f_i与g_i进行相位一致特征检测，分别得到f_i与g_i的主方向相位一致特征PC(f_i)与PC(g_i)。

常用的空域检测算法主要采用亮度梯度值来表征边缘的强度，因此，检测出的边缘对图像亮度和对比度的变化非常敏感。1987年，Morrone等人在研究马赫带现象时发现，人类视觉感知的图像特征出现在图像Fourier谐波分量叠合最大的相位处，而特征的类型由相位的值决定。据此他提出了相位一致的理论和计算方法，并通过实验和理论证明了相位一致与人类视觉系统对图像特征的认知相符合。

Kovesi提出了改进相位一致计算方法(参见：Peter Kovesi，“Image FeaturesFrom Phase Congruency”)，但是需要多个方向(如6、8个方向)的滤波器与图像进行卷积，使得特征提取的计算量增加，对于大幅面的遥感图像提取相位一致特征耗费的内存大、时间长。本发明在Kovesi改进方法的基础上，只利用方向为Φ_i的二维Log-Gabor滤波器计算相位一致特征，极大地减少了计算量。

所述图像块f_i的梯度主方向Φ_i的计算步骤如下：

(1)将[0,π)区间划分为n个方向，记为集合n为正整数，一般取值4～10。

(2)计算图像块f_i中每个像素的梯度幅值mag(x,y)和梯度方向θ(x,y)。

$\mathrm{mag}(x,y)=\sqrt{{\mathrm{\&Delta;x}}^2+{\mathrm{\&Delta;y}}^2}$

θ(x,y)＝tg^-1|Δy/Δx|

其中，Δx＝f_i(x+1,y)-f_i(x-1,y)，Δy＝f_i(x,y+1)-f_i(x,y-1)；当Δx＝0时，θ(x,y)＝π/2，θ(x,y)的范围为[0,π)。

(3)根据划分的n个方向将步骤(2)中的得到的该像素的梯度幅值进行分解，具体步骤如下：

计算d＝θ(x,y)×n/π，则该像素的梯度方向介于第个方向和第个方向之间，表示下取整，mod表示求余运算；将该像素的梯度幅值分解到第个方向和第个方向上，得到在第个方向的梯度幅值为在第个方向的梯度幅值为

(4)对图像块f_i中的所有像素重复步骤(2)与(3)，将所有像素的梯度幅值分解到n个方向上。

(5)分别计算n个方向上的幅值和，最大的幅值和记为MAG_max，对应的方向记为属于 $\{0,\frac{\mathrm{\&pi;}}{n},\frac{2\mathrm{\&pi;}}{n},\&CenterDot;\&CenterDot;\&CenterDot;,\frac{(n-1)\mathrm{\&pi;}}{n}\},$ 与相邻的两个方向表示为 ${\mathrm{\&Phi;}}_L=\left\{\begin{array}{cc}{\mathrm{\&Phi;}}_{{\mathrm{MAG}}_{\max }}-\frac{\mathrm{\&pi;}}{n}& {\mathrm{\&Phi;}}_{{\mathrm{MAG}}_{\max }}\&GreaterEqual;\frac{\mathrm{\&pi;}}{n}\\ \frac{(n-1)\mathrm{\&pi;}}{n}& {\mathrm{\&Phi;}}_{{\mathrm{MAG}}_{\max }}=0\end{array}\right.,{\mathrm{\&Phi;}}_R=\left\{\begin{array}{cc}{\mathrm{\&Phi;}}_{{\mathrm{MAG}}_{\max }}+\frac{\mathrm{\&pi;}}{n}& {\mathrm{\&Phi;}}_{{\mathrm{MAG}}_{\max }}<\frac{(n-1)\mathrm{\&pi;}}{n}\\ 0& {\mathrm{\&Phi;}}_{{\mathrm{MAG}}_{\max }}=\frac{(n-1)\mathrm{\&pi;}}{n}\end{array}\right.,$ 对应的幅值和分别为MAG_L、MAG_R，利用MAG_max、MAG_L、MAG_R及其对应的方向Φ_L、Φ_R拟合一元二次曲线，曲线方程为MAG＝r·Φ²+s·Φ+t，MAG、Φ表示幅值和方向大小，r、s、t表示拟合得到的曲线系数，该曲线的最大峰值所对应的方向即为图像块f_i的梯度主方向Φ_i。

所述构造方向为Φ_i的二维Log-Gabor滤波器对f_i与g_i进行相位一致特征检测，步骤如下：

(1)构造方向为Φ_i的二维Log-Gabor滤波器。

为减少计算量，仅利用计算出的主方向角Φ_i构造方向为Φ_i的二维Log-Gabor滤波器G(ω,θ)，该滤波器定义在频域的极坐标上，表示形式为：

$G(\mathrm{\&omega;},\mathrm{\&theta;})=G\left(\mathrm{\&omega;}\right)\&CenterDot;G\left(\mathrm{\&theta;}\right)=\exp \left(\frac{{-\left(\log (\mathrm{\&omega;}/{\mathrm{\&omega;}}_0)\right)}^2}{{2\left(\log ({\mathrm{\&sigma;}}_{\mathrm{\&omega;}}/{\mathrm{\&omega;}}_0)\right)}^2}\right)\exp \left(\frac{{-\left({\mathrm{\&theta;}-\mathrm{\&Phi;}}_i\right)}^2}{{{2\mathrm{\&sigma;}}_{\mathrm{\&theta;}}}^2}\right)$

其中G(ω)为径向成分，ω₀为中心频率，σ_ω用于确定径向带宽；G(θ)为角度成分，Φ_i为滤波器的方向角，σ_θ是角度方向高斯函数的标准差，用于确定方向带宽。

可以看出，二维Log-Gabor滤波器在对数坐标系下为高斯函数，没有直流分量，并且包含更多的高频成分。本发明只利用计算出的主方向角建立滤波器，较传统相位一致特征提取的多个方向滤波器参与计算，在最大可能保留信息的同时，能够极大地减少计算量。

(2)利用构造的主方向滤波器对f_i与g_i进行相位一致特征检测。

Kovesi提出的改进相位一致表示为(参见：Peter Kovesi，“Image FeaturesFrom Phase Congruency”)：

式中，Φ_i表示二维滤波器的方向，即计算得到的梯度主方向，n表示滤波器的尺度；为相位一致的权值函数；为相位偏离函数；为经过计算估计出的在方向Φ_i上的噪声干扰阈值；为在方向Φ_i、尺度n的幅值。

首先，将频域的二维滤波器G(ω,θ)变换到实域，则二维滤波器包括偶对称滤波器和奇对称滤波器两部分，它们与图像块f_i或g_i的卷积表示为：

$[e_{{\mathrm{n\&Phi;}}_i}\left(f_i\right),o_{{\mathrm{n\&Phi;}}_i}\left(f_i\right)]=[f_i*M_{{\mathrm{n\&Phi;}}_i}^{\mathrm{even}},f_i*M_{{\mathrm{n\&Phi;}}_i}^{\mathrm{odd}}]$

$[e_{{\mathrm{n\&Phi;}}_i}\left(g_i\right),o_{{\mathrm{n\&Phi;}}_i}\left(g_i\right)]=[g_i*M_{{\mathrm{n\&Phi;}}_i}^{\mathrm{even}},g_i*M_{{\mathrm{n\&Phi;}}_i}^{\mathrm{odd}}]$

计算幅值：

$A_{{\mathrm{n\&Phi;}}_i}\left(f_i\right)=\sqrt{e_{{\mathrm{n\&Phi;}}_i}{\left(f_i\right)}^2+o_{{\mathrm{n\&Phi;}}_i}{\left(f_i\right)}^2},A_{{\mathrm{n\&Phi;}}_i}\left(g_i\right)=\sqrt{e_{{\mathrm{n\&Phi;}}_i}{\left(g_i\right)}^2+o_{{\mathrm{n\&Phi;}}_i}{\left(g_i\right)}^2}$

幅值与相位偏离函数相乘，分别得到：

$A_{n{\mathrm{\&Phi;}}_i}\left(f_i\right){\mathrm{\&Delta;\&Phi;}}_{{\mathrm{n\&Phi;}}_i}\left(f_i\right)=(e_{n{\mathrm{\&Phi;}}_i}\left(f_i\right)\&CenterDot;{\stackrel{\&OverBar;}{\mathrm{\&Phi;}}}_{{\mathrm{\&Phi;}}_i}^{\mathrm{even}}\left(f_i\right)+o_{{\mathrm{n\&Phi;}}_i}\left(f_i\right)\&CenterDot;{\stackrel{\&OverBar;}{\mathrm{\&Phi;}}}_{{\mathrm{\&Phi;}}_i}^{\mathrm{odd}}\left(f_i\right))-|e_{{\mathrm{n\&Phi;}}_i}\left(f_i\right)\&CenterDot;{\stackrel{\&OverBar;}{\mathrm{\&Phi;}}}_{{\mathrm{\&Phi;}}_i}^{\mathrm{odd}}\left(f_i\right){-o}_{{\mathrm{n\&Phi;}}_i}\left(f_i\right)\&CenterDot;{\stackrel{\&OverBar;}{\mathrm{\&Phi;}}}_{{\mathrm{\&Phi;}}_i}^{\mathrm{even}}\left(f_i\right)|$

$A_{n{\mathrm{\&Phi;}}_i}\left(g_i\right){\mathrm{\&Delta;\&Phi;}}_{{\mathrm{n\&Phi;}}_i}\left(g_i\right)=(e_{n{\mathrm{\&Phi;}}_i}\left(g_i\right)\&CenterDot;{\stackrel{\&OverBar;}{\mathrm{\&Phi;}}}_{{\mathrm{\&Phi;}}_i}^{\mathrm{even}}\left(g_i\right)+o_{{\mathrm{n\&Phi;}}_i}\left(g_i\right)\&CenterDot;{\stackrel{\&OverBar;}{\mathrm{\&Phi;}}}_{{\mathrm{\&Phi;}}_i}^{\mathrm{odd}}\left(g_i\right))-|e_{{\mathrm{n\&Phi;}}_i}\left(g_i\right)\&CenterDot;{\stackrel{\&OverBar;}{\mathrm{\&Phi;}}}_{{\mathrm{\&Phi;}}_i}^{\mathrm{odd}}\left(g_i\right){-o}_{{\mathrm{n\&Phi;}}_i}\left(g_i\right)\&CenterDot;{\stackrel{\&OverBar;}{\mathrm{\&Phi;}}}_{{\mathrm{\&Phi;}}_i}^{\mathrm{even}}\left(g_i\right)|$

其中：

$({\stackrel{\&OverBar;}{\mathrm{\&Phi;}}}_{{\mathrm{\&Phi;}}_i}^{\mathrm{even}}\left(f_i\right),{\stackrel{\&OverBar;}{\mathrm{\&Phi;}}}_{{\mathrm{\&Phi;}}_i}^{\mathrm{odd}}\left(f_i\right))=\frac{{\mathrm{\&Sigma;}}_n{e}_{{\mathrm{n\&Phi;}}_i}\left(f_i\right),{\mathrm{\&Sigma;}}_n{o}_{n{\mathrm{\&Phi;}}_i}\left(f_i\right)}{\sqrt{{\left({\mathrm{\&Sigma;}}_n{e}_{{\mathrm{n\&Phi;}}_i}\left(f_i\right)\right)}^2+{\left({\mathrm{\&Sigma;}}_n{o}_{{\mathrm{n\&Phi;}}_i}\left(f_i\right)\right)}^2}}$

$({\stackrel{\&OverBar;}{\mathrm{\&Phi;}}}_{{\mathrm{\&Phi;}}_i}^{\mathrm{even}}\left(g_i\right),{\stackrel{\&OverBar;}{\mathrm{\&Phi;}}}_{{\mathrm{\&Phi;}}_i}^{\mathrm{odd}}\left(g_i\right))=\frac{{\mathrm{\&Sigma;}}_n{e}_{{\mathrm{n\&Phi;}}_i}\left(g_i\right),{\mathrm{\&Sigma;}}_n{o}_{n{\mathrm{\&Phi;}}_i}\left(g_i\right)}{\sqrt{{\left({\mathrm{\&Sigma;}}_n{e}_{{\mathrm{n\&Phi;}}_i}\left(g_i\right)\right)}^2+{\left({\mathrm{\&Sigma;}}_n{o}_{{\mathrm{n\&Phi;}}_i}\left(g_i\right)\right)}^2}}$

根据公式(1)得到f_i与g_i在方向Φ_i上的主方向相位一致特征PC(f_i)与PC(g_i)。

3.3、将图像块f_i的正向角点集中的每个角点作为待匹配点，在图像块g_i中基于相位相关法和二元二次曲面拟合计算对应的亚像素匹配点，则得到图像块g_i中的正向亚像素匹配角点集；正向角点集和正向亚像素匹配角点集记为正向匹配角点对集Γ_i。

特别是由于蓝、绿、红谱段与近红外谱段在灰度上有很大的差异，导致两幅图像提取的角点不一致，从而在不一致角点的基础上进行匹配得不到正确的匹配点，因此，直接用角点匹配得不到正确的匹配结果。图3为红谱段图像(图3a)和近红外谱段图像(图3b)的基于区域相关的角点匹配结果，图像大小1000×1000个像素，因此对于灰度差异较大的两幅图像基于区域相关的方法很难获得正确的匹配结果。

本发明提出的提取两幅图像中的Harris角点，不是建立基准图像的角点与待配准图像的角点的匹配关系，而是通过相位相关建立基准图像的角点到待配准图像的点、以及从待配准图像的角点到基准图像的点的双向匹配策略，一方面有效避免了两幅灰度差异大的图像中提取的角点不一致，同时充分利用了两幅图像的角点信息。1975年，Kuglin和Hines首次提出了基于傅里叶变换的相位相关算法用于检测两幅图像的平移关系(参见：The phase correlation imagealignment method)，假设f₂(x,y)为图像f₁(x,y)在x和y方向分别平移x₀和y₀后的图像，即：

f₂(x,y)＝f₁(x-x₀,y-y₀)

令F₁(u,v)和F₂(u,v)分别为f₁(x,y)和f₂(x,y)的傅里叶变换，由上式得：

$F_2(u,v){=F}_1(u,v)e^{-j2\mathrm{\&pi;}({\mathrm{ux}}_0+{\mathrm{vy}}_0)}$

图像f₁(x,y)和f₂(x,y)的互功率谱为：

$\frac{F_1^{*}(u,v)F_2(u,v)}{\left|F_1^{*}(u,v)F_2(u,v)\right|}=e^{-j2\mathrm{\&pi;}({\mathrm{ux}}_0+{\mathrm{vy}}_0)}$

式中，F₁ ^*(u,v)为F₁(u,v)的复共轭。互功率谱的傅里叶逆变换为一个二维脉冲函数δ(x-x₀,y-y₀)，相位相关算法就是找到该脉冲函数的位置来确定平移参数。

所述基于相位相关法和二元二次曲面拟合计算对应的亚像素匹配点，主要步骤包括：

(1)当待配准图像是蓝谱段或绿谱段图像时，以待匹配点为中心，在图像块f_i、g_i中分别截取大小为m₁×m₁个像素的块区域；当待配准图像是近红外谱段图像时，以待匹配点为中心，在相位一致特征图PC(f_i)与PC(g_i)中分别截取大小为m₁×m₁个像素的块区域，m₁取128或256；

(2)基于相位相关计算两个块区域的相位相关系数，最大相位相关系数对应的平移量即为两个块区域的平移量，得到像素级匹配点；

(3)以像素级匹配点为中心的局部区域内的相位相关系数拟合二元二次曲面，所述局部区域的大小为3×3个像素，曲面方程为z＝ax²+by²+cxy+dx+ey+f，(x,y)表示点的坐标，z表示点(x,y)对应的相位相关系数，a～f为拟合得到的曲面系数，计算曲面最大值所对应的位置即为亚像素匹配点。

3.4、计算从g_i提取的每个角点与正向亚像素匹配角点集中所有角点的最小距离，如果最小距离小于设定阈值则剔除g_i中的该角点，设定阈值一般取5～10个像素，否则保留g_i中的该角点，所有符合上述条件的角点构成g_i中的反向角点集。

3.5、将g_i中的反向角点集中的每个角点作为待匹配点，在图像块f_i中基于相位相关法和二元二次曲面拟合计算其对应的亚像素匹配点，则得到在f_i中的反向亚像素匹配角点集，反向角点集和反向亚像素匹配角点集构成反向匹配角点对集Γ_i；正反向匹配角点对集表示为Γ_i+Γ'_i。

3.6、将图像块f_i划分若干个大小为K×K个像素的块区域，K一般为200～500范围的正整数，遍历所有块区域，判断块区域内是否有匹配的角点对，即是否有集合Γ_i+Γ'_i中的元素，若在区域内没有匹配的角点对则将该块区域的中心点作为构造点，否则继续遍历各个块区域，直到所有的块区域都处理完毕，则得到的所有构造点构成构造点集。

3.7、将构造点集中的每个点作为待匹配点，在图像块g_i中基于相位相关法和二元二次曲面拟合计算其对应的亚像素匹配点，则得到在g_i中的亚像素匹配点集，构造点集及其亚像素匹配点集构成匹配构造点对集Γ”_i。因此，完整的匹配点集合记为Γ_i+Γ'_i+Γ”_i。

3.8、基于RANSAC法剔除误匹配点对，误匹配点集记为Γ”'_i，则图像块f_i与g_i的匹配点对集Γⁱ _total＝Γ_i+Γ'_i+Γ”_i-Γ”'_i。

RANSAC法剔除误匹配点对的具体步骤如下：

(1)将图像块f_i与g_i都划分为大小为r×r个像素的块，r取值2048～4196的整数；

(2)针对每一块区域，统计最终匹配点对集合Γⁱ _total的匹配点对落入该区域的所有匹配点对；

(3)任取3对匹配点对计算从f_i到g_i的变换参数：

$\left\{\begin{array}{c}{x}^{\&prime;}=a_{1}x+b_{1}y+c_1\\ y^{\&prime;}=a_{2}x+b_{2}y+c_2\end{array}\right.$

(4)对落入该区域的匹配点对集合中寻找所有符合距离小于阈值的点对，将它们作为内点，并记录满足条件的内点数量，如果内点数量大于事先设定的阈值则保留，否则舍弃；

(5)对步骤(3)和(4)重复L次，记录每一次的内点数量，迭代次数L的计算方法为：

$L=\frac{\log (1-\mathrm{\&zeta;})}{\log \left({1-e}^R\right)}$

其中，R为计算约束模型所需的最小匹配点对数量，在这里取3对匹配点对计算变换参数；ζ是期望达到的概率。

(6)选取对应内点数最大的变换参数，对应的内点作为最终的内点，也就是正确的匹配点对，不符合距离小于阈值的点对作为误匹配点对，予以剔除。

3.9、重复执行步骤3.1～3.8，直到f(x,y)和g(x,y)的对图像块都处理完毕，得到所有图像块的匹配点对集，即整幅图像的匹配点对集Γ_total。

如图4所示是本发明方法的红外谱段图像与近红外谱段图像的最终匹配点对示意图，图4a是红外谱段图像的结果，图4b是近红外谱段图像的结果。

4、在整幅图像的匹配点对集Γ_total的基础上，基于Delaunay三角网格划分基准图像和待配准图像，分别计算基准图像和待配准图像中对应三角网格的变换系数，将变换系数分别作用于待配准图像对应的三角网格内部的每个点，得到配准图像。在基准图像和待配准图像中划分的Delaunay三角网格分别如图4a和图4b所示。

5、重复步骤2～4，直到蓝、绿、近红外谱段图像都完成配准，分别得到它们的配准图像。

本发明未详细说明部分属本领域技术人员公知常识。

Claims (3)

1.一种遥感卫星多光谱图像配准方法，其特征在于步骤如下：

1)根据多光谱相机蓝、绿、红、近红外四个谱段探测器间的距离，对获取的四个谱段图像分别进行裁剪，得到粗配准的多光谱图像，具体步骤如下：

多光谱相机焦平面上的蓝谱段探测器到绿、红、近红外谱段探测器的距离分别为d₁、d₂、d₃，则从绿、红、近红外谱段图像的顶部分别裁剪δ₁、δ₂、δ₃个像素，从蓝、绿、红谱段图像的底部分别裁剪δ₃、δ₂、δ₁个像素，其中， τ为多光谱相机的单个探测器尺寸，表示下取整，裁剪后的图像即为粗配准的多光谱图像；

2)将粗配准的红谱段图像作为基准图像f(x,y)，粗配准的蓝、绿、近红外图像中的任意一幅图像作为待配准图像，建立基准图像到待配准图像的直方图映射关系，以建立的直方图映射关系对待配准图像进行映射，得到映射图像g(x,y)，所述f(x,y)与g(x,y)的大小为M×N个像素；M、N为正整数；

3)以大小为m×m个像素的图像块分别划分f(x,y)与g(x,y)，得到的每对图像块表示为f_i与g_i，共得到对图像块，其中，m为正整数，i＝1,2,…,表示上取整；针对每一对图像块f_i与g_i进行如下步骤：

3.1)分别提取图像块f_i与g_i的Harris角点，图像块f_i的所有角点构成f_i的正向角点集；

3.2)当待配准图像为近红外谱段图像时，计算图像块f_i的梯度主方向Φ_i；构造方向为Φ_i的二维Log-Gabor滤波器对f_i与g_i进行相位一致特征检测，分别得到f_i与g_i的主方向相位一致特征PC(f_i)与PC(g_i)；

3.3)将图像块f_i的正向角点集中的每个角点作为待匹配点，在图像块g_i中基于相位相关法和二元二次曲面拟合计算对应的亚像素匹配点，则得到在图像块g_i中的正向亚像素匹配角点集；正向角点集和正向亚像素匹配角点集记为正向匹配角点对集Γ_i；

3.4)计算从g_i提取的每个角点与正向亚像素匹配角点集中所有角点的最小距离，如果最小距离小于设定阈值则剔除g_i中的该角点，否则保留g_i中的该角点，所有符合上述条件的角点构成g_i中的反向角点集；

3.5)将g_i中的反向角点集中的每个角点作为待匹配点，在图像块f_i中基于相位相关法和二元二次曲面拟合计算对应的亚像素匹配点，则得到在f_i中的反向亚像素匹配角点集；反向角点集和反向亚像素匹配角点集构成反向匹配角点对集Γ′_i；正反向匹配角点对集表示为Γ_i+Γ′_i；

3.6)将图像块f_i划分若干个大小为K×K个像素的块区域，K为正整数；遍历所有块区域，判断块区域内是否有匹配的角点对，即是否有集合Γ_i+Γ′_i中的元素，若没有则将该块区域的中心点作为构造点，所有的构造点构成构造点集；

3.7)将构造点集中的每个构造点作为待匹配点，在图像块g_i中基于相位相关法和二元二次曲面拟合计算对应的亚像素匹配点，则得到在g_i中的亚像素匹配点集；构造点集及其亚像素匹配点集构成匹配构造点对集Γ″_i；完整的匹配点集合记为Γ_i+Γ′_i+Γ″_i；

3.8)基于RANSAC法剔除误匹配点对，误匹配点集记为Γ″′_i，则图像块f_i与g_i的匹配点对集 ${\mathrm{\&Gamma;}}_{\mathrm{total}}^i={\mathrm{\&Gamma;}}_i+{{\mathrm{\&Gamma;}}^{\&prime;}}_i+{{\mathrm{\&Gamma;}}^{\&prime;\&prime;}}_i-{{\mathrm{\&Gamma;}}^{\&prime;\&prime;\&prime;}}_i;$

3.9)重复执行步骤3.1)～步骤3.8)，直到f(x,y)和g(x,y)的对图像块都处理完毕，得到所有图像块的匹配点对集，即整幅图像的匹配点对集Γ_total；

4)在整幅图像的匹配点对集Γ_total的基础上，基于Delaunay三角网格划分基准图像和待配准图像，分别计算基准图像和待配准图像中对应三角网格的变换系数，将变换系数分别作用于待配准图像对应的三角网格内部的每个点，得到配准图像；

5)重复步骤2)～步骤4)，直到蓝、绿、近红外谱段图像都完成配准，分别得到它们的配准图像。

2.根据权利要求1所述的一种遥感卫星多光谱图像配准方法，其特征在于：所述步骤3.2)中计算图像块f_i的梯度主方向Φ_i，通过以下步骤实现：

3.2.1)将[0,π)区间划分为n个方向，记为集合n为正整数；

3.2.2)计算图像块f_i中每个像素的梯度幅值mag(x,y)和梯度方向θ(x,y)；

$\mathrm{mag}(x,y)=\sqrt{{\mathrm{\&Delta;x}}^2+{\mathrm{\&Delta;y}}^2},$ θ(x,y)＝tg^-1|Δy/Δx|；

其中，Δx＝f_i(x+1,y)-f_i(x-1,y)，Δy＝f_i(x,y+1)-f_i(x,y-1)；当Δx＝0时，θ(x,y)＝π/2，θ(x,y)的范围为[0,π)；

3.2.3)根据划分的n个方向将步骤3.2.2)中得到的该像素的梯度幅值进行分解，具体步骤如下：

计算d＝θ(x,y)×n/π，则该像素的梯度方向介于第个方向和第个方向之间，mod表示求余运算；将该像素的梯度幅值分解到第个方向和第个方向上，得到在第个方向的梯度幅值为在第个方向的梯度幅值为

3.2.4)对图像块f_i中的所有像素重复步骤3.2.2)、步骤3.2.3)，将所有像素的梯度幅值分解到n个方向上；

3.2.5)分别计算n个方向上的幅值和，最大的幅值和记为MAG_max，对应的方向记为属于与相邻的两个方向表示为 ${\mathrm{\&Phi;}}_L=\left\{\begin{array}{cc}{\mathrm{\&Phi;}}_{{\mathrm{MAG}}_{\max }}-\frac{\mathrm{\&pi;}}{n}& {\mathrm{\&Phi;}}_{{\mathrm{MAG}}_{\max }}\&GreaterEqual;\frac{\mathrm{\&pi;}}{n}\\ \frac{(n-1)\mathrm{\&pi;}}{n}& {\mathrm{\&Phi;}}_{{\mathrm{MAG}}_{\max }}=0\end{array}\right.,{\mathrm{\&Phi;}}_R=\left\{\begin{array}{cc}{\mathrm{\&Phi;}}_{{\mathrm{MAG}}_{m\mathrm{ax}}}+\frac{\mathrm{\&pi;}}{n}& {\mathrm{\&Phi;}}_{{\mathrm{MAG}}_{\max }}<\frac{(n-1)\mathrm{\&pi;}}{n}\\ 0& {\mathrm{\&Phi;}}_{{\mathrm{MAG}}_{\max }}=\frac{(n-1)\mathrm{\&pi;}}{n}\end{array}\right.,$ 对应的幅值和分别为MAG_L、MAG_R，利用MAG_max、MAG_L、MAG_R及其对应的方向Φ_L、Φ_R拟合一元二次曲线，则曲线最大峰值所对应的方向即为图像块f_i的梯度主方向Φ_i。

3.根据权利要求1所述的一种遥感卫星多光谱图像配准方法，其特征在于：所述步骤3.3)、步骤3.5)、步骤3.7)中任一项基于相位相关法和二元二次曲面拟合计算对应的亚像素匹配点，通过以下步骤实现：

(a)当待配准图像是蓝谱段或绿谱段图像时，以待匹配点为中心，在图像块f_i、g_i中分别截取大小为m₁×m₁个像素的块区域；当待配准图像是近红外谱段图像时，以待匹配点为中心，在主方向相位一致特征PC(f_i)与PC(g_i)中分别截取大小为m₁×m₁个像素的块区域，m₁取128或256；

(b)基于相位相关计算两个块区域的相位相关系数，最大相位相关系数对应的平移量即为两个块区域的平移量，得到像素级匹配点；

(c)以像素级匹配点为中心的局部区域内的相位相关系数拟合二元二次曲面，所述局部区域的大小为3×3个像素；计算曲面最大值所对应的位置即为亚像素匹配点。