CN105303545A - 一种多光谱影像和全色卫星影像融合方法 - Google Patents

Abstract

本发明公开了一种多光谱影像和全色卫星影像融合方法，其包括：分别提取多光谱影像与全色卫星影像的SIFT特征点，并做降维处理；采用BBF算法对降维处理后的SIFT特征点进行粗匹配，并进行误匹配消除，再采用仿射变换进行精确匹配；对精确匹配后的多光谱影像进行HIS变换获得亮度、饱和度与色度分量，将其中的亮度分量以及全色卫星影像的直方图分别进行小波变换，并将二者的小波变化结果依次进行高频替换、小波逆变换得到新的亮度分量；通过预设的参数对所述新的亮度分量进行更新，再将更新后的亮度分量与所述饱和度与色度分量进行HIS逆变换，从而实现影像的融合。本方法可以对全色卫星影像和多光谱影像进行准确、有效配准，并有效获取融合图像。

Description

一种多光谱影像和全色卫星影像融合方法

技术领域

本发明涉及卫星遥感图像处理技术领域，尤其涉及一种多光谱影像和全色卫星影像融合方法。

背景技术

目前，在遥感影像的融合技术中，多光谱影像MS和高空间分辨率全色影像PAN的融合是最重要的一类。通过多光谱影像和全色影像融合后，获得了同时具有高空间分辨率和光谱质量的遥感图像是当前研究热点。过去常用的一些遥感影像处理方法，如HIS(亮度、饱和度、色度)变换和PCA(主成分分析)变换融合方法，其融合后影像，在保持较高空间分辨率的同时，总是发生严重的光谱失真现象。如何在提高空间分辨率的同时，能够保持融合影像的光谱质量，这依旧是当前遥感影像融合研究中的难点。

国内现有的遥感融合技术主要包括代数运算法、高通滤波法、HIS变换法、主分量分析法、分量替换法以及较热门的小波分析融合法等等。影像融合方法从传统的HIS变换、PCA变换等方法，经过金字塔分解融合、小波变换融合等方法的洗礼，逐步进入到新的应用Contourlet变换融合阶段。与此同时，在融合算法的改造上，越来越多的传统算法相结合的方式被使用，在一定程度上增进了融合效果，提升了融合效率。对于HIS变换融合和PCA变换融合等方法，他们对于多光谱及全色影像的融合结果都能有不错的空间几何信息保有度，能保持较高的空间分辨率，但总是会发生较大的光谱扭曲。小波变换影像融合尽管同传统的像素级影像融合方法相比，具有很好的光谱质量，但其对空间细节信息的保持则显得不足，融合影像常常出现模糊甚至“重影”现象。

鉴于此，有必要研究一种可以对多光谱影像和全色卫星影像进行有效融合的方案。

发明内容

本发明的目的是提供一种多光谱影像和全色卫星影像融合方法，可以对全色卫星影像和多光谱影像进行准确、有效配准，并有效获取融合图像。

本发明的目的是通过以下技术方案实现的：

一种多光谱影像和全色卫星影像融合方法，包括：

分别提取多光谱影像与全色卫星影像的SIFT特征点，并对提取到的SIFT特征点做降维处理；

采用BBF算法对降维处理后的SIFT特征点进行粗匹配，并对粗匹配后的结果进行误匹配消除，再采用仿射变换进行精确匹配；

对精确匹配后的多光谱影像进行HIS变换获得亮度、饱和度与色度分量，将其中的亮度分量以及精确匹配后的全色卫星影像的直方图分别进行小波变换，并将二者的小波变化结果依次进行高频替换、小波逆变换得到新的亮度分量；

通过预设的参数对所述新的亮度分量进行更新，再将更新后的亮度分量与所述饱和度与色度分量进行HIS逆变换，从而实现多光谱影像和全色卫星影像的融合。

所述分别提取多光谱影像与全色卫星影像的SIFT特征点，并对提取到的SIFT特征点做降维处理包括：

分别构建光谱影像与全色卫星影像的SIFT空间金字塔，利用SIFT空间金字塔求得极值点，并对极值点进行筛选，以去除不稳定极值点，再根据筛选后的极值点进行方向和尺度的计算，从而获得多光谱影像的SIFT特征点以及全色卫星影像的SIFT特征点；

采用SIFT特征描述与PCA算法相结合的方式分别对多光谱影像的SIFT特征点以及全色卫星影像的SIFT特征点做降维处理。

所述采用BBF算法对降维处理后的SIFT特征点进行粗匹配，并对粗匹配后的结果进行误匹配消除后，再采用仿射变换进行精确匹配包括：

对降维处理后的多光谱影像的SIFT特征点以及全色卫星影像的SIFT特征点进行基于堆的优先队列方法搜索匹配；

然后，通过RANSAC算法，对匹配结果进行消除误匹配；

再基于仿射变换模型，对误匹配消除后的多光谱影像的SIFT特征点以及全色卫星影像的SIFT特征点进行平移、缩放和旋转处理，从而完成影像的精确匹配。

所述对精确匹配后的多光谱影像进行HIS变换获得亮度、饱和度与色度分量，将其中的亮度分量以及精确匹配后的全色卫星影像的直方图分别进行小波变换，并将二者的小波变化结果依次进行高频替换、小波逆变换得到新的亮度分量包括：

将多光谱影像从RGB颜色空间转换到HIS颜色空间，获得亮度、饱和度与色度分量，并对全色卫星影像进行直方图匹配；

分别对多光谱影像的亮度分量及全色卫星影像匹配后的直方图进行小波变换，并将亮度分量小波变换后的高频部分移去，用全色卫星影像匹配后的直方图小波变换后的高频部分进行替换，再对替换结果进行小波逆变换，从而得到新的亮度分量。

所述通过预设的参数对所述新的亮度分量进行更新，再将更新后的亮度分量与所述饱和度与色度分量进行HIS逆变换，从而实现多光谱影像和全色卫星影像的融合包括：

利用预设的参数a、b、c与新的亮度分量进行下述计算，从而获得更新后的亮度分量：

I″＝a·I′+b·I_PAN+c·I_MS

其中，参数a、b、c均为常数；I′为所述新的亮度分量，I_PAN为全色卫星影像的亮度分量，I_MS为多光谱影像通过HIS变换获得亮度分量，I″为更新后的亮度分量；

再将更新后的亮度分量I″与多光谱影像通过HIS变换获得的饱和度与色度分量HIS逆变换，变换后获得的融合图像即为多光谱影像和全色卫星影像的融合图像。

由上述本发明提供的技术方案可以看出，利用SIFT与PCA-SIFT结合的方法，找到全色图像与多光谱图像的特征点，基于这些特征点依次使用BBF、RANSAC和仿射空间变换进行精确匹配；然后，根据HIS变换和小波变换的改进方法进行融合；再引入预设的参数，结合全色图像的亮度和HIS变换得到的亮度得到新的亮度分量，最后通过HIS逆变换得到融合图像；本方案在一定程度上解决了目前融合方法所存在的部分问题，使融合后图像质量有明显改进。

附图说明

为了更清楚地说明本发明实施例的技术方案，下面将对实施例描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域的普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他附图。

图1为本发明实施例提供的一种多光谱影像和全色卫星影像融合方法的流程图；

图2为本发明实施例提供的详细的多光谱影像和全色卫星影像融合方法的流程图；

图3为本发明实施例提供的二维图像f(x,y)分解示意图；

图4为本发明实施例提供的图像的Mallat小波分解过程示意图；

图5为本发明实施例提供的图像的重构过程示意图；

图6为本发明实施例提供的对遥感影使用小波变换进行融合的原理示意图。

具体实施方式

下面结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明的保护范围。

图1为本发明实施例提供的一种多光谱影像和全色卫星影像融合方法的流程图。如图1所示，其主要包括如下步骤：

步骤11、分别提取多光谱影像与全色卫星影像的SIFT(尺度不变特征变换)特征点，并对提取到的SIFT特征点做降维处理。

步骤12、采用BBF算法(最优节点优先算法)对降维处理后的SIFT特征点进行粗匹配，并对粗匹配后的结果进行误匹配消除，再采用仿射变换进行精确匹配。

步骤13、对精确匹配后的多光谱影像进行HIS变换获得亮度、饱和度与色度分量，将其中的亮度分量以及精确匹配后的全色卫星影像的直方图分别进行小波变换，并将二者的小波变化结果依次进行高频替换、小波逆变换得到新的亮度分量。

步骤14、通过预设的参数对所述新的亮度分量进行更新，再将更新后的亮度分量与所述饱和度与色度分量进行HIS逆变换，从而实现多光谱影像和全色卫星影像的融合。

本发明实施例的上述方案可以对全色遥感图像和多光谱图像进行准确、有效配准，并有效获取融合图像。

为了便于理解，下面结合附图2对上述过程做进一步的介绍。

如图2所示，为本发明实施例提供的详细的多光谱影像和全色卫星影像融合方法的流程图，其主要包括如下步骤：

1、分别提取多光谱影像与全色卫星影像的SIFT(尺度不变特征变换)特征点，并做降维处理。

其具体过程如下：

1)分别构建光谱影像与全色卫星影像的SIFT空间金字塔，利用SIFT空间金字塔求得极值点，并对极值点进行筛选，以去除不稳定极值点，再根据筛选后的极值点进行方向和尺度的计算，从而获得多光谱影像的SIFT特征点以及全色卫星影像的SIFT特征点；

示例性的，对于影像I(x,y)的尺度空间函数S(x,y,σ)，可由该影像I(x,y)与高斯函数G(x,y,σ)通过卷积运算求得，其中x_i,y_i为空间坐标：

$G(x_i,y_i,\mathrm{\σ})=\frac{1}{2{\mathrm{\π\σ}}^2}\exp (-\frac{{(x-x_i)}^2+{(y-y_i)}^2}{2{\mathrm{\σ}}^2});$

S(x,y,σ)＝G(x,y,σ)·I(x,y)；

将高斯函数进行差分运算，运算结果同图像卷积，从而得到空间函数D(x,y,σ)，并在其中寻找极值点；公式如下：

D(x,y,σ)＝(G(x,y,kσ)-G(x,y,σ))·I(x,y)＝S(x,y,kσ)-S(x,y,σ)

其中k为SIFT空间金字塔层数相关系数，为寻找尺度空间中的特征点，需求得影像的尺度L；通过高斯的差分运算方式求其特征点，则仅需计算相邻的尺度函数之差，这一过程大大减少了计算量。同时，空间函数D的表达形式与性质同尺度归一化的Laplace高斯函数极为近似，该函数表达形式如下：

${\mathrm{\σ}}^2{\▿}^{2}G=\frac{\∂G}{\∂\mathrm{\σ}}\≈\frac{G(x,y,k\mathrm{\σ})-G(x,y,\mathrm{\σ})}{k\mathrm{\σ}-\mathrm{\σ}};$

$G(x,y,k\mathrm{\σ})-G(x,y,\mathrm{\σ})\≈(k-1){\mathrm{\σ}}^2{\▿}^{2}G;$

然后，利用SIFT尺度空间金字塔，可以将原始图像生成高斯差分图像，通过对高斯差分图像的像素点遍历可以求出极值点，由于建立尺度空间过程中，金字塔处理存在降采样，高斯差分图像提取出的极值点在原始图像中的位置会发生变化。当我们需求出某个任意尺度上的极值点的原始位置，可以通过一个3维的2次函数实现；该函数公式如下：

$D\left(\hat{X}\right)=D+\frac{1}{2}\frac{\∂D^T}{\∂X}\hat{X}$

对于选定的极值点，还需进行筛选，以去除不稳定极值点，增强特征点匹配精度，提升其抗噪声能力。对于不稳定极值点进行主曲率壁纸，需要采用近似的方法，先计算待测极值点的海森矩阵H。

$H=\left[\begin{array}{cc}{D}_{xx}& D_{xy}\\ D_{yx}& D_{yy}\end{array}\right]$

计算出H的较大特征值与较小特征值之比即为所求。设α是较大的特征值，β是较小的特征值，由矩阵性质知，矩阵的迹Tr(H)与行列式Det(H)求解如下：

Tr(H)＝D_xx+D_yy＝α+β；

Det(H)＝D_xxD_yy-(D_xy)²＝αβ；

在上式结果中不会为负。若出现负值，即两个主曲率不同号的情况，则视其为非极值点，并丢弃这个点。设r＝α/β，即r为D的主曲率比值。此时计算如下除法运算。

$\frac{Tr{\left(H\right)}^2}{Det\left(H\right)}=\frac{{(\mathrm{\α}+\mathrm{\β})}^2}{\mathrm{\α}\mathrm{\β}}=\frac{{(r\mathrm{\β}+\mathrm{\β})}^2}{{\mathrm{r\β}}^2}=\frac{{(r+1)}^2}{r};$

选取阈值r＝10时，可利用该主曲率比值去除边缘的点。

进一步进行特征点方向和尺度计算，根据SIFT旋转不变的特性，选取关键点赋上方向，定义关键点描述子服从此方向即实现匹配时图像的旋转不变。同时为实现与尺度无关，首先选择与该尺度最接近的高斯平滑图像L，对于L上的每个点L(x,y)，计算其梯度(gradI(x,y))和方向(m(x,y)和θ(x,y))：

$gradI(x,y)=(\frac{\∂I}{\∂x},\frac{\∂I}{\∂y});$

$m(x,y)=\sqrt{{\left(L\right(x+1,y)-L(x-1,y\left)\right)}^2+{\left(L\right(x,y+1)-L(x,y-1\left)\right)}^2};$

$\mathrm{\θ}(x,y)={\tan }^{-1}\[\frac{L(x,y+1)-L(x,y-1)}{L(x+1,y)-L(x-1,y)}\];$

同时出于对稳定性的考虑，需对关键点的领域点梯度大小做高斯加权处理，梯度加权系数w可采用如下公式确定，i、j为像素点在以关键点为原点的坐标系中所处的坐标值，exp_denom为关键点尺度系数：

w＝exp(-(i·i+j·j)/exp_denom)。

2)采用SIFT特征描述与PCA算法相结合的方式分别对多光谱影像的SIFT特征点以及全色卫星影像的SIFT特征点做降维处理。

本发明实施例中，进一步运用SIFT特征描述与PCA相结合的方式(PCA-SIFT特征点描述方法)，对SIFT特征点进行描述，从而改进SIFT特征的影像匹配过程。

基于PCA-SIFT描述符与传统的标准SIFT描述符的相同点如：有相同的尺度(Scale)，相同的亚像素位置(Sub－pixel)和相同的主方向(DominantOrientations)以及PCA分析法对传统SIFT的128维特征向量做降维处理，得到更加精确的特征描述。降维具体方法如下：

a)输入两幅待匹配图像中所有关键点(设为n个)的128维SIFT特征描述符，将这n个特征描述符作为样本，得到样本矩阵[x₁,x₂,...,x_n]·T，其中x_i表示第i个特征点的128维特征向量。

b)计算n个样本的平均特征向量

c)将所有样本点的特征向量与平均特征向量求差，得到差值向量 $d_i=x_i-\stackrel{\‾}{x},i=1,2,...,n.$

d)构建协方差矩阵 $C=\frac{1}{n}\underset{i=1}{\overset{n}{\Σ}}{d}_i{d}_i^T=\frac{1}{n}{\mathrm{QQ}}^T,$ 其中Q＝[d₁,d₂,...,d_n]。

e)求协方差矩阵的128个特征值λ_i和128个特征向量e_i。

f)将求出的128个特征值按从小到大的顺序进行排列λ₁≥λ₂≥…≥λ₁₂₈和对应的特征向量[e₁,e₂,…,e₁₂₈]。

g)选取对应t个最大特征值的特征向量作为主成分的方向。示例性的，可设置t＝20。

h)构造一个128×t的矩阵A，它的列由t个特征向量组成。

i)把原始的128维SIFT描述符依据y_i＝x_i·A投影到所计算出的n维子空间M中，从而得到PCA－SIFT的描述符y₁,y₂,…,y_n。

2、对降维后的多光谱影像的SIFT特征点以及全色卫星影像的SIFT特征点依次进行粗匹配、误匹配消除及精确匹配处理。具体过程如下：

1)采用BBF算法(最优节点优先算法)对降维处理后的SIFT特征点进行粗匹配，即对降维处理后的多光谱影像的SIFT特征点以及全色卫星影像的SIFT特征点进行基于堆的优先队列方法搜索匹配。具体来说：

BBF算法在KD-tree基础上使用优先队列方法。该方法先计算结点与被查询结点间距离，并对其进行排序。而后依据该距离由小到大的方式进行搜索；在搜索到一分支上的独立结点时，优先队列记忆该结点。这一记忆过程包括记录结点成员在KD-tree中的位置及其与查询结点间距离。搜索完成后，从优先队列中删掉该结点信息。随后，以该节点的最近邻结点为中心，朝其他分支方向搜索。当我们设定最近邻距离/次最近邻距离(例如，Ratio＜0.8)时，即可找到所需要的最近邻，又保持了较高效率。

2)通过RANSAC算法，对匹配结果进行消除误匹配；

本发明实施例中，基于随机抽样一致(RANSAC)算法对这些匹配点进行分类，筛选出正确的匹配点，从而消除无匹配，减少误差。具体过程如下：

a)考虑一个最小抽样集的势为n'的模型(n'为初始化模型参数所需的最小样本数)和一个样本集P，集合P的样本数#(P)＞n'，从P中随机抽取包含n'个样本的P的子集S初始化模型M'；

b)余集SC＝P/S中与模型M的误差小于某一设定阈值c的样本集以及S构成S*。S*认为是内点集，它们构成S的一致集；

c)若#(S*)≥N，认为得到正确的模型参数，并利用集S*，采用最小二乘等方法重新计算新的模型M*；重新随机抽取新的S，重复以上过程。

d)在完成一定的抽样次数后，若未找到一致集则算法失败，否则选取抽样后得到的最大一致集判断内外点，算法结束。

由上可知存在两个可能的算法优化策略。优化算法一是，如果在选取子集S时可以根据某些已知的样本特性等采用特定的选取方案或有约束的随机选取来代替原来的完全随机选取。另一种算法则是通过一致集S*计算出模型M*后，将P中所有与模型M*的误差小于c的样本加入S*，然后重新计算M*。

3)基于仿射变换模型，对误匹配消除后的多光谱影像的SIFT特征点以及全色卫星影像的SIFT特征点进行平移、缩放和旋转处理，从而完成影像的精确匹配。

仿射变换是通过平移、旋转、缩放等变换，将二维平面的平行线映射到新的二维空间。具体变换公式如下：

1)平移变换

$\left[\begin{array}{c}{x}^{\′}\\ y^{\′}\\ 1\end{array}\right]=\left[\begin{array}{ccc}1& 0& t_x\\ 0& 1& t_y\\ 0& 0& 1\end{array}\right]\left[\begin{array}{c}x\\ y\\ 1\end{array}\right]=\left[\begin{array}{c}x+t_x\\ y+t_y\\ 1\end{array}\right]=T(t_x,t_y)\left[\begin{array}{c}x\\ y\\ 1\end{array}\right];$

2)缩放变换

$\left[\begin{array}{c}{x}^{\′}\\ y^{\′}\\ 1\end{array}\right]=\left[\begin{array}{ccc}{S}_x& 0& 0\\ 0& S_y& 0\\ 0& 0& 1\end{array}\right]\left[\begin{array}{c}x\\ y\\ 1\end{array}\right]=\left[\begin{array}{c}{S}_x\·x\\ S_y\·y\\ 1\end{array}\right]=S(S_x,S_y)\left[\begin{array}{c}x\\ y\\ 1\end{array}\right]$

3)旋转变换

围绕原点进行θ角度旋转，变换形式如下。

$\left[\begin{array}{c}{x}^{\′}\\ y^{\′}\\ 1\end{array}\right]=\left[\begin{array}{ccc}cos\mathrm{\θ}& -sin\mathrm{\θ}& 0\\ sin\mathrm{\θ}& \cos \mathrm{\θ}& 0\\ 0& 0& 1\end{array}\right]\left[\begin{array}{c}x\\ y\\ 1\end{array}\right]=\left[\begin{array}{c}xcos\mathrm{\θ}-ysin\mathrm{\θ}\\ xsin\mathrm{\θ}+ycos\mathrm{\θ}\\ 1\end{array}\right]=R\left(\mathrm{\θ}\right)\left[\begin{array}{c}x\\ y\\ 1\end{array}\right]$

逆时针旋转θ角取正值，顺时针旋转θ角取负值。

在进行图像的空间几何变换时，图形要做一次以上的几何变换，可以将各个变换矩阵综合起来进行一步到位的变换，即复合变换。复合变换包括复合平移、复合缩放和复合旋转。复合平移顾名思义就是两次或多次平移，它就是将原图像平移两次并叠加。其变换公式如下：

$\begin{array}{c}T(t_{x2},t_{y2})\·T(t_{x1},t_{y1})=\left[\begin{array}{ccc}1& 0& t_{x2}\\ 0& 1& t_{y2}\\ 0& 0& 1\end{array}\right]\left[\begin{array}{ccc}1& 0& t_{x1}\\ 0& 1& t_{y1}\\ 0& 0& 1\end{array}\right]\\ \left[\begin{array}{ccc}1& 0& t_{x2}+t_{x1}\\ 0& 1& t_{y2}+t_{y1}\\ 0& 0& 1\end{array}\right]=T(t_{x2}+t_{x1},t_{y2}+t_{y1})\end{array}$

其中，(t_x1,t_y1)为第一次的平移量，(t_x2,t_y2)为第二次的平移量。

复合缩放是指两次连续的缩放。其实质相当于将缩放操作进行相乘。其变换公式如下：

$\begin{array}{c}S(s_{x2},s_{y2})\·S(s_{x1},s_{y1})=\left[\begin{array}{ccc}{s}_{x2}& 0& 0\\ 0& s_{x2}& 0\\ 0& 0& 1\end{array}\right]\left[\begin{array}{ccc}{s}_{x2}& 0& 0\\ 0& s_{x2}& 0\\ 0& 0& 1\end{array}\right]\\ =\left[\begin{array}{ccc}{s}_{x2}\·s_{x1}& 0& 0\\ 0& s_{y2}\·s_{y1}& 0\\ 0& 0& 1\end{array}\right]=S(s_{x2}\·s_{x1},s_{y2}\·s_{y1})\end{array}$

其中，(S_x1,S_y1)为第一次的缩放量，(S_x2,S_y2)为第二次的缩放量。

复合旋转是指两次连续的旋转，即对原图像两次旋转并叠加，求得总的旋转角度。其变换公式如下：

$\begin{array}{c}R\left({\mathrm{\θ}}_2\right)\·R\left({\mathrm{\θ}}_1\right)=\left[\begin{array}{ccc}{\mathrm{cos\θ}}_2& -{\mathrm{sin\θ}}_2& 0\\ {\mathrm{sin\θ}}_2& {\mathrm{cos\θ}}_2& 0\\ 0& 0& 1\end{array}\right]\left[\begin{array}{ccc}{\mathrm{cos\θ}}_1& -{\mathrm{sin\θ}}_1& 0\\ {\mathrm{sin\θ}}_1& {\mathrm{cos\θ}}_1& 0\\ 0& 0& 1\end{array}\right]\\ =\left[\begin{array}{ccc}\cos ({\mathrm{\θ}}_2+{\mathrm{\θ}}_1)& -\sin ({\mathrm{\θ}}_2+{\mathrm{\θ}}_1)& 0\\ \sin ({\mathrm{\θ}}_2+{\mathrm{\θ}}_1)& \cos ({\mathrm{\θ}}_2+{\mathrm{\θ}}_1)& 0\\ 0& 0& 1\end{array}\right]=R({\mathrm{\θ}}_2+{\mathrm{\θ}}_1)\end{array}$

其中，θ₁为第一次的平移量，θ₂为第二次的平移量。

对于图像的缩放和旋转变换，必须考虑参照点。前文提及的变换都以原点做为参考点的。当我们需要对任一点(x_f,y_f)作上述变换时，需将该点移到原点，变换后再将(x_f,y_f)移回。在完成上述复合变换时，先作用的变换矩阵在等式右侧变换矩阵先进行变换，然后才是等式左侧矩阵变换。对于图像上点(x_f,y_f)的旋转与缩放变换，公式如下：

1)关于(x_f,y_f)点的缩放变换

$\begin{array}{c}S(x_f,y_f;s_x,s_y)=T(x_f,y_f)\·S(s_x,s_y)\·T(-x_f,-y_f)\\ =\left[\begin{array}{ccc}1& 0& x_f\\ 0& 1& y_f\\ 0& 0& 1\end{array}\right]\left[\begin{array}{ccc}{s}_x& 0& 0\\ 0& s_y& 0\\ 0& 0& 1\end{array}\right]\left[\begin{array}{ccc}1& 0& -x_f\\ 0& 1& -y_f\\ 0& 0& 1\end{array}\right]\\ =\left[\begin{array}{ccc}{s}_x& 0& x_f(1-s_x)\\ 0& s_y& y_f(1-s_y)\\ 0& 0& 1\end{array}\right]\end{array}$

2)绕(x_f,y_f)点的旋转变换

$\begin{array}{c}R(x_f,y_f;\mathrm{\θ})=T(x_f,y_f)\·R\left(\mathrm{\θ}\right)\·T(-x_f,-y_f)\\ =\left[\begin{array}{ccc}1& 0& x_f\\ 0& 1& y_f\\ 0& 0& 1\end{array}\right]\left[\begin{array}{ccc}\cos \mathrm{\θ}& -\sin \mathrm{\θ}& 0\\ \sin \mathrm{\θ}& \cos \mathrm{\θ}& 0\\ 0& 0& 1\end{array}\right]\left[\begin{array}{ccc}1& 0& -x_f\\ 0& 1& -y_f\\ 0& 0& 1\end{array}\right]\\ =\left[\begin{array}{ccc}\cos \mathrm{\θ}& -\sin \mathrm{\θ}& x_f(1-\cos \mathrm{\θ})+y_f\sin \mathrm{\θ}\\ \sin \mathrm{\θ}& \cos \mathrm{\θ}& y_f(1-\cos \mathrm{\θ})+x_f\sin \mathrm{\θ}\\ 0& 0& 1\end{array}\right]\end{array}$

根据仿射变换的基础变换，两幅卫星遥感影像的匹配可以由单应性矩阵解决，单应性矩阵的形式如下：

$\left[\begin{array}{ccc}a& b& c\\ d& e& f\\ g& h& i\end{array}\right];$

在上述公式中， $\left[\begin{array}{cc}a& b\\ d& e\end{array}\right]$ 用来完成缩放、旋转和对称等变换； $\left[\begin{array}{c}c\\ f\end{array}\right]$ 是对图形进行平移变换；[gh]是对图形作投影变换；[i]则是对图形整体进行缩放变换。

通过对正确的匹配点迭代计算，找到并返回源影像和目的影像之间的平面透视变换矩阵H'，要求矩阵误差最小。具体公式如下：

$H^{\′}=\left[\begin{array}{ccc}{h}_{11}& h_{12}& h_{13}\\ h_{21}& h_{22}& h_{23}\\ h_{31}& h_{32}& h_{33}\end{array}\right];$

$S_i\left[\begin{array}{c}{x}_i^{\′}\\ y_i^{\′}\\ 1\end{array}\right]~H^{\′}\left[\begin{array}{c}{x}_i\\ y_i\\ 1\end{array}\right];$

$\underset{i}{\Σ}{(x_i^{\′}-\frac{h_{11}{x}_i+h_{12}{y}_i+h_{13}}{h_{31}{x}_i+h_{32}{y}_i+h_{33}})}^2+{(y_i^{\′}-\frac{h_{21}{x}_i+h_{22}{y}_i+h_{23}}{h_{31}{x}_i+h_{32}{y}_i+h_{33}})}^2;$

精确匹配图像的具体步骤如下：假定以图像1(I1)作为参考系，将图像2(I2)变换到图像1所在的坐标系，变换后的图像为I'2，则其对应关系为

$\left[\begin{array}{c}{x}^{\′}\\ y^{\′}\\ 1\end{array}\right]=\left(\begin{array}{ccc}{h}_{11}& h_{12}& h_{13}\\ h_{21}& h_{22}& h_{23}\\ h_{31}& h_{32}& 1\end{array}\right)\left[\begin{array}{c}x\\ y\\ 1\end{array}\right]=H\left[\begin{array}{c}x\\ y\\ 1\end{array}\right];$

式中:(x,y)和(x',y')是一对匹配点，则只需4对匹配点即可算出H'，通过前问RANSAC可计算得到精确的平面透视变换矩阵H'，将图像2(I2)中的每一个像素和H'矩阵相乘，可以得到精确的图像匹配结果。公式为：

$\left[\begin{array}{cccc}{x}_1^{\′}& x_2^{\′}& & x_n^{\′}\\ y_1^{\′}& y_2^{\′}& ...& y_n^{\′}\\ 1& 1& & 1\end{array}\right]=\left(\begin{array}{ccc}{h}_{11}& h_{12}& h_{13}\\ h_{21}& h_{22}& h_{23}\\ h_{31}& h_{32}& 1\end{array}\right)\left[\begin{array}{cccc}{x}_1& x_2& & x_n\\ y_1& y_2& ...& y_n\\ 1& 1& & 1\end{array}\right]=H\left[\begin{array}{cccc}{x}_1& x_2& & x_n\\ y_1& y_2& ...& y_n\\ 1& 1& & 1\end{array}\right];$

式中，(x,y)和(x',y')为图像像素位置，此时图像2(I'2)和图像1是一一对应的图像，即完成了影像的匹配。

3、基于改进的HIS变换和小波变换相结合的参数可调融合方法。多光谱HIS变换后的亮度I与全色图像进行小波变换后进行高频替换处理，进而小波逆变换得出亮度分量I'。具体过程如下：

将多光谱影像从RGB颜色空间转换到HIS颜色空间，获得亮度I、饱和度S与色度H分量，并对全色卫星影像进行直方图匹配；

分别对多光谱影像的亮度分量及全色卫星影像匹配后的直方图进行小波变换，并将亮度分量小波变换后的高频部分移去，用全色卫星影像匹配后的直方图小波变换后的高频部分进行替换，再对替换结果进行小波逆变换，从而得到新的亮度分量。

其中，对多光谱影像和高分辨率全色影像进行精确匹配，并把多光谱影响的3个波段分别赋予R、G、B3种颜色，把R、G、B彩色影像转换到I、H、S空间。公式如下：

$\left[\begin{array}{c}I\\ V_1\\ V_2\end{array}\right]=\left[\begin{array}{ccc}\frac{1}{\sqrt{3}}& \frac{1}{\sqrt{3}}& \frac{1}{\sqrt{3}}\\ \frac{1}{\sqrt{6}}& \frac{1}{\sqrt{6}}& -\frac{1}{\sqrt{6}}\\ \frac{1}{\sqrt{2}}& -\frac{1}{\sqrt{2}}& 0\end{array}\right]\×\left[\begin{array}{c}R\\ G\\ B\end{array}\right];$

$H=\arctan \left(\frac{V_1}{V_2}\right);S=\sqrt{V_1^2+V_2^2};$

然后遥感图像的Mallats小波分解与重构，通过Mallats小波的二维分解，可以知道，当将图像作为二维信号时，任意二维图像f(x,y)∈L²(R×R)可小波分解，为二维空间尺度函数，ψ_j,m,n(x,y)二维中的小波函数：

$D_j^{1}f=<f(x,y),{\mathrm{\&psi;}}_{j,m,n}^1(x,y)>,(m,n)\&Element;Z^2$

$D_j^{2}f=<f(x,y),{\mathrm{\&psi;}}_{j,m,n}^2(x,y)>,(m,n)\&Element;Z^2$

$D_j^{3}f=<f(x,y),{\mathrm{\&psi;}}_{j,m,n}^3(x,y)>,(m,n)\&Element;Z^2$

如果在2^j分辨率下，f(x,y)可以被分解为A_jf、这4个子图。其中A_jf为图像的低频部分，包含着原图像主要的内容，承载了大部分能量，是原图像的细节信息，用LL表示。为原图的水平方向低频部分与垂直方向高频部分信息，可以反映出水平方向的主要边缘细节信息，用LH表示。承载了原图的水平方向高频部分与垂直方向低频部分信息，可以反映出垂直方向的主要边缘细节信息，用HL表示。承载了原图的对角线方向高频信息，用HH表示。其分解变化如图3所示，其中，LL、LH、HH、HL后面的0、1、2表示分割的层数。

图像的Mallat小波分解与重构过程图如下4-5所示。

当对遥感影使用小波变换进行融合时，设为N层分解，可以得到3N+1个不同频带。这些频带中包含有3N个高频带和1个低频带。其融合步骤为：

(1)对原始影像进行小波滤波，将其分解为一个低频部分和包含水平、垂直及对角线方向的3个高频部分，根据分解层数对低频子图像不断重复小波滤波，建立塔形分解。

(2)对各个分解层分别进行融合处理，融合采用的融合算子可根据不同的融合目的和需求进行选取，并最终得到融合后的小波金字塔。

(3)对上一步融合后得到的小波金字塔进行小波逆变换，完成重构，即得到融合后影像。融合原理如图6所示。

4、通过预设的参数对所述新的亮度分量进行更新，再将更新后的亮度分量与所述饱和度与色度分量进行HIS逆变换，从而实现多光谱影像和全色卫星影像的融合。

当多光谱图像HIS变换后，分解的I分量包含了空间信息，H和S分量包含了光谱信息。不对H和S分量进行考虑，而将亮度I部分和全色图像做小波融合，得到新的亮度分量I’。此时引入参数a、b、c，将参数a和亮度分量I’相乘，之后，将其与全色图像的的亮度分量I_PAN与参数b的乘积、多光谱影像的亮度分量I_MS与参数c乘积三者一起求和，获得了最终的亮度分量I″。这一相加求和过程如下式所示：

I″＝a·I′+b·I_PAN+c·I_MS

其中，参数a、b、c均为常数；I′为所述新的亮度分量，I_PAN为全色卫星影像的亮度分量，I_MS为多光谱影像通过HIS变换获得亮度分量，I″为更新后的亮度分量(即最终的亮度分量)。

由于多光谱影像中富含光谱信息，但相应的缺少高分辨率的空间细节信息。而全色影像是高分辨率影像，其具有丰富的空间几何细节信息，但光谱质量较差。对于基于HIS变换和小波变换相结合的参数可调融合算法而言，由上节中的公式为入口进行分析。当参数a＝0，b＝0，c＝1时，仅仅进行了多光谱影像的HIS变换与逆变换，并未进行融合处理，故不存在分析意义。当参数a＝0，b＝1，c＝0时，即为使用了全色影像直接替代多光谱影像HIS变换后的亮度I部分，为HIS变换融合算法。当参数a＝1，b＝0，c＝0时，则该算法即为HIS变换与小波变换相结合的传统算法。

本发明实施例所提供的上述方法，实际是在传统算法上通过对参数b和c的引入，将原全色影像和多光谱影像再次引入到融合算法中，通过在HIS逆变换完成前引入原有多光谱及全色影像，来对小波变换过程中损失的空间细节信息及光谱信息进行补偿，从而通过简单易行的参数引入方法，进一步提升融合影像的光谱及空间细节质量。同时，因参数的引用，我们也可以根据对融合影像应用场景不同的需要，对融合影像的光谱及细节质量进行动态的调节。当我们需要较高的分辨率时，可通过增加b系数来提高全色影像在新的亮度分量I″中的占比，从而实现目的。同理，当我们需要较高的光谱质量的融合影像时，可以通过增加c系数，来提高多光谱影像在新的亮度分量I″中的占比，从而提升融合后影像的光谱质量。而对于承载HIS变换与小波变换相结合的融合算法结果的a系数，它同时保有光谱信息和空间细节信息，在融合算法中的主体，也是新的亮度分量I″构成中的主体。

最后，将更新后的亮度分量I″与多光谱影像通过HIS变换获得的饱和度与色度分量HIS逆变换，变换后获得的融合图像即为多光谱影像和全色卫星影像的融合图像。

HIS逆变换公式如下：

$\left[\begin{array}{c}R\\ G\\ B\end{array}\right]=\left[\begin{array}{ccc}\frac{1}{\sqrt{3}}& \frac{1}{\sqrt{6}}& \frac{1}{\sqrt{2}}\\ \frac{1}{\sqrt{3}}& \frac{1}{\sqrt{6}}& -\frac{1}{\sqrt{2}}\\ \frac{1}{\sqrt{3}}& -\frac{2}{\sqrt{6}}& 0\end{array}\right]\&times;\left[\begin{array}{c}{I}^{\&prime;\&prime;}\\ V_1\\ V_2\end{array}\right]$

通过以上的实施方式的描述，本领域的技术人员可以清楚地了解到上述实施例可以通过软件实现，也可以借助软件加必要的通用硬件平台的方式来实现。基于这样的理解，上述实施例的技术方案可以以软件产品的形式体现出来，该软件产品可以存储在一个非易失性存储介质(可以是CD-ROM，U盘，移动硬盘等)中，包括若干指令用以使得一台计算机设备(可以是个人计算机，服务器，或者网络设备等)执行本发明各个实施例所述的方法。

以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明披露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应该以权利要求书的保护范围为准。

Claims (5)

1.一种多光谱影像和全色卫星影像融合方法，其特征在于，包括：

分别提取多光谱影像与全色卫星影像的SIFT特征点，并对提取到的SIFT特征点做降维处理；

采用BBF算法对降维处理后的SIFT特征点进行粗匹配，并对粗匹配后的结果进行误匹配消除，再采用仿射变换进行精确匹配；

对精确匹配后的多光谱影像进行HIS变换获得亮度、饱和度与色度分量，将其中的亮度分量以及精确匹配后的全色卫星影像的直方图分别进行小波变换，并将二者的小波变化结果依次进行高频替换、小波逆变换得到新的亮度分量；

通过预设的参数对所述新的亮度分量进行更新，再将更新后的亮度分量与所述饱和度与色度分量进行HIS逆变换，从而实现多光谱影像和全色卫星影像的融合。

2.根据权利要求1所述的方法，其特征在于，所述分别提取多光谱影像与全色卫星影像的SIFT特征点，并对提取到的SIFT特征点做降维处理包括：

分别构建光谱影像与全色卫星影像的SIFT空间金字塔，利用SIFT空间金字塔求得极值点，并对极值点进行筛选，以去除不稳定极值点，再根据筛选后的极值点进行方向和尺度的计算，从而获得多光谱影像的SIFT特征点以及全色卫星影像的SIFT特征点；

采用SIFT特征描述与PCA算法相结合的方式分别对多光谱影像的SIFT特征点以及全色卫星影像的SIFT特征点做降维处理。

3.根据权利要求1所述的方法，其特征在于，所述采用BBF算法对降维处理后的SIFT特征点进行粗匹配，并对粗匹配后的结果进行误匹配消除后，再采用仿射变换进行精确匹配包括：

对降维处理后的多光谱影像的SIFT特征点以及全色卫星影像的SIFT特征点进行基于堆的优先队列方法搜索匹配；

然后，通过RANSAC算法，对匹配结果进行消除误匹配；

再基于仿射变换模型，对误匹配消除后的多光谱影像的SIFT特征点以及全色卫星影像的SIFT特征点进行平移、缩放和旋转处理，从而完成影像的精确匹配。

4.根据权利要求1所述的方法，其特征在于，所述对精确匹配后的多光谱影像进行HIS变换获得亮度、饱和度与色度分量，将其中的亮度分量以及精确匹配后的全色卫星影像的直方图分别进行小波变换，并将二者的小波变化结果依次进行高频替换、小波逆变换得到新的亮度分量包括：

将多光谱影像从RGB颜色空间转换到HIS颜色空间，获得亮度、饱和度与色度分量，并对全色卫星影像进行直方图匹配；

分别对多光谱影像的亮度分量及全色卫星影像匹配后的直方图进行小波变换，并将亮度分量小波变换后的高频部分移去，用全色卫星影像匹配后的直方图小波变换后的高频部分进行替换，再对替换结果进行小波逆变换，从而得到新的亮度分量。

5.根据权利要求1所述的方法，其特征在于，所述通过预设的参数对所述新的亮度分量进行更新，再将更新后的亮度分量与所述饱和度与色度分量进行HIS逆变换，从而实现多光谱影像和全色卫星影像的融合包括：

利用预设的参数a、b、c与新的亮度分量进行下述计算，从而获得更新后的亮度分量：

I″＝a·I′+b·I_PAN+c·I_MS

其中，参数a、b、c均为常数；I′为所述新的亮度分量，I_PAN为全色卫星影像的亮度分量，I_MS为多光谱影像通过HIS变换获得亮度分量，I″为更新后的亮度分量；

再将更新后的亮度分量I″与多光谱影像通过HIS变换获得的饱和度与色度分量HIS逆变换，变换后获得的融合图像即为多光谱影像和全色卫星影像的融合图像。