CN107169946B - 基于非负稀疏矩阵与超球面彩色变换的图像融合方法 - Google Patents

Abstract

本发明提供了一种基于非负稀疏矩阵与超球面彩色变换的图像融合方法，涉及图像融合领域，本发明通过使用一种非负矩阵分解算法对全色和多光谱图像共同提取亮度分量，然后对亮度分量进行建模计算并调整，再使用超球面彩色变换对图像进行融合，得到融合结果图像，本发明利用NMF算法将全色图像和多光谱图像结合起来提取I分量，在解决光谱匹配性差的问题的同时，提高了亮度分量的提取精度，通过利用全色图像及其滤波后的图像对I分量进行调整，更好的融入了空间细节信息，最大限度地防止了光谱畸变，使得新型卫星的融合结果图像在空间细节信息融入方面和光谱特性保持方面都较现有算法有很大的提高，主观评价与客观分析结果能够达到一致。

Description

基于非负稀疏矩阵与超球面彩色变换的图像融合方法

技术领域

本发明涉及图像融合领域，尤其是一种星载全色与多光谱图像融合方法。

背景技术

近年来，搭载全色与多光谱图像传感器的卫星相继升空，全色图像与多光谱图像的融合研究随之成为卫星遥感领域一个重要而活跃的研究方向，在遥感领域里，往往需要同时具有高空间分辨率和高光谱分辨率的多光谱图像，图像融合技术就是综合全色与多光谱图像的特点与优势，利用高空间分辨率的全色图像去提高多光谱图像的空间分辨率，同时尽量保持多光谱图像的光谱特性，从而克服单一卫星传感器获取信息的局限性，获得更加详细的场景信息。

2009年10月，新型地球观测卫星WorldView-2由美国DigitalGlobe公司发射升空，该新型卫星在高分辨图像上的表现非常突出：能够同时提供8波段1.84m分辨率的多光谱图像和单波段0.46m分辨率的全色图像。与前期卫星相比，如，Landsat、SPOT、IKONOS、QuickBird、THEOS、GeoEye-1等，WorldView-2卫星在高分辨图像上的突出表现集中在：图像波段数增多、分辨率更高，且光谱范围更宽、光谱划分更细，这也导致全色图像与多光谱图像的光谱覆盖范围变窄，光谱匹配性差。随着美国DigitalGlobe与GeoEye两公司于2013年1月完成合并，后续研制的新型卫星表现在遥感图像中的上述特点将会更加突出。例如，2014年8月发射的WorldView-3卫星能够提供0.31m分辨率的全色图像和1.24m分辨率的8波段多光谱图像；2016年11月11日发射的WorldView-4卫星是美国DigitalGlobe与GeoEye两公司合并之后的又一重磅产品，其传感器与WorldView-3具有相同的分辨率，即0.31m的全色波段分辨率和1.24m的多光谱分辨率，但WorldView-4的高分辨率影像能支持用户进行前所未有的精确的2D或3D地图制图、变化监测和图像分析。由此可以看出，WorldView系列的新型卫星图像代表了下一代超高分辨率遥感图像的发展趋势，但同时也使现有的融合方法适应性变差，难以解决新型卫星“光谱匹配变差”的问题，因此融合效果不佳。

关于全色与多光谱图像的融合方法，国内外学者进行了广泛深入的研究，其中大都具有提取亮度分量的过程。例如，PCA变换算法能够实现空间细节信息的有效融入，但却由于丢失部分多光谱图像的光谱信息，而造成光谱失真问题。MSFIM算法是一种基于亮度平滑滤波调制(Smoothing Filter-based Intensity Modulation，SFIM)算法的改进算法，这种改进算法虽然改善了细节信息的融入效果，却产生了较SFIM算法更大的光谱畸变。为降低MSFIM方法中的光谱畸变，则需要改变I分量与全色图像的比值，使其比值更接近于1，也就是使I分量与全色图像的光谱响应特性更相似。HCS变换算法虽然具有良好的光谱保持特性，但却在提取亮度分量时丢失了部分空间细节信息，导致融合效果在细节信息融入方面不够理想。综上所述，目前关于全色与多光谱图像的融合中大都具有提取亮度分量的过程，该过程所用方法主要包括均值法、加权均值法和各波段的算术平均法等，这些方法主要针对传统的IKONOS、QuickBird等4波段星载多光谱与全色图像为研究对象展开研究，但是此类方法应用在WorldView系列新型卫星的遥感数据时，由于存在“光谱匹配差”的问题，所以提取的亮度分量不准确，从而导致空间信息融入时融合结果的光谱畸变明显。

如上所述，现有的全色与多光谱图像融合方法针对的是传统卫星遥感图像，对于WorldView-2/WorldView-3新型卫星的全色与多光谱图像而言，现有的亮度分量提取方法使得提取的亮度分量不够精确，并不能较好的解决空间细节信息融入和光谱畸变的问题，融合效果欠佳。

发明内容

为了克服现有技术的不足，本发明通过使用一种非负矩阵分解(Non-negativeMatrix Factorization，NMF)的算法对全色和多光谱图像共同提取亮度分量，然后对亮度分量进行建模计算并调整，再使用超球面彩色变换(Hyperspherical Color Sharpening，HCS)变换对图像进行融合，得到融合结果图像。本发明方法在空间细节信息融入和光谱特性保持方面都较现有方法有着明显提升，最终得到了高质量的融合结果图像。

本发明解决其技术问题所采用的技术方案的详细融合步骤如下：

步骤1.使用NMF算法提取亮度分量I

首先，将全色图像PAN和八波段的多光谱图像X₁,X₂,...,X₈按行拉直成列向量，再按公式(1)组成待分解矩阵V，即

V＝[P,M₁,M₂,...,M₈] (1)

其中，P,M₁,M₂,...,M₈分别为全色图像PAN和多光谱图像的八个波段X₁,X₂,...,X₈的图像矩阵按行拉直所得的对应列向量；

令

[P,M₁,M₂,...,M₈]＝WH (2)

其中，W为n×r矩阵，n为矩阵W的行数，r为矩阵W的列数，H为r×9矩阵，取r＝1，则分解所得的W为一个列向量，将分解得到的列向量W中的值按行依次恢复到原图像中的对应位置，构成新的图像矩阵作为I分量；

步骤2.使用分量PAN及其滤波后分量PAN_smooth调整I分量

首先利用滤波器对全色图像PAN进行平滑滤波，得到滤波后的全色图像PAN_smooth：

PAN_smooth＝SMOOTH(PAN) (3)

其中，SMOOTH(PAN)表示对原全色图像进行平滑滤波；

令

P'²＝(PAN_smooth)² (4)

P"²＝(PAN)² (5)

即用P'变量代替PAN_smooth，用P"变量代替PAN；

利用上述分量P'²、分量P"²对步骤1所得的分量I²进行匹配，匹配公式如下：

其中，μ₀、σ₀分别为分量I²的均值和标准方差，μ₁、σ₁分别为分量P'²的均值和标准方差，则调整后的分量如下：

其中，I_adj表示对I分量进行调整后得到的新的分量；

步骤3.使用HCS变换最终得到融合结果图像

首先，对八波段的多光谱图像X₁,X₂,...,X₈进行HCS正变换得到相应的I分量和角度分量

HCS正变换如下：

其次，对步骤2求得的分量I_adj和角度分量

做HCS反变换得到新的八波段分量X₁',X₂',...,X₈'，HCS反变换如下：

波段分量X₁',X₂',...,X₈'对应八波段的高空间分辨率多光谱图像，即为得到的融合结果图像。

本发明的有益效果在于由于利用NMF算法将全色图像和多光谱图像结合起来提取I分量，在解决光谱匹配性差的问题的同时，提高了亮度分量的提取精度，较传统算法中I分量的提取方法更合理，通过利用全色图像及其滤波后的图像对I分量进行调整，更好的融入了空间细节信息，利用HCS变换将亮度信息和角度信息分开处理，最大限度地防止了光谱畸变，使得WorldView-2/WorldView-3新型卫星的融合结果图像在空间细节信息融入方面和光谱特性保持方面都较现有算法有很大的提高，主观评价与客观分析结果能够达到一致，另外，本发明不存在频率分解与重建过程，计算简洁、实时性好，对于需要快速交互处理和实时可视化的融合系统具有重要的实用参考价值。

附图说明

图1为本发明WorldView-2卫星原始图像和不同方法融合结果，其中图1(a)是全色图像，图1(b)是多光谱图像，图1(c)是参考图像，图1(d)是PCA融合结果图像，图1(e)是WT+PCA融合结果图像，图1(f)是

融合结果图像，图1(g)是HCS-smart融合结果图像，图1(h)是NMF+HCS融合结果图像。

图2为本发明WorldView-3卫星原始图像和不同方法融合结果，其中图2(a)是全色图像，图2(b)是多光谱图像，图2(c)是参考图像，图2(d)是PCA融合结果图像，图2(e)是WT+PCA融合结果图像，图2(f)是

融合结果图像，图2(g)是HCS-smart融合结果图像，图2(h)是NMF+HCS融合结果图像。

其中，PCA为使用PCA融合方法得到的融合结果图像，WT+PCA为小波变换+PCA融合方法得到的融合结果图像，

为直接提取I分量的

模型方法得到的融合结果图像，HCS-smart为直接提取I分量的smart模型方法得到的融合结果图像，NMF+HCS为使用本发明方法得到的融合结果图像。

具体实施方式

下面结合附图和实施例对本发明进一步说明。

非负矩阵分解作为一种相对成熟的数据分析手段，广泛应用于图像分析、数据挖掘、语音处理、人脸识别等诸多领域。非负矩阵分解所得基矩阵具有一定的线性无关性和稀疏性，能有力表达原始数据的特征及结构，能自动发现数据背后的隐含模式和趋势。考虑到从多光谱图像中提取的I分量要尽可能的与全色图像的光谱响应相一致，因此在利用NMF算法提取I分量时，用全色图像和八波段的多光谱图像共同组成待分解矩阵。

HCS融合方法对图像波段数没有限制，因此适合新型WorldView-2/WorldView-3卫星多波段遥感图像的融合。由于在HCS变换中，角度变量决定图像的光谱信息，I分量决定图像的亮度信息，所以I分量的变化对图像的光谱信息没有影响，这里重点研究I分量的替换策略。

步骤1.使用NMF算法提取亮度分量I

首先，将全色图像PAN和八波段的多光谱图像X₁,X₂,...,X₈按行拉直成列向量，再按公式(1)组成待分解矩阵V，即

V＝[P,M₁,M₂,...,M₈] (1)

其中，P,M₁,M₂,...,M₈分别为全色图像PAN和多光谱图像的八个波段X₁,X₂,...,X₈的图像矩阵按行拉直所得的对应列向量；

令

[P,M₁,M₂,...,M₈]＝WH (2)

其中，W为n×r矩阵，n为矩阵W的行数，r为矩阵W的列数，H为r×9矩阵，取r＝1，则分解所得的W为一个列向量，将分解得到的列向量W中的值按行依次恢复到原图像中的对应位置，构成新的图像矩阵作为I分量；

步骤2.使用分量PAN及其滤波后分量PAN_smooth调整I分量

首先利用滤波器对全色图像PAN进行平滑滤波，得到滤波后的全色图像PAN_smooth：

PAN_smooth＝SMOOTH(PAN) (3)

SMOOTH(PAN)表示对原全色图像进行平滑滤波，本发明使用均值滤波器，其中，滤波器模板的大小可以进行适当调整，但经过多次实验发现，将滤波器模板设置为7*7大小时，融合结果的空间伪影最小，使用比7*7更小或更大的模板都可能会导致融合结果出现重影和模糊，所以此处滤波器选择模板大小为7*7的均值滤波器。

令

P'²＝(PAN_smooth)² (4)

P"²＝(PAN)² (5)

即用P'变量代替PAN_smooth，用P"变量代替PAN；

利用上述分量P'²、分量P"²对步骤1所得的分量I²进行匹配，匹配公式如下：

其中，μ₀、σ₀分别为分量I²的均值和标准方差，μ₁、σ₁分别为分量P'²的均值和标准方差，则调整后的分量如下：

其中，I_adj表示对I分量进行调整后得到的新的分量；

步骤3.使用HCS变换最终得到融合结果图像

首先，对八波段的多光谱图像X₁,X₂,...,X₈进行HCS正变换得到相应的I分量和角度分量

HCS正变换如下：

其次，对步骤2求得的分量I_adj和角度分量

做HCS反变换得到新的八波段分量X₁',X₂',...,X₈'，HCS反变换如下：

波段分量X₁',X₂',...,X₈'对应八波段的高空间分辨率多光谱图像，即为得到的融合结果图像。

为验证基于非负稀疏矩阵与超球面彩色变换的WorldView-2/WorldView-3全色与多光谱图像融合新方法的有效性与优越性，进行以下实验：将本发明方法与不同的融合方法进行对比分析，对比方法分别为：PCA融合方法记作PCA，小波变换+PCA融合方法记作WT+PCA，直接提取I分量的

模型方法记作

直接提取I分量的smart模型方法记作HCS-smart，本发明方法记作NMF+HCS，实验数据采用了两组真实的新型卫星图像数据，第一组为WorldView-2图像，2011年4月拍摄的California USA，大小为300×300，第二组为WorldView-3图像，2014年10月拍摄的Sydney Australia，大小为512×512，为验证算法并清楚的显示原图像和融合结果图像用于主观评价，实验截取了原图像的一部分，为更加方便的对融合结果图像进行主客观评价，在实验中，先对全色与多光谱图像进行下采样，将全色与多光谱图像的空间分辨率同时降低4倍，作为待融合的模拟全色与多光谱图像数据，再将原始的多光谱图像作为标准融合结果进行参考。

图1和图2中(a)(b)分别为原始全色图像与多光谱图像下采样后得到的模拟全色图像和模拟多光谱图像；(c)为原始多光谱图像，作为参考图像用于实验结果对比；(d)-(g)为上述对比方法对应的融合结果；(h)为本发明方法的融合结果。通过观察图1和图2，从道路及建筑物的边缘细节处以及植被颜色上可以看出，图(h)的细节信息更加丰富，且光谱特征与(c)更为相近，与其他对比方法相比较，显然本发明方法的融合结果与参考图像差异最小，主观视觉效果最好，由此可见，本发明方法融合结果在空间细节信息的融入和光谱特性的保持上都较对比方法有不同程度的提高。

本发明选取空间相关系数(Spatial Correlate Coefficient，sCC)、平均梯度(Average Gradient，AG)、信息熵(InformationEntropy，IE)、交叉熵(Cross Entropy，CE)和边缘强度(Edge Intensity，EI)五种常用的客观评价指标，对不同方法的融合结果进行客观评价。其中，sCC反应细节融合的好坏，其值越大，融合效果越好；AG表示相邻像元的灰度值变化程度，在理想状况下，梯度值越大，图像越清晰，空间细节信息保持越好；IE反映了图像的信息含量，融合图像的熵值越大，说明信息量越丰富；CE用来度量融合结果图像与原始多光谱图像对应像素的差异，交叉熵越小，融合图像与参考图像间差异越小，光谱信息保持越好；EI表示边缘与邻近像素的对比强度，其值越大，对比越强烈，边缘细节信息保持越好。两组卫星图像数据的不同融合结果对应的各项指标大小分别如表1和表2所示。

表1

表2

表1为WorldView-2图像融合结果客观评价结果，表2为WorldView-3图像融合结果客观评价结果，由表1和表2可以看出，五项评价指标中，本发明方法的各项评价指标均为最优的，即对于sCC、AG、IE、EI四项指标，本发明方法的值均是最大的，说明本方法对空间细节信息融入效果最好，而对于CE指标而言，本发明方法的值是最小的，说明本方法对光谱特性保持最优。综上所述，本发明方法无论在空间细节增强还是光谱保真方面，都较其他对比方法有不同程度的提高，即本发明方法较对比方法具有更好的融合性能。

Claims (1)

1.一种基于非负稀疏矩阵与超球面彩色变换的图像融合方法，其特征在于包括下述步骤：

步骤1.使用NMF算法提取亮度分量I

首先，将全色图像PAN和八波段的多光谱图像X₁,X₂,...,X₈按行拉直成列向量，再按公式(1)组成待分解矩阵V，即

V＝[P,M₁,M₂,...,M₈] (1)

其中，P,M₁,M₂,...,M₈分别为全色图像PAN和多光谱图像的八个波段X₁,X₂,...,X₈的图像矩阵按行拉直所得的对应列向量；

令

[P,M₁,M₂,...,M₈]＝WH (2)

其中，W为n×r矩阵，n为矩阵W的行数，r为矩阵W的列数，H为r×9矩阵，取r＝1，则分解所得的W为一个列向量，将分解得到的列向量W中的值按行依次恢复到原图像中的对应位置，构成新的图像矩阵作为I分量；

步骤2.使用全色图像PAN及其滤波后的全色图像PAN_smooth调整I分量

首先利用滤波器对全色图像PAN进行平滑滤波，得到滤波后的全色图像PAN_smooth：

PAN_smooth＝SMOOTH(PAN) (3)

其中，SMOOTH(PAN)表示对原全色图像进行平滑滤波；

令

P'²＝(PAN_smooth)² (4)

P"²＝(PAN)² (5)

即用P'变量代替PAN_smooth，用P"变量代替PAN；

利用上述分量P'²、分量P"²对步骤1所得的分量I²进行匹配，匹配公式如下：

其中，μ₀、σ₀分别为分量I²的均值和标准方差，μ₁、σ₁分别为分量P'²的均值和标准方差，则调整后的分量如下：

其中，I_adj表示对I分量进行调整后得到的新的分量；

步骤3.使用HCS变换最终得到融合结果图像

首先，对八波段的多光谱图像X₁,X₂,...,X₈进行HCS正变换得到相应的I分量和角度分量

HCS正变换如下：

其次，对步骤2求得的分量I_adj和角度分量

做HCS反变换得到新的八波段分量X₁',X₂',...,X₈'，HCS反变换如下：

波段分量X₁',X₂',...,X₈'对应八波段的高空间分辨率多光谱图像，即为得到的融合结果图像。

Images