CN103186893A - 一种普适的高分辨率遥感图像融合方法 - Google Patents

Abstract

一种普适的高分辨率遥感图像融合方法。它在利用多光谱影像基于辐射亮度对全色高分辨率影像模拟的基础上，通过线性回归分析，建立模拟全色影像与原始全色影像的信息量特征之间的联系。同时针对全色与多光谱波段波谱范围不一致的情况提出了解决方法。

Description

一种普适的高分辨率遥感图像融合方法

技术领域

本发明属遥感技术领域，涉及遥感数据的处理及应用。 

背景技术

遥感影像融合能够集成或者整合多个源图像中的互补信息，利用优势互补的数据来提高图像的信息可用程度，同时增加对研究对象解译的可靠性。遥感影像数据融合作为一种遥感信息综合处理与分析技术，是遥感领域的热点研究课题之一。数据融合可以在以下三个层次进行：1)像素层融合(Pixel-Level Fusion)；2)特征层融合(Feature-level Fusion)；3)决策层融合(Decision-level Fusion)。像素层融合是目前三个层次中研究最为成熟的一级，已经形成了各种各样的融合算法。 

合成变量比方法是一种基于像素层的遥感影像融合算法。合成变量比(SVR)是由Munechika et al.(1993)(C.K.Munechika，J.S.Warnick，C.Salvaggio and J.R.Schott.Resolution enhancement of multispectral image data to improve classificationaccuracy[J].Photogrammetric Engineering and Remote Sensing，1993(59)：67-72.)提出的。SVR融合方法一般用于同一传感器的多光谱与全色数据融合，其计算公式为： 

${\mathrm{XSP}}_i={\mathrm{Pan}}_H\frac{{\mathrm{XS}}_{\mathrm{Li}}}{{\mathrm{Pan}}_{\mathrm{LSyn}}}---\left(1\right)$

其中XSP_i是融合后的高分辨率图像第i个波段的灰度值，Pan_H是原始全色图像的灰度值，XS_Li是原始多光谱图像中第i个波段的灰度值，Pan_LSyn是由低空间分辨率的多光谱波段模拟合成的低空间分辨率全色图像的灰度值。 

在Munechika et al.(1993)的研究中，模拟合成的全色图像Pan_LSyn是由公式(2)获得 

参数

是通过一种改进的大气模型来对接收到的目标反射和相对光谱响应进行修正，并测量了5种具有代表性的土地覆盖反射光谱，包括城市、土壤、水体、树和草。通过对使用大气模型模拟得到的值与实际测到的五类值的线性回归分析计算出了

在建立了Pan_LSyn之后，将原始全色图像相对于Pan_LSyn进行了直方图匹配，以消除大气和辐射角度的差异。 

Zhang(1999)和Zhang(2001)(Y.Zhang.A new merging method and its spectral andspatial effects[J].International Journal of Remote Sensing，1999(20)：2003～2014.Y.Zhang.Detection of urban housing development by fusing multi-sensor satellite data and performing spatial feature post-classification[J].International Journalof Remote Sensing，2001(22)：3339～3355.)认为，Munechika et al(1993)的SVR方法的计算中，只使用到几种特定的土地覆盖类型。而且系数

很难确定，不适用于实际应用中的大面积区域融合，同时在一个较大的城市区域中，5种土地覆盖类型远远不足以具有代表性，不能保证其它土地覆盖类型的光谱属性在融合之后得到保留。此外还要将原始全色影像相对与模拟的全色影像进行直方图匹配。 

为了解决以上问题，Zhang(1999)和Zhang(2001)将SVR方法的公式修改为： 

${\mathrm{XSP}}_i={\mathrm{Pan}}_H\frac{{\mathrm{XS}}_{\mathrm{Hi}}}{{\mathrm{Pan}}_{\mathrm{HSyn}}}---\left(3\right)$

其中，XSP_i和Pan_H不变，XS_Hi是将多光谱图像重采样为全色图像相同像元大小后，第i个波段的灰度值。Pan_HSyn是由下式模拟合成的高空间分辨率全色图像的灰度值： 

参数

是直接利用原始全色图像Pan_H与原始多光谱图像XS_Hi(具有与Pan_H相同的像元大小)进行多元回归分析计算得到。 

这里对SVR方法所做的修改包括： 

第一，参数

是直接从采样到和全色图像相同分辨率的多光谱图像和全色图像Pan_H中计算得到的，并未测量各种土地覆盖类型。模拟的全色图像的Pan_HSyn具有全色影像的高分辨率，同时与全色影像具有很高的相关性。节省了大量的测量时间和复杂的计算过程，因为不需要使用者测量不同的土地覆盖类型，参数更稳定，且不需要对全色影像与Pan_HSyn图像进行直方图调整。 

第二，模拟合成全色影像Pan_HSyn是从融合过程中将要使用的多光谱波段中计算出来的，而不是固定的某些波段。因而，任何一个波段都可以进行对全色图像Pan_HSyn的模拟合成，从而可以降低多光谱波段的颜色扭曲。 

曹鑫(2005)(曹鑫.ISVR遥感数据融合方法及其应用[D].北京师范大学，2005.)认为Zhang(1999)和Zhang(2001)对SVR的改进存在一定问题，即： 

(1)Zhang(1999)和Zhang(2001)对SVR的改进，只有当相关性越弱的时候，才能保留更多的光谱信息，如此根据回归分析方法计算出来的回归系数值得怀疑。 

(2)可扩展性差，若进行多幅图像的融合时，要对所有图像分别进行回归分析，计算模 拟将使用到的转换系数。 

因此，曹鑫对Zhang(1999)和Zhang(2001)的SVR方法进行了改进，提出了改进的合成变量比方法(ISVR)。SVR方法和Zhang(1999)和Zhang(2001)的SVR方法的核心问题是由低空间分辨率多光谱波段模拟合成低空间分辨率的全色波段，即Pan_HSyn。ISVR方法认为，在辐射亮度(radiance)条件下，模拟合成低空间分辨率的全色波段更容易实现。 

辐射亮度Radiance是表示传感器接收到的反射能量的物理量。ISVR最基本的假设是：模拟的全色波段的radiance可以通过多光谱波段(MS)来获得，MS波段和全色波段(PAN)覆盖的光谱范围基本一致。基于以上假设，模拟的全色波段的radiance可以认为是在相应波段上radiance的曲线积分。若将radiance曲线的积分部分分解开，那么模拟的全色波段的radiance将等于所有MS波段的radiance之和，再加上MS波段之间的间隔区域的面积。如图1所示，灰色区域为已知MS的radiance，如果近似地将空缺波长范围内的radiance曲线看作直线的话，那么白色梯形区域的面积可以近似认为是空缺波段的radiance。所有矩形和梯形的面积之和就是PAN波段的模拟radiance值。 

对于多光谱的1、2、3、4波段，其radiance分别可以看作图中矩形的面积(或者看作是radiance曲线在相应波长范围内的积分)，则在全色波段范围内，radiance对应的就是整个PAN波段radiance曲线的积分。如果近似地将空缺波长范围内的radiance曲线看作直线的话，那么只需要计算空白部分梯形的面积即可推算出整个PAN波段的radiance。 

图1ISVR方法中模拟全色波段radiance的原理图示 

ISVR方法的主要计算方法与Zhang(1999)相同，但是其中XSP_i、Pan_H和XS_Hi都是radiance而不是DN值。依图1所示，XS_H1，XS_H2，XS_H3和XS_H4分别表示多光谱影像四个波段的辐射亮度，即各波谱段所对应的矩形面积；Pan_HSyn的计算方法由公式(5)给出： 

${\mathrm{Pan}}_{\mathrm{HSyn}}={\mathrm{XS}}_{H1}+{\mathrm{XS}}_{H2}+...+{\mathrm{XS}}_{\mathrm{Hn}}+(\frac{{\mathrm{XS}}_{H1}}{{\mathrm{\λ}}_{\mathrm{1,2}}-{\mathrm{\λ}}_{1,1}}+\frac{{\mathrm{XS}}_{H2}}{{\mathrm{\λ}}_{\mathrm{2,2}}-{\mathrm{\λ}}_{\mathrm{2,1}}})\×\frac{{\mathrm{\λ}}_{\mathrm{2,1}}-{\mathrm{\λ}}_{\mathrm{1,2}}}{2}---\left(5\right)$

$+(\frac{{\mathrm{XS}}_{H2}}{{\mathrm{\λ}}_{\mathrm{2,2}}-{\mathrm{\λ}}_{\mathrm{2,1}}}+\frac{{\mathrm{XS}}_{H3}}{{\mathrm{\λ}}_{\mathrm{3,2}}-{\mathrm{\λ}}_{\mathrm{3,1}}})\×\frac{{\mathrm{\λ}}_{\mathrm{3,1}}-{\mathrm{\λ}}_{\mathrm{2,2}}}{2}+...+(\frac{{\mathrm{XS}}_{H(n-1)}}{{\mathrm{\λ}}_{(n-1),2}-{\mathrm{\λ}}_{(n-1),1}}+\frac{{\mathrm{XS}}_{\mathrm{Hn}}}{{\mathrm{\λ}}_{n,2}-{\mathrm{\λ}}_{n,1}})\×\frac{{\mathrm{\λ}}_{n,1}-{\mathrm{\λ}}_{(n-1),2}}{2}$

也可以写成： 

其中参数

由下式获得： 

其中，λ_i，1和λ_i，2是多光谱图像的第i个波段的两端波长数值，n是多光谱总的波段数目。 

相对于Zhang(1999)的SVR方法，ISVR的计算比较简单，模拟全色图像时所使用到的系数只与可见光/近红外各个通道的波长值有关，只要知道各个通道的波长范围，就可以利用公式简单地计算出模拟系数。这样的模拟系数具有普适性。对于相同的传感器，不需要在多次应用中反复计算系数，也不需要用户具备相当的经验。系数同时具有可解释性，具有一定的物理意义。 

曹鑫(2005)提出的ISVR融合方法在一定程度上解决了Zhang提出的改进的SVR融合方法中的问题，并通过实验，验证了ISVR有着较好的融合效果，唯一不足之处在于其理论机制限制了ISVR在使用MS波段模拟PAN波段图像时，只能局限于使用PAN波段波长内的MS波段。一般情况下，PAN图像的波段都是从可见光到近红外(VNIR)，因此，ISVR融合算法只适合于多光谱波段的波长范围在可见光和近红外的图像。 

与Zhang(1999)提出的对SVR方法的改进相比，ISVR方法赋予了SVR方法中的系数一定的物理意义。但是在ISVR方法中，只是基于辐射亮度，简单地通过多光谱影像对全色影像进行模拟，而且模拟的假设前提是，全色影像的波谱范围与多光谱的波谱范围基本一致，模拟的全色波段的辐射亮度可以认为是相应波段上辐射亮度的曲线积分，同时将空缺波长范围内的辐射亮度曲线近似看作是直线，用所有矩形和梯形的面积之和来模拟PAN波段的模拟radiance值。 

融合技术存在明显颜色差异，不能保证较好光谱保持性的主要原因是PAN波段波谱范围的改变，利用多光谱低空间分辨率影像对高分辨率波段进行模拟，应当尽量保证模拟的高分辨率波段影像信息量特性与高分辨率全色波段影像的信息量特性比较接近。在ISVR的方法中提到的全色波段的模拟方法中，仅仅是利用了多光谱影像和全色影像光谱范围之间的相关性，二者之间并未建立任何影像信息特性的联系， 

发明内容

本发明通过对多种融合方法原理的综合理解和融合试验，提出一种基于SVR的普适高分辨率遥感图像融合方法。它在利用多光谱影像基于辐射亮度radiance对全色高分辨率影像模拟的基础上，通过线性回归分析，建立模拟全色影像Pan_HSyn与原始PAN影像的信息量特征之间的联系。同时针对全色与多光谱波段波谱范围不一致的情况提出了解决方法。 

根据ISVR算法的假设前提，该方法物理意义上的可行性要求全色PAN波段的波谱范围和多光谱波段的波谱范围相一致。然而对于一些新型高分辨率遥感卫星影像，比如ALOS，全色波段的波谱范围是0.52-0.77μm，多光谱波段的波谱范围分别是(0.42-0.50μm、0.52-0.60μm、0.61-0.69μm、0.76-0.89μm)，显然全色影像的波谱范围和多光谱波谱范围并不一致，只是多光谱范围的一部分。因此，针对不同的新型高分辨率遥感影像的波谱设置特点，如何弥补全色影像和多光谱影像光谱范围的差异和缺失，改善影像的融合效果，使得融合方法具有更准确的物理意义，值得进一步思考和研究。 

针对这个问题，本发明尝试用波谱范围比值来修正最后得到的Pan_HSyn，以弥补全色与多光谱数据波谱范围的差异。 

具体实施方式

本发明融合算法的步骤如下： 

第一步、由全色影像波谱范围与多光谱影像波谱范围之间的关系和公式5-7计算出基于辐射亮度radiance的、由多光谱波段影像模拟的全色影像Pan_HSyn； 

第二步、用波谱范围比值来修正Pan_HSyn，以弥补全色与多光谱数据波谱范围的差异，即： 

设： 

$p=\frac{d_{{\mathrm{\λ}}_{\mathrm{pan}}}}{d_{\mathrm{\λmul}}}---\left(8\right)$

其中，p为影像波谱范围比，

为全色波段的波谱范围，d_λmul为多光谱波段的波谱范围； 当全色影像的波谱范围和多光谱影像的波谱范围相一致的时候，p＝1； 

令 

${\mathrm{Pan}}_{\mathrm{HSyn}}^{*}={\mathrm{Pan}}_{\mathrm{HSyn}}*p---\left(9\right)$

第三步、在原始全色波段影像和

之间寻找一一对应的像素点，这些点分别代表图像上的几种不同类型地物，例如：植被、水体、建筑物、裸地等；每一种地物类别分别取5-10个点作为样本点； 

第四步、计算原始全色波段影像和

样本点之间的线性回归系数k、b； 

第五步、由公式 

${\mathrm{Pan}}_{\mathrm{HSyn}}^{**}={\mathrm{kPan}}_{\mathrm{HSyn}}^{*}+b---\left(10\right)$

计算得到新的模拟全色的值

第六步、将

代入到式11中 

${\mathrm{XSP}}_i={\mathrm{Pan}}_H\frac{{\mathrm{XS}}_{\mathrm{Hi}}}{{\mathrm{Pan}}_{\mathrm{HSyn}}^{**}}---\left(11\right)$

其中XSP_i为融合后的高分辨率图像第i个波段的灰度值；Pan_H是全色图像，XS_Hi是将多光谱图像重采样为全色图像相同像元大小后，第i个波段的灰度值。 

从数学和物理角度分析，本发明融合方法一方面通过比值运算突出了高分辨率影像的结构纹理信息，同时与原始高分辨率全色波段相乘，在原始多光谱影像中引入了高分辨率影像的结构特征，大大提高了多光谱波段数据的空间分辨率；另一方面，由于

与原始PAN的相关性，光谱信息主要是多光谱波段变换得来的，融合之后的影像具有较好的光谱保持性。另外，该算法也适合于全色影像的波谱范围与多光谱的波谱范围不一致的情况，是一种普适性强的融合技术。 

融合实验与比较 

以ALOS高分辨率卫星影像为数据源，进行原始全色、多光谱、传统Brovey融合、ISVR融合以及本发明方法(简称IPS融合)的比较。 

以下列出上述几种融合结果的主要客观评价指标统计值，包括：灰度均值、标准差、信息熵和梯度值。 

表1ALOS不同融合算法灰度均值对比 

均值	蓝波段	绿波段	红波段	近红外波段
					多光谱	74.387	59.500	47.946	42.752
Brovey		21.376	17.239	15.249
					ISVR	17.842	14.282	11.524	10.216
IPS	73.628	58.906	47.496	42.134

表2ALOS不同融合算法标准差对比 

标准差	蓝波段	绿波段	红波段	近红外波段
					多光谱	4.012	5.025	6.174	7.513
Brovey		2.475	2.676	2.322
					ISVR	1.619	1.703	1.862	1.722
IPS	5.914	6.222	6.973	6.811

表3ALOS不同融合算法熵值对比 

熵	蓝波段	绿波段	红波段	近红外波段
					多光谱	3.903	4.244	4.604	4.844
Brovey		3.240	3.394	3.205
					ISVR	2.586	2.699	2.874	2.792
IPS	4.404	4.530	4.755	4.741

表4ALOS不同融合算法梯度值对比 

从表1-表4列出的各个客观评价指标统计值都可以看出，IPS融合法取得了很好的融合效果，IPS融合方法的灰度均值与原始多光谱影像最为接近，标准差和熵值都得到了很大的提高，某种程度上说明了其融合影像信息含量得到极大的丰富，并具有良好的光谱保持性。同时，每一种算法的融合影像的梯度值都显著提升，说明融合影像的纹理细节信息较原始多光谱影像大大增加，空间分辨率明显提高。 

ISVR融合方法，除了一些波段的梯度指标有所提高外，波段的灰度均值较原始多光谱影像差距很大，说明了ISVR融合方法虽然提高了影像的空间分辨率，但光谱保真性较差。 

从以上试验结果可以看出，经过线性回归分析和波谱范围修正的IPS方法，很好地改善了原ISVR融合算法光谱保持性差的问题，不仅显著提高了影像的空间分辨率，增强了纹理细节信息，还很好的保持了原始多光谱影像的光谱信息。 

Claims (1)

1.一种普适的高分辨率遥感图像融合方法，其步骤为：

第一步、由全色影像波谱范围与多光谱影像波谱范围之间的关系和公式1-3计算出基于辐射亮度radiance的、由多光谱波段影像模拟的全色影像Pan_HSyn；

${\mathrm{Pan}}_{\mathrm{HSyn}}={\mathrm{XS}}_{H1}+{\mathrm{XS}}_{H2}+...+{\mathrm{XS}}_{\mathrm{Hn}}+(\frac{{\mathrm{XS}}_{H1}}{{\mathrm{\λ}}_{\mathrm{1,2}}-{\mathrm{\λ}}_{1,1}}+\frac{{\mathrm{XS}}_{H2}}{{\mathrm{\λ}}_{\mathrm{2,2}}-{\mathrm{\λ}}_{\mathrm{2,1}}})\×\frac{{\mathrm{\λ}}_{\mathrm{2,1}}-{\mathrm{\λ}}_{\mathrm{1,2}}}{2}---\left(1\right)$

$+(\frac{{\mathrm{XS}}_{H2}}{{\mathrm{\λ}}_{\mathrm{2,2}}-{\mathrm{\λ}}_{\mathrm{2,1}}}+\frac{{\mathrm{XS}}_{H3}}{{\mathrm{\λ}}_{\mathrm{3,2}}-{\mathrm{\λ}}_{\mathrm{3,1}}})\×\frac{{\mathrm{\λ}}_{\mathrm{3,1}}-{\mathrm{\λ}}_{\mathrm{2,2}}}{2}+...+(\frac{{\mathrm{XS}}_{H(n-1)}}{{\mathrm{\λ}}_{(n-1),2}-{\mathrm{\λ}}_{(n-1),1}}+\frac{{\mathrm{XS}}_{\mathrm{Hn}}}{{\mathrm{\λ}}_{n,2}-{\mathrm{\λ}}_{n,1}})\×\frac{{\mathrm{\λ}}_{n,1}-{\mathrm{\λ}}_{(n-1),2}}{2}$

也可以写成：

其中参数由下式获得：

其中，XS_Hi表示多光谱影像第i个波段的辐射亮度；λ_i，1和λ_i，2是多光谱影像的第i个波段的两端波长数值，n是多光谱总的波段数目；

第二步、用波谱范围比值来修正Pan_HSyn，以弥补全色与多光谱数据波谱范围的差异，即：

设：

$p=\frac{d_{{\mathrm{\λ}}_{\mathrm{pan}}}}{d_{\mathrm{\λmul}}}---\left(4\right)$

其中，p为影像波谱范围比，

为全色波段的波谱范围，d_λmul为多光谱波段的波谱范围；

当全色影像的波谱范围和多光谱影像的波谱范围相一致的时候，p＝1；

令

${\mathrm{Pan}}_{\mathrm{HSyn}}^{*}={\mathrm{Pan}}_{\mathrm{HSyn}}*p---\left(5\right)$

第三步、在原始全色波段影像和之间寻找一一对应的像素点，这些点分别代表图像上的几种不同类型地物，例如：植被、水体、建筑物、裸地等；每一种地物类别分别取5-10个点作为样本点；

第四步、计算原始全色波段影像和

样本点之间的线性回归系数k、b；

第五步、由公式

${\mathrm{Pan}}_{\mathrm{HSyn}}^{**}={\mathrm{kPan}}_{\mathrm{HSyn}}^{*}+b---\left(6\right)$

计算得到新的模拟全色的值

第六步、将

代入到式7中

${\mathrm{XSP}}_i={\mathrm{Pan}}_H\frac{{\mathrm{XS}}_{\mathrm{Hi}}}{{\mathrm{Pan}}_{\mathrm{HSyn}}^{**}}---\left(7\right)$

其中XSP_i为融合后的高分辨率图像第i个波段的灰度值；Pan_H是全色图像，XS_Hi是将多光谱图像重采样为全色图像相同像元大小后，第i个波段的灰度值。

Images