A Deep Unfolding Network for Multispectral and Hyperspectral Image Fusion

Abstract

Multispectral and hyperspectral image fusion (MS/HS fusion) aims to generate a high-resolution hyperspectral (HRHS) image by fusing a high-resolution multispectral (HRMS) and a low-resolution hyperspectral (LRHS) images. The deep unfolding-based MS/HS fusion method is a representative deep learning paradigm due to its excellent performance and sufficient interpretability. However, existing deep unfolding-based MS/HS fusion methods only rely on a fixed linear degradation model, which focuses on modeling the relationships between HRHS and HRMS, as well as HRHS and LRHS. In this paper, we break free from this observation model framework and propose a new observation model. Firstly, the proposed observation model is built based on the convolutional sparse coding (CSC) technique, and then a proximal gradient algorithm is designed to solve this model. Secondly, we unfold the iterative algorithm into a deep network, dubbed as MHF-CSCNet, where the proximal operators are learned using convolutional neural networks. Finally, all trainable parameters can be automatically learned end-to-end from the training pairs. Experimental evaluations conducted on various benchmark datasets demonstrate the superiority of our method both quantitatively and qualitatively compared to other state-of-the-art methods.

1. Introduction

Hyperspectral imagery (HSI) with high spatial resolution has the potential to greatly enhance our understanding and capabilities across a wide range of applications, including the monitoring and management of natural resources, ecosystems, biodiversity, and disasters [1]. However, due to the inherent limitations of imaging sensors, there exists a trade-off between spatial resolution and spectral resolution. Consequently, directly obtaining an HSI with high spatial resolution is a challenging task. Therefore, multispectral and hyperspectral image fusion (MS/HS fusion), i.e., generating a high-resolution hyperspectral image (HRHS) by combining a high-resolution multispectral image (HRMS) and a low-resolution hyperspectral image (LRHS), is an important task for remote sensing image processing. In recent decades, a large number of MS/HS fusion methods have been developed [1,2,3,4,5,6], and these approaches can be roughly classified into four categories, i.e., component substitution (CS) approaches [7], multiresolution analysis (MRA) methods [8,9], subspace-based methods [10,11,12,13,14,15,16,17], and deep learning (DL) approaches.

The CS-based MS/HS fusion method mainly originates from the CS-based pansharpening approaches [7,18,19], e.g., intensity–hue–saturation (IHS) [19], principal component analysis (PCA), and Gram–Schmidt adaptive (GSA) [7]. The general CS-based scheme sharpens the low-resolution images by multiplying the difference between a high-resolution image and a synthetic intensity component by a band-wise modulation coefficient. For example, in GSA, a representative CS-based pansharpening method obtains the intensity component by the Gram–Schmidt transform [18]. A direct approach to using GSA for MS/HS fusion is to create multiple sets of images for pansharpening subproblems. Each set consists of an MS band and its associated HS bands, which are then grouped together based on correlation-based clustering. The CS-based MS/HS fusion method can enhance the spatial details well, but spectral distortion may occur, especially when there is a significant difference in the spectral range between the MS and HS images [1,7].

The MRA-based approach differs from the CS-based approach in how spatial details are obtained from high-resolution images. In the conventional pansharpening procedure, the MRA method [8,9] sharpens the entire MSI by individually sharpening the characteristics of each MS band. Specifically, it supplements the spatial details of each interpolated MS band, which is then upsampled to the same dimensions as the Pan image. Two representative methods are smoothing filtered-based intensity modulation (SFIM) [8] and generalized Laplacian pyramid (GLP) [9]. When applying such methods to MS/HS fusion task, Selva et al. proposed a novel framework called hyper-sharpening [20], where the source of high spatial resolution details is no longer a specific MS band, but rather a high-resolution image synthesized using linear regression with the MS images as input. The SFIM hyper-sharpening method (SFIM-HS) is an extension of the traditional SFIM [8] method, where a synthesized image is generated by applying a linear regression of MS bands using the least square. Similar to SFIM, GLP [9] has also been applied to MS/HS fusion effectively, i.e., GLP-HS. Compared to CS-based methods, MRA-based approaches exhibit superior spectral quality but worse spatial sharpening effects.

The subspace-based methods can be subdivided mainly into two categories, i.e., unmixing [11,12,13,16,17,21,22] and Bayesian-based approaches [10,14,15]. Most subspace-based methods involve decomposing the HS and MS images to obtain spectral signals of an endmember and their corresponding fractional abundances. Constrained Nonnegative Matrix Factorization (CNMF) [11] alternately unmixes the HS and MS images by NMF [23] to estimate the spectral signatures of endmembers and the high-resolution abundance maps, respectively. Similarly to CNMF, the Lanaras algorithm [13] jointly unmixes the two input images into the spectral signatures of the endmembers and the associated fractional abundances. Akhtar et al. [12] employed dictionary learning and sparse coding techniques to obtain the endmembers and high-resolution abundances, respectively. To preserve edges and smooth out noise in homogeneous regions, HySure [10] introduced total variation regularization into the subspace-based HSI super-resolution framework. The CSTF method [16] formulated the MS/HS fusion task as a tensor factorization problem. Bayesian-based methods first model the MS/HS fusion problem as a probabilistic optimization problem, and then solve it using an iterative algorithm. The MAP-SMM [14] algorithm uses the stochastic mixing model (SMM) to estimate the underlying spectral scene statistics, as well as the local statistics of the spectral mixing model. FUSE [15] utilizes a Sylvester equation to solve the maximization problem of the likelihoods obtained from the forward observation models. However, the subspace-based methods, i.e., unmixing or Bayesian-based, always require assuming subjective prior information and high computational complexity.

In recent years, deep learning-based approaches [24,25,26,27,28] have been extensively studied for MS/HS fusion. This type of method directly learns the mapping from HRMS and LRHS to HRHS using some network architectures. For example, Han et al. [29] presented an MS/HS fusion method using cluster-based multi-branch BP neural networks. Dian et al. [30] proposed an MS/HS fusion method based on the subspace representation and convolutional neural network (CNN) denoiser. Dian et al. [31] designed an HSRNet by restoring spatial and spectral details using a spectral preservation module and a spatial preservation module, respectively. Palsson et al. [32] introduced the deep residual network for the MS/HS fusion task. There are also some other methods to address the challenge that multimodal remote sensing images cannot be directly compared due to their modal heterogeneity. For example, Chen et al. [33] take advantage of two types of modality-independent structural relationships in multi-modal images. In particular, they present a structural relationship graph representation learning framework for measuring the similarity of the two structural relationships. Additionally, they also designed an adaptive fusion method [34] based on frequency-decoupling to effectively fuse the local and nonlocal structural difference maps. However, these deep learning-based methods are black-box and lack interpretation. To deal with the limitations of the black-box networks, researchers proposed various model-based deep learning methods by unfolding the solution of the variational model into CNN. Xie et al. [35] first proposed a deep unfolding MS/HS fusion network based on the linear observation model, i.e., there are linear relationships between HRHS and HRMS, as well as between HRMS and LRHS. Dong et al. [36] proposed an iterative Hyperspectral Image Super-Resolution (HSISR) algorithm based on a deep HSI denoiser to leverage both domain knowledge likelihood and deep image prior and then unfold the iterative HSISR algorithm into a novel model-guided deep convolutional network. Ma et al. [37] proposed a model-based deep learning network named unfolding spatiospectral super-resolution network (US3RN) to reformulate the image degradation and incorporate the spatiospectral super-resolution (SSSR) model, which takes the observation models of SISR and SSR into consideration and then solve the model-based energy function via the alternative direction multiplier method (ADMM) technique. Sun et al. [38] propose a new image fusion network that uses deep prior information as spatial guide information. The fusion model uses spatial information as a regularization term and is solved using half splitting quadratic method. These deep unfolding methods can obtain excellent performance and sufficient interpretability by relying on a fixed linear observation model. However, the fixed linear model only considers the interaction between HRHS, HRMS, and LRHS, while the characteristics of HRHS, HRMS, and LRHS are ignored. In this paper, we continue our research along the deep unfolding method and attempt to jump out of this observation model framework by considering the image features. Specifically, we first propose a new observation model based on the convolutional sparse coding (CSC) [39,40] technique and its corresponding proximal gradient algorithm. Secondly, we unfold the iterative algorithm into a deep network (MHF-CSCNet), where the proximal operators are learned using convolutional neural networks. Finally, all trainable parameters can be learned end-to-end from the training pairs.

In summary, the contributions of our work are three-fold:

• We propose a new observation model for the MS/HS fusion task using the convolutional sparse coding (CSC) technique. Specifically, we separately model the HS and MS images using CSC, incorporating two types of features: common features and unique features. These features are derived from two key observations. First, the HS and MS images capture the same scene, indicating the presence of common features. Second, the images are generated by distinct sensors on the satellite, implying specific unique features specific to each image.
• We reformulate the observation model as an optimization problem and then place implicit prior on both the common and unique features. This allows us to leverage additional information and improve the quality of the fused image. Then, we develop a proximal gradient algorithm to solve the optimization problem efficiently, which is then unfolded into a deep network. Each module of the network corresponds to a specific operation in the iterative algorithm and, thus, the interpretability of the network can be significantly improved.
• Experimental results on benchmark datasets demonstrate the superiority of our network both quantitatively and qualitatively compared with other state-of-the-art methods.

2. Proposed Method

In this section, we introduce the proposed methodology. Initially, we formalize the established model as an optimization problem and decompose it into three distinct subproblems, utilizing the proximal gradient algorithm to derive the model’s solution. Next, we describe the solution steps by presenting them within a network framework and provide a clear overview of the network’s structure. Finally, we explain the basic settings used during the training of the network.

2.1. Problem Formulation

The MS/HS fusion model is proposed on the basis of the following observations. First, since both the LRHS image $H\in {\mathbb{R}}^{m\times n\times B}$ and the HRMS image $M\in {\mathbb{R}}^{M\times N\times b}$ are obtained from the same scenario, they should have some common features. Second, since the two images are acquired by different sensors in the satellite, thus each of them has its own unique features. Furthermore, we also assume that the information underlying both the LRHS image and the HRMS image is beneficial to estimate the HRHS image $O\in {\mathbb{R}}^{M\times N\times B}$. Therefore, we first extract these effective features from both images and then fuse these features to generate the HRHS image O. To acquire the features from both images, we adopt the CSC technique to model the LRHS and HRMS images separately.

Specifically, our new model is defined as follows:

$$M=\sum _{k=1}^K{D}_k^c\otimes C_k+\sum _{k=1}^K{D}_k^u\otimes U_k,$$

(1)

$$\tilde{H}=\sum _{k=1}^K{P}_k^c\otimes C_k+\sum _{k=1}^K{P}_k^v\otimes V_k,$$

(2)

where $\tilde{H}\in {\mathbb{R}}^{M\times N\times B}$ is the upsampled version of LRHS image H using the polynomial kernel technique; ⊗ is the convolutional operation; $D_k^c\in {\mathbb{R}}^{s\times s\times b}$ and $P_k^c\in {\mathbb{R}}^{s\times s\times B}$ denote the common filter for both images; $D_k^u\in {\mathbb{R}}^{s\times s\times b}$ and $P_k^v\in {\mathbb{R}}^{s\times s\times B}$ represent the filters of both images, respectively; $C_k\in {\mathbb{R}}^{M\times N}$ denotes the shared common features; and $U_k\in {\mathbb{R}}^{M\times N}$ and $V_k\in {\mathbb{R}}^{M\times N}$ are unique features of HRMS image and LRHS image, respectively. Then, the HRHS image $O\in {\mathbb{R}}^{M\times N\times B}$ can be generated by fusing these features as follows:

$$O=\sum _{k=1}^K{G}_k^c\otimes C_k+\sum _{k=1}^K{G}_k^u\otimes U_k+\sum _{k=1}^K{G}_k^v\otimes V_k,$$

(3)

where $G_k^c\in {\mathbb{R}}^{s\times s\times B}$ represents the common filter, $G_k^u\in {\mathbb{R}}^{s\times s\times B}$ and $G_k^v\in {\mathbb{R}}^{s\times s\times B}$ denote the unique filters concerning the features of both images. To make the notations clear, we denote $C={\left\{C_k\right\}}_{k=1}^K\in {\mathbb{R}}^{M\times N\times K}$, $U={\left\{U_k\right\}}_{k=1}^K\in {\mathbb{R}}^{M\times N\times K}$, and $V={\left\{V_k\right\}}_{k=1}^K\in {\mathbb{R}}^{M\times N\times K}$.

In the following, we assume that all the filters are known. Therefore, to obtain the fused HRHS image O, we need to calculate C, U, and V as follows:

$$\begin{array}{cc}\hfill mi{n}_{U,V,C}& {\displaystyle \frac{1}{2}}\Vert M-\sum _{k=1}^K{D}_k^c\otimes C_k-\sum _{k=1}^K{D}_k^u\otimes U_k{\Vert }_F^2\hfill \\ & +{\displaystyle \frac{1}{2}}\Vert \tilde{H}-\sum _{k=1}^K{P}_k^c\otimes C_k-\sum _{k=1}^K{P}_k^v\otimes V_k{\Vert }_F^2\hfill \\ & +{\lambda }_1{f}_1\left(U\right)+{\lambda }_2{f}_2\left(V\right)+{\lambda }_3{f}_3\left(C\right),\hfill \end{array}$$

(4)

where ${\lambda }_1,{\lambda }_2$, and ${\lambda }_3$ are the trade-off parameters. $f_1(·)$, $f_2(·)$, and $f_3(·)$ are the regularization terms modeling the priors on $U,V$, and C, respectively.

2.2. Model Optimization

To solve the problem, we adopt the alternate optimization algorithm to alternately update each variable with other variables fixed. First, the problem is broken down into three sub-problems:

$$U^{\left(t\right)}=argmi{n}_U\Vert M-\sum _{k=1}^K{D}_k^c\otimes C_k^{(t-1)}-\sum _{k=1}^K{D}_k^u\otimes U_k{\Vert }_F^2+{\lambda }_1{f}_1\left(U\right),$$

(5)

$$V^{\left(t\right)}=argmi{n}_V\Vert \tilde{H}-\sum _{k=1}^K{P}_k^c\otimes C_k^{(t-1)}-\sum _{k=1}^K{P}_k^v\otimes V_k{\Vert }_F^2+{\lambda }_2{f}_2\left(V\right),$$

(6)

$$\begin{array}{cc}\hfill C^{\left(t\right)}=argmi{n}_C& \Vert M-\sum _{k=1}^K{D}_k^c\otimes C_k-\sum _{k=1}^K{D}_k^u\otimes U_k^{\left(t\right)}{\Vert }_F^2\hfill \\ & +\Vert \tilde{H}-\sum _{k=1}^K{P}_k^c\otimes C_k-\sum _{k=1}^K{P}_k^v\otimes V_k^{\left(t\right)}{\Vert }_F^2+{\lambda }_3{f}_3\left(C\right).\hfill \end{array}$$

(7)

Next, we introduce the optimization for each sub-problem.

2.2.1. U-Subproblem

Instead of directly optimizing the objective (5), we first adopt the quadratic function to approximate the objective function, which is expressed as:

$$\begin{array}{cc}\hfill U^{\left(t\right)}=argmi{n}_U& g\left(U^{(t-1)}\right)+{\displaystyle \frac{1}{2{\eta }_1}}{\Vert U-U^{(t-1)}\Vert }_F^2\hfill \\ & +\langle \nabla g\left(U^{(t-1)}\right),U-U^{(t-1)}\rangle +{\lambda }_1{f}_1\left(U\right),\hfill \end{array}$$

(8)

where

$$g\left(U^{(t-1)}\right)=\Vert M-\sum _{k=1}^K{D}_k^c\otimes C_k^{(t-1)}-\sum _{k=1}^K{D}_k^u\otimes U_k^{(t-1)}{\Vert }_F^2,$$

(9)

$$\nabla g\left(U^{(t-1)}\right)=D^u{\otimes }^{†}(\sum _{k=1}^K{D}_k^c\otimes C_k^{(t-1)}+\sum _{k=1}^K{D}_k^u\otimes U_k^{(t-1)}-M).$$

(10)

∇ is gradient operator, ${\eta }_1$ is the step size, $D^u\in {\mathbb{R}}^{s\times s\times K\times b}$ is a 4-D tensor stacked by all $D_k^u$, and ${\otimes }^{†}$ denotes the transposed convolution. By a simple mathematical derivation, Equation (8) can be simplified to

$$U^{\left(t\right)}=argmi{n}_U{\displaystyle \frac{1}{2}}\left|\right|U-(U^{(t-1)}-{\eta }_1\nabla g\left(U^{(t-1)}\right)){\left|\right|}_F^2+{\eta }_1{\lambda }_1{f}_1\left(U\right).$$

(11)

Equation (11) can be optimized by the general proximal operator, which is defined as

$$U^{\left(t\right)}=pro{x}_{{\eta }_1{\lambda }_1}(U^{(t-1)}-{\eta }_1\nabla g\left(U^{(t-1)}\right)),$$

(12)

where $pro{x}_{{\eta }_1{\lambda }_1}(·)$ is the proximal operator related to $f_1(·)$. To overcome the constraint of human-designed regularizers, we employ deep CNNs to autonomously acquire prior knowledge from the data.

2.2.2. V-Subproblem

To update V, we employ a similar approach. Firstly, we define the quadratic approximation of Equation (6) with respect to V as follows:

$$\begin{array}{cc}\hfill V^{\left(t\right)}=argmi{n}_V& h\left(V^{(t-1)}\right)+{\displaystyle \frac{1}{2{\eta }_2}}{\Vert V-V^{(t-1)}\Vert }_F^2\hfill \\ & +\langle \nabla h\left(V^{(t-1)}\right),V-V^{(t-1)}\rangle +{\lambda }_2{f}_2\left(V\right),\hfill \end{array}$$

(13)

where

$$h\left(V^{(t-1)}\right)=\Vert \tilde{H}-\sum _{k=1}^K{P}_k^c\otimes C_k^{(t-1)}-\sum _{k=1}^K{P}_k^v\otimes V_k^{(t-1)}{\Vert }_F^2,$$

(14)

$$\nabla h\left(V^{(t-1)}\right)=P^v{\otimes }^{†}(\sum _{k=1}^K{P}_k^c\otimes C_k^{(t-1)}+\sum _{k=1}^K{P}_k^v\otimes V_k^{(t-1)}-\tilde{H}).$$

(15)

∇ is gradient operator, ${\eta }_2$ is the step size, $P^v\in {\mathbb{R}}^{s\times s\times K\times B}$ is a 4-D tensor stacked by all $P_k^v$, and ${\otimes }^{†}$ denotes the transposed convolution. Then, Equation (13) can be equivalently to

$$V^{\left(t\right)}=argmi{n}_V{\displaystyle \frac{1}{2}}\left|\right|V-(V^{(t-1)}-{\eta }_2\nabla h\left(V^{(t-1)}\right)){\left|\right|}_F^2+{\eta }_2{\lambda }_2{f}_2\left(V\right).$$

(16)

Equation (16) can also be optimized by the general proximal operator, which is referred to as

$$V^{\left(t\right)}=pro{x}_{{\eta }_2{\lambda }_2}(V^{(t-1)}-{\eta }_2\nabla h\left(V^{(t-1)}\right)),$$

(17)

where $pro{x}_{{\eta }_2{\lambda }_2}(·)$ is the proximal operator related to $f_2(·)$.

2.2.3. C-Subproblem

Equation (7) can be re-written as:

$$\begin{array}{cc}\hfill C^{\left(t\right)}=argmi{n}_C& \Vert \widehat{M}-\sum _{k=1}^K{D}_k^c\otimes C_k{\Vert }_F^2\hfill \\ & +\Vert \widehat{H}-\sum _{k=1}^K{P}_k^c\otimes C_k{\Vert }_F^2+{\lambda }_3{f}_3\left(C\right),\hfill \end{array}$$

(18)

where $\widehat{M}=M-{\sum }_{k=1}^K{D}_k^u\otimes U_k^{\left(t\right)}$, $\widehat{H}=\tilde{H}-{\sum }_{k=1}^K{P}_k^v\otimes V_k^{\left(t\right)}$. Furthermore, Equation (18) can be expressed as

$$C^{\left(t\right)}=argmin\left|\right|N-\sum _{k=1}^K{L}_k^c\otimes C_k{\left|\right|}_F^2+{\lambda }_3{f}_3\left(C\right),$$

(19)

where $N\in {\mathbb{R}}_{M\times N\times (B+b)}$ is the concatenation of $\widehat{H}$ and $\widehat{M}$, and $L_k^c\in {\mathbb{R}}^{s\times s\times (B+b)}$ is the concatenation of $D_k^c$ and $P_k^c$. Equation (19) can then be solved by the same way as optimizing U and V. As before, we first derive the quadratic approximation of Equation (19) with respect to C as follows:

$$\begin{array}{cc}\hfill C^{\left(t\right)}=argmi{n}_C& l\left(C^{(t-1)}\right)+{\displaystyle \frac{1}{2{\eta }_3}}{\Vert C-C^{(t-1)}\Vert }_F^2\hfill \\ & +\langle \nabla l\left(C^{(t-1)}\right),C-C^{(t-1)}\rangle +{\lambda }_3{f}_3\left(C\right).\hfill \end{array}$$

(20)

where

$$l\left(C^{(t-1)}\right)=\Vert N-\sum _{k=1}^K{L}_k^c\otimes C_k^{(t-1)}{\Vert }_F^2,$$

(21)

$$\nabla l\left(C^{(t-1)}\right)=L^c{\otimes }^{†}(\sum _{k=1}^K{L}_k^c\otimes C_k^{(t-1)}-N).$$

(22)

∇ is gradient operator, ${\eta }_3$ is the step size, $L^c\in {\mathbb{R}}^{s\times s\times K\times B}$ is a 4-D tensor stacked by all $L_k^c$, and ${\otimes }^{†}$ denotes the transposed convolution. Then, Equation (20) can be equivalently to

$$C^{\left(t\right)}=argmi{n}_C{\displaystyle \frac{1}{2}}\left|\right|C-(C^{(t-1)}-{\eta }_3\nabla l\left(C^{(t-1)}\right)){\left|\right|}_F^2+{\eta }_3{\lambda }_3{f}_3\left(C\right).$$

(23)

Equation (23) can also be optimized by the general proximal operator, which is referred to as

$$C^{\left(t\right)}=pro{x}_{{\eta }_3{\lambda }_3}(C^{(t-1)}-{\eta }_3\nabla l\left(C^{(t-1)}\right)),$$

(24)

where $pro{x}_{{\eta }_3{\lambda }_3}(·)$ is the proximal operator related to $f_3(·)$.

2.3. Network Architecture

Based on the model and solution proposed in the above two sections, we construct a deep neural network, which is shown in Figure 1. This network consists of T stages, which correspond to T iterations of the iterative algorithm used to solve Equation (4). These stages are connected sequentially and each phase completes an update of U, V, and C. And finally, the output of the last stage is the ultimate result of the three features.

As shown in Equations (12), (17), and (24), the three proximal operators $pro{x}_{{\eta }_1{\lambda }_1}(·)$, $pro{x}_{{\eta }_2{\lambda }_2}(·)$, and $pro{x}_{{\eta }_2{\lambda }_2}(·)$ play a important role in the algorithm. Usually, various hand-crafted priors are designed to make the proximal operators better and be closed form. However, due to the limited representation abilities, human-designed priors can not always achieve the expected effect. Therefore, we utilize convolutional neural networks to explore prior for U, V, and C, namely learning their corresponding proximal operators $pro{x}_{{\eta }_1{\lambda }_1}(·)$, $pro{x}_{{\eta }_2{\lambda }_2}(·)$, and $pro{x}_{{\eta }_2{\lambda }_2}(·)$ for updating $U^{\left(t\right)}$, $V^{\left(t\right)}$, and $C^{\left(t\right)}$ in each stage t. Specifically, we use ResNet to approximate the three proximal operators. Our network is designed as several stages, and at stage t, $U^{\left(t\right)}$ is updated using M, $U^{(t-1)}$, and $C^{(t-1)}$ by the following sub-steps:

$$U-Net:\left\{\begin{array}{c}{M}_c^{(t-1)}=\sum _{k=1}^K{D}_k^c\otimes C_k^{(t-1)},\hfill \\ M_u^{(t-1)}=\sum _{k=1}^K{D}_k^u\otimes U_k^{(t-1)},\hfill \\ {\epsilon }_M^{(t-1)}=M_c^{(t-1)}+M_u^{(t-1)}-M,\hfill \\ \nabla g^{(t-1)}=D^u{\otimes }^{†}{\epsilon }_M^{(t-1)},\hfill \\ U^{(t-0.5)}=U^{(t-1)}-{\eta }_1\nabla g^{(t-1)},\hfill \\ U^{\left(t\right)}=proxNe{t}_{{\Theta }_U^t}\left(U^{(t-0.5)}\right).\hfill \end{array}\right.$$

(25)

Each step is unfolded into one module of the U-Net shown in Figure 1. Then, $V^{\left(t\right)}$ is updated by $\tilde{H}$, $V^{(t-1)}$, and $C^{(t-1)}$:

$$V-Net:\left\{\begin{array}{c}{\tilde{H}}_c^{(t-1)}=\sum _{k=1}^K{P}_k^c\otimes C_k^{(t-1)},\hfill \\ {\tilde{H}}_u^{(t-1)}=\sum _{k=1}^K{P}_k^v\otimes V_k^{(t-1)},\hfill \\ {\epsilon }_{\tilde{H}}^{(t-1)}={\tilde{H}}_c^{(t-1)}+{\tilde{H}}_v^{(t-1)}-\tilde{H},\hfill \\ \nabla h^{(t-1)}=P^v{\otimes }^{†}{\epsilon }_{\tilde{H}}^{(t-1)},\hfill \\ V^{(t-0.5)}=V^{(t-1)}-{\eta }_2\nabla h^{(t-1)},\hfill \\ V^{\left(t\right)}=proxNe{t}_{{\Theta }_V^t}\left(V^{(t-0.5)}\right).\hfill \end{array}\right.$$

(26)

Similarly, each step is unfolded into one module of the V-Net as shown in Figure 1. Finally, $C^{\left(t\right)}$ is updated using M, $\tilde{H}$, $U^{\left(t\right)}$, and $V^{\left(t\right)}$

$$C-Net:\left\{\begin{array}{c}{F}_c^{(t-1)}=\sum _{k=1}^K{L}_k^c\otimes C_k^{(t-1)},\hfill \\ {\epsilon }_C^{(t-1)}=F_c^{(t-1)}-N,\hfill \\ \nabla l^{(t-1)}=L^c{\otimes }^{†}{\epsilon }_C^{(t-1)},\hfill \\ \nabla l^{(t-0.5)}=C^{(t-1)}-{\eta }_3\nabla l^{(t-1)},\hfill \\ C^{\left(t\right)}=proxNe{t}_{{\Theta }_C^{\left(t\right)}}\left(l^{(t-0.5)}\right).\hfill \end{array}\right.$$

(27)

Similarly, each step is unfolded into one module of the C-Net, as shown in Figure 1. In the formula above, ${\Theta }_U^{\left(t\right)}$, ${\Theta }_V^{\left(t\right)}$, and ${\Theta }_C^{\left(t\right)}$ are the parameters of the three ResNets at the $t^{th}$ stage. Note that all of the network parameters can be end-to-end learned from training data, including ${\{{\Theta }_U^{\left(t\right)},{\Theta }_V^{\left(t\right)},{\Theta }_C^{\left(t\right)}\}}_{t=1}^T, D^c={\left\{D_k^c\right\}}_{k=1}^K, D^u={\left\{D_k^u\right\}}_{k=1}^K, P^c={\left\{P_k^c\right\}}_{k=1}^K,$ $P^v={\left\{P_k^v\right\}}_{k=1}^K,$ $G^c={\left\{G_k^c\right\}}_{k=1}^K, G^u={\left\{G_k^u\right\}}_{k=1}^K, G^v={\left\{G_k^v\right\}}_{k=1}^K,{\eta }_1,{\eta }_2,{\eta }_3$. Note that our network is interpretable. For example, to update $U^{\left(t\right)}$ by Equation (25), computing $U^{(t-0.5)}$ means the gradient descent, and the $proxNe{t}_{{\Theta }_U^{\left(t\right)}}(·)$ performs nonlinear operation on $U^{(t-0.5)}$.

2.4. Network Training

The loss function of training process is defined as

$$L=\sum _{j=1}^J\Vert O_p^j-O_g^j{\Vert }_2^2,$$

(28)

where $O_p^j$ is the predicted HRHS image of the $j^{th}$ training sample, $O_g^j$ represents the reference HRHS image of the $j^{th}$ training sample, and J is the number of training pairs. The proposed network is trained on Python 3.9.0 with Tensorflow 2.5.0 and Linux operating system with NVIDIA GPU GeForce GTX 2080Ti. In the network, all the convolutional kernels are set as 8 × 8, and the stage number of our network is set as 3. and the C-Net contains 3 ResBlocks, while the U-Net and V-Net each containing 5 ResBlocks. During optimizing the model, we chose the Adam algorithm and a decay technique to set the learning rate, i.e., decayed by 0.95 every 20,000 iterations with a fixed initial learning rate 0.00005. The epoch number is 95, and the batch size is 3.

3. Experimental Results

To verify the efficiency of the proposed model, we compare the proposed network with several state-of-the-art methods on three datasets. The performance assessment contains two cases. In the first case where the ground-truth HRMS is known, we utilize popular reference-based indexes. In the second case where the ground-truth HRMS is not available, we use $D_{lambda}$, $D_s$, and the quality without reference ($QNR$) criteria [41]. Specifically, the competing methods consist of GSA [7], SFIM-HS [8], GLP-HS [9], CNMF [11], S3RHSR [12], HySure [10], FUSE [15], UTV [17], CSTF [16], MHFnet [35], and HSRnet [31].

The three selected datasets are the CAVE database [42], Harvard database [43], and Chikusei database [44], respectively. Specifically, the CAVE database consists of 32 HSIs with 512 × 512 and 31 spectral bands between 400 nm to 700 nm at 10 nm steps. Meanwhile, every HSI has a corresponding RGB image with 512 × 512 and three spectral bands. The Harvard dataset is a public dataset that contains 77 HSIs of indoor and outdoor scenes covering various kinds of objects and buildings. Each HSI has a spatial size of 1392 × 1040 with 31 spectral bands acquired at an interval of 10 nm in the range of 420–720 nm. The Chikusei dataset is an airborne hyperspectral dataset, which was taken by Headwall Hyperspec-VNIR-C imaging sensor over agricultural and urban areas in Chikusei, Ibaraki, Japan. It only has one image with 128 bands in the spectral range from 363 nm to 1018 nm and the scene consists of 2517 × 2335 pixels.

We use the following six complementary and widely used quality measures for quantitative fusion assessment, i.e., the peak signal-to-noise ratio (PSNR), the spectral angle mapper (SAM [45]), the erreur relative globale asimensionnelle de synthèse (ERGAS [46]), the structure similarity (SSIM [47]), spatial correlation coefficient (SCC), and $Q{2}^n$ [48]. Let $x\in R^{B\times P}$ be the reference image with B bands and P pixels. $x={[x_1,\dots ,x_B]}^T=[x_1,\dots ,x_j]$, where $x_i\in R^{P\times 1}$ is the ith band ($i=1,\dots ,B$) and $x_j\in R^{B\times 1}$ is the spectral signature of the jth pixel ($j=1,\dots ,P$). $\widehat{x}$ denotes the estimated HS image. The PSNR measures the average qualities of spatial reconstruction between all bands of the target and of the reference image. It is the ratio of the maximum square of a signal to the square of relative residual errors. The PSNR is defined as

$$PSNR(x,\widehat{x})={\displaystyle \frac{1}{B}}\sum _{i}10·lo{g}_{10}{\displaystyle \frac{max{\left(x_i\right)}^2}{\Vert x_i-{\widehat{x}}_i{\Vert }_2^2/P}}.$$

(29)

The SAM index [45] represents the spectral similarity between reference data and reconstructed hyperspectral data. More precisely, The spectral distance is determined by calculating the angle between two vectors of the estimated and reference spectra. The SAM index is defined as

$$SAM(x,\widehat{x})={\displaystyle \frac{1}{P}}\sum _{i}arccos\left({\displaystyle \frac{x_j^T{\widehat{x}}_j}{\Vert x_j{\Vert }_2{\Vert {\widehat{x}}_j\Vert }_2}}\right).$$

(30)

The SSIM [47] provides a criterion to evaluate the similarity of images from three aspects: brightness, contrast, and structure. Usually, its calculation formula is described as

$$SSIM(x,\widehat{x})={\displaystyle \frac{1}{B}}\sum _i{\displaystyle \frac{({\mu }_{x_i}{\mu }_{{\widehat{x}}_i}+c_1)(2{\sigma }_{x_i{\widehat{x}}_i}+c_2)}{({\mu }_{x_i}^2+{\mu }_{{\widehat{x}}_i}^2+c_1)({\sigma }_{x_i}^2+{\sigma }_{{\widehat{x}}_i}^2+c_2)}}.$$

(31)

The ERGAS [46] conducted an overall evaluation of space and spectrum. ERGAS calculates the band-wise normalized root-mean-square error(RMSE) and multiplies it with the ratio of the ground sampling distance, which is defined as

$$ERGAS(x,\widehat{x})=100d\sqrt{{\displaystyle \frac{1}{B}}\sum _{i=1}^B{\displaystyle \frac{\Vert x_i-{\widehat{x}}_i{\Vert }_2^2}{{\left(\frac{1}{P}{1}_P^T{x}_i\right)}^2}}}.$$

(32)

The SCC is a method of assessing the similarity between the estimated and reference data. It can be achieved by calculating the correlation coefficient of the corresponding pixel values between the two images. Specifically, the calculation formula is as follows

$$SCC(x,\widehat{x})={\displaystyle \frac{1}{P}}\sum _{j=1}^P{\displaystyle \frac{cov(x_j,{\widehat{x}}_j)}{{\sigma }_{x_j}{\sigma }_{{\widehat{x}}_j}}}.$$

(33)

$Q{2}^n$ [48] is an extension of the universal image quality index ($UIQI$ [49]) to HS image depended on hypercomplex numbers. The $UIQI$ is designed for monochromatic images to measure any distortion from three factors, i.e., loss of correlation, luminance distortion, and contrast distortion, and by the product method of the three terms. The $UIQI$ between one reference image band (z) and its corresponding target image band (y) can be written as

$$UIQI(z,y)={\displaystyle \frac{4{\sigma }_{zy}\overline{z}\overline{y}}{({\sigma }_z^2+{\sigma }_y^2)({\overline{z}}^2+{\overline{y}}^2)}}.$$

(34)

To overcome this limitation and additionally consider the spectral distortion, $Q{2}^n$ extends the $UIQI$ based on modeling each pixel spectrum ($x_j$) as a hypercomplex number [48], expressed as

$$x_j=x_{j,0}+x_{j,1}{i}_1+x_{j,2}{i}_2+\dots +x_{j,2^{n-1}}{i}_{2^n-1}.$$

(35)

The smaller ERGAS and SAM are, and the larger PSNR, SSIM, SCC and $Q{2}^n$ are, the better the fusion result is. The ideal values for SSIM and SCC are one.

Additionally, we also adopt $D_{\lambda },D_s$, and $QNR$ [50,51] as the indexes to measure the qualities of data with no reference at full resolutions. The spectral distortion, $D_m$, is determined by comparing the LRHS images with the fused HS images [51]. Interband $UIQI$ values are computed in two distinct sets for low and high resolutions. By subtracting the corresponding $UIQI$ values obtained at the two scales, the spectral distortion caused by the pansharpening process is derived. The spatial distortion, $D_s$, is determined by calculating the $UIQI$ values between each HS band and the MS band that has been degraded to the resolution of the HS images, as well as between the fused HS images and the full-resolution MS image. The absolute difference between the corresponding $UIQI$ values, averaged over all bands, gives the spatial distortion $D_s$ [51]. The two distortions are then combined to provide a unique quality index, referred to as $QNR\in [0,1]$, with 1 being the best attainable value:

$$QNR={(1-D_{\lambda })}^{\alpha }{(1-D_s)}^{\beta },$$

(36)

where usually $\alpha =\beta =1$.

3.1. Evaluation at Reduced Resolution

We measure the performance of our network in the reduced resolution case. We conduct experiments on two datasets, CAVE and Chikusei. In the experiment on CAVE dataset, we did not use all the data, but chose 20 images as the training set and randomly selected 10 patches of size 200 × 200 from the other 11 images for testing. The 11 CAVE test images are shown in Figure 2.

We extracted 18,000 overlapped patches with a size of 32 × 32 × 31 from the 20 images of the CAVE dataset used as GT of training dataset, thus forming the HR-HSI patches. Accordingly, the LR-HSI patches are generated starting from the HR-HSI by applying a Gaussian blur with kernel size equal to 3 × 3 and standard deviation equal to 0.5 and then downsampling the blurred patches to the size of 8 × 8, i.e., with a downsampling factor of 4. The corresponding RGB patches are generated based on the spectral response function of the Nikon D700 camera. Thus, 18,000 patches of size of 32 × 32 × 3 are available to represent the HR-MSI of training dataset. We perform the experiments on all the 10 testing reduced images.

Table 1. Average $QIs$ of the results on 10 patches of size 200 × 200 on the CAVE datasets. The best values are highlighted in boldface and the second place is underlined.

Method	PSNR	SSIM	SAM	ERGAS	SCC	Q2n
GSA [7]	37.578	0.976	6.621	3.483	0.992	0.717
SFIM-HS [8]	30.179	0.941	10.140	22.596	0.967	0.645
GLP-HS [9]	32.956	0.962	6.106	5.136	0.985	0.670
CNMF [11]	36.647	0.971	8.251	3.847	0.985	0.708
S3RHSR [12]	31.564	0.922	19.770	7.500	0.977	0.657
HySure [10]	35.879	0.968	7.274	4.056	0.989	0.742
FUSE [15]	27.247	0.934	5.035	9.248	0.955	0.601
UTV [17]	29.952	0.815	12.209	7.988	0.959	0.441
CSTF [16]	33.694	0.954	7.326	4.948	0.985	0.621
HSRnet [31]	42.796	0.994	2.814	2.147	0.995	0.817
MHFnet [35]	41.308	0.987	5.573	2.802	0.989	0.732
Ours	43.485	0.993	3.319	1.986	0.995	0.799

presents the average $QIs$ on the 10 images. It can be seen that, in Table 1, the DL-based methods (i.e., HSRnet, MHFnet and our network) have better performance than other traditional methods, which may be due to the powerful learning capabilities of DL methods. Additionally, our method performs best in terms of the PSNR, ERGAS, and SCC indexes, and is comparable to the best value on the remaining metrics. The visual comparison shown in Figure 3, implying our network performs almost the best in spatial domain, further verifies the numerical results. For better illustration, we enlarge two regions of the fused image. Meanwhile, the selected spectral vectors of the sample which is shown in Figure 3 are also presented for the fused results coming from the different fusion methods and the GT. It is worth to be remarked that the spectral vectors estimated by our method and the GT ones are very close to each other.

To demonstrate the performance of our MHF-CSCNet on remote sensing HSIs, we conducted an experiment on the Chikusei dataset, which is taken over agricultural and urban areas in Chikusei, Ibaraki, Japan. We take the original data as the GT HRHS and simulate the LRHS in the same way as the previous experiments. As for the HRMS, we use the corresponding RGB image obtained by Canon EOS 5D Mark II together with the HRHS. Then, we selected the 1000 × 2200 area in the upper left corner of the image for training, and extracted 64 × 64 overlapping patches from it as training GT HRHS patches. Moreover, the input HRMS and LRHS patches are of size 64 × 64 × 3 and 16 × 16 × 128, respectively. Additionally, we crop six nonoverlap 680 × 680 × 128 area from the remaining part of the Chikusei dataset for testing.

Table 2. Average $QIs$ of the results for six testing images on the Chikusei dataset. The best values are highlighted in boldface and the second place is underlined.

Method	PSNR	SSIM	SAM	ERGAS	SCC	Q2n
GSA [7]	32.519	0.881	3.745	5.511	0.698	0.630
SFIM-HS [8]	28.468	0.885	4.392	60.326	0.782	0.617
GLP-HS [9]	32.418	0.917	3.288	5.717	0.922	0.647
CNMF [11]	31.163	0.883	3.457	6.059	0.876	0.560
S3RHSR [12]	25.936	0.751	10.115	12.367	0.749	0.441
HySure [10]	30.464	0.842	3.716	7.020	0.606	0.555
FUSE [15]	31.871	0.906	3.251	6.170	0.759	0.627
UTV	25.085	0.652	9.130	12.253	0.712	0.261
CSTF [16]	29.701	0.829	4.853	7.376	0.871	0.460
HSRnet [31]	35.685	0.948	2.131	3.668	0.955	0.778
MHFnet [35]	33.611	0.938	3.587	5.345	0.936	0.659
Ours	36.462	0.955	2.242	3.174	0.959	0.753

shows the average $QIs$ for the six testing images on all methods. It can clearly be seen that our MHF-CSCNet has the best performance on almost all indices compared to all other methods. These results indicate that our method has excellent performances on reconstruction of spatial and spectral senses. Although our approach is slightly worse than HSRnet (which directly preserves the spectrum from the hyperspectral images and then supplements it with spectral information from the high-resolution images) on SAM metrics, it is still significantly better than other approaches. It shows that our method not only has a good protection for the reconstruction of spatial details, but also for the spectrum. We also exhibit the pseudocolor images of the obtained outcomes for a visual comparison in Figure 4. Obviously, the fused results of our MHF-CSCNet either exceeds some methods in spatial feature recovery, or outperforms in spectral protection. At the same time, we show the spectral vector at a certain position in Figure 4. It can be observed from Figure 4a,b that the spectral vectors obtained by the SFIMHS, HySure, FUSE, UTV, and CSTF methods have a large degree of error, While HSRnet, MHFnet and our methods have a slight spectral distortion. Specifically, our method and HSRnet both protect the spectrum well. The effectiveness of our method is checked again.

3.2. Evaluation at Full Resolution

This section evaluates the performance of our network in the full resolution case. We conduct this experiment on one of complete test image pairs from the CAVE test set. The image pair contains one HS image ($128\times 128\times 31$), and one MS image ($512\times 512\times 3$). Since the ground-truth HRHS image is not available, we thus adopt $D_{\lambda },D_s$, and $QNR$ as the indexes. The quantitative results of all competing methods. are reported in

Table 3. Average $QIs$ of the results on one image (clay) of CAVE test set. The best values are highlighted in boldface and the second place is underlined.

Method	${\mathit{D}}_{\mathit{\lambda }}$	${\mathit{D}}_{\mathit{s}}$	$\mathit{QNR}$
GSA [7]	0.1096	0.0948	0.806
SFIM-HS [8]	0.1712	0.0602	0.7789
GLP-HS [9]	0.0778	0.0891	0.8401
CNMF [11]	0.0704	0.0786	0.8565
S3RHSR [12]	0.1512	0.1484	0.7229
HySure [10]	30.0789	0.0885	0.8396
FUSE [15]	0.0223	0.1022	0.8778
UTV [17]	0.1046	0.0545	0.8466
CSTF [16]	30.0741	0.0631	0.8675
HSRnet [31]	0.0171	0.1027	0.882
MHFnet [35]	0.1284	0.1124	0.7736
Ours	0.0222	0.0967	0.8832

. As can be seen, compared with DL-based methods, our MHF-CSCNet has the best evaluation indexes. However, some traditional methods are superior than DL-based methods on indices $D_{\lambda }$ and $D_s$. This phenomenon may be caused by the characteristics of different traditional methods. Figure 5 shows the composite images of a test sample obtained by the competing methods, with bands 15-8-29 as R-G-B. It is seen that the composite image obtained by MHF-net is clearer than the fused results of the other methods, while the results of other methods usually contain blurred boundary.

3.3. Generalization to New Satellites

In this section, we randomly choose the 512 × 512 size in the lower right part of the image, then twenty images are randomly determined for testing from the Harvard dataset. As in the previous settings, the original data are regarded as the GT HRHS. The LRHS are simulated by the same method as in the CAVE dataset. The HRMS is also obtained by applying the spectral response of the Nikon D700 camera. In the experiment on the Harvard dataset, we would like to remark that our method is trained on the CAVE dataset, and directly test them without any retraining or fine-tuning. Thus, the performance on the Harvard dataset of these methods could reflect the model’s generalization abilities.

Table 4. Average $QIs$ of the results for twenty testing images on the Harvard dataset. The best values are highlighted in boldface and the second place is underlined.

Method	PSNR	SSIM	SAM	ERGAS	SCC	Q2n
GSA [7]	31.926	0.952	3.350	5.796	0.988	0.738
SFIM-HS [8]	28.701	0.916	3.446	7.200	0.971	0.665
GLP-HS [9]	29.558	0.928	3.630	6.671	0.977	0.678
CNMF [11]	31.440	0.949	3.065	6.027	0.984	0.734
S3RHSR [12]	30.420	0.920	5.778	7.001	0.977	0.683
HySure [10]	31.240	0.948	3.341	6.141	0.986	0.730
FUSE [15]	27.803	0.876	3.242	7.636	0.952	0.516
UTV [17]	33.170	0.931	4.722	4.351	0.975	0.647
CSTF [16]	33.927	0.930	3.660	4.042	0.978	0.654
HSRnet [31]	38.634	0.974	3.159	2.828	0.985	0.752
MHFnet [35]	37.820	0.969	4.412	4.747	0.978	0.736
Ours	38.930	0.975	3.405	3.585	0.986	0.755

records the average $QIs$ over 20 images of all competing methods. In the four indicators, PSNR, SSIM, ERGAS and Q2n, the method based on deep learning is significantly better than the traditional method as shown in Table 4. Figure 6 shows the pseudo-graph composed of bands 8-15-29 on all the completing methods. It is easy to observe that the proposed method performs better than other competing ones, in the better recovery of both finer-grained textures and spectral distribution. Additionally, we plot the spectral vectors for two exemplary cases in Figure 7. It is worth to be remarked that the spectral vectors estimated by our method and the GT ones are very close to each other. Similar to the results in the CAVE dataset, we can observe the superiority of our method on almost all evaluation measures, which implies that our network can generalize well to new satellites.

4. Parameters and Time Comparisons

Table 5. Comparison of running time, parameter number, and FLOPs.

Method	Image Size	HSRnet	MHFnet	Our
Parameters ($\times {10}^6$)	128 × 128 × 31,512 × 512	0.63	1.21	5.52
FLOPs ($\times {10}^{11}$)	128 × 128 × 31,512 × 512	2.57	5.10	52.33
Running time (s)	128 × 128 × 31,512 × 512	0.47	12.91	23.24

shows the parameter number and floating-point operations per second (FLOPs) of our network with the other two DL-based methods. All the methods used 20 images of Harvard dataset as same as the test data the last subsection for inference with NVIDIA GPU GeForce GTX 2080Ti. As Table 5 shows the parameter number, running time and FLOPs of our model are larger than other methods, although our method has better results, and we will try to improve this problem in our future work.

5. Conclusions

This paper proposes a novel model-based deep neural network for the MS/HS fusion task. Firstly, a new MS/HS fusion observation model based on convolutional sparse coding (CSC) is proposed. Then, we design a proximal gradient algorithm to solve this model, and the algorithm is unfolded into a network. Our network has a clear physical interpretation because each module corresponds to a specific operation of the algorithm. Experimental results conducted on some benchmark datasets demonstrate that our network can obtain comparable or better quantitative and qualitative performance compared to other methods.