CN104809691A - 一种基于滑动窗口混合像元分解的图像融合方法 - Google Patents

Abstract

本发明公开了一种基于滑动窗口混合像元分解的图像融合方法，包括如下步骤：A：对Hspec影像和Hspa影像进行预处理；B:对Hspec和Hspa影像，进行波段的相关性分析；C:对分析后的Hspa影像在每个滑动窗口中进行非监督分类；D:对预处理的Hspec影像进行混合像元分解；E:对预处理的Hspec影像进行空间插值；F:对所述D、E处理后的影像，生成Hspec和Hspa融合的遥感影像。本发明针对内陆浑浊水体，通过滑动窗口分类，进而建立线性混合像元分解的图像融合方法，将低空间分辨率高光谱分辨率数据与高空间分辨率低光谱分辨率数据融合，生成具有较高空间分辨率和光谱分辨率的数据，并在融合结果中，尽可能多地保留高光谱信息。

Description

一种基于滑动窗口混合像元分解的图像融合方法

技术领域

本发明涉及的是遥感技术领域，主要针对内陆水色遥感进行图像的预处理，具体涉及一种基于滑动窗口混合像元分解的图像融合方法。

背景技术

遥感技术可以快速地对水体的整体状况进行监测，具有传统人工实地监测无法比拟的优势，是水环境监测的重要技术手段。为了提高遥感数据对水环境监测的有效性，国际上已经发射了多颗水色遥感传感器，如SeaWiFS(Sea-Viewing Wide Field-of-View Sensor)、MERIS(MediumResolution Imaging Spectrometer)、GOCI(Geostationary Ocean ColorImager)、MODIS(Moderate Resolution Imaging Spectroradiometer)等。但是，这些传感器都是针对海洋水色遥感设计的，其空间分辨率往往从几百米到1千米，对于水域面积较小而光学特性复杂的内陆水体，其适用性受到了较大的限制；目前内陆水环境遥感使用的卫星数据源主要是空间分辨率相对较高的宽波段遥感数据，如搭载于美国陆地资源卫星的Thematic Mapper(TM)、法国的Systeme Probatoried’Observationdela Tarre(SPOT)、印度的Indian Space Research(IRS)-1C、AdvancedSpaceborn Thermal Emission and Reflection Radiometer(ASTER)、中巴资源卫星(China-brazil Earth Resource Satellite,CBERS)、我国的环境1号小卫星(HJ1)等。但是，这些数据较宽的波段宽度往往影响了水色参数反演的精度。因此，针对内陆水色遥感的特点，研究适宜的数据融合方法，生成既具有高光谱或水色遥感传感器的波段优势，又具有较高空间分辨率的遥感数据，对内陆水环境遥感具有重要的意义。

发明内容

发明目的：针对现有技术上存在的不足，本发明目的是在于提供一种基于滑动窗口混合像元分解的图像融合方法，针对内陆浑浊水体，通过滑动窗口分类，进而建立线性混合像元分解的图像融合方法，将低空间分辨率高光谱分辨率数据与高空间分辨率低光谱分辨率数据融合，生成具有较高空间分辨率和光谱分辨率的数据，并在融合结果中，尽可能多地保留高光谱信息。

技术方案：为了实现上述目的，本发明是通过如下的技术方案来实现：一种基于滑动窗口混合像元分解的图像融合方法，包括如下步骤：

1)对获取的低空间分辨率高光谱分辨率影像和高空间分辨率低光谱分辨率影像进行预处理；

2)将处理后的低空间分辨率高光谱分辨率影像和高空间分辨率低光谱分辨率影像进行波段的相关性分析，挑选出后续处理所需要的高空间分辨率低光谱分辨率影像数据；

3)对经过步骤2处理后的高空间分辨率低光谱分辨率影像在每个滑动窗口中进行非监督分类；

4)利用步骤3得到的分类结果，对步骤1处理后的低空间分辨率高光谱分辨率影像进行混合像元分解；

5)对步骤1处理后的低空间分辨率高光谱分辨率影像进行空间插值；

6)将步骤4处理后的影像与步骤5处理后的影像融合，生成最终的遥感影像。

优选的，所述的步骤1具体包括如下步骤：

1-1)获取成像时间相近的低空间分辨率高光谱分辨率影像和高空间分辨率低光谱分辨率影像；

1-2)对获取的影像进行辐射校正，将DN值转化为大气顶层辐亮度(L(λ))，L(λ)＝Gain×DN+Bias，其中，DN为图像的灰度值，无量纲；Gain和Bias分别为增益和偏置，单位为Wm^-2μm^-1sr^-1，增益和偏置的取值由传感器官方给定；

1-3)对步骤1-2处理的数据进行几何校正，几何校正采用常规的方法进行；

1-4)对步骤1-3处理的数据进行大气校正，大气校正采用6S大气校正模型、ENVI遥感软件嵌入的FLAASH大气校正模块进行。

优选的，所述的步骤2具体包括如下步骤：

2-1)针对步骤1中预处理后的数据，将低空间分辨率高光谱分辨率影像和高空间分辨率低光谱分辨率影像的波段进行两两匹配的相关性分析；

2-2)对比步骤2-1中的相关性分析，针对低空间分辨率高光谱分辨率影像的每个波段，挑选与之相关性最大的高空间分辨率低光谱分辨率影像波段，并将其作为后续处理的影像数据。

优选的，所述的步骤3具体包括如下步骤：

3-1)率定滑动窗口大小w：首先设置多个滑动窗口尺度，然后将每个窗口尺度方案一一代入算法中，利用相对尺度误差ERGAS对结果进行评价，以ERGAS逐步稳定达到最小值时的w值为最终确定的滑动窗口尺度大小值。ERGAS的计算公式如下：

${\mathrm{ERGAS}}_{\mathrm{scale}}=100\frac{S_h}{S_l}\sqrt{\frac{1}{N}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}({\mathrm{RMSE}}_i^2/M_i^2)}$

其中S_h代表高空间分辨率影像的空间分辨率；S_l代表低空间分辨率影像的空间分辨率；N为参与融合的波段数量；当计算低空间分辨率尺度的ERGAS时，RMSE_i代表在空间尺度退化之后的融合影像与低分辨率影像之间的均方根误差，M_i代表低分辨率影像的均值；当计算高空间分辨率尺度的ERGAS时，RMSE_i代表在光谱尺度退化之后的融合影像与高分辨率影像之间的均方根误差，M_i代表高分辨率影像的均值。

3-2)率定归一化因子α：首先设置多个w与α的组合，然后，将每种组合方案一一代入算法中，利用相对尺度误差ERGAS对结果进行评价，以ERGAS逐步稳定达到最小值时的w与α的组合为最终确定的α值的大小。

3-3)根据步骤3-1确定的w值，取窗口大小为w×w，设定分类类别数为w×w，利用步骤2挑选出的高空间分辨率低光谱分辨率影像波段，在每个滑动窗口中进行非监督分类。

优选的，所述的步骤4具体包括如下步骤：

4-1)利用所述步骤3的分类结果，计算步骤1处理的低空间分辨率高光谱分辨率影像的像元内每个类别的面积比例，即该类别的像元丰度，记为c_i(l,s；k)，其中i表示低空间分辨率高光谱分辨率影像第i个波段，k表示第k类别，(l,s)为低空间分辨率高光谱分辨率影像像元行列号；

4-2)利用计算出的丰度，求算低空间分辨率高光谱分辨率影像像元的最佳端元值，最佳端元值的计算通过优化目标函数J完成：

$J=\mathrm{\Σ}{(S_i(l,s)-\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{c}_i(l,s;k){\stackrel{\‾}{S}}_i\left(k\right))}^2+\mathrm{\α}\frac{N_{\mathrm{Hspec}}-1}{K}\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{({\stackrel{\‾}{S}}_i\left(k\right)-{{\stackrel{\‾}{S}}_i}^{\′}\left(k\right))}^2$

其中i为当前进行分解的第i个波段；S_i(l,s)为低空间分辨率高光谱分辨率影像在(l,s)位置的像元值；c_i(l,s；k)为类别k的像元丰度；α为归一化参数；K和N_Hspec分别是窗口中的类别数和低空间分辨率高光谱分辨率影像像元的个数；为窗口中的特定类别中低空间分辨率高光谱分辨率影像像元的均值；为未知的端元亮度值。

优选的，在步骤5中对步骤1处理的低空间分辨率高光谱分辨率影像进行重采样，使得其空间分辨率与高空间分辨率低光谱分辨率影像空间分辨率一致，重采样后的影像记为I(m,n)。

优选的，所述的步骤6具体包括如下步骤：

6-1)将所述步骤4计算的端元值与所述步骤5得到的重采样计算的影像融合计算生成最终的融合影像，表达式如下所示：

${F_i}^{*}(m,n)=\frac{K_c}{K}\·{\stackrel{\‾}{S}}_i^{*}\left(k(m,n)\right)+(1-\frac{K_c}{K})\·I(m,n)$

F_i ^*(m,n)为最终生成的窗口中心像元的像素值；和I(m,n)分别是所述步骤4计算得到的端元值与所述步骤5得到的插值影像，Kc为低空间分辨率高光谱分辨率影像窗口中心像元中在所述步骤3中分类的类别数；

6-2)以w×w窗口遍历所有的低空间分辨率高光谱分辨率影像和高空间分辨率低光谱分辨率影像，完成波段i中所有像元的计算。

6-3)对低空间分辨率高光谱分辨率影像上所有波段重复上述步骤进行计算，完成所有波段的图像融合计算。

有益效果：本发明的数据融合方法，为每个低空间分辨率影像波段，选取对应的高空间分辨率波段进行分类，提高了水体的分类精度；提出在每个滑动窗口中进行非监督分类的方法，避免了通常分类中容易造成水体“唯一端元效应”的缺陷；对低空间分辨率影像进行空间插值重采样，体现了水质参数连续变化的特点。综上，本发明特别适合于水体遥感影像的图像融合计算。

附图说明

图1为本发明实施例所述的基于滑动窗口混合像元分解的图像融合的方法流程图。

具体实施方式

为使本发明实现的技术手段、创作特征、达成目的与功效易于明白了解，下面结合具体实施方式，进一步阐述本发明。

本具体实施方式采用以下技术方案：基于滑动窗口混合像元分解的图像融合方法，其图像融合方法步骤：

1)对低空间分辨率高光谱分辨率影像(以下记为Hspec)和高空间分辨率低光谱分辨率影像(以下记为Hspa)进行预处理；

2)对所述步骤1处理后的Hspec和Hspa影像，进行波段的相关性分析；

3)对所述步骤2处理后的Hspa影像，在每个滑动窗口中进行非监督分类；

4)利用步骤3得到的分类结果，对步骤1处理后的Hspec影像进行混合像元分解；

5)对所述步骤1处理后的Hspec影像进行空间插值；

6)将步骤4处理后的影像与步骤5处理后的影像融合，生成Hspec和Hspa融合的遥感影像(记为Cspec-spa)。

下面利用一个具体实施例进行说明：

本实施例以2009年4月25日太湖HJ1-CCD数据、MERIS数据为示例，其中HJ1-CCD数据代表高空间分辨率低光谱分辨率数据(Hspec)，MERIS数据代表高光谱分辨率低空间分辨率数据(Hspa)；HJ1-CCD数据的空间分辨率为30m，光谱通道为4个，MERIS数据的空间分辨率为300m，光谱通道为15个，图1是本发明实施例所述的基于滑动窗口混合像元分解的图像融合的方法流程图，如图1所述，所述方法包括如下的步骤：

A：通过网站获取了4月25日太湖HJ1-CCD数据、MERIS数据,并对数据进行预处理，包括：辐射定标、几何校正、大气校正；

B:对所述A处理后的HJ1影像4个波段反射率与MERIS影像的15个波段反射率，进行两两相关分析，将对应MERIS影像某波段(如第i波段)相关性最大的波段(如第j波段)HJ1影像作为i波段MERIS影像分类的基准图像；

C:对所述A处理后的HJ1影像，经过所述B分析后，在j波段影像每个滑动窗口中进行非监督分类，类别数是窗口大小的平方。如窗口大小为w×w的话，就将窗口中高分辨率影像分类类别个数的上限设定为w×w。其中，需要率定参数w和α，具体步骤如下：

C1：率定滑动窗口大小w:首先设置多个滑动窗口尺度，然后，将不同方案代入算法中，利用相对尺度误差ERGAS对结果进行评价，以ERGAS逐步稳定达到最小值时的w值为最终确定的滑动窗口尺度大小值。ERGAS的计算公式如下：

${\mathrm{ERGAS}}_{\mathrm{scale}}=100\frac{S_h}{S_l}\sqrt{\frac{1}{N}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}({\mathrm{RMSE}}_i^2/M_i^2)}$

其中S_h代表HJ1影像的空间分辨率；S_l代表MERIS影像的空间分辨率；N为参与融合的波段数量；当计算低空间分辨率尺度的ERGAS时，RMSE_i代表在空间尺度退化之后的融合影像与MERIS影像之间的均方根误差，M_i代表MERIS影像的均值；当计算高空间分辨率尺度的ERGAS时，RMSE_i代表在光谱尺度退化之后的融合影像与HJ1影像之间的均方根误差，M_i代表HJ1影像的均值。

C2：率定归一化因子α:首先设置多个w与α的组合，然后，将不同方案代入算法中，利用相对尺度误差ERGAS对结果进行评价，以ERGAS逐步稳定达到最小值时的w与α的组合为最终确定的α值的大小。

C3：根据所述C1确定的w值，取窗口大小为w×w，设定分类类别数为w×w，利用所述B的HJ影像波段，在每个滑动窗口中进行非监督分类。

本实施例率定后的参数如下：w＝7，α＝0.001。

D:对所述A处理后的MERIS影像进行混合像元分解，具体步骤如下：

D1：利用所述C3的分类结果，计算所述A处理的MERIS影像像元内每个类别的面积比例，称为该类别的像元丰度，记为c_i(l,s；k)，其中i表示MERIS影像第i个波段，k表示第k类别，(l,s)为MERIS影像像元行列号；

D2：利用D1所述计算的丰度，求算MERIS影像像元的最佳端元值。最佳端元值的计算，通过优化目标函数J完成：

$J=\mathrm{\Σ}{(S_i(l,s)-\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{c}_i(l,s;k){\stackrel{\‾}{S}}_i\left(k\right))}^2+\mathrm{\α}\frac{N_{\mathrm{Hspec}}-1}{K}\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{({\stackrel{\‾}{S}}_i\left(k\right)-{{\stackrel{\‾}{S}}_i}^{\′}\left(k\right))}^2$

其中i为当前进行分解的第i个波段；S_i(l,s)为MERIS影像在(l,s)位置的像元值；c_i(l,s；k)为类别k的像元丰度；α为归一化参数；K和N_Hspec分别是窗口中的类别数和MERIS像元的个数；为窗口中的特定类别中MERIS像元的均值；为未知的端元亮度值。

E:对所述A处理后的MERIS影像进行空间插值；

F:对所述D、E处理后的影像，生成MERIS和HJ1融合的遥感影像，具体步骤如下。

F1：对所述D计算的端元值与所述E重采样计算的影像融合计算生成最终的融合影像。表达式如下所示：

${F_i}^{*}(m,n)=\frac{K_c}{K}\·{\stackrel{\‾}{S}}_i^{*}\left(k(m,n)\right)+(1-\frac{K_c}{K})\·I(m,n)$

F_i ^*(m,n)为最终生成的窗口中心像元的像素值；和I(m,n)分别是所述D计算得到的端元值与所述E得到的插值影像，Kc为Hspec影像窗口中心像元中在所述C中分类的类别数。

F2：以7×7窗口遍历MERIS和HJ1影像，完成波段i中所有像元的计算。

F3：对MERIS影像所有波段重复上述步骤进行计算，完成所有波段的图像融合计算。

本发明的数据融合方法，为每个低空间分辨率影像波段，选取对应的高空间分辨率波段进行分类，提高了水体的分类精度；发明中提出的在每个滑动窗口中进行非监督分类的方法，避免了通常分类中容易造成水体“唯一端元效应”的缺陷；对低空间分辨率影像进行空间插值重采样，体现了水质参数连续变化的特点。综上，本发明适合于水体遥感影像的图像融合计算，特别是针对内陆混浊水体，一方面，由于水质的空间差异较大，低空间分辨率遥感数据难以体现其空间差异；另一方面，由于水体是弱信号体，导致其光谱信息的细微差异又难以从低光谱分辨率遥感影像中反映出来，因此，本发明尽可能多地保留了高光谱信息，同时提高了空间分辨率，在水环境遥感的数据处理中具有广泛的应用前景。

Claims (7)

1.一种基于滑动窗口混合像元分解的图像融合方法，其特征在于：包括如下步骤：

1)对获取的低空间分辨率高光谱分辨率影像和高空间分辨率低光谱分辨率影像进行预处理；

2)将处理后的低空间分辨率高光谱分辨率影像和高空间分辨率低光谱分辨率影像进行波段的相关性分析，挑选出后续处理所需要的高空间分辨率低光谱分辨率影像数据；

3)对经过步骤2处理后的高空间分辨率低光谱分辨率影像在每个滑动窗口中进行非监督分类；

4)利用步骤3得到的分类结果，对步骤1处理后的低空间分辨率高光谱分辨率影像进行混合像元分解；

5)对步骤1处理后的低空间分辨率高光谱分辨率影像进行空间插值；

6)将步骤4处理后的影像与步骤5处理后的影像融合，生成最终的遥感影像。

2.根据权利要求1所述的基于滑动窗口混合像元分解的图像融合方法，其特征在于：所述的步骤1具体包括如下步骤：

1-1)获取成像时间相近的低空间分辨率高光谱分辨率影像和高空间分辨率低光谱分辨率影像；

1-2)对获取的影像进行辐射校正，将DN值转化为大气顶层辐亮度(L(λ))，L(λ)＝Gain×DN+Bias，其中，DN为图像的灰度值，无量纲；Gain和Bias分别为增益和偏置，单位为Wm^-2μm^-1sr^-1，增益和偏置的取值由传感器官方给定；

1-3)对步骤1-2处理的数据进行几何校正，几何校正采用常规的方法进行；

1-4)对步骤1-3处理的数据进行大气校正，大气校正采用6S大气校正模型、ENVI遥感软件嵌入的FLAASH大气校正模块进行。

3.根据权利要求1所述的基于滑动窗口混合像元分解的图像融合方法，其特征在于：所述的步骤2具体包括如下步骤：

2-1)针对步骤1中预处理后的数据，将低空间分辨率高光谱分辨率影像和高空间分辨率低光谱分辨率影像的波段进行两两匹配的相关性分析；

2-2)对比步骤2-1中的相关性分析，针对低空间分辨率高光谱分辨率影像的每个波段，挑选与之相关性最大的高空间分辨率低光谱分辨率影像波段，并将其作为后续处理的影像数据。

4.根据权利要求1所述的基于滑动窗口混合像元分解的图像融合方法，其特征在于：所述的步骤3具体包括如下步骤：

3-1)率定滑动窗口大小w：首先设置多个滑动窗口尺度，然后将每个窗口尺度方案一一代入算法中，利用相对尺度误差ERGAS对结果进行评价，以ERGAS逐步稳定达到最小值时的w值为最终确定的滑动窗口尺度大小值。ERGAS的计算公式如下：

${\mathrm{ERGAS}}_{\mathrm{scale}}=100\frac{S_h}{S_l}\sqrt{\frac{1}{N}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}(\mathrm{RMS}{E_i}^2/{M_i}^2)}$

其中S_h代表高空间分辨率影像的空间分辨率；S_l代表低空间分辨率影像的空间分辨率；N为参与融合的波段数量；当计算低空间分辨率尺度的ERGAS时，RMSE_i代表在空间尺度退化之后的融合影像与低分辨率影像之间的均方根误差，M_i代表低分辨率影像的均值；当计算高空间分辨率尺度的ERGAS时，RMSE_i代表在光谱尺度退化之后的融合影像与高分辨率影像之间的均方根误差，M_i代表高分辨率影像的均值。

3-2)率定归一化因子α：首先设置多个w与α的组合，然后，将每种组合方案一一代入算法中，利用相对尺度误差ERGAS对结果进行评价，以ERGAS逐步稳定达到最小值时的w与α的组合为最终确定的α值的大小。

3-3)根据步骤3-1确定的w值，取窗口大小为w×w，设定分类类别数为w×w，利用步骤2挑选出的高空间分辨率低光谱分辨率影像波段，在每个滑动窗口中进行非监督分类。

5.根据权利要求1所述的基于滑动窗口混合像元分解的图像融合方法，其特征在于：所述的步骤4具体包括如下步骤：

4-1)利用所述步骤3的分类结果，计算步骤1处理的低空间分辨率高光谱分辨率影像的像元内每个类别的面积比例，即该类别的像元丰度，记为c_i(l,s；k)，其中i表示低空间分辨率高光谱分辨率影像第i个波段，k表示第k类别，(l,s)为低空间分辨率高光谱分辨率影像像元行列号；

4-2)利用计算出的丰度，求算低空间分辨率高光谱分辨率影像像元的最佳端元值，最佳端元值的计算通过优化目标函数J完成：

$J=\mathrm{\Σ}{(S_i(l,s)-\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{c}_i(l,s;k){\stackrel{\text{-}}{S}}_i\left(k\right))}^2+\mathrm{\α}\frac{N_{\mathrm{Hspec}}-1}{K}\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{({\stackrel{\‾}{S}}_i\left(k\right)-{{\stackrel{\‾}{S}}_i}^{\′}\left(k\right))}^2$

其中i为当前进行分解的第i个波段；S_i(l,s)为低空间分辨率高光谱分辨率影像在(l,s)位置的像元值；c_i(l,s；k)为类别k的像元丰度；α为归一化参数；K和N_Hspec分别是窗口中的类别数和低空间分辨率高光谱分辨率影像像元的个数；为窗口中的特定类别中低空间分辨率高光谱分辨率影像像元的均值；为未知的端元亮度值。

6.根据权利要求1所述的基于滑动窗口混合像元分解的图像融合方法，其特征在于：在步骤5中对步骤1处理的低空间分辨率高光谱分辨率影像进行重采样，使得其空间分辨率与高空间分辨率低光谱分辨率影像空间分辨率一致，重采样后的影像记为I(m,n)。

7.根据权利要求1所述的基于滑动窗口混合像元分解的图像融合方法，其特征在于：所述的步骤6具体包括如下步骤：

6-1)将所述步骤4计算的端元值与所述步骤5得到的重采样计算的影像融合计算生成最终的融合影像，表达式如下所示：

${F_i}^{*}(m,n)=\frac{K_c}{K}\·{\stackrel{\‾}{S}}_i^{*}\left(k(m,n)\right)+(1-\frac{K_c}{K})\·I(m,n)$

F_i ^*(m,n)为最终生成的窗口中心像元的像素值；和I(m,n)分别是所述步骤4计算得到的端元值与所述步骤5得到的插值影像，Kc为低空间分辨率高光谱分辨率影像窗口中心像元中在所述步骤3中分类的类别数；

6-2)以w×w窗口遍历所有的低空间分辨率高光谱分辨率影像和高空间分辨率低光谱分辨率影像，完成波段i中所有像元的计算。

6-3)对低空间分辨率高光谱分辨率影像上所有波段重复上述步骤进行计算，完成所有波段的图像融合计算。