一种基于增强反投影的全色锐化后处理方法与流程

本发明属于遥感图像处理、图像融合以及全色锐化技术领域，特别涉及一种基于增强反投影的全色锐化后处理方法。

背景技术：

遥感图像的全色锐化就是在已知低空间分辨率的多光谱图像和高空间分辨率、低光谱分辨率的全色图像的基础上，恢复其兼具高空间和高光谱分辨率的多光谱图像，进而提升多光谱图像对地物的精细表达能力。全色锐化问题的研究源于当前的高分卫星受技术所限，对同一场景仅能产生一幅高空间分辨率的全色图像(panchromaticimage)和具有多个波段的低空间分辨率多光谱图像(multispectralimage)，而不能产生同时具备高空间分辨率和高光谱分辨率的图像。高空间分辨率有利于对地物的精准定位和精确描述，高光谱分辨率能刻画物质属性，有利于对地物的分类和识别。结合二者优势的高空间分辨率多光谱图像在目视判读、变化监测、提高分类和目标识别精度等诸多方面具有重要的应用价值。为了满足这一需求，遥感图像全色锐化概念被提出并迅速成为遥感图像处理领域的研究热点。

全色锐化问题属于图像反问题，是数学与信息科学的交叉研究领域。从摄影成像的角度出发，全色锐化数学上可归结为如下典型的不适定反问题模型：

其中，为观测得到的低空间分辨率多光谱图像的第k个波段图，为理想的高空间分辨率多光谱图像的第k个波段图，b为下采样矩阵，hk为第k个波段所对应的模糊算子，e为噪声或模型误差，n为多光谱图像的波段数。

为了利用观测得到的高空间分辨率全色图像x^pan，全色锐化问题从统计学的角度也可以看成为贝叶斯统计推断问题：

其中，为所有高空间分辨率多光谱波段图，为所有低空间分辨率多光谱波段图，p(·)为概率密度函数。

基于上述两种模型，学者们提出了一些典型的全色锐化方法，例如：贝叶斯数据融合法、基于全变差正则化的方法和基于稀疏表示的方法，等等。这些方法源于机器学习和图像处理领域的最新研究成果，虽然取得了相对较好的效果，但是，其求解复杂度较高，且进一步提升性能的空间有限。从信号处理的角度出发，低分辨率图像相对于其高分辨率图像而言，丢失了部分高频细节信息，即：

其中，为第k个波段图所丢失的空间细节。因此，为恢复高分辨率多光谱图像，需要对其低分辨率多光谱图像注入一些高频空间细节

基于细节注入的思想，全色锐化方法大致可以分为两类。一类是成分替代(componentsubstitution:cs)法，其代表性方法有主成分分析(principalcomponentanalysis:pca)，ihs变换,波段依赖空间细节法等等；另一类是多分辨率分析(multi-resolutionanalysis:mra)法，其典型方法有：mtf法，hpf法，基于光滑滤波的mtf法，上下文相关决策的广义laplace金字塔，等等。cs类方法可以极大提高多光谱图像的空间分辨率但常常伴随光谱扭曲，而mra类方法能够较好的保持多光谱图像的光谱信息但其重构图像的空间细节注入往往不足。为了进一步提升两类方法的效果，一些学者尝试自适应的调节空间细节的注入量，例如：自适应部分成分替代法，自适应正交变换法，基于二叉树分割的上下文自适应法，等等；还有一些学者尝试把两类方法结合。这些改进在性能上虽有所提高但很难找到提高空间分辨率与保持光谱信息的平衡点。总之，这些方法仍存在光谱扭曲和空间扭曲问题。

为了解决全色图像与多光谱图像光谱不匹配导致的光谱扭曲问题，一方面，学者利用直方图匹配(histogrammatching：hm)技术对全色锐化融合过程进行预处理，进而一定程度上改善了光谱扭曲问题。另一方面，实验发现高通调制(high-passmodulation：hpm)技术可以进一步提升融合图像的空间细节信息。但是，这两项技术往往导致其融合结果不满足空间一致性，即对融合后的图像进行滤波下采样后与原始观测到的低分辨率图像不一致。因此，为满足一致性要求，并考虑真实传感器的物理参数信息，即真实传感器的奈奎斯特截断频率，可以引入调制传递函数匹配滤波器(modulationtransferfunction(mtf)-matchedfilter)，进而可以在提升融合图像空间分辨率的同时，尽可能多的保留原始多光谱图像的光谱特征。

尽管直方图匹配hm、高通调制hpm以及mtf-匹配滤波器在提升全色锐化的结果方面非常有效，但是，它们仅仅是作为算法的预处理或中间步骤而没有用于后处理。此外，这三项技术在全色锐化方法中通常只是利用其中的一种，它们的结合利用至今是一个空白。

技术实现要素：

本发明的目的在于提供一种基于增强反投影的全色锐化后处理方法，以解决上述存在的技术问题。本发明可应用于全色锐化的后处理当中，能在保持光谱图像光谱信息的同时极大提升融合结果的空间分辨率。

为达到上述目的，本发明采用以下技术方案：

一种基于增强反投影的全色锐化后处理方法，包括以下步骤：

步骤1，对全色锐化处理获得的融合结果进行增强处理；增强处理包括调整光谱信息和提升空间分辨率；

步骤2，对步骤1增强处理后的融合结果进行重构误差反投影；包括以迭代的方式计算重构误差，直至重构误差达到预设收敛条件，获得后处理后的全色锐化融合结果。

进一步的，步骤1中调整光谱信息具体包括：通过直方图匹配技术调整融合结果的光谱信息。

进一步的，通过直方图匹配技术调整融合结果的光谱信息的具体步骤包括：

将高分辨率全色图像x^pan与第k个低分辨率多光谱波段进行直方图匹配，得到匹配的全色图像计算公式为：

其中，和分别表示全色图像和第k个低分辨率多光谱波段的均值，和分别表示全色图像和第k个低分辨率多光谱波段的标准差。

进一步的，步骤1中提升空间分辨率的过程包括：通过高通调制技术弥补注入空间细节信息；其中，提升空间分辨率的过程设置于调整光谱信息过程之后。

进一步的，步骤1中提升空间分辨率的过程具体包括：根据直方图匹配后的高分辨率全色图像对每个初始hr多光谱波段进行调制，第1次迭代后得到高分辨率多光谱波段计算公式为：

其中，是第k次迭代由的低通滤波获得的低分辨率全色图像。

进一步的，步骤2中，利用mtf匹配滤波器，对增强处理后的融合结果进行重构误差反投影处理。

进一步的，步骤2中的重构误差反投影处理具体包括：

步骤2.1，以成像模型为基础计算低分辨率下的重构误差；

步骤2.2，将重构误差反投影给融合结果，调整空间分辨率和光谱信息；

步骤2.3，重复步骤2.1和2.2，直至重构误差达到预设收敛阈值，输出反投影处理后的融合结果。

进一步的，步骤2利用mtf匹配滤波器，对增强的融合结果进行重构误差反投影处理；对每一个多光谱波段迭代计算以下两个步骤：

步骤2.1：计算第t次迭代原始低分辨率多光谱波段与反投影的低分辨率多光谱波段的重建误差计算式为：

其中，是第t次迭代第k个多光谱波段相应的mtf匹配滤波函数，*表示卷积算子，↓r为下采样算子；

步骤2.2：反投影重构误差给融合结果以调整空间分辨率和光谱信息，调整后的多光谱波段为：

其中，是相应mtf匹配滤波器的反投影滤波器，*表示卷积算子，↑r为上采样算子；

步骤2.3，重复步骤2.1和2.2，直至重构误差达到预设收敛阈值，输出反投影处理后的融合结果。

进一步的，重构误差的预设收敛阈值为0.005。

进一步的，步骤1中，全色锐化处理方法为bdsd、gfpca、gsa、mf、nlihs或sfim。

与现有技术相比，本发明具有以下有益效果：

本发明基于融合结果的增强和反投影处理，可应用于全色锐化处理融合结果的后处理，能在保持光谱图像光谱信息的同时极大提升融合结果的空间分辨率；后处理增强的融合图像可以更加精确的识别图像中的地物信息，可以更好得应用于目视判读、变化监测、提高分类和目标识别精度等；本发明是遥感图像全色锐化中独立、互补的后处理模块，不需要对现有的全色锐化方法做任何改进，可以有效提高融合结果的质量，可进一步开拓了基于遥感光谱图像的应用。

进一步的，直方图匹配hm技术可以有效减少光谱扭曲；高通调制hpm技术可以提高图像空间质量，尤其是当要匹配或者调制的图像与原始图像很相似时；对真实成像传感器物理参数(与调制传递函数相匹配的滤波器，即mtf匹配滤波器)的利用也必将提升算法的性能。本发明将hm、hpm、mtf匹配滤波器有效地结合在反投影过程当中，并应用于全色锐化的后处理当中，充分利用这三项技术的各自特性使得融合图像的质量得到极大的提升。

附图说明

图1是针对ikonos数据集通过八种全色锐化方法得到的融合图像主观对比采用(结果位于子图右侧)和不采用(结果位于子图左侧)本发明的ebp后处理方法两种结果的差异；图1中，(a)图为bdsd方法，(b)图为gfpca方法，(c)图为gsa方法，(d)图为mf方法，(e)图为nlihs方法，(f)图为pracs方法，(g)图为sfim方法，(h)图为pnn方法；

图2是针对quickbird数据集通过八种全色锐化方法得到的融合图像主观对比采用(结果位于子图右侧)和不采用(结果位于子图左侧)本发明的ebp后处理方法两种结果的差异；图2中，(a)图为bdsd方法，(b)图为gfpca方法，(c)图为gsa方法，(d)图为mf方法，(e)图为nlihs方法，(f)图为pracs方法，(g)图为sfim方法，(h)图为pnn方法；

图3是针对worldview-2数据集通过八种全色锐化方法得到的融合图像主观对比采用(结果位于子图右侧)和不采用(结果位于子图左侧)本发明的ebp后处理方法两种结果的差异，图3中，(a)图为bdsd方法，(b)图为gfpca方法，(c)图为gsa方法，(d)图为mf方法，(e)图为nlihs方法，(f)图为pracs方法，(g)图为sfim方法，(h)图为pnn方法；

图4是针对geoeye-1数据集通过八种全色锐化方法得到的融合图像主观对比采用(结果位于子图右侧)和不采用(结果位于子图左侧)本发明的ebp后处理方法两种结果的差异；图4中，(a)图为bdsd方法，(b)图为gfpca方法，(c)图为gsa方法，(d)图为mf方法，(e)图为nlihs方法，(f)图为pracs方法，(g)图为sfim方法，(h)图为pnn方法。

具体实施方式

下面结合附图和具体实施例对本发明作进一步详细说明。

本发明的一种基于增强反投影的全色锐化后处理方法，具体包括以下阶段和步骤：

阶段1-增强阶段：此阶段为了进一步调整光谱信息和提升空间分辨率，具体包括以下两个步骤，

步骤1.1：通过直方图匹配(hm)技术进一步调整融合结果的光谱信息以减轻光谱扭曲；

步骤1.2：通过高通调制(hpm)技术进一步提升融合结果的空间分辨率以弥补注入的空间细节信息；

阶段2-反投影阶段：阶段1的增强使得空间分辨率得到了极大提升，但也带来了光谱和空间扭曲以及融合结果不满足空间一致性等问题，为此，需要进一步调节融合结果。此外，大多数全色锐化算法在其设计当中都未考虑真实传感器的物理参数信息，特别地，如真实传感器的奈奎斯特截断频率。综上，我们以成像模型为基础，利用与真实调制传递函数相匹配的滤波器(mtf匹配滤波器)，对增强阶段的融合结果进行重构误差反投影过程。此过程以迭代的方式进行以下两个步骤：

步骤2.1，以成像模型为基础计算低分辨率下的重构误差；

步骤2.2，将重构误差反投影给融合结果以调整空间分辨率和光谱信息；

直到重建误差足够小为止，并在最后输出反投影阶段的结果。

直方图匹配hm技术可以有效减少光谱扭曲；高通调制hpm技术可以提高图像空间质量，尤其是当要匹配或者调制的图像与原始图像很相似时；对真实成像传感器物理参数(与调制传递函数相匹配的滤波器，即mtf匹配滤波器)的利用也必将提升算法的性能。本发明的基于增强反投影(enhancedback-projection：ebp)的全色锐化后处理算法将hm、hpm、mtf匹配滤波器有效地结合在反投影过程当中，并应用于全色锐化的后处理当中，从而充分利用这三项技术的各自特性使得融合图像的质量得到极大的提升，能在保持光谱图像光谱信息的同时极大提升融合结果的空间分辨率。本发明是遥感图像全色锐化中独立、互补的后处理模块，不需要对现有的全色锐化方法做任何改进，易于实现和应用。本发明的ebp首先利用hm和hpm技术增强其它算法的融合结果，然后以一种迭代的方式将重构误差反投影回去进而调整融合图像的质量。为了保持融合结果的空间一致性，不同于传统的反投影技术，ebp利用了mtf匹配滤波器来引导反投影过程。针对八种经典的全色锐化算法，在四种卫星数据上进行实验，结果表明我们提出的ebp作为全色锐化问题的后处理算法在不需要改动传统算法的基础上，可以极大提升融合图像的质量。对几种经典全色锐化方法，例如：bdsd、gfpca、gsa、mf、nlihs、sfim，利用本后处理方法进行后处理在多个数据集上其融合效果均有不同程度的提升。

实施例

数据集的选取：

请参阅表1，实施例中选择ikonos、quickbird、worldview-2、geoeye-1四个卫星数据集作为实验数据。

ikonos数据集的全色图像分辨率为1米，多光谱图像的分辨率为4米，多光谱图像由蓝波段、绿波段、红波段、近红外1波段四个波段构成，主要包含大块物体，如山、农场、马路以及地震过后的房屋等；quickbird数据集的全色图像分辨率为0.7米，多光谱图像的分辨率为2.8米，多光谱图像由蓝波段、绿波段、红波段、近红外波段四个波段构成，主要包含巨大的建筑，如摩天大厦、商业工业建筑，以及许多小物体，如汽车，操场等；worldview-2数据集的全色图像分辨率为0.5米，多光谱图像的分辨率为2米，多光谱图像由蓝波段、绿波段、红波段、近红外1波段、沿海波段、黄波段、红边波段、近红外2波段八个波段构成，图像主要为一些有阴影的建筑和树木；geoeye-1数据集全色图像分辨率为0.46米，多光谱图像分辨率为1.84米，多光谱图像由蓝波段、绿波段、红波段、近红外1波段、沿海波段、黄波段、红边波段、近红外2波段八个波段构成，主要包含一些同质和异质区域，并具有大量精细的空间细节信息。

本实施例中，分别从以上数据集中截取大小为512×512的全色图像，大小为128×128的多光谱图像作为实验数据。

表1四种数据集的数据信息

全色锐化：

本实施例中选择八种经典的全色锐化方法：bdsd、gfpca、gsa、mf、nlihs、pracs、sfim、pnn，分别对四个数据集的高空间分辨率全色图像与低分辨率多光谱图像进行全色锐化，得到高空间、高光谱分辨率的多光谱图像。由于真实的高分辨率多光谱图像无法通过卫星获取，为了对比直接全色锐化与进行后处理融合结果的差异，分别在全尺度下和降尺度下进行两组实验。

本发明的ebp算法进行后处理：

对于全色锐化所得融合结果采用本发明提出的ebp算法进行后处理，包括以下两个阶段。

阶段1-增强阶段：此阶段进一步调整光谱信息和提升空间分辨率；

已知全色图像x^pan、低分辨率多光谱波段以及由全色锐化方法得到的原始hr多光谱波段增强阶段进行以下两个步骤：

步骤1.1：通过直方图匹配技术进一步调整融合结果的光谱信息以减轻光谱扭曲；

将高分辨率全色图像x^pan与第k个低分辨率多光谱波段进行直方图匹配，得到匹配的全色图像公式如下：

其中，分别表示全色图像x^pan、第k个低分辨率多光谱波段的均值，分别表示全色图像x^pan、第k个低分辨率多光谱波段的标准差。

步骤1.2：通过高通调制技术弥补注入的空间细节信息进一步提升融合结果空间分辨率；

根据直方图匹配后的高分辨率全色图像对每个初始hr多光谱波段进行调制，得到1次迭代后的高分辨率多光谱波段调制方法如下：

其中，是由的低通滤波获得的低分辨率全色图像。

阶段2-反投影阶段：利用mtf匹配滤波器，对增强阶段的融合结果进行重构误差反投影过程。对每一个多光谱波段迭代计算以下两个步骤：

步骤2.1：计算第t次迭代原始低分辨率多光谱波段与反投影的低分辨率多光谱波段的重建误差计算式为：

其中，是第t次迭代的第k个多光谱波段相应的mtf匹配滤波函数，*表示卷积算子，↓r为下采样算子，下采样算子在全色锐化中通常取为4。

步骤2.2：反投影重构误差给融合结果以调整空间分辨率和光谱信息，调整后的多光谱波段为：

其中，是相应mtf匹配滤波器的反投影滤波器，*表示卷积算子，↑r为上采样算子。

重复以上两步直到重建误差达到规定阈值为止，此时原始低分辨率多光谱图像与反投影的低分辨率多光谱图像差异也充分小，输出反投影阶段的结果，得到由ebp算法后处理的融合图像。

评价指标：

一般来说，全色锐化方法的性能好坏有两种评价方式：主观评价和客观评价。主观评价主要是观察融合后的多光谱图像的彩色显示，即rgb三个波段的真彩色显示，分析其色彩是否自然，是否和原始多光谱rgb显示图像相似，融合后的物体轮廓是否更加清晰等等；客观评价需要一些量化指标来评价。这些量化指标中，有的是将真实高分辨率多光谱图像作为参照，根据一些量化指标，如相关系数(correlationcoefficient：cc)、空间相关系数(spatialcorrelationcoefficient：scc)、均值平方根误差(rootmeansquarederror：rmse)、ergas、光谱角度映射器(spectralanglemapper：sam)、q4/q8等进行评价，但是由于实际中无法获得真实高分辨率多光谱图像，这使得这些需要参照物的指标只能在降尺度下进行评估，即将从卫星获得的高分辨率全色图像与低分辨率多光谱图像分别降至低尺度下(一般尺度缩放因子为4)，再进行全色锐化，得到的融合结果与原始多光谱图像进行对比。有的指标无需参照物，如无参考质量指标(qualitywithnoreference：qnr)。qnr由光谱扭曲指标dλ和空间扭曲指标ds构成，它可以在全尺度下对融合图像进行量化评估。

本实施例选取cc、rmse、ergas、sam、q4/q8这五个指标在降尺度下对融合结果的性能进行评价，选取qnr指标在全尺度下评价融合结果。

主观评价结果分析：

四个数据集原始融合结果与ebp算法后处理过的融合结果的对比分别见图1、图2、图3和图4。

参阅图1至图4可知，ebp后处理方法可以显著提高全色锐化所得融合图像的空间细节信息，尤其是对数据集ikonos、quickbird、worldview-2，在采用gfpca和nlihs两种全色锐化方法时，对融合结果后处理后，大大增强了原始融合图像模糊区域的空间细节信息，新的融合图像显示出更多的建筑、房屋等信息。对比其他六种全色锐化方法得到的融合图像，本发明提出的ebp算法的确会增加或保持图像空间细节信息，并有效减少了光谱扭曲问题。

客观评价结果分析：

全尺度下：

根据qnr无参照指标对融合结果在全尺度下进行评价，计算四个数据集八种全色锐化方法融合结果的qnr值，结果见表2。表2表明，对于ikonos、quickbird、worldview-2这三个数据集，将全色锐化方法得到的融合结果应用本发明提出的ebp方法后，qnr指标值都高于或者等于未进行后处理时的指标值(见表2加粗指标值)，有效提高了融合图像的空间质量并减少了光谱扭曲。

表2、在全尺度下分别对比四个数据集(ikonos、quickbird、worldview-2、geoeye-1)全色锐化结果采用(√)与不采用(×)ebp方法所得qnr无参照指标值

*thesymbol-denotesnoresults.

降尺度下：

因为理想的高分辨率多光谱图像不能得到，故实验在降尺度下进行对比。首先分别将四个数据集中的全色图像分辨率降低4倍，此时图像大小为128×128，多光谱图像分辨率也降低4倍，图像大小变为32×32。然后将大小为128×128的原始多光谱图像作为参考图像，分别与八种全色锐化方法得到融合图像进行对比，最后根据两组融合结果的五个指标值进行评价，结果见表3。

表3、在降尺度下分别对比四个数据集(ikonos、quickbird、worldview-2、geoeye-1)全色锐化结果采用(√)与不采用(×)ebp方法所得cc、ergas、sam、rmse、q4指标值

*thesymbol-denotesnoresults.

表3表明，对于worldview-2、geoeye-1数据集，由ebp算法对八种全色锐化方法的融合结果后处理，以上五个指标值的结果都优于或等于未进行后处理时的指标值(表3加粗数据)。因此本发明提出的ebp方法对提高worldview-2、geoeye-1数据集融合图像的空间质量与细节信息很有优势。对于ikonos、quickbird数据集，本发明除了pracs和pnn全色锐化方法后处理的图像效果不好，而对于其余方法的后处理，可有效提高数据集融合图像的空间信息。

实验结论：

第一，对于数据集worldview-2和geoeye-1，本发明提出的ebp后处理算法可以显著增强八种全色锐化方法(bdsd、gfpca、gsa、mf、nlihs、pracs、sfim、pnn)所得融合图像的空间质量，并减少光谱扭曲问题。对于ikonos和quickbird数据集，ebp算法可以增强六种全色锐化方法(bdsd、gfpca、gsa、mf、nlihs、sfim)的融合结果，得到具有更多空间细节信息的融合图像。第二，后处理增强的融合图像可以更加精确的识别图像中的地物信息，因此可以更好的应用于目视判读、变化监测、提高分类和目标识别精度等。第三，全色锐化，即把高分辨率的全色图像和低分辨率的多光谱图像融合，是得到高分辨率ms图像的有效手段，而本发明提出的ebp算法作为全色锐化独立、互补的后处理模块，不需要对现有的全色锐化方法做任何改进，可以有效提高融合结果的质量，因此进一步开拓了基于遥感光谱图像的应用。