本发明涉及遥感影像处理领域，尤其是一种能有效提高多光谱影像空间分辨率和光谱分辨率的、计算复杂度适中、方向自适应性强、光谱特性保真度高、基于log-gabor变换与方向区域熵的引导滤波遥感影像融合方法。
背景技术：
遥感对地观测的根本任务是获取地理空间信息。地理空间信息包括几何信息和属性信息两方面，其中，几何信息包括观测对象的空间位置、大小和形状等，属性信息则包括观测对象的种类、对不同电磁波谱段的反射、透射和吸收特性。只有将几何信息和属性信息有机结合才能实现全面的对地观测。然而，在实际的对地观测中，由于传感器成像机制和参数的差异性以及性能局限，单一传感器尚无法获得空间分辨率和光谱分辨率均较高的遥感影像。其中，多光谱影像(multispectral,ms)具有较高的光谱分辨率和较低的空间分辨率，并且具有多个波段；全色影像(panchromatic,pan)具有较高的空间分辨率和较低的光谱分辨率，并且只有一个波段。于是，遥感图像融合便应运而生，其目的是在保持光谱分辨率的前提下，尽可能地增加影像的空间分辨率。融合后的遥感影像既具有较高的空间分辨率，又具有较高的光谱分辨率，进而利用源图像之间的优势互补实现有效的对地观测。作为图像融合领域中的重要分支之一，遥感影像融合技术目前已被广泛应用在军事战略防御、目标侦测跟踪、测绘制图、森林覆盖率分析、自然灾害检测和天气预报等诸多领域，具有非常重要的研究价值与现实意义。
早期的遥感影像融合技术主要采用空间域像素级融合策略，提出了加权法、主成分分析法、亮度色度饱和度变换方法、gram-schmidt方法等融合技术。在多光谱图像出现以后，为了将原始的ihs变换应用到3个以上波段的影像融合中，leung等人提出了广义ihs算法和自适应的ihs算法。但是，这些方法会导致融合后的影像出现严重的光谱扭曲现象或视觉质量降质，具有明显的局限性。为此，burt等人提出了基于拉普拉斯金字塔变换的图像融合算法，开创了变换域图像融合技术的先河，使图像融合技术从空间域跨越到了变换域。其基本思想是将遥感影像从空间域映射到变换域，得到一个近似于原始影像的低频子带和多个反映影像细节信息的不同尺度、不同方向的高频子带，然后根据一定的融合规则对得到的低频和高频子带分别进行融合处理，从而在变换域内整合不同图像之间的冗余与互补信息。此类方法的关键在于，选取一种恰当的变换，有效地从具有较高空间分辨率、较低光谱分辨率的影像中提取出清晰的纹理细节信息，并将其注入空间分辨率较低、光谱分辨率较高的影像中。为了达到这一目的，pajares等人将小波变换引入融合过程，提出了一种基于小波变换的遥感影像融合方法。但是，由于小波变换不具备平移不变性，该方法产生的融合影像中往往呈现吉布斯效应。于是，伴随着多尺度几何分析领域的快速发展，ridgelet、curvelet、contourlet、shearlet等多尺度变换逐渐被应用于图像融合领域。尽管这些变换在一定程度上克服了融合过程引起的吉布斯现象，可是其计算效率均较低，无法满足实际应用的需求。同时，现有的遥感影像融合技术在纹理细节注入过程中，大多仅依据变换系数的能量对比关系，简单地选择能量较大的变换系数进行注入，却忽视了待融合图像的边缘纹理特征的保持，以致往往在结果图像中引入伪边缘并干扰遥感影像的光谱特性。这将不可避免地影响人类对地观测的精度，甚至产生错误的决策。
总体来看，现有技术方案仍存在计算复杂度高、方向自适应性弱、光谱特征保持不佳的缺点。并且，由于遥感影像的统计特性不同于自然图像，这些方案尚缺乏对遥感影像光谱特征的有效度量和约束。
技术实现要素：
本发明是为了解决现有技术所存在的上述技术问题，提供一种计算复杂度适中、方向自适应性强、光谱特性保真度高、基于log-gabor变换与方向区域熵的引导滤波遥感影像融合方法。
本发明的技术解决方案是：一种基于log-gabor变换与方向区域熵的引导滤波遥感影像融合方法，其特征在于按如下步骤进行：
步骤1.输入待融合的多光谱影像m和全色影像p，其空间分辨率分别为wm×hm、wp×hp；
步骤2.利用双三次插值对多光谱影像m进行上采样，使其空间分辨率达到wp×hp；
步骤3.对多光谱影像m进行ihs变换，提取出变换后的强度分量i，并保持h分量与s分量不变；
步骤4.利用log-gabor变换对强度分量i进行方向数量为2、尺度为2的分解，从而得到了4个子带，令其分别为l1、l2、l3、l4；
步骤5.构建决策模板mapdecision，用来监督细节注入的过程；
步骤5.1根据公式(1)和公式(2)的定义，计算预决策模板mapa和mapb；
mapb(i,j)＝ω×sumh(i,j)+(1-ω)×suml(i,j)(2)
所述lk∈{l1,l2,l3,l4}，i∈{1,2,3，…，wp}，j∈{1,2，3,…,hp}，lk(i,j)表示第k个子带lk中位于坐标(i,j)处的变换系数值，mapa(i,j)表示预决策模板mapa中位于坐标(i,j)处的系数值，ω表示预设的权重，mapb(i,j)表示预决策模板mapb中位于坐标(i，j)处的系数值，sumh(i，j)表示每个尺度下第一层子带中位于坐标(i，j)处的系数之和，suml(i,j)表示每个尺度下、除第一层子带以外的所有其他子带中位于坐标(i，j)处的系数之和；
步骤5.2根据公式(3)的定义，利用预决策模板mapa和mapb计算得到决策模板mapdecision；
所述mapdecision(i，j)表示模板mapdecision中位于坐标(i,j)处的系数值；
步骤6.根据公式(4)的定义，对全色影像p进行局部均值滤波，得到平滑后的全色影像
所述p(i，j)表示全色影像p中位于坐标(i，j)处的像素值，
步骤7.根据公式(5)的定义，计算强度分量i的每个log-gabor变换子带的方向区域熵值矩阵；
所述k∈{1,2,3,4}，lk(i+l,j+c)表示第k个子带lk中位于坐标(i+l,j+c)处的变换系数值，enk表示强度分量i的第k个log-gabor变换子带的方向区域熵值矩阵，enk(i,j)表示enk中位于坐标(i,j)处的元素值；
步骤8.根据公式(6)的定义，计算强度分量i的边缘引导矩阵ζ；
所述ζ(i,j)表示边缘引导矩阵ζ中位于坐标(i,j)处的元素值；
步骤9.根据公式(7)和公式(8)的定义，计算基于方向区域熵的边缘引导模型的系数矩阵a和b；
所述a(i,j)表示系数矩阵a中位于坐标(i,j)处的元素值，b(i,j)表示系数矩阵b中位于坐标(i,j)处的元素值，i(i,j)表示强度分量i中位于坐标(i,j)处的像素值，μ(i,j)和δ2(i,j)分别表示强度分量i的处于正方形窗口ωi,j中的像素均值和方差，η表示预设的正则化参数；
步骤10.根据公式(9)和公式(10)的定义，计算系数矩阵a和b进行局部平均滤波，得到平滑后的系数矩阵
所述
步骤11.根据公式(11)的定义，计算基于方向区域熵的边缘引导模型的估计矩阵
所述
步骤12.根据公式(12)的定义，计算待注入的边缘纹理细节特征图d；
所述d(i，j)表示d中位于坐标(i,j)处的元素值；
步骤13.利用前向差分法，计算强度分量i在每个像素处的梯度，得到每个像素处的梯度模长矩阵ig；
步骤14.根据公式(13)的定义，计算梯度模长矩阵ig的梯度能量熵eu；
所述maxθig(i,j)和minθig(i,j)分别表示ig中以坐标(i,j)为中心、半径为7的正方形窗口内的最大梯度模长和最小梯度模长，pg(h)表示梯度模长为h的像素所占的比重，eu(i，j)表示位于坐标(i，j)处的梯度能量熵值；
步骤15.根据公式(14)和公式(15)的定义，计算自适应放大系数模板α和自适应缩小系数模板β；
β(i,j)＝α(i，j)-1(15)
所述mid(eu)表示梯度能量熵eu的中值，σ表示梯度能量熵eu的标准差，α(i,j)表示自适应放大系数模板α中位于坐标(i,j)处的元素值，β(i,j)表示自适应缩小系数模板β中位于坐标(i,j)处的元素值；
步骤16.根据公式(16)的定义，计算融合后的i分量ifused；
所述ifused(i,j)表示ifused中位于坐标(i,j)处的像素值，δ表示预设的经验阈值；
步骤17.将h分量、s分量与ifused分量进行ihs逆变换，得到融合后的多光谱影像mfused，输出mfused。
与现有技术相比，本发明具有三个方面优点：第一，较之传统的多尺度几何分析方法，log-gabor变换更加符合人类的视觉认知特性，能够提取出更完整的高频信息，且其计算复杂度适中；第二，利用方向区域熵充分发掘了log-gabor变换系数所具有的多尺度、多方向性特点，构建了新的边缘引导矩阵及其自适应引导滤波模型，更加充分地提取出全色影像pan中丰富的边缘和纹理；第三，引进基于监督决策模板的纹理注入监督机制，在兼顾考虑pan影像和ms影像空间信息的同时，既有效地保持了ms影像的光谱特征，避免了光谱信息的扭曲和噪声斑块等信息的干扰，又能保证融合过程向低分辨率影像中注入更丰富的纹理细节信息。因此，本发明具有计算复杂度适中、方向自适应性强、光谱特性保真度高的优点。
具体实施方式
本发明的一种基于log-gabor变换与方向区域熵的引导滤波遥感影像融合方法，按照如下步骤进行：
步骤1.输入待融合的多光谱影像m和全色影像p，其空间分辨率分别为wm×hm、wp×hp；
步骤2.利用双三次插值对多光谱影像m进行上采样，使其空间分辨率达到wp×hp；
步骤3.对多光谱影像m进行ihs变换，提取出变换后的强度分量i，并保持h分量与s分量不变；
步骤4.利用log-gabor变换对强度分量i进行方向数量为2、尺度为2的分解，从而得到了4个子带，令其分别为l1、l2、l3、l4；
步骤5.构建决策模板mapdecision，用来监督细节注入的过程；
步骤5.1根据公式(1)和公式(2)的定义，计算预决策模板mapa和mapb；
mapb(i,j)＝ω×sumh(i,j)+(1-ω)×suml(i,j)(2)
所述lk∈{l1,l2,l3,l4}，i∈{1,2,3,…,wp}，j∈{1,2,3,…,hp}，lk(i,j)表示第k个子带lk中位于坐标(i,j)处的变换系数值，mapa(i,j)表示预决策模板mapa中位于坐标(i,j)处的系数值，ω表示预设的权重，mapb(i,j)表示预决策模板mapb中位于坐标(i,j)处的系数值，sumh(i，j)表示每个尺度下第一层子带中位于坐标(i，j)处的系数之和，suml(i，j)表示每个尺度下、除第一层子带以外的所有其他子带中位于坐标(i，j)处的系数之和，本实施例中令ω＝0.7；
步骤5.2根据公式(3)的定义，利用预决策模板mapa和mapb计算得到决策模板mapdecision；
所述mapdecision(i，j)表示模板mapdecision中位于坐标(i,j)处的系数值；
步骤6.根据公式(4)的定义，对全色影像p进行局部均值滤波，得到平滑后的全色影像
所述p(i,j)表示全色影像p中位于坐标(i，j)处的像素值，
步骤7.根据公式(5)的定义，计算强度分量i的每个log-gabor变换子带的方向区域熵值矩阵；
所述k∈{1,2，3,4}，lk(i+l,j+c)表示第k个子带lk中位于坐标(i+l,j+c)处的变换系数值，enk表示强度分量i的第k个log-gabor变换子带的方向区域熵值矩阵，enk(i,j)表示enk中位于坐标(i,j)处的元素值；
步骤8.根据公式(6)的定义，计算强度分量i的边缘引导矩阵ζ；
所述ζ(i,j)表示边缘引导矩阵ζ中位于坐标(i,j)处的元素值；
步骤9.根据公式(7)和公式(8)的定义，计算基于方向区域熵的边缘引导模型的系数矩阵a和b；
所述a(i，j)表示系数矩阵a中位于坐标(i，j)处的元素值，b(i，j)表示系数矩阵b中位于坐标(i，j)处的元素值，i(i,j)表示强度分量i中位于坐标(i,j)处的像素值，μ(i,j)和δ2(i,j)分别表示强度分量i的处于正方形窗口ωi,j中的像素均值和方差，η表示预设的正则化参数，本实施例中令η＝0.64；
步骤10.根据公式(9)和公式(10)的定义，计算系数矩阵a和b进行局部平均滤波，得到平滑后的系数矩阵
所述
步骤11.根据公式(11)的定义，计算基于方向区域熵的边缘引导模型的估计矩阵
所述
步骤12.根据公式(12)的定义，计算待注入的边缘纹理细节特征图d；
所述d(i,j)表示d中位于坐标(i,j)处的元素值；
步骤13.利用前向差分法，计算强度分量i在每个像素处的梯度，得到每个像素处的梯度模长矩阵ig；
步骤14.根据公式(13)的定义，计算梯度模长矩阵ig的梯度能量熵eu；
所述maxθig(i,j)和minθig(i,j)分别表示ig中以坐标(i,j)为中心、半径为7的正方形窗口内的最大梯度模长和最小梯度模长，pg(h)表示梯度模长为h的像素所占的比重，eu(i,j)表示位于坐标(i，j)处的梯度能量熵值；
步骤15.根据公式(14)和公式(15)的定义，计算自适应放大系数模板α和自适应缩小系数模板β；
β(i，j)＝α(i，j)-1(15)
所述mid(eu)表示梯度能量熵eu的中值，σ表示梯度能量熵eu的标准差，α(i，j)表示自适应放大系数模板α中位于坐标(i,j)处的元素值，β(i，j)表示自适应缩小系数模板β中位于坐标(i，j)处的元素值；
步骤16.根据公式(16)的定义，计算融合后的i分量ifused；
所述ifused(i，j)表示ifused中位于坐标(i,j)处的像素值，δ表示预设的经验阈值，本实施例中令δ＝9；
步骤17.将h分量、s分量与ifused分量进行ihs逆变换，得到融合后的多光谱影像mfused，输出mfused。
为验证本发明的有效性，以test-rs-1、test-rs-2和test-rs-3三个遥感影像数据集为例进行实验，采用ssim、sam、rase、rmse、ergas、cc、snr作为客观指标对融合结果进行评价，并将本发明的评价结果与ihs算法、传统基于引导滤波的融合算法(gif)、基于方向区域能量和非下采样剪切波的融合算法(nsst+dmd)、基于梯度域加权的融合算法(ctmr)、wica算法、基于自适应字典学习的融合算法(asr)和基于深度学习的融合算法(dlfs)进行比较。
从表1可见，本发明对test-rs-1遥感影像的融合结果取得了最优的ssim、rase、ergas、rmse、cc与snr值。尽管本发明的sam评价结果略逊于asr，但是本发明在其他指标上评价结果明显高于asr算法。这表明，本发明在test-rs-1遥感影像数据集上得到的融合影像含有最丰富的纹理信息与光谱信息。从表2可见，本发明对test-rs-2遥感影像的融合结果取得了最优的ssim、rase、ergas、sam、cc与snr值，只有其rmse评价结果略低于nsst+dmd方法。这表明，本发明有效地提高了test-rs-2遥感影像的空间分辨率，并且能够充分保持遥感影像的光谱信息。从表3可见，本发明在test-rs-3遥感影像上取得的客观评价结果整体优于其他现有技术。综合表1至表3的对比结果可知，通过有效提取遥感影像的方向纹理信息并保持其光谱特性，本发明获得了最佳的遥感影像融合性能，有效提高了融合后影像的客观质量。
表1不同技术在test-rs-1遥感影像数据集上的融合结果定量分析统计表
表2test-rs-2遥感影像数据集上的融合结果定量分析统计表
表3test-rs-3遥感影像数据集上的融合结果定量分析统计表