基于自适应注入模型的遥感图像融合方法

2020-01-02杨勇卢航远黄淑英涂伟李露奕

北京航空航天大学学报 2019年12期

杨勇，卢航远，黄淑英，涂伟，李露奕

（1.江西财经大学信息管理学院，南昌330032； 2.江西财经大学软件与物联网工程学院，南昌330032）

遥感的目的是通过获取卫星的光谱测量，提取有关地球表面结构和内容的信息［1］。由于传感器技术的限制，同时在空间域和光谱域上获取高分辨率难以实现，因此通常利用图像处理技术融合现有的2类遥感图像，即高空间分辨率、低光谱分辨率的全色（Panchromatic，PAN）图像和低空间分辨率、高光谱分辨率的多光谱（Multispectral，MS）图像。对于在同一区域同时获得的PAN图像和MS图像，通过适当的算法，将这些数据结合起来，生成具有高空间分辨率和高光谱分辨率的图像，这一过程也称为PAN锐化［2］。遥感图像融合广泛应用于地理、军事、农业、生态环境等领域。

现有的传统PAN锐化方法主要分为2类：一类是成分替代法（Component Substitute，CS），即将MS图像的成分替换成PAN图像成分，如亮度-色度-饱和度变换法（Intensity-Hue-Saturation，IHS）［3］、主成分分析法（Principle Component Analysis，PCA）［4］、基于施密特正交化方法（Gram-Schmidt，GS）［5］等。基于CS方法的融合结果具有较高的空间分辨率，但其光谱容易失真。另一类是多尺度分析法（Multi-Resolution Analysis，MRA），即对PAN图像进行多尺度分解后将提取的细节加入到MS图像中。基于MRA的方法和基于CS的方法的主要区别在于如何提取空间细节［6］。对于MRA方法，通过PAN图像与其低通滤波之后的图像之间的差异来获得空间结构信息。传统的MRA方法主要包括基于金字塔分解或者小波变换的方法，如离散小波变换［7］、àtrous小波变换［8］、非抽样的轮廓波变换［9］、非抽样的剪切波变换［10］等。MRA方法可以较好地保存光谱信息，但可能会损失部分空间信息，因此文献［11］提出了一种基于引导滤波多尺度分解的方法，效果有所改进。最新的理论研究提出了稀疏表示（Sparse Representation，SR）的方法［12-13］，以及非线性的途径如卷积神经网络［14］等，取得了较好的效果。但是，这些算法需要大量的计算，运行时间较长，且需要大量的样本。由于遥感图像样本数量的限制，以及算法效率的要求，本文综合考虑了上述几种方法，对基于MRA方法的注入模型进行优化，以实现融合时对光谱信息与细节信息的双保真。

基于MRA方法中关键的一步在于如何提取PAN图像细节。一种较为先进的方法是利用线性时不变的滤波器去匹配MS传感器的点传播函数（Point Spread Function，PSF）。Aiazzi等［15］提出使用高斯滤波器去匹配MS传感器的调制传输函数，但需要使用MS传感器的出厂信息，而该信息不易获得。另一种较为可靠的方法是从现有图像中估计滤波器以模拟MS传感器的PSF。在此框架下，一幅低空间分辨率的MS（Low Resolution MS，LRMS）图像，可以通过滤波器对一幅高空间分辨率的MS（High Resolution MS，HRMS）图像进行空间滤波得到。该滤波器的冲击响应函数可模拟MS传感器的PSF，且该滤波器通常有与高斯模型相似的形状［2，15］。使用该滤波器卷积PAN图像，所得的高频分量即为HRMS图像中缺失的细节成分，且与MS图像高线性相关。此外，综合考虑光谱的保真和细节的注入，本文提出了一种基于相关性分析的高斯滤波估计和自适应细节注入的遥感图像融合方法。首先，通过引导滤波多尺度分解的注入模型得到初始融合图像；其次，用高斯滤波模拟MS传感器的PSF，并用估计的高斯滤波器卷积PAN图像得到优化的注入细节；然后，递归计算自适应的注入系数，使光谱保真与细节注入达到联合最佳；最后，通过一种新的边缘保持策略，得到融合图像。

本文主要创新点如下：①提出了一种基于高斯滤波估计的细节提取方法，使注入的细节得到优化；②提出了一种综合考虑光谱信息与细节信息的自适应注入系数模型，实现融合图像对光谱信息与细节信息的双保真；③提出了一种减少光谱信息损失的边缘保持模型，在增强边缘信息的同时保留光谱信息。

1 相关工作

1.1 注入模型的一般框架

基于MRA方法的注入模型是指通过滤波提取PAN图像细节，再利用注入系数将其注入到MS图像中的一种融合方法。该方法可根据实际需要，组合不同的融合技术，以发挥不同融合技术的优点。基于传统的MRA方法注入模型框架如图1所示，该注入模型的一般表示形式如下［16］：

图1 注入模型一般框架Fig.1 Basic framework of injection model

对于注入系数Gk，文献［16］通过迭代算法，提出基于回归模型的通道间权重系数；文献［17］提出基于像素通道间比例的权重系数；文献［18］通过高通调制（High Pass Modulation，HPM）方法，用MS图像和PAN图像的乘法组合作为注入系数。其中，基于通道间比例的注入系数能更好地保持光谱信息，融合效果较好，通道间注入系数定义如下：

1.2 线性加权法提取I分量

由于注入细节的同时不希望改变色度与饱和度，因此图像融合过程中往往仅对图像的亮度分量进行处理。MS图像的亮度分量I可由如下公式计算：

式中：αk为通道的组合系数。

文献［19-20］提出I分量可自适应地由各通道线性组合表示，即通过求解以下优化问题得到组合系数：

式中：P为PAN图像。

在得到MS图像的I分量之后，将P与式（3）的I进行直方图匹配，得到的PI可由如下公式计算：

式中：mean（P）为P的均值；mean（I）为I分量的均值；std（I）为I分量的标准差；std（P）为P的标准差。

1.3 引导滤波

引导滤波最早由文献［21］提出，其可以保存输入图像的主要信息，同时获得引导图像的变化趋势。将PI作为输入图像，I作为引导图像，PI通过I引导滤波产生PL，获得了与MS图像线性相关的细节信息。引导滤波可由如下简化公式而得：

式中：GF（·）表示引导滤波函数。

多尺度引导滤波［22］可表示为

2 图像融合方法

目前，遥感图像融合注入模型主要存在2类问题：因PAN图像与MS图像相关性不高导致的光谱失真；细节注入过多或者注入不足。针对上述问题，本文提出了一种基于高斯滤波估计的细节提取及自适应的注入系数优化方法。本文融合方法框架如图2所示。

首先，根据式（1），将上采样后的MS图像与PAN图像通过直方图匹配、引导滤波等传统注入模型得到初始融合图像F（1）；其次，进行细节优化，即用多尺度高斯滤波器滤波该初始融合图像以模拟MS传感器的特性，将得到的高斯滤波估计滤波PAN图像，得到优化后的细节；然后，优化注入效益，即综合考虑光谱信息与细节信息计算自适应的注入量系数，具体流程如图2上虚线框所示，详细算法见2.2节；最后，通过一种自适应的边缘保持权重矩阵保持边缘信息，具体流程如图2下虚线框所示，详细算法见2.3节。将优化后的细节与注入系数及边缘权重系数相乘以后，注入到上采样后的MS图像中，得到最终的融合图像。

2.1 高斯滤波估计及细节提取

文献［15］提出使用高斯滤波器匹配MS传感器的调制传输函数，即HRMS图像经过MS传感器之后成为LRMS图像，使用高斯滤波模拟该过程，所得的高斯滤波器即为符合MS传感器PSF的滤波估计。使用估计的滤波器卷积PAN图像，提取的细节成分即为低分辨率MS图像所缺失的细节。文献［15］通过实验证明了该滤波估计有较强的鲁棒性，即使在初始融合图像中加入白噪声，也不影响细节的提取。另外，文献［17］提出了二步骤多尺度分解的方法来精细化PAN锐化。受此启发，本文提出使用多尺度高斯滤波模拟MS传感器的PSF，即对初始融合图像F（1）提取亮度分量，记为I1，并迭代地进行高斯滤波，具体可表示为

图2 本文融合方法框架Fig.2 Framework of proposed fusion method

为了获取最佳高斯滤波估计，计算式（8）所得的结果Ii1与上采样后的MS图像的亮度分量I的相关系数，相关系数最高时，估计的高斯滤波最准确。相关系数计算如下：

式中：cov（·）表示协方差函数。

随着滤波层次的增加，滤波后的图像与上采样的MS图像越来越接近，达到极大值之后，随着滤波层次继续增加，滤波后图像与上采样的MS图像偏离越来越大。通过迭代计算式（8）和式（9）的值，取得式（9）最大时的迭代次数m，即为所估计的高斯滤波器Hm，表示如下：

Hm即表示一幅锐化图像经MS传感器模糊化的PSF最佳匹配，将滤波器Hm作用于PAN图像，所提取的细节成分也应符合MS传感器PSF，这样就减少了光谱失真。细节提取D定义如下：

结合式（1）、式（7）和式（11），细节优化后，第k通道融合图像可更新为

2.2 自适应注入量系数

取得优化的注入细节之后，对注入量进行优化。不同的图像有不同的结构特征和光谱信息，从而对细节注入量有不同的要求，直接等比例注入细节容易导致注入过多从而引起光谱失真，或者注入不足导致图像模糊，因此需要通过一个注入量系数g来平衡光谱信息与细节信息。加入该系数后，结合式（12），新的融合图像可以定义为

需要对不同结构特征的遥感图像确定不同的注入量系数g。文献［23］认为注入模型中光谱的保真与空间分辨率要求需要保持平衡，可以将两者的评价指标进行加权求和作为融合图像指标，权重表示重视程度。设定细节注入量系数g的初值g和步长r，对于g0≤g≤1，得到相应的融合后图像并计算加权评价指标，取得使该评价指标最高的gmax作为最终的注入量系数。对于评价指标，本文提出了一种基于相关性分析的光谱信息与空间信息综合评价体系，将光谱信息与空间信息的综合评价指标Q定义为

式中：ESP和EHF分别为光谱信息评价指标和空间信息评价指标；α为两者的权重。

对于一幅由式（13）融合后的图像F（3）

k，其亮度分量记为I3，则ESP用融合后的图像与上采样的MS图像各通道k之间相关系数的平均值来表示，EHF则用I3与PAN图像的线性相关性表示，定义为

文献［6］表明，在注入成分与注入量系数g的优化确定上，遵循的前提是MS图像的亮度成分与PAN图像相关性越大，融合效果越好。即若两者相关性高，注入细节的同时能较好地保存光谱信息，此时可以适当增加空间信息的权重，减少光谱信息的权重；反之，若相关性低，则注入的细节容易引起光谱失真，应当增加光谱信息权重而减少空间信息权重。设初值g0产生的融合图像对应的亮度成分记为I0，则式（14）中的权重α可由下式确定：

平方项的作用为：①避免I0与PI高相关时光谱权重过小；②放大不同图像空间信息之间的区别。

结合式（13）和式（18），使用优化后的注入量系数，新的融合图像可更新为

2.3 边缘信息保护

为保持边缘纹理信息，Leung等［20］提出利用边缘信息作为各通道融合权重，Yang等［11，24］提出利用改进的边缘信息作为注入各像素的权重。其中，改进的边缘信息计算PAN图像的边缘检测矩阵WP，以及上采样后的MS图像各通道的边缘检测矩阵W，并将两者进行加权和作为融合图像的边缘信息权重。该方法在边缘信息保护上起到了较好的效果。第k通道的边缘保持权重矩阵Ek定义如下［20］：

式中：βk为各通道的加权系数，可由式（21）计算得到［11］；边缘检测矩阵［19］可由式（22）得到。

式中：A为某图像；ΔA为图像A的梯度；λ和ε为调制参数。

式（21）获取各通道的边缘信息权重，在细节上有更好的表现，但作为注入系数时改变了各通道间的比例，光谱信息略有损失。为了尽量保持原有光谱信息，同时实现边缘信息的保持，本文提出一种新的边缘保持权重矩阵。针对式（21），将求取各通道的边缘检测矩阵Wk替换为求取MS图像亮度成分I的边缘检测矩阵WI，并计算PAN图像的边缘检测矩阵WP与WI加权和，以综合考虑PAN图像和MS图像的边缘信息。边缘保持权重矩阵E可表示为

式中：WP和WI计算方法同式（22）；自适应的权重因子βnew由式（24）计算，并由式（25）归一化得到。

根据式（23）得到的边缘保持权重矩阵E，再结合式（19），可得到最终融合图像，表示如下：

与式（20）的边缘保持权重矩阵相比，本文得到的边缘保持权重矩阵E仅对细节成分D产生权重变化，以尽可能增强边缘信息，并不改变通道间的比例，从而减少光谱损失。

3 实验结果分析

3.1 实验设置

为评估本文方法的有效性，使用来自不同遥感卫星的四大数据集，包括Quick Bird、IKONOS、WorldView-2及pleiades，这些遥感图像在空间分辨率、光谱波长等方面都有不同的特征，且都包括红、绿、蓝、近红外4个通道。

参数设置中，设置注入量系数初值g0=0.1，步长r=0.05，以综合考虑计算效率与计算精度；高斯滤波器窗口设置为5×5，参考文献［25］中的默认设置；式（22）中设置参数λ=10-9，ε=10-10，同样参考文献［19］中的默认设置。

本文进行了2种类型的实验。一种是对真实图像进行下采样之后的仿真图像，也称为有参考遥感图像融合实验，即遵循Wald协议［26］，原图像作为参考图像，仿真图像包含四大数据集共计80组图片。另一种是真实图像的实验，也称为无参考遥感图像融合实验，真实图像包括四大数据集共计50组图片。

用于对比的方法包括最新的几种方法，其中基于MRA模型的方法有基于调制传递函数模型的多尺度分析方法（CBD）［15］、加性小波亮度比的方法（Additive Wavelet Luminance Proportion，AWLP）［27］、基于双边滤波亮度比的方法（Bilateral Filter Luminance Proportional，BFLP）［28］；将MRA和CS两种模型相结合的方法有基于多尺度引导滤波和改进的成分替代法（Improved IHS and Multi-Scale Guided Filter，IMG）［11］；将注入模型与其他算法相结合的方法有基于抠图模型和多尺度变换的方法（Matting Model and Multi-scale Transform，MMMT）［24］、基于自适应光谱与亮度调制的方法（Adaptive Spectral Intensity Modulation Pansharpening，ASIMP）［25］。对这些方法分别进行主观评价与客观评价，也即定性和定量评价。

3.2 质量评价

评价包括肉眼视觉的主观评价和客观评价指标的定量评价。常用的有参考图像的客观评价指标如下：

1）相关系数（Correlation Coefficient，CC）［12］。主要用于评价融合图像与参考图像的空间相似度，其值区间为［0，1］，越接近1表示2幅图像越相近，融合效果越好。

2）通用图像质量指数（Universal Image Quality Indexes，UIQI）［29］。UIQI联合相关度损耗、亮度失真和对比度失真3个因子来评价融合图像的效果，值区间为［0，1］。UIQI的值越接近1，说明结构失真越小。

3）均方根误差（Root Mean Square Error，RMSE）［30］。RMSE计算融合图像与参考图像各像素之间均方根误差，用于表示两者的差异，其理想值为0。

4）相对平均光谱误差（Relative Average Spectral Error，RASE）［31］。RASE反映融合图像在光谱方面的误差，其理想值为0。

5）光谱角映射（Spectral Angle Mapper，SAM）［32］。SAM用于度量2幅图像之间的光谱信息接近程度，其值越小，光谱失真越小，理想值为0。

6）全局相对光谱损失（ERGAS）［33］。ERGAS用于评估融合图像的空间和光谱质量，其值越小说明总体性能越好，理想值为0。

常用的无参考图像的客观评价指标包括Ds、Dλ和QNR［34］。其中，QNR值由Ds和Dλ共同计算而得，并基于UIQI测度评价融合图像和参考图像间的局部相关性、亮度和对比度。Ds和Dλ越小越好，理想值为0，而QNR值越大越好，理想值为1。

3.3 算法性能分析

为了验证2.1节～2.3节提出的优化方法的有效性，定量说明本文方法所带来的性能提升，分别计算各优化步骤的平均客观评价指标。4个优化步骤对应的融合图像表示为：F1表示初始融合图像；F2表示经高斯滤波优化细节后图像；F3表示经自适应优化注入量系数后的图像；F4表示加入自适应边缘保持后的图像。并分别计算80组图像融合结果的平均客观评价指标。这4个步骤是依次递进的过程。为了更直观地展示各指标的相对变化，各客观评价指标分别进行最大值归一化处理，结果如图3所示。可知，对于初始融合图像F1，注入细节是通过PAN图像与MS图像的I分量进行直方图匹配再由引导滤波得到，细节存在与MS图像线性非相关的问题，且未考虑注入量，从而导致融合结果与参考图像的光谱相关度和空间相关度较低；通过估计的高斯滤波器优化细节后，所得的细节与MS图像有更高的相关性，各项指标明显提升；而优化注入量系数后，有效防止了注入过多或注入不足，所有的指标进一步提升；进行边缘保护后，改进的边缘保持方法在保护边缘信息的同时，减少了光谱失真。随着F2、F3、F4优化步骤的依次进行，总体指标也越来越优。其中SAM指标基本保持平稳，随着优化的进行，SAM指标先小幅上升再小幅降低，说明本文方法对光谱信息损失的影响很小。

图3 优化过程平均客观评价指标对比Fig.3 Average objective evaluation indicator comparison of optimization steps

如2.3节所述，本文提出的一种新的自适应边缘保持方法是对文献［11］的改进。为了验证本文方法的有效性，在保持前3个步骤不变的前提下，对改进前［11］与改进后的方法在80组遥感图像上进行平均客观评价指标对比，指标同样进行最大值归一化处理，结果如图4所示。可知，改进后的边缘保持方法在所有指标上都有所提升，其中，在RASE、RMSE及ERGAS等指标上提升效果明显，可知本文方法能有效地提升图像融合的质量，减少光谱损失。

图4 边缘保持方法改进前后平均客观评价指标对比Fig.4 Comparison of average objective evaluation indicators before and after improvement of edge preserving method

3.4 仿真实验

本文采用3组来自不同卫星的遥感图像作为展示，分别进行主观效果和客观评价指标的评价，然后对80组图像测量平均客观评价指标。

第1组图像来自QuickBird数据集，参考图像大小为256×256，经不同方法融合后主观效果如图5所示。小红色矩形框区域经放大后在大矩形框中显示。可知，AWLP方法、CBD方法、IMG方法森林区域存在明显光谱失真，与参考图像相比，颜色偏亮；BFLP方法、MMMT方法有稍微的光谱失真，且存在过度注入问题，引入了不必要的噪声。ASIMP方法与本文方法在海滩细节与森林的光谱信息上与参考图像较为接近，难以用肉眼区别，用表1所示的客观评价指标加以区分。表中：黑体数据表示最优值，下划线数据表示次优值。可知，本文方法除了在SAM 指标上比BFLP方法略低以外，在所有其他指标上都达到最优。

图5 QuickBird数据集遥感图像融合结果Fig.5 Fusion results of remote sensing images from QuickBird dataset

第2组图像来自World View-2数据集，参考图像大小为256×256，经不同方法融合后主观效果如图6所示，客观评价指标如表2所示。由图中红色放大矩形框可知，IMG方法有较为明显的光谱真失，亮度偏亮；AWLP方法、BFLP方法和CBD方法有轻微的光谱失真，表现为屋顶区域颜色更浅，且在屋顶、中间白色地面有较为明显的噪声；MMMT方法、ASIMP方法与本文方法能较好地保持光谱信息，但MMMT方法在边缘上更加模糊，而本文方法与参考图像不论在空间信息上还是光谱信息上都更为接近。为更好地区分融合效果，可通过表2查看融合结果客观评价指标，可知，本文方法在所有指标上都达到最优。

第3组图像来自pleiades数据集，经不同方法融合后主观效果如图7所示，客观评价指标如表3所示。从放大的红色矩形区域可知，AWLP方法、BFLP方法和CBD方法存在较为明显的光谱失真，且引入了较多的噪声；IMG方法、MMMT方法和ASIMP方法有较轻微的光谱失真，引入少量的噪声；而本文方法能在保持空间细节的同时，较好地保留光谱信息，色彩更加真实。对于难以用肉眼明显区分的方法，可通过表3查看融合结果客观评价指标，可知，本文方法除了在UIQI指标上比ASIMP略低以外，在所有其他指标上都达到最优。

表1 图5对应的遥感图像融合结果客观评价指标Table 1 Objective evaluation indicators of remote sensing image fusion results in Fig.5

图6 WorldView-2数据集遥感图像融合结果Fig.6 Fusion results of remote sensing images from WorldView-2 dataset

表2 图6对应的遥感图像融合结果客观评价指标Table 2 Objective evaluation indicators of remote sensing image fusion results in Fig.6

为了更客观地评价各方法的性能，对来自QuickBird、IKONOS、WorldView-2及pleiades四大数据集的共计80组有参考遥感图像融合后计算平均客观评价指标，结果如表4所示。可知，本文方法在所有的平均客观评价指标上达到了最优。

图7 pleiades数据集遥感图像融合结果Fig.7 Fusion results of remote sensing images from pleiades dataset

表3 图7对应的遥感图像融合结果客观评价指标Table 3 Objective evaluation indicators of remote sensing image fusion results in Fig.7

表4 80组遥感图像仿真实验的平均客观评价指标Table 4 Average objective evaluation indicators of simulation experiments from 80 groups of remote sensing images

3.5 真实图像实验

真实图像的实验采用2组图像进行主观和客观评价，并对所有来自四大数据集的50组真实图像计算平均客观评价指标。

第1组图像来自IKONOS数据集，图像大小为780×608。不同方法融合后的主观效果如图8所示，相应的客观评价指标如表5所示。由图8可知，BFLP方法和MMMT方法在橙色屋顶区域有明显光谱失真，其他方法与本文方法难以用肉眼直接判别，可通过表5中图8的客观评价指标来评价，可知，本文方法除了在Ds指标上略低于ASIMP方法以外，其他指标都达到最优。

第2组图像来自QuickBird数据集，图像大小为512×512，融合后视觉效果如图9所示。可知，AWLP方法、CBD方法和IMG方法亮度偏亮，存在过度注入的问题，其他方法难以用肉眼直接观察，可由表6中图9的客观评价指标测量，可知，本文方法在所有客观评价指标上都达到最优。

为更加客观评价不同方法的有效性，本文对来自四大数据集的50组真实图像进行实验，并测算平均客观评价指标，如表7所示。可知，本文方法所有指标都达到最优。