CN108921809B - 整体原则下基于空间频率的多光谱和全色图像融合方法 - Google Patents

Abstract

本发明提供了整体原则下基于空间频率的多光谱和全色图像融合方法，设A为多光谱图像，B为全色图像；首先对A进行IHS变换，得到亮度IA、色度HA、饱和度SA；然后对IA和B分别进行NSCT分解；再对IA和B的NSCT系数进行融合；最后根据确定的低通子带系数和带通子带系数，进行NCST逆变换，重构生成融合新的IA分量，称之为Inew，将Inew和HA、SA进行IHS逆变换，转换到RGB空间，得到融合后的新图像。本发明在整体原则下基于空间频率将NSCT多光谱和全色图像相融合，融合图像的质量好，光谱信息分辨率高，信息量丰富，层次分明，具有更好的视觉效果和统计指标，有利于遥感图像提取、分类等后续处理。

Description

整体原则下基于空间频率的多光谱和全色图像融合方法

技术领域

本发明属于遥感科学技术领域，具体的涉及整体原则下基于空间频率的多光谱和全色图像融合方法。

背景技术

遥感是利用任何物体都有不同的电磁波反射或辐射特征，通过获取地面目标反射或辐射出来的电磁波来得到目标信息，并对所获得的信息进行处理，实现目标的定位、识别、定性或定量的描述等。

遥感科学是一门综合性科学，它以空间信息科学、电子信息科学、光学、计算机科学、生物、化学等学科为基础，通过对空间对象的观测，为国民经济和科学研究等人类活动，提供了丰富的遥感影像数据。

Contourlet变换和NSCT变换(非下采样Contourlet变换)在遥感图像处理领域应用的很多，包括去噪、压缩、编码、融合等。因为Contourlet和NSCT变换都能有效的捕捉图像的边缘特征，而图像的边缘部分集中了图像的大部分信息，且是图像的重要特征，这在处理图像的时候就等于抓住了图像的“重心”。

高光谱图像超多的波段数目形成了海量的数据，处理高光谱图像就面临着维数过多带来的波段如何选择、运算效率低等困难。因此在高光谱图像的处理中，一直存在着压缩、波段选择、融合的研究，许多研究成果能起到降低高光谱维数、提高效率的目的。

图像融合的常规方法非常多，已经很成熟，但是并不代表这些方法都很“优良”。目前，融合领域的热点集中在基于多尺度几何分析工具方法上。这是因为NSCT具有平移不变性，分解得到的所有子图具有与源图像相同的像素空间，可以充分利用源图像分割的区域信息指导融合过程；相比小波、拉普拉斯金字塔方法等，非下采样的NSCT具有对图像几何特征更强的表达能力，其中小波变换方法只具有有限的方向性，非下采样的NSCT具有任意角度的方向性，更能“捕捉”图像中“沿”边缘信息。一些学者在基于NCST的图像融合方面做了大量的工作，主要有针对医学图像的融合、红外图像和可见光图像融合、全色图像、高光谱和多光谱图像之间的融合、极化SAR图像的融合、SAR图像和多光谱图像的融合等等。其中大量研究工作都集中在医学图像融合或红外图像和可见光图像的融合，而其它方面涉及较少。

图像融合规则可以分为三类：基于像素的融合、基于邻域窗口的融合和基于区域的融合。基于上述三种基本方法的也有加权的、自适应加权的，也有按照某种标准择优选取的。基于单个像素的融合是比较常见的一种方法，但是这种方法等于隔裂了整个图像的背景，往往需要做一致性的检测。而基于邻域窗口的融合规则以邻域窗口为考察对象，利用当前像素为中心的掩模窗口计算邻域内的统计特性，如方差、灰度梯度、局部灰度能量等。在一个邻域窗口内，统计特性越明显，说明图像灰度级变化越大，细节越丰富。基于区域的融合规则以图像分割后的区域为研究对象，每个区域对应一个目标或目标的一部分，部分分割区域是融合中富含信息的部分，要求两幅融合图像的配准精度必须很高，另外基于区域的融合增加了融合步骤，同时也增加了计算量。

对于低频信息的融合可以是简单的使用一些加权方法，或者自适应的加权方法。这些方案都没有充分考虑图像的边缘、纹理等图像特征，其结果是降低图像的对比度。

综上所述，目前的图像融合方法存在诸多缺陷，在图像融合质量上还有待进一步提高。

发明内容

本发明要解决的技术问题是如何提高图像融合的质量，使融合图像具有更好的视觉效果和统计指标。

为了解决上述技术问题，本发明提供了一种整体原则下基于空间频率的多光谱和全色图像融合方法，其采用了如下的技术方案：

一种整体原则下基于空间频率的多光谱和全色图像融合方法，其特征在于：具体步骤如下：

步骤1：设A代表多光谱图像，B代表全色图像；步骤2：对多光谱图像A进行IHS变换，得到亮度I_A、色度H_A、饱和度S_A三个分量；

步骤3：对I_A和B分别进行NSCT分解；

步骤4：对I_A和B的NSCT系数进行融合；

步骤5：NSCT重构；

根据步骤4中确定的低通子带系数和带通子带系数，进行NCST逆变换，重构生成融合新的I_A分量，称之为I_new；

步骤6：IHS逆变换；

将得到新的分量I_new和步骤2中的色度H_A、饱和度S_A三个分量进行IHS逆变换，转换到RGB空间，得到融合后的新图像。

优选地，所述步骤3中，NSCT分解的尺度层数为3，对应的方向个数是[2,4,8]；I_A图像的NSCT分解系数为L_A和

B图像的NSCT分解系数为L_B和

其中，L_A、L_B分别代表I_A和B的低通子带系数；

分别代表I_A和B的带通子带系数；k＝1,2,3，是多尺度分解层数；相应的j＝2^k＝[1,2；1,2,3,4；1,2,…,8]，是分解层数j下的非下采样分解方向个数，此处的j循环定义。

优选地，所述步骤4中，低通子带系数的融合规则：低频融合采用多光谱图像的低频子带系数L_F＝L_A,L_A是A图像的低通子带系数。

优选地，所述步骤4中，带通子带系数的融合规则为：带通子带将依据整体原则下的区域空间频率为衡量准则，进行融合；图像A和B都按照步骤2中的设置进行NSCT分解，得到带通子带

和

最终是由非下采样方向滤波器分解得到的，但是其大小与图像A、B一样，各个方向子带分别保存了该方向分解的系数，也保存了非支撑区域的信息；除了计算k＝1,2,3、对应的j＝[1,2；1,2,3,4；1,2,…,8]时的

方向子带系数，还计算k＝1,2,3、j＝[1；1；1]时的方向子带系数，这时每层带通子带的方向个数都是1，记为：

的大小也和图像A和B是一样的，，这样的结构可以生成多个，目的是用于保存

方向子带系数。

采用基于像素的绝对值最大的高频系数融合规则，其实质是比较高频子带单个像素的绝对值大小，从中选取绝对值较大的作为融合高频子带的融合结果，其模型用公式(1.1)表示：

式(1.1)中，C_new为融合得到的最大高频系数,C_A为A图像的高频系数,C_B为B图像的高频系数；

基于领域窗口统计特性的自适应加权，是基于像素的一种扩展方法；将当前元素作为中心，划定一个大小为3×3或者5×5的窗口，计算该窗口的统计特性；然后将统计特性作为指标，而后或者进行大小比较，直接取舍，或者计算权重进行加权融合；其模型用公式(1.2)表示：

C_new＝qC_A+pC_B,p+q＝1 (1.2)

C_new为融合得到的最大高频系数,C_A为A图像的高频系数,C_B为B图像的高频系数，p为B图像的高频系数的权，q为A图像的高频系数的权。图像的空间频率反映了图像在空域上的总体活跃程度，空间频率越大，图像越清晰，反之，图像越模糊；大小为M×N图像的空间频率SF的定义是：

其中

M、N为不小于2的正整数；f(m，n)是图像m行n列的像素值。

取大小为3×3或5×5的窗口，按照公式(1.3)、(1.4)和(1.5)计算带通子带

的区域空间频率，相应的结果保存在

基于带通子带依据整体原则下的区域空间频率，计算带通子带系数

如下：

本发明提供的方法克服了现有技术的不足，在整体原则下基于空间频率将NSCT多光谱和全色图像相融合，融合图像的质量好，光谱信息分辨率高，图像清晰，信息量丰富，层次分明，具有更好的视觉效果和统计指标，有利于遥感图像的提取、分类等后续处理；本发明提出的计算方案计算量减小，计算效率更高。

附图说明

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

图1为3×3的邻域窗口示意图；

图2为5×5的邻域窗口示意图；

图3为IHS变换得到的融合图像；

图4为小波变换融合得到的融合图像；

图5为Contourlet变换融合得到的融合图像；

图6为传统的NSCT融合方法得到的融合图像；

图7为本发明方法得到的融合图像。

具体实施方式

以下将参照附图，通过实施例方式详细地描述本发明提供的方法。在此需要说明的是，对于这些实施例方式的说明用于帮助理解本发明，但并不构成对本发明的限定。

本文中术语“和/或”，仅仅是一种描述关联对象的关联关系，表示可以存在三种关系，例如，A和/或B，可以表示：单独存在A，单独存在B，同时存在A和B三种情况，本文中术语“/和”是描述另一种关联对象关系，表示可以存在两种关系，例如，A/和B，可以表示：单独存在A，单独存在A和B两种情况，另外，本文中字符“/”，一般表示前后关联对象是一种“或”关系。

实施例一

本实施例主要介绍整体原则下基于空间频率的多光谱和全色图像融合方法的基本原理和步骤。

一、基本定义

高频系数代表是图像的“明显”特征，反应了图像的细节，如边缘、纹理等特征；NSCT变换又是善于捕捉图像“边缘”的，因此，高频系数的融合很关键。本发明采用基于像素的绝对值最大的高频系数融合规则，其实质是比较高频子带单个像素的绝对值大小，从中选取绝对值较大的作为融合高频子带的融合结果，其模型可以用公式(1.1)表示。

式(1.1)中，C_new为融合得到的最大高频系数,C_A为A图像的高频系数,C_B为B图像的高频系数。

基于领域窗口统计特性的自适应加权，可以说是基于像素的一种扩展方法。它是将当前元素作为中心，通常划定一个大小为3×3或者5×5的窗口，如附图1-2所示(附图1为3×3的邻域窗口示意图，附图2为5×5的邻域窗口示意图)，计算该窗口的统计特性；然后将统计特性作为指标，而后或者进行大小比较，直接取舍，或者计算权重进行加权融合。用公式(1.2)描述加权融合的模型。本实施例选择基于3×3的领域窗口，计算整体下的3×3窗口的空间频率，作为指标，比较大小，取绝对值大的空间频率，作为融合的高频系数。公式(1.2)如下：

C_new＝qC_A+pC_B,p+q＝1 (1.2)

式(1.2)中，C_new为融合得到的最大高频系数,C_A为A图像的高频系数,C_B为B图像的高频系数，p为为B图像的高频系数的权，q为为A图像的高频系数的权；

图像的空间频率反映了图像在空域上的总体活跃程度，空间频率越大，图像越清晰，反之，图像越模糊。

大小为M×N图像的空间频率SF的定义是：

其中

M、N为不小于2的正整数，f(m，n)是图像m行n列的像素值。

NSCT变换的带通子带中绝对值较大的系数对应着一些突变，如图像的边缘、纹理等重要特征信息。选择以空间频率为度量标准进行基于区域的融合，可以反映出图像的较为真实的高频信息。

二、整体原则下基于空间频率的多光谱和全色图像融合方法

本发明的一种整体原则下基于空间频率的多光谱和全色图像融合方法的具体步骤包括：

步骤1：设A代表多光谱图像，B代表全色图像；

步骤2：对多光谱图像A进行IHS变换，得到亮度I_A、色度H_A、饱和度S_A三个分量，然后I_A分量参与下面步骤的NSCT分解、融合和NSCT逆换变换；

步骤3：对I_A和B分别进行NSCT分解；

具体的，本融合方法中，NSCT分解的尺度层数为3，对应的方向个数是[2,4,8]；I_A图像的NSCT分解系数为L_A和

B图像的NSCT分解系数为L_B和

其中，L_A、L_B分别代表I_A和B的低通子带系数；

分别代表I_A和B的带通子带系数；k＝1,2,3，是多尺度分解层数；相应的j＝2^k＝[1,2；1,2,3,4；1,2,…,8]，是分解层数下对应的分解的非下采样分解方向数，此处的j循环定义。

步骤4：对I_A和B的NSCT系数进行融合；

具体的，低通子带系数的融合规则：多尺度分解的时候，发现PAN图像(高分辨率全色图像)空间信息较弱，其低频能量很低，所以低频融合采用多光谱图像的低频子带L_F＝L_A,L_A是A图像的低通子带系数。

具体的，带通子带系数的融合规则：带通子带将依据整体原则下的区域空间频率为衡量准则，进行融合。图像A和B都按照步骤2中的设置进行NSCT分解，得到带通子带

和

最终是由非下采样方向滤波器分解得到的，但是其大小与图像A、B一样，各个方向子带分别保存了该方向分解的系数，也保存了非支撑区域的信息，而这些非支撑区域信息如果参与空间频率的计算将会削弱该方向系数。因此，提出在整体原则下计算区域空间频率，作为融合的衡量标准。除了计算k＝1,2,3、对应的j＝[1,2；1,2,3,4；1,2,…,8]时的

方向子带系数，还计算k＝1,2,3、j＝[1；1；1]时的方向子带系数，这时每层带通子带的方向个数都是1，记为：

的大小也与图像A、B是一样的，保存了第k层的方向系数。取大小为3×3或5×5的窗口Q(兼顾计算量，本实施例取3×3)，按照公式(1.3)、(1.4)和(1.5)计算带通子带

的区域空间频率，相应的保存在

至此，可以基于带通子带将依据整体原则下的区域空间频率，计算带通子带方向融合系数

如下：

步骤5：NSCT重构：根据上述步骤中确定的低通子带L_F和带通子带

进行NCST逆变换，重构生成融合新的I_A分量，称之为I_new。

步骤6：IHS逆变换：将得到新的量度分量I_new和步骤2中的色度H_A、饱和度S_A三个分量进行IHS逆变换，转换到RGB空间，得到融合后的新图像。

本实施例中，将多光谱图像A的亮度分量I_A和全色图像B都分别进行NSCT两种变换(一种是分层数k＝3，l＝[1,2,3]，方向数为[2¹,2²,2³]；另外一种是分解层数k＝3，而l＝[0,0,0]，方向数为[2⁰,2⁰,2⁰])(l是对应的分解层的分解方向的指数，比如对应分解层数的方向数为2⁰，即分解方向数是1)，其中后面的这种分解方法是用来计算整体原则下窗口的空间频率，并作为高频融合的依据；低频子带的融合直接采用了全色图像的低频子带，高频带通子带的融合是依据空间频率的绝对值最大原则；将融合后的各个子带按照NSCT逆变换重构得到新的亮度分量I_new；将I_new、H_A、S_A三个分量进行IHS逆变换得到融合新图像，如附图7所示。

本发明融合图像低频系数直接采用了全色图像的低频系数，其原由是多光谱图像包含的更多的光谱信息；而全色图像具有较高的分辨率，其包含的信息量大，边缘、轮廓和纹理等图像特征更明显深刻，相比于多光谱图像其空间信息更准确，分辨率较高的全色图像更能近似地反映图像。

本发明将NSCT和PCA(主成分分析法)变换结合起来，融合高光谱遥感图像，经过对比分析，获取了质量更好的遥感图像，也将会更有利于遥感图像的提取、分类等后续处理。

实施例二

为了进一步的说明本发明技术方案的技术效果，本实施例将本发明技术方案的处理结果与传统方法的处理结果进行比较分析。IHS变换得到的融合图像如附图3所示。

小波变换融合原理首先是多光谱图像A的IHS变换，然后对其亮度分量I_A和全色遥感图像B进行小波三层分解，小波函数使用的是“db45”，将低频数据进行平均融合，高频数据按照绝对值最大的原则，进行融合，然后小波逆变换得到新的亮度分量I_new，将I_new、H_A、S_A进行IHS逆变换，得到融合新图像，如附图4所示。

Contourlet变换融合方法是经典传统的融合规则，同样将多光谱图像A进行IHS逆变换，将变换得到的亮度分量I_c和全色遥感图像进行Contourlet变换，k＝3的三层分解，其带通子带的方向数目为l＝[1,2,3]；将二者按照低频子带平均值融合，高频子带其绝对值最大的原则；然后按照Contourlet逆变换重构I_c分量，将得到的I_c、H、S三个分量进行IHS逆变换，得到Contourlet变换融合的新图像，如附图5所示。

使用NSCT结合传统的融合规则做了融合试验，比较起来更有说服力。传统的NSCT融合方法其流程和本发明提出的基于整体原则下的空间频率的NSCT融合方案相似，分解层数和带通子带的方向数都是相同的，只是传统的NSCT融合使用较为常用的融合规则，即低频子带去平均值，高频子带按照绝对值最大的原则选取，处理结果如附图6所示。

采用本发明的一种整体原则下基于空间频率的多光谱和全色图像融合方法，其处理结果如附图7所示。

融合图像信息量的丰富和图像质量的改善始终是图像融合的根本目的，也是衡量各种融合方法效果的基本准则。理想的融合过程既对有用信息的摄入，也应有对原有有用信息的保留，融合效果评价应该包括创新性和继承性。

通关观察分析图3-7，进行主观视觉评价。图3是基于彩色空间的融合算法IHS变换，其融合质量是几种方法中最差的，其图像较模糊，空间分辨率相对于其它融合图像来说较低，其原因在于亮度分量发生了变化。整体视觉观察：本发明提出的方法融合图像(附图7)处理结果是最好的，其在光谱信息，分辨率上都较好，图像清晰，层次分明，优于其它算法的融合结果。总体来讲，本发明方法的融合图像、传统NSCT融合图像、Contourlet变换融合图像都优于小波变换融合图像、IHS变换融合图像，质量依次渐差。

为了进一步比较几种融合方法的融合质量，本实施例计算了融合图像的熵、平均梯度、标准差，具体结果见表1。

通过分析表1，可以发现客观评价指标中，本发明提出的算法的融合结果较好，这与上述主观视觉评价是一致的。方法IHS变换、小波变换WT、Contourlet变换、传统的NSCT、本发明方法的融合图像的质量依次变高，融合图像的RGB三个波段的指标依次提高了，表明信息量依次变得更丰富。和WT和HIS变换相比，基于Contourlet和NSCT变换的融合，在光谱信息和空间信息上都得到了提高，效果较优；基于Contourlet和NSCT变换的融合图像，其边缘、细节信息也更加明显，这反映了多尺度几何分析工具在图像处理中的优势；由于NSCT的平移不变性，其结果也由于基于Contourlet变换的融合结果。NSCT变换得到的融合图像解译判读效果更佳，对于后续工作分类，提取处理也更为有利。

表1不同方法的统计评价

注：TM代表的是多光谱图像，包含了RGB三个波段。

以上所述，仅为本发明的较佳实施例，并非对本发明任何形式上和实质上的限制，应当指出，对于本技术领域的普通技术人员，在不脱离本发明方法的前提下，还将可以做出若干改进和补充，这些改进和补充也应视为本发明的保护范围。凡熟悉本专业的技术人员，在不脱离本发明的精神和范围的情况下，当可利用以上所揭示的技术内容而做出的些许更动、修饰与演变的等同变化，均为本发明的等效实施例；同时，凡依据本发明的实质技术对上述实施例所作的任何等同变化的更动、修饰与演变，均仍属于本发明的技术方案的范围内。

Claims (1)

1.一种整体原则下基于空间频率的多光谱和全色图像融合方法，其特征在于：该方法的具体步骤包括：

步骤1：设A代表多光谱图像，B代表全色图像；

步骤2：对多光谱图像A进行IHS变换，得到亮度I_A、色度H_A、饱和度S_A三个分量；

步骤3：对I_A和B分别进行NSCT分解；

步骤4：对I_A和B的NSCT系数进行融合；

步骤5：NSCT重构；根据步骤4中确定的低通子带系数和带通子带系数，进行NCST逆变换，重构生成融合新的I_A分量，称之为I_new；

步骤6：IHS逆变换；将得到新的分量I_new和步骤2中的色度H_A、饱和度S_A三个分量进行IHS逆变换，转换到RGB空间，得到融合后的新图像；

所述步骤3中，NSCT分解的多尺度分解层数为3，对应的方向个数是[2,4,8]；I_A图像的NSCT分解系数为L_A和

B图像的NSCT分解系数为L_B和

其中，L_A、L_B分别代表I_A和B的低通子带系数；

分别代表I_A和B的带通子带系数；k＝1,2,3，是多尺度分解层数；相应的j＝2^k＝[1,2；1,2,3,4；1,2,…,8]，是分解层数下的非下采样分解方向个数；

所述步骤4中包括低通子带系数的融合，其规则为：低频融合采用全色图像的低频子带系数L_F＝L_B，L_B是B图像的低通子带系数；

所述步骤4中还包括带通子带系数的融合，其规则为：带通子带将依据整体原则下的区域空间频率为衡量准则，进行融合；图像A和B都按照步骤3中的设置进行NSCT分解，得到带通子带分解系数

和

最终是由非下采样方向滤波器分解得到的，但是其大小与图像A、B一样，各个方向子带分别保存了该方向分解的系数，也保存了非支撑区域的信息，而这些非支撑区域信息如果参与空间频率的计算将会削弱该方向系数；计算k＝1,2,3,对应的j＝[1,2；1,2,3,4；1,2,…,8]时的方向子带系数

计算k＝1,2,3,j＝[1；1；1]时的方向子带系数

这时每层带通子带的方向个数都是1；

图像的空间频率反映了图像在空域上的总体活跃程度，空间频率越大，图像越清晰，反之，图像越模糊；大小为M×N图像的空间频率SF的定义是：

其中

式中，M、N为不小于2的正整数，f(m,n)是图像m行n列的像素值；取大小为3×3或5×5的窗口，按照公式(1.3)、(1.4)和(1.5)计算

的区域空间频率，相应的结果保存在

依据整体原则下的区域空间频率，计算带通子带系数

如下：

Images