CN104933690B - 基于混合亚像元解混的遥感多光谱与全色图像融合方法 - Google Patents

Abstract

本发明公开了一种基于混合亚像元解混的遥感多光谱与全色图像融合方法，包括将原始低分辨率多光谱图像MS ₀ 采用cubic方式上采样到原始全色图像PAN的分辨率，得到MS；并基于该MS计算归一化植被指数NDVI；利用PAN和NDVI，对MS中的亚像元进行分类，得到植被‑非植被混合亚像元C _M 和其他亚像元C _O ；利用PAN和NDVI，判断植被‑非植被混合亚像元的地物类型类别；对于MS中每个亚像元t，根据其对应的亚像元类别和地物类型类别进行融合。本发明的有益效果：通过对多光谱图像中的混合亚像元，在融合过程中进行光谱解混处理，从而降低融合图像中混合像元数量，锐化融合图像边界和降低光谱失真。

Description

基于混合亚像元解混的遥感多光谱与全色图像融合方法

技术领域

本发明涉及遥感图像处理技术领域，具体来说，涉及一种基于混合亚像元解混的遥感多光谱与全色图像融合方法。

背景技术

遥感数据融合是将不同传感器数据所提供的信息加以综合，以获得更全面、准确的图像，它能实现优势互补和消除冗余，克服单一传感器在几何、光谱、空间和时间分辨率等方面的局限，从而有利于影像解译、分类和目标识别。随着航天、传感器和通信等技术的发展，遥感数据在空间、光谱、时间分辨率上快速提高。由于传感器物理特性和数据传输能力等的限制，高空间分辨率和高光谱空间分辨率影像难以同时获取。因此，大量目前在轨的卫星(例如Landsat 7ETM+、QuickBird、SPOT-5和WorldView-2/3等)会同时提供高空间分辨率的全色波段(PAN)和低空间分辨率的多光谱波段(MS)。由于大量的应用需要使用高空间分辨率的多光谱图像，近年来国内外研究人员发展了大量遥感融合技术来融合MS和PAN图像以得到高空间分辨率的MS图像。

现有的MS与PAN融合方法可以分为三类：基于成分替代的方法、基于调制的方法和基于多尺度分析的方法。基于成分替代融合的典型代表算法有IHS变换、PCA,Gram-Schmidt等，而基于调制技术的融合算法主要有Brovey变换、Pradines’、合成变化比率(SyntheticVariable Ratio)、Smoothing Filter-based Intensity Modulation、PANSHARP(PS)等。基于成分替代和调制技术的融合算法的特点是快速且易于实现，但其融合结果会导致一定程度的光谱失真。基于多尺度变换的融合方法，有效地保持了多光谱图像的光谱信息。然而这类方法的融合图像可能会出现空间畸变，典型的现象有振铃效应、虚景混淆、边缘以及纹理模糊。为综合成分替代和多分辨率分析等方法分别在空间和光谱信息保留方面的优势，一些研究人员提出了成分替代(PCA、HIS、ICA等)或调制技术与多分辨率分析相结合的融合算法，这类方法能获得优于成分替代和标准多分辨率分析融合的结果，但增加了计算复杂度。有一些研究利用正则化优化的方法对融合结果进行迭代。尽管国内外研究人员提出了大量新的融合方法来降低融合图像光谱失真，如何进一步降低光谱失真仍然是遥感融合研究面临的主要挑战。

由于MS和PAN图像空间分辨率的差异，低分辨率MS图像的一个混合像元在重采样到PAN分辨率的MS图像中对应多个亚像元。这些亚像元中部分像元对应PAN图像中的纯像元，另一部分则对应PAN图像中的边界像元。在现有融合方法中，对应PAN纯像元的MS混合亚像元在融合图像中仍然为混合像元，这使得融合图像中植被和非植被像元的边界模糊，存在较大的光谱失真。因此，为了改善融合图像质量，需要在融合过程中将MS混合亚像元解混到对应PAN纯像元的地物类别。然而，现有融合方法鲜有考虑混合亚像元的解混。在少量考虑了混合亚像元解混的融合方法中，主要使用面对对象分类的方法来获取地物类别信息。但面向对象分类在参数设置、样本训练等方面的复杂性限制了这些方法的实际应用。

针对相关技术中的问题，目前尚未提出有效的解决方案。

发明内容

针对相关技术中的上述技术问题，本发明提出一种基于混合亚像元解混的遥感多光谱与全色图像融合方法，该方法对MS图像中的混合亚像元(即对应PAN纯像元的混合MS像元的亚像元)，在融合过程中进行光谱解混处理，从而降低融合图像中混合像元数量，锐化融合图像边界和降低光谱失真。

本发明的技术方案是这样实现的：

一种基于混合亚像元解混的遥感多光谱与全色图像融合方法，包括以下步骤：

将原始低分辨率多光谱图像(MS₀)采用cubic方式上采样到原始全色图像PAN的分辨率，得到MS；并基于该MS计算归一化植被指数NDVI；

利用PAN和NDVI，对MS中的亚像元进行分类，得到植被-非植被混合亚像元C_M和其他亚像元C_O；

利用PAN和NDVI，判断植被-非植被混合亚像元的地物类型类别，其中，所述地物类型类别包括植被类、非植被类；

对于MS中每个亚像元t，根据其对应的亚像元类别和地物类型类别进行融合。

此外，上述的基于混合亚像元解混的遥感多光谱与全色图像融合方法还包括：将PAN采用平均法下采样到R倍空间分辨率，再采用立方卷积cubic方式上采用到原始分辨率，得到PAN_S；其中，R为MS₀与PAN的空间分辨比率。

另外，在利用PAN和NDVI，对MS中的亚像元进行分类，得到植被-非植被混合亚像元C_M和其他亚像元C_O时，可利用边缘检测Canny算子分别计算PAN_S和NDVI的边界E_PAN和E_NDVI；并利用直径为L_SE1的圆盘结构元素，对E_NDVI进行形态学膨胀得到宽的NDVI边界图像E^W _NDVI；同时将在E_PAN中值为1但在E^W _NDVI中值为0的像元设置为0，得到PAN边界图像E′_PAN；再采用直径为L_SE2的圆盘结构元素对E′_PAN进行形态学膨胀得到宽的PAN边界图像E^W _PAN；利用E^W _PAN，对MS中的亚像元进行分类，得到植被-非植被混合亚像元C_M和其他亚像元C_O。

此外，在利用E^W _PAN，对MS中的亚像元进行分类，得到植被-非植被混合亚像元C_M和其他亚像元C_O时，可利用E^W _PAN，将MS中的亚像元分为两类，得到亚像元分类图，并根据亚像元分类图，得到植被-非植被混合亚像元C_M和其他亚像元C_O；其中，亚像元分类图M_PI的分类规则如下：

其中，如果亚像元t在E^W _PAN中对应值为1，则被分为植被-非植被混合亚像元C_M；如果t在E^W _PAN中对应值为0，则被分为其它亚像元类C_O。

另外，在利用PAN和NDVI，判断植被-非植被混合亚像元的地物类型类别时，可根据PAN边界图像E′_PAN，计算植被边界像元图像和非植被像元图像并利用直径为3的圆盘结构元素，对进行形态学膨胀，得到膨胀后的植被边界像元图像M_N，并将M_N中值为1且在中值为1的像元设置为0，并重复该步骤R-1次，使得M_N中值为1的像元均为植被像元；同时利用直径为3的圆盘结构元素，对进行形态学膨胀，得到M_NV，并将M_NV中值为1且在中值为1的像元设置为0，并重复该步骤R-1次，使得M_NV中值为1的像元均为非植被像元；而对于和中每个值为1的像元t₁，分别在以t₁为中心、大小为S_T×S_T的邻域窗口W(t₁,S_T)内计算所有值为1的像元的NDVI值，并根据该NDVI值，采用以下公式，分别得到图像和

并对于M_PI中每个属于混合亚像元类C_M的像元t₂，根据M_N、M_NV、和采用以下公式判断其属于植被C_V或非植被C_NV类别，得到混合亚像元分类图M_LC：

此外，在根据PAN边界图像E′_PAN，计算植被边界像元图像和非植被像元图像时，可对于每个亚像元t，在其8邻域内查找与其在PAN图像中灰度差异最大的邻居像元t’，得到像元对{t，t’}；其中，如果NDVI(t)＞NDVI(t’)，则且如果NDVI(t)＜NDVI(t’)，则且

另外，在对于MS中每个亚像元t，根据其对应的亚像元类别和地物类型类别进行融合时，可采用暗像元法计算MS和PAN图像的雾气值H_i和H_p，其中，i等于1，……N，N为MS图像的波段数量；并判断亚像元t是否属于其它亚像元类C_O，并在判断结果为是的情况下，利用该亚像元t的原始光谱，采用以下公式进行融合：其中，MS_i为MS的第i波段，为融合图像MS^f的第i波段；同时判断亚像元t是否属于植被-非植被混合亚像元C_M，并在判断结果为是的情况下，利用其邻域内的纯植被或纯非植被像元n的光谱，采用以下公式进行融合：

其中，邻域内的纯植被或纯非植被像元n在以t为中心、大小为S_N×S_N的窗口W(t,S_N)内根据以下公式查找：

其中，S_N的值小于等于L_SE2。

其中，在采用暗像元法计算MS和PAN图像的雾气值H_i和H_p时，可确定PAN_S中最小值对应的位置IDX_min，并利用该IDX_min和以下公式计算H_i和H_p：H_i＝MS_i(IDX_min)，H_p＝PAN_s(IDX_min)。

本发明的有益效果：通过对植被-非植被混合亚像元进行识别并判断其属于植被或非植被类型，在融合过程中根据其类型将其解混为植被或非植被像元，从而降低融合图像中混合亚像元的数量，能有效锐化植被-非植被边界和显著降低光谱失真。

附图说明

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

图1是根据本发明实施例的基于混合亚像元解混的遥感多光谱与全色图像融合方法的流程示意图；

图2是根据本发明实施例的基于混合亚像元解混的遥感多光谱与全色图像融合方法中MS亚像元分类流程图。

具体实施方式

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

根据本发明的实施例，提供了一种基于混合亚像元解混的遥感多光谱与全色图像融合方法。

如图1-2所示，根据本发明实施例的基于混合亚像元解混的遥感多光谱与全色图像融合方法包括以下步骤：

一、将原始低分辨率MS图像(MS₀)采用cubic(立方卷积)方式上采样到原始PAN图像(PAN)的分辨率，得到MS；将PAN采用平均法下采样到R倍空间分辨率，然后再采用cubic方法上采样到原始分辨率，得到PAN_S，其中，R为MS₀与PAN的空间分辨比率，基于MS图像计算归一化植被指数(NDVI)。

二、将MS亚像元分类为植被-非植被边界及其附近混合亚像元和其它亚像元两类。具体实施步骤如下：

步骤1：利用Canny(边缘检测)算子分别计算PAN和NDVI的边界E_PAN和E_NDVI，阈值分别设置为T_NV和T_PAN，通常T_PAN小于T_NV。在E_PAN和E_NDVI中，边界像元和非边界像元的值分别为1(True)和0(False)；

步骤2：利用直径为L_SE1的圆盘结构元素对E_NDVI进行形态学膨胀得到E^W _NDVI，将在E_PAN中为1但在E^W _NDVI中为0的像元设置为0得到PAN边界图像E′_PAN。E′_PAN中值为1的像元主要包含植被和非植被边界像元，以及少量非植被-非植被边界像元。利用NDVI将E′_PAN中值为1的非植被-非植被边界像元的值设置为0。最后得到的E′_PAN中值为1的像元主要包含植被和非植被边界像元；其中L_SE1的值设为2R-1；步骤3：采用直径为L_SE2的圆盘结构元素对E′_PAN进行形态学膨胀得到E^W _PAN。E^W _PAN中值为1的像元包含植被和非植被边界像元(即E′_PAN中为1的像元)，和植被-非植被边界附近的像元(即E^W _PAN中值为1但在E′_PAN中为0的像元)；其中L_SE2的值设为2R-3(L_SE2<L_SE1)；

步骤4：根据E^W _PAN，将MS亚像元分为C_O和C_M二类，采用以下公式：

即如果亚像元t在E^W _PAN中对应值为1，则被分为植被边界及其附近的混合亚像元类C_M；如果t在E^W _PAN中值为0，则被分为其它亚像元类C_O，从而得到亚像元分类图M_PI。

三、利用PAN和NDVI判断植被-非植被混合亚像元类C_M中亚像元的为植被或非植被地物类型。具体实施步骤如下：

1)根据PAN边界图像E′_PAN，计算植被边界像元图像和非植被边界像元图像首先对于每个植被-非植被边界像元t(即满足E′_PAN(t)＝1)，在其8邻域内查找与其在PAN图像中灰度差异最大的邻居像元t’，得到像元对{t，t’}。如果NDVI(t)＞NDVI(t’)，则且如果NDVI(t)＜NDVI(t’)，则且

2)使用直径为3的圆盘结构元素，对进行形态学膨胀，得到M_V，并将M_V中值为1且在中为1的像元设置为0；

3)重复步骤2共R-1次，最后得到的M_V中值为1的像元均为植被像元；

4)使用直径为3的圆盘结构元素，对行形态学膨胀，得到M_NV，并将M_NV中值为1且在中为1的像元设置为0；

5)重复步骤4共R-1次，最后得到的M_NV中值为1的像元均为非植被像元；

6)对于和中每个值为1的像元t₁，分别在以t₁为中心、大小为S_T×S_T的邻域窗口W(t₁,S_T)内计算所有值为1的像元的NDVI均值，采用以下公式：

其中,S_T的值设为2R-1。

7)对于M_PI中每个属于混合亚像元类C_M的像元t₂，根据M_V、M_NV、和采用以下公式判定其属于植被(C_V)或非植被(C_NV)类别，从而得到混合亚像元分类图M_LC。

四、对于MS中每个亚像元t，根据其对应的亚像元类别和地物类型类别进行融合。具体实施步骤如下：

1)采用暗像元法计算MS和PAN图像的雾气值H_i(i＝1，...，N)和H_p，其中N为MS图像的波段数量。具体为，首先搜索PAN_S中最小值对应的位置IDX_min，然后分别采用以下公式计算H_i和H_p：

H_i＝MS_i(IDX_min)

H_p＝PAN_S(ID X_min)。

2)如果t属于其它亚像元类C_O，利用该亚像元t的原始光谱，采用以下公式进行融合：

其中，MS_i为MS的第i波段，为融合图像MS^f的第i波段。

3)如果t属于混合亚像元类C_M，利用其邻域内的纯植被或纯非植被像元n的光谱进行融合，采用以下公式进行融合：

其中，邻居像元n在以t为中心、大小为S_N×S_N的窗口W(t,S_N)内根据如下公式查找：

其中，S_N的值小于等于L_SE2；从而实现植被非植被混合亚像元融合为纯植被或纯非植被像元。

为了更好的了解本发明上述技术方案的性能，以下通过对比试验(选取目前常用的两种图像融合方法PANSHARP方法和Gram-Schmidt方法与本发明方法进行性能对比)，对本发明的上述技术方案进行详细说明。

1.实验数据

使用数据为已配准的512×512多光谱图像和2048×2048全色遥感图像。

2.对比实验及实验评价指标

采用相对平均误差(EASE)、相对全局维数综合指标(ERGAS)、光谱角(SAM)和综合质量指数Q4等常用定量指标从客观上评价了这几个融合方法，如表1所示。其中EASE反映了融合图像跟参考图像的偏差，值越小融合效果越好；ERGAS反映了融合图像跟参考图像的全局光谱辐射变形误差，越小越好；SAM反映融合图像光谱跟参考图像光谱的差异，值越小融合效果越好；Q4是同时考虑了融合图像跟参考图像的局部均值偏差、对比度变化以及相关性丢失情况的综合质量指标，值越大越好。

表1融合图像的质量评价指标统计

通过比较融合图像以及表中的统计指标可见，跟PANSHARP和Gram-Schmidt方法的融合图像相比，本发明方法的融合图像锐化了植被-非植被边界，其光谱失真最小。

由此可见，借助于本发明的上述技术方案，通过对植被-非植被混合亚像元进行识别并判断其属于植被或非植被类型，在融合过程中根据其类型将其解混为植被或非植被像元，从而降低融合图像中混合亚像元的数量，能有效锐化植被-非植被边界和显著降低光谱失真。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (7)

1.一种基于混合亚像元解混的遥感多光谱与全色图像融合方法，其特征在于，包括以下步骤：

将原始低分辨率多光谱图像MS₀采用立方卷积cubic方式上采样到原始全色图像PAN的分辨率，得到MS；并基于该MS计算归一化植被指数NDVI；

利用PAN和NDVI，对MS中的亚像元进行分类，得到植被-非植被混合亚像元C_M和其他亚像元C_O；

利用PAN和NDVI，判断植被-非植被混合亚像元的地物类型类别，其中，所述地物类型类别包括植被类、非植被类；

对于MS中每个亚像元t，根据其对应的亚像元类别和地物类型类别进行融合。

2.根据权利要求1所述的基于混合亚像元解混的遥感多光谱与全色图像融合方法，其特征在于，所述将原始低分辨率多光谱图像MS₀采用立方卷积cubic方式上采样到原始全色图像PAN的分辨率，得到MS；并基于该MS计算归一化植被指数NDVI还包括：

将PAN采用平均法下采样到R倍空间分辨率，再采用立方卷积cubic方式上采用到原始分辨率，得到PAN_S；其中，R为MS₀与PAN的空间分辨比率。

3.根据权利要求2所述的基于混合亚像元解混的遥感多光谱与全色图像融合方法，其特征在于，利用PAN和NDVI，对MS中的亚像元进行分类，得到植被-非植被混合亚像元C_M和其他亚像元C_O包括：

利用边缘检测Canny算子分别计算PAN_S和NDVI的边界E_PAN和E_NDVI；

利用直径为L_SE1的圆盘结构元素，对E_NDVI进行形态学膨胀得到宽的NDVI边界图像E^W _NDVI；

将在E_PAN中值为1但在E^W _NDVI中值为0的像元设置为0，得到PAN边界图像E′_PAN；

采用直径为L_SE2的圆盘结构元素对E′_PAN进行形态学膨胀得到宽的PAN边界图像E^W _PAN；

利用E^W _PAN，对MS中的亚像元进行分类，得到植被-非植被混合亚像元C_M和其他亚像元C_O。

4.根据权利要求3所述的基于混合亚像元解混的遥感多光谱与全色图像融合方法，其特征在于，利用E^W _PAN，对MS中的亚像元进行分类，得到植被-非植被混合亚像元C_M和其他亚像元C_O包括：

利用E^W _PAN，将MS中的亚像元分为两类，得到亚像元分类图，并根据亚像元分类图，得到植被-非植被混合亚像元C_M和其他亚像元C_O；

其中，亚像元分类图M_PI的分类规则如下：

其中，如果亚像元t在E^W _PAN中对应值为1，则被分为植被-非植被混合亚像元C_M；如果t在E^W _PAN中对应值为0，则被分为其它亚像元类C_O。

5.根据权利要求4所述的基于混合亚像元解混的遥感多光谱与全色图像融合方法，其特征在于，利用PAN和NDVI，判断植被-非植被混合亚像元的地物类型类别包括：

根据PAN边界图像E′_PAN，计算植被边界像元图像和非植被像元图像

利用直径为3的圆盘结构元素，对进行形态学膨胀，得到膨胀后的植被边界像元图像M_N，并将M_N中值为1且在中值为1的像元设置为0，并重复该步骤R-1次，使得M_N中值为1的像元均为植被像元；

利用直径为3的圆盘结构元素，对进行形态学膨胀，得到膨胀后的非植被边界像元图像M_NV，并将M_NV中值为1且在中值为1的像元设置为0，并重复该步骤R-1次，使得M_NV中值为1的像元均为非植被像元；

对于和中每个值为1的像元t₁，分别在以t₁为中心、大小为S_T×S_T的邻域窗口W(t₁,S_T)内计算所有值为1的像元的NDVI值，并根据该NDVI值，采用以下公式，分别得到图像和

对于M_PI中每个属于混合亚像元类C_M的像元t₂，根据M_N、M_NV、和采用以下公式判断其属于植被C_V或非植被C_NV类别，得到混合亚像元分类图M_LC：

6.根据权利要求5所述的混合亚像元解混的遥感多光谱与全色图像融合方法，其特征在于，根据PAN边界图像E′_PAN，计算植被边界像元图像和非植被像元图像包括：

对于每个亚像元t，在其8邻域内查找与其在PAN图像中灰度差异最大的邻居像元t’，得到像元对{t，t’}；

如果NDVI(t)＞NDVI(t’)，则且如果NDVI(t)＜NDVI(t’)，则且

7.根据权利要求6所述的混合亚像元解混的遥感多光谱与全色图像融合方法，其特征在于，对于MS中每个亚像元t，根据其对应的亚像元类别和地物类型类别进行融合包括：

采用暗像元法计算MS和PAN图像的雾气值H_i和H_p，

确定PAN_S中最小值对应的位置IDX_min，并利用该IDX_min和以下公式计算H_i和H_p：

H_i＝MS_i(IDX_min)

H_p＝PAN_S(IDX_min)

其中，i等于1，……N，N为MS图像的波段数量；

判断亚像元t是否属于其它亚像元类C_O，并在判断结果为是的情况下，利用该亚像元t的原始光谱，采用以下公式进行融合：

其中，MS_i为MS的第i波段，为融合图像MS^f的第i波段；

判断亚像元t是否属于植被-非植被混合亚像元C_M，并在判断结果为是的情况下，利用其邻域内的纯植被或纯非植被像元n的光谱，采用以下公式进行融合：

其中，邻域内的纯植被或纯非植被像元n在以t为中心、大小为S_N×S_N的窗口W(t,S_N)内根据以下公式查找：

其中，S_N的值小于等于L_SE2。