CN102842124A - 基于矩阵低秩分解的多光谱图像与全色图像融合方法 - Google Patents

Abstract

本发明公开了一种基于矩阵低秩分解的多光谱图像和全色图像融合方法，主要解决现有多光谱图像和全色图像融合中高光谱分辨率信息丢失的问题。其实现步骤为：（1）对输入已配准的4幅多光谱图像进行插值使其与全色图像具有相同的像素;(2)将插值后多光谱图像全部拉成列按照顺序依次堆叠构成大数据矩阵；（3）利用矩阵低秩分解算法对大数据矩阵分解得到低秩矩阵和稀疏矩阵；（4）利用标准主成分分析融合算法对低秩矩阵和全色图像进行融合得到初步的粗略融合结果；（5）将稀疏矩阵加至粗略融合结果得到最终的多光谱融合图像。本发明能有效解决高光谱分辨率信息丢失的问题，获取清晰的图像，可用于多光谱图像的预处理。

Description

基于矩阵低秩分解的多光谱图像与全色图像融合方法

技术领域

本发明属于图像处理技术领域，涉及图像融合技术，具体地说是一种结合了矩阵低秩分解理论的多光谱图像和全色图像融合方法，可用于多光谱图像的预处理。

背景技术

多光谱图像与全色图像的融合是利用它们在时空上的相关性及信息上的互补性把具有低空间分辨率、高光谱分辨率的多光谱图像和具有高空间分辨率、低光谱分辨率的全色图像进行综合，使融合后多光谱图像具有较高空间分辨率，又同时保留高光谱分辨率，从而得到对景物更全面、清晰的描述。

传统的多光谱图像与全色图像融合算法包括格拉姆-施密特正交化方法、亮度-色调-色饱和度方法等。这些方法计算量小，原理简单，已经被广泛地应用到多光谱图像与全色图像融合上，但是这些方法没有利用多光谱图像自身特点，无法有效保留高光谱分辨率信息。

针对上述传统的多光谱图像与全色图像融合算法效果较差，在实际应用中不能很好地实现的问题，目前国际上提出了一些改进上述缺点的算法。如，P.S.Chavez,J.S.Sides和J.A.Anderson提出一种主成分分析的方法，是多光谱图像与全色图像融合中非常经典的算法，具体参见《Comparison of Three Different Methods to MergeMultiresolution and Multispectral Data:Landsat TM and SPOT Panchromatic》PHOTOGRAMMETRIC ENGINEERING & REMOTE SENSING，Vol.57,No.3,March1991,pp.295-303。这种方法将全色图像直接替换掉多光谱图像的低空间分辨率信息从而得到融合的多光谱图像。这种方法速度快，比较灵活，但由于缺乏自适应性，而会造成高光谱分辨率信息的丢失；此后，T.M.Tu等人提出一种广义亮度-色调-色饱和度的方法，具体参见《A Fast Intensity-Hue–Saturation Fusion Technique With SpectralAdjustment for IKONOS Imagery》IEEE GEOSCIENCE AND REMOTE SENSINGLETTERS,VOL.1,NO.4,pp.309-312.OCTOBER2004.这种方法将传统只适用于三波段多光谱图像的亮度-色调-色饱和度方法推广到多波段的多光谱图像上，并采用固定的模型参数，由于模型参数被固定而造成光谱信息失真。

发明内容

本发明的目的在于针对上述已有技术的不足，提出了一种基于矩阵低秩分解的多光谱图像与全色图像融合方法，以保留较多的高光谱分辨率信息，提高融合图像的效果。

实现本发明目的的技术方案包括如下步骤：

(1)输入已配准的大小为m×n的4幅多光谱图像MS_l1，MS_l2，MS_l3，MS_l4，输入1幅大小为4m×4n全色图像PAN；

(2)利用matlab软件中的imresize函数对4幅多光谱图像MS_l1，MS_l2，MS_l3，MS_l4进行插值，得到与全色图像PAN相同像素的多光谱图像MS_l1′，MS_l2′，MS_l3′，MS_l4′；

(3)将插值后的4幅多光谱图像MS_l1′，MS_l2′，MS_l3′，MS_l4′全部拉成列，按顺序依次排列构成大数据矩阵X∈R^(4m×4n)×4，其中R^(4m×4n)×4表示行数为4m×4n，列数为4的2维整数型矩阵；

(4)通过低秩分解算法对大数据矩阵X进行低秩分解，得到低秩矩阵L∈R^(4m×4n)×4和稀疏矩阵SG∈R^(4m×4n)×4，其中低秩矩阵L表示4幅多光谱图像中的低空间分辨率信息，稀疏矩阵SG表示4幅多光谱图像中的高光谱分辨率信息；

(5)利用标准主成分分析融合算法对低秩矩阵L和全色图像PAN进行融合，得到4幅初步的粗略融合结果MS_h1′，MS_h2′，MS_h3′，MS_h4′;

(6)利用matlab软件中的reshape函数，将步骤(4)得到的稀疏矩阵SG的每一列还原成大小为4m×4n的图像，得到4幅高光谱分辨率图像sg_i，i＝1,2,...,4；

(7)将步骤(6)得到的4幅高光谱分辨率图像sg_i，i＝1,2,...,4,与步骤(5)得到的粗略融合结果MS_h1′，MS_h2′，MS_h3′，MS_h4′进行相加，得到最终的4幅融合多光谱图像MS_h1，MS_h2，MS_h3，MS_h4。

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

本发明由于在标准主成分分析融合方法的基础上引入了矩阵低秩分解，与传统的融合方法相比，充分利用了多光谱图像时空上的相关性及信息上的互补性，保留了更多的高光谱分辨率信息。仿真实验表明，本发明能有效保留多光谱图像中的高光谱分辨率信息，同时获得较好的视觉效果。

附图说明

图1是本发明的流程图；

图2是本发明在仿真试验中应用的4幅多光谱图像；

图3是本发明在仿真试验中应用的1幅全色图像；

图4是本发明在仿真试验中应用的4幅参考多光谱图像；

图5是用本发明对图2和图3进行融合的结果图；

图6是用现有的广义亮度-色调-色饱和度方法对图2和图3进行融合的结果图；

图7是用现有的标准主成分分析方法对图2和图3进行融合的结果图。

具体实施方式

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

步骤1、利用matlab软件中的imresize命令对输入已配准的4幅多光谱图像MS_l1，MS_l2，MS_l3，MS_l4进行插值，得到与全色图像PAN相同像素的多光谱图像MS_l1′，MS_l2′，MS_l3′，MS_l4′。

上述输入已配准的4幅多光谱图像MS_l1，MS_l2，MS_l3，MS_l4对应多光谱图像4个不同的波段，其中，图像MS_l1对应多光谱图像的红色波段，图像MS_l2对应多光谱图像的绿色波段,图像MS_l3对应多光谱图像的蓝色波段，图像MS_l4对应多光谱图像的近红外波段。

步骤2、将插值后多光谱图像MS_l1′，MS_l2′，MS_l3′，MS_l4′全部拉成列并按顺序依次排列构成大数据矩阵X∈R^(4m×4n)×4，其中R^(4m×4n)×4表示行数为4m×4n，列数为4的2维整数型矩阵。

步骤3、通过低秩分解算法对大数据矩阵X进行低秩分解，得到低秩矩阵L∈R^(4m×4n)×4和稀疏矩阵SG∈R^(4m×4n)×4。

上述的低秩分解方法，是由T.Y.Zhou和D.C.Tao提出的，参见文献《GoDec:Randomized Low-rank & Sparse Matrix Decomposition in Noisy Case》Proceedingsof the 28th International Conference on Machine Learning，2011，具体操作如下：

3a)初始化迭代次数t=0，迭代误差ε为0.0001；

3b)设t=t+1，利用matlab软件中的randn函数生成随机高斯矩阵A，按照如下公式得到3个中间变量矩阵：

P₁=X×A，

P₂=X^T×P₁，

P₃=X×P₂，

其中，(·)^T操作表示矩阵转置操作；

3c)计算第t次迭代中的低秩矩阵L_t和中间变量矩阵S_t：

L_t=P₃×(P₁ ^T×P₃)^-1P₂ ^T，

S_t=P_Ω|X-L_t|，

其中，(·)^-1表示矩阵求逆操作，P_Ω(·)表示取(·)中最大的前Ω个数值，Ω取值为262144；

3d)根据步骤3c)得到的低秩矩阵L_t和中间变量矩阵S_t计算分解误差

其中

表示矩阵2范数的平方；

3e)将分解误差ω与迭代误差ε进行比较，如果ω<ε，则停止迭代，并将矩阵L_t设为所求的低秩矩阵L，大数据矩阵X与低秩矩阵L的差设为所求的稀疏矩阵SG，否则返回步骤3b)，其中低秩矩阵L表示4幅多光谱图像中的低分辨空间信息，稀疏矩阵SG表示4幅多光谱图像中的高光谱分辨率信息。

步骤4、利用标准主成分分析融合算法对低秩矩阵L和全色图像PAN进行融合得到4幅初步的粗略融合结果MS_h1′，MS_h2′，MS_h3′，MS_h4′。

上述标准主成分分析融合算法，是由P.S.Chavez,J.S.Sides和J.A.Anderson提出的，参见文献《Comparison of Three Different Methods to Merge Multiresolution andMultispectral Data:Landsat TM and SPOT Panchromatic》PHOTOGRAMMETRICENGINEERING & REMOTE SENSING，Vol.57,No.3,March 1991,pp.295-303，具体操作如下：

4a)按照如下公式计算低秩矩阵L的相关矩阵C：

$C=\frac{L^T\&times;L}{4m\&times;4n\&times;4},$

4b)求相关矩阵C满足如下条件的特征值λ：

|C-λ×I|=0，

其中，|·|表示求矩阵的特征多项式，I是单位矩阵；

4c)求相关矩阵C满足如下条件的特征向量矩阵V：

(C-λ×I)×V=0，

其中，λ是步骤4b)中得到的特征值；

4d)按照如下公式计算中间变量矩阵PC和Y：

PC=L×V，

$Y=\frac{(\mathrm{PAN}-{\mathrm{\&mu;}}_{\mathrm{PAN}})\&times;D_{p{c}_4}}{D_{\mathrm{PAN}}}+{\mathrm{\&mu;}}_{{\mathrm{pc}}_4},$

其中，μ_PAN表示利用matlab软件的mean命令求出的全色图像PAN的平均值，表示利用matlab软件的mean命令求出的中间变量矩阵PC第4列pc₄的平均值，D_PAN表示利用matlab软件的std命令求出的全色图像PAN的标准方差，

表示利用matlab软件的std命令求出的中间变量矩阵PC第4列pc₄的标准方差；

4e)将中间变量矩阵PC的第4列pc₄用中间变量矩阵Y代替，得到新的中间变量矩阵PC'=[pc₁,pc₂,pc₃,Y]，并按照如下公式计算粗略融合结果矩阵MS′_h：

MS′_h＝PC'×V^T，

其中，(·)^T表示矩阵转置操作；

4f)利用matlab软件的reshape命令将粗略融合结果矩阵MS′_h变成4幅大小为4m×4n的粗略融合多光谱图像MS_h1′，MS_h2′，MS_h3′，MS_h4′。

步骤5、利用matlab软件中的reshape函数，将步骤(3)得到的稀疏矩阵SG的每一列还原成大小为4m×4n的图像，得到4幅高光谱分辨率图像sg_i，i＝1,2,...,4。

步骤6、将步骤5得到的4幅高光谱分辨率图像sg_i，i＝1,2,...,4，与步骤4得到的粗略融合结果MS_h1′，MS_h2′，MS_h3′，MS_h4′进行相加，得到最终的4幅融合多光谱图像MS_h1，MS_h2，MS_h3，MS_h4。其中，图像MS_h1对应融合多光谱图像的红色波段，图像MS_h2对应融合多光谱图像的绿色波段,图像MS_h3对应融合多光谱图像的蓝色波段，图像MS_h4对应融合多光谱图像的近红外波段。

本发明的效果可以通过仿真实验具体说明：

1.实验条件

实验所用微机CPU为Intel Core(TM)2Duo 3GHz内存2GB，编程平台是MatlabR2012a。实验采用4波段QuickBird卫星已配准多光谱图像和一幅全色图像和4幅参考多光谱图像，图像来源于图像融合网站http://www.glcf.umd.edu/data/。4幅已配准多光谱图像其大小均为64×64，如图2所示，每幅多光谱图像对应一个波段，其中图2(a)为多光谱图像的红色波段,图2(b)为多光谱图像的绿色波段,图2(c)为多光谱图像的蓝色波段,图2(d)为多光谱图像的近红外波段；1幅全色图像大小为256×256，如图3所示；图4为4幅参考多光谱图像其大小均为256×256，其中图4(a)为参考多光谱图像的红色波段,图4(b)为参考多光谱图像的绿色波段,图4(c)为参考多光谱图像的蓝色波段,图4(d)为参考多光谱图像的近红外波段。

2.实验内容与结果

实验一：利用本发明对图2和图3进行融合，结果如图5所示，其中图5(a)为融合多光谱图像的红色波段,图5(b)为融合多光谱图像的绿色波段,图5(c)为融合多光谱图像的蓝色波段,图5(d)为融合多光谱图像的近红外波段。

从图5可以看出，本发明得到的融合多光谱图像边缘清晰，具有较好的视觉效果。

实验二：利用现有的广义亮度-色调-色饱和度方法对图2和图3进行融合，结果如图6所示，其中图6(a)为融合多光谱图像的红色波段,图6(b)为融合多光谱图像的绿色波段,图6(c)为融合多光谱图像的蓝色波段,图6(d)为融合多光谱图像的近红外波段。

从图6可以看出，现有的基于广义亮度-色调-色饱和度方法得到的融合多光谱图像过于明亮，没有很好保持原有多光谱图像的高光谱分辨率信息。

实验三：利用现有的标准主成分分析方法对图2和图3进行融合，结果如图7所示，其中图7(a)为融合多光谱图像的红色波段,图7(b)为融合多光谱图像的绿色波段,图7(c)为融合多光谱图像的蓝色波段,图7(d)为融合多光谱图像的近红外波段。

从图7可以看出，现有的基于标准主成分分析方法得到的融合多光谱图像平滑，图像较模糊。

仿真实验中，应用相关系数CC和光谱角映射SAM指标来评价融合多光谱图像光谱保持度的优劣。其中，相关系数CC的定义如下：

$\mathrm{CC}=\frac{\underset{i=1}{\overset{M}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}[F(i,j)-\stackrel{\&OverBar;}{F}][X(i,j)-\stackrel{\&OverBar;}{X}]}{\sqrt{\underset{i=1}{\overset{M}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}{[F(i,j)-\stackrel{\&OverBar;}{F}]}^2\underset{i=1}{\overset{M}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}{[X(i,j)-\stackrel{\&OverBar;}{X}]}^2}},$

其中，X表示大小为M×N的参考多光谱图像，F表示融合后的大小为M×N的多光谱图像，

表示参考多光谱图像像素的平均值，

表示融合后的多光谱图像像素的平均值。

相关系数CC越接近1表明高光谱分辨率信息保持的越多。

光谱角映射SAM指标的定义如下：

$\mathrm{SAM}=\arccos \left(\frac{<u_X,u_F>}{{\left|\right|u_X\left|\right|}_2\&CenterDot;{\left|\right|u_F\left|\right|}_2}\right),$

其中，arccos(·)表示求反余弦操作，<·,·>表示向量的点乘操作，||·||₂表示矩阵2范数，u_X表示参考4波段多光谱图像X中对应像素点构成的谱向量，u_F表示融合4波段多光谱图像F中对应像素点构成的谱向量。

光谱角映射SAM指标越小表明高光谱分辨率信息保持的越多。

三种融合方法对实验所用多光谱图像和全色图像融合结果的定性评价指标如表1所示，其中，GIHS表示现有的广义亮度-色调-色饱和度融合方法，PCA表示现有的标准主成分分析融合方法。R表示参考多光谱图像与融合多光谱图像红色波段间的相关系数，G表示参考多光谱图像与融合多光谱图像绿色波段间的相关系数，B表示参考多光谱图像与融合多光谱图像蓝色波段间的相关系数，N表示参考多光谱图像与融合多光谱图像近红外波段间的相关系数，avg表示4个波段相关系数的平均值。

表1、多光谱图像和全色图像融合定性评价指标

从表1可以看出，本发明融合结果得到的相关系数CC更接近于1，光谱角映射SAM指标最小，因此，本发明能有效保持高光谱分辨率信息。

综上，本发明提出的基于矩阵低秩分解的多光谱图像和全色图像融合方法对多光谱图像和全色图像融合问题能有效保持高光谱分辨率信息并取得较好的视觉效果。

Claims (4)

1.一种基于矩阵低秩分解的多光谱图像和全色图像融合方法，包括如下步骤：

(1)输入已配准的大小为m×n的4幅多光谱图像MS_l1，MS_l2，MS_l3，MS_l4，输入1幅大小为4m×4n全色图像PAN；

(2)利用matlab软件中的imresize函数对4幅多光谱图像MS_l1，MS_l2，MS_l3，MS_l4进行插值，得到与全色图像PAN相同像素的多光谱图像MS_l1′，MS_l2′，MS_l3′，MS_l4′；

(3)将插值后的4幅多光谱图像MS_l1′，MS_l2′，MS_l3′，MS_l4′全部拉成列，按顺序依次排列构成大数据矩阵X∈R^(4m×4n)×4，其中R^(4m×4n)×4表示行数为4m×4n，列数为4的2维整数型矩阵；

(4)通过低秩分解算法对大数据矩阵X进行低秩分解，得到低秩矩阵L∈R^(4m×4n)×4和稀疏矩阵SG∈R^(4m×4n)×4，其中低秩矩阵L表示4幅多光谱图像中的低空间分辨率信息，稀疏矩阵SG表示4幅多光谱图像中的高光谱分辨率信息；

(5)利用标准主成分分析融合算法对低秩矩阵L和全色图像PAN进行融合，得到4幅初步的粗略融合结果MS_h1′，MS_h2′，MS_h3′，MS_h4′；

(6)利用matlab软件中的reshape函数，将步骤(4)得到的稀疏矩阵SG的每一列还原成大小为4m×4n的图像，得到4幅高光谱分辨率图像sg_i，i＝1,2,...,4；

(7)将步骤(6)得到的4幅高光谱分辨率图像sg_i，i＝1,2,...,4,与步骤(5)得到的粗略融合结果MS_h1′，MS_h2′，MS_h3′，MS_h4′进行相加，得到最终的4幅融合多光谱图像MS_h1，MS_h2，MS_h3，MS_h4。

2.根据权利要求1所述的基于矩阵低秩分解的多光谱图像和全色图像融合方法，其中所述步骤（1）中输入的已配准的大小为m×n的4幅多光谱图像MS_l1，MS_l2，MS_l3，MS_l4，对应不同的波段，即MS_l1对应多光谱图像的红色波段，MS_l2对应多光谱图像的绿色波段,MS_l3对应多光谱图像的蓝色波段，MS_l4对应多光谱图像的近红外波段。

3.根据权利要求1所述的基于矩阵低秩分解的多光谱图像和全色图像融合方法，其中步骤（4）所述通过低秩分解算法对大数据矩阵X进行低秩分解，得到低秩矩阵L和稀疏矩阵SG，具体步骤如下：

4a)初始化迭代次数t=0，迭代误差ε为0.0001；

4b)设t＝t+1，利用matlab软件中的randn函数生成随机高斯矩阵A，按照如下公式得到3个中间变量矩阵：

P₁=X×A，

P₂=X^T×P₁，

P₃=X×P₂，

其中，(·)^T操作表示矩阵转置操作；

4c)计算第t次迭代中的低秩矩阵L_t和中间变量矩阵S_t：

L_t=P₃×(P₁ ^T×P₃)^-1P₂ ^T，

S_t=P_Ω|X-L_t|，

其中，(·)^-1表示矩阵求逆操作，P_Ω(·)表示取(·)中最大的前Ω个数值，Ω取值为262144；

4d)根据步骤4c)得到的低秩矩阵L_t和中间变量矩阵S_t计算分解误差其中

表示矩阵2范数的平方；

4e)将分解误差ω与迭代误差ε进行比较，如果ω<ε，则停止迭代，并将矩阵L_t设为所求的低秩矩阵L，矩阵X与低秩矩阵L的差设为所求的稀疏矩阵SG，否则返回步骤4b)。

4.根据权利要求1所述的基于矩阵低秩分解的多光谱图像和全色图像融合方法，其中步骤（5）所述利用标准主成分分析融合算法对低秩矩阵L和全色图像PAN进行融合，具体步骤如下：

5a)按照如下公式计算低秩矩阵L的相关矩阵C：

$C=\frac{L^T\&times;L}{4m\&times;4n\&times;4},$

5b)求相关矩阵C满足如下条件的特征值λ：

|C-λ×I|=0，

其中，|·|表示求矩阵的特征多项式，I是单位矩阵；

5c)求相关矩阵C满足如下条件的特征向量矩阵V：

(C-λ×I)×V=0，

其中，λ是步骤5b)中得到的特征值；

5d)按照如下公式计算中间变量矩阵PC和Y：

PC＝L×V，

$Y=\frac{(\mathrm{PAN}-{\mathrm{\&mu;}}_{\mathrm{PAN}})\&times;D_{p{c}_4}}{D_{\mathrm{PAN}}}+{\mathrm{\&mu;}}_{{\mathrm{pc}}_4},$

其中，μ_PAN表示利用matlab软件的mean命令求出的全色图像PAN的平均值，

表示利用matlab软件的mean命令求出的中间变量矩阵PC第4列pc₄的平均值，D_PAN表示利用matlab软件的std命令求出的全色图像PAN的标准方差，表示利用matlab软件的std命令求出的中间变量矩阵PC第4列pc₄的标准方差；

5e)将中间变量矩阵PC的第4列pc₄用中间变量矩阵Y代替，得到新的中间变量矩阵PC'=[pc₁,pc₂,pc₃,Y]，按照如下公式计算粗略融合结果矩阵MS′_h：

MS′_h＝PC'×V^T，

其中，(·)^T表示矩阵转置操作；

5f)利用matlab软件的reshape命令将粗略融合结果矩阵MS′_h变成4幅大小为4m×4n的粗略融合多光谱图像MS_h1′，MS_h2′，MS_h3′，MS_h4′。

Images