CN110942451B - 一种无参考图像的遥感图像融合性能评价方法 - Google Patents

Abstract

本发明提供一种无参考图像的遥感图像融合性能评价方法，主要解决在没有高分辨率参考图像的情况下，如何有效地对遥感图像融合性能进行有效评价的问题，以及评价指标的有效性问题。本发明的核心机制是，将图像融合性能分解为光谱保持度、空间平滑度和细节保持度三个维度，通过对融合图像进行子带分解，并对每一子带与参与融合的多光谱图像和全色图像分别进行比对，计算以相关系数为基础的指标，最后进行指标综合，获得最终评价指标。本发明的方法对融合图像性能评价结果与专家经验一致，可以很好地表现融合算法的不同性能方面，运算速度快，适用于遥感图像融合评价。

Description

一种无参考图像的遥感图像融合性能评价方法

技术领域

本发明涉及属于航天光学遥感成像及地面处理系统和数字图像处理技术领域，特别是涉及一种无参考图像的遥感图像融合性能评价方法。

背景技术

卫星遥感技术可以对地球表面进行大规模全覆盖的观测，在测绘、气象、海洋、农业、自然资源调查、灾害监测、国防安全等诸多方面均发挥巨大作用，世界各国都十分重视。由于传感器信噪比的限制，光学遥感卫星所获取图像的光谱分辨率和空间分辨率之间存在矛盾，因此目前光学遥感卫星一般配置一个高分辨率全色波段和多个低分辨率多光谱波段。

遥感图像融合技术的目标是，通过软件算法，将高空间分辨率低光谱分辨率的全色图像与低空间分辨率高光谱分辨率的多光谱图像，合成为一幅可由更高性能的传感器，或同等传感器在较低轨道上获取的图像，使其具有与多光谱图像相同的光谱分辨率，和与全色图像相同的空间分辨率。

根据人们对遥感数据利用的方式，对融合图像评价主要包括两个方面：空间性能和光谱性能。空间性能主要指融合图像在视觉效果上，与原始的多光谱图像相比，是否得到充分的改进；光谱性能则是面向后续的计算机定量处理，表征的是融合图像是否能够完整地保留多光谱图像中包含的光谱信息。因此，融合性能评价包括主观评价和客观评价/定量评价。随着定量遥感的发展，定量评价逐渐成为学界的主流，发展现状是，虽然有多种定量评价方案，但主要侧重于光谱性能的评价，综合评价结果往往与主观评价之间存在较大的冲突。问题的核心在于现有的定量评价方法不能很好地对空间性能进行定量表达。

定量评价方法可分为有参考图像的评价和无参考图像的评价两种。有参考的图像融合算法评价的基本原理，就是将经某种方法得到的融合图像与这一理想参考图像进行比较，进而作为对融合算法的评价。由于理想参考图像并不存在，因此目前主要采用两种方法进行模拟。第一种方法采用其他来源的更高分辨率数据，主要是机载数据，根据目标传感器的特性进行模拟。这种方法成本高，适用范围小，同时，模拟方法也有争议，因此学界主要采用第二种方法——降分辨率方案，即将全色和多光谱图像均按照二者分辨率比例进行退化，对降分辨率的全色和多光谱图像进行融合，从而可将原始多光谱图像作为参考图像。这一方法虽然得到普遍应用，但有研究表明，退化方案会严重影响评价结果；在实践中也发现，在降分辨率方案中得分很高的方法，在原始尺度下融合性能很差。

无参考图像的评价方法则是根据融合图像与参与融合的全色和多光谱图像之间的定量关系，在一定的理论假设下，对融合性能进行定量评价。无参考图像的评价适用性更广泛，稳定性更好，是目前研究的重点，但现有无参考图像在指标构建方面还不太成功，对空间性能和光谱性能两方面的表征都有较大偏差。如何构造与主观评价一致性高的无参考图像融合评价指标和方法，是一个需要解决的问题。

发明内容

为了在没有高分辨率参考图像的情况下，如何有效地对遥感图像融合性能进行有效评价的问题，以及评价指标的有效性问题，本发明提供一种无参考图像的遥感图像融合性能评价方法，本发明将图像融合性能分解为光谱保持度、空间平滑度和细节保持度三个维度，通过对融合图像进行子带分解，并对每一子带与参与融合的多光谱图像和全色图像分别进行比对，计算以相关系数为基础的指标，最后进行指标综合，获得最终评价指标，为达此目的，本发明提供一种无参考图像的遥感图像融合性能评价方法，包括下述步骤：

步骤1：参数初始化，设置低分辨率多光谱图像M的光学MTF值m、多光谱图像M与全色图像P之间的名义分辨率比值r；

步骤2：对F各波段进行线性组合，构造亮度分量I；

步骤3：图像频率成分提取，包括以下步骤：

步骤3.1：根据m和r构造频率域低通滤波器G⁽⁰⁾和相应的高通滤波器H⁽⁰⁾；

步骤3.2：利用G⁽⁰⁾分别对融合图像F的每一波段进行低通滤波，并使用盒采样方法向下r倍重采样至多光谱图像M尺寸，获得小尺寸的退化图像F⁽¹⁾；

步骤3.3：根据m构造频率域低通滤波器G⁽¹⁾，对F⁽¹⁾和M各波段分别进行低通滤波，获得各自低频成分F^(1)，l和M^l，以及相应的高频成分F^(1)，h和M^h；

步骤3.4：利用H⁽⁰⁾对I和P进行高通滤波，获取各自高频成分I^h和P^h，然后利用H⁽⁰⁾对I^h和P^h再进行一次频率成分分离，获得高频成分I^hh和P^hh，以及低频成分I^hl和P^hl；

步骤4：子指标计算：根据F^(1)，l和M^l计算光谱保持度Q_λ；根据F^(1)，h和M^h计算空间平滑度第一组分

根据I^hl和P^hl计算空间平滑度第二组分

根据I^hh和p^hh计算空间保持度Q_σ；

步骤5：根据各子指标构造综合评价指标。

作为本发明进一步改进，步骤1中采用的多光谱图像M的光学MTF值m据下式估算m＝m_tot/0.6366

其中m_tot为多光谱图像的MTF值，由数据提供商提供，或自行进行测量，m缺省值为0.5。

作为本发明进一步改进，步骤2中对F各波段进行线性组合，构造亮度分量I，线性组合方法不限，权重系数之和也无需为1.各波段权重可由F对P作线性回归取得，缺省值为1。

作为本发明进一步改进，步骤3对融合图像F进行低通滤波后，盒采样方法向下r倍重采样至多光谱图像M尺寸，获得小尺寸的退化图像F⁽¹⁾.低通滤波器G⁽⁰⁾构造方法如下：

其中σ₀为空间域高斯滤波器的尺度参数；

对融合图像的第i波段滤波为

对应的高通滤波器H⁽⁰⁾为G⁽⁰⁾取反后中心元素+1。

作为本发明进一步改进，步骤3退化后的融合图像F⁽¹⁾和M各波段分别进行低通滤波，获得各自低频成分F^(1)，l和M^l，以及相应的高频成分F^(1)，h和M^h，所采用的低通滤波器G⁽¹⁾构造方法如下：

其中σ₁为空间域高斯滤波器的尺度参数；

对融合退化图像的第i波段滤波为

对多光谱图像的第i波段滤波为

各波段高频成分直接由该波段减去低频成分获得，即

F^(1)，h＝F⁽¹⁾-F^(1)，l

M^h＝M-M^l。

作为本发明进一步改进，步骤3利用H⁽⁰⁾对I和P进行高通滤波，获取各自高频成分

和

然后利用H⁽⁰⁾对I^h和p^h再进行一次频率成分分离，获得高频成分

和

以及低频成分I^hl＝I^h-I^hh和p^hl＝p^h＝p^hh。

作为本发明进一步改进，步骤4计算F^(1)，l和M^l对应波段间的UIQI值，将所有波段UIQI值的均值作为融合图像的光谱保持度指标Q_λ如下：

其中n为波段数，对于给定图像A和B，UIQI计算如下：

其中，μ_A，μ_B为图像A和B的均值，σ_A，σ_B为图像A和B的标准差，σ_AB为图像A和B的协方差。

作为本发明进一步改进，步骤4计算F^(1)，h和M^h对应波段间的相关系数C^MH，将所有波段相关系数的均值作为融合图像的空间平滑度指标第一组分

如下：

其中n为波段数。

作为本发明进一步改进，步骤4计算I^hl和P^hl之间的相关系数C^PL，作为空间平滑度指标第二组分

计算I^hh和P^hh之间的相关系数C^PH，作为空间保持度指标Q_σ。

作为本发明进一步改进，根据各子指标构造综合评价指标Q如下：

Q＝Q_λQ_sQ_σ

其中，Q_s定义为

本发明的方法与现有方法相比，改进主要体现在：1、适用性更广，无需参考图像；2、采取的理论假设较少，避免了在现有无参考图像和降分辨率方案中存在的问题，对融合图像性能评价结果更稳定可靠，评价结果与专家经验一致；3、可以很好地表现融合算法的不同性能方面，可以帮助开发人员有针对性地对融合算法进行改进；4、运算速度快，适用于遥感图像融合评价。

具体实施方式

下面结合具体实施方式对本发明作进一步详细描述：

本发明提供一种无参考图像的遥感图像融合性能评价方法，本发明将图像融合性能分解为光谱保持度、空间平滑度和细节保持度三个维度，通过对融合图像进行子带分解，并对每一子带与参与融合的多光谱图像和全色图像分别进行比对，计算以相关系数为基础的指标，最后进行指标综合，获得最终评价指标。

作为本发明具体实施例，本发明采用真实WorldView 2星载遥感多光谱和全色光图像，分别采用IHS,Brovey,PCA,GS,GSA,HPF,ATWT等方法进行融合，并采用本发明方法对上述融合结果，以及通过立方卷积上采样至全色尺寸的多光谱图像进行评价。实施本发明包括以下步骤：

步骤1：参数初始化。依据卫星数据，设置低分辨率多光谱图像imgM的光学MTF值m＝0.6、多光谱图像imgM与全色图像imgP之间的名义分辨率比值r＝4.

步骤2：为待测试图像imgF计算波段间均值，构造亮度分量imgI.MATLAB代码示例如下：

imgI＝mean(imgF,3)；

步骤3：构造低通滤波器G0,G1和相应的高通滤波器H0.MATLAB代码示例如下：

N＝41；

alpha＝sqrt((N*(0.5/r))^2/(-2*log(m)))；

H＝fspecial('gaussian',N,alpha)；

Hd＝H./max(H(:))；

h＝fwind1(Hd,kaiser(N))；

G0＝real(h)；

H0＝padarray(1,[(N-1)/2(N-1)/2],0,'both')-G0；

alpha＝sqrt((N*0.5)^2/(-2*log(m)))；

H＝fspecial('gaussian',N,alpha)；

Hd＝H./max(H(:))；

h＝fwind1(Hd,kaiser(N))；

G1＝real(h)；

步骤4：利用G0分别对待测试图像imgF的每一波段进行低通滤波，并使用盒采样方法向下r倍重采样至多光谱图像尺寸.MATLAB示例代码如下：

imgF_L＝double(imfilter(imgF,G1,'circular'))；

imgF_L＝imresize(imgF_L,1/r,'box')；

步骤5：利用G1分别对imgF_L和imgM各波段进行滤波，提取各波段的低频和高频成分。MATLAB示例代码如下：

imgF_LL＝double(imfilter(imgF_L,G1,'circular'))；

imgF_LH＝imgF_L-imgF_LL；

imgM_L＝double(imfilter(imgM,G1,'circular'))；

imgM_H＝imgM-imgM_L；

步骤6：计算图像的光谱保持度指标Q_Lamb和空间平滑度第一组分Q_Sig1.MATLAB示例代码如下：

for ii＝1:nbands

Q_ML(ii)＝img_qi(imgF_LL(:,:,ii),imgM_LL(:,:,ii))；

Q_MH(ii)＝corr2(imgF_LH(:,:,ii),imgM_LH(:,:,ii))；

end

Q_Lamb＝mean(Q_ML)；

Q_Sig1＝mean(Q_MH)；

步骤7：利用G0和H0分别对imgI和imgP进行滤波，提取各波段的低频和高频成分。MATLAB示例代码如下：

imgI_H＝double(imfilter(imgI,H0,'circular'))；

imgI_HL＝double(imfilter(imgI_H,G0,'circular'))；

imgI_HH＝imgI_H-imgI_HL；

imgP_H＝double(imfilter(imgP,H0,'circular'))；

imgP_HL＝double(imfilter(imgP_H,G0,'circular'))；

imgP_HH＝imgP_H-imgP_HL；

步骤8：计算图像的空间平滑度第二组分Q_Sig2和空间保持度指标Q_Spa和.MATLAB示例代码如下：

Q_Sig2＝corr2(imgI_HL,imgP_HL)；

Q_Spa＝corr2(imgI_HH,imgP_HH)；

步骤9：计算并输出评价结果Q。MATLAB示例代码如下：

Q_Sig＝sqrt(Q_Sig1*Q_Sig2)；

Q＝Q_Lamb*Q_Sig*Q_Spa；

融合图像评价结果如下：

上采样MS在光谱保持方面较好，但空间平滑度不够，细节方面则完全没有得到增强，因此评分很低，仅为0.155；IHS和Brovey同为基于色彩空间的分量替换方法，得分非常接近，符合理论预期；随后发展出的PCA方法主要改善了边缘处的空间平滑度，GS方法则进一步提高了PCA的光谱保持能力(已有理论证明，PCA为GS方法的特例)；GSA方法在亮度分量构造方面对GS方法进行了改进，因此光谱保持能力大幅度提高；HPF为代表的多尺度分析方法天然具有较好的光谱保持能力，因此光谱保持度高，但空间保持度下降，表现为融合图像不够锐利；ATWT则通过滤波器设计，提高了空间细节保持能力，总体得分最高。从评价结果可以看出，三个子指标可以很好的反映图像在这三个方面的表现，结果与已有研究和主观评价结论一致。

以上所述，仅是本发明的较佳实施例而已，并非是对本发明作任何其他形式的限制，而依据本发明的技术实质所作的任何修改或等同变化，仍属于本发明所要求保护的范围。

Claims (10)

1.一种无参考图像的遥感图像融合性能评价方法，包括下述步骤，其特征在于：

步骤1：参数初始化，设置低分辨率多光谱图像M的光学MTF值m、多光谱图像M与全色图像P之间的名义分辨率比值r；

步骤2：对F各波段进行线性组合，构造亮度分量I；

步骤3：图像频率成分提取，包括以下步骤：

步骤3.1：根据m和r构造频率域低通滤波器G⁽⁰⁾和相应的高通滤波器H⁽⁰⁾；

步骤3.2：利用G⁽⁰⁾分别对融合图像F的每一波段进行低通滤波，并使用盒采样方法向下r倍重采样至多光谱图像M尺寸，获得小尺寸的退化图像F⁽¹⁾；

步骤3.3：根据m构造频率域低通滤波器G⁽¹⁾，对F⁽¹⁾和M各波段分别进行低通滤波，获得各自低频成分F^(1)，l和M^l，以及相应的高频成分F^(1)，h和M^h；

步骤3.4：利用H⁽⁰⁾对I和P进行高通滤波，获取各自高频成分I^h和P^h，然后利用H⁽⁰⁾对I^h和P^h再进行一次频率成分分离，获得高频成分I^hh和P^hh，以及低频成分I^hl和P^hl；

步骤4：子指标计算：根据F^(1)，l和M^l计算光谱保持度Q_λ；根据F^(1)，h和M^h计算空间平滑度第一组分

根据I^hl和P^hl计算空间平滑度第二组分

根据I^hh和P^hh计算空间保持度Q_σ；

步骤5：根据各子指标构造综合评价指标。

2.根据权利要求1所述的一种无参考图像的遥感图像融合性能评价方法，其特征在于：步骤1中采用的多光谱图像M的光学MTF值m据下式估算

m＝m_tot/0.6366

其中m_tot为多光谱图像的MTF值，由数据提供商提供，或自行进行测量，m缺省值为0.5。

3.根据权利要求1所述的一种无参考图像的遥感图像融合性能评价方法，其特征在于：步骤2中对F各波段进行线性组合，构造亮度分量I，线性组合方法不限，权重系数之和也无需为1.各波段权重可由F对P作线性回归取得，缺省值为1。

4.根据权利要求1所述的一种无参考图像的遥感图像融合性能评价方法，其特征在于：步骤3对融合图像F进行低通滤波后，盒采样方法向下r倍重采样至多光谱图像M尺寸，获得小尺寸的退化图像F⁽¹⁾.低通滤波器G⁽⁰⁾构造方法如下：

其中σ₀为空间域高斯滤波器的尺度参数；

对融合图像的第i波段滤波为

对应的高通滤波器H⁽⁰⁾为G⁽⁰⁾取反后中心元素+1。

5.根据权利要求1所述的一种无参考图像的遥感图像融合性能评价方法，其特征在于：步骤3退化后的融合图像F⁽¹⁾和M各波段分别进行低通滤波，获得各自低频成分F^(1)，l和M^l，以及相应的高频成分F^(1)，h和M^h，所采用的低通滤波器G⁽¹⁾构造方法如下：

其中σ₁为空间域高斯滤波器的尺度参数；

对融合退化图像的第i波段滤波为

对多光谱图像的第i波段滤波为

各波段高频成分直接由该波段减去低频成分获得，即

F^(1)，h＝F⁽¹⁾-F^(1)，l

M^h＝M-M^l。

6.根据权利要求1所述的一种无参考图像的遥感图像融合性能评价方法，其特征在于：步骤3利用H⁽⁰⁾对I和P进行高通滤波，获取各自高频成分

和

然后利用H⁽⁰⁾对I^h和P^h再进行一次频率成分分离，获得高频成分

和

以及低频成分I^hl＝I^h-I^hh和P^hl＝P^h-P^hh。

7.根据权利要求1所述的一种无参考图像的遥感图像融合性能评价方法，其特征在于：步骤4计算F^(1)，l和M^l对应波段间的UIQI值，将所有波段UIQI值的均值作为融合图像的光谱保持度指标Q_λ如下：

其中n为波段数，对于给定图像A和B，UIQI计算如下：

其中，μ_A，μ_B为图像A和B的均值，σ_A，σ_B为图像A和B的标准差，σ_AB为图像A和B的协方差。

8.根据权利要求1所述的一种无参考图像的遥感图像融合性能评价方法，其特征在于：步骤4计算F^(1)，h和Mh对应波段间的相关系数C^MH，将所有波段相关系数的均值作为融合图像的空间平滑度指标第一组分

如下：

其中n为波段数。

9.根据权利要求1所述的一种无参考图像的遥感图像融合性能评价方法，其特征在于：步骤4计算I^hl和P^hl之间的相关系数C^PL，作为空间平滑度指标第二组分

计算I^hh和P^hh之间的相关系数C^PH，作为空间保持度指标Q_σ

10.根据权利要求1和6所述的一种无参考图像的遥感图像融合性能评价方法，其特征在于：根据各子指标构造综合评价指标Q如下：

Q＝Q_λQ_sQ_σ

其中，Q_s定义为