CN1296871C - 基于残差的遥感图像融合方法 - Google Patents

Abstract

本发明为一种基于残差的遥感图像融合方法。该方法首先使用残差抽取算法获取多光谱图像的残差图像和全色图像的残差图像；然后将多光谱图像的残差图像插值到和全色图像的残差图像同样的尺寸，通过主成分分析对插值后的多光谱残差图像和全色图像的残差图像进行融合，恢复出多光谱图像的高分辨率残差；最后将多光谱图像的插值图像和多光谱图像的残差图像相加得到多光谱图像的高分辨率图像估计。本发明在较大地增强融合图像的空间细节表现能力的同时，很好地保留了原多光谱图像的光谱信息，其物理意义明确、实现结构简单、融合效果好，可为改善遥感图像的目视判读精度，提高信息清晰度和可靠性方面提供有效的技术支持。

Description

基于残差的遥感图像融合方法

技术领域

本发明属于遥感图像处理技术领域，具体涉及一种基于残差的遥感图像融合方法。

背景技术

遥感图像融合就是将不同类型的图像数据进行空间配准后，采用一定方法将各图像的优点结合起来产生新图像的技术。例如，通常人们希望将多光谱图像和全色图像进行融合，使融合后的多光谱图像具有较高的空间细节表现能力，同时保持多光谱图像的光谱特性，以获得对地表景物更全面、清晰的描述，从而提高遥感图像的判读与分析精度。

遥感图像的融合方法目前主要有HIS(Hue-Intensity-Saturation)变换融合方法^[1，2]、PCA(Principal Component Analysis)融合方法^[3，4]和小波变换融合方法^[5-7]等，但这些方法各自有以下主要缺点。

HIS变换融合方法在图像处理中已经成为一个标准方法，由于采用全色图像替换多光谱图像的亮度分量，从而在增强多光谱图像的空间细节表现能力的同时，也带来了较大的光谱失真。PCA融合方法适用于多光谱图像的所有波段，由于采用全色图像替换多光谱图像的第一主成分，在增强了融合图像的空间细节表现能力的同时，也使得光谱分辨率受到很大的影响。这两种“成分替换”方法的共同缺点在于多光谱图像的光谱信息会遭受比较严重的破坏^[1-4]。

小波变换融合方法由于直接舍弃了全色图像的低频分量，因此容易出现分块效应。同时，小波变换的分解阶数的大小也对融合效果产生很大的影响。如果小波分解阶数选得较低，虽然光谱特性保留较好，但增强后多光谱图像的空间细节表现能力较差，容易出现地物纹理模糊；若小波分解的阶数选得较高，虽然增强后多光谱图像的空间细节表现能力较好，但光谱特性保留较差^[5-7]。

针对以上问题，在遥感图像融合的研究中，如何在较大地增强融合图像的空间细节表现能力的同时，很好地保留原多光谱图像的光谱信息成为目前研究的热点。

                                 参考文献

[1]Carper W.J.，Lillesand T.M.，Kiefer R.W.，The use of the intensity-hue-saturation transformations formerging SPOT panchromatic and multispectral image data[J]，Photogrammetric Engineering and RemoteSensing，1990，56：459-467

[2]R.Welch，M.Ehlers，Merging multi-resolution SPOT HRV and Landsat TM data，PhotogrammetricEngineering and Remote Sensing，1987，53：301-303

[3]P.S.Chavez，J.Stuart，C.Slides，J.A.Anderson，Comparison of three different methods to mergemulti-resolution and multi-spectral data：Landsat TM and SPOT panchromatic[J]，PhotogrammetricEngineering and Remote Sensing，1991，57：295-303

[4]V.K.Shettigara，A generalized component substitution technique for spatial enhancement of multispectralimages using a higher resolution data set[J]，Photogrammetric Engineering and Remote Sensing，1992，58：561-567

[5]Nunez J.，Otazu X.，Fors O.，Prades A.，Pala V.，Arbiol R.，Multiresolution-based image fusion withadditive wavelet decomposition[J]，IEEE Transaction on Geoscience and Remote Sensing，1999，37(3)：1204-1211

[6]Gonzalez-Audicana M.，Saleta J.L.，Catalan R.G.，Garcia R.，Fusion of multi-spectral and panchromaticimages using improved IHS and PCA mergers based on wavelet decomposition[J]，IEEE Transaction onGeoscience and Remote Sensing，2004，42(6)：1291-1299

[7]Li Jun，Zhou YueQin，Li DeRen，Fusion of high-resolution panchromatic and multispectral images by usingwavelet transform[J]，Journal of remote sensing，1999，3(2)：116-122

[8]Li BiCheng，Wei Jun，Peng TianQiang，Objective analysis and evaluation of remote sensing image fusioneffect[J]，Computer engineering & science(李弼程，魏俊，彭天强，遥感图像融合效果的客观分析和评价，计算机工程与科学)，2004，26(1)：42-46

发明内容

本发明的目的是提出一种基于残差的遥感图像融合方法，使融合图像具有较高的空间细节表现能力，又尽可能多的保持多光谱图像的光谱信息。

本发明提出的基于残差的遥感图像融合方法，具体步骤如下：首先使用残差抽取算法抽取多光谱图像的残差图像和全色图像的残差图像；然后将多光谱图像的残差图像插值到和全色图像的残差图像同样的尺寸，通过主成分分析对插值后的多光谱残差图像和全色图像的残差图像进行融合，恢复出多光谱图像的高分辨率残差；最后将多光谱图像的插值图像和多光谱图像的高分辨率残差图像相加，得到多光谱图像的高分辨率图像估计。

本发明中，所述的残差抽取算法E(·)如附图2所示。残差抽取算法首先通过图像插值算法获得其插值图像，然后对该插值图像通过低通滤波和降采样过程获得原图像的低分辨率图像估计，最后用原图像与低分辨率的图像估计相减得到图像的残差图像。通过上述残差抽取算法E(□)可以分别获取多光谱图像和全色图像的残差图像。

本发明中，所述的插值方法Z(□)选用双三次插值方法，通过待插值点附近的4×4个像素来确定待插值像素的灰度。通过该插值方法将多光谱图像的残差图像放大到和全色图像的残差图像同样的大小。

本发明中，所述的基于主成分分析的残差融合方法，首先进行如下所示的主成分分析：

(1)计算参与融合的多光谱图像n个波段灰度图像的协方差矩阵S_x。

(2)算出S_x的全部特征值λ₁，λ₂，…，λ_n和对应的特征向量u¹，u²，…，uⁿ，其中λ₁≥λ₂≥…≥λ_n。

(3)求出n个新特征向量Y＝(PC¹，PC²，…，PCⁿ)^T，它们满足Y＝U^TX，其中U＝(u¹，u²，…，uⁿ)^T，而且S_y＝Λ＝diag{λ₁，λ₂，…，λ _n}。

此时的PC¹，PC²，…，PCⁿ分别被称为第1，2，…，n个主成分，它的逆变换为X＝UY。

这里我们用全色图像的残差图像e_f ^p经过灰度调整后替换第一主成分PC¹，以此来增强残差图像e_f的空间细节。由于残差图像e_f ^p与第一主成分图像PC¹的均值均为零，其灰度调整可如式(2)进行，

        z(i，j)＝x(i，j)×σ_y/σ_x，

其中，x(i，j)表示在图像e_f ^p在点(i，j)处的灰度值，σ_x ²和σ_y ²分别表示图像e_f ^p和第一主成分PC¹的方差，z(i，j)是替换后的第一主成分图像PC^1′在点(i，j)处的灰度值。通过对替换后的第一主成分图像PC^1′和其他主成分进行PCA的逆变换可以获得高分辨率的残差图像估计e′_f。

本发明中，所述的叠加方法，是把多光谱图像的插值图像I和多光谱图像的高分辨率残差图像e′_f逐像素相加得到多光谱图像的高分辨率图像估计MS′，如式(2)所示：

        MS′＝I+e′_f。

在以上的图像融合方法中，我们通过多光谱图像的插值图像I使得融合后的图像保留了原多光谱图像中主要的光谱信息。

本发明所提出的基于残差的遥感图像融合方法，有效地解决了遥感图像融合中增强空间细节信息和保持光谱信息的矛盾。为了提高融合的综合效果，可以在残差的层次上进行主成分分析的融合。

附图说明

图1为基于残差的遥感图像融合方法的流程图示。

图2为残差抽取的流程图示。

图3为不同融合方法结果的比较。其中，图3(a)为多光谱彩色图像，图3(b)为全色图像，图3(c)为基于残差的融合结果，图3(d)为HIS变换的融合结果，图3(e)为PCA的融合结果，图3(f)为小波变换的融合结果。

具体实施方式

以下通过实例，进一步对发明中的各个组成加以描述

1.残差的抽取

残差抽取算法E(□)如图2所示。首先通过图像插值算法Z(□)获得其插值图像I，然后通过低通滤波和降采样过程L(□)获得原多光谱图像的低分辨率图像估计MS_L，最后用原多光谱图像MS与低分辨率的图像估计MS_L相减得到多光谱图像的残差图像e_g。这里的插值算法Z(□)选用双三次插值方法，通过待插值点附近的4×4个像素来确定待插值像素的灰度。采用的低通滤波器是11×11的均值低通滤波器。对于全色图像，可以采用同样的残差抽取算法E(□)去获取全色图像的残差图像e_f ^p。

为了将多光谱图像的残差图像e_g放大到与全色图像的残差图像e_f ^p一样大，这里我们利用双三次插值方法来放大多光谱图像的残差图像，如下式所示

e_f＝Z(e_g).                                           (3)

2.残差的融合

经过插值后得到的残差图像e_f需要增强空间细节信息，本发明采用主成分分析的方法在多光谱图像的残差图像中融入全色图像的残差图像的空间细节信息。

主成分分析的步骤如下：

(1)计算参与融合的多光谱图像3个波段灰度图像的协方差矩阵S_x。

(2)算出S_x的全部特征值λ₁，λ₂，λ₃和对应的特征向量u¹，u²，u³，其中λ₁≥λ₂≥λ₃。

(3)求出3个新特征向量Y＝(PC¹，PC²，PC³)^T，它们满足Y＝U^TX，其中U＝(u¹，u²，u³)^T，而

且S_y＝Λ＝diag{λ₁，λ₂，λ₃}。

此时的PC¹，PC²，PC³分别被称为第1，2，3个主成分，它的逆变换为X＝UY。其中PC¹表示e_f的第一主成分。

用全色图像的残差图像e_f ^p经过灰度调整后替换第一主成分PC¹，以此来增强多光谱残差图像的空间细节。由于残差图像e_f ^p与第一主成分图像PC¹的均值均为零，其灰度调整可如式(1)进行，

z(i，j)＝x(i，j)×σ_y/σ_x

其中，x(i，j)表示在图像e_f ^p在点(i，j)处的灰度值，σ_x ²和σ_y ²分别表示图像e_f ^p和第一主成分PC¹的方差，z(i，j)是替换后的第一主成分图像PC^1′在点(i，j)处的灰度值。通过对替换后的第一主成分PC^1′和其他主成分进行PCA的逆变换得到融合后的多光谱图像的残差图像

最后，把多光谱图像的插值图像I和多光谱图像的高分辨率残差图像 相叠加得到多光谱图像的高分辨率图像估计MS′，如式(2)所示

${\mathrm{MS}}^{\′}=I+{\tilde{e}}_f.$

对本发明方法，进行了仿真计算。具体的仿真条件如下：

以2000年7月14日Landsat ETM+传感器拍摄到的上海地区多光谱图像和全色图像的融合为例进行说明。以Landsat ETM+第3、2和1波段分别作为Red、Green和Blue通道进行彩色图像的合成，其结果如图3(a)所示，其空间分辨率为30m，图像尺寸是256×256像素。另外，对应于同一区域的全色图像如图3(b)所示，其空间分辨率为15m，图像的尺寸是512×512像素。

为了融合图3(a)和(b)所示的多光谱图像和全色图像，采用了本发明的基于残差的图像融合方法，其融合结果如图3(c)所示。这里采用了双三次插值方法来放大多光谱图像的残差图像，通过待插值点附近的4×4个像素来确定待插值像素的灰度。采用的低通滤波器是11×11的均值低通滤波器。

为了与本发明方法的结果进行比较，这里，我们同时给出了标准的HIS变换融合方法、PCA融合方法、小波变换的融合方法的融合结果。这里的HIS变换融合方法采用圆柱模型把图像从RGB彩色空间变换为HIS模型，其融合结果如图2(d)所示。PCA融合方法则使用主成分分析方法，并用全色图像替换多光谱图像的第一主成分，其融合结果如图3(e)所示。小波变换的融合方法，首先使用HIS变换获得多光谱图像的亮度分量，然后通过小波变换方法融合全色图像和多光谱图像的亮度分量，采用了Daubechies的3阶小波系数，小波分解的层数选用3层，其结果如图3(f)所示。

实验结果如下：

首先，从视觉效果上进行分析和比较。对比图3(a)和图3(c)，本发明方法的融合图像纹理特征清晰，细节信息丰富。例如，在图3(a)所示图像的中部的栅格和左下角的机场跑道比较模糊，而在图3(c)的融合图像中可以清晰的观察到纵横交错的栅格和跑道的星簇图案。由图3(c)可看到，本发明方法较大地改善了图3(a)的空间分辨能力。

同时，从视觉效果上可以明显看到，其他的融合方法会引起多光谱图像光谱信息的严重失真。以HIS融合方法为例，在图3(d)中右上角的河流色彩接近于全色图像的黑色，而不是原始多光谱图像中的淡色。另外，图3(a)所示图像的中间区域呈现墨绿色，而在HIS融合结果中该区域的色彩改变较大。相反地，本发明方法的融合结果基本上没有上述的现象，图3(c)的融合图像的色彩信息与图3(a)的原始图像基本保持一致。

对于遥感图像融合效果的定量评价，应综合考虑空间细节的增强与光谱信息的保持两个方面。一般应综合考虑两类统计参数^[8]：一类反映空间细节信息，如方差、信息熵、清晰度；另一类反映光谱信息，如扭曲程度、偏差指数、相关系数。这里以多光谱图像第2波段(Green，0.52-0.60um)为例分析比较以上各种方法的融合效果，其结果如表1所示。

                       表1  基于统计参数的融合性能比较

图像	方差	信息熵	清晰度	扭曲程度	偏差指数	相关系数
							Green	35.461	4.833	4.535
PAN	49.008	5.207	12.669
				Residual	39.313	5.019	10.876	10.428	0.106	0.931
HIS	40.628	5.011	10.199								35.903	0.376	0.094
				PCA	35.585	4.929	9.105	33.374	0.342	0.145
WT	36.399	4.857	10.424								24.376	0.255	0.541

表1中的“Green”表示第2波段原始图像，“PAN”表示全色图像。“Residual”是本发明方法的融合结果，“HIS”表示HIS变换融合结果，“PCA”是PCA融合结果，“WT”代表小波变换融合结果。

表格中的方差、信息熵和清晰度用来衡量图像的空间细节信息，其定义如下。

(1)方差 ${\widehat{\mathrm{\σ}}}^2=\frac{1}{M\×N}\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}(f(i,j)-\stackrel{\‾}{f})---\left(4\right)$

其中，f(i，j)是图像上(i，j)点的灰度值， f是图像的均值，M和N是图像中行和列的像素的个数。

(2)信息熵 $H=-\underset{i=0}{\overset{v}{\mathrm{\Σ}}}{p}_i\×\log p_i---\left(5\right)$

对于[0，v]的灰度范围，p_i为灰度值等于i的像素数目与图像总的像素数目之比，v＝2ⁿ-1，n一般为6～10。

(3)清晰度 $\stackrel{\‾}{g}=\frac{1}{M\×N}\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}\sqrt{{\mathrm{\ΔI}}_x^2(i,j)+{\mathrm{\ΔI}}_y^2(i,j)},---\left(6\right)$

其中，ΔI_x与ΔI_y分别为x与y方向上的差分，M和N分别表示图像中行和列的像素的个数。

通常，较大的方差、信息熵和清晰度说明图像的灰度级分布比较均衡，图像携带的信息量较大，图像的空间细节比较丰富。从表1中可以看到，在本发明的方法中，融合图像的方差、信息熵和清晰度比原始图像要大，与其他融合方法的参数大致相当，这说明新方法成功地将全色图像的空间细节信息融入到多光谱图像中。

扭曲程度、偏差指数和相关系数反映了融合图像和原始图像之间光谱信息的保持情况，其定义如下。

(1)扭曲程度 $D=\frac{1}{M\×N}\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}|f(i,j)-g(i,j)|,---\left(7\right)$

(2)偏差指数 $D_{\mathrm{mdex}}=\frac{1}{M\×N}\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}\frac{|f(i,j)-g(i,j)}{f(i,j)},$

(3)相关系数 $C(f,g)=\frac{\mathrm{\Σ}[(f(i,j)-\stackrel{\‾}{f})\×(g(i,j)-\stackrel{\‾}{g})]}{\sqrt{\mathrm{\Σ}\left[{(f(i,j)-\stackrel{\‾}{f})}^2]\×\mathrm{\Σ}[{(g(i,j)-\stackrel{\‾}{g})}^2\right]}},---\left(9\right)$

其中，f(i，j)和g(i，j)分别为融合前和融合后图像(i，j)点的灰度值， f和 g分别为两幅图像的均值，M和N分别表示图像中行和列的像素的个数。

通常，较小的扭曲程度、偏差指数以及较大的相关系数，说明原始图像的光谱信息在融合的过程中得到了较好的保持。表1所示的结果表明，本发明方法将扭曲程度和偏差指数降低到了很小的程度，而相关系数是各种方法中最大的。这说明新方法保留了多光谱图像的绝大部分光谱信息，与现有的融合方法相比具有明显的优势。

综合以上，本发明方法在较大地增强融合图像的空间细节表现能力的同时，很好地保留了原多光谱图像的光谱信息，并且具有物理意义明确、实现结构简单、融合效果好的优点。实验结果表明，本方法的性能优于目前为止所广泛采用的基于HIS变换的融合方法、基于PCA的融合方法以及基于HIS变换和小波变换相结合的遥感图像融合方法等。新方法可为改善遥感图像的目视判读精度，提高信息清晰度和可靠性上提供新的有效的技术支持。

Claims (2)

1、一种基于残差的遥感图像融合方法，其特征在于使用残差抽取算法获取多光谱图像的残差图像和全色图像的残差图像；然后将多光谱图像的残差图像使用插值方法插值到和全色图像的残差图像同样的尺寸，通过主成分分析方法对插值后的多光谱残差图像和全色图像的残差图像进行融合，恢复出多光谱图像的高分辨率残差；最后将多光谱图像的插值图像和多光谱图像的高分辨率残差图像迭加，得到多光谱图像的高分辨率图像估计；其中：

所述残差抽取算法首先通过图像插值算法获得其插值图像，然后对该插值图像通过低通滤波和降采样过程获得原图像的低分辨率图像估计，最后用原图像与低分辨率的图像估计相减得到图像的残差图像；

所述插值方法选用双三次插值方法，通过待插值点附近的4×4个像素来确定待插值像素的灰度；

所述主成分分析方法是用全色图像的残差图像经过灰度调整后替换第一主成分，然后和其他的主成分一起经过主成分分析的逆变换得到多光谱图像的高分辨率残差图像估计，其步骤如下：

(1)计算参与融合的多光谱图像n个波段灰度图像的协方差矩阵S_x；

(2)算出S_x的全部特征值λ₁，λ₂，…λ_n和对应的特征向量u¹，u²，…，uⁿ，其中λ₁≥λ₂≥…≥λ_n；

(3)求出n个新特征向量Y＝(PC¹，PC²，…，PCⁿ)^T，它们满足Y＝U^TX，其中U＝(u¹，u²，…，uⁿ)^T，

而且S_y＝Λ＝diag{λ₁，λ₂，…，λ_n}；

此时的PC¹，PC²，…，PCⁿ分别被称为第1，2，…，n个主成分，它的逆变换为X＝UY。

2、根据权利要求1所述的基于残差的遥感图像融合方法，其特征在于把多光谱图像的插值图像I和多光谱图像的高分辨率残差图像e′_f逐像素相加得到多光谱图像的高分辨率图像估计MS′，

                    MS′＝I+e′_f。

Images