基于稀疏表示和光谱自回归的多光谱和高光谱图像融合

Multispectral and Hyperspectral Image Fusion Based on Sparse Representation and Spectral Autoregression

摘要: 现有的高光谱图像融合算法采用将空间维度与光谱维度进行分开融合重建的方式，而高光谱图像光谱维度存在大量的光谱信息，这些像素信息比空间上的像素更加接近目标像素，对于空间维度的重建非常关键。因此，本文提出一种将光谱像素用于对空间维度进行空间信息修复的融合方法，称作基于稀疏表示和光谱自回归的多光谱和高光谱图像融合。该方法通过将光谱维度上的像素，通过自回归模型将其用于空间维度的信息修复，自回归模型由低分辨率的高光谱图像中学习得到，该模型自然保持了光谱间的关联性，通过在公开数据集上的验证了提出模型的有效性。

Abstract: Existing hyperspectral image fusion algorithms use separate fusion reconstruction of spatial and spectral dimensions, whereas there is a large amount of spectral information in the spectral dimension of hyperspectral images, and the information of these pixels is much closer to the target pixels than spatially similar pixels, which is very critical for the reconstruction of the spatial dimension. Therefore, in this paper, we propose a method to use spectral pixels for spatial information restoration on the spatial dimension, called multispectral and hyperspectral image fusion based on sparse representation and spectral autoregression. The method works by taking the pixels on the spectral dimension and using them for information restoration on the spatial dimension through an autoregressive model, the autoregressive model is learned from low-resolution hyperspectral images, which naturally maintains the correlation between the spectra, and the validity of the proposed model is verified by the validation on the publicly available dataset.

1. 引言

图像在现实生活中随处可见，其作为日常生活中不可或缺的一部分，扮演着传递信息、表达情感的重要角色。在实际的生活里，常见的数字图像是RGB图像，该图像记录了红、绿、蓝三通道的可见光谱波段的部分区域。在一些类似于图像识别任务中，RGB图像所蕴含的光谱信息太少，不足以精确地将物体进行区分，我们就需要采集更多不同波长的光线反射所形成的光谱信息来对目标物体进行甄别。于是，高光谱成像技术应运而生。高光谱图像具有几十甚至上千个光谱波段，对地物特征的捕捉更加的精细，在遥感任务 [1] 、目标分类 [2] 、目标检测 [3] 等领域有着广泛的应用。不幸地是，由于当前的高光谱成像设备的物理条件限制，想要同时获取一副在空间和光谱维度上均有较高分辨率的高光谱图像是十分困难的。这也就意味着在获取高光谱图像的时候，人们需要在光谱分辨率和空间分辨率之间有一个权衡，但是一幅在空间以及光谱上具有超高分辨率的高光谱图像却在生产活动中的应用具有重要意义。多光谱图像具有较高的空间分辨率，但其光谱分辨力较低。于是学者们自然地想到将一幅低分辨率的高光谱图像(Low-Resolution Hyperspectral Image, LR-HSI)和一幅高分辨率的多光谱图像(High-Resolution Multispectral Image, HR-MSI)进行融合来获得高分辨率的高光谱图像成为了该领域内的一种主流方法。

在上述背景下，大量学者做了很多相关工作。早期的高光谱图像融合算法是融合一张高分辨率的全色图像(Panchromatic Image, PAN)和一副低分辨率的多光谱图像 [4] [5] [6] 。随后人们将该融合方法扩展到高光谱图像融合。Chen等人 [7] 首次将低分辨的高光谱图像按照其光谱的分布范围将其分成几组，在每组数据里面，分别使用基于全色图像的方法进行高光谱图像融合。Selva等人 [8] 通过对高分辨率多光谱图像进行线性回归，合成高光谱图像的每个谱带的高分辨率图像，再通过基于全色图像的方法将低分辨率的高光谱图像的每个谱带与合成的高分辨图像进行融合。基于全色图像的高光谱图像融合方法，虽然具有较低的计算复杂度，可以被快速实现。但是当高分辨率的多光谱图像与低分辨率的高光谱图像之间存在较大的像素差异时，会造成比较严重的光谱扭曲。

随后，矩阵分解方法在图像处理中取得了较好的结果，于是开始涌现出了大量基于矩阵分解方法。基于矩阵分解的方法主要做法是将三维结构高光谱图像转化更容易处理的二维空间。该方法通常会将高光谱图像分解为光谱字典乘以空间系数的形式。其中，光谱字典包含了高光谱图像的光谱特征，而空间系数则包含了大量的空间信息。基于矩阵分解的高光谱图像融合方法根据其类型又可被分为稀疏表示方法 [9] [10] [11] 和低秩方法 [12] [13] [14] 。在稀疏表示模型中，光谱字典被视为一个过完备的字典，为了去获得其系数的稀疏性，通常来说字典所包含的原子数是要比光谱数要更大的，该方法假定高光谱图像的光谱信号是由几个光谱原子进行线性表示形成的。在稀疏表示模型中，光谱字典的学习通常来自于低分辨的高光谱图像，因为其包含了最多的光谱信息。Kawakami等人 [15] 首次提出了稀疏矩阵分解的方法来解决高光谱和多光谱图像的融合问题，他们首先通过稀疏字典学习算法从低分辨的高光谱图像中学习光谱字典，再通过使用高分辨率的多光谱图像进行稀疏编码算法来求解稀疏系数。Akhtar等人 [9] 充分利用了高分辨率的多光谱图形中的非局部自相似性，同时设计了一个耦合的贪心算法来求解稀疏系数。Dong等人 [11] 首次提出了一种过完备字典的学习算法，通过稀疏正则来确立稀疏系数，为了进一步利用空间信息存在的非局部自相似性，他们使用高分辨的多光谱图像来设计非局部自相似先验。而对于基于低秩的方法来说，因为光谱信号具有一定的冗余性，所以该方法假定光谱信号可以用低维的子空间中进行表示，子空间的维度是小于光谱维数的。Zhou等人 [16] 和Veganzones等人 [17] 使用了不同的局部低秩正则从高分辨的多光谱图像中探究局部信息，同时采用独立重建每一块局部区域的方法获得融合结果。

由于高光谱图像和多光谱图像本身就是三维的结构，所以可以自然地想到使用张量对高光谱图像以及多光谱图像进行表示。基于张量表示的方法通常聚焦于分解方式和构建的正则项。Dian等人 [18] 从低分辨率的高光谱图像中学习张量子空间，并通过按照非局部自相似性分组的张量块施加基于张量多秩的低秩正则来探索低秩性和非局部自相似性。同时，Dian等人 [19] 又提出基于张量列秩的分解，通过非局部自相似性对低分辨率的高光谱图像进行分解，形成四维的张量组，在其上施加低列秩正则，但由于张量列秩无法直接求解，于是他们找到近似的下届去逼近模型的最小值。Li等人 [20] 提出耦合稀疏张量分解，该方法通过使用近端交替优化方法去确立核张量以及三个维度上的因子矩阵，并对它们都施加稀疏正则来求解融合模型。最近，深度学习技术在各个研究领域都有着广泛应用。随着卷积神经网络的兴起，深度学习在图像处理领域大获成功，在图像超分辨率、去噪、识别等各种任务上表现出色，因此在多源图像融合领域上也涌现了大量基于深度学习的方法。Xie等人 [21] 通过使用卷积神经网络学习图像先验，然后将其加入到融合任务公式中组合成一个新的优化问题。Wang等人 [22] 设计了一个迭代增强的单元来充分利用两个待融合图像的信息，将其融入到融合网络中。

高光谱图像拥有大量的光谱信息，这些信息对于重建空间像素也具有指导意义。在光谱的可搜寻到的相似像素要多于在空间中搜索到的光谱像素。所以利用光谱像素可以更好地对空间上的进行重建。且利用光谱像素构建空间像素时候，可以在一定程度上保持住光谱的关联性，然而现存的方法往往会忽略这一属性。为了解决这一问题，本文设计了一种基于稀疏表示和光谱自回归高光谱图像融合方法。具体地说，为了可以充分构建空间上信息，本文设计了光谱自回归模型将其加入稀疏表示模型的正则项中，该模型通过在光谱维度上对空间像素进行重建，从另一维度对空间像素进行线性拟合。由于该拟合关系来源于低分辨率的高光谱图像，在一定程度上保证了光谱关联性。

2. 问题定义

2.1. 符号定义

本文用大写字母来表示矩阵，比如 $A\in {ℝ}^{n_1\times n_2}$ 。 $A^i\in {ℝ}^{1\times n_2}$ 表示矩阵 $A$ 第 $i$ 行。同时，本文使用花体字母来定义张量，例如 $\mathcal{A}$ 。对于一个三维的张量 $\mathcal{A}\in {ℝ}^{n_1\times n_2\times n_3}$ ，其在 $\left(i_1,i_2,i_3\right)$ 上的元素被定义为 $a_{i_1{i}_2{i}_3}$ 。 $\mathcal{A}\left(i_1,:,:\right)$ 、 $\mathcal{A}\left(:,i_1,:\right)$ 和 $\mathcal{A}\left(:,:,i_3\right)$ 分别定义的张量 $\mathcal{A}$ 的水平切面、侧切面和正面切面，其中我们特别指出将正面切面 $\mathcal{A}\left(:,:,i_3\right)$ 定义为 $A_{i_3}$ 。与此同时，本文定义 $\mathcal{A}\left(:,i_2,i_3\right)$ 、 $\mathcal{A}\left(i_1,:,i_3\right)$ 和 $\mathcal{A}\left(i_1,i_2,:\right)$ 分别表示张量 $\mathcal{A}$ 的一模纤维、二模纤维和三模纤维。张量的三模展开通过将三模纤维作为矩阵的列，定义为 $A_{(3)}={\text{unfold}}_{\text{3}}\left(\mathcal{A}\right)\in {ℝ}^{n_3\times n_1{n}_2}$ 。本文也需要使用到一些向量、矩阵以及张量的范数，下面给出这些范数的定义形式。对于一个向量 $v\in {ℝ}^{n_1}$ ，

其 $l_2$ 范数定义为 ${\Vert v\Vert }_2=\sqrt{{\displaystyle {\sum }_{i_1}{\left|v_{i_1}\right|}^2}}$ 。对于一个矩阵 $A$ ，其 $l_1$ 范数定义为 ${\Vert A\Vert }_1={\text{max}}_{i_2}{\displaystyle {\sum }_{i_1}\left|a_{i_1{i}_2}\right|}$ ，其F范数定义为 ${\Vert A\Vert }_F=\sqrt{{\displaystyle {\sum }_{i_1{i}_2}{\left|a_{i_1{i}_2}\right|}^2}}$ 。对于一个张量 $\mathcal{A}$ ，本文定义其F范数为 ${\Vert \mathcal{A}\Vert }_F=\sqrt{{\displaystyle {\sum }_{i_1{i}_2{i}_3}{\left|a_{i_1{i}_2{i}_3}\right|}^2}}$ 。

2.2. 问题定义

高光谱图像融合是对同一场景下的高分辨率的多光谱图像(HR-MSI)和低分辨率的高光光谱图像(LR-HSI)进行融合，如图1所示。本文定义高分辨的高光谱图像为 $\mathcal{Z}\in {ℝ}^{W\times H\times S}$ ，包含了 $W\times H$ 个空间像素和 $S$ 个光谱维度。 $\mathcal{X}\in {ℝ}^{w\times h\times S}$ 表示低分辨率的高光谱图像，它包含 $w\times h$ 个空间像素和 $S$ 个光谱维度。低分辨率的高光谱图像可以使用高分辨率的高光谱图像通过空间下采样的方式获：

$X_{(3)}=Z_{(3)}B$ (1)

这里的 $Z_{(3)}\in {ℝ}^{S\times WH}$ 和 $X_{(3)}\in {ℝ}^{S\times wh}$ 分别 $\mathcal{Z}$ 是 $\mathcal{X}$ 和三模展开。 $B\in {ℝ}^{WH\times wh}$ 表示点扩散函数和下采样矩阵的结合。

Figure 1. Schematic of hyperspectral image fusion

图1. 高光谱图像融合示意图

高分辨率的多光谱图像定义为，它包含 $W\times H$ 个空间像素点和 $s$ 个光谱维度。高分辨率的多光谱图像可以通过高分辨率的高光谱图像通过光谱下采样获得：

$Y_{(3)}=R{Z}_{(3)}$ (2)

这里的 $Y_{(3)}\in {ℝ}^{s\times WH}$ 定义的是 $\mathcal{Y}$ 的三模展开。 $R\in {ℝ}^{s\times S}$ 表示光谱响应函数。在多源图像融合中，通常需要解决一个如下的最优化问题：

$\underset{\mathcal{Z}}{\min }{\Vert X_{(3)}-Z_{(3)}B\Vert }_F^2+{\Vert Y_{(3)}-R{Z}_{(3)}\Vert }_F^2+\varphi \left(\mathcal{Z}\right)$ (3)

该最优化问题包含两个保真项，此处的定义的 $\varphi \left(\mathcal{Z}\right)$ 基于 $\mathcal{Z}$ 的某种先验，通常来说，多源图像融合是一个不可逆的反问题，需要设计合理的先验来对该最优化问题进行正则，才能获得较好的融合结果。

3. 提出模型

3.1. 稀疏表示模型

稀疏表示是针对不确定方程组提出的，对于不确定性方程组，在方程组系数满秩的条件下，方程组具有无数组解，稀疏表示通过将解的稀疏性作为方程组的约束，使得方程组具有唯一解。稀疏表示模型在信号处理等领域有着广泛的应用，它通过从字典中寻找最小原子的线性组合来表示信号。通过这种方式，即可将高维的数据进行线性分解，能够更好地压缩信号，抑制噪声。

在高光谱图像处理领域，稀疏表示模型已经通过大量的工作证明了其编码信号的有效性。该模型通过假定一个一幅清晰的高光图像可以通过以下公式进行获得：

$Z_{(3)}=DA$ (4)

这里 $D\in {ℝ}^{S\times d}$ 是光谱字典， $A\in {ℝ}^{d\times WH}$ 是其系数。在稀疏表示模型中通常假定该系数是稀疏的，该模型的可视化示意图如图2所示。但是在多源图像融合问题中，高分辨率的高光谱图像 $\mathcal{Z}$ 通常是未知的，所以我们无法直接从 $\mathcal{Z}$ 中直接求得光谱字典 $D$ 。幸运地是在多源图像融合问题中，同一场景下的低分辨率的高光谱图像 $\mathcal{X}$ 很好地保存了 $\mathcal{Z}$ 的光谱信息，因此我们可以从 $\mathcal{X}$ 中学习出光谱字典。

Figure 2. Schematic of sparse representation model

图2. 稀疏表示模型示意图

根据公式，低分辨率的高光谱图像 $\mathcal{X}$ 的每一像素也能够被光谱字典中的原子线性表示，即：

$X_{(3)}=D\alpha$ (5)

这里的 $\alpha$ 是低分辨率的高光谱图像对应的系数。该系数通常也被假定是稀疏的，因此字典的求解通过解下面这个优化问题：

$\left(D,\alpha \right)=\underset{D,\alpha }{\text{argmin}}\frac{1}{2}{\Vert X_{(3)}-D\alpha \Vert }_F^2+\lambda {\Vert \alpha \Vert }_0$ (6)

此处的 ${\Vert \alpha \Vert }_0$ 表示的是 $\alpha$ 的零范数，该范数统计的 $\alpha$ 中的非零个数，因此该范数通常用来表示目标的稀疏性。然而零范数并不是一个真正意义上的范数，因此其不满足范数的三角不等式性质，在实际应用中，其往往是难以优化的。相关研究学者已经证明在一定系数条件下，零范数可以转换成 $l_1$ 范数进行估计。有很多经典的字典学习算法已经被提出，本文将采用非负矩阵分解，因此公式(6)可以转换成如下形式：

$\begin{array}{c}\left(D,\alpha \right)=\underset{D,\alpha }{\text{argmin}}\frac{1}{2}{\Vert X_{(3)}-D\alpha \Vert }_F^2+\lambda {\Vert \alpha \Vert }_1\\ s.t.D\ge 0,\alpha \ge 0.\end{array}$ (7)

此处的 $\lambda$ 是一个正则项系数，公式(7)求解需要使用到非负字典学习算法，本文将使用文献 [11] 中的非负光谱字典学习算法来求解光谱字典D。我们将稀疏表示模型加入到公式(3)中，可以的得到如下公式：

$\begin{array}{c}{\min }_{\mathcal{Z},A}{\Vert X_{(3)}-DAB\Vert }_F^2+{\Vert Y_{(3)}-RDA\Vert }_F^2\\ +\varphi \left(\mathcal{Z}\right)+\eta {\Vert A\Vert }_1,\\ s.t.Z_{(3)}=DA\end{array}$ (8)

3.2. 光谱自回归模型

高光谱图像含有丰富的光谱信息，这些光谱维度上的信息要比空间上的搜索到局部或者是非局部信更加接近于真实的像素值，但是现存的很多方法更多的是在于将光谱和空间的信息独立处理，比如利用多光谱图像空间信息来重建空间信息，利用低分辨率的高光谱图像重建光谱信息。而本文提出的光谱自回归模型，假定光谱维度上的像素点可以由其前后n个范围内的点进行线性拟合，根据光谱局部关联性的特点，在本文中我们设定 $n=2$ 。这种交叉视角的融合方法不仅能对空间像素进行更加精确的拟合，同时由于像素的关系来源于LR-HSI，因此其在一定程度上保留的光谱间的关联性信息，其处理流程图如图3所示。

Figure 3. Spectral autoregressive model flowchart

图3. 光谱自回归模型流程图

高光谱图像的光谱信息主要保存在低分辨率的高光谱图像 $\mathcal{X}\in {ℝ}^{h\times w\times S}$ 中，因此我们将从LR-HSI中学习光谱自回归模型。首先需要将LR-HSI图像分割成T个光谱子条 ${\hat{x}}_1,\cdots ,{\hat{x}}_T$ ，其中 ${\hat{x}}_t\in {ℝ}^{2n+1}$ ，然后将它们使用k-means++算法进行光谱相似性聚类，将其聚成J类 $\left\{G_1,\cdots ,G_J\right\}$ ，其中第j个类 $G_j$ 的质心定义为 ${\beta }_j$ 。对于每一类 $G_j$ ，我们都可以学习一个光谱自回归模型。我们定义其第j类的光谱自回归模型为 $b_j$ ，它可以被通过如下的公式获得：

$b_j=\underset{\text{b}}{\text{argmin}}{\displaystyle \underset{\hat{x}\in G_j}{\sum }\left({\hat{x}}_0-b^{\top }{u}_{\hat{x}}\right)}$ (9)

这里的 ${\hat{x}}_0\in ℝ$ 是光谱子条上的中心像素， $u_{\hat{x}}$ 是中心点像素在光谱字条 ${\hat{x}}_t$ 前后n个范围内的像素点。我们可以获得 $J$ 个光谱自回归模型 $\left\{b_1,\cdots ,b_J\right\}$ 。在获得光谱自回归模型之后，我们就需要利用该模型完成图像融合。

然而，通常一幅高光谱图像具有广袤的信息，其空间范围巨大，因此逐个像素去进行处理，计算复杂度太大，不具备现实意义，因而需要找到一个合适的方式对空间像素进行分组并且在光谱上也要保证其相似性，这也就意味着，我们需要按照图像的色彩、亮度等因素对图像完成分组。在图像分割领域的超像素分割技术恰巧可以解决这一问题，本文使用线性迭代聚类算法(Simple Linear Iterative Clustering, SCIL) [23] 进行超像素分割。高分辨率的多光谱图像 $\mathcal{Y}\in {ℝ}^{H\times W\times s}$ 含有大量的空间信息，因此通过对其进行超像素分割，找到空间中色彩、亮度、纹理一样的 $K$ 个图像块 ${\bar{Y}}_1,\cdots ,{\bar{Y}}_K$ ，需要注意的是这些图像块并不是规则的。我们将多光谱图像的分组直接应用到 $\mathcal{Z}$ ，获得 $K$ 个不规则张量立方体 $\left\{{\hat{\mathcal{Z}}}_1,\cdots ,{\hat{\mathcal{Z}}}_K\right\}$ 。接着，需要为每个张量立方体 ${\hat{\mathcal{Z}}}_k$ 计算出其平均立方体，它可以通过下面公式获得：

${\alpha }_k=\frac{1}{N_k}{\displaystyle \underset{d_n\in \widehat{{\mathcal{Z}}_k}}{\sum }{d}_n}$ (10)

这里的 $d_n$ 是张量立方体 ${\hat{\mathcal{Z}}}_k$ 中的光谱条， $N_k$ 是其拥有的光谱条数量。根据公式(10)我们可以得到K个平均立方体 ${\alpha }_k$ 。对于需要获得高分辨率的高光谱图像 $\mathcal{Z}$ ，其空间中心像素 $z_0$ 需要和其在光谱维度上n范围内的像素向量 $u_z\in {ℝ}^{2n}$ 应该为 $\min {\Vert z_0-b_{jz}^{\top }{u}_z\Vert }_2^2$ ， $b_{jz}$ 通过 ${\text{argmin}}_j{\Vert {\alpha }_k-{\beta }_j\Vert }_2^2$ 获得。我们将光谱自回归模型作为正则项加入公式(8)得到如下形式：

$\begin{array}{c}{\min }_{\mathcal{Z},A}{\Vert X_{(3)}-DAB\Vert }_F^2+{\Vert Y_{(3)}-RDA\Vert }_F^2\\ +\lambda {\displaystyle \underset{k=1}{\overset{K}{\sum }}{\displaystyle \underset{z\in {\hat{\mathcal{Z}}}_k}{\sum }{\Vert z_0-b_{jz}^{\top }{u}_z\Vert }_2^2}}+\eta {\Vert A\Vert }_1,\\ s.t.Z_{(3)}=DA\end{array}$ (11)

这里的 $\lambda$ 是一个用来平衡光谱自回归模型的参数。为了简化计算，我们重写光谱自回归项，将其变为 ${\displaystyle \underset{i=1}{\overset{K}{\sum }}{\Vert \left(I-H_k\right){\widehat{Z_{(3)}}}^k\Vert }_2^2}$ ，这里的I是一个单位矩阵， ${\widehat{Z_{(3)}}}_k$ 是 ${\hat{\mathcal{Z}}}_k$ 的三模展开矩阵， $H_i\in {ℝ}^{(2\text{n+1)}\times N_k}$ 是光谱自回归矩阵，该矩阵的第i个元素的可以通过以下公式获得：

$H\left(l,m\right)=\{\begin{array}{l}{\beta }_l\text{if}{z}_m\in u_z,\beta \in b_{jz}\\ 0\text{otherwise}\end{array}$ (12)

与此同时，本文提出的模型变成：

$\begin{array}{c}{\min }_{\mathcal{Z},A}{\Vert X_{(3)}-DAB\Vert }_F^2+{\Vert Y_{(3)}-RDA\Vert }_F^2\\ +\lambda {\displaystyle \underset{k=1}{\overset{K}{\sum }}{\Vert \left(I-H_k\right){\widehat{Z_{(3)}}}^k\Vert }_2^2}+\eta {\Vert A\Vert }_1,\\ s.t.Z_{(3)}=DA\end{array}$ (13)

3.3. 优化

从上一小节提出的模型可以看出，该模型是一个非凸优化问题，因此我们需要使用交替方向乘子法(Alternating Direction Method of Multipliers, ADMM)迭代的方法来解决这个问题。我们通过引入变量 $Q=Z_{(3)}$ 以及 $E=A$ ，然后将公式(13)写成其增广拉个朗日的形式如下：

$\begin{array}{c}L\left(A,Z_{(3)},E,V_4\right)={\Vert X_{(3)}-Z_{(3)}B\Vert }_F^2+{\Vert Y_{(3)}-RDA\Vert }_F^2\\ +\lambda {\displaystyle \underset{k=1}{\overset{K}{\sum }}{\Vert \left(I-H_k\right){\hat{Q}}^k\Vert }_F^2}+\mu {\Vert E-A+\frac{V_2}{2\mu }\Vert }_F^2\\ +\mu {\Vert DE-Z_{(3)}+\frac{V_1}{2\mu }\Vert }_F^2+\mu {\Vert Q-Z_{(3)}+\frac{V_4}{2\mu }\Vert }_F^2\\ +{\eta }_1{\Vert E\Vert }_1\end{array}$ (14)

这里 $V_1\in {ℝ}^{S\times HW}$ 、 $V_2\in {ℝ}^{S\times d}$ 和 $V_4\in {ℝ}^{S\times HW}$ ，是拉格朗日乘子， $\mu$ 是惩罚项。ADMM通过固定一项求解另外一项的方法求解公式(14)，下面将介绍每个子问题的具体解法。

3.3.1. 求解Z₍₃₎

为了优化 $Z_{(3)}$ ，公式(14)被重写成关于 $Z_{(3)}$ 的子问题：

$\underset{Z_{(3)}}{\text{argmin}}{\Vert X_{(3)}-Z_{(3)}B\Vert }_F^2+\mu {\Vert DE-Z_{(3)}+\frac{V_1}{2\mu }\Vert }_F^2+\mu {\Vert Q-Z_{(3)}+\frac{V_4}{2\mu }\Vert }_F^2$ (15)

公式(15)有如下闭解：

$Z_{(3)}={\left[\left(B\right){\left(B\right)}^{\top }+2\mu I\right]}^{-1}\left[B{X}_{(3)}+\mu \left(DE+\frac{V_1}{2\mu }\right)+\mu \left(Q+\frac{V_4}{2\mu }\right)\right]$ (16)

3.3.2. 求解Q

求解 $Q$ ，公式(14)被重写成：

$\underset{Q}{\text{argmin}}\mu {\Vert Q-Z_{(3)}+\frac{V_4}{2\mu }\Vert }_F^2+\lambda {\displaystyle \underset{k=1}{\overset{K}{\sum }}{\Vert \left(I-H_k\right){\hat{Q}}^k\Vert }_F^2}$ (17)

由于 $K$ 个类别是按超像素分割进行分组的，所以每个 $k$ 之间是相互独立的。因此我们将采用子类之间独立优化的方法，将等式(17)重写成：

$\underset{Q}{\text{argmin}}{\displaystyle \underset{k=1}{\overset{K}{\sum }}\mu {\Vert {\hat{Q}}^k-{\widehat{Z_{(3)}}}^k+\frac{{\widehat{V_4}}^k}{2\mu }\Vert }_F^2+\lambda {\Vert \left(I-H_k\right){\hat{Q}}^k\Vert }_F^2}$ (18)

这里的 ${\widehat{V_4}}^k$ 、 ${\widehat{Z_{(3)}}}^k$ 分别是 $V_4$ 和 $Z_{(3)}$ 按照超像素分割后的结果。公式(18)可以通过下面的公式求解：

${\hat{Q}}^k={\left[\mu I+\lambda {\left(I-H_k\right)}^{\top }\left(I-H_k\right)\right]}^{-1}\mu \left({\widehat{Z_{(3)}}}^k-\frac{{\widehat{V_4}}^k}{2\mu }\right)$ (19)

我们将求解完成的 ${\hat{Q}}^k$ ，通过超像素分组位置进行拼接，还原成 $Q=\left\{{\hat{Q}}^1,\cdots ,{\hat{Q}}^K\right\}$ 。

3.3.3. 求解A

为求解 $A$ ，需要将公式(14)转化为：

$\underset{A}{\text{argmin}}{\Vert Y_{(3)}-RDA\Vert }_F^2+\mu {\Vert E-A+\frac{V_2}{2\mu }\Vert }_F^2$ (20)

问题(20)可通过如下公式获得闭解：

$A={\left[{\left(RD\right)}^{\top }\left(RD\right)+\mu I\right]}^{-1}\left[{\left(RD\right)}^{\top }{Y}_{(3)}+\mu \left(E+\frac{V_2}{2\mu }\right)\right].$ (21)

3.3.4. 求解E

对于参数E来说，我们优化E通过重写公式(14)为下面的式子：

$\underset{E}{\text{argmin}}\mu {\Vert E-A+\frac{V_2}{2\mu }\Vert }_F^2+{\eta }_1{\Vert E\Vert }_1$ (22)

公式(22)可以通过软阈值收缩算子进行解决，其解的形式为：

$E={\left[\text{Soft}\left(A-\frac{V_2}{2\mu },\frac{{\eta }_1}{2\mu }\right)\right]}_{+}.$ (23)

3.3.5. 求解拉格朗日乘子

拉格朗日乘子 $V_1$ 、 $V_2$ 、 $V_4$ 通过如下公式进行更新：

$\begin{array}{l}{V}_1=V_1+\mu \left(DA-Z_{(3)}\right)\\ V_2=V_2+\mu \left(E-A\right)\\ V_4=V_4+\mu \left(Q-Z_{(3)}\right)\end{array}$ (24)

以上就是本文提出模型的全部解法。

4. 实验

4.1. 实验设置

4.1.1. 数据集

为了验证提出方法的有效性，我们将对比方法和本文提出的方法在三个公开数据集上进行对比实验。第一个数据集是哥伦比亚大学计算机视觉实验室提供的CAVE数据集。该数据集有32张室内的高质量图片。在CAVE数据集中的每个数据集有512 × 512个空间像素和31个光谱维度，光谱波长间隔10 nm，范围从400 nm到700 nm。CAVE数据集中的高光谱图像就是原始数据，我们使用其对LR-HSI进行仿真，用平均模糊核通过平均大小为s × s图像对原始数据进行模糊，这里的s是下采样因子。HR-MSI的生成使用尼康D700相机的光谱响应，大小为512 × 512 × 3。

第二个数据集是帕维亚大学(Pavia University)数据集，Pavia University是通过使用rosis传感器在位于意大利的帕维亚大学收集的。该传感器收集了115个光谱带的数据，在移除一些异常和包含噪声的数据之后，剩下了103维光谱维的数据用于实验。该数据集的包含了610 × 340个空间像素。我们截取了其中320 × 320个空间像素，最终形成用于实验的数据集是320 × 320 × 103的真实图像。同样，LR-HIS的生成使用大小为s的平均模糊核进行生成。HR-MSI的光谱生成依靠伊科诺斯卫星(IKONOS)去生成一个320 × 320 × 4的图像。

4.1.2. 对比方法

本文是面向传统优化方法的高光谱图像和多光谱图像融合，所以和同类型的方法做对比，而非深度学习方法。我们选择目前在该领域内较为先进的方法：HySure [24] (A Convex Formulation for Hyperspectral Image Super-resolutionvia Subspace-Based Regularization)，NSSR [11] (Non-negative Structure Sparse Representation)，CSTF [20] (Coupled Sparse Tensor Factorization, CSTF)，LTTR [19] (Learning a Low Tensor-train rank Representation for Hyperspectral Image super-resolution)，IR-TenSR [25] (Iterative Regularization Method Based on Tensor Subspace Representation)。比较方法均使用各个方法中的最优参数。

4.1.3. 评价指标

多源图像融合任务，是需要对同一场景的地物来源于不同传感器的图像进行融合，以获得对地物更加精准的观测。因此需要对融合后的结果使用评价指标来衡量融合的算法模型性能，目前主要包括主观评价指标和客观评价指标。

主观评价方法主要是从视觉的角度主观地对融合后图像的纹理信息和边缘清晰度等进行评判。客观评价指标依靠数据通过某种方法对重建后的图像 $\mathcal{Z}$ 与对应的真实图像进行对比，目前常用的高光谱图像质量评价指标主要有以下几种：

均方根误差(Root Mean Squared Error, RMSE)：该指标通过融合后的图像和真实图像间的差值来衡量图像质量。其值越小，表示融合结果与真实图像越接近。该指标无单位，公式表示如下：

$\text{RMSE}=\sqrt{\frac{1}{WHS}{\displaystyle \underset{j=1}{\overset{S}{\sum }}{\displaystyle \underset{i=1}{\overset{WH}{\sum }}{\left(Z^{<i,j>}-Z_{\text{ori}}{}^{<i,j>}\right)}^2}}}$ (25)

峰值信噪比(PeekSignal-to-Noise Ratio, PSNR)：信号的最大功率和信号噪声功率之比，该指标反映融合后的图像与真实图像相比是否失真，通常以分贝(dB)来表示。PSNR指标越高，说明图像质量越好。

$\text{PSNR}=-10{\text{log}}_{10}\frac{1}{WH}{\displaystyle \underset{i=1}{\overset{W}{\sum }}{\displaystyle \underset{i=1}{\overset{H}{\sum }}{\left(Z^{<i,j>}-Z_{\text{ori}}{}^{<i,j>}\right)}^2}}$ (26)

结构相似性(Structure Similarity, SSIM)：该指标通过亮度、对比度和结构这三个因素来衡量两幅图像的相似度，其取值范围为[0,1]，SSIM的值越大，表示图像失真程度越小，图像的结构与真实图像的结构相似度越高，也反映出融合图像质量越好

(27)

其中，和 ${\mu }_{\mathcal{Z}}$ 、和 ${\sigma }_{\mathcal{Z}}$ 分别为图像信号和 $\mathcal{Z}$ 的平均值和标准差，为 $\mathcal{Y}$ 和 $\mathcal{Z}$ 的相关系数， $C_1$ 和 $C_2$ 为常数。

光谱角匹配(Spectral Angle Mapping, SAM)：评估每个像素点的光谱信息保存，通过计算重建光谱和地况光谱的两个向量之间的角度来确定图像的光谱失真，用于衡量图像光谱信息的完整程度。SAM数值越小，表示重建效果越好，该指标单位为度。公式如下：

$\text{SAM}\left({\mathcal{Z}}_{\text{ori}}\text{,}\mathcal{Z}\right)=\frac{1}{HW}{\displaystyle \underset{i=1}{\overset{HW}{\sum }}\arccos \left(\frac{\langle {\mathcal{Z}}_{\text{ori}}\text{,}\mathcal{Z}\rangle }{{\Vert {\mathcal{Z}}_{\text{ori}}\Vert }_2{\Vert \mathcal{Z}\Vert }_2}\right)}$ (28)

相对无量纲全局误差(Erreur Relative Globale Adimensionnelle De Synthsès, ERGAS)：可以衡量重建图像的全局质量，ERGAS数值越低，重建效果越好，该指标无单位。

$\text{ERGAS}\left({\mathcal{Z}}_{\text{ori}}\text{,}\mathcal{Z}\right)=100r\sqrt{\frac{1}{C}{\displaystyle \underset{i=1}{\overset{C}{\sum }}{\left(\frac{{\text{RMSE}}^i}{\text{Mean}\left({\mathcal{Z}}_{\text{ori}}^i\right)}\right)}^2}}$ (29)

其中， ${\text{RMSE}}^i$ 表示在第i通道上的均方误差， $\text{Mean}(\cdot )$ 表示图像的通道均值，r为超分辨率的比例因子。

通用图像质量指数(Universal Image Quality Index, UIQI)：该指标从相关损失、亮度畸变和对比度畸变这三个方面对融合结果和真实图像进行计算，数值越高，表示融合图像的效果越好，该指标无单位，公式如下：

$\text{UIQI}\left({\mathcal{Z}}_{\text{ori}}\text{,}\mathcal{Z}\right)=\frac{4{\mu }_{{\mathcal{Z}}_{\text{ori}}}{\mu }_{\mathcal{Z}}{\sigma }_{\mathcal{Z}{\mathcal{Z}}_{\text{ori}}}}{\left({\mu }_{{\mathcal{Z}}_{\text{ori}}}^2+{\mu }_{\mathcal{Z}}^2\right)\left({\sigma }_{{\mathcal{Z}}_{\text{ori}}}^2+{\sigma }_{\mathcal{Z}}^2\right)}$ (30)

4.2. 量化分析

4.2.1. 参数设置

对于提出的模型来说，有三个变量最为重要，它们分别是 $\lambda$ 、 $\eta$ 、 $K$ ，CAVE数据集具有较少的光谱维度，而Pavia University和Washington DC Mall都具有较高的数据集，因此我们选择对CAVE这类光谱较短的数据集，Pavia University及Washington DC Mall这一类光谱维度较长的数据集分开讨论。

Figure 4. Transformation curves of PSNR for each parameter

图4. 各参数的PSNR的变换曲线

图4中显示的是对于CAVE数据集和Pavia University数据集上的参数变化示意图。红色代表CAVE数据集的，蓝色代表Pavia University数据集。从图中的变化曲线可以看出，对于CAVE数据集来说，λ从 ${10}^{-6}$ 左右开始成上升趋势，直到 $0.5\times {10}^{-4}$ 附近达到最大值，之后开始呈现下降的趋势，因此在实验中，我们设置 $\lambda =0.5\times {10}^{-4}$ 。从关于η的变化曲线上可以看出，其峰值集中在 $0.1\times {10}^{-5}$ 到 $0.6\times {10}^{-5}$ ，所以在实验中我们设置 ${\eta }_1=0.6\times {10}^{-5}$ 。另一个较为关键的参数是K，通过观察其变化曲线可以发现，该曲线在200时候到达峰值，因此对于CAVE数据集我们设置 $K=200$ 。蓝色曲线是Pavia University数据集参数变化趋势，从图中我们可以看出，对于参数λ其变化趋势呈现先上升后下降的趋势，并在 $0.5\times {10}^{-3}$ 处达到峰值，因而对于Pavia University数据集来说，我们设置 $\lambda =0.5\times {10}^{-3}$ 。参数η其峰值位于 $0.8\times {10}^{-7}$ 到 $0.2\times {10}^{-6}$ 之间，因此在我们实验中，设置参数 $\eta =0.1\times {10}^{-6}$ 。K的峰值位于300左右，我们设置 $K=300$ 。Pavia University数据集各项参数是与Washington DC Mall数据集进行共享的，因为这两个数据集都属于光谱维度较长，所以我们使用共同的参数。

4.2.2. 对比实验

本节为了验证本文提出方法的有效性在三个数据集的实验结果。表1显示的是在CAVE数据集上不同方法在不同评价指标的实验结果。最好的实验结果已经加粗字体标出。

Table 1. Average experimental results for each comparison method on CAVE data

表1. 各个比较方法在CAVE数据集上的平均实验结果

从表1中够可以看出HySure方法在CAVE这种较高分辨率上的结果表现并不是很好，因为其探究的变分正则对于CAVE这类空间分辨率较高且其高频信息含量较少的数据集来说，并不是一个好的正则项。IR-TenSR方法主要面向的是对不同噪声的数据集有更加优秀的表现，而对于此类噪声单一，更加关注图像本身融合质量的任务来说，表现不够亮眼。相较于其他对比方法，NSSR取得了更加满意的结果，这依赖于其充分探究了图像的非局部相似性。但是空间上的像素的相似程度并没有光谱上的相似程度高，因此本文提出的方法在各个不同的下采样因子下面都获得了较好的实验结果。

表2中列出的是在Pavia University数据集上的结果。该光谱数据集相较于CAVE数据集，其拥有更多的光谱维度，但空间分辨率也更低一点，且包含较多的地物信息。因此HySure方法在该数据集上表现更加出色，在不同下采样因子下都取得了较好的实验结果。NSSR及LTTR在CAVE数据集上表现都较好，但对于Pavia University这类数据集自身含有一定噪声，且高频信息较多的，其设置的先验并不能够得到充分地学习，因此其效果产生下降。而反观本文提出的方法，在不同的下采样因子下都取得了最好的结果，这主要由于我们使用了光谱自回归利用光谱维度的信息对空间进行重建，且在一定程度上保存了光谱信息。

Table 2. Average experimental results for each comparison method on Pavia University data

表2. 各个比较方法在Pavia University数据集上的平均实验结果

4.3. 视觉效果展示

本章节将对提出的模型进行图像视觉上进行比较，我们通过制作了三个数据集的彩色图像及其误差图像。首先是CAVE数据集的视觉展示在图5，通过图中所展示的彩色图像分别是取了×16的结果图像10、15、25维组成的，误差图像为第15维与真实图像之间的差值。从视觉图像上可以看出来HySure、NSSR、IR-TenSR在插值图像上有较多的误差，CSTF在平滑区域有误差波纹，而本文提出的方法更加接近于真实图像这与我们的实验结果是一致的。

Figure 5. Full-color and error images of the CAVE dataset

图5. CAVE数据集的全彩图像和误差图像

而对于光谱维度更高的Pavia University数据集来说其视觉效果展示如图6所示，HySure和CSTF取的了较好的视觉效果，并且误差矩阵更小，LTTR以及IR-TenSR中心高频区域误差较大，且其全色图像在相同的光谱上形成的色彩效果也有差异，而相较于比较方法，本文提出的方法的无论是全彩图像还是误差图像都更加解决于真实图像。

Figure 6. Full-color and error images of the Pavia University dataset

图6. Pavia University数据集的全彩图像和误差图像

5. 总结与展望

本文提出了一种基于稀疏表示和光谱自回归的高光谱图像融合模型，该模型通过使用稀疏表示模型作为基本模型，在一定程度上抑制了空间噪声。同时，为了获得更好的融合结果，本文通过将光谱像素引入空间维度进行空间像素的重建，构建出光谱自回归模型对模型进行正则。光谱自回归模型不仅可以将光谱维度上相似的像素引入空间，对空间上的信息更加精确的重建，同时可以自然保留光谱关联性质。通过在两个公开数据集上的实验结果证明该方法的有效性。但该模型对空间维度信息利用不够充分，将在未来的研究中加以改进。