CN103236047A - 一种基于替换分量拟合的全色与多光谱图像融合方法 - Google Patents

Abstract

本发明公开了一种基于替换分量拟合的全色与多光谱图像融合方法，该方法首先对全色图像进行亮度拉伸，并将多光谱图像插值成与全色图像相同分辨率的插值图像；对多光谱插值图像进行矩阵变换，得到待替换分量和保留分量，并将全色图像与待替换分量做减法，得到差值图像；然后，对差值图像进行高斯滤波，得到差值趋势图像，并将全色图像减去差值趋势图像，得到拟合的全色图像；最后，将拟合的全色图像替换待替换分量，并进行矩阵逆变换得到融合图像。与现有技术相比较，本融合方法避免现有分量替换分量的光谱失真现象，适应性好。

Description

一种基于替换分量拟合的全色与多光谱图像融合方法

技术领域

本发明涉及一种遥感图像融合方法，尤其涉及一种基于替换分量拟合的全色与多光谱图像融合方法，属于数字图像处理技术领域。

背景技术

随着遥感成像技术的不断发展，世界各科技强国发射了越来越多的多传感器卫星，例如IKONOS、QuickBird、ALOS、WordView-2等卫星同时搭载了全色与多光谱等传感器。为了满足我国遥感应用的需求，在2012年发射了资源三号光学卫星，首次实现了中国民用卫星的多角度、多光谱综合立体成像，获取了大量的全色和多光谱遥感图像。随着我国高分辨率对地观测系统国家重大科技专项的推进，届时将拥有数量更多、分辨率更高的多光谱和全色成像卫星。

由于传感器成像原理不同，单一全色图像或单一多光谱图像难以全面反映图像中目标的特性，具有一定的应用局限制，因此需要对全色与多光谱图像进行融合操作。全色与多光谱图像融合是将同一地区的两种图像数据进行像素级合成，得到更加适合图像解译人员分析或机器处理的复合图像。一般而言，全色与多光谱图像融合必须满足以下要求：融合图像的光谱信息必须与多光谱图像的光谱信息保持一致，即光谱保真；融合图像的空间细节信息必须与全色图像的空间细节信息保持一致，即空间细节保真。

迄今为止，研究人员已经提出了大量的图像融合方法，其中研究最为广泛且应用效果最好的方法是分量替换融合方法。它首先利用矩阵变换生成待替换分量和保留分量；然后对全色图像和待替换分量进行匹配，使全色图像的整体结构信息与待替换分量的整体结构信息保持一致；最后，利用匹配的全色图像替换待替换分量，再进行矩阵逆变换得到融合图像。根据矩阵变换类型的不同，分量替换融合方法可以分为基于IHS变换、基于PCA变换和基于Gram-Schmidt变换（GS变换）的分量替换融合方法。这些方法的优点在于融合图像的空间细节保真度高，但是同时存在着光谱失真的严重不足。因此，解决分量替换融合方法的光谱失真问题已成为当前全色与多光谱图像融合研究的热点与难点。目前，研究人员侧重于改进矩阵变换方法，以生成更好的待替换分量，即待替换分量的整体结构信息与全色图像的整体结构信息处处相等。尽管如此，在部分图像区域内，二者的整体结构信息仍然存在较大的差异，导致融合图像的光谱出现严重失真。

在此背景下，研究一种可以有效避免融合图像光谱和空间细节失真的高保真融合方法，有效地避免现有分量替换融合方法的光谱失真现象显得尤为重要。

发明内容

本发明所要解决的技术问题在于提供一种基于替换分量拟合的全色与多光谱图像融合方法。该方法通过矩阵变换、图像做差和图像平滑滤波操作实现多光谱与全色图像的高保真融合，适应性好。

为实现上述的发明目的，本发明采用下述的技术方案：

一种基于替换分量拟合的全色与多光谱图像融合方法，包括如下步骤：

步骤（1）对全色图像进行亮度线性拉伸，使全色图像的方差σ^P与多光谱图像的方差σ^MS相等；

步骤（2）将多光谱图像插值成与全色图像同样分辨率的插值图像MS_r；

步骤（3）对多光谱插值图像MS_r进行矩阵变换，得到待替换分量与保留分量；

步骤（4）利用所述线性拉伸后的全色图像减去所述的待替换分量，得到差值图像I_d，并对差值图像I_d进行高斯滤波，得到差值趋势图像I_s；

步骤（5）将所述线性拉伸后的全色图像与差值趋势图像I_s相减，得到拟合的全色图像；

步骤（6）利用所述拟合的全色图像替换所述的待替换分量，并与保留分量一起进行矩阵逆变换，得到全色与多光谱融合图像。

其中，在所述步骤（1）中，将全色图像的像素值乘以

进行线性亮度拉伸。

所述步骤（1）中，随机选取一定数量的像素点来估计多光谱图像的方差以及全色图像的方差。

所述步骤（2）中，利用双线性插值法进行图像插值放大。

所述步骤（3）中，所述的矩阵变换为IHS变换、PCI变换、GS变换中的一种。

所述步骤（4）中，利用一维高斯滤波器G先对图像按行进行平滑滤波，再利用所述一维高斯滤波器G的转置G^T对图像按列进行平滑滤波。

所述步骤（4）中，一维高斯滤波器G的长度为13，G的尺度参数为3.33。

本发明所提供的基于替换分量拟合的全色与多光谱图像融合方法具有如下的优点：

1.本发明较为有效地避免了融合图像的光谱和细节失真，融合图像的视觉效果好；

2.本发明方法简单易行，计算复杂度低，可快速完成宽幅遥感图像的融合；

3.本发明方法的参数少，普适性好，可适用于不同卫星拍摄的全色和多光谱图像。

附图说明

下面结合附图和具体实施方式对本发明作进一步的详细说明。

图1为本发明所述融合方法的流程示意图；

图2显示了本发明所述融合方法与现有典型融合方法在QuickBird卫星图像上的实验结果对比；其中，(a)为全色图像，(b)为多光谱图像，(c)为Massip的方法，(d)为ENVI-GS变换法，(e)为PCI锐化法，(f)为本发明的方法；

图3显示了本发明所述融合方法与现有典型融合方法在IKONOS卫星图像上的实验结果对比；其中，(a)为全色图像，(b)为多光谱图像，(c)为Massip的方法，(d)为ENVI-GS变换法，(e)为PCI锐化法，(f)为本发明的方法。

具体实施方式

本发明主要利用数据拟合方法对矩阵变换生成的待替换分量与全色图像进行拟合，使二者的整体结构信息保持一致，从而实现全色与多光谱图像的高保真融合。图1展示了基于替换分量拟合的全色与多光谱图像融合方法的流程图。本发明的基本思路如下：首先对全色图像进行亮度线性拉伸，并将多光谱图像插值至全色图像相同的分辨率，并对其进行矩阵变换，得到待替换分量与保留分量；接着，将全色图像分别减去待替换分量，得到差值图像，并对差值图像进行高斯平滑，得到差值趋势图像；然后，将全色图像减去差值趋势图像，得到拟合的全色图像；最后，对拟合的全色图像和保留分量进行矩阵逆变换，得到融合图像。

下面具体介绍本发明所提供的基于替换分量拟合的全色与多光谱图像融合方法的具体步骤：

（1）对全色图像进行亮度线性拉伸，使全色图像的方差σ^P与多光谱图像的方差σ^MS相等；

首先分别计算多光谱图像和全色图像的方差。为了提高方差的计算速度，随机选取一定数量的像素点来估计多光谱图像和全色图像的方差，以减少方差计算的时间消耗。然后将全色图像的像素值乘以

进行线性亮度拉伸。

（2）将多光谱图像插值成与全色图像同样分辨率的插值图像MS_r；

多光谱图像的分辨率低于全色图像，图像融合利用插值算法将多光谱图像插值至全色图像相同的分辨率。本发明利用双线性插值法对多光谱图像进行插值。假设像素点I(i+u,j+v)为像素点{I(i,j),I(i+1,j),I(i,j+1),I(i+1,j+1)}之间的一个插值点，其中0<u<1，0<v<1，i，j为像素点的行坐标和列坐标，则插值的计算方法具体如下：

I(i+u,j+v)=(1-u)(1-v)I(i,j)+uvI(i+1,j+1)

+v(1-u)I(i,j+1)+u(1-v)I(i+1,j)

（3）对多光谱插值图像MS_r进行矩阵变换，得到待替换分量与保留分量；

在该步骤中，可以选择IHS变换、PCA变换和GS变换等对多光谱插值图像进行矩阵变换，且融合图像的保真效果基本相同。考虑到算法的计算复杂度，优选IHS变换进行矩阵变换。下面以IHS变换为例介绍矩阵变换的计算过程。在图像MS_r中，像素点(i,j)的光谱值记为

经IHS变换得到

具体如下：

$\left[\begin{array}{c}{I}_{i,j}\\ H_{i,j}\\ S_{i,j}\end{array}\right]=\left[\begin{array}{ccc}\frac{1}{3}& \frac{1}{3}& \frac{1}{3}\\ \frac{-\sqrt{2}}{6}& \frac{-\sqrt{2}}{6}& \frac{2\sqrt{2}}{6}\\ \frac{1}{\sqrt{2}}& \frac{-1}{\sqrt{2}}& 0\end{array}\right]\left[\begin{array}{c}{\mathrm{MS}}_{i,j}^1\\ {\mathrm{MS}}_{i,j}^2\\ {\mathrm{MS}}_{i,j}^3\end{array}\right]$

（4）利用所述线性拉伸后的全色图像减去所述的待替换分量，得到差值图像I_d，并对差值图像I_d进行高斯滤波，得到差值趋势图像I_s；

将全色图像与待替换分量作矩阵减法，得到差值图像I_d。差值图像包含两部分信息：一部分是全色图像的空间细节信息，另一部分是全色与待替换分量的整体结构信息差异。

对差值图像进行高斯平滑是为了滤除其中的空间细节信息，得到全色与待替换分量的整体结构信息差异。为了快速地完成高斯滤波操作，利用一维高斯滤波器G先对图像按行进行平滑滤波，再利用所述一维高斯滤波器G的转置G^T对图像按列进行平滑滤波。其中，一维高斯滤波器G的长度为13，G的尺度参数为3.33。

（5）将所述线性拉伸后的全色图像与差值趋势图像I_s相减，得到拟合的全色图像；

在该步骤中，全色图像与差值趋势图像I_s作矩阵减法即可。此时，拟合的全色图像与待替换分量之间的整体结构信息基本相同，可以较好地避免融合图像出现光谱失真。

（6）利用所述拟合的全色图像替换所述的待替换分量，并与保留分量一起进行矩阵逆变换，得到全色与多光谱融合图像。

在该步骤中，首先计算变换矩阵的逆矩阵，然后将该逆矩阵与“拟合的全色图像和保留分量构成的向量”进行矩阵相乘，得到全色与多光谱融合图像。下面以IHS逆变换为例介绍矩阵逆变换的计算过程。在像素点(i,j)处，拟合的全色图像记为P_i,j，保留分量记为[H_i，j,S_i，j],“拟合的全色图像和保留分量构成向量”记为[P_i,j，H_i,j,S_i,j],经IHS逆变换得到

（F即为融合图像），具体如下：

$\left[\begin{array}{c}{F}_{i,j}^1\\ F_{i,j}^2\\ F_{i,j}^3\end{array}\right]={\left[\begin{array}{ccc}\frac{1}{3}& \frac{1}{3}& \frac{1}{3}\\ \frac{-\sqrt{2}}{6}& \frac{-\sqrt{2}}{6}& \frac{2\sqrt{2}}{6}\\ \frac{1}{\sqrt{2}}& \frac{-1}{\sqrt{2}}& 0\end{array}\right]}^{-1}\left[\begin{array}{c}{\mathrm{MS}}_{i,j}^1\\ {\mathrm{MS}}_{i,j}^2\\ {\mathrm{MS}}_{i,j}^3\end{array}\right]$

本发明方法与主流遥感图像处理软件中最好的融合方法，即ENVI软件的Gram-Schmidt变换融合法（简称ENVI-GS变换法）、PCI软件的全色锐化融合法（简称PCI锐化法），以及最新的融合方法，即“Massip P.，Blanc P.，Wald L.A method to better account for modulationtransfer functions in ARSIS-based pansharpening methods[J].IEEE Trans.Geosci.Remote Sens.，2012，V50(3):800-808”（简称Massip的方法）进行了实验对比。实验数据为QuickBird卫星和IKONOS卫星等拍摄的全色与多光谱图像，共计19景。其中，全色图像平均每景尺寸约为12000×12000像素，多光谱图像每景尺寸约为3000×3000像素。

图2和图3分别展示QuickBird卫星和IKONOS卫星的全色和多光谱图像以及融合图像。由于实验图像的尺寸较大，为了清晰地展示图中的地物，图2和图3仅给出了实验图像的局部区域。在主观视觉效果上，本发明的方法对所有的实验数据均有较好的光谱和空间细节保真效果见图2(f)和图3(f)。但是，对比方法的均存在失真现象：对于Massip的方法，在部分区域融合图像的空间细节变模糊，见图2(c)和图3(c)；对于ENVI-GS变换法，在某些区域内融合图像的光谱严重失真，见图2(d)和图3(d)；对于PCI锐化法，在某些区域内融合图像的光谱出现失真，见图2(e)。

下面利用在光谱扭曲度、光谱角和全局相对误差来评价融合图像的光谱保真效果，同时利用细节扭曲度来评价融合图像的空间细节保真效果。这些指标的取值越小，融合图像的保真效果越小，其定义如下：

（1）光谱扭曲度（Spectral Distortion，SD）

光谱扭曲度反映了融合图像[F_i,j,k]_m×n×p与多光谱图像[M_i,j,k]_m×n×p的光谱失真程度，其中m和n为图像的高度与宽度，p为图像的波段数量，G为高斯滤波器。它的取值越小，说明融合图像的光谱失真越小，计算公式如下：

$\mathrm{SD}=\frac{1}{m\&times;n\&times;p}\underset{k=1}{\overset{p}{\mathrm{\&Sigma;}}}\underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{n}{\mathrm{\&Sigma;}}}|{(F*G)}_{i,j,k}-M_{i,j,k}|$

（2）细节扭曲度（Detail Distortion，DD）

细节扭曲度是指融合图像[F_i,j,k]_m×n×p与全色图像[P_i,j]_m×n空间细节信息差值的平均值，其中m和n为图像的高度与宽度，p为图像的波段数量，G为高斯滤波器。它的取值越小，融合图像的细节保真程度越好，计算公式如下：

$\mathrm{DD}=\frac{1}{m\&times;n\&times;p}\underset{k=1}{\overset{p}{\mathrm{\&Sigma;}}}\underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{n}{\mathrm{\&Sigma;}}}|{(F-F*G)}_{i,j,k}-{(P-P*G)}_{i,j,}|$

表1列出了本发明方法与对比方法的客观评价结果。观察表1的客观指标可知，在本发明融合图像的光谱扭曲度和细节扭曲度均比对比方法小。这说明本发明方法的光谱和空间细节保真效果优于对比方法。在配置为3.2GHz的4核CPU，2GB的内存，Windows XP操作系统的计算机上，ENVI-GS变换法每景实验数据融合平均耗时170.73秒，PCI锐化法每景实验数据融合平均耗时163.26秒，本发明方法每景实验数据融合平均耗时24.17秒。实验结果表明本文融合方法的光谱保真和细节保真效果好、计算速度快，优于对比方法。

表1  图像融合质量的客观评价表

	光谱扭曲度	细节扭曲度
			Massip的方法	33.29	6.37
ENVI-GS变换法	28.24	5.11
			PCI锐化法	9.62	6.13
本发明的方法	5.36	4.85

本发明所提供的基于替换分量拟合的全色与多光谱图像融合方法主要是为提高全色与多光谱图像融合的保真效果专门提出的。但显然，本融合方法也适用于全色与高光谱图像，以及普通成像设备如数码相机所拍摄图像的融合操作，所取得的有益效果也是相似的。

以上对本发明所提供的基于替换分量拟合的全色与多光谱图像融合方法进行了详细的说明，但显然本发明的具体实现形式并不局限于此。对于本技术领域的一般技术人员来说，在不背离本发明的权利要求范围的情况下对它进行的各种显而易见的改变都在本发明的保护范围之内。

Claims (7)

1.一种基于替换分量拟合的全色与多光谱图像融合方法，其特征在于，包括如下步骤：

步骤（1）、对全色图像进行亮度线性拉伸，使全色图像的方差σ^P与多光谱图像的方差σ^MS相等；

步骤（2）、将多光谱图像插值成与全色图像同样分辨率的插值图像MS_r；

步骤（3）、对多光谱插值图像MS_r进行矩阵变换，得到待替换分量与保留分量；

步骤（4）、利用所述线性拉伸后的全色图像减去所述的待替换分量，得到差值图像I_d，并对差值图像I_d进行高斯滤波，得到差值趋势图像I_s；

步骤（5）、将所述线性拉伸后的全色图像与差值趋势图像I_s相减，得到拟合的全色图像；

步骤（6）、利用所述拟合的全色图像替换所述的待替换分量，并与保留分量一起进行矩阵逆变换，得到全色与多光谱融合图像。

2.如权利要求1所述的基于替换分量拟合的全色与多光谱图像融合方法，其特征在于：在所述步骤（1）中，将全色图像的像素值乘以

进行线性亮度拉伸。

3.如权利要求2所述的基于替换分量拟合的全色与多光谱图像融合方法，其特征在于：所述步骤（1）中，随机选取一定数量的像素点来估计多光谱图像的方差以及全色图像的方差。

4.如权利要求1所述的基于替换分量拟合的全色与多光谱图像融合方法，其特征在于：在所述步骤（2）中，利用双线性插值法对多光谱图像进行插值。

5.如权利要求1所述的基于替换分量拟合的全色与多光谱图像融合方法，其特征在于：所述步骤（3）中，所述的矩阵变换为IHS变换、PCA变换、GS变换中的一种。

6.如权利要求1所述的基于替换分量拟合的全色与多光谱图像融合方法，其特征在于：所述步骤（4）中，利用一维高斯滤波器G先对图像按行进行平滑滤波，再利用所述一维高斯滤波器G的转置G^T对图像按列进行平滑滤波。

7.如权利要求5所述的基于替换分量拟合的全色与多光谱图像融合方法，其特征在于：所述步骤（4）中，一维高斯滤波器G的长度为13，G的尺度参数为3.33。

Images