CN104867124A - 基于对偶稀疏非负矩阵分解的多光谱与全色图像融合方法 - Google Patents

Abstract

本发明公开了一种基于对偶稀疏非负矩阵分解的多光谱与全色图像融合方法，主要解决现有技术中空间信息模糊和光谱扭曲的问题。其步骤为：(1)分别输入低空间分辨率的多光谱图像、高空间分辨率的全色图像及光谱退化矩阵、空间退化矩阵；(2)分别将低空间分辨率的多光谱图像和高空间分辨率的全色图像分块列向量化；(3)对列向量化后的图像进行对偶稀疏非负矩阵分解，得到高空间分辨率的字典和具有光谱信息的系数矩阵；(4)将高空间分辨率的字典和具有光谱信息的系数矩阵相乘得到列向量化的高空间分辨率的多光谱图像，并将其恢复为高空间分辨率多光谱图像。本发明能获得准确的空间和光谱信息，可用于目标识别、地物分类和环境监测等遥感领域。

Description

基于对偶稀疏非负矩阵分解的多光谱与全色图像融合方法

技术领域

本发明属于图像处理技术领域，更进一步涉及遥感图像中的图像融合方法，可用于目标识别、地物分类和环境监测等遥感领域。

背景技术

由于光学传感器成像的物理限制，遥感卫星通常只能获得低空间分辨率的多光谱图像和高空间分辨率的全色图像。多光谱图像具有丰富的光谱信息，但空间分辨率较低；全色图像空间分辨率较高，具有丰富的空间信息，但光谱信息较少。为了得到高空间分辨率的多光谱图像，人们采用图像融合技术，综合利用全色图像的空间信息与多光谱图像的光谱信息，最终得到更准确、更全面的地物信息，从而获得更好的视觉效果。

目前，多光谱图像与全色图像融合技术主要可以分为两类：一类是基于成分替代的方法，例如IHS色彩变换法、主成分分析法PCA。由于基于成分替代的方法原理简单，易于操作，因此得到广泛应用。但是这类方法会产生严重的光谱扭曲。另一类是基于注入的方法，例如小波变换法、Laplacian金字塔分解法。这类方法大大减小了图像的光谱扭曲，但是容易产生振铃效应。

中国科学院自动化研究所和北京市信息遥感所在“基于非负矩阵分解的高光谱和可见光图像融合方法”(申请号201210245835.4，申请公布号CN102789639A)的专利申请中，提出将非负矩阵分解技术运用到高光谱图像和可见光图像的融合的方案是：先利用基于稀疏约束的非负矩阵分解方法对高光谱图像进行解混,得到端元光谱曲线；再根据可见光图像的光谱吸收曲线和解混得到的端元光谱曲线,利用有约束的二次规划来求解图像中每个像元的混合系数；最后将解混得到的端元光谱曲线和每个像元的混合系数进行线性混合,得到高空间分辨率和高光谱分辨率的融合图像。该方法虽然能够较好的提高高光谱图像的空间分辨率，但是主要针对高光谱图像与可见光图像的融合，而不能够实现多光谱图像与全色图像的融合，如果将该方案直接用于多光谱图像与全色图像的融合，将产生严重的光谱扭曲。

发明内容

本发明的目的在于针对上述已有技术的不足，提出一种基于对偶稀疏非负矩阵分解的多光谱与全色图像融合方法，以提高空间分辨率，减少光谱扭曲。

为实现上述目的技术方案是：将全色图像与多光谱图像分别作为最终得到的高空间分辨率的多光谱图像光谱退化和空间退化的结果，获取光谱退化矩阵和空间退化矩阵；利用对偶稀疏非负矩阵分解方法的更新规则求解目标函数，得到高空间分辨率的字典和具有光谱信息的系数矩阵；根据高空间分辨率的字典和具有光谱信息的系数矩阵得到高空间分辨率的多光谱图像。具体步骤包括如下：

(1)输入图像和矩阵：

(1a)分别输入低空间分辨率的多光谱图像和高空间分辨率的全色图像；

(1b)分别输入空间退化矩阵D和光谱退化矩阵S；

(2)分别将低空间分辨率的多光谱图像和高空间分辨率的全色图像分块列向量化，获得分块列向量化的低空间分辨率的多光谱图像V_lms和分块列向量化的高空间分辨率的全色图像V_pan；

(3)将分块列向量化的低空间分辨率的多光谱图像V_lms和分块列向量化的高空间分辨率的全色图像V_pan同时进行对偶稀疏非负矩阵分解，得到高空间分辨率的字典W_h和具有光谱信息的系数矩阵H_h：

(3a)构建目标方程：

$\left\{\begin{array}{c}\underset{W_D,H_h,W_h,H_S}{\min }\frac{1}{2}[{\left|\right|V_{\mathrm{lms}}-W_D{H}_h\left|\right|}_2^2+{\left|\right|V_{\mathrm{pan}}-W_h{H}_S\left|\right|}_2^2+\mathrm{\β}[{\left|\right|H_h\left|\right|}_{1/2}+{\left|\right|H_S\left|\right|}_{1/2}\left]\right]\\ \mathrm{subject\; to}{W}_D,H_h,W_h,H_S\≥0;W_D={\mathrm{DW}}_h;H_S=H_{h}S\end{array}\right.$

其中，min[·]为最小化，subject to表示W_D,H_h,W_h,H_S≥0；W_D＝DW_h；H_S＝H_hS为 $\frac{1}{2}[{\left|\right|V_{\mathrm{lms}}-W_D{H}_h\left|\right|}_2^2+{\left|\right|V_{\mathrm{pan}}-W_h{H}_S\left|\right|}_2^2+\mathrm{\β}[{\left|\right|H_h\left|\right|}_{1/2}+{\left|\right|H_S\left|\right|}_{1/2}\left]\right]$ 的约束条件，W_D为低空间分辨率的字典，H_S为无光谱信息的系数矩阵，和为数据的保真项，||H_h||_1/2和||H_S||_1/2为稀疏项，分别表示具有光谱信息的系数矩阵H_h和无光谱信息的系数矩阵H_S的范数，β为正则参数，权衡保真项与稀疏项之间的关系；

(3b)用0到1之间的任意数值对高空间分辨率的字典W_h和具有光谱信息的系数矩阵H_h的每一个元素进行初始化；

(3c)设置最大循环次数n＝1000；

(3d)固定W_D和H_h，利用下式更新W_h和H_S：

$W_h\←W_h.*V_{\mathrm{pan}}{H}_S^T./W_h{H}_S{H}_S^T$

$H_S\←H_S.*W_h^T{V}_{\mathrm{pan}}./(W_h^T{W}_h{H}_S+\frac{\mathrm{\β}}{2}{H_S}^{-\frac{1}{2}})$

其中，T为矩阵转置操作，.*为点乘操作，./为点除操作，←为结果更新操作；

(3e)固定W_h和H_S，利用下式更新W_D和H_h：

$W_D\←W_D.*V_{\mathrm{lms}}{H}_h^T./W_D{H}_h{H}_h^T$

$H_h\←H_h.*W_D^T{V}_{\mathrm{lms}}./(W_D^T{W}_D{H}_h+\frac{\mathrm{\β}}{2}{H_g}^{-\frac{1}{2}});$

(3f)循环(3d)和(3e)，直到循环次数达到预先设定的最大循环次数n；

(4)将高空间分辨率的字典W_h和具有光谱信息的系数矩阵H_h相乘得到分块列向量化的高空间分辨率的多光谱图像V_hms：

V_hms＝W_hH_h；

(5)将分块列向量化的高空间分辨率的多光谱图像V_hms转化为图像块，恢复到图像的对应位置，得到高空间分辨率的多光谱图像。

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

第一，本发明将非负矩阵分解引入到多光谱图像与全色图像融合领域来，克服了现有技术中空间信息模糊和光谱扭曲的不足，使得本发明可以获得准确的空间和光谱信息，有利于后续图像处理技术的使用。

第二，本发明将对偶稀疏正则引入到非负矩阵分解中，克服了现有技术空间信息模糊的不足，使得本发明可以获得边缘信息更加丰富的图像。

第三，本发明采用对偶稀疏的更新规则求解非负矩阵的分解，减小了矩阵初始化对融合结果的影响，使得本发明比现有技术更加鲁棒，融合结果更加稳定。

本发明的技术过程和效果可结合以下附图详细说明。

附图说明

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

图2为用本发明对低空间分辨率多光谱图像和高空间分辨率全色图像融合的结果对比图。

具体实施方法

参照图1，本发明具体实现步骤如下：

步骤1.输入图像和矩阵。

(1a)分别输入低空间分辨率的多光谱图像和高空间分辨率的全色图像；

本发明实施例中，采用Quickbrid卫星于2002年9月在印度孙德尔本斯地区拍摄的低空间分辨率的多光谱和高空间分辨率的全色图像，低空间分辨率的多光谱图像大小为64×64×4，空间分辨率为9.6m，高空间分辨率的全色图像大小为256×256，空间分辨率为2.4m。

(1b)分别输入空间退化矩阵D和光谱退化矩阵S；

空间退化矩阵D，表示最终得到的高空间分辨率的多光谱图像与低空间分辨率的多光谱图像之间的空间退化关系；

光谱退化矩阵S为谱间下采样矩阵，表示最终得到的高空间分辨率的多光谱图像与高空间分辨率的全色图像之间的光谱退化关系。

空间退化矩阵D的表示式为：

$D=\frac{I_{m\×m}\⊗(e_{1\×4}\⊗(I_{m\×m}\⊗e_{1\×4}))}{16}$

其中，I_m×m为单位矩阵，m为小块的尺寸，e_1×4为元素全为1的行向量，为克罗内克积。本发明实施例中，小块的尺寸m＝8，则空间退化矩阵D的大小为64×1024。

光谱退化矩阵S的表示式为：

其中，A_N×N为基矩阵，该矩阵对角线元素为1，其余元素为0，N为低空间分辨率的多光谱图像每个颜色通道中小块的数目，w＝[w₁,w₂,w₃,w₄]为行向量，w₁，w₂，w₃，w₄分别为红，绿，蓝，近红外通道的权值系数，本发明实施例中，w₁＝0.2308，w₂＝0.2315，w₃＝0.1139，w₄＝0.4239，低空间分辨率的多光谱图像每个颜色通道中小块的数目N＝81，则空间退化矩阵S的大小为324×81。

步骤2.分别将低空间分辨率的多光谱图像和高空间分辨率的全色图像分块列向量化，获得分块列向量化的低空间分辨率的多光谱图像V_lms和分块列向量化的高空间分辨率的全色图像V_pan。

(2a)将低空间分辨率的多光谱图像中红颜色通道的图像分为m×m的小块，且每个小块之间重叠1个像素；

(2b)将每个小块的第2列串联到第1列，第3列串联到第2列，以此类推，完成块的列向量化，构成红色通道矩阵；

(2c)重复步骤(2a)和(2b)，构成低空间分辨率的多光谱图像中绿色，蓝色和近红外通道的矩阵；

(2d)将绿色通道矩阵连接到红色通道矩阵，将蓝色通道矩阵连接到绿色通道矩阵，以此类推，构成列向量化的低空间分辨率的多光谱图像V_lms；

(2e)将高空间分辨率的全色图像分为4m×4m的小块，且每个小块之间重叠4个像素；

(2f)将每个小块的第2列串联到第1列，第3列串联到第2列，以此类推，构成列向量化的高空间分辨率的全色图像V_pan。

本发明实施例中，低空间分辨率的多光谱图像被分成8×8的小块，则分块列向化的低空间分辨率的多光谱图像V_lms大小为64×324，高空间分辨率的全色图像被分成32×32的小块，则分块列向量化的高空间分辨率的全色图像V_pan大小为1024×81。

步骤3.将分块列向量化的低空间分辨率的多光谱图像V_lms和分块列向量化的高空间分辨率的全色图像V_pan同时进行对偶稀疏非负矩阵分解，得到高空间分辨率的字典W_h和具有光谱信息的系数矩阵H_h。

(3a)构建目标方程：

$\left\{\begin{array}{c}\underset{W_D,H_h,W_h,H_S}{\min }\frac{1}{2}[{\left|\right|V_{\mathrm{lms}}-W_D{H}_h\left|\right|}_2^2+{\left|\right|V_{\mathrm{pan}}-W_h{H}_S\left|\right|}_2^2+\mathrm{\β}[{\left|\right|H_h\left|\right|}_{1/2}+{\left|\right|H_S\left|\right|}_{1/2}\left]\right]\\ \mathrm{subject\; to}{W}_D,H_h,W_h,H_S\≥0;W_D={\mathrm{DW}}_h;H_S=H_{h}S\end{array}\right.$

其中，min[·]为最小化，subject to表示W_D,H_h,W_h,H_S≥0；W_D＝DW_h；H_S＝H_hS为 $\frac{1}{2}[{\left|\right|V_{\mathrm{lms}}-W_D{H}_h\left|\right|}_2^2+{\left|\right|V_{\mathrm{pan}}-W_h{H}_S\left|\right|}_2^2+\mathrm{\β}[{\left|\right|H_h\left|\right|}_{1/2}+{\left|\right|H_S\left|\right|}_{1/2}\left]\right]$ 的约束条件，W_D为低空间分辨率的字典，H_S为无光谱信息的系数矩阵，和为数据的保真项，||H_h||_1/2和||H_S||_1/2为稀疏项，分别表示具有光谱信息的系数矩阵H_h和无光谱信息的系数矩阵H_S的范数，β为正则参数，权衡保真项与稀疏项之间的关系；

(3b)用0到1之间的任意数值对高空间分辨率的字典W_h和具有光谱信息的系数矩阵H_h的每一个元素进行初始化；

(3c)设置最大循环次数n＝1000；

(3d)固定W_D和H_h，利用下式更新W_h和H_S：

$W_h\←W_h.*V_{\mathrm{pan}}{H}_S^T./W_h{H}_S{H}_S^T$

$H_S\←H_S.*W_h^T{V}_{\mathrm{pan}}./(W_h^T{W}_h{H}_S+\frac{\mathrm{\β}}{2}{H_S}^{-\frac{1}{2}})$

其中，T为矩阵转置操作，.*为点乘操作，./为点除操作，←为结果更新操作；

(3e)固定W_h和H_S，利用下式更新W_D和H_h：

$W_D\←W_D.*V_{\mathrm{lms}}{H}_h^T./W_D{H}_h{H}_h^T$

$H_h\←H_h.*W_D^T{V}_{\mathrm{lms}}./(W_D^T{W}_D{H}_h+\frac{\mathrm{\β}}{2}{H_g}^{-\frac{1}{2}});$

(3f)循环(3d)和(3e)，直到循环次数达到预先设定的最大循环次数n；

本发明实施例中，设β为0.03，最终得到的高空间分辨率的字典W_h大小为1024×256，具有光谱信息的系数矩阵H_h大小为256×324，低空间分辨率的字典W_D大小为64×256，无光谱信息的系数矩阵H_S大小为256×81。

步骤4.将高空间分辨率的字典W_h和具有光谱信息的系数矩阵H_h相乘得到分块列向量化的高空间分辨率的多光谱图像：V_hms＝W_hH_h。

本发明实施例中，分块列向量化的高空间分辨率的多光谱图像V_hms大小为1024×324。

步骤5.将分块列向量化的高空间分辨率的多光谱图像V_hms转化为图像块，恢复到图像的对应位置，得到高空间分辨率的多光谱图像。

本发明实施例中，分块列向量化的高空间分辨率的多光谱图像V_hms每一列被转化成一个32×32的小块，然后恢复到图像的对应位置，得到高空间分辨率的多光谱图像大小256×256×4，空间分辨率为2.4m。

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

1.仿真环境：

MATLAB R2013a,Intel(R)Core(TM)i5-3470CPU 3.00GHz,Window 7。

2.仿真内容：

仿真1，用本发明对Quickbrid卫星于2002年9月在印度孙德尔本斯地区拍摄的低空间分辨率的多光谱和高空间分辨率的全色图像进行融合，结果如图2所示，其中：

图2(a)是低空间分辨率的多光谱图像，大小为64×64×4，

图2(b)是高空间分辨率的全色图像，大小为256×256，

图2(c)是一幅参考图像，大小为256×256×4，

图2(d)为采用本发明对图2(a)和图2(b)进行融合后获得的高空间分辨率的多光谱图像，大小为256×256×4。

由图2可见，图2(d)的空间纹理信息相比于图2(a)有明显提高，尤其在道路与田地等细节部分，边缘更加清晰，与图2(b)相比，图2(d)的光谱信息更加丰富，由此可见，本发明能够对图2(a)和图2(b)进行较好的融合。

仿真2，为了证明本发明的效果，分别用本发明的方法和现有技术GIHS变换法、PCA变换法、AWLP变换法、压缩感知方法对图2(a)和图2(b)的待融合图像进行融合，并对融合结果进行客观指标评价，评价指标如下：

1)相关系数CC，表示光谱信息的保留程度，结果在区间[0,1]，相关系数越接近1，表示融合结果与参考图像越相似。

2)整体图像质量指数UIQI，表示融合图像与参考图像的接近程度，越接近1，融合结果越好。

3)光谱弧度SAM，表示光谱的扭曲程度，越接近0，融合结果越好。

4)全局质量评估Q4，表示图像在空间和光谱上的整体相似度，结果范围在区间[0,1]，全局质量评估指数越大，融合图像越与参考图像相似。

5)平均梯度AG，表示图像的清晰度，平均梯度值越大，图像越清晰，融合结果越好

根据上述评价指标对本发明与现有技术的融合结果从客观评价指标上进行评价，结果如表1。

表1.各类方法融合结果的客观评价

由表1可见，本发明的整体图像质量指数UIQI、全局质量评估Q4和平均梯度AG均大于现有技术的评价值，相关系数CC在红色、蓝色、近红外波段和平均值上均大于现有技术的评价值，相关系数CC在绿色波段上小于现有技术的评价值，虽然光谱弧度SAM小于PCA方法的评价值，但是本发明的光谱弧度SAM均大于其他现有技术，由此可以看出本发明的大部分客观评价指标均优于现有技术的客观评价指标。

Claims (4)

1.一种基于对偶稀疏非负矩阵分解的多光谱与全色图像融合方法，包括以下步骤：

(1)输入图像和矩阵：

(1a)分别输入低空间分辨率的多光谱图像和高空间分辨率的全色图像；

(1b)分别输入空间退化矩阵D和光谱退化矩阵S；

(2)分别将低空间分辨率的多光谱图像和高空间分辨率的全色图像分块列向量化，获得分块列向量化的低空间分辨率的多光谱图像V_lms和分块列向量化的高空间分辨率的全色图像V_pan；

(3)将分块列向量化的低空间分辨率的多光谱图像V_lms和分块列向量化的高空间分辨率的全色图像V_pan同时进行对偶稀疏非负矩阵分解，得到高空间分辨率的字典W_h和具有光谱信息的系数矩阵H_h：

(3a)构建目标方程：

$\left\{\begin{array}{c}\underset{W_D,H_h,W_h,H_s}{\min }\frac{1}{2}[{\left|\right|V_{\mathrm{lms}}-W_D{H}_h\left|\right|}_2^2+{\left|\right|V_{\mathrm{pan}}-W_h{H}_S\left|\right|}_2^2+\mathrm{\β}[{\left|\right|H_h\left|\right|}_{1/2}+{\left|\right|H_S\left|\right|}_{1/2}\left]\right]\\ \mathrm{subject\; to}{W}_D,H_h,W_h,H_S\≥0;W_D={\mathrm{DW}}_h;H_S=H_{h}S\end{array}\right.$

其中，min[·]为最小化，subjectto表示W_D,H_h,W_h,H_S≥0；W_D＝DW_h；H_S＝H_hS为 $\frac{1}{2}[{\left|\right|V_{\mathrm{lms}}-W_D{H}_h\left|\right|}_2^2+{\left|\right|V_{\mathrm{pan}}-W_h{H}_S\left|\right|}_2^2+\mathrm{\β}[{\left|\right|H_h\left|\right|}_{1/2}+{\left|\right|H_S\left|\right|}_{1/2}\left]\right]$ 的约束条件，W_D为低空间分辨率的字典，H_S为无光谱信息的系数矩阵，和为数据的保真项，||H_h||_1/2和||H_S||_1/2为稀疏项，分别表示具有光谱信息的系数矩阵H_h和无光谱信息的系数矩阵H_S的21范数，β为正则参数，权衡保真项与稀疏项之间的关系；

(3b)用0到1之间的任意数值对高空间分辨率的字典W_h和具有光谱信息的系数矩阵H_h的每一个元素进行初始化；

(3c)设置最大循环次数n＝1000；

(3d)固定W_D和H_h，利用下式更新W_h和H_S：

$W_h\←W_h\·*V_{\mathrm{pan}}{H}_S^T\·{/W}_h{H}_S{H}_S^T$

$H_S\←H_S\·*W_h^T{V}_{\mathrm{pan}}\·/(W_h^T{W}_h{H}_S+\frac{\mathrm{\β}}{2}{H}_S^{-\frac{1}{2}})$

其中，T为矩阵转置操作，.*为点乘操作，./为点除操作，←为结果更新操作；

(3e)固定W_h和H_S，利用下式更新W_D和H_h：

$W_D\←W_D\·*V_{\mathrm{lms}}{H}_h^T\·/W_D{H}_h{H}_h^T$

$H_h\←H_h\·*W_D^T{V}_{\mathrm{lms}}\·/(W_D^T{W}_D{H}_h+\frac{\mathrm{\β}}{2}{H}_h^{-\frac{1}{2}});$

(3f)循环(3d)和(3e)，直到循环次数达到预先设定的最大循环次数n；

(4)将高空间分辨率的字典W_h和具有光谱信息的系数矩阵H_h相乘得到分块列向量化的高空间分辨率的多光谱图像V_hms：

V_hms＝W_hH_h；

(5)将分块列向量化的高空间分辨率的多光谱图像V_hms转化为图像块，恢复到图像的对应位置，得到高空间分辨率的多光谱图像。

2.根据权利要求1所述的基于对偶稀疏非负矩阵分解的多光谱与全色图像融合方法，其特征在于：所述步骤(1b)中获取空间退化矩阵D，按如下公式进行：

$D=\frac{I_{m\×m}\⊗(e_{1\×4}\⊗(I_{m\×m}\⊗e_{1\×4}))}{16}$

其中，I_m×m为单位矩阵，e_1×4为元素全为1的行向量，为克罗内克积。

3.根据权利要求1所述的基于对偶稀疏非负矩阵分解的多光谱与全色图像融合方法，其特征在于：所述步骤(1b)中获取光谱退化矩阵S，按如下公式进行：

$S=w^T\⊗A_{N\×N}$

其中，A_N×N为基矩阵，该矩阵对角线元素为1，其余元素为0，N为低空间分辨率的多光谱图像每个颜色通道中小块的数目，w＝[w₁,w₂,w₃,w₄]为行向量，w₁，w₂，w₃，w₄分别为红，绿，蓝，近红外通道的权值系数。

4.根据权利要求1所述的基于对偶稀疏非负矩阵分解的多光谱与全色图像融合方法，其特征在于：所述步骤(2)中低空间分辨率的多光谱图像和高空间分辨率的全色图像分块列向量化，按如下步骤进行：

(2a)将低空间分辨率的多光谱图像中红颜色通道的图像分为m×m的小块，且每个小块之间重叠1个像素，m为小块的尺寸；

(2b)将每个小块的第2列串联到第1列，第3列串联到第2列，以此类推，完成块的列向量化，构成红色通道矩阵；

(2c)重复步骤(2a)和(2b)，构成低空间分辨率的多光谱图像中绿色，蓝色和近红外通道的矩阵；

(2d)将绿色通道矩阵连接到红色通道矩阵，将蓝色通道矩阵连接到绿色通道矩阵，以此类推，构成列向量化的低空间分辨率的多光谱图像V_lms；

(2e)将高空间分辨率的全色图像分为4m×4m的小块，且每个小块之间重叠4个像素；

(2f)将每个小块的第2列串联到第1列，第3列串联到第2列，以此类推，构成列向量化的高空间分辨率的全色图像V_pan。