大多数遥感卫星，如SPOT、Landsat 7、IKONOS和QuickBird，能提供低分辨的多光谱图像(multispectral,MS)和高分辨率的全色图像(panchromatic,PAN)，在目标检测、地物分类等遥感图像的应用中，往往需要具有高空间和高光谱分辨率的多光谱图像，因此，多光谱图像和全色图像融合已成为当前重要的研究课题[1, 2]。
目前，MS与PAN图像的融合算法大致分为两类：基于成分替代的融合算法和基于多分辨率分析的融合算法。基于成分替代的融合算法首先对MS图像进行颜色空间变换，然后用PAN图像代替亮度分量，主要包括HSI变换法[3]、主成份分析[4](principal component analysis,PCA)法，该类算法在有效提高融合图像空间分辨率的同时却引入了严重的频谱失真。基于多分辨率分析的融合算法首先对MS图像和PAN图像进行多分辨率分解，然后对分解系数采用不同融合规则，再进行反重构得到融合图像。多分辨率分析的工具包主要包括：小波变换[4]、Contourlet变换[5]、非下采样的Contourlet变换(NSCT)[6]等，NSCT变换是一种有效的二维图像多尺度几何分析方法，与小波变换相比，它具有多个方向的基，不局限于小波变换的3个方向，同时基的支撑区间具有大小变化的长方形结构，满足各向异性尺度关系；与Contourlet变换相比，它取消了图像在分解和重构中的上采样及下采样操作，以数据的冗余度换取了平移不变性，避免了Contourlet变换的混频现象。NSCT变换为图像融合提供了有效的工具，对于融合系数的选择，通常是将2幅源图像低频系数通过加权平均法融合，而高频系数按照一定的策略从2幅源图像对应的高频系数选择最大值融合，这种融合规则忽略了图像特征和区域的相关性，因此，融合后的图像对比度下降并且边缘模糊。
针对上述问题，本文提出一种基于多特征的MS与PAN图像的融合方法，首先，将MS图像转换至HSI空间，以亮度分量(I)为标准对PAN图像进行直方图匹配，获得新的全色图分量；然后，对新的PAN图像分量和亮度分量分别进行NSCT分解；接着，对变换后得到的低频系数采用系数绝对值取大的融合规则，对高频系数提出了基于多特征的融合规则。最后通过逆NSCT变换和逆HSI变换得到融合图像。实验结果证明，本文所提出融合算法能够在保持光谱特性的前提下，尽可能提高融合图像的空间细节信息，具有较好的融合效果。
1 非下采样的Contourlet变换Contourlet变换[5]具多尺度，多分辨率和各向异性等特点。然而，由于下采样和上采样，Contourlet变换不具有移动不变性，此外，根据多样本率理论，采样结果在低频和高频部分会出现混频现象。以上提及的方法是Contourlet变换的最弱性能。
为了消除Contourlet变换混频现象和增强它的移动不变性，Cunha等人[7]提出了基于非下采样的金子塔分解和非下采样滤波器的NSCT变换。NSCT是Contourlet变换的平移不变版本，并通过不可分离双通道非下采样滤波器(NSFB)来获得平移不变性。与Contourlet变换不同的是，NSCT多分辨分解步骤是通过移动不变滤波器满足Bezout等式，而不是拉普拉斯金子塔滤波器来完成重构的。
由于在金子塔分解过程中没有采样，即使低通滤波器的带宽大于 
，低通滤波器也并不会带来频率混频现象。因此，NSCT变换比Contourlet变换具有更好的频率特性。图 1给出了NSCT分解的框架。
 |
| 图 1 NSCT分解框架 |
为了得到更好的空间分辨率和光谱信息，本文提出的MS与PAN图像融合算法的框架如图 2所示。
 |
| 图 2 本文融合算法的流程图 |
1)首先对MS和PAN图像进行配准，得到配准后的MS和PAN图，然后采用文献[8]提出的HSI空间转换方法，将低分辨率多光谱图像从RGB彩色空间转换到HSI空间，并得到色调(H)和饱和度(S)、亮度(I)；
2)对得到的亮度分量和全色图像采用文献[9]提供的方法进行直方图匹配，得到新的全色图分量Pnew；
3)对亮度分量和Pnew分别进行NSCT变换分解，得到低频系数和高频系数；
4)对低频系数和高频系数分别使用不同的融合规则，以此获得融合后的NSCT变换系数；
5)对融合后的系数进行逆NSCT变换，得到融合后的亮度分量；
6)对融合后的亮度分量、H和S进行逆HSI变换，得到最终的融合图像。
2.2 低频系数的融合规则经NSCT变换分解得到的低频部分包含了原始图像的概貌特征，为了更好地保留多光谱图像的光谱信息，并将全色图像的更多细节信息注入到多光谱图像中去，本文采用了低频系数像素绝对值选大的融合规则，融合后的低频系数为：
式中,LPnew(i,j)、LI(i,j)分别为图像经NSCT变换后在位置(i,j)处的低频系数,LF(i,j)为融合后的图像F对应的低频系数。 2.3 高频系数的融合规则
经NSCT变换分解获得的高频系数包含了图像边缘、区域边界等大量的细节信息,为了更好地在融合后的图像中反映这些信息,本文提出了一种多特征的高频子带系数融合方法并将该方法描述如下:
高频系数融合规则的选择将影响融合后图像的清晰度和光谱失真程度,其特征在标准差、平均梯度和能量对应有最大值,并且传统的方法以其不同的方式用来区分图像的清晰度。标准差是由均值间接求得的,反映了图像当前的像素值与整体像素均值的离散程度。标准差越大,整体像素灰度级分布就越分散,图像显示就越清晰。大小为M×N的图像F,其标准差可定义为:
式中,F(i,j)为像素在(i,j)位置像素的灰度值,表示图像的平均灰度值: 1) 基于标准方差的融合规则 式中,CPnewj(x,y)、CjI和CjF(x,y)分别是新全色图分量Pnew,MS图像的亮度分量I和融合图像F在(x,y)位置的高频系数。
平均梯度反映了融合图像的清晰度和微小细节反差的变化程度,其值越大,图像将更清晰,大小为M×N的图像F,其平均梯度定义为:
式中,F(i,j)是像素(i,j)的灰度值。
2) 基于平均梯度的融合规则
式中,CPnewj(x,y)、CjI和CjF(x,y)分别是新全色图分量Pnew,MS图像的亮度分量I和融合图像F在(x,y)位置的高频系数。
能量反映了图像的活跃水平和纹理的均匀程度,其值越大,图像包含的纹理信息分布越均匀,并且图像越清晰。图像F的能量定义为:
3) 基于能量的融合规则
式中,CPnewj(x,y)、CjI和CjF(x,y)分别是新全色图分量Pnew,MS图像的亮度分量I和融合图像F在(x,y)位置的高频系数。
在大多数现有MS与PAN图像的融合中,通常只考虑以上3个特征之一用来构造高频系数的融合规则。然而,以上提到的3个特征并不能表示图像的综合特征。甚至在一点上,不同的特征证明了不同表示能力。因此,在融合图像时,基于3种特征的单一融合规则会导致融合图像时丢失细节息,无法实现全局显示和质量评价综合指标的最优。为了处理传统方法用单一特征进行高频融合的问题,本文用以上的3种特征混合来进行高频子带系数的融合,其算法实现过程描述如下:
a) 在点(i,j)位置用方程(3)、(4)、(6)和(8)分别计算新全色图像分量Pnew和MS图像亮度分量I的高频系数的区域标准差DPnew(i,j)和DI(i,j),区域平均梯度GPnew(i,j)和GI(i,j),区域能量EPnew(i,j)和EI(i,j)。在这一步邻近区域的大小选择为(3×3);
b) 计算整幅图像的标准差DPnew和DI,平均梯度GPnew和GI,能量EPnew和EI;
c) 如果DPnew≠0,DI≠0,GPnew≠0,GI≠0,EPnew≠0,EI≠0,计算2幅图像的相对标准差,相对平均梯度和相对能量之间的比值如下:
式中,K1(i,j)是源图像在点(i,j)之间的相对标准差。如果,K1(i,j)>1新全色图分量Pnew在点(i,j)的相对标准差大于图像I,并且在图像Pnew内分散的灰度值就比在图像I内多。K2(i,j)和K3(i,j)做法与K1(i,j)相似,转到步骤e);
d) 如果全局标准差(DPnew或DI),平局梯度(GPnew或GI),能量(EPnew或EI)中任何一个等于0,它的特征表示在以上3个特征(方差、梯度、能量)之中是最差的,然后比较方程(13)~(15)的其他2种比率。对应率最大的值将被选择为融合;如果不同类型的2种特征都为0,融合规则将由全色图像的特征决定而不是以上提到的2种特征;如果不同类型的3种特征都为0,融合的系数将由2幅源图像上在这一点的最大系数决定。转到步骤f);
e) 通过比较以下的3种比率决定高频系数的融合规则:
对于任何一点,对应率的最大值规则将选择为在这一点的最终融合规则。例如,对于点(i,j),如果R1(i,j)拥有以上3种比率的最大值,基于方差的融合规则,将选择为这一点的融合规则。R2(i,j)和R3(i,j)的评价规则和R1(i,j)相似。
f) 对在源图像上的点重复步骤a)~e)。
3 实验结果和分析 3.1 实验实验选用Landsat-TM多光谱图像作为图像源，其中MS图像的地面分辨率是2.5 m，大小是256×256像素。SPOT全色图像的地面分辨率是10 m，大小为256×256像素。
为了验证推荐算法性能，将本文算法与传统的HSI变换、Contourlet变换、NSCT变换(融合规则均为低频加权平均、高频取局部方差选大规则)进行对比。Contourlet变换的分解参数为：4、4、8、8；NSCT变换与本文算法的分解参数均为1、3、4、4。
实验中采用的定量指标包括综合相对误差(relative global error in synthesis,ERGAS)[9]、相对平均光谱误差(relative average spectral error,RASE)[9]、均方根误差(root mean squared error,RMSE)[10]、平均梯度(average gradient,AG)、信息熵(information entropy,IE)和相关系数(correlation coefficient,CC)，其中ERGAS、RASE和RMSE几个指标的值越小表示融合后图像的光谱质量就越高；IE是衡量图像信息丰富程度的最重要指标。IE越大，图像就具有越多的信息；AG是图像中细节对比、纹理变化和清晰度表达的影响指标，AG越大，图像就具有越多的层次结构，融合图像就越清晰；CC反映了融合图像与多光谱图像的相关程度，相关系数值越大，则融合图像的光谱信息改变程度就越小。
3.2 实验结果分析实验采用的MS图像和PAN图像以及各种算法融合后的图像如图 3所示。从图 3可以看出，多光谱图(图 3a))包含有丰富的光谱信息，全色图(图 3b))包含有丰富的、显著的细节特征。从图 3c)可以看出，HSI变换法融合图像的结果光谱扭曲十分严重；从图 3d)对比可以看出Contourlet变换的融合方法虽能较好保留图像的多光谱信息，提高融合图像的空间分辨率，但是色彩失真明显。从图 3e)看出NSCT变换方法与图 3d)相比，色彩变淡。本文提出的融合方法(图 3f))与前3种方法相比能够更好地保留多光谱图像的光谱信息，从右下角部分放大图像可以看出，本文方法融合后的图像纹理更加清晰，细节轮廓更加丰富，色彩更明亮。
 |
| 图 3 源图像及融合实验结果 |
除了主观视觉效果外，实验还采用了客观评价指标对其进行了进一步评价。客观指标对比结果见表 1所示，从表 1中可以看出，本文方法所得的融合结果在指标ERGAS、RASE、RMSE、AG、IE和CC等都是最好的。其中，ERGAS表示融合后图像的光谱保持程度，ERGAS值越接近0表示光谱相关度信息保持的就越高。因此本文算法在光谱保持度方面是最优。而IE和AG表示空间信息的保留程度，其值越大，表示融合后图像的空间信息越接近于PAN图像。总体而言，本文算法在空间细节和光谱保持程度方面相对平衡，总体性能最好。
| 方法 | Best value | HSI | Contourlet | NSCT | Our |
| ERGAS | 0 | 11.373 3 | 16.115 7 | 8.660 4 | 7.029 3 |
| RASE | 0 | 39.354 5 | 46.868 2 | 33.702 0 | 29.337 9 |
| RMSE | 0 | 37.716 1 | 44.742 9 | 32.299 2 | 28.116 7 |
| AG | ∞ | 13.228 4 | 13.838 5 | 13.906 4 | 14.124 1 |
| IE | ∞ | 7.391 0 | 7.256 7 | 7.501 6 | 7.562 3 |
| CC | 1 | 0.713 9 | 0.667 8 | 0.786 7 | 0.857 1 |
由以上实验分析可知，本文提出的融合方法能够有效将全色图像的空间信息注入到多光谱图像中，并且能够很好保留多光谱图像的光谱信息，使融合后的图像光谱损失达到最小。
4 结 论本文提出了一种基于多特征的遥感图像融合算法，该方法对经NSCT变换分解得到的低频和高频系数的不同，分别制定了不同的融合规则，对于低频系数采用像素绝对值选大的融合规则；对于高频系数根据MS和PAN图融合选取不同区域特征方式，对融合图像的细节信息表示不一样，提出了多特征的融合规则。实验结果表明，文中算法能够很好保留多MS图像的光谱信息，有效将PAN图像的空间细节信息注入到多光谱图像中，取得更好的融合效果。
| [1] | Wald L. Some Terms of Reference in Data Fusion[J]. IEEE Trans on Geoscience and Remote Sensing, 1999, 37(3): 1190-1193 |
| Click to display the text | |
| [2] | Thomas C, Ranchin T, Wald L, Chanussot J. Synthesis of Multispectral Images to High Spatial Resolution: A Critical Review of Fusion Methods Based on Remote Sensing Physics[J]. IEEE Trans on Geoscience and Remote Sensing, 2008, 46(5): 1301-1312 |
| Click to display the text | |
| [3] | Shettigara V K . A Generalized Component Substitution Technique for Spatial Enhancement of Multispectral Images Using a Higher Resolution Data Set[J]. Photogrammetric Engineering and Remote Sensing, 1992, 58(5): 561-567 |
| Click to display the text | |
| [4] | Maria G A, Jose L S, Raquel G C, et al. Fusion of Multispectral and Panchromatic Images Using Improved HIS and PCA Mergers Based on Wavelet Decomposition[J]. IEEE Trans on Geoscience and Remote Sensing, 2004, 42(6):1291-1299 |
| Click to display the text | |
| [5] | 陈志刚,尹福昌. 基于Contourlet变换的遥感图像增强算法[J]. 光学精密工程,2008, 16(10): 2031-2037 Chen Z G, Yin F C. Enhancement of Remote Sensing Imgae Based on Contourlet Transform[J]. Opt Precision Eng, 2008, 16(10): 2031-2037 (in Chinese) |
| Cited By in Cnki (44) | Click to display the text | |
| [6] | 冯鑫,王小明,党建武,等. 基于改进非下采样轮廓波的融合算法[J]. 农业机械学报,2012,43(12):192-196 Feng Xin, Wang Xiaoming, Dang Jianwu, et al. Image Fusion Algorithm Based on Improved Non-Subsampled Contourlet[J]. Transactions of the Chinese Society for Agricultural Machinery, 2012, 43(12): 192-196 (in Chinese) |
| Cited By in Cnki (8) | |
| [7] | Cunha A L, Zhou J P, Donoho M N. The Non-Subsampled Contourlet Transforms: Theory, Design, and Applications[J]. IEEE Trans on Image Processing, 2006, 15(10): 3089-3101 |
| Click to display the text | |
| [8] | Tu T, Su S, Shyn H, Huang P. A New Look at IHS Like Image Fusion Methods[J]. Information Fusion, 2001(2): 177-186 |
| Click to display the text | |
| [9] | Tania Stathaki. Image Fusion: Algorithms and Applications[M]. 1st ed, Elsevier, 2008 |
| Click to display the text | |
| [10] | Du Q, Younan N, King R, Shah V. On the Performance Evaluation of Pan-Sharpening Techniques[J]. IEEE Geoscience and Remote Sensing Letters, 2007, 4(4): 518-522 |
| Click to display the text |