基于稀疏张量近邻嵌入的遥感图像融合方法与流程

本发明属于图像处理技术领域，特别涉及一种遥感图像融合方法，可用于遥感探测、安全导航、医学图像分析、反恐检查和环境保护。

背景技术：

遥感图像融合是对多遥感器获得的图像数据和其他信息处理的过程。其中研究最广的是多光谱图像与全色图像融合。多光谱图像，色彩信息丰富，但空间分辨率较低，一般为四个波段。全色图像，空间分辨率很高，细节清晰，但缺少颜色信息，一般指使用0.5μm到0.75μm左右的单波段。在目标识别时，色彩信息和目标的空间分辨率对于分辨地物特征都具有非常重要的作用。

遥感图像融合是对多遥感器获得的图像数据处理的过程。目的是将两类图像通过一定的规则或方法有机的结合起来，得到一幅具有高光谱分辨率和高空间分辨率的图像。多光谱与全色图像融合方法可以被归结为三大类：(1)分量替换CS_Based融合框架；(2)基于多尺度分析工具MRA_Based的融合框架；(3)基于模型Model_Based的融合框架。

基于模型Model_Based的融合算法在光谱信息的保持方面仍然存在不足之处：1)忽略了波段之间的联系导致了融合结果的光谱扭曲；2)忽略了对高阶统计特性的应用。

近年来，压缩感知和稀疏表示的理论的发展，使得基于退化模型为理论的模型复原的图像融合方法得以发展。参见Jiang C,Zhang H,Shen H,et al.Two-step sparse coding for the pan-sharpening of remote sensing images[J].Selected Topics in Applied Earth Observations and Remote Sensing,IEEE Journal of,2014,7(5):1792-1805。该算法为两步稀疏编码多光谱图像锐化方法。第一步稀疏编码过程是低分辨多光谱图像的图像块由对应位置的低分辨全色图像块线性表示，最终求得一个系数。第二步稀疏编码是通过低分辨字典来稀疏表示残差，然后得到残差的稀疏向量。两步稀疏编码结束后，两个系数被整合成一个向量。融合结果由最终整合的稀疏向量和字典原子相乘得到。该算法由于将多光谱图像分波段融合，没能充分考虑对波段间的信息充分利用，使得融合结果产生了色彩失真及光谱扭曲。

技术实现要素：

本发明目的在于针对上述已有技术的不足，提出一种基于稀疏张量近邻嵌入的遥感图像融合方法，以充分利用多光谱图像波段间的信息，保持融合结果的光谱信息，减少融合后的多光谱图像的颜色失真及光谱扭曲现象。

本发明的技术方案是这样实现的：

一.技术原理

“张量”是高阶数据处理的有效模型。张量作为矢量与矩阵的高维推广，在多维信号处理时张量形式能够很好地保持多维信号的空间结构。与基于向量的方法相比，基于张量的方法具有如下优点：(a)能够保持更多的张量数据的结构信息，从而使得数据蕴含的信息得到充分的利用，充分利用多光谱波段间的信息；(b)通过设计与张量分解有关的模型和算法，可以改变以往融合算法中多图像量化操作的过程，可以在融合的过程中有效的利用多光谱数据的高阶统计特性。

本发明结合张量高阶数据处理的特点，克服基于模型的融合框架中的量化的操作，利用张量能够很好地保持多维信号的空间结构的特点充分利用多光谱图像各个波段间的联系，将低分辨多光谱图像在张量的形式下构建张量形式的近邻嵌入模型，将低分辨的多光谱图像在低分辨的近邻嵌入字典下稀疏表示后，将稀疏系数张量映射到高维流形结构上得到高分辨的多光谱图像。

二.实现方案

基于上述原理，本发明的实现方案包括如下：

(1)分别输入低分辨率多光谱图像M和高分辨率全色图像P；

(2)对低分辨多光谱图像M进行上采样，得到上采样多光谱图像M1，该M1与全色图像P的分辨率尺寸相同；

(3)对高分辨全色图像P进行下采样，得到下采样全色图像P1，该P1与低分辨多光谱图像M的分辨率尺寸相同；

(4)对下采样全色图像P1进行上采样，得到上采样全色图像P2，该P2与高分辨全色图像P的分辨率尺寸相同；

(5)对上采样多光谱图像M1进行分块，得到三维的多光谱张量块M；

(6)根据上采样全色图像P2，获得低分辨模1矩阵D_L1：

6a)将上采样全色图像P2进行分块，得到上采样全色图像小块P3，对P3进行近邻扩充操作，得到扩充后的上采样全色图像I_L；

6b)对扩充后的上采样全色图像I_L划窗取块，并对其进行列向量化，得到低分辨模1矩阵D_L1；

(7)根据高分辨全色图像P，获得高分辨模1矩阵D_H1：

7a)将高分辨全色图像P进行分块，得到全色图像小块P4，对P4进行近邻扩充操作，得到扩充后的全色图像I_H；

7b)将扩充后的全色图像I_H划窗取块，并经过列向量化操作，得到高分辨模1矩阵D_H1；

(8)分别优化低分辨模1矩阵D_L1中的原子、高分辨模1矩阵D_H1中的原子：

8a)根据多流形理论，分别对低分辨模1矩阵D_L1和高分辨模1矩阵D_H1中的原子进行单位化欧氏距离操作，使得这两个矩阵D_L1和D_H1中原子在同一流形结构上；

8b)根据多流形理论，分别对低分辨模1矩阵D_L1和高分辨模1矩阵D_H1中的原子进行添加惩罚正则项的优化求解操作，使得这两个矩阵D_L1和D_H1中的原子在同一切空间平面上；

(9)对多光谱张量块M进行量化操作，得到低分辨模2矩阵D_L2和低分辨模3矩阵D_L3，及高分辨模2矩阵D_H2和高分辨模3矩阵D_H3；

(10)用低分辨模1矩阵D_L1、低分辨模2矩阵D_L2和低分辨模3矩阵D_L3进行组合，得到低分辨多模字典LD；

(11)将高分辨模1矩阵D_H1、高分辨模2矩阵D_H2和高分辨模3矩阵D_H3进行组合，得到高分辨多模字典HD；

(12)利用张量基追踪算法，求解多光谱张量块M在低分辨多模字典LD下的稀疏张量系数A；

(13)将高分辨多模字典HD与稀疏张量系数A进行模乘，获得融合的多光谱图像F1；

(14)利用空域残差补偿法更新融合的多光谱图像F1，获得高分辨的多光谱图像F。

本发明与现有的技术相比具有以下优点：

第一，本发明将低分辨多光谱图像以张量的形式进行表示，保证了在计算时各个波段间的信息的完整性，克服了量化操作对光谱信息所造成的影响，改善了光谱扭曲及色彩失真的现象；

第二，本发明构建了稀疏张量近邻嵌入模型，在张量的形式下，利用了流形假设中的局部相似性特征，使得融合图像的空间信息更加精确；

第三，本发明对近邻字典进行优化，优化后的字典原子与对应的多光谱图像块不仅在同一流形结构上，且在同一个切空间上，提高了张量的近邻字典的性能，增强了融合结果的空间分辨率。

仿真结果表明，本发明得到的融合后多光谱图像比现有技术的融合图像能更好的改善光谱扭曲及颜色失真的现象，保持了融合图像的光谱分辨率，同时提高了空间分辨率。

附图说明

图1为本发明的实现流程图；

图2为本发明的仿真结果图。

具体实施方法

以下结合附图为本发明做进一步详细描述：

参照图1，本发明的具体实施方式如下：

步骤1，分别输入低分辨率多光谱图像M和高分辨率全色图像P，

(1.1)本发明实施例中输入的低分辨多光谱图像M大小为64×64×4，分辨率为2m；高分辨率的全色图像P大小为256×256，分辨率为0.5m；

(1.2)将低分辨的多光谱图像M上采样至与高分辨率全色图像P的相同分辨率尺寸，得到上采样的多光谱图像M1，本实例中M1的大小为256×256×4。

步骤2，输入高分辨率全色图像P，得到下采样全色图像P1和上采样全色图像P2。

(2.1)输入高分辨率全色图像P，并对其下采样得到大小为64×64的下采样全色图像P1；

(2.2)对下采样全色图像P1上采样，得到大小为256×256的上采样全色图像P2。

步骤3，对上采样的多光谱图像M1分块，构成三维的多光谱张量块M。

(3.1)上采样多光谱图像M1包含四个波段，对M1的每一个波段按5×5的大小划窗取块，得到四个波段的图像块分别为B1、B2、B3、B4；

(3.2)对四个波段的图像块B1、B2、B3、B4按照第三维方向组合，得到三维的多光谱张量块M。

步骤4，输入上采样全色图像P2，得到低分辨模1矩阵D_L1。

(4.1)对上采样全色图像P2分块，得到上采样全色图像小块P3；

(4.2)对上采样全色图像小块P3以大小为21进行邻域扩充，得到大小为21×21的扩充后上采样全色图像I_L，且I_L的中心为P3；

(4.3)对每一个扩充后的上采样全色图像I_L用5×5大小的窗口，按“Z”字形划窗取块，将每一个块列向量化后，得到大小为25×1的上采样全色图像小块列向量；

(4.4)对上采样全色图像小块列向量进行组合，得到大小为25×289的低分辨模1矩阵D_L1。

步骤5，根据高分辨全色图像P，得到高分辨模1矩阵D_H1。

(5.1)对高分辨全色图像P按5×5的大小划窗分块，得到全色图像小块P4；

(5.2)对全色图像小块P4以大小为21进行邻域扩充，得到大小为21×21的扩充后全色图像I_H，且I_H的中心为P4；

(5.3)对每一个扩充后的全色图像I_H用5×5大小的窗口，以1为步长，按“Z”字形划窗取块，并对每一个块列进行向量化，得到大小为25×1的高分辨全色图像小块列向量；

(5.4)对高分辨全色图像小块列向量进行组合，得到大小为25×289的高分辨模1矩阵D_H1。

步骤6，计算低分辨模1矩阵D_L1中原子的单位化欧氏距离。

(6.1)对上采样全色图像小块P3列向量化，得到上采样全色图像小块向量t_P3；

(6.2)使用公式||D_L1-t_P3||₂计算低分辨模1矩阵D_L1的原子与上采样全色图像小块向量t_P3的模值；

(6.3)使用公式计算低分辨模1矩阵D_L1中的原子的单位化欧氏距离。

步骤7，计算高分辨模1矩阵D_H1中原子的单位化欧氏距离。

(7.1)对全色图像小块P4列向量化，得到全色图像小块向量t_P4；

(7.2)使用公式||D_H1-t_P4||₂计算高分辨模1矩阵D_H1的原子与全色图像小块向量t_P4的模值；

(7.3)使用公式计算高分辨模1矩阵D_H1的原子的单位化欧氏距离。

步骤8，优化低分辨模1矩阵D_L1中的原子。

(8.1)使用公式t_P3＝D_L1C_L计算低分辨权值C_L，其中，t_P3为上采样全色图像小块向量，D_L1为低分辨模1矩阵；

(8.2)使用公式：计算低分辨惩罚正则项Q_L，

其中，D_L1为低分辨模1矩阵，s.t表示求解的限制条件，||Q_LC_L||₁为Q_LC_L的一范数，||D_L1C_L||₂为D_L1C_L的二范数，T表示矩阵的转置，ε为重构误差；

(8.3)将高分辨惩罚正则项Q_L中元素编号映射到低分辨模1矩阵D_L1；

(8.4)在矩阵D_L1中删除零元素编号的原子，将剩下非零元素编号的原子作为优化后的低分辨模1矩阵D_L1原子。

步骤9，优化高分辨模1矩阵D_H1中的原子。

(9.1)使用公式t_P4＝D_H1C_H计算高分辨权值C_H，其中，t_P4为全色图像小块向量，D_H1为高分辨模1矩阵；

(9.2)使用公式：计算高分辨惩罚正则项Q_H，

其中，D_H1为高分辨模1矩阵，||Q_HC_H||₁为Q_HC_H的一范数，||D_H1C_H||₂为D_H1C_H的二范数；

(9.3)将高分辨惩罚正则项Q_H中元素编号映射到高分辨模1矩阵D_H1；

(9.4)在矩阵D_H1中删除零元素编号的原子，将剩下非零元素编号的近邻原子作为优化后的高分辨模1矩阵D_H1原子。

步骤10，根据多光谱张量块M得到模矩阵。

(10.1)对多光谱张量块M按照模2方向进行量化，得到量化矩阵，将该量化矩阵同时作为低分辨模2矩阵D_L2和高分辨模2矩阵D_H2；

(10.2)对多光谱张量块M按照模3方向进行量化，得到量化矩阵，将该量化矩阵同时作为低分辨模3矩阵D_L3和高分辨模3矩阵D_H3。

步骤11，利用模矩阵得到多模字典。

(11.1)组合低分辨模1矩阵D_L1、低分辨模2矩阵D_L2和低分辨模3矩阵D_L3，得到低分辨多模字典LD；

(11.2)组合高分辨模1矩阵D_H1、高分辨模2矩阵D_H2和高分辨模3矩阵D_H3，得到高分辨多模字典HD。

步骤12，利用多光谱张量块M、低分辨多模字典LD和高分辨多模字典HD，得到融合的多光谱图像F1。

(12.1)使用公式：计算多光谱张量块M在低分辨多模字典LD下的稀疏张量系数A，

其中，A为待求解的高分辨率稀疏张量系数，||A||₁为稀疏张量系数A的一范数，LD为低分辨率多模字典，为A×₁LD-M的二范数平方；

(12.2)使用公式F1＝A×HD计算高分辨多模字典HD与稀疏张量系数A模乘得的融合的多光谱图像F1。

步骤13，更新融合的多光谱图像F1，获得高分辨的多光谱图像F。

(13.1)对融合的多光谱图像F1分块，得到融合多光谱图像小块M_F；

(13.2)使用公式R＝M-M_F计算融合多光谱图像小块M_F与多光谱张量块M的残差R；

(13.3)使用公式计算残差R在低分辨多模字典LD下的稀疏系数A_R；

(13.4)根据稀疏系数A_R与高分辨多模字典HD，得到融合后残差R_F＝A_R×HD；

(13.5)根据融合后残差R_F与融合多光谱图像小块M_F，得到高分辨的多光谱图像F＝R_F+M_F。

本发明的效果可以通过以下仿真实验作进一步说明。

1.实验仿真环境为MATLAB R2014a on PC with Intel(R)core(TM)/2.50G/2.50G；

2.仿真内容：

仿真1，使用本发明方法，对低分辨的多光谱图像和高分辨的全色图像融合，得到高分辨的多光谱图像，结果如图2。其中：

图2(a)是低分辨的多光谱图像，图2(b)是高分辨的全色图像，图2(c)是参考图像，这三幅图像选自GeoEye，图2(d)是采用本发明对图2(a)和图2(b)融合后获得的高分辨多光谱图像。

由图2(d)与图2(a)的比较可见，图2(a)的细节不清晰，如道路的清晰度，房屋的色彩等均不清楚，而图2(d)在空间分辨率、光谱信息、视觉效果上明显优于图2(a)。

仿真2，为了证明本发明的效果，用本发明的方法和现有的代表性融合算法，分别对图2(a)和图2(b)的图像融合，融合结果的评价指标对比，如表1。

现有的代表性融合算法包括：基于主成分分析变换的融合算法PCA、基于小波变换的融合算法AWLP、基于小波变换的融合算法CT、基于压缩感知融合算法CS。

表1各方法融合结果的评价指标对比

表1中的常用评价指标包括：相关系数CC、通用图像质量指标UIQI、均方根误差RMSE、整体图像质量指数Q4、光谱弧度SAM、相对整体维数综合误差ERGAS。

从表1数据对比表明，本发明融合结果的相关系数CC、通用图像质量指标UIQI和整体图像质量指数Q4均大于现有代表性融合算法的指标；其均方根误差RMSE、光谱弧度SAM和相对整体维数综合误差ERGAS均小于代表性融合算法的指标；其中，本发明融合结果得到的相关系数CC更接近1，SAM指标最小，表明本发明能更好的保持融合图像的光谱信息，且整体图像质量指数Q4最高，其客观评价结果优于现有技术的客观评价结果。

综上可见，本发明的融合图像优于现有代表性融合算法的结果，且与参考图像最接近，色彩失真程度最小，图像整体质量最优，融合结果最好。