CN102800099A - 多特征多级别的可见光与高光谱图像高精度配准方法 - Google Patents

Abstract

本发明公开了一种可见光图像与高光谱图像的配准方法，包括对可见光图像进行多尺度分解，形成低分辨率的可见光图像；根据高光谱图像生成高光谱图像显著波段图像；在低分辨率可见光图像和高光谱图像显著波段图像上提取SIFT特征、多尺度角点特征和面点特征，匹配SIFT特征并去除外点，利用匹配的SIFT特征对求取变换模型；在每层可见光图像和高光谱图像显著波段图像上利用上一层配准的变换模型作为该层的初始变换模型，进行基于图像块对的多尺度角点特征、面点特征提取，根据初始变换和多尺度角点特征、面点特征集合利用迭代重加权最小二乘法选择变换类型并求取变换参数；根据所述变换模型对高光谱图像进行变换，得到变换后的高光谱图像。

Description

多特征多级别的可见光与高光谱图像高精度配准方法

技术领域

本发明属于图像处理技术领域，具体涉及遥感图像的配准方法，特别是一种用于可见光与高光谱图像的配准方法。本发明可以广泛适用于航天、航空传感器平台所获取的遥感图像的配准。

背景技术

高光谱图像将物体的光谱频段切分的更细，连续的光谱曲线反映了目标的材质信息，这对于识别伪装目标具有重要的意义。例如，处于同一环境下的真坦克和伪装橡胶坦克，在可见光图像上无法区分；虽然借助可见光图像的温度图像可以区分工作状态的真坦克和伪装坦克，但是如果真坦克也处于休息状态或伪装橡胶坦克被置入了发热物体，则仍无法区分坦克的真伪；而通过高光谱图像的光谱曲线可以正确其材质，从而准确将真坦克和伪装橡胶坦克区分开。因此，综合利用可见光与高光谱图像进行融合对情报获取和分析具有重要的作用。作为可见光与高光谱图像融合的前提，可见光与高光谱图像必须进行自动、高精度配准。可见光与高光谱图像高精度自动配准技术是制约可见光与高光谱图像融合应用的关键技术，具有重要的应用价值。

然而，由于空间分辨率差异很大(如40倍)，同一物体在可见光与高光谱图像上的表现差距很大，这种差距给传统的特征匹配方法带来了很大的挑战。实际上，传统的尺度不变特征只能在空间分辨率差异在不超过4倍时保持尺度不变，空间分辨率差异超过4倍时的图像的缩放都是一个科学难题，直接对空间分辨率差异40倍的可见光图像和高光谱图像配准是不可能的。正是因为可见光与高光谱图像的差异较大，目前可见光与高光谱图像配准基本上仍处于研究阶段，尚未有专门的可见光-高光谱图像配准软件。尽管商业软件如Erdas和Envi都提供了配准功能，但这些配准功能都是为可见光图像设计的。在这些软件上对可见光图像和高光谱图像进行配准，成功的概率很小。目前，对可见光与高光谱图像配准的研究大都沿用可见光图像配准的技术路线；然而，由于空间分辨率差异巨大以及高光谱图像质量较差(噪声、漏缝、扭曲)等因素的影响，目前尚无通用、自动、实用的可见光-高光谱图像配准算法。在实际应用中，往往靠人工标定控制点对，费时费力；而且，在空间分辨率差异40倍的情况下，手工标定的误差很大，配准精度无法保证。在这种情况下，必须另辟思路研究新的配准方法，在保证配准精度的同时，尽可能提供其自动化程度和计算速度。

发明内容

(一)要解决的技术问题

本发明所要解决的技术问题是针对可见光图像与高光谱图像的互补性而提出一种特别适合于可见光与高光谱图像的自动、高效、高精度的配准方法，以推动遥感识别技术的发展。

(二)技术方案

本发明提出一种可见光图像与高光谱图像的配准方法，包括如下步骤：

步骤S1、对可见光图像进行多尺度分解，形成低分辨率的可见光图像；

步骤S2、根据高光谱图像生成高光谱图像显著波段图像；

步骤S3、在低分辨率可见光图像和高光谱图像显著波段图像上提取SIFT特征、多尺度角点特征和面点特征，匹配SIFT特征并去除外点，利用匹配的SIFT特征对求取变换模型；

步骤S4、在每层可见光图像和高光谱图像显著波段图像上利用上一层配准的变换模型作为该层的初始变换模型，进行基于图像块对的多尺度角点特征、面点特征提取，根据初始变换和多尺度角点特征、面点特征集合利用迭代重加权最小二乘法选择变换类型并求取变换参数；

步骤S5、根据所述变换模型对高光谱图像进行变换，得到变换后的高光谱图像。

(三)有益效果

本发明的可见光与高光谱图像高精度自动配准方法基于多特征多级别的配准，综合利用了模式识别、计算机视觉领域中最新的特征提取和特征匹配技术及并行计算技术，有效地解决了可见光与高光谱图像的配准问题，弥补了现有配准软件对可见光与高光谱图像配准方面的缺陷。

本发明的配准方法对可见光与高光谱图像配准具有很好的通用性和实用性，能够大大推动可见光与高光谱图像的广泛应用，具有很好的经济效益。

附图说明

图1是本发明的可见光与高光谱图像高精度自动配准方法的系统框图

图2是本发明的方法的一个具体实施例的原始可见光图像和生成的低分辨率图像的例图

图3是本发明的方法的一个具体实施例的原始高光谱图像与生成的显著波段图像的例图

图4是本发明的方法的一个具体实施例的提取SIFT特征的例图

图5是本发明的方法的一个具体实施例的SIFT特征提取中DOG构建图例

图6是本发明的方法的一个具体实施例的SIFT特征提取中极值点提取图例

图7是本发明的方法的一个具体实施例的SIFT特征提取中主方向检测图例；

图8是本发明的一个具体实施例的为SIFT特征构造描述向量的一个示意图；

图9是本发明的方法的一个具体实施例的粗配准流程；

图10是本发明的方法的一个具体实施例的GDBICP核心流程；

图11是本发明的方法的一个具体实施例的逐级配准流程图。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚明白，以下结合具体实施例，并参照附图，对本发明作进一步的详细说明。

本发明的多特征多级别的可见光与高光谱图像的高精度配准方法既可以硬件方式实现，也可以软件方式实现。例如在个人计算机、工控机及服务器上以软件的形式安装并执行，也可本发明的方法做成嵌入式芯片以硬件的形式来体现。下面参照附图对本发明的具体实施方式进行说明。

在下面的描述中，所述的“图像”特指通过遥感设备获得的遥感图像，并且是已进行了数字化的数字图像。然而，本发明并不限于遥感图像，对于其它技术领域中获得的需要配准的可见光图像和高光谱图像，本发明也可适用。

图1给出了一个实施例系统框图，如图1所示，总的来说，本发明的方法包括如下步骤：

步骤S1：多尺度分解。对可见光图像进行多尺度分解，形成低分辨率的可见光图像；

根据本发明的一种具体实施方式，该实施例中的多尺度分解通过降采样实现。虽然多尺度分解可以由小波金字塔、高斯金字塔等多种方法实现，但计算量比降采样大的多；而在精配准阶段，SIFT(Scale Invariant FeatureTransform，尺度不变特征变换)特征还要在每层可见光图像上提取。因此，基于将采样的多尺度分解方式在保持配准精度的同时，还大大减少了计算量。

步骤S2：根据高光谱图像生成高光谱图像显著波段图像。

根据本发明的一种具体实施方式，将高光谱图像的各波段的均值图像作为显著波段图像。

步骤S3：粗配准。在低分辨率可见光图像和高光谱图像显著波段图像上提取SIFT特征、多尺度角点特征和面点特征，匹配SIFT特征并去除外点，利用匹配的SIFT特征对求取变换模型，变换模型求取包括变换类型选择和变换参数求解两部分；

根据本发明的一种具体实施方式，角点特征分为主动匹配角点特征和被动匹配角点特征，面点特征分为主动匹配面点特征和被动匹配面点特征。主动匹配角(面)点特征是参考图像(如红外图像)上需要匹配的角(面)点特征，被动匹配角(面)点特征是待配准图像(如高光谱图像)上可能匹配的角(面)点特征。主动匹配角(面)点特征是被动匹配角(面)点特征的子集。

根据本发明的一种具体实施方式，采用最近邻比值法来匹配SIFT特征，并去除外点，外点指不满足上述变换模型的SIFT特征对；

根据本发明的一种具体实施方式，利用GDBICP(Generalized DualBootstrap Iterative Closest Point)法得到初始变换模型。

根据本发明的一种具体实施方式，变换模型包含变换类型和变换参数-两层含义。同样，变换模型的确定也包括变换类型选择和变换参数求解。上述变换模型是后续配准变换模型的初始值即初始变换模型，在精配准阶段变换模型的类型及其参数都将不断调整。

步骤S4：精配准。在每层可见光图像和高光谱图像显著波段图像上利用上一层配准的变换模型作为该层的初始变换模型，进行基于图像块对的多尺度角点特征、面点特征提取，根据初始变换和多尺度角点特征、面点特征集合利用迭代重加权最小二乘法选择变换类型并求取变换参数。上述过程迭代进行，在原始可见光图像和高光谱图像显著波段图像上得到的变换模型作为最终的变换模型。

根据本发明的一种具体实施方式，利用迭代重加权最小二乘法选择每一层的变换类型并求取变换参数。

步骤S5：图像变换。根据所述变换模型对高光谱图像进行变换，得到变换后的高光谱图像。

根据本发明的一种具体实施方式，根据变换模型对高光谱图像进行双三次样条插值，然后根据变换模型及原始可见光图像和高光谱图像的尺寸计算参考图像和变换后的高光谱图像的最大重叠区域的坐标，取出、保存所述最大重叠区域的坐标内的可见光图像和高光谱图像。

为使本发明的目的、技术方案和优点更加清楚明白，以下结合具体实施例，对本发明作进一步的详细说明。

在该实施例中，步骤S1利用双线性插值方法对可见光图像进行逐尺度降采样，生成多尺度可见光图像。设原始可见光图像为I，则多尺度可见光图像序列中第k层图像I_k是由第k-1层图像I_k-1经双线性插值方法降采样得到的，降采样的倍数为2，I₀＝I，0≤k＜n。可见光图像的多尺度分解层数n由如下公式确定：n＝log2(N/M)，其中，N为可见光图像长度和宽度的最小像素值个数，M为高光谱图像长度和宽度的最小像素值个数。图2为原始可见光图像和生成的多尺度可见光图像的例图，其中上图为原始可见光图像，下图为低分辨率图像。

显著波段图像也可通过波段选择、波段融合等方式生成，但计算量很大。为了减少计算量并克服噪声的影响，在该实施例中，将高光谱图像的各波段的均值图像作为显著波段图像。图3为原始高光谱图像与生成的显著波段图像的例图，其中上图为原始高光谱图像，下图为显著波段图像。

图4是本发明的方法的具体实施例中步骤S3的提取SIFT特征的流程图。在该实施例中，SIFT特征的提取是按照如下流程实现的：

步骤S31、构建图像的DOG(Difference of Gaussian，高斯差分)金字塔。

设图像为I(x，y)，则DOG金字塔中第k层高斯差分图像为D(x，y，σ)＝L(x，y，kσ)-L(x，y，σ)，其中，L(x，y，σ)＝G(x，y，σ)*I(x，y)，

*表示卷积运算。图5是构建DOG金字塔的示意图，如图5所示，图像在5个尺度上进行高斯卷积，得到的DOG金字塔中有4个图像。

步骤S32、在DOG金字塔的每层上提取极值点，所谓极值点是指在局部邻域内D(x，y，kσ)取值最大的点。

提取极值点的过程为，选择DOG金字塔上的任一点，若该点在该层及上、下相邻两层的26个邻域内不是极值点，则将该点去除，否则将该点作为极值点。图6是提取极值点的示意图，如图6所示，DOG金字塔中的第k层图像的标“×”的点与第k-1层、第k+1层的26个标

的点进行比较，若标“×”的点对应的D(x，y，kσ)是这26个邻域点中的最大值，则标“×”的点为极值点。

步骤S33、对于所提取的极值点，去除局部曲率非常不对称的极值点。在该实施例中，计算差分图像D的局部Hessian矩阵H，去除满足下面条件的极值点：tr(H)²/det(H)＞10，其中det(H)表示矩阵H的行列式值，tr(H)表示矩阵H的迹；

步骤S34、计算SIFT特征的亚像素级别的空间位置、尺度，所述SIFT特征指的是保留下来的极值点。

假设SIFT特征X＝(x，y，σ)，其中x，y，σ分别为极值点X的x、y方向坐标及尺度参数，此时三个方向的坐标都为正数。根据高斯差分图像D(x，y，σ)泰勒展开公式 $D\left(X\right)=D+\frac{{\∂D}^T}{\∂X}X+\frac{1}{2}{X}^T\frac{{\∂}^{2}D}{\∂X^2}X$ 计算出以SIFT特征X为原点相对于X在行、列及层数上的偏移量，即

其中D表示泰勒展开公式的一次项，D、

和

可由SIFT特征X及其邻域的像素根据差分法计算。然后按照如下规则对极值点进行亚像素级插值：如果三个方向上的偏移量都小于0.5像素，则该点就是极值点

为所求的亚像素级极值点坐标；如果某一方向上的偏移量大于等于0.5像素，如假设x方向的偏移量大于0.5像素，将偏移量四舍五入到一整数值a，将a与x相加得到新的SIFT特征X2＝(x+a，y，s)，然后将极值点X2按上述步骤进行操作，直到三个方向上的偏移量都小于0.5。

步骤S35、确定所述SIFT特征主方向，所谓主方向是指与以SIFT特征为中心的邻域内的梯度方向直方图的峰值对应的梯度方向。

在该实施例中，在以SIFT特征X＝(x，y，σ)为中心、以1.5σ为半径的邻域窗口内采样，计算高斯平滑图像L(x，y，σ)在上述邻域窗口内的每一像素的梯度方向并用直方图统计邻域像素的梯度方向，得到一个梯度方向直方图。梯度方向直方图是一种关于梯度方向θ的统计图，其范围是0～360度，其中每10度一个柱，总共36个柱。该梯度方向直方图的峰值代表了该极值点处邻域梯度的主方向，即作为该SIFT特征的方向。图7是本发明的该实施例的梯度直方图的一个示例图。在该图中显示了采用7个柱为SIFT特征确定主方向的示例。在该梯度方向直方图中，当存在另一个相当于主峰值80％能量的峰值时，则将这个方向认为是该SIFT特征的辅方向。一个SIFT特征可能会被指定具有多个方向(一个主方向，一个以上辅方向)。

步骤S36、为SIFT特征构造描述向量，所谓描述向量是指用于刻画该SIFT特征周围图像块统计特征的、由梯度方向直方图构成的向量。

在该实施例中，首先将SIFT特征周围图像块的坐标轴旋转为所述SIFT特征的方向，以确保旋转不变性；然后将SIFT特征周围16×16像素的窗口内分成4×4像素的小块，计算每个小块的8个方向的梯度方向直方图，将每个小块的梯度方向直方图连接起来形成128维的特征向量；最后，将其归一化到单位长度。图8是本发明的为SIFT特征构造描述向量的一个示意图。

本发明的步骤S3还涉及到在多个尺度上提取角点特征和面点特征。对图像I，在

多个尺度上提取角点特征和面点特征。k＝0，…，N，N控制尺度的个数，本实例N＝4。在每个尺度上，角点特征和面点特征提取方法如下：

步骤S31’、对于图像的每个像素点x计算互相关矩阵M(x)。互相关矩阵是用来描述该像素点与其邻域像素的关系。计算方式如下式，其中

表示图像I在水平、竖直方向的梯度图像，w()表示高斯权重函数，其标准差为σ，像素点x邻域N(x)的尺寸为3σ。

$M\left(x\right)={\mathrm{\Σ}}_{y\∈N\left(x\right)}w(x-y)(\▿I\left(y\right)){(\▿I\left(y\right))}^T.$

步骤S32’、根据互相关矩阵的特征值选取角点特征和面点特征，得到角点特征集合和面点特征集合。

首先对互相关矩阵M(x)进行特征值分解，即M(x)＝∑_i＝1，2λ_i(x)Γ_i(x)Γ_i(x)^T，λ₁(x)≤λ₂(x)，λ₁(x)和λ₂(x)表示x处M(x)的特征值，Γ_i(x)表示与λ_i(x)对应的特征向量，Γ_i(x)^T表示Γ_i(x)的转置。

然后，根据特征值来选取角点特征和面点特征。若λ₁(x)/λ₂(x)＞0.1且m(x)＝λ₁(x)+λ₂(x)＞t_m，则x为候选角点特征；若λ₁(x)/λ₂(x)≤0.1且m(x)＝λ₁(x)+λ₂(x)≤t_m，x为候选面点特征。对每个像素点进行上述操作，得到角点特征集合C和面点特征集合F。t_m是一阈值，表示特征的强度。t_m应取一正数，t_m越大得到的候选角点特征和面点特征越少，不同的t_m取值对后续的特征匹配影响不大，在该实施例中，t_m＝1。

步骤S33’、去除角点特征集合和面点特征集合中的噪声点。对于角点特征集合C或面点特征集合F中的某特征点x，若m(x)＜μ_m+0.5σ_m，将x从C或F中去除。其中，μ_m和σ_m分别为半径r＝30的邻域内的m(x)的均值和方差。

步骤S34’、对角点特征和面点特征进行非极值点抑制。对于角点特征集合C中的角点特征x，沿特征值λ₁(x)和λ₂(x)对应的特征向量Γ₁(x)和Γ₂(x)的方向进行非极大值点抑制；对于面点特征集合F中的某面点特征x，沿λ₂(x)对应的特征向量的方向Γ₂(x)进行非极大值点抑制。具体方法如下：对以x为中心、沿Γ_i(x)方向的半径r＝30的邻域内，若存在x′满足m(x′)＞m(x)，则x不是好的特征点，将其从候选特征点集合中去掉。

步骤S35’、根据特征强度的大小来选择被动角点特征和被动面点特征。

在该实施例中，对角点特征集合C或面点特征集合F中的特征按特征强度m(x)从大到小排序，选择前n个角点特征和面点特征作为被动匹配角点特征和面点特征的集合。n根据图像大小确定，对于1024×1024像素的图像，n一般取1000。

步骤S36’、根据特征强度的大小来选择主动角点特征和主动面点特征。

改变特征强度阈值，令t_m＝2t_m，r＝2r，重复步骤S32’-S35’。对角点特征集合C或面点特征集合F中的特征点按特征强度m(x)从大到小排序，选择前m个角点特征和m个面点特征作为主动匹配角点特征、面点特征集合。m根据图像大小确定，对于1024×1024像素的图像，m一般取500。同一图像的主动匹配角(面)点特征是被动匹配角(面)点特征的子集，在匹配时，参考图像的主动匹配角(面)点特征与待配准图像的被动匹配角(面)点特征(而不是待配准图像的主动匹配角(面)点特征)进行匹配，这样做的好处是可以得到更稠密的角(面)点特征匹配点对，从而有利用提高配准精度。

按照上述步骤在每个尺度σ上提取被动匹配角点特征、面点特征和主动匹配角点特征、面点特征，将各尺度、各类型特征点合并在一起，组成被动匹配角点特征集合C_m、被动匹配面点特征集合F_m、主动匹配角点特征集合C_d、主动匹配面点特征集合F_d。每一个角点特征c_i由三维向量组成，即c_i＝(x_i，y_i，s_i)，s_i等于它所在的尺度参数σ_k。每一个面点特征f_i由四维向量组成，即f_i＝(x_i，y_i，s_i，η_i)，η_i表示方向，η_i对应步骤S26’中在(x_i，y_i)处的Γ₂(x)。

图9是本发明的方法的步骤S3的流程图。该步骤为粗配准步骤。如图9所示，本发明中的步骤S3是按如下流程实现的：

选择置信度高且未尝试过的SIFT特征对，由该SIFT特征对及被动匹配角点特征集合C_m、被动匹配面点特征集合F_m、主动匹配角点特征集合C_d、主动匹配面点特征集合F_d通过GDBICP方法进行匹配，若该次匹配成功，粗配准终止；否则，选择下一个置信度高的SIFT特征对通过GDBICP方法进行匹配。

图10是本发明的方法的步骤S3中的GDBICP方法的流程图。如图10所示，GDBICP方法包括如下步骤：

步骤S31”、根据选择的SIFT特征对确定初始变换模型(相似变换)，SIFT特征对所在的邻域作为初始匹配区域；

设当前选择的SIFT特征对为A＝(x₁，y₁，s₁，θ₁)和B＝(x₂，y₂，s₂，θ₂)，其中x_i和y_i是坐标值，s_i和θ₂是对应的尺度参数和主方向参数。该步骤为粗配准，变换类型为相似变换，SIFT特征对之间应满足 $\left[\begin{array}{c}{x}_2\\ y_2\end{array}\right]=\left[\begin{array}{ccc}s\·\cos \mathrm{\θ}& s\·\sin \mathrm{\θ}& t_x\\ -\sin \mathrm{\θ}& s\·\cos \mathrm{\θ}& t_y\end{array}\right]\left[\begin{array}{c}{x}_1\\ y_1\\ 1\end{array}\right],$ s、θ、t_x、t_y为所求的变换参数。根据A和B的尺度参数和主方向参数的对应关系，可求得

θ＝θ₂-θ₁，x_t＝x₂-s·cosθ·x₁-s·sinθ·y₁，y_t＝y₂+s·sinθ·x₁-s·cosθ·y₁。以(x_i，y_i)为中心、以h+3·s_i为半径的区域分别作为A和B的初始匹配区域。

步骤S32”、根据当前迭代步骤的变换模型(注：求变换模型是一个迭代的过程，S4的变换模型是最终的变换模型，当前迭代步骤的变换模型是一个中间结果)和当前匹配区域，分别从前向和后向利用ICP(IterativeClosest Point)方法匹配当前匹配区域内的角点特征和面点特征；

前向匹配以可见光图像为参考图像，为高光谱图像上的每一个角(面)点特征寻找其在可见光图像上对应的角(面)点特征。假设当前变换模型为

(k表示变换类型，k＝0表示相似变换、k＝1表示仿射变换、k＝2表示投影变换)，对于高光谱图像的任一角点特征c_i＝(x_i，y_i，s_i)，它对应可见光图像的坐标记为在可见光图像上距离

最近的三个角点特征为i＝1，2，3.然后计算

中每个元素与c_i的匹配度

将与最大的

(i＝1，2，3)对应的

作为c_i的匹配角点特征，记为e_i。类似的，为高光谱图像的每一面点特征f_i＝(x_i，y_i，s_i，η_i)寻找可见光图像对应的匹配面点特征h_i，面点特征的匹配度计算公式为 $w_i^f=\min (s_i^f/s_i,s_i/s_i^f)\·|{\mathrm{\η}}_i^T\·{\mathrm{\η}}_i^f|,$ $g_i^f=(x_i^g,y_i^g,s_i^g,{\mathrm{\η}}_i^g)$ (i＝1，2，3)为可见光图像上距离

最近的三个面点特征。

后向匹配以高光谱图像为参考图像，为可见光图像上的每一个角点特征和面点特征寻找其在高光谱图像上对应的角点特征和面点特征，操作过程与前向匹配类似。

步骤S33”、根据当前匹配的角点特征和面点特征利用迭代重加权最小二乘法重新估计前向和后向变换参数；通过ICP和参数估计的交替迭代，确定最优变换参数和最优的匹配；

当前变换模型与当前匹配的角点特征和面点特征紧密相关。为此，在给定变换类型的前提下采用迭代重加权最小二乘法不断调整变换参数及匹配的角点特征和面点特征集合。为方便叙述，在步骤S33中，当变换类型给定后将变换模型

记为θ^pq。迭代重加权最小二乘法的目标函数为：

$E({\mathrm{\θ}}^{\mathrm{pq}};C_c^{\mathrm{pq}},C_f^{\mathrm{pq}})=\underset{(p_i,q_i)\∈C_c^{\mathrm{pq}}}{\mathrm{\Σ}}{w}_{c;i}\mathrm{\ρ}(d_c(p_i,q_i;{\mathrm{\θ}}^{\mathrm{pq}})/{\mathrm{\σ}}_c)+\underset{(p_i,q_i)\∈C_f^{\mathrm{pq}}}{\mathrm{\Σ}}{w}_{f;i}\mathrm{\ρ}(d_f(p_i,q_i;{\mathrm{\θ}}_{\mathrm{pq}})/{\mathrm{\σ}}_f)$

上式中， $\mathrm{\ρ}\left(u\right)=\left\{\begin{array}{cc}\frac{a^2}{6}[1-(1-{\left(\frac{u^2}{a}\right)}^3]& \left|u\right|\≤a\\ \frac{a^2}{6}& \left|u\right|>a\end{array}\right.,$ a＝4，d_c(p，q；θ^pq)＝||T(p；θ^pq)-q||/s_q

d_f(p，q；θ^pq)＝|(T(p；θ^pq)-q)^Tη_q|/s_q，w_c；i和w_f；i为每一对角点特征和面点特征的权重。为简化计算，实际求解时将上述目标函数拆分为两个步骤：求解每一对角(面)点对匹配权重w_d，i和更新变换参数-，这两个步骤循环交叉迭代，直到算法收敛。匹配权重w_d，i根据当前变换模型θ^pq下匹配点对(p_i，q_i)的匹配误差进行计算，对于角点，

对于面点， $w_{d,i}=w(d_f(p_i,q_i;{\mathrm{\θ}}^{\mathrm{pq}})/{\mathrm{\σ}}_f)/{\mathrm{\σ}}_f^2,$ 其中， $w\left(u\right)=\left\{\begin{array}{cc}1-{\left(\frac{u}{a}\right)}^2& \left|u\right|\≤a\\ 0& \left|u\right|>a\end{array}\right..$ 在当前匹配点对权重w_d，i的基础上，利用加权最小二乘法估计新通过优化 $F({\mathrm{\θ}}^{\mathrm{pq}};C_c^{\mathrm{pq}},C_f^{\mathrm{pq}})=\underset{(p_i,q_i)\∈C_c^{\mathrm{pq}}}{\mathrm{\Σ}}{w}_{d;i}{w}_{c;i}{d}_c^2(p_i,q_i;{\mathrm{\θ}}^{\mathrm{pq}})+\mathrm{\ρ}\underset{(p_i,q_i)\∈C_f^{\mathrm{pq}}}{\mathrm{\Σ}}{w}_{d;i}{w}_{f;i}{d}_f^2(p_i,q_i;{\mathrm{\ρ}}^{\mathrm{pq}})$ 求解新的变换参数。σ_c和σ_f为角点特征和面点特征的尺度因子，计算公式如下：

${\mathrm{\σ}}_c=\sqrt{\underset{(p_i,q_i)\∈C_c^{\mathrm{pq}}}{\mathrm{\Σ}}{w}_{d;i}{w}_{c;i}{d}_c^2(p_i,q_i;{\mathrm{\θ}}^{\mathrm{pq}})/\underset{(p_i,q_i)\∈C_c^{\mathrm{pq}}}{\mathrm{\Σ}}{w}_{d;i}{w}_{c;i}},$

${\mathrm{\σ}}_f=\sqrt{\underset{(p_i,q_i)\∈C_f^{\mathrm{pq}}}{\mathrm{\Σ}}{w}_{d;i}{w}_{f;i}{d}_f^2(p_i,q_i;{\mathrm{\θ}}^{\mathrm{pq}})/\underset{(p_i,q_i)\∈C_f^{\mathrm{pq}}}{\mathrm{\Σ}}{w}_{d;i}{w}_{f;i}}.$

步骤S34”、利用当前最优匹配的角点特征和面点特征集合，根据Akaike Information准则进行变换类型选择；

令变换类型k＝0，1，2，并分别进行步骤S33”操作，得到对应的前向变换模型和后向变换分别为

和

当前最优匹配的前向匹配角点特征和面点特征集合为和

后向匹配角点特征和面点特征集合为

和

在Akaike Information准则下的不适度为：

$I_i=-2\left[\right|C_c^{\mathrm{pq}}|\log \left({\mathrm{\σ}}_c^{\mathrm{pq}}\right)+|C_f^{\mathrm{pq}}|\log \left({\mathrm{\σ}}_f^{\mathrm{pq}}\right)+E({\mathrm{\θ}}^{\mathrm{pq}};C_c^{\mathrm{pq}},C_f^{\mathrm{pq}})]$

$-2\left[\right|C_c^{\mathrm{qp}}|\log \left({\mathrm{\σ}}_c^{\mathrm{qp}}\right)+|C_f^{\mathrm{qp}}|\log \left({\mathrm{\σ}}_f^{\mathrm{qp}}\right)+E({\mathrm{\θ}}^{\mathrm{qp}};C_c^{\mathrm{qp}},C_f^{\mathrm{qp}})]+2\mathrm{nl}/(n-l-1).$ 上式中，l表示变换类型k的自由度，k＝0时，l＝4；k＝1时，l＝6；k＝2时，l＝9；|A|表示集合A中元素个数。与最小的I_i对应的

和

为当前匹配区域的最优前向、后向变换模型。

步骤S35”、根据当前变换模型的不确定性(用变换参数的方差矩阵刻画)进行匹配区域增长，区域增长的速度与当前匹配区域的前向匹配和后向匹配的角点特征和面点特征的匹配误差成反比。当匹配区域增长到角点特征和面点特征的最大外接矩形时，GDBICP收敛，此时的前向匹配变换模型

作为步骤S4的初始变换模型。

本发明的该实施例的步骤S4是按如下流程实现的，该步骤S4是精配准步骤：

S4在第i层可见光图像I_i(n-2≤i≤0)和高光谱图像显著波段图像上利用第i层可见光图像和高光谱图像显著波段图像的变换模型

作为初始变换模型，然后逐级求取第i层的变换模型。上述过程是一个循环的迭代过程，在第0层的变换模型即为最终的变换模型。图11是本发明的方法的步骤S4中的逐级配准的流程图，逐级配准是按如下流程实现的。

步骤S41：根据初始变换模型

进行基于图像块对的多尺度角点特征、面点特征特征提取；

首先将高光谱图像进行均匀分块，每块大小为256×256像素。假设某子块C1的左上角、右下角坐标分别为lu＝(x0，y0)，rd＝(x1，y1)，根据变换模型

可得到对应的可见光图像I_i的子块C2，在C1和C2上按照步骤S31”至S36”的方法提取多尺度角点特征和面点特征。之所以将图像分块，是为了对大尺寸图像进行并行计算，即并行控制端将图像块对发送到多个并行终端，每个并行终端可以为CPU或计算机，在每个并行终端分块特征提取完毕后，只需将特征发送到并行控制端并对特征的坐标值进行修正后即可得到特征的全局坐标。例如，在子块C1上提取的某角点特征或面点特征k的坐标为(x，y)，则修正后的k的坐标为(x0+x，y0+y)。将在原始可见光图像和高光谱图像显著波段图像上按照上述方法可以得到主动匹配角点特征、主动匹配面点特征、被动匹配角点特征、被动匹配面点特征。

步骤S42：根据初始变换模型

和多尺度角点特征、面点特征集合利用迭代重加权最小二乘法选择变换类型并求取变换参数。

角点特征和面点特征的最大外接矩形作为匹配区域，利用初始变换模型

和步骤S33”至步骤S34”的方法选择变换类型并求取相应的变换参数。

本发明的该实施例的步骤S5是按如下流程实现的：根据步骤S4求得的变换模型，按照高光谱图像的波段排列顺序对高光谱图像的每一波段图像分别进行双三次样条插值得到变换后的高光谱图像，变换后的高光谱图像与可见光图像已经配准。

以上所述的具体实施例，对本发明的目的、技术方案和有益效果进行了进一步详细说明，应理解的是，以上所述仅为本发明的具体实施例而已，并不用于限制本发明，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (10)

1.一种可见光图像与高光谱图像的配准方法，其特征在于，包括如下步骤：

步骤S1、对可见光图像进行多尺度分解，形成低分辨率的可见光图像；

步骤S2、根据高光谱图像生成高光谱图像显著波段图像；

步骤S3、在低分辨率可见光图像和高光谱图像显著波段图像上提取SIFT特征、多尺度角点特征和面点特征，匹配SIFT特征并去除外点，利用匹配的SIFT特征对求取变换模型；

步骤S4、在每层可见光图像和高光谱图像显著波段图像上利用上一层配准的变换模型作为该层的初始变换模型，进行基于图像块对的多尺度角点特征、面点特征提取，根据初始变换和多尺度角点特征、面点特征集合利用迭代重加权最小二乘法选择变换类型并求取变换参数；

步骤S5、根据所述变换模型对高光谱图像进行变换，得到变换后的高光谱图像。

2.如权利要求1所述的可见光图像与高光谱图像的配准方法，其特征在于，在步骤S1中，对可见光图像进行多尺度分解，形成低分辨率的可见光图像。

3.如权利要求1所述的可见光图像与高光谱图像的配准方法，其特征在于，在步骤S2中，将高光谱图像的各波段的均值图像作为显著波段图像。

4.如权利要求1所述的可见光图像与高光谱图像的配准方法，其特征在于，在步骤S3中，采用最近邻比值法来匹配SIFT特征，并去除外点，外点指不满足上述变换模型的SIFT特征对。

5.如权利要求1所述的可见光图像与高光谱图像的配准方法，其特征在于，在步骤S3中，利用GDBICP法得到初始变换模型。

6.如权利要求1所述的可见光图像与高光谱图像的配准方法，其特征在于，在步骤S3中，提取SIFT特征的步骤包括：

步骤S31、构建图像的DOG金字塔；

步骤S32、在DOG金字塔的每层上提取极值点，所谓极值点是指在局部邻域内取值最大的点；

步骤33、对于所提取的极值点，去除局部曲率不对称的极值点；

步骤S34、计算SIFT特征的亚像素级别的空间位置、尺度，所述SIFT特征指的是保留下来的极值点；

步骤S35、确定所述SIFT特征主方向，所谓主方向是指与以SIFT特征为中心的邻域内的梯度方向直方图的峰值对应的梯度方向；

步骤S36、为SIFT特征构造描述向量，所谓描述向量是指用于刻画该SIFT特征周围图像块统计特征的、由梯度方向直方图构成的向量。

7.如权利要求1所述的可见光图像与高光谱图像的配准方法，其特征在于，在步骤S3中，提取角点特征和面点特征的步骤包括：

步骤S31’、对于图像的每个像素点计算互相关矩阵，所述互相关矩阵是用来描述该像素点与其邻域像素的关系；

步骤S32’、根据互相关矩阵的特征值选取角点特征和面点特征，得到角点特征集合和面点特征集合；

步骤S33’、去除角点特征集合和面点特征集合中的噪声点；

步骤S34’、对角点特征和面点特征进行非极值点抑制；

步骤S35’、根据特征强度的大小来选择被动角点特征和被动面点特征；

步骤S36’、根据特征强度的大小来选择主动角点特征和主动面点特征。

8.如权利要求1所述的可见光图像与高光谱图像的配准方法，其特征在于，在步骤S3中，选择置信度高且未尝试过的SIFT特征对，由该SIFT特征对及被动匹配角点特征集合Cm、被动匹配面点特征集合Fm、主动匹配角点特征集合Cd、主动匹配面点特征集合Fd通过GDBICP方法进行匹配，若该次匹配成功，粗配准终止；否则，选择下一个置信度高的SIFT特征对通过GDBICP方法进行匹配。

9.如权利要求8所述的可见光图像与高光谱图像的配准方法，其特征在于，所述GDBICP方法包括如下步骤：

步骤S31”、根据选择的SIFT特征对确定初始变换模型(相似变换)，SIFT特征对所在的邻域作为初始匹配区域；

步骤S32”、根据当前迭代步骤的变换模型和当前匹配区域，分别从前向和后向利用ICP方法匹配当前匹配区域内的角点特征和面点特征；

步骤S33”、根据当前匹配的角点特征和面点特征利用迭代重加权最小二乘法重新估计前向和后向变换参数；通过ICP和参数估计的交替迭代，确定最优变换参数和最优的匹配；

步骤S34”、利用当前最优匹配的角点特征和面点特征集合，根据Akaike Information准则进行变换类型选择；

步骤S35”、根据当前变换模型的不确定性(进行匹配区域增长，区域增长的速度与当前匹配区域的前向匹配和后向匹配的角点特征和面点特征的匹配误差成反比。当匹配区域增长到角点特征和面点特征的最大外接矩形时，GDBICP收敛，此时的前向匹配变换模型

作为步骤S4的初始变换模型。

10.如权利要求1所述的可见光图像与高光谱图像的配准方法，其特征在于，所述步骤S4包括：

步骤S41：根据初始变换模型

进行基于图像块对的多尺度角点特征、面点特征特征提取；

步骤S42：根据初始变换模型

和多尺度角点特征、面点特征集合利用迭代重加权最小二乘法选择变换类型并求取变换参数。

Images