CN109191450A - 一种遥感图像融合质量评价方法 - Google Patents

Abstract

本发明公开了一种遥感图像融合质量评价方法，其考虑到光谱域失真和空间域失真对遥感融合影像质量的影响，提取了红波段图像、绿波段图像、蓝波段图像、近红外波段图像、归一化植被指数图像、归一化差分水体指数图像、道路和建筑物指数图像，并通过提取所有子块的统计特征矢量，构造出训练样本集合的原始多元高斯模型；在测试阶段，通过获取测试遥感融合影像构成的测试样本集合的原始多元高斯模型，并根据训练阶段构造的原始多元高斯模型，预测得到测试遥感融合影像的客观质量评价预测值，由于获得的统计特征矢量信息能够较好地反映遥感融合影像的质量变化情况，因此有效地提高了客观评价结果与主观感知之间的相关性。

Description

一种遥感图像融合质量评价方法

技术领域

本发明涉及一种图像质量评价方法，尤其是涉及一种遥感图像融合质量评价方法。

背景技术

全色多光谱融合是目前遥感图像融合领域的一个研究热点。多光谱图像通常空间分辨率较低，而全色图像通常具有较高的空间分辨率，但缺乏光谱信息。全色多光谱融合的目的就是将全色图像的空间信息融合到多光谱图像中，得到具有高空间分辨率的多光谱图像。但是，不同融合方法对不同的全色图像和多光谱图像有着不同的融合效果，如何评价融合后图像质量问题是目前需要解决的一个问题。

目前遥感全色多光谱融合图像质量评价领域中，主要存在两方面问题：(1)缺少全色多光谱融合图像的主观评价库，现有的客观质量评价方法无法通过人类主观评价的对比来衡量效果。(2)光谱域失真问题难以用无参考评价方法来衡量，而光谱域失真过于严重可能导致地物识别困难，因此度量光谱域失真是评价融合图像质量必须要解决的问题。

发明内容

本发明所要解决的技术问题是提供一种遥感图像融合质量评价方法，其能够有效地提高客观评价结果与主观感知之间的相关性。

本发明解决上述技术问题所采用的技术方案为：一种遥感图像融合质量评价方法，其特征在于包括训练阶段和测试阶段两个过程；

所述的训练阶段过程的具体步骤如下：

①_1、选取N幅宽度为W且高度为H的原始的无失真多光谱影像，每幅原始的无失真多光谱影像包括红波段图像、绿波段图像、蓝波段图像、近红外波段图像，将第u幅原始的无失真多光谱影像的红波段图像、绿波段图像、蓝波段图像、近红外波段图像对应记为R_u、G_u、B_u、NIR_u；其中，N＞1，u为正整数，1≤u≤N；

①_2、计算每幅原始的无失真多光谱影像的归一化植被指数图像，将第u幅原始的无失真多光谱影像的归一化植被指数图像记为NDVI_u，将NDVI_u中坐标位置为(x,y)的像素点的像素值记为NDVI_u(x,y)，其中，1≤x≤W,1≤y≤H，NIR_u(x,y)表示NIR_u中坐标位置为(x,y)的像素点的像素值，R_u(x,y)表示R_u中坐标位置为(x,y)的像素点的像素值；

计算每幅原始的无失真多光谱影像的归一化差分水体指数图像，将第u幅原始的无失真多光谱影像的归一化差分水体指数图像记为NDWI_u，将NDWI_u中坐标位置为(x,y)的像素点的像素值记为NDWI_u(x,y)，其中，G_u(x,y)表示G_u中坐标位置为(x,y)的像素点的像素值；

计算每幅原始的无失真多光谱影像的道路和建筑物指数图像，将第u幅原始的无失真多光谱影像的道路和建筑物指数图像记为RatioG_u，将RatioG_u中坐标位置为(x,y)的像素点的像素值记为RatioG_u(x,y)，其中，B_u(x,y)表示B_u中坐标位置为(x,y)的像素点的像素值；

①_3、从每幅原始的无失真多光谱影像的红波段图像、绿波段图像、蓝波段图像、近红外波段图像、归一化植被指数图像、归一化差分水体指数图像、道路和建筑物指数图像各自中同坐标位置截取一个尺寸大小为96×96的子块，并对截取的每个子块进行标记以指示其来源，以随机截取方式，重复截取的过程，共截取M次后从每幅原始的无失真多光谱影像的红波段图像、绿波段图像、蓝波段图像、近红外波段图像、归一化植被指数图像、归一化差分水体指数图像、道路和建筑物指数图像各自中截取了M个子块；然后将截取的每个子块中的所有像素点的像素值按序排列构成该子块对应的列向量；再将每种标记的M'个子块对应的列向量构成一个训练样本集合，将标记为γ的M'个子块对应的列向量构成的训练样本集合记为其中，M≥7，M'表示标记相同的子块的总个数，M'＝N×M，γ∈{R,G,B,NIR,NDVI,NDWI,RatioG}，γ＝R表示该子块截取自红波段图像，γ＝G表示该子块截取自绿波段图像，γ＝B表示该子块截取自蓝波段图像，γ＝NIR表示该子块截取自近红外波段图像，γ＝NDVI表示该子块截取自归一化植被指数图像，γ＝NDWI表示该子块截取自归一化差分水体指数图像，γ＝RatioG表示该子块截取自道路和建筑物指数图像，m为正整数，1≤m≤M'，表示标记为γ的第m个子块对应的列向量，的维数为(96×96)×1；

①_4、计算7个训练样本集合中同索引号的列向量对应的子块的统计特征矢量，共有M'个统计特征矢量，将7个训练样本集合中索引号同为m的列向量对应的子块的统计特征矢量记为F_m；其中，F_m的维数为42×1；

①_5、根据7个训练样本集合对应的M'个统计特征矢量，提取出7个训练样本集合的原始多元高斯模型，记为(μ,C)；其中，μ表示(μ,C)的均值向量，C表示(μ,C)的协方差矩阵；

所述的测试阶段过程的具体步骤如下：

②_1、对于任意一幅宽度为W'且高度为H'的测试遥感融合影像，其包括红波段图像、绿波段图像、蓝波段图像、近红外波段图像，对应记为R_test,G_test,B_test,NIR_test；

②_2、按照步骤①_2的过程，以相同的操作方式，得到测试遥感融合影像的归一化植被指数图像、归一化差分水体指数图像、道路和建筑物指数图像，对应记为NDVI_test,NDWI_test,RatioG_test；

②_3、采用窗口尺寸大小为96×96且窗口移动步长为1像素的移动窗口，在R_test,G_test,B_test,NIR_test和NDVI_test,NDWI_test,RatioG_test各自中移动以截取同坐标位置的尺寸大小为96×96的子块，并对从R_test,G_test,B_test,NIR_test和NDVI_test,NDWI_test,RatioG_test中截取的每个子块进行标记以指示其来源，共移动M”次后从R_test,G_test,B_test,NIR_test和NDVI_test,NDWI_test,RatioG_test各自中截取了M”个子块；然后将从R_test,G_test,B_test,NIR_test和NDVI_test,NDWI_test,RatioG_test中截取的每个子块中的所有像素点的像素值按序排列构成该子块对应的列向量；再将每种标记的M”个子块对应的列向量构成一个测试样本集合，将标记为γ的M”个子块对应的列向量构成的测试样本集合记为其中，M”表示R_test,G_test,B_test,NIR_test和NDVI_test,NDWI_test,RatioG_test中标记相同的子块的总个数，1＜M”≤(W'-96)×(H'-96)，m”为正整数，1≤m”≤M”，表示标记为γ的第m”个子块对应的列向量，的维数为(96×96)×1；

②_4、按照步骤①_4的过程，以相同的方式计算7个测试样本集合中同索引号的列向量对应的子块的统计特征矢量，共有M”个统计特征矢量，将7个测试样本集合中索引号同为m”的列向量对应的子块的统计特征矢量记为F_test,m”；其中，F_test,m”的维数为42×1；

②_5、根据7个测试样本集合对应的M”个统计特征矢量，按照步骤①_5的过程，以相同的方式提取出7个测试样本集合的原始多元高斯模型，记为(μ_test,C_test)；其中，μ_test表示(μ_test,C_test)的均值向量，C_test表示(μ_test,C_test)的协方差矩阵；

②_6、根据(μ,C)和(μ_test,C_test)，计算测试遥感融合影像的客观质量评价预测值，记为Q_test，其中，(μ-μ_test)^T为(μ-μ_test)的转置，为的逆。

所述的步骤①_4中的F_m的获取过程为：

①_4a、对于γ＝R的子块，对中的所有元素的值进行高斯和分布拟合，拟合得到中的所有元素的值的高斯和分布拟合曲线，记为 其中，表示输入变量，表示通过拟合得到的中的6个高斯和分布参数，exp()表示以自然基数e为底的指数函数；

对于γ＝G的子块，对中的所有元素的值进行高斯和分布拟合，拟合得到中的所有元素的值的高斯和分布拟合曲线，记为 其中，表示通过拟合得到的中的6个高斯和分布参数；

对于γ＝B的子块，对中的所有元素的值进行高斯和分布拟合，拟合得到中的所有元素的值的高斯和分布拟合曲线，记为 其中，表示通过拟合得到的中的6个高斯和分布参数；

对于γ＝NIR的子块，对中的所有元素的值进行高斯和分布拟合，拟合得到中的所有元素的值的高斯和分布拟合曲线，记为 其中，表示通过拟合得到的中的6个高斯和分布参数；

对于γ＝NDVI的子块，对中的所有元素的值进行高斯和分布拟合，拟合得到中的所有元素的值的高斯和分布拟合曲线，记为 其中，表示通过拟合得到的中的6个高斯和分布参数；

对于γ＝NDWI的子块，对中的所有元素的值进行高斯和分布拟合，拟合得到中的所有元素的值的高斯和分布拟合曲线，记为 其中，表示通过拟合得到的中的6个高斯和分布参数；

对于γ＝RatioG的子块，对中的所有元素的值进行高斯和分布拟合，拟合得到中的所有元素的值的高斯和分布拟合曲线，记为 其中，表示通过拟合得到的中的6个高斯和分布参数；

①_4b、将 和按序排列构成F_m。

所述的步骤①_5中的(μ,C)的获取过程为：

①_5a、根据7个训练样本集合对应的M'个统计特征矢量，获取截取的7M'个子块的特征矩阵，记为X，其中，X的维数为M'×42，(F₁)^T为F₁的转置，F₁表示7个训练样本集合中索引号同为1的列向量对应的子块的统计特征矢量，(F₂)^T为F₂的转置，F₂表示7个训练样本集合中索引号同为2的列向量对应的子块的统计特征矢量，(F_m)^T为F_m的转置，(F_M')^T为F_M'的转置，F_M'表示7个训练样本集合中索引号同为M'的列向量对应的子块的统计特征矢量；

①_5b、计算X中的列向量的均值向量，记为μ，μ＝[μ₁ μ₂ … μ₄₂]，其中，μ的维数为1×42，X(m,1)表示X中第m行第1列的元素的值，X(m,2)表示X中第m行第2列的元素的值，X(m,42)表示X中第m行第42列的元素的值；

计算X中的列向量的协方差矩阵，记为C，σ_1,1＝cov(X(m,1),X(m,1))，σ_1,2＝cov(X(m,1),X(m,2))，σ_1,42＝cov(X(m,1),X(m,42))，σ_2,1＝cov(X(m,2),X(m,1))，σ_2,2＝cov(X(m,2),X(m,2))，σ_2,42＝cov(X(m,2),X(m,42))，σ_42,1＝cov(X(m,42),X(m,1))，σ_42,2＝cov(X(m,42),X(m,2))，σ_42,42＝cov(X(m,42),X(m,42))；其中，C的维数为42×42，cov(a,b)表示计算向量a和b的协方差操作；

①_5c、根据μ和C，构造得到(μ,C)。

与现有技术相比，本发明的优点在于：

本发明方法考虑到光谱域失真和空间域失真对遥感融合影像质量的影响，提取了红波段图像、绿波段图像、蓝波段图像、近红外波段图像、归一化植被指数图像、归一化差分水体指数图像、道路和建筑物指数图像，并通过提取所有子块的统计特征矢量，构造出训练样本集合的原始多元高斯模型；在测试阶段，通过获取测试遥感融合影像构成的测试样本集合的原始多元高斯模型，并根据训练阶段构造的原始多元高斯模型，预测得到测试遥感融合影像的客观质量评价预测值，由于获得的统计特征矢量信息能够较好地反映遥感融合影像的质量变化情况，因此有效地提高了客观评价结果与主观感知之间的相关性。

附图说明

图1为本发明方法的总体实现框图。

具体实施方式

以下结合附图实施例对本发明作进一步详细描述。

本发明提出的一种遥感图像融合质量评价方法，其总体实现框图如图1所示，其包括训练阶段和测试阶段两个过程；

所述的训练阶段过程的具体步骤如下：

①_1、选取N幅宽度为W且高度为H的原始的无失真多光谱影像，每幅原始的无失真多光谱影像包括红波段图像、绿波段图像、蓝波段图像、近红外波段图像，将第u幅原始的无失真多光谱影像的红波段图像、绿波段图像、蓝波段图像、近红外波段图像对应记为R_u、G_u、B_u、NIR_u；其中，N＞1，在本实施里中取N＝100，u为正整数，1≤u≤N。

①_2、计算每幅原始的无失真多光谱影像的归一化植被指数图像，将第u幅原始的无失真多光谱影像的归一化植被指数图像记为NDVI_u，将NDVI_u中坐标位置为(x,y)的像素点的像素值记为NDVI_u(x,y)，其中，1≤x≤W,1≤y≤H，NIR_u(x,y)表示NIR_u中坐标位置为(x,y)的像素点的像素值，R_u(x,y)表示R_u中坐标位置为(x,y)的像素点的像素值。

计算每幅原始的无失真多光谱影像的归一化差分水体指数图像，将第u幅原始的无失真多光谱影像的归一化差分水体指数图像记为NDWI_u，将NDWI_u中坐标位置为(x,y)的像素点的像素值记为NDWI_u(x,y)，其中，G_u(x,y)表示G_u中坐标位置为(x,y)的像素点的像素值。

计算每幅原始的无失真多光谱影像的道路和建筑物指数图像，将第u幅原始的无失真多光谱影像的道路和建筑物指数图像记为RatioG_u，将RatioG_u中坐标位置为(x,y)的像素点的像素值记为RatioG_u(x,y)，其中，B_u(x,y)表示B_u中坐标位置为(x,y)的像素点的像素值。

①_3、从每幅原始的无失真多光谱影像的红波段图像、绿波段图像、蓝波段图像、近红外波段图像、归一化植被指数图像、归一化差分水体指数图像、道路和建筑物指数图像各自中同坐标位置截取一个尺寸大小为96×96的子块，并对截取的每个子块进行标记以指示其来源，以随机截取方式，重复截取的过程，共截取M次后从每幅原始的无失真多光谱影像的红波段图像、绿波段图像、蓝波段图像、近红外波段图像、归一化植被指数图像、归一化差分水体指数图像、道路和建筑物指数图像各自中截取了M个子块，即针对N幅红波段图像共有N×M个子块、针对N幅绿波段图像共有N×M个子块、针对N幅蓝波段图像共有N×M个子块、针对N幅近红外波段图像共有N×M个子块、针对N幅归一化植被指数图像共有N×M个子块、针对N幅归一化差分水体指数图像共有N×M个子块、针对N幅道路和建筑物指数图像共有N×M个子块；然后将截取的每个子块中的所有像素点的像素值按序排列构成该子块对应的列向量；再将每种标记的M'个子块对应的列向量构成一个训练样本集合，将标记为γ的M'个子块对应的列向量构成的训练样本集合记为其中，M≥7，在本实施例中取M＝100，M'表示标记相同的子块的总个数，M'＝N×M，γ∈{R,G,B,NIR,NDVI,NDWI,RatioG}，γ＝R表示该子块截取自红波段图像，γ＝G表示该子块截取自绿波段图像，γ＝B表示该子块截取自蓝波段图像，γ＝NIR表示该子块截取自近红外波段图像，γ＝NDVI表示该子块截取自归一化植被指数图像，γ＝NDWI表示该子块截取自归一化差分水体指数图像，γ＝RatioG表示该子块截取自道路和建筑物指数图像，m为正整数，1≤m≤M'，表示标记为γ的第m个子块对应的列向量，的维数为(96×96)×1。

①_4、计算7个训练样本集合中同索引号的列向量对应的子块的统计特征矢量，共有M'个统计特征矢量，将7个训练样本集合中索引号同为m的列向量对应的子块的统计特征矢量记为F_m，将F_m作为原始的无失真多光谱影像中与7个训练样本集合中索引号同为m的列向量对应的子块同坐标位置的尺寸大小为96×96的区域的统计特征矢量；其中，F_m的维数为42×1。

在本实施例中，步骤①_4中的F_m的获取过程为：

①_4a、对于γ＝R的子块，对中的所有元素的值进行高斯和(Sum-of-Gaussians)分布拟合，拟合得到中的所有元素的值的高斯和分布拟合曲线，记为 其中，表示输入变量，表示通过拟合得到的中的6个高斯和分布参数，exp()表示以自然基数e为底的指数函数。

对于γ＝G的子块，对中的所有元素的值进行高斯和(Sum-of-Gaussians)分布拟合，拟合得到中的所有元素的值的高斯和分布拟合曲线，记为 其中，表示通过拟合得到的中的6个高斯和分布参数。

对于γ＝B的子块，对中的所有元素的值进行高斯和(Sum-of-Gaussians)分布拟合，拟合得到中的所有元素的值的高斯和分布拟合曲线，记为 其中，表示通过拟合得到的中的6个高斯和分布参数。

对于γ＝NIR的子块，对中的所有元素的值进行高斯和(Sum-of-Gaussians)分布拟合，拟合得到中的所有元素的值的高斯和分布拟合曲线，记为 其中，表示通过拟合得到的中的6个高斯和分布参数。

对于γ＝NDVI的子块，对中的所有元素的值进行高斯和(Sum-of-Gaussians)分布拟合，拟合得到中的所有元素的值的高斯和分布拟合曲线，记为 其中，表示通过拟合得到的中的6个高斯和分布参数。

对于γ＝NDWI的子块，对中的所有元素的值进行高斯和(Sum-of-Gaussians)分布拟合，拟合得到中的所有元素的值的高斯和分布拟合曲线，记为 其中，表示通过拟合得到的中的6个高斯和分布参数。

对于γ＝RatioG的子块，对中的所有元素的值进行高斯和(Sum-of-Gaussians)分布拟合，拟合得到中的所有元素的值的高斯和分布拟合曲线，记为 其中，表示通过拟合得到的中的6个高斯和分布参数。

①_4b、将 和按序排列构成F_m。

①_5、根据7个训练样本集合对应的M'个统计特征矢量，提取出7个训练样本集合的原始多元高斯(Pristine Multivariate Gaussian，MVG)模型，记为(μ,C)；其中，μ表示(μ,C)的均值向量，C表示(μ,C)的协方差矩阵。

在本实施例中，步骤①_5中的(μ,C)的获取过程为：

①_5a、根据7个训练样本集合对应的M'个统计特征矢量，获取截取的7M'个子块的特征矩阵，记为X，其中，X的维数为M'×42，(F₁)^T为F₁的转置，F₁表示7个训练样本集合中索引号同为1的列向量对应的子块的统计特征矢量，(F₂)^T为F₂的转置，F₂表示7个训练样本集合中索引号同为2的列向量对应的子块的统计特征矢量，(F_m)^T为F_m的转置，(F_M')^T为F_M'的转置，F_M'表示7个训练样本集合中索引号同为M'的列向量对应的子块的统计特征矢量。

①_5b、计算X中的列向量的均值向量，记为μ，μ＝[μ₁ μ₂ … μ₄₂]，其中，μ的维数为1×42，X(m,1)表示X中第m行第1列的元素的值，X(m,2)表示X中第m行第2列的元素的值，X(m,42)表示X中第m行第42列的元素的值。

计算X中的列向量的协方差矩阵，记为C，σ_1,1＝cov(X(m,1),X(m,1))，σ_1,2＝cov(X(m,1),X(m,2))，σ_1,42＝cov(X(m,1),X(m,42))，σ_2,1＝cov(X(m,2),X(m,1))，σ_2,2＝cov(X(m,2),X(m,2))，σ_2,42＝cov(X(m,2),X(m,42))，σ_42,1＝cov(X(m,42),X(m,1))，σ_42,2＝cov(X(m,42),X(m,2))，σ_42,42＝cov(X(m,42),X(m,42))；其中，C的维数为42×42，cov(a,b)表示计算向量a和b的协方差操作。

①_5c、根据μ和C，构造得到(μ,C)。

所述的测试阶段过程的具体步骤如下：

②_1、对于任意一幅宽度为W'且高度为H'的测试遥感融合影像，其包括红波段图像、绿波段图像、蓝波段图像、近红外波段图像，对应记为R_test,G_test,B_test,NIR_test。

②_2、按照步骤①_2的过程，以相同的操作方式，得到测试遥感融合影像的归一化植被指数图像、归一化差分水体指数图像、道路和建筑物指数图像，对应记为NDVI_test,NDWI_test,RatioG_test。

②_3、采用窗口尺寸大小为96×96且窗口移动步长为1像素的移动窗口，在R_test,G_test,B_test,NIR_test和NDVI_test,NDWI_test,RatioG_test各自中移动以截取同坐标位置的尺寸大小为96×96的子块，并对从R_test,G_test,B_test,NIR_test和NDVI_test,NDWI_test,RatioG_test中截取的每个子块进行标记以指示其来源，共移动M”次后从R_test,G_test,B_test,NIR_test和NDVI_test,NDWI_test,RatioG_test各自中截取了M”个子块；然后将从R_test,G_test,B_test,NIR_test和NDVI_test,NDWI_test,RatioG_test中截取的每个子块中的所有像素点的像素值按序排列构成该子块对应的列向量；再将每种标记的M”个子块对应的列向量构成一个测试样本集合，将标记为γ的M”个子块对应的列向量构成的测试样本集合记为其中，M”表示R_test,G_test,B_test,NIR_test和NDVI_test,NDWI_test,RatioG_test中标记相同的子块的总个数，1＜M”≤(W'-96)×(H'-96)，m”为正整数，1≤m”≤M”，表示标记为γ的第m”个子块对应的列向量，的维数为(96×96)×1。

②_4、按照步骤①_4的过程，以相同的方式计算7个测试样本集合中同索引号的列向量对应的子块的统计特征矢量，共有M”个统计特征矢量，将7个测试样本集合中索引号同为m”的列向量对应的子块的统计特征矢量记为F_test,m”，再将F_test,m”作为测试遥感融合影像中与7个测试样本集合中索引号同为m”的列向量对应的子块同坐标位置的尺寸大小为96×96的区域的统计特征矢量；其中，F_test,m”的维数为42×1。

②_5、根据7个测试样本集合对应的M”个统计特征矢量，按照步骤①_5的过程，以相同的方式提取出7个测试样本集合的原始多元高斯(Pristine Multivariate Gaussian，MVG)模型，记为(μ_test,C_test)；其中，μ_test表示(μ_test,C_test)的均值向量，C_test表示(μ_test,C_test)的协方差矩阵。

②_6、根据(μ,C)和(μ_test,C_test)，计算测试遥感融合影像的客观质量评价预测值，记为Q_test，其中，(μ-μ_test)^T为(μ-μ_test)的转置，为的逆。

为了进一步说明本发明方法的可行性和有效性，对本发明方法进行试验。

在本实施例中，采用本发明方法对宁波大学建立的遥感融合影像数据库中的遥感融合影像进行测试。该遥感融合影像数据库收集了20个由Quickbird卫星拍摄的全色空间分辨率为0.6米的Quickbird影像，然后对多光谱影像和全色影像采用18种融合方法进行融合，得到对应的遥感融合影像，这样该遥感融合影像数据库共有360幅遥感融合影像，并给出了每幅遥感融合影像的平均主观评分均值。在本实施例中，利用评估图像质量评价方法的3个常用客观参量作为评价指标，即非线性回归条件下的Pearson相关系数(Pearsonlinear correlation coefficient，PLCC)、Spearman相关系数(Spearman rank ordercorrelation coefficient，SROCC)、均方误差(root mean squared error，RMSE)，PLCC和RMSE反映客观质量评价预测值的准确性，SROCC反映其单调性。表1给出了采用本发明方法得到的客观质量评价预测值与平均主观评分均值之间的相关性，从表1中可以看出，采用本发明方法得到的客观质量评价预测值与平均主观评分均值之间的相关性是很高的，表明客观评价结果与人眼主观感知的结果较为一致，足以说明本发明方法的有效性。

表1采用本发明方法得到的客观质量评价预测值与平均主观评分均值之间的相关性

方法	PLCC	SROCC	RMSE
				本发明方法	0.7410	0.6913	8.8399

Claims (3)

1.一种遥感图像融合质量评价方法，其特征在于包括训练阶段和测试阶段两个过程；

所述的训练阶段过程的具体步骤如下：

①_1、选取N幅宽度为W且高度为H的原始的无失真多光谱影像，每幅原始的无失真多光谱影像包括红波段图像、绿波段图像、蓝波段图像、近红外波段图像，将第u幅原始的无失真多光谱影像的红波段图像、绿波段图像、蓝波段图像、近红外波段图像对应记为R_u、G_u、B_u、NIR_u；其中，N＞1，u为正整数，1≤u≤N；

①_2、计算每幅原始的无失真多光谱影像的归一化植被指数图像，将第u幅原始的无失真多光谱影像的归一化植被指数图像记为NDVI_u，将NDVI_u中坐标位置为(x,y)的像素点的像素值记为NDVI_u(x,y)，其中，1≤x≤W,1≤y≤H，NIR_u(x,y)表示NIR_u中坐标位置为(x,y)的像素点的像素值，R_u(x,y)表示R_u中坐标位置为(x,y)的像素点的像素值；

计算每幅原始的无失真多光谱影像的归一化差分水体指数图像，将第u幅原始的无失真多光谱影像的归一化差分水体指数图像记为NDWI_u，将NDWI_u中坐标位置为(x,y)的像素点的像素值记为NDWI_u(x,y)，其中，G_u(x,y)表示G_u中坐标位置为(x,y)的像素点的像素值；

计算每幅原始的无失真多光谱影像的道路和建筑物指数图像，将第u幅原始的无失真多光谱影像的道路和建筑物指数图像记为RatioG_u，将RatioG_u中坐标位置为(x,y)的像素点的像素值记为RatioG_u(x,y)，其中，B_u(x,y)表示B_u中坐标位置为(x,y)的像素点的像素值；

①_3、从每幅原始的无失真多光谱影像的红波段图像、绿波段图像、蓝波段图像、近红外波段图像、归一化植被指数图像、归一化差分水体指数图像、道路和建筑物指数图像各自中同坐标位置截取一个尺寸大小为96×96的子块，并对截取的每个子块进行标记以指示其来源，以随机截取方式，重复截取的过程，共截取M次后从每幅原始的无失真多光谱影像的红波段图像、绿波段图像、蓝波段图像、近红外波段图像、归一化植被指数图像、归一化差分水体指数图像、道路和建筑物指数图像各自中截取了M个子块；然后将截取的每个子块中的所有像素点的像素值按序排列构成该子块对应的列向量；再将每种标记的M'个子块对应的列向量构成一个训练样本集合，将标记为γ的M'个子块对应的列向量构成的训练样本集合记为其中，M≥7，M'表示标记相同的子块的总个数，M'＝N×M，γ∈{R,G,B,NIR,NDVI,NDWI,RatioG}，γ＝R表示该子块截取自红波段图像，γ＝G表示该子块截取自绿波段图像，γ＝B表示该子块截取自蓝波段图像，γ＝NIR表示该子块截取自近红外波段图像，γ＝NDVI表示该子块截取自归一化植被指数图像，γ＝NDWI表示该子块截取自归一化差分水体指数图像，γ＝RatioG表示该子块截取自道路和建筑物指数图像，m为正整数，1≤m≤M'，表示标记为γ的第m个子块对应的列向量，的维数为(96×96)×1；

①_4、计算7个训练样本集合中同索引号的列向量对应的子块的统计特征矢量，共有M'个统计特征矢量，将7个训练样本集合中索引号同为m的列向量对应的子块的统计特征矢量记为F_m；其中，F_m的维数为42×1；

①_5、根据7个训练样本集合对应的M'个统计特征矢量，提取出7个训练样本集合的原始多元高斯模型，记为(μ,C)；其中，μ表示(μ,C)的均值向量，C表示(μ,C)的协方差矩阵；

所述的测试阶段过程的具体步骤如下：

②_1、对于任意一幅宽度为W'且高度为H'的测试遥感融合影像，其包括红波段图像、绿波段图像、蓝波段图像、近红外波段图像，对应记为R_test,G_test,B_test,NIR_test；

②_2、按照步骤①_2的过程，以相同的操作方式，得到测试遥感融合影像的归一化植被指数图像、归一化差分水体指数图像、道路和建筑物指数图像，对应记为NDVI_test,NDWI_test,RatioG_test；

②_3、采用窗口尺寸大小为96×96且窗口移动步长为1像素的移动窗口，在R_test,G_test,B_test,NIR_test和NDVI_test,NDWI_test,RatioG_test各自中移动以截取同坐标位置的尺寸大小为96×96的子块，并对从R_test,G_test,B_test,NIR_test和NDVI_test,NDWI_test,RatioG_test中截取的每个子块进行标记以指示其来源，共移动M”次后从R_test,G_test,B_test,NIR_test和NDVI_test,NDWI_test,RatioG_test各自中截取了M”个子块；然后将从R_test,G_test,B_test,NIR_test和NDVI_test,NDWI_test,RatioG_test中截取的每个子块中的所有像素点的像素值按序排列构成该子块对应的列向量；再将每种标记的M”个子块对应的列向量构成一个测试样本集合，将标记为γ的M”个子块对应的列向量构成的测试样本集合记为其中，M”表示R_test,G_test,B_test,NIR_test和NDVI_test,NDWI_test,RatioG_test中标记相同的子块的总个数，1＜M”≤(W'-96)×(H'-96)，m”为正整数，1≤m”≤M”，表示标记为γ的第m”个子块对应的列向量，的维数为(96×96)×1；

②_4、按照步骤①_4的过程，以相同的方式计算7个测试样本集合中同索引号的列向量对应的子块的统计特征矢量，共有M”个统计特征矢量，将7个测试样本集合中索引号同为m”的列向量对应的子块的统计特征矢量记为F_test,m”；其中，F_test,m”的维数为42×1；

②_5、根据7个测试样本集合对应的M”个统计特征矢量，按照步骤①_5的过程，以相同的方式提取出7个测试样本集合的原始多元高斯模型，记为(μ_test,C_test)；其中，μ_test表示(μ_test,C_test)的均值向量，C_test表示(μ_test,C_test)的协方差矩阵；

②_6、根据(μ,C)和(μ_test,C_test)，计算测试遥感融合影像的客观质量评价预测值，记为Q_test，其中，(μ-μ_test)^T为(μ-μ_test)的转置，为的逆。

2.根据权利要求1所述的一种遥感图像融合质量评价方法，其特征在于所述的步骤①_4中的F_m的获取过程为：

①_4a、对于γ＝R的子块，对中的所有元素的值进行高斯和分布拟合，拟合得到中的所有元素的值的高斯和分布拟合曲线，记为 其中，表示输入变量，表示通过拟合得到的中的6个高斯和分布参数，exp()表示以自然基数e为底的指数函数；

对于γ＝G的子块，对中的所有元素的值进行高斯和分布拟合，拟合得到中的所有元素的值的高斯和分布拟合曲线，记为 其中，表示通过拟合得到的中的6个高斯和分布参数；

对于γ＝B的子块，对中的所有元素的值进行高斯和分布拟合，拟合得到中的所有元素的值的高斯和分布拟合曲线，记为 其中，表示通过拟合得到的中的6个高斯和分布参数；

对于γ＝NIR的子块，对中的所有元素的值进行高斯和分布拟合，拟合得到中的所有元素的值的高斯和分布拟合曲线，记为 其中，表示通过拟合得到的中的6个高斯和分布参数；

对于γ＝NDVI的子块，对中的所有元素的值进行高斯和分布拟合，拟合得到中的所有元素的值的高斯和分布拟合曲线，记为 其中，表示通过拟合得到的中的6个高斯和分布参数；

对于γ＝NDWI的子块，对中的所有元素的值进行高斯和分布拟合，拟合得到中的所有元素的值的高斯和分布拟合曲线，记为 其中，表示通过拟合得到的中的6个高斯和分布参数；

对于γ＝RatioG的子块，对中的所有元素的值进行高斯和分布拟合，拟合得到中的所有元素的值的高斯和分布拟合曲线，记为 其中，表示通过拟合得到的中的6个高斯和分布参数；

①_4b、将 和按序排列构成F_m。

3.根据权利要求1或2所述的一种遥感图像融合质量评价方法，其特征在于所述的步骤①_5中的(μ,C)的获取过程为：

①_5a、根据7个训练样本集合对应的M'个统计特征矢量，获取截取的7M'个子块的特征矩阵，记为X，其中，X的维数为M'×42，(F₁)^T为F₁的转置，F₁表示7个训练样本集合中索引号同为1的列向量对应的子块的统计特征矢量，(F₂)^T为F₂的转置，F₂表示7个训练样本集合中索引号同为2的列向量对应的子块的统计特征矢量，(F_m)^T为F_m的转置，(F_M')^T为F_M'的转置，F_M'表示7个训练样本集合中索引号同为M'的列向量对应的子块的统计特征矢量；

①_5b、计算X中的列向量的均值向量，记为μ，μ＝[μ₁ μ₂ … μ₄₂]，其中，μ的维数为1×42，X(m,1)表示X中第m行第1列的元素的值，X(m,2)表示X中第m行第2列的元素的值，X(m,42)表示X中第m行第42列的元素的值；

计算X中的列向量的协方差矩阵，记为C，σ_1,1＝cov(X(m,1),X(m,1))，σ_1,2＝cov(X(m,1),X(m,2))，σ_1,42＝cov(X(m,1),X(m,42))，σ_2,1＝cov(X(m,2),X(m,1))，σ_2,2＝cov(X(m,2),X(m,2))，σ_2,42＝cov(X(m,2),X(m,42))，σ_42,1＝cov(X(m,42),X(m,1))，σ_42,2＝cov(X(m,42),X(m,2))，σ_42,42＝cov(X(m,42),X(m,42))；其中，C的维数为42×42，cov(a,b)表示计算向量a和b的协方差操作；

①_5c、根据μ和C，构造得到(μ,C)。