CN108257107A - 基于局部低秩矩阵和光谱分解的图像融合算法 - Google Patents

Abstract

本发明公开了一种基于局部低秩矩阵和光谱分解的图像融合算法，该方法包括以下步骤：步骤(1)、局部化多光谱图像和高光谱图像，包括多光谱图像局部化处理和高光谱图像局部化处理；步骤(2)、进行多光谱丰度矩阵估计，求解丰度信息；步骤(3)、初始化高光谱图像局部块的端元矩阵E_p；步骤(4)、利用将多光谱图像的空间信息提取与高光谱图像的光谱信息提取两个过程通过交替迭代的方法计算E_m,p和X_h,p。与现有技术相比，本发明能够有效融合高光谱图像和多光谱图像，有着广泛的应用前景。

Description

基于局部低秩矩阵和光谱分解的图像融合算法

技术领域

本发明涉及图像处理领域，更具体地，涉及一种基于局部低秩矩阵和光谱分解的高 光谱与多光谱图像融合算法。

背景技术

根据光谱分辨率的大小，可以将光谱图像分为多光谱图像(MultispectralImage，SI) 和高光谱图像(Hyperspectral Image，HSI)。常见的多光谱图像包含具有7个波段的 Landsat-TM5，其波段范围分别是450-520nm、520-600nm、630-690nm、760-900nm、1550-1750nm、1040-1250nm2080-2350nm。常见的高光谱图像包含具有224个波段的AVIRIS，其光谱覆盖范围为400-2500nm。高光谱图像比多光谱图像具有更高的光谱分 辨率。然而这种光谱分辨率的提升带来的副作用就是其空间分辨率要低于多光谱图像。 由于高光谱图像较低的空间分辨率，其成像时往往造成大量的混合像元，极大的影响了 高光谱图像在各个领域中的应用。如何获得同时具有高的空间分辨率的高光谱图像 (Highspatial-resolutionHSI，HHSI)成为亟需解决的问题。

如果目标场景存在一幅高空间分辨率的多光谱数据，可以通过数据融合的方式来生 成高空间分辨率的高光谱图像。融合得到的图像可以看成是多光谱图像中的高空间分辨 率信息与高光谱图像中的高光谱分辨率信息的结合物。融合的结果能够使得以往只能在 机载高空间分辨率光谱成像系统上进行的应用扩展到全球尺度。例如：高空间分辨率生态系统监测，高空间分辨率专题制图(矿物质、城市地标物质、植物种类等)。尽管在 目前位置，同时具备高光谱图像与多光谱图像传感器的卫星很有限，但是随着可用的高 分辨率的光学卫星的增加和回归周期的改善，使得相同季节下和相同的大气状况和照明 条件下互补的高光谱与多光谱图像的获取成为可能。

迄今为止，在国内外公开发表的论文和文献中尚未见到基于局部低秩和耦合光谱分 解的高光谱和多光谱图像融合算法。因此，本发明具有独创性，该算法能够有效融合高光谱和多光谱图像，有着广泛的应用前景。

发明内容

基于对现有的高光谱图像和多光谱图像融合算法总结和分析的基础上，本发明提出 一种基于局部低秩矩阵和光谱分解的图像融合算法，利用高光谱图像局部低秩的性质，首先将高光谱图像与多光谱图像进行分块处理，然后，针对高光谱图像与多光谱图像之 间的每一对应块进行融合；对于每一个对应块，利用光谱图像线性分解模型来实现融合。

本发明的一种基于局部低秩矩阵和光谱分解的图像融合算法，该方法包括以下步骤：

步骤1、局部化多光谱图像和高光谱图像，包括多光谱图像局部化处理和高光谱图像局部化处理，具体如下：

多光谱图像局部化处理：将多光谱图像分割成ω×ω，固定大小的图像块，其 中ω是表示图像块的长度，表示多光谱图像所在的空间)；将多光谱图像所有的图像块构成的集合表示为其中，n表示块的数量，|Y_m,p|表示每一 个多光谱图像局部块中像素点的个数，λ_m表示多光谱图像的波段数；

高光谱图像局部化处理：将高光谱图像分割成(ω/d)×(ω/d)，固定大小的图像块；d表示空间下采样因子，即高光谱图像与多光谱图像在空间上的分辨率差异； 将高光谱图像所有的图像块构成的集合表示为其中|Y_h,p|表示每一个高光谱图像局部块中像素点的个数，λ_h表示高光谱图像的波段数，n表 示块的数量；

步骤2、进行多光谱丰度矩阵估计，对于每一个多光谱图像局部块，其对应的局部丰度矩阵定义为n_s表示光谱图像的子空间维度，|Y_m,p|表示每一个多光 谱图像局部块中像素点的个数；求解丰度信息，公式如下：

其中表示每一个局部多光谱图像块中的局部端元矩阵，p表示图像块 的索引，X_p≥0表示的是图像块中的每一个元素都符合非负的性质，Y_m,p表示每一个多光谱图像局部块中像素点的个数，表示Frobenius范数；

步骤3、初始化高光谱图像局部块的端元矩阵E_p，即将每一个高光谱图像局部块中的局部端元矩阵转换为下面的二次线性优化问 题的形式：

其中，X_h,p表示第p个高光谱图像块对应的丰度信息；

将正则项加到高光谱图像端元信息提取的二次线性优化式中得到端元提取的最优 化代价函数：

上式中的是数据拟合项，用来使估计出来的高光谱图像块端元矩阵 和块丰度矩阵能够准确的恢复出原始的高光谱图像块；μ·|Y_h,p|·V(E_p)是用来约束单纯形 体积的正则项；{μ·|Y_h,p}∈[0,1]是正则化参数，用来控制拟合项和正则项之间的相对重要 性；|Y_h,p|是使得正则参数会对不同的块大小做出不同的平衡的参数，μ是上式中唯一的 正则化参数；；

对代价函数进行求导，公式如下：

其中，n_s是子空间维数，tr(·)表示求矩阵的主对角线上元素之和，表示元素全为 1的矩阵；

令上述导数等于零，即得到最优的高光谱图像块端元矩阵：

步骤4、利用将多光谱图像的空间信息提取与高光谱图像的光谱信息提取两个过程 通过交替迭代的方法计算E_m,p和X_h,p，首先，利用VCA算法对高光谱图像块中的端元 矩阵E_p进行初始化，获得多光谱图像局部块的端元矩阵E_m,p和高光谱图像的丰度矩阵 X_h,p；

在之后的每一次迭代中，都通过对高光谱图像的端元矩阵E_p进行光谱下采样来获得E_m,p；然后在固定E_m,p的情况下，利用求解多光谱图像的 丰度矩阵X_p；利用该多光谱图像的丰度矩阵，通过空间模糊和下采样来得到高光谱图 像的X_h,p；然后在固定X_h,p的情况下，利用 求解多光谱图像块的端元矩阵 E_p；最后，通过E_p与X_p相乘就可以得到待重建的高空间分辨率的高光谱图像块Z_p：

Z_p＝E_pX_p

对多光谱图像空间信息提取和高光谱图像光谱信息提取这两个过程交替更新，直到 算法收敛；算法的收敛条件为：

其中，σ是一个接近零的常量；

输出HHSIZ＝{Z₁,Z₂,...,Z_n}，将得到的高光谱局部图像块Z_p按照分解时的方法重新 组合成一个高光谱图像。

与现有技术相比，本发明能够有效融合高光谱图像和多光谱图像，有着广泛的应用 前景。

附图标记

图1为本发明的一种基于局部低秩和光谱分解的光谱图像融合算法整体模型图；

图2为交替迭代优化模型示意图；

图3为利用MoffettField数据集在SAM算法下获得的融合过程的伪彩色图；(a)CNMF算法，(b)CLS(5*5)算法，(c)BSR算法，(d)本发明算法；

图4为利用Salinas数据集在SAM算法下获得的融合过程的伪彩色图；(a)CNMF 算法，(b)CLS(5*5)算法，(c)BSR算法，(d)本发明算法；

图5为利用Washington DC的SAM数据集在SAM算法下获得的融合过程的伪彩 色图；(a)CNMF算法，(b)CLS(5*5)算法，(c)BSR算法，(d)本发明算法。

具体实施方式

下面将结合附图对本发明的实施方式作进一步的详细描述。

考虑到光谱图像退化模型和光谱分解模型，基于光谱图像退化和光谱分解的高光谱 和多光谱图像融合模型(问题模型)可以表示为如下数学模型：

Y_m＝REX+N_m≈E_mX (1)

Y_h＝EXBM+N_h≈EX_h (2)

其中，Y_m表示多光谱图像模型，Y_h表示Y_m高光谱图像模型，表示与传感器光谱响应函数有关的光谱下采样矩阵。是一个点扩散函数，用来表示与波 段无关的传感器空间模糊运算。表示在横纵两个空间维度进行下采样系数为 d的均匀下采样矩阵。和表示高斯噪声，表示光谱退 化端元矩阵，表示空间退化丰度矩阵。

则有

X_h＝XBM (3)

E_m＝RE (4)

待重建的光谱图像Z是光谱信息的端元矩阵E和空间信息的丰度矩阵X的线性组合。又因为观察到的多光谱图像具有高的空间分辨率，观察到的高光谱图像具有高的光 谱分辨率。因此，重构高的空间和高光谱分辨率的光谱图像时采用两个线性优化模型分 别提取多光谱图像的丰度矩阵和高光谱图像的端元矩阵。

对于提取多光谱图像Y_m的丰度矩阵X，建立如下线性优化模型：

其中，||·||_F表示Frobenius范数。相似的，提取高光谱图像Y_h的端元矩阵E时，建立 光谱信息的端元矩阵E的线性优化模型如下：

因此，重构高的空间分辨率和高光谱分辨率的光谱图像Z转换为求解上述两个线性 优化模型。

实施例1

建立基于局部低秩矩阵和光谱分解的高光谱与多光谱图像融合模型；过程如下：

步骤1、局部化多光谱图像(MSI)和高光谱图像(HSI)，包括：

多光谱图像局部化处理：将多光谱图像分割成ω×ω固定大小的图像块，其 中ω是图像块的长度，是多光谱图像所在的空间，将多光谱图像所有的图像块构成的集合表示为其中n表示块的数量，|Y_m,p|表示每一个多光谱 图像局部块中像素点的个数，λ_m表示多光谱图像的波段数。

高光谱图像局部化处理：与多光谱图像局部化过程相同，将高光谱图像分割成(ω/d)×(ω/d)，固定大小的图像块；d表示空间下采样因子，即高光谱图像 与多光谱图像在空间上的分辨率差异；将高光谱图像所有的图像块构成的集合表示为其中|Y_h,p|表示每一个高光谱图像局部块中像素点的个数，λ_h表示高光谱图像的波段数，n表示块的数量。因为观察到的高光谱图像与多光谱 图像是经过配准的图像，因此高光谱图像中的每一个(ω/d)×(ω/d)的图像块都能在多 光谱图像中找到与之相对应的ω×ω的高光谱图像块；

步骤2、进行多光谱丰度矩阵估计，多光谱图像丰度矩阵估计问题可以转换为线性优化问题，即对于每一个多光谱图像局部块，其对应的局部丰度矩阵定义为 n_s表示光谱图像的子空间维度，|Y_m,p|表示每一个多光谱图像局部块中 像素点的个数。同时，因为丰度矩阵表示物质的光谱特性，因此X_p要保证非负。为了 求解丰度信息，将式(6)转换为约束最小二乘形式如下：

其中表示每一个局部多光谱图像块中的局部端元矩阵，p表示图像块 的索引，X_p≥0表示的是图像块中的每一个元素都符合非负的性质，Y_m,p表示每一个多光谱图像局部块中像素点的个数，表示Frobenius范数；

步骤3、通过顶点成分分析(Vertex ComponentAnalysis，VCA)算法初始化高光谱图像局部块的端元矩阵E_p

为了求解每一个高光谱图像局部块中的局部端元矩阵

问题转换为下面的二次线性优化问题的形式：

其中，X_h,p表示第p个高光谱图像块对应的丰度信息。

将正则项加到高光谱图像端元信息提取的二次线性优化式中得到端元提取的最优 化代价函数：

上式中的是数据拟合项，用来使估计出来的高光谱图像块端元矩阵 和块丰度矩阵能够准确的恢复出原始的高光谱图像块。μ·|Y_h,p|·V(E_p)是一项用来约束单 纯形体积的正则项。{μ·|Y_h,p|}∈[0,1]是正则化参数，用来控制拟合项和正则项之间的相对 重要性。|Y_h,p|的作用是使得正则参数会对不同的块大小做出不同的平衡，μ是上式中唯 一的正则化参数；。

上式的求解可以看成一个二次规划问题，为了求得最优的块端元矩阵，先对代价函 数进行求导：

其中，n_s是子空间维数，tr(·)表示求矩阵的主对角线上元素之和，表示元素全为 1的矩阵。

然后，令上述导数等于零，即得到最优的高光谱图像块端元矩阵：

步骤4、利用交替迭代更新的方式来计算E_m,p和X_h,p，提取多光谱图像块的X_p和高光谱图像块的E_p时，E_m,p和X_h,p是作为已知量的。因此，本专利采用图2所示的交替 迭代更新的方法来计算E_m,p和X_h,p，即：

首先，利用VCA算法对高光谱图像块中的端元矩阵E_p进行初始化。因为高光谱图像和多光谱图像拍摄的是同一场景，它们的潜在端元矩阵应该高度相关。因此，可以用 来获得多光谱图像局部块的端元矩阵E_m,p和高光谱图像的丰度矩阵X_h,p。

根据上面的原理，在之后的每一次循环中，都通过对高光谱图像的端元矩阵E_p进行光谱下采样来获得E_m,p。然后在固定E_m,p的情况下，利用式(7)求解多光谱图像的 丰度矩阵X_p。一旦获得了多光谱图像的丰度矩阵，通过空间模糊和下采样来得到高光 谱图像的X_h,p。然后在固定X_h,p的情况下，利用式(11)求解多光谱图像块的端元矩阵 E_p。这样就将多光谱图像的空间信息提取与高光谱图像的光谱信息提取两个过程通过 交替迭代的方法联系到了一起。最后，通过E_p与X_p相乘就可以得到待重建的高空间分 辨率的高光谱图像块Z_p：

Z_p＝E_pX_p (12)

对多光谱图像空间信息提取和高光谱图像光谱信息提取这两个过程交替更新，直到 算法收敛。算法的收敛条件为：

其中，σ是一个接近零的常量。

输出HHSIZ＝{Z₁,Z₂,...,Z_n}，将得到的高光谱局部图像块Z_p按照分解时的方法重新 组合成一个高光谱图像。

采用三幅高光谱图像来验证本发明的算法的有效性，通过这三幅图像合成出作为融 合算法输入的高光谱与多光谱图像。对于作为参考的高光谱图像来说，我们将其中对应水蒸气吸收波段和一些噪声波段进行去除。如表1所示，为本发明实验数据集。

表1

为了验证本文所提出算法的有效性，本文与现有的主流算法：CNMF，CLS和BSR 进行了对比。因为CLS算法同样需要局部块，为了方便对比，这里选取能在CLS下取 得较好效果的局部块大小来进行实验。对于每一个数据本文选取三种不同的局部块大 小，观察到的高光谱图像MoffettField、Salinas和WashingtonDC的局部块大小分别为 {5×5，9×9，15×15}，{5×5，10×10，25×25}和{20×20，30×30，60×60}。对于CNMF算法， 其内循环和外循环的最大迭代次数分别设为100和1。对于BSR算法，这里将参数a₀、 b₀、c₀、d₀和e₀设置为10^-6，Q＝35，迭代次数设为100。MoffettField、Salinas和Washington DC三个数据集在不同算法下获得的融合结果的伪彩色图分别显示在图3中(a)～(d)， 图4中(a)～(d)和图5中(a)～(d)。如表2所示，为客观指标对比。

表2

Claims (1)

1.一种基于局部低秩矩阵和光谱分解的图像融合算法，其特征在于，该方法包括以下步骤：

步骤(1)、局部化多光谱图像和高光谱图像，包括多光谱图像局部化处理和高光谱图像局部化处理，具体如下：

多光谱图像局部化处理：将多光谱图像分割成ω×ω固定大小的图像块，其中ω表示图像块的长度，表示多光谱图像所在的空间；将多光谱图像所有的图像块构成的集合表示为其中n表示块的数量，|Y_m,p|表示每一个多光谱图像局部块中像素点的个数，λ_m表示多光谱图像的波段数；

高光谱图像局部化处理：将高光谱图像分割成(ω/d)×(ω/d)，固定大小的图像块；d表示空间下采样因子，即高光谱图像与多光谱图像在空间上的分辨率差异；将高光谱图像所有的图像块构成的集合表示为其中|Y_h,p|表示每一个高光谱图像局部块中像素点的个数，λ_h表示高光谱图像的波段数，n表示块的数量；

步骤(2)、进行多光谱丰度矩阵估计，对于每一个多光谱图像局部块，其对应的局部丰度矩阵定义为n_s表示光谱图像的子空间维度，|Y_m,p|表示每一个多光谱图像局部块中像素点的个数)，求解丰度信息，公式如下：

其中表示每一个局部多光谱图像块中的局部端元矩阵，p表示图像块的索引，X_p≥0表示的是图像块中的每一个元素都符合非负的性质，Y_m,p表示每一个多光谱图像局部块中像素点的个数，表示Frobenius范数；

步骤(3)、初始化高光谱图像局部块的端元矩阵E_p，即将每一个高光谱图像局部块中的局部端元矩阵转换为下面的二次线性优化问题的形式：

其中，X_h,p表示第p个高光谱图像块对应的丰度信息；

将正则项加到高光谱图像端元信息提取的二次线性优化式中得到端元提取的最优化代价函数：

上式中的是数据拟合项，用来使估计出来的高光谱图像块端元矩阵和块丰度矩阵能够准确的恢复出原始的高光谱图像块；μ·|Y_h,p|·V(E_p)是用来约束单纯形体积的正则项；{μ·|Y_h,p|}∈[0,1]是正则化参数，用来控制拟合项和正则项之间的相对重要性；|Y_h,p|是使得正则参数会对不同的块大小做出不同的平衡的参数，μ是上式中唯一的正则化参数；

对代价函数进行求导，公式如下：

其中，n_s是子空间维数，tr(·)表示求矩阵的主对角线上元素之和，表示元素全为1的矩阵；

令上述导数等于零，即得到最优的高光谱图像块端元矩阵：

步骤(4)、利用将多光谱图像的空间信息提取与高光谱图像的光谱信息提取两个过程通过交替迭代的方法计算E_m,p和X_h,p，首先，利用VCA算法对高光谱图像块中的端元矩阵E_p进行初始化，获得多光谱图像局部块的端元矩阵E_m,p和高光谱图像的丰度矩阵X_h,p；

在之后的每一次迭代中，都通过对高光谱图像的端元矩阵E_p进行光谱下采样来获得E_m,p；然后在固定E_m,p的情况下，利用求解多光谱图像的丰度矩阵X_p；

利用该多光谱图像的丰度矩阵，通过空间模糊和下采样来得到高光谱图像的X_h,p；然后在固定X_h,p的情况下，利用求解多光谱图像块的端元矩阵E_p；最后，通过E_p与X_p相乘就可以得到待重建的高空间分辨率的高光谱图像块Z_p：

Z_p＝E_pX_p

对多光谱图像空间信息提取和高光谱图像光谱信息提取这两个过程交替更新，直到算法收敛；算法的收敛条件为：

其中，σ是一个接近零的常量；

输出HHSIZ＝{Z₁,Z₂,...,Z_n}，将得到的高光谱局部图像块Z_p按照分解时的方法重新组合成一个高光谱图像。