遥感图像融合是采用合适的算法将全色图像和多光谱图像进行融合，目的是获得一幅光谱信息丰富并且空间分辨率高的遥感图像。在多分辨率分析的方法中，研究人员最先提出了小波变换，虽然其可以获得3个方向的高频信息，但是由于其不具备各向异性，无法得到图像多方向的特征，因此近年来，更多的研究人员将Contourlet变换和Shearlet变换应用到图像融合领域。TORO-GARAY et al[1]将Contourlet变换应用到了遥感图像进行边缘增强，取得了不错的效果。但Contourlet变换会对图像进行下采样操作，不具备平移不变性，并且由于其自身原因，容易引入频谱混叠等现象，影响融合结果的细节和光谱信息。针对此问题，2016年，DACUNHA et al[2]提出了NSCT变换，将图像在非下采样的情况下完成了多尺度和多方向的分解，其局部性在时间和频域上表现优秀。在使用NSCT变换对图像进行分解重构时，具有多尺度和多方向的优势，可以更加准确地表达图像的细节信息，在遥感图像融合领域有很好的应用前景[3]。
随着机器学习的飞速发展，各行各业的研究人员都将其引入到了各自的领域，遥感图像处理领域也不例外。2010年，ZEILER et al[4]提出了卷积稀疏表示，采用卷积的形式来表示图像或信息，考虑到了像素间的语义关系，并且从本质上来讲该模型具有平移不变性，可以更好地保留图像的细节信息。ZHANG et al[5]利用卷积稀疏表示建立融合模型，来获取融合结果，使得空间分辨率和光谱信息都得到了提升，但是该方法没有将细节信息和光谱信息分开处理，在融合过程中不可避免地产生了光谱损失。WU et al[6]在融合框架中用到了卷积稀疏表示，提高了融合图像的性能，但是其分解高低频子带的方式过于简单，无法获得图像多方向的特征信息，在融合过程中会造成细节信息的缺失。本文结合NSCT变换和卷积稀疏表示的优势，提出了基于卷积稀疏表示和NSCT变换的遥感图像融合，通过实验，并结合主观和客观两方面指标评价，证明了本文方法的有效性。
1 相关理论介绍
1.1 卷积稀疏表示
卷积稀疏表示[4]是传统稀疏表示的一种替代形式，采用字典滤波器和特征映射的卷积形式来替代字典和稀疏系数的乘积形式，它通过应用卷积形式来实现整个信号或图像的稀疏表示，充分考虑到了图像内的语义关系，有效解决了稀疏表示中因为“分块”导致的细节缺失问题。具体来讲，卷积稀疏表示是将整个信号或图像建模为特征映射和字典滤波器的一组卷积的和，其中特征映射的大小与源图像或信号相同。具体公式如下：
(1)
式中：s是待稀疏表示的信号或图像；{dm}是m维的字典滤波器；{xm}是一组特征映射，它与待表示的信号或图像具有相同的尺寸；λ表示模型的准确性和稀疏性的平衡；*定义为卷积运算。
1.2 NSCT变换
NSCT变换[2]的提出是在Coutourlet变换的基础上进行了改变，去掉了Coutourlet变换中采样步骤，使图像处理过程不再出现伪吉布斯现象，减小了融合结果中的光谱扭曲现象。NSCT变换的最大优势在于变换后的图像具有平移不变性，可以保留原图像更多的细节信息。NSCT变换之所以有此优势，根本原因在于其将多尺度分解和方向滤波器结合使用。
NSCT变换[7]主要由非下采样金字塔滤波器和非下采样方向滤波器组成。变换过程中先采用非下采样金字塔滤波器进行多尺度分解，对于分解结果中的高频子带使用非下采样方向滤波器，从而得到图像各个方向中更详细的细节特征，其中所有的子带图像都具有和源图像相同的尺寸大小，低频部分继续进行上述分解。具体框架图如图1所示。
图1 NSCT框架[2]
Fig.1 NSCT framework[2]
2 基于卷积稀疏表示和NSCT的遥感图像融合
本文算法结合两种算法的优势，提出了卷积稀疏表示和NSCT相结合的遥感图像融合算法。将融合过程分为两部分，分别为图像的超分辨率分析和图像融合，融合流程如图2所示。先使用卷积稀疏表示建立模型完成图像的超分辨率分析。之后再利用NSCT完成融合步骤，先将超分辨后的图像和原始的Pan图像进行NSCT分解，由于高频子带和低频子带所含有的信息不同，为了更好地在融合结果中保留空间信息和光谱信息，针对高频和低频子带的不同特点分别采用了不同的融合策略来获得新的子带信息，最后进行NSCT逆变换得到融合结果。
图2 NSCT和卷积稀疏表示结合的融合算法流程图
Fig.2 Flow chart of fusion algorithm based on NSCT and
convolutional sparse representation
2.1 基于卷积稀疏表示的超分辨率重构
通过卷积稀疏表示来完成图像的超分辨率(super-resolution，SR)，增强图像的细节信息。首先，训练得到图像高分辨的字典滤波器和低分辨率的字典滤波器，使用低分辨率的字典滤波器表示超分辨率之前的图像得到特征映射，最后用特征映射和高分辨率的卷积完成图像的超分辨率，得到细节增强的图像。由于获取的原始Pan图像和MS图像大小并不匹配，因此在实验之前对数据进行了统一的配准和插值处理，获得大小相同的Pan图像和MS图像。超分辨率阶段用到的字典滤波器在实验中大小统一为8×8×32[8].
具体步骤如下：
1) 对MS图像进行IHS变换，将MS图像转换到IHS空间，得到以下3个分量，分别为亮度(Intensity，I)、色调(Hue，H)和饱和度(Saturation，S)三个独立分量，后续的步骤主要在I分量上进行；
2) 采用形态学滤波器对Pan图像进行滤波，得到低分辨率的Pan图像，即LPan.为了获得更加准确的低频信息，采用了由RESTAINO et al[9]提出的由半梯度形态学算子来组成形态学滤波器，对应的滤波器如式(2)所示：
0.5(δB[I]-I) .
(2)
式中：εB[I]为腐蚀算子；δB[I]为膨胀算子；I表示图像。
3) 对Pan图像和LPan图像进行训练，得到高低分辨率的字典滤波器DH和DL.根据文献[8]将字典学习定义为式(3)的优化问题。
s.t.‖dm‖2=1 .
(3)
式中：dm为对应的字典滤波器；xm为特征映射；s为图像。
4) 用低分辨率的字典滤波器DL对I分量进行卷积稀疏表示，得到对应的特征映射X，式(4)所示：
(4)
5) 用4)得到的特征映射X和3)中的高频字典滤波器DH结合，对I分量进行重构，得到重构后的I分量，即ISR.式(5)所示：
ISR=DH*X.
(5)
对应的流程图如图3所示：
图3 图像超分辨率算法框架
Fig.3 Super resolution algorithm framework of image
2.2 低频融合规则
将超分辨后的I分量和Pan图像进行NSCT变换，得到各自的高、低频子带图像，分别为PL、PH和IH、IL.由于低频信息中主要是图像的概要信息，包含的细节信息少，因此为了将更多的细节信息保留到融合结果中，这里采用了绝对值取大的融合规则。其次，由于低频子带具有较低的稀疏性，为了更好地保留细节信息，在融合规则中引入了稀疏表示，采用稀疏系数绝对值取大的规则。先利用图像库中的图像训练得到字典，之后用字典表示低频子带，得到对应的稀疏系数，即xI和xP，同一位置对应的更大的稀疏系数即为融合结果的稀疏系数。具体的公式如下：
FL(x,y)=
(6)
2.3 高频融合规则
源图像大量的细节信息和纹理信息都存在于高频子带图像中，同时，也有少量的光谱信息存在于高频子带中，为了确保在融合步骤中尽可能少的损失光谱信息和空间信息，故将局部方差和加权平均结合的融合规则应用到高频子带融合。
首先，计算以像素{x,y}为中心，窗口大小为n*n的局部方差。相应的计算公式如下：
(7)
(8)
然后根据式(7)和式(8)来计算各自的权重系数，公式如下所示：
(9)
(10)
使用式(9)和式(10)计算得到的权重，并与式(11)结合即可得到融合后的高频子带。公式如下所示：
FH=WIH(x,y)×IH(x,y)+
WPH(x,y)×PH(x,y) .
(11)
2.4 融合算法描述
对重构后的I分量和Pan图像进行NSCT变换，然后针对变换得到的子带信息进行融合得到新的子带信息，最后进行NSCT逆变换得到融合结果，其中NSCT变换具有平移不变性，因此分解前后图像的尺寸不会发生变换，融合图像的尺寸与Pan图像尺寸相同。流程图如图2所示，具体的步骤如下所示。
1) 将I分量和Pan图像进行NSCT变换，得到I分量和Pan图像的高频子带图像和低频子带图像，即PL、PH和IH、IL.
2) 根据高频和低频子带各自的特点采取不同的融合规则，得到新的高频和低频子带图像，即FL和FH.
3) 将2)中的子带图像进行NSCT逆变换得到新的I分量。
4) 将新的I分量和H、S分量进行IHS逆变换，得到融合图像。
3 实验
3.1 数据来源和预处理
本文所采用的实验数据来源于QuickBird 卫星和法国国家太空研究中心(centre national d’etudes spatiales，CNES).其中QuickBird卫星由美国DigitalGlobe公司负责发射，获取的MS图像分辨率和Pan图像分辨率分别为2.44 m和0.61 m，此处为了实验方便将MS图像插值放大至与Pan图像相同大小，图像尺寸为512像素×512像素。CNES获取的图像全色图像分辨率为0.6 m，此处为了实验方便将MS图像插值放大至与Pan图像相同大小，图像尺寸为1 024像素×1 024像素，处理好的图像如图4和图5所示。
图4 QuickBird卫星图像
Fig.4 QuickBird satellite image
图5 CNES卫星图像
Fig.5 CNES satellite image
3.2 实验结果和分析
为证明本文算法的优势，将其与5种算法进行对比，分别为IHS变换[10]的融合算法(IHS)，基于形态学滤波器的融合算法[9](HG)，基于稀疏表示的融合算法[11](SR)，基于乘法注入细节融合算法(HPM)[12]，基于NSCT的融合算法[13](NSCT).为了对融合结果进行更加客观的评价，本文选取了以下5个评价指标，分别为相对全局融合误差(erreur relative global adimensionnelle de synthese，ERGAS)[12]、光谱角(spectral angle mapper，SAM)[13]、通用图像质量评价指标(universal image quality index，UIQI)[14]、相关系数(correlation coefficient，CC)、标准差(standard deviation，STD).
图6和图7的源图像均来自于QuickBird 卫星，其中图6为印度某公路段(Road)，图7为印度某国家公园湖泊(Lake)附近景区。为了便于观测分析，选取了图中有代表性的区域进行标记并放大，放大区域位于图像的右下角或左上角。图6(c)和图7(c)为使用形态学滤波器获得的融合结果，观测放大区域可以看出，其细节信息丰富，但是出现了光谱畸变，与源多光谱图像相比，这两幅图像草坪的颜色都要更深一点。图6(d)和图7(d)为IHS变换得到的融合结果，由于其是用全色图像直接替代I分量得到融合结果，因此虽然其细节信息也较为丰富，但是出现了严重的光谱损失，原本绿色的草坪现在颜色已经开始发白。
图6 Road源图像及不同方法获得的融合结果
Fig.6 Fusion results of original Road image and different methods
图7 Lake源图像及不同方法获得的融合结果
Fig.7 Fusion results of original Lake image and different methods
图6(e)和图7(e)为稀疏表示的融合方法，由于卷积稀疏表示自身存在的“滑窗”技术，破坏了图像块之间的相关性，会导致融合结果的光谱失真和细节损失，观察放大区域会发现道路的边缘不是很清晰，而且颜色也没有很明显的分界线。图6(f)和图7(f)为HPM的融合算法，相比于加法注入细节，该算法可以使得融合结果中注入更多的空间细节信息，图中可以看出其细节信息还是很清晰的，但是也存在一些光谱畸变。图6(g)和图7(g)为NSCT变换的融合结果，NSCT变换具有平移不变性，在保持图像信息方面表现优异，因此可以看出图像中细节信息丰富，但是存在轻微的光谱扭曲现象。图6(h)和图7(h)为本文提出的算法，先对多光谱图像进行超分辨率的细节增强，然后再进行后续的融合步骤。从放大区域可以看出，图像具有更加丰富的细节信息，且光谱畸变也较小。对上述两组图像性能的比较结果如表1和表2所示。
表1 Road图像融合结果的性能比较
Table 1 Performance comparison of Road
image fusion results
MethodCCERGASSAMSTDUIQIHG0.802 05.887 91.673 20.089 80.763 5IHS0.346 99.744 61.164 20.167 10.346 9SR0.561 67.705 41.054 70.127 90.543 4HPM0.834 913.416 50.807 60.087 70.607 7NSCT0.815 57.047 61.156 80.084 80.784 7proposed0.891 94.019 40.557 20.066 30.888 2
表2 Lake图像融合结果的性能比较
Table 2 Performance comparison of Lake
image fusion results
MethodCCERGASSAMSTDUIQIHG0.748 414.510 20.166 70.116 70.552 5IHS0.486 714.230 72.885 40.161 80.389 6SR0.630 78.269 41.486 80.141 20.594 7HPM0.845 214.456 01.298 40.105 70.576 7NSCT0.834 37.883 62.093 80.098 80.800 2proposed0.912 84.162 80.743 00.074 00.912 5
从上述表中可以看出，HG算法获得的融合结果虽然CC和UIQI两个指标值较好，甚至在表2中SAM的值最佳，但是ERGAS的值均偏高，说明该方法存在较为严重的光谱失真。IHS算法获得的融合结果同样存在ERGAS和SAM两个指标值偏大的现象，说明该算法同样存在较为严重的光谱失真。SR算法获得的融合结果的指标与IHS相比，虽然各项指标的值都有所好转，但是仍有可以改进的空间，说明其仍然存在光谱失真和细节缺失问题。HPM获得的融合结果指标CC和STD都较好，但是指标ERGAS偏高，说明其空间分辨率有所提升，但是仍存在光谱畸变现象。NSCT获得的融合结果中CC、STD和UIQI指标都较为不错，证明了NSCT应用于图像融合领域的有效性，但是指标ERGAS和SAM偏高，说明同样存在光谱畸变的问题。本文所提出的算法在各个指标上均优于其余算法，尽管在表2中指标SAM的值没有达到最优，但仅次于最优指标。因此本文所提算法在主观视觉效果和客观指标评价两方面都取得了不错的结果，说明本文算法所得的融合结果在增强空间分辨率的同时保持了较高的光谱分辨率。
图8中的源图像是CNES提供的遥感图像Pléiades.图8(c)和(d)分别为HG算法和IHS算法获得的融合结果，可以看出这两种算法的融合结果细节信息都较为丰富，但是都存在较为严重的光谱失真，屋顶颜色与源多光谱图像相比颜色较浅。图8(e)为SR算法获得的融合结果，和源多光谱图像相比其细节信息有所增加，但是还有提升的空间，而且存在光谱畸变。图8(f)为HPM算法获取的融合结果，可以看出其空间和光谱分辨率都有所提升，但是仍存在光谱失真。图8(g)和(h)分别为NSCT算法和本文算法获得的融合结果，可以看出两幅图像的光谱分辨率和空间分辨率都得到了提升，但单从视觉上观察无法很好地区分哪个融合结果更优，需要借助评价指标。
图8 Pléiades源图像及不同方法获得的融合结果
Fig.8 Fusion results of original Pléiades image and different methods
表3为上述图像的性能比较结果。HG算法获得融合结果的CC和STD指标都较优，但是指标ERGAS和SAM都偏高，说明其空间分辨率高但是图像光谱分辨率低。IHS算法与其相比虽然指标UIQI有所提升，指标ERGAS和SAM的值都有所好转，但是仍偏高，说明其光谱信息损失较大。SR算法获得的融合结果中CC和UIQI指标都优于前两种算法，但ERGAS值偏高，说明该方法虽然整体质量不错但是存在光谱扭曲。HPM算法获取的融合结果的CC和SAM两项指标优于其他算法，证明其细节信息丰富，但是ERGAS值仍然偏高，有很大的提升空间。观察NSCT算法和本文算法获得融合结果的指标可以发现两者很接近，但是本文算法更胜一筹，而且SAM指标明显低于NSCT算法，说明融合结果中保留了更多的光谱信息。整体来看，本文算法的CC和SAM两个指标较好，且ERGAS、STD和UIQI这3个指标明显优于其他算法，说明本文算法可以同时拥有丰富的光谱信息和较高的空间分辨率。
表3 Pléiades源图像融合结果的性能比较
Table 3 Performance comparison of the original image
fusion results of PLéiades
MethodCCERGASSAMSTDUIQIHG0.874 419.682 99.164 30.045 20.702 7IHS0.802 315.571 87.668 40.051 90.704 4SR0.877 418.954 66.402 10.041 80.765 9HPM0.913 412.322 41.811 10.040 30.732 9NSCT0.899 76.926 62.496 90.039 60.892 1proposed0.897 96.561 61.885 90.038 30.892 5
4 结束语
为了解决遥感图像融合步骤中细节信息缺失和光谱畸变问题，提出了基于卷积稀疏表示和NSCT变换的遥感图像融合。首先利用卷积稀疏表示训练得到高分辨率的字典滤波器和低分辨率的字典滤波器，之后用低分辨的字典滤波器对低分辨的图像进行稀疏表示，获得特征映射，最后用特征映射和高分辨率的字典滤波器完成图像的超分辨率。在融合步骤中，先对超分辨率的图像和全色图像进行NSCT变换，得到各自的高低频子带信息。为了在融合结果中保留更多细节信息的同时减少光谱畸变，因此根据高频和低频子带不同的特点，在高频部分采用了局部方差和加权平均结合的融合规则，低频部分采用了稀疏系数取大的融合规则，最后将得到的新高低频子带进行逆变换得到融合结果。实验证明，该方法在提高融合结果空间分辨率的同时也保留了更多的光谱信息。
[1] TORO-GARAY G H,MEDINA-DAZA R J,et al.Fusion of worldview2 images using contourlet,curvelet and ridgelet transforms for edge enhancement[J].Revista Facultad de Ingeniería Universidad de Antioquia,2017(85):8-17.
[2] DA CUNHA A L,ZHOU J P,DO M N.The nonsubsampled contourlet transform:theory,design,and applications[J].IEEE Transactions on Image Processing,2006,15(10):3089-3101.
[3] XUE X R,XIANG F,WANG H F.A parallel fusion method of remote sensing image based on NSCT[J].Компьютерная оптика,2019,43(1):123-131.
[4] ZEILER M D,KRISHNAN D,TAYLOR G W,et al.Deconvolutional networks[J].Computer Vision & Pattern Recognition,2010,2528-2535.
[5] ZHANG K,WANG M,YANG S Y,et al.Convolution structure sparse coding for fusion of panchromatic and multispectral images[J].Journal of Technology & Science,2019,57(2):1117-1130.
[6] WU H L,ZHAO S Z,ZHANG J M,et al.Remote sensing image sharpening by integrating multispectral image super-resolution and convolutional sparse representation fusion[J].IEEE Access,2019(7):46562-46574.
[7] LI H,QIU H,YU Z,et al.Infrared and visible image fusion scheme based on NSCT and low-level visual features[J].Infrared Physics & Technology,2016,76(5):174-184.
[8] WOHLBERG B.Efficient algorithms for convolutional sparse representations[J].IEEE Transactions on Image Processing,2016,25(1):301-315.
[9] RESTAINO R,VIVONE G,DALLA MURA M,et al.Fusion of multi spectral and panchromatic images based on morphological operators[J].IEEE Transactions on Image Processing,2016,25(6):2882-2895.
[10] MUFIT C,ABDULKADIR T.Intensity-hue-saturationbased-based image fusion using iterative linear regression[J].Journal of Applied Remote Sensing,2016,10(4):045019.
[11] VICINANZA M R,RESTAINO R,VIVONE G,et al.A pansharpening method based on the sparse representation of injected details[J].IEEE Geoscience and Remote Sensing Letters,2015,12(1):180-184.
[12] VIVONE G,ALPARONE L,CHANUSSOT J,et al.A critical comparison among Pan sharpening algorithms[J].IEEE Transactions on Geoscience and Remote Sensing,2015,53(5):2565-2586.
[13] YUHAS R,GOETZ A,BOARDMAN J.Discrimination among semi-arid landscape endmembers using the Spectral Angle Mapper (SAM) algorithm[J].JPL Airborne Geosci,1992,1:147-149.
[14] WANG Z,BOVIK A C.A universal image quality index[J].IEEE Signal Processing Letters,2002,9(3):81-84.