Unsupervised Hyperspectral and Multispectral Images Fusion Based on the Cycle Consistency

Conventional HSI and MSI fusion methods generally assume the imaging model is known, where both the parameters of point spread function (PSF) in the spatial degradation model from HrHSI to LrHSI, and spectral response function (SRF) in the spectral merging processes from HrHSI to HrMSI, are known or at least one of them is known. In practice, however, the spatial and spectral downsampling processes are complex and it is non-trivial to accurately know these parameters [9]. Therefore this type of traditional model has limited performance when downsampling parameters mismatch the actual system.

Deep learning-based approaches are introduced to fusion processes to recover the spatial and spectral details of HrHSI and obtain promising results, which often use HrHSI as ground truth to supervise the training process. However, HrHSI is usually inaccessible in practical fusion scenarios. Moreover, these methods require a number of paired images as training data. Likewise, strictly paired training data is hard to obtain, i.e., LrHSI and HrMSI may be non-rigidly aligned in practice. The above two facts limit the practical generalization ability of deep learning-based models for the HSI and MSI fusion.

Recently, cycle-consistent generative adversarial networks (CycleGAN) [10] has been successfully used for cross-domain transformation tasks, e.g., image style transfer [11], unsupervised language translation[12], image superresolution[13], near-infrared (NIR) image to RGB image transformation [14]. This method, also known as unsupervised distribution alignment, works in an unsupervised manner and does not impose the constraint of using paired training data. Inspired by CycleGAN, we propose an HSI and MSI fusion model based on cycle consistency during the transformation between domain LrHSI and domain HrMSI, named CycFusion. Instead of using known PSF and SRF as priors in traditional methods, or known HrHSI as the ground truth to supervisely train deep fusion networks, our proposed CycFusion can estimate the image degradation processes and then implement the fusion task in an unsupervised manner. Specifically, the proposed model consists of two downsampling modules and two upsampling modules. The downsampling modules contain a shared weight depthwise separable convolution network to estimate the kernel of PSF in the spatial degradation process, and another $1\times 1$ convolution module to estimate the SRF matrix in the spectral degradation model. Inspired by the PSF and SRF estimation processes[15], we draw above downsampling modules to respectively degrade HrMSI in spatial dimension and LrHSI in spectral dimension to obtain two low resolution MSI (LrMSI), then constrain the consistency of these two images to learn the downsampling processes. The estimated results are used as pre-training parameters of CycFusion. Besides, the other two upsampling modules, namely, the spatial superresolution module and spectral superresolution module, which restore the spatial details and spectral signatures of desired results through the mapping from LrHSI to HrHSI and from HrMSI to HrHSI, respectively. After building the above modules, we can transfer an image in the domain LrHSI to the image that contains same scene in the domain HrMSI through spatial upsampling and spectral downsampling modules and the same is true for the reverse conversion via another two modules. To clarify, the proposed CycFusion is shown in Fig. 1. After perform the training strategy that will be detailed in the Section III, the fused results will be obtained by extracting the feature maps in the intermediate transformation.

The main contributions of our proposed model are summarized as follows.

1. 1.

   Inspired by the advanced progress of the domain transformation, a novel unsupervised hyperspectral and multispectral fusion model is proposed.

2. 2.

   The observation model is embedded in the CycFusion, and the estimated results make it possible to implement the blind fusion task using deep neural networks with high model capacity.

3. 3.

   The proposed model consists of cross-domain transformations between LrHSI and HrMSI. The experiment results show that the fused images benefit from the objective function with cycle consistency between the input LrHSI and HrMSI, and themselves through twice domain transformations.

The rest of this paper is organized as follows. Section II introduces related work. The proposed CycFusion is presented in Section III. In Section IV, the effectiveness of CycFusion is demonstrated through experiments on three publicly available datasets. Section V concludes this paper and gives future directions.

Pansharpening algorithms [16] aim at injecting the spatial details of high-resolution panchromatic images into LrHSI, including component substitution (CS) approaches and multi-resolution analysis (MRA) approaches. These methods have been extended to the HSI-MSI fusion. Gram–Schmidt adaptive (GSA) [17] is a representative CS-based pansharpening algorithm, which uses Gram–Schmidt transformation to separate the spatial component of the LrHSI that needs to be substituted by the HrMSI. The MRA-based algorithms [18] produce the desired HrHSI by injecting the spatial details obtained by the subtraction of HrMSI and its low-pass version. The generalized Laplacian pyramid-based hyper-sharpening (GLP-HS) [19] is the representative method that uses the pyramidal decompositions to obtain the high spatial resolution structures. Pansharpening-based HSI-MSI fusion methods are computationally efficient and do not rely on known PSF and SRF, however, these approaches compute the fused results band-by-band, which may introduce the spectral distortion.

Subspace-based HSI-MSI fusion methods leverage the low-rankness and sparsity of the desired HrHSI, which generally assume there is a subspace that can represent all spectral signatures of HrHSI. Among the representative approaches are unmixing-based methods and orthogonal subspace-based methods. Coupled non-negative matrix factorization (CNMF) [20] pioneered the unmixing-based fusion, which embeds the linear mixing model (LMM) in the fusion problem and iteratively optimizes the endmembers and abundance maps of the HrHSI. Orthogonal subspace-based fusion methods assume the HrHSI underlies a low dimensional and orthogonal subspace. These methods generally require proper regularization on the coefficients of HrHSI on the subspace due to the under-estimate problem [21] caused by the dimensions of the orthogonal subspace may large than the bands of HrMSI. Commonly used regularization imposes smoothness and sparsity priors on the coefficients, e.g., hyperspectral superresolution (HySure) [15] constrains the smoothness of HrHSI in the subspace using the band-by-band total variation (TV), the work in [22] using dictionary learning with sparse coefficients, and a non-negative structured sparse representation (NSSR) [23] method models the spatial correlations via clustering. Furthermore, the convolutional neural networks (CNNs) denoisers have been introduced to produce reasonable results which substitute using the regularized term, named CNN-Fus [24]. It is worth mentioning that the work in [25] implements the fusion process by solving the Sylvester equation (FUSE), which greatly improved the fusion efficiency. Beyond solving the 2D matrix-based optimization problems, tensor-based methods direct process the 3D HSI data cube. The representative approach of this class is the coupled sparse tensor factorization (CSTF) [26] method based on Tucker decomposition. Subspace-based HSI-MSI fusion methods have interpretable model design and physical meaningful regularizers, leading to better superior performance over pansharpening-based methods. However, these models generally need PSF and SRF as priors, which limits their application in real scenarios. Besides, these methods only adopt the linear mixing model and linear subspace resulting in restrictive models.

Recently, deep learning as an expressive model has been introduced to HSI-MSI fusion [27]. The deep learning-based fusion models can be divided into two parts, namely, supervised methods and unsupervised methods. The supervised fusion methods use an amount of paired LrHSI and HrMSI as inputs and HrHSI as ground truth. Once properly trained, the fusion models can fuse two degraded images in the test dataset and output the desired HrHSI. Most deep learning-based methods focus on designing the expressive modules to extract the high-resolution spatial information and spectral signatures and then fuse them into output images. Wang et al.[28] introduced the deep prior into HSI super-resolution to learn the spatial-spectral priors automatically. Xie et al. [29] proposed a model-inspired MSI-HSI fusion network (MHF-net) which enhanced the interpretability in estimating the observation model from the training data. Hu et al. [30] designed a deep spatial and spectral attention CNN to implement the hyperspectral image super-resolution, named HSRnet. Wang et al.[31] proposed an alternately and iteratively optimization algorithm to estimate both the degradation model and fusion model, called enhanced deep blind hyperspectral image fusion network (EDBIN). Instead of using the known PSF and SRF as in subspace-based models, deep learning-based models leverage the high capacity of deep neural networks to recovery the HrHSI from low-resolution images and obtain state-of-the-art performance. However, these models need large numbers of HrHSI as training data that cannot be acquired in the real world. Fortunately, unsupervised fusion methods have been proposed to reduce the reliance on HrHSI as training data. In [32], the authors assume LrHSI and HrHSI have similar abundances and propose a two-branch unsupervised sparse Dirichlet-Net (uSDN) that iteratively learns the shared features of abundance. Then, a nonlinear variational probabilistic generative model (NVPGM) [33] extends the fusion model that globally trains the model parameters and obtains more accurate fusion results. Moreover, the CNNs have been introduced to unsupervised fusion models to learn the spatially correlated structures and further improve the fusion quality, like FusionNet [34]. Guide deep decoder (GDD) [35] based on the deep image prior was proposed to produce the fused images from a noise input with the degraded images. A coupled unmixing model with a cross-attention module is embedded into the fusion network (CUCaNet) [36], which can learn the unknown SRF and PSF. These unsupervised fusion models are independent of the ground truth of HrHSI as training data, which have promising applications in practice.

To summarize, the prior information of representative fusion methods are listed in Table I.

Given a pair of observed LrHSI cube $\mathcal{Y}\in\mathbb{R}^{L\times w\times h}$ and HrMSI cube $\mathcal{Z}\in\mathbb{R}^{l\times W\times H}$, the objective of HSI-MSI fusion is to aggregate the spectral and spatial information of these two data cubes and then produce the desired HrHSI cube $\mathcal{X}\in\mathbb{R}^{L\times W\times H}$, where $\{W,H,L\}$ and $\{w,h,l\}$ denote the widths, heights and band numbers of high-resolution and low-resolution images. Generally, $w\ll W,h\ll H$ and $l\ll L$ and thus the fusion is ill-posed. When unfolding the 3D tensor to a 2D matrix, most of the existing literature assumes the observation model is linear as follows

where $\mathbf{X}\in\mathbb{R}^{L\times WH}$, $\mathbf{Y}\in\mathbb{R}^{L\times wh}$ and $\mathbf{Z}\in\mathbb{R}^{l\times WH}$ are HrHSI, LrHSI and HrMSI matrices, respectively. $\mathbf{B}\in\mathbb{R}^{WH\times WH}$ is the bluring matrix constructed by the PSF kernel. $\mathbf{S}\in\mathbb{R}^{WH\times wh}$ and $\mathbf{R}\in\mathbb{R}^{l\times L}$ are the spatial downsampling matrix and the SRF matrix. $\mathbf{N}_{y}$ and $\mathbf{N}_{z}$ are the additive noise matrices. Note that the above observation model (1) is adopted in the most of existing HSI-MSI fusion literature [26, 23]. The task of HSI-MSI fusion is to recover the matrix, $\mathbf{X}$, using $\mathbf{Y},\mathbf{Z}$ as inputs, which may use known $\mathbf{B,S}$ and $\mathbf{R}$ in the non-blind fusion or use estimated values of them in the blind fusion.

Based on the above observation model, we propose two CNN-based modules, namely, the spatial downsampling module and spectral downsampling module, to estimate the parameters in the observation model for performing the blind fusion task.

First, considering the spatial degradation process is imposed only in the spatial dimension and is independent of the spectral dimension. Inspired by it, we introduce a depthwise separable CNN [37] with a shared kernel among all channels of size equal to the stride size and also to the spatial downsampling ratio, $S$, to model the PSF kernel. The spatial degradation process is then formulated by a convolutional operator as

where $\mathcal{X}_{i,:,:}$ means to take the $i$-th slice from the HrHSI tensor and also the $i$-th channel image in the data cube and $[\dots,\dots,\dots]$ means concatenating matrices into a tensor. The convolutional operator is simplified by $\mathcal{F}_{y}(\cdot)$ parameterized by a shared convolutional kernel $\mathbf{K}\in\mathbb{R}^{S\times S}$. In addition, the kernel parameter should ensure the sum-to-one constraint to satisfy the energy conservation in the spatial degradation model as

and the Softmax function is added to the model parameter $\boldsymbol{\Theta}^{(PSF)}\in\mathbf{R}^{S\times S}$ to ensure the above constraint as shown in the top left of Fig. 1.

Second, the opposite of the spatial degradation process is the spectral degradation process only performing on the spectral dimension leading to no spatial correction information bothered here. So we use a $1\times 1$ CNN to model the SRF between HrHSI and HrMSI, where the spectral downsampling module in (1) can be formulated by a pointwise convolutional operator as

where the convolutional operator is simplified by $\mathcal{F}_{z}(\cdot)$, the SRF matrix $\mathbf{R}$ is parameterized by $\mathbf{\Theta}^{(SRF)}\in\mathbb{R}^{l\times L}$ and the Softmax function is also used here to ensure the energy conservation as

The pointwise spectral downsampling module is depicted in the bottom right of Fig. 1.

To summarize, inspired by the observation model, we build two single-layer CNNs to estimate the downsampling domain transformation from HrHSI to LrHSI and HrMSI with the operators and kernel sizes listed in Table II.

Furthermore, we propose two upsampling modules to implement the inverse transformation from LrHSI to HrHSI and HrMSI to HrHSI, namely, spatial upsampling module and spectral upsampling module. The former consists of an interpolation block and a refining block. The interpolation block alternately and iteratively increases the spatial resolution of the input LrHSI correspond to HrHSI and aggregate spatially correlated information via bicubic interpolation and $3\times 3$ convolution operation, respectively. The input LrHSI will pass through several interpolation blocks with different learnable parameters during the forward computation. Then followed by a refining block with two $1\times 1$ convolutional layers further preserves the high-fidelity details of HrHSI. Additionally, the instance normalization (IN) layers[38] are plugged into the spatial upsampling module to accelerate the network training process and the rectified linear unit (ReLU) is adopted as the activation function in our proposed model provided as

The Sigmoid function, $\mathbf{\sigma}(x)=1/(1+\exp{({-x})})$, is at the end of this module to ensure that the outputs are in [0,1] where data are normalized in the range when the network is trained. To clarify, the forward computation of the spatial upsampling module can be formulated as

where $\mathbf{Y}^{(0)}=\mathbf{Y}$, $\uparrow_{2}$ represents the spatial resolution is doubled by each interpolation step, $i=0,1,\dots,4$ due to we conducted on data with a spatial downsampling ratio of 32 in the next section and $\mathbf{X}_{Spa}$ is the spatial superresolution result. We denote this module as $\mathbf{X}_{Spa}=\mathcal{G}_{y}(\mathbf{Y})$ for short and the diagram is depicted in the bottom left of Fig. 1.

Besides, we propose another spectral upsampling module which is a 5-layer $1\times 1$ CNNs mapping HrMSI to HrHSI given as

where $\mathbf{Z}^{(0)}=\mathbf{Z}$, $i=0,1,\dots,4$ and $\mathbf{X}_{Spe}$ is the spectral superresolution result. We denote this module as $\mathbf{X}_{Spe}=\mathcal{G}_{z}(\mathbf{Z})$ for short and the diagram is depicted in the top right of Fig. 1.

These two upsampling modules are parameterized by the neural networks with parameters in (9) and (10) denoted by $\boldsymbol{\Theta}^{(SpaUp)}$ and $\boldsymbol{\Theta}^{(SpeUp)}$, respectively. To clarify, the network architecture of CycFusion for the CAVE dataset is listed in Table II.

In this part, we construct the objective function based on learning the forward and reverse domain transformations between LrHSI and HrHSI. Specifically, the combination of $\mathcal{G}_{y}(\cdot)$ and $\mathcal{F}_{z}(\cdot)$ completes the forward mapping, while the reverse transformation consists of $\mathcal{G}_{z}(\cdot)$ and $\mathcal{F}_{y}(\cdot)$. The objective function of CycFusion contains three parts, namely, marginal matching, cycle consistency and fusion identity.

First, single domain transformation between LrHSI and HrMSI is supposed to satisfy the marginal matching, e.g., an image from domain LrHSI, $p(\mathbf{Y})$, through the forward mapping is subject to the distribution of HrMSI, $p(\mathbf{Z})$. And it holds for images in the domain HrMSI via the reverse transformation. These relationships are given as

where $q_{1}(\mathbf{Z})$ and $q_{1}(\mathbf{Y})$ are learned distributions of HrMSI and LrHSI via the single domain transformation, respectively. In the HSI-MSI fusion task, paired $\mathbf{Y}$ and $\mathbf{Z}$ are generally used and (11) can be rephrased to the first term of the loss function as

where $\mathbf{\Theta}=\{\boldsymbol{\Theta}^{(PSF)},\boldsymbol{\Theta}^{(SRF)},\boldsymbol{\Theta}^{(SpaUp)},\boldsymbol{\Theta}^{(SpeUp)}\}$ contains all parameters that need to be optimized, and $\|\cdot\|_{1}$ denotes $\ell_{1}$-norm of a matrix/vector.

Second, dual-domain transformations between LrHSI and HrMSI are supposed to satisfy the cycle consistency, e.g., an image from domain LrHSI transfers to domain HrMSI and then transfers back close to itself. And it also holds for images in the domain HrMSI. The cycle consistency can be formulated as

where $q_{2}(\mathbf{Z})$ and $q_{2}(\mathbf{Y})$ are learned distributions of HrMSI and LrHSI via dual-domain transformations, respectively. We also use the $\ell_{1}$ loss to construct the second term of the objective function and (13) can be rephrased as

Third, the corresponding images in the domain HrHSI can be obtained from single-side superresolution. We constrain the identity for fused results via the following loss function

To sum up, the objective function of CycFusion is

When the CycFusion training is complete, we simply obtain the fused result by

The parameters of the PSF and SRF may be hard-acquired in real applications. An estimation method for the PSF kernel and the SRF matrix was proposed in Hysure [15], which further downsample the LrHSI and the HrMSI in the spectral and spatial dimensions, respectively, to produce the same low spatial resolution MSI (LrMSI). A similar strategy is used here to learn the downsampling parameters, $\boldsymbol{\Theta}^{(PSF)}$ and $\boldsymbol{\Theta}^{(SRF)}$. Specifically, we use the following loss function to pre-train the downsampling modules and then load the pretraining weights when training the whole network.

To specify, the training process of CycFusion is summarized in Algorithm 1. When we know the PSF kernel and SRF matrix as priors, the proposed model can skip the pretraining phase and implements the fusion task in a no-blind manner, which we refer to as CycFusion-noblind.

The CAVE dataset contains 32 indoor HSI with the size of $512\times 512$ and 31 spectral channels covering wavelengths ranging from 0.4 $\mu m$ to 0.7 $\mu m$. The original Chikusei dataset was acquired over Chikusei, Japan, by the Headwall Hyperspec-VNIR-C imaging sensor in 2014. The scene is with the size of $2517\times 2335$ and has 128 bands covering the spectral range from 0.363 $\mu m$ to 1.018 $\mu m$ with the ground sampling distance (GSD) of 2.5m. We select 16 non-overlapped regions with the size of $512\times 512$ from the raw data to evaluate the performance of all comparison methods. The last HSI we used in the experiment is the Pavia University dataset collected by the reflective optics system imaging spectrometer (ROSIS) covering the University of Pavia, Italy. The image contains $610\times 340$ pixels and has 103 bands ranging from 0.43 $\mu m$ to 0.86 $\mu m$ with the GSD of 1.3m. We crop the top-left area with the size of $608\times 320$ for the convenience of data processing. The above data we refer to as the ground truth (GT) of HrHSI to evaluate the fused results. Some benchmark images in the three datasets are depicted in Fig. 2.

Following [30, 33] and [34], we directly average the $32\times 32$ spatially disjoint block in the HrHSI to simulate the observed LrHSI and use the Nikon D700 camera¹²¹²12https://www.maxmax.com/spectral_response.htm, LANDSAT-8¹³¹³13https://landsat.gsfc.nasa.gov/article/preliminary-spectral-response-of-the-operational-land-imager-in-band-band-average-relative-spectral-response and an IKONOS-like SRF to generate the HrMSI via HrHSI for the CAVE, Chikusei and Pavia University dataset, respectively. The HrMSI of the CAVE dataset have 3 bands coresponding to red, green, and blue channels, while HrMSI of other two remote sensing dataset have 4 bands with one more NIR band than the former.

The average performance metrics on the CAVE dataset for all comparison methods are listed in Table III. Overall, GSA and GLP-HS as pansharpening-based methods with parameters of SRF and PSF absent yield poor fusion results. Compared with them, the subspace-based methods combine the observation model and give competitive fusion performance, especially NSSR, bringing 2.8dB improvement in PSNR compared to GLP-HS. GDD attempts to learn the prior information on HrHSI from a noise input with a deep CNN, however, GDD may not learn enough prior to the indoor dataset, resulting in poor results. CUCaNet couples the spectral unmixing model and utilizes the multiple consistency loss to learn the parameters of the observation model and fusion network, which works in a blind fusion manner and brings 0.45dB increment in PSNR compared to GDD. Unsurprisingly, the proposed CycFusion outperforms all other fusion algorithms with performance metrics as listed in Table III. Specifically, CycFusion and CycFusion-noblind increase 2.41dB and 0.26dB in terms of PSNR compared to the second-best blind and no-blind fusion methods, GLP-HS and NSSR, respectively.

Fig. 4 and Fig. 5 show the fused HrHSI in the 21st band (600nm) for balloons and CD, respectively. Two regions of interest (ROIs) are highlighted for comparing detailed differences of all testing algorithms. It can be observed that there are prominent bumps in the fused image provided by the pansharpening-based GLP-HS as shown in Fig. 4 (c) and Fig. 5 (c). Besides, unusual ripples are present in the result of GDD especially in the Fig. 5 (h) that may be caused by insufficient prior learning. Overall, the proposed CycFusion and CycFusion-noblind approaches give better fusion results of the HrHSI and corresponding lower error maps than other comparison methods. To compare the reconstruction in each band obtained by the testing algorithms, the band-by-band root mean squared error (RMSE) of the first five best methods on these two images are shown in Fig. 6. These RMSE results further illustrate the superiority of our proposed methods in spectral reconstruction.

To further evaluate the effectiveness, we conduct two remote sensing datasets. Different from the indoor images, the remote sensing data have lower spatial resolution and generally contain several spectral signatures in one image. Table IV shows the quantitative metrics of all testing methods over the Chikusei and Pavia University datasets. As one can see, the SAM results evaluated on these two datasets are much smaller than the same metrics on the CAVE dataset in Table III. In the no-blind fusion methods, our proposed CycFusion-noblind also obtains the best results among all quantitative metrics. We also depict one band in each fused image for visual comparison as shown in Fig. 7 and Fig. 8. Our proposed model also produces the best fusion results and well spectral reconstruction as shown in Fig. 9 .