# Coupled Convolutional Neural Network With Adaptive Response Function Learning for Unsupervised Hyperspectral Super Resolution

Ke Zheng, Lianru Gao, Wenzhi Liao, Danfeng Hong, Bing Zhang, Ximin Cui, Jocelyn Chanussot

► To cite this version:

Ke Zheng, Lianru Gao, Wenzhi Liao, Danfeng Hong, Bing Zhang, et al. Coupled Convolutional Neural Network With Adaptive Response Function Learning for Unsupervised Hyperspectral Super Resolution. IEEE Transactions on Geoscience and Remote Sensing, 2021, 59 (3), pp.2487-2502. 10.1109/TGRS.2020.3006534. hal-03429693

HAL Id: hal-03429693

https://hal.science/hal-03429693v1

Submitted on 25 Sep 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Coupled Convolutional Neural Network with Adaptive Response Function Learning for Unsupervised Hyperspectral Super-Resolution

Ke Zheng, Lianru Gao, *Senior Member, IEEE*, Wenzhi Liao, *Senior Member, IEEE*, Danfeng Hong, *Member, IEEE*, Bing Zhang, *Fellow, IEEE*, Ximin Cui, and Jocelyn Chanussot, *Fellow, IEEE*

**Abstract**—Due to the limitations of hyperspectral imaging systems, hyperspectral imagery (HSI) often suffers from poor spatial resolution, thus hampering many applications of the imagery. Hyperspectral super-resolution refers to fusing HSI and MSI to generate an image with both high spatial and high spectral resolutions. Recently, several new methods have been proposed to solve this fusion problem, and most of these methods assume that the prior information of the Point Spread Function (PSF) and Spectral Response Function (SRF) are known. However, in practice, this information is often limited or unavailable. In this work, an unsupervised deep learning based fusion method HyCoNet that can solve the problems in HSIMSI fusion without the prior PSF and SRF information is proposed. HyCoNet consists of three coupled autoencoder nets in which the HSI and MSI are unmixed into endmembers and abundances based on the linear unmixing model. Two special convolutional layers are designed to act as a bridge that coordinates with the three autoencoder nets, and the PSF and SRF parameters are learned adaptively in the two convolution layers during the training process. Furthermore, driven by the joint loss function, the proposed method is straightforward and easily implemented in an end-to-end training manner. The experiments performed in the study demonstrate that the proposed method performs well and produces robust results for different datasets and arbitrary PSFs and SRFs.

**Index Terms**—Hyperspectral Image, Super-Resolution, Coupled Convolutional Neural Network, Autoencoder, Adaptive Learning

## I. INTRODUCTION

A Hyperspectral image (HSI) is a data cube containing hundreds of contiguous narrow-bandwidth images covering a large wavelength range [1]. Because of its high spectral resolution, HSI is very important in many applications such as land cover classification [2]–[5], target detection [6]–[9], feature extraction and dimensionality reduction [10]–[13], data fusion [14]–[16], and spectral unmixing [17]–[19]. However, the hyperspectral imaging system often has a trade-off between spectral resolution and spatial resolution, due to hardware restrictions. This causes the spatial resolution of HSI to usually be coarser than that of multispectral imagery (MSI). The limited spatial resolution restricts further applications of HSI. To enhance the spatial resolution of HSI, a natural solution is to fuse it with higher-resolution MSI. This approach is called hyperspectral and multispectral image fusion (HSI-MSI Fusion).

HSI-MSI Fusion is similar to the MSI pansharpening process in which a low spatial resolution MSI is fused with a high-resolution panchromatic (PAN) image. However, applying these pansharpening methods to directly fuse HSI and PAN imagery remains challenging as PAN images contain limited spectral information, and, as a result, spectral distortion can easily occur [20]. Recently, HSI-MSI Fusion has attracted a lot of attention because the result preserves more accurate spectral information with a high-spatial resolution. The existing fusion methods can be categorized as: (1) extensions of pansharpening methods; (2) Bayesian-based approaches; and (3) matrix factorization-based methods [21].

In the first category, Richard et al. first attempted to use a pansharpening-based method to fuse HSI and MSI using a wavelet technique [22]. However, the results were highly dependent on spectral resampling, which made it difficult to enhance the spatial resolution. Zhang et al. proposed a pansharpening-based fusion method that used a 3-D wavelet transform [23]. Chen et al. proposed a framework for fusing HSI and MSI by dividing the HSI into several regions and fusing the HSI and MSI in each region using the pansharpening method [24]. Aiazzi et al. proposed a component substitution fusion method that took the spectral response function (SRF) as part of the model [25]. Liu et al. proposed a spectral

This work was supported by the National Natural Science Foundation of China under Grant No. 41722108 and No. 91638201, as well as with the support of the AXA Research Fund. (Corresponding author: Lianru Gao.)

K. Zheng is with the Key Laboratory of Digital Earth Science, Aerospace Information Research Institute, Chinese Academy of Sciences, Beijing 100094, China, and also with the College of Geoscience and Surveying Engineering, China University of Mining and Technology (Bei Jing), Beijing 100083, China. (e-mail: zhengkevic@aircas.ac.cn)

L. Gao is with the Key Laboratory of Digital Earth Science, Aerospace Information Research Institute, Chinese Academy of Sciences, Beijing 100094, China. (e-mail: gaolr@aircas.ac.cn)

W. Liao is with the Sustainable Materials Management, Flemish Institute for Technological Research (VITO), 2400 Mol, Belgium, and also with the Image Processing and Interpretation, IMEC Research Group, Ghent University, 9000 Ghent, Belgium. (e-mail: wenzhi.liao@vito.be)

D. Hong is with the Remote Sensing Technology Institute (IMF), German Aerospace Center (DLR), 82234 Weling, Germany. (e-mail: danfeng.hong@dlr.de)

B. Zhang is with the Key Laboratory of Digital Earth Science, Aerospace Information Research Institute, Chinese Academy of Sciences, Beijing 100094, China, and also with the College of Resources and Environment, University of Chinese Academy of Sciences, Beijing 100049, China. (e-mail: zb@radi.ac.cn)

X. Cui is with the College of Geoscience and Surveying Engineering, China University of Mining and Technology (Bei Jing), Beijing 100083, China. (e-mail: cxm@cumtb.edu.cn)

J. Chanussot is with the Univ. Grenoble Alpes, CNRS, Grenoble INP, GIPSA-lab, F-38000 Grenoble, France, also with Aerospace Information Research Institute, Chinese Academy of Sciences, Beijing 100094. (e-mail: jocelyn.chanussot@grenoble-inp.fr)
preservation fusion method that was based on a simplified solar radiation and land-surface reflection model [26].

In the second category, Eismann et al. proposed a Bayesian-based fusion method that used a stochastic mixing model of the underlying spectral content to achieve resolution enhancement [27]. Qi et al. proposed a variational-based fusion method that assumed the target image in a low-dimensional subspace and that solved the fusion problem by alternating optimization with respect to the coding coefficients and the target image [28]. Simes et al. formulated the fusion problem as a minimization of a convex objection containing two quadratic terms and an edge-preserving term [29]. Akhtar et al. proposed a non-parametric Bayesian sparse coding strategy which first inferred the probability distributions of the material spectra and then computed the sparse codes of the high-resolution image [30].

Methods in the third category the matrix factorization methods usually assume that the HSI is composed of a series of pure spectral vectors and that the matrix HSI can be decomposed into abundances and endmembers. The fusion problem becomes an estimation problem for the endmembers of the low-resolution HSI and abundances of the high-resolution MSI. Kawakami et al. proposed an unmixing approach to fusing a low-resolution HSI with a high-resolution RGB image. Firstly, the unmixing algorithm was employed to estimate the basis endmembers of the low-resolution HSI; this was then combined with a high-resolution RGB image to produce the final result [31]. Instead of keeping the estimated endmembers of the low-resolution HSI fixed, Yokoya and Lanaras presented a coupled NMF (CNMF) to estimate endmembers and abundances using an alternating unmixing approach [32], [33]. Wycoff et al. restricted sparse regularization to abundances which assume that each pixel is composed of only a small number of endmembers [34]. Akhtar et al. proposed a sparse representation-based approach with a local spatial structure constraint whose main feature was the exploitation of the local patch prior knowledge using a greedy pursuit algorithm [35]. Yi et al. proposed an interactive feedback strategy fusion method with spectral unmixing and spatial constraints [36]. Tensor-based fusion methods utilizing tensor factorization with sparse constraint or subspace projection have also been proposed [14], [37].

In recent years, deep learning has been successfully applied in many spectral tasks [38]–[41]. To deal with the HSIMSI fusion problem with deep learning, Wang et al. proposed a convolutional neural network to fuse HSI and MSI using a residual network and a preserved spectral loss function [42]. Han et al. proposed a partial densely connected network to fuse MSI and HSI spatial and spectral information [43]. Palsson et al. proposed a 3-D convolutional neural network together with reducing the dimensionality of the HSI to make the fusion more computationally efficient [44]. Dian et al. initialized the HSI by solving a Sylvester equation and then implementing a neural network to learn the mapping between the initialized HSI and target HSI [45]. Xie et al. proposed a fusion network that took the observation models of low-resolution images and the low-rank knowledge into consideration [46]. Wang et al. proposed a deep learning-based blind hyperspectral image fusion method with iterative

and alternating optimization strategy [47]. Han et al. presented a multi-scale spatial and spectral fusion convolutional neural network (CNN) for HSI-MSI Fusion [48]. However, all of the deep learning-based methods mentioned above are supervised learning methods that are difficult to apply in practice because the high-resolution HSI needed for training is unavailable. Qu et al. proposed an unsupervised deep learning-based fusion method with a sparse Dirichlet network [49]. Zhou et al. proposed a registration algorithm and a fusion algorithm to handle HSI and MSI image with significant scale difference and nonrigid distortion [50]. Fu et al. proposed a camera spectral response (CSR) optimization layer to learn the spectral response with an unsupervised way [51].

Some of fusion methods assume that the prior information of the SRF or point spread function (PSF) is known. However, in practice, this information is often difficult to obtain [29]. In this paper, an unsupervised deep learning-based fusion network that can handle situations where the PSF and SRF are unknown is proposed. The only information our proposed method requires is the spectral coverage of the MSI and HSI, which is easy to obtain from the data provider. To our knowledge, this is the first time which the unsupervised coupled CNN was developed with learnable PSFs for the HSI-MSI Fusion task. The main contributions of this study can be summarized as follows:

* A novel unsupervised network called HyCoNet is proposed to solve the HSIMSI fusion problem for an unknown SRF and PSF. The results show that the proposed method can deal well with arbitrary SRFs and PSFs in comparison with nine state-of-the-art HSIMSI fusion methods, as applied to four remote sensing datasets;

* Based on the linear unmixing theory, three autoencoder networks are jointly coupled in the proposed method. Within these networks, the endmembers comprise the parameters of one convolution layer, which is shared between two autoencoder networks. Also, in order to improve the consistency of these networks, a learned PSF layer acts as a bridge connecting the low- and high-resolution abundances;

* A joint-loss function that drives the model using unsupervised learning and an end-to-end training manner, thus providing a simple and direct training strategy for obtaining the final result, is introduced.

This paper is organized as follows. Section II describes the basic formulation of the HSI and MSI fusion problem. Section III introduces the proposed fusion model, including the network architecture and the joint-loss function. Section IV presents the experimental results and discussion; Section V is the conclusion.

## II. PROBLEM FORMULATION

The HSI-MSI fusion problem requires the estimation of the HSI, which has both high spectral and high spatial resolution and is denoted as $X \in \mathbb{R}^{M \times N \times L}$, where $M$, $N$ and $L$ are the width, height, and number of spectral bands. The input images include a high spatial resolution MSI denoted as $Y \in \mathbb{R}^{M \times N \times l}$ and a low spatial resolution HSI denoted as $Z \in \mathbb{R}^{M \times N \times l}$.
<figure><img src="image_2.png" /><figcaption>Fig. 1. The illustration of the relationship between the HrMSI, the LrHSI and the target HrHSI based on the linear unmixing model, where $Y_{lr}$ represents low spatial resolution multispectral image (LrMSI).</figcaption></figure>

$\mathbb{R}^{m \times n \times L}$, where $l$ is the number of spectral bands in $Y$, and $m$ and $n$ are the width and height of $Z$. In particular, $n \le N$, $m \le M$ and $l \le L$. For convenience, we call the low spatial resolution HSI, high spatial resolution MSI and target HSI and low spatial resolution MSI as LrHSI, HrMSI and HrHSI, respectively. To simplify the notation in this section, we unfold the 3D data cube to form a 2D matrix: the 2D matrices of the three types of image are denoted as $X \in \mathbb{R}^{MN \times L}$, $Y \in \mathbb{R}^{MN \times l}$ and $Z \in \mathbb{R}^{mn \times L}$, respectively.

The relationship between these images $X$, $Y$ and $Z$ according to the hyperspectral linear unmixing model [52] is shown in Fig. 1. In the linear spectral unmixing model, each pixel of the HSI is assumed to be a linear combination of a set of pure spectral bases (called endmembers) and the coefficients of each pure spectral basis (called abundances). According to the linear spectral unmixing model, the target HrHSI can be described as:

$$X = AE \quad (1)$$

where the matrix $A \in \mathbb{R}^{MN \times p}$ is formed from the abundances, the matrix $E \in \mathbb{R}^{p \times L}$ is made up of the endmembers, and $p$ is the number of pure spectral bases. This equation describes the degree of mixing for each pixel in the image $X$.

Similarly, the input $Z$ can also be expressed as a linear combination of the same endmembers $E$:

$$Z = A_h E = S * X = S * AE \quad (2)$$

where the matrix $A_h \in \mathbb{R}^{mn \times p}$ represents the abundances, the matrix $S \in \mathbb{R}^{L_m \times L_n}$ is the point spread function (PSF), which describes the spatial degradation function, and $*$ denotes the convolution operator [53]. $L_m$ and $L_n$ are the spatial size of the convolution filter. $Z$ also denotes the spatially degraded version of image $X$. The input $Y$ is the spectrally degraded version of $X$:

$$Y = XR = AER \quad (3)$$

where the matrix $R \in \mathbb{R}^{L \times l}$ is the spectral response function (SRF), which describes the spectral degradation process.

Moreover, the spectral degraded version of $Z$ should approximate to the spatially degraded version of $Y$:

$$Y_{lr} = ZR = S * Y \quad (4)$$

where the matrix $Y_{lr} \in \mathbb{R}^{m \times n \times l}$ represents the low spatial resolution multispectral image (LrMSI).

In addition, our goal is to estimate $X$ using the inputs $Y$ and $Z$ with the following constraints also satisfied:

$$\begin{aligned} \sum_{j=1}^{p} a_{ij} &= 1 \quad \forall i, j \\ a_{ij} &\ge 0 \quad \forall i, j \\ 1 &\ge e_{ij} \ge 0 \quad \forall i, j \end{aligned} \quad (5)$$

where $a_{ij}$ is a component unit of $A$ and $e_{ij}$ is a component unit of $E$. These constraints relate to the sum-to-one property of the abundance, the non-negative property of the abundance and the bounded non-negative property of the endmembers, respectively [28]. In addition, the abundances should be sparse, meaning that each HSI pixel is composed of only a few pure spectral bases.

### III. PROPOSED METHOD

According to Eqs. (1)-(3), to solve the HSIMSI fusion problem, the key point is to estimate the high spatial resolution abundance matrix $A$ and the spectral bases matrix $E$. The HrMSI provides detailed spatial contextual information that is highly correlated with $A$. Also, the LrHSI preserves the spectral information, which is highly consistent with the target spectral endmembers matrix $E$. The basic idea of the proposed method is based on matrix factorization. The proposed Hy-CoNet is an unsupervised network that includes three coupled autoencoder networks. The target HrHSI is embedded in one of the networks this will be elaborated on in the part A of this section, and the Part B introduces the joint objective function used in the training.

#### A. Coupled Autoencoder Network for Image Fusion

The proposed network is composed of three autoencoder nets, as shown in Fig. 2. The upper one called LrHSI autoencoder, and the lower one called HrMSI autoencoder. The estimated target HrHSI $X$ is embedded in the HrMSI autoencoder. The LrMSI autoencoder can be seen at the upper right.

Since the 2D convolution layer, which has a kernel size of $1 \times 1$, is equivalent to the fully connected layer when it applied to a spectral vector, all the fully connected layers in the traditional autoencoder network are replaced by the convolution layers, as shown in Fig. 2. In addition, all the convolution kernel sizes are set to be $1 \times 1$, except for the PSF convolution layer. Instead of using the fully connected layer, the convolution layer is used to preserve the spatial structure of the input image cube for easy implementation of the PSF operation. Further details of the PSF operation are discussed below.

In the LrHSI autoencoder, the network tries to learn an approximation to the identity function $f(Z) \approx Z$. Since the input LrHSI contains sufficient spectral information, the
<figure><img src="image_3.png" /><figcaption>Fig. 2. The generator networks of proposed HyCoNet method which contains three coupled autoencoder nets: LrHSI-Autoencoder, HrMSI-Autoencoder and LrMSI autoencoder.</figcaption></figure>

endmembers **E** and abundances **A**<sub>h</sub> can be extracted during the reconstruction process for this autoencoder.

The module before the latent variable **A**<sub>h</sub><sup>a</sup> is the encoder:

$$ \mathbf{A}_h^a = f_{en}(\mathbf{Z}) \qquad (6) $$

where the $f_{en}()$ try to learn a nonlinear mapping which transforms the input LrHSI to its abundances $\mathbf{A}_h^a$; $\mathbf{A}_h^a \in \mathbb{R}^{m \times n \times p}$ represents the abundances of **Z**, and *p* is the number of endmembers.

The module after the latent variable is the decoder function:

$$ \tilde{\mathbf{Z}}_a = f_{de}(\mathbf{A}_h^a) \qquad (7) $$

where $\tilde{\mathbf{Z}}_a$ is the output of this upper autoencoder, which represents the reconstructed input image cube. The decoder function $f_{de}()$ is a convolution layer without bias and is shown as the green convolution layer 1x1 Conv in Fig. 2. This layer is also used in the HrMSI autoencoder, meaning that these two autoencoders share the same parameters as this convolution layer. The parameters of the shared convolution layer are denoted as the endmembers **E** ∈ $\mathbb{R}^{p \times 1 \times 1 \times L}$, where 1 × 1 represents the spatial size of this convolution kernel. The size of each convolution kernel is *p* × 1 × 1 and the number of kernels is *L*.

The reason, that we call the latent variables $\mathbf{A}_h^a$ and parameters **E** as the abundances and endmembers, is that they are trying to reconstruct **Z** by matrix multiplication: $\tilde{\mathbf{Z}}_a = \mathbf{A}_h^a \mathbf{E}$. This is the same idea as matrix decomposition except that the process is optimized by gradient descent in neural network. In order to satisfy the non-negativity restrictions for $\mathbf{A}_h$ and **E**, several tricks are applied, as described in more detail in Part III-B.

The structure of the HrMSI autoencoder is similar to that

of the LrHSI autoencoder and includes an encoder function $h_{en}()$ and a decoder function $h_{de}()$. The encoder function can be expressed as:

$$ \mathbf{A} = h_{en}(\mathbf{Y}) \qquad (8) $$

where **A** is the high-resolution abundance, $h_{en}()$ is the encoder function and **Y** is the input HrMSI.

The decoder for the HrMSI autoencoder consists of two parts a shared convolution layer and the SRF:

$$ \tilde{\mathbf{Y}} = h_{de}(\mathbf{A}) = \mathrm{SRF}(f_{de}(\mathbf{A})) \qquad (9) $$

where $\tilde{\mathbf{Y}}$ is the reconstructed HrMSI **Y**, $h_{de}()$ is the decoder function, $f_{de}()$ is the shared convolution layer containing the parameters of the endmember matrix **E**, and $\mathrm{SRF}()$ is the spectral resampling operation. The HrHSI is the output of the shared convolution layer:

$$ \tilde{\mathbf{X}} = f_{de}(\mathbf{A}) \qquad (10) $$

where $\tilde{\mathbf{X}}$ is the estimated target image.

To handle the situation where the SRF parameters are unknown, a convolution layer and a normalization layer are placed after the target $\tilde{\mathbf{X}}$ to learn the unknown parameters of the SRF. The SRF consists of spectral resampling from HSI to MSI and the process can be defined as:

$$ \varphi_i = \frac{\int_{\lambda_{i,L}}^{\lambda_{i,U}} R(\lambda)\varepsilon(\lambda)d\lambda}{\int_{\lambda_{i,L}}^{\lambda_{i,U}} R(\lambda)d\lambda} \qquad (11) $$

where $\varphi_i$ is the spectral radiance of band *i* of the HrMSI, $\lambda$ is the wavelength, $\lambda_{i,U}$ and $\lambda_{i,L}$ are the wavelength bounds of band *i* of the HrMSI, $R()$ is the spectral response function, and $\varepsilon$ is the spectral radiance of the HrHSI. To implement the spectral resampling in the neural network, a convolution layer
with kernel size 1 × 1 (shown in red and labeled 1x1Conv in Fig. 2) is added after the target image $\tilde{\mathbf{X}}$ to simulate the numerator of Eq. 11. A normalization layer (labeled Norm in Fig. 2) follows the convolution layer and simulates the denominator of Eq. 11.

Therefore, the SRF process within our network can be expressed as:

$$ \varphi_i = \text{SRF}(\varepsilon_\lambda) = \frac{\sum_{\lambda=\lambda_{i,L}}^{\lambda_{i,U}} w_{i,\lambda} \varepsilon_\lambda}{\sum_{\lambda=\lambda_{i,L}}^{\lambda_{i,U}} w_{i,\lambda}} \quad (12) $$

where $\varphi_i$ is the band i image in $\tilde{\mathbf{Y}}$, $w_{i,\lambda}$ is the weight of the SRF convolution layer, and $\varepsilon_\lambda$ is the band with wavelength $\lambda$ in $\tilde{\mathbf{X}}$. The use of this function means that the convolution layer and the normalization layer integrate the HrHSI $\tilde{\mathbf{X}}$ with the weights of the convolution layer between the upper and lower spectral bounds. It also means that our network assumes the spectral coverage of HrMSI bands is known and that each convolution kernel of the convolution layer only covers the spectral range corresponding to the bands in $\tilde{\mathbf{X}}$. The number of convolution kernels is equals to the number of bands in $\tilde{\mathbf{Y}}$.

The PSF means that a given pixel is a weighted combination of contributions from the pixel and its neighboring pixels [54], [55]. In the fusion problem, this means that a pixel in the LrHSI is a weighted combination of local pixels from the HrHSI. Therefore, the PSF is a convolution process and a convolution operation can easily be implemented as part of the convolutional neural network. In our network, to simulate the PSF, a convolution layer with 1 input channel and 1 output channel is implemented for every band of the abundance A. According to the definition of the PSF, the relationship between X and Z can be expressed as $Z = PSF(X)$. Also, X and Z are composed of the same linear unmixing endmembers. Therefore, another low-resolution abundance $A_h^b$ can be characterized as:

$$ A_h^b = PSF(\mathbf{A}) \quad (13) $$

where PSF() is the convolution layer labeled with the blue arrow and 1LayerConv in Fig. 2. The spatial size of this PSF convolution kernel is the same as the ratio of the Ground Sampling Distance (GSD) of the LrHSI to that of the HrMSI, and the stride of the convolution layer equals the kernel size. The PSF process acts as a bridge between the LrHSI autoencoder and the HrMSI autoencoder and forces the reconstructed image to be spectrally consistent $\tilde{\mathbf{X}}$. Therefore, another LrHSI $\tilde{\mathbf{Z}}_b$ can be reconstructed using $A_h^b$ and E:

$$ \tilde{\mathbf{Z}}_b = f_{de}(A_h^b) \quad (14) $$

In addition, the spatial degraded version of the HrMSI is equivalent to the spectrally degraded version of the LrHSI shown as LrMSI- Autoencoder in Fig. 2. This relation can be expressed as:

$$ PSF(\mathbf{Y}) = \tilde{\mathbf{Y}}_{lr}^a \approx \tilde{\mathbf{Y}}_{lr}^b = SRF(\mathbf{Z}) \quad (15) $$

where $\tilde{\mathbf{Y}}_{lr}$ is the estimated LrMSI.

<figure><img src="image_4.png" /><figcaption>Fig. 3. The illustration of flow charts for the inputs Y, Z and estimated target $\tilde{\mathbf{X}}$.</figcaption></figure>

## B. Joint Loss Function

Fig. 3 shows flow charts for the proposed method, including input images, output images and corresponding target images. The first row of Fig. 3 indicates the inference of the input Z in the LrHSI-Autoencoder and is also represented by the yellow arrows in Fig. 2. The inference of HrMSI Y includes two parts. The first one is shown as the second row of Fig. 3, which includes the PSF operation to generate the LrHSI this process is represented by the blue arrows Fig. 2. The second part can be seen in the third row of Fig. 3 and consists of the SRF operation; this process is represented by the red arrows in Fig. 2. The last row of Fig. 3 shows the relationship between the two LrMSIs. Therefore, the objective function for reconstruction can be expressed as:

$$ L_{base}(\mathbf{Y}, \mathbf{Z}) = \left\| \mathbf{Z} - \tilde{\mathbf{Z}}_{a} \right\|_{1} + \alpha \left\| \mathbf{Z} - \tilde{\mathbf{Z}}_{b} \right\|_{1} + \beta \left\| \mathbf{Y} - \tilde{\mathbf{Y}} \right\|_{1} + \gamma \left\| \tilde{\mathbf{Y}}_{lr}^{a} - \tilde{\mathbf{Y}}_{lr}^{b} \right\|_{1} \quad (16) $$

where $\alpha$, $\beta$ and $\gamma$ are trade-off parameters that tune the weights between these reconstruction errors.

The sum-to-one and non-negative properties given in Eq. 5 also need to be satisfied. First, we constrain the sum of the abundances in the channel dimension to meet the sum-to-one property:

$$ L_{sum2one}(\mathbf{Y}, \mathbf{Z}) = \left\| 1 - \sum_{i=1}^{p} A_i \right\|_{1} + \left\| 1 - \sum_{i=1}^{p} A_{h,i}^{a} \right\|_{1} + \left\| 1 - \sum_{i=1}^{p} A_{h,i}^{b} \right\|_{1} \quad (17) $$

where i indicates the ith band of the abundance matrix A. Although the softmax function can be used to strictly enforce the sum-to-one property for the abundances, the resulting convergence accuracy is lower than for the proposed method. Part IV-D will explain this in detail.

Secondly, to enforce the non-negative property, several tricks are applied during training. The clamp function is applied to the output of the last convolution layer of both the encoder nets and decoder nets to force all the elements of the abundances and reconstruction images into the range [0,1]. Although we tried to use the sigmoid activation layer for this purpose, we found that it was difficult to make the network converge using this function. In addition, the weights of the shared convolution layer (containing the endmember parameter matrix, E), PSF layer and SRF layer should also meet the non-negative property. Since the weights of these
<table><caption>TABLE I
THE USED DATA IN THE EXPERIMENT.</caption><thead><tr><th></th><th>Pavia University</th><th>Indian Pines</th><th>Washington DC</th><th>University of Houston</th></tr></thead><tbody><tr><td>Spatial size of HrHSI</td><td>336x336</td><td>144x144</td><td>304x304</td><td>320x320</td></tr><tr><td>Spectral range of HrHSI</td><td>466-834nm</td><td>400-2500nm</td><td>400-2500nm</td><td>403-1047nm</td></tr><tr><td>Number bands of HrHSI</td><td>103</td><td>191</td><td>191</td><td>46</td></tr><tr><td>GSD ratio</td><td>4</td><td>4</td><td>8</td><td>8</td></tr><tr><td>Spatial size of LrHSI</td><td>84x84</td><td>36x36</td><td>38x38</td><td>40x40</td></tr><tr><td>Bands of HrMSI</td><td>Blue-Green-Red</td><td>Blue to SWIR2</td><td>Blue to SWIR2</td><td>Blue-Green-Red</td></tr><tr><td>Number bands of HrMSI</td><td>3</td><td>6</td><td>6</td><td>3</td></tr></tbody></table>

layers may be updated to a negative value after the back-propagation, we applied the clamp function to these layers after the weights were updated to force these weights into the non-negative range. The range of the clamp function was [0, 1]. As a result, the weights satisfied the non-negative and bounding constraints for each forward propagation, except for the first time.

Since each pixel of the HSI is composed of a small number of pure spectral bases, the abundance matrix **A** should be sparse. To guarantee the sparsity of the abundance, the KL divergence is used to ensure that most of the elements in the abundance are close to a small number:

$$ 
\begin{aligned}
L_{sparse}(\mathbf{Y}, \mathbf{Z}) &= \sum_{i=1}^{s} \sum_{j=1}^{p} KL(a \| \tilde{a}_{i,j}) \\
&= \sum_{i=1}^{s} \sum_{j=1}^{p} \left( a \log\left(\frac{a}{\tilde{a}_{i,j}}\right) + (1-a) \log\left(\frac{1-a}{1-\tilde{a}_{i,j}}\right) \right)
\end{aligned} 
\quad (18) 
$$

where *s* is the number of pixels, *p* is the number of convolution kernels and also the number of endmembers, *a* is a sparsity parameter which is set to a small value close to zero (0.0001 in our network), and $\tilde{a}_{i,j}$ is an element of the abundance matrix **A**. To satisfy this constraint, the elements of **A** must mostly be near zero. The sparsity constraint is also applied to the abundance matrix $\mathbf{A}_h^a$.

Ultimately, our aim is to solve the fusion problem in accordance with the optimization problem:

$$ L(\mathbf{Y}, \mathbf{Z}) = L_{base}(\mathbf{Y}, \mathbf{Z}) + \mu L_{sum2one}(\mathbf{Y}, \mathbf{Z}) + \nu L_{sparse}(\mathbf{Y}, \mathbf{Z}) \quad (19) $$

where $\mu$ and $\nu$ are the trade-off parameters used to balance the errors. This loss function can be directly used in the optimizer, thus providing a simple and direct solution.

## IV. EXPERIMENTS

To get an accurate assessment of the fusion quality and evaluate the performances of different fusion methods, simulation experiments [56] are used in the experiment. The proposed HyCoNet was implemented using four different simulated data sets. Firstly, the sensitivity of the trade-off parameters $\alpha, \beta, \gamma, \mu$ and $\nu$ was evaluated. Secondly, the constraints on the abundances were investigated. Thirdly, we compared the effectiveness of the method for different numbers of endmembers. Fourthly, the learned PSF kernels were investigated using different spatially down-sampled kernels. The character of the

estimated abundances was then explored. Finally, the fused images obtained using the different methods were evaluated using both visual and quantitative measures.

### A. Experimental Dataset

The proposed HyCoNet was evaluated using four widely used HSI datasets: Pavia University, Indian Pines, Washington DC, and University of Houston. The Pavia University data were acquired by the ROSIS-3 optical airborne sensor in 2003. This image consists of 610 × 340 pixels with a GSD of 1.3m and spectral range of 430 nm840 nm in 115 bands. Due to the effects of noise and water vapor absorption, 12 bands have been removed. An area covering 366 × 366 pixels in the lower-left corner of the image and containing 103 bands was selected for use in this experiment. The Indian Pines data were acquired by the AVIRIS sensor in 1992. This image consists of 145 × 145 pixels with a 20m GSD; the spectral range is 400 nm2500 nm covering 224 bands. After removing 33 noisy bands, we selected a 144 × 144-pixel area with 191 bands as experimental data. The Washington DC data were acquired by the HYDICE sensor in 1995. This image has an area of 1280 × 307 pixels and a GSD of 2.5m. The spectral range is 400 nm2500 nm consisting of 210 bands. After removing 19 noisy bands, we selected 191 bands covering 304 × 304 pixels for use. The University of Houston data were used in the 2018 IEEE GRSS Data Fusion Contest [57], and consist of 601 × 2384 pixels with a 1-m GSD. The data covers the spectral range 380 nm1050nm with 48 bands. We selected 46 bands consisting of 320 × 320 pixels from this imagery for use as experimental data.

### B. Implementation Details

In the experiment, simulated spatial downsampling was used to generate the LrHSI using a Gaussian filter, as is widely used in remote sensing [58]. In the experiment, the width and height of the Gaussian filter was set equal to the ratio between the high-resolution GSD and the low-resolution GSD. The standard deviation of all the Gaussian filters was set to 0.5, except as described in part IV-F and IV-H where different standard deviations were used to evaluate the robustness of the fusion model. To simulate the HrMSI, the SRF for the blue to SWIR2 bands of the Landsat 8 were used [59]. To verify the stability of the model, different GSD ratios and number bands were used to simulate the LrHSI and HrMSI. For the Pavia University and Indian Pines data, the GSD ratio was set
<figure><img src="image_5.png" /><figcaption>Fig. 4. Parameters discussion. (a): Fusion accuracy of using different α and β in Pavia University data, where α and β mean the weights of reconstruction errors for LrHSI and HrMSI, respectively; (b): PSNR results of the different γ in Pavia University data, where γ represents the weight of LrMSI reconstruction error; (c): Accuracy comparison between μ and ν in the Pavia University data. μ and ν are the weights of sum-to-one error and sparsity error.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 5. Ablation study for the proposed network. ≪ means removing part of the network/losses and studying its performances.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 6. Convergence curve when using different constrained functions for abundances. (a) Convergence accuracy for comparing clamp function and softmax function, respectively. (b) Function curves of clamp and softmax function for a vector evenly spaced 0.01 over the range [-5, 5].</figcaption></figure>

to 4; it was set to 8 for the Washington DC and University of Houston imagery. The bluegreenred bands of the Landsat 8 SRF were used for Pavia University and University of Houston data, and the blue to SWIR2 part of the Landsat 8 SRF was

used for Indian Pines and Washington DC. Table I summarizes the simulated parameters for all of the datasets used in this experiment.

To evaluate the performances of different fusion methods, simulation experiments are used in this experiment. Simulation experiments refer to that the spatial and spectral down-sampling are implemented on the original HrHSI, and this one is the truth target image to evaluate the performance of the estimated HrHSI. Five different quality measures were used to evaluate the performance of the fusion results: the root mean square error (RMSE), peak SNR (PSNR), spectral angle mapper (SAM), relative global dimension error (ERGAS), and mean relative absolute error (MRAE) [56], [60]. Of these measures, the RMSE, MRAE, and SAM were used in the visual evaluation, and the PSNR, SAM, and ERGAS were used in the quantitative evaluation.

The proposed model was trained using an Adam optimizer [61] with the default parameters $β_1 = 0.9$, $β_2 = 0.999$ and $ε = 10^{-8}$; the initial learning rate was set to $5×10^{-3}$. Learning rate schedules seek to adjust the learning rate during training by reducing the learning rate according to a linear decay. After a total of 10000 epochs, the learning rate drops to 0. In our experiment, the number of input images was one; therefore, 1 epoch was equal to 1 iteration and the batch size was also 1. The Pytorch deep learning framework was used to train the proposed model [62]. The training environment consisted of an Intel i7-6850K CPU, 128-GB RAM, and 4×NVIDIA TITAN Xp 12G GPU.

### C. Parameters Discussion

In the proposed method, the parameters α, β, γ, μ and ν in the loss function Eqs. 16 and 19 need to be set. The parameters α, β and γ are used to balance the weights of the different reconstruction errors. The parameters μ and ν are the trade-off parameters for the sum-to-one loss and sparsity loss.

Due to the fact that, in our method, the HSIMSI fusion is driven by the autoencoder network, the reconstruction errors are the main factor that influence this process. Therefore, firstly, we set the parameters γ, μ and ν to a fixed value
<figure><img src="image_8.png" /><figcaption>Fig. 7. The effect of endmember number on the fusion quality for different datasets. (a) Pavia University, (b) Indian Pines, (c) Washington DC, (d) University of Houston.</figcaption></figure>

of 1 to evaluate the effect of using different values of the parameters $\alpha$ and $\beta$. In this part of the experiment, we used the Pavia University data to investigate the performance. Fig. 4(a) shows the PSNR of the target image $X$ for different values of the trade-off parameters $\alpha$ and $\beta$. It can be seen that both parameters have a considerable effect on the fusion quality with the results being more sensitive to $\beta$ than to $\alpha$. This is because $\beta$ controls the weight of the reconstructed HrMSI, and therefore affects the spatial quality of the fused image. $\alpha$ also has an important effect on the quality of the fused target image. It is the weight of the $\tilde{Z}_b$ reconstruction error and the process of reconstructing image $\tilde{Z}_b$ is constrained by both the LrHSI and HrMSI autoencoder networks. Overall, the target image has a high PSNR when the two parameters $\alpha$ and $\beta$ are equal: these PSNR values have been marked in Fig. 4(a). Clearly the PSNR increases as these two parameters increase and tends to be smoother when $\alpha = \beta = 1$. The best result is achieved for $\alpha = \beta = 10$. Therefore, we set $\alpha$ and $\beta$ to 10 in the later experiments.

In the second experiment, we tested the effect of parameter $\gamma$, which controls the weight for the reconstruction of the LrMSI images. Fig. 4(b) shows the experiment results. It can be seen that $\gamma$ has only a slight effect on the fusion quality. Since the reconstruction accuracy is higher and more stable when $\gamma = 100$, we set $\gamma$ to 100 in the later experiments.

In the third experiment, the effects of $\mu$ and $\nu$ were investigated and the results are shown in Fig. 4(c): $\mu$ is the weight of the sum-to-one loss and $\nu$ is the weight of the sparsity loss. It can clearly be seen that the performance is sensitive to $\mu$ as this parameter affects the precision of the fusion reconstruction. In contrast, the performance is not sensitive to $\nu$ since we set the number of endmembers $p$ to be 100. This is because the larger the number of endmembers, with the sum-to-one property, the more likely the endmembers are to be sparse. More details about the experiment carried out using different numbers of endmembers will be given in part IV-E. Accordingly, we set $\mu$ and $\nu$ to 0.001 in the subsequent experiments.

To investigate the essentiality of the proposed network, as shown in Fig. 5, ablation study was implemented in the case of the hyperparameters setting mentioned above. In this experiment, we can clearly see the performances when missing certain parts of the network or losses. As shown in Fig. 5, $\zeta$ means removing certain parts of the network or losses. It can be seen that the performances suddenly drop when removing $\tilde{Z}_b$. Compared to this branch, removing one of the other branches interacts only with less effect. This indicates the branch of $\tilde{Z}_b$ strongly ablation affects the overall fusion performance. Moreover, the branch $\tilde{Z}_b$ means that the core advantage of the learnable PSF layer play an important role in improving the fusion performance.

### D. Constraint Function for the Abundance

Eq. (17) indicates that the sum-to-one property is restricted in the loss function; furthermore, a clamp function is used before the abundance to produce non-negative abundances. [63] reported that combing a softmax function before the abundances can also perfectly restrict the property of sum-to-one and nonnegative. However, in contrast to the clamp function, we found that using the softmax function leads to a lower reconstruction accuracy. Fig. 6(a) shows the fusion accuracy obtained when using different constraint functions for the abundances. It can be observed that the fusion accuracy obviously improves when the clamp function is used. Although the softmax function has a smoother curve, the convergence accuracy and convergence speed are slightly lower.

Fig. 6(b) can explain the results of this experiment: the figures shows the curves for the clamp and softmax functions. It can be seen that, due to the characteristics of the softmax function, it does not converge uniformly, meaning that different points converge at different rates and may converge arbitrarily slowly. This causes the gradient to be smaller for points with smaller values, which is equivalent to the gradient vanishing. In contrast, the gradient of the clamp function is updated faster in the range [0, 1].
<figure><img src="image_9.png" /><figcaption>Fig. 8. The visualization of original and estimated PSF kernel. (a) Pavia University, (b) Washington DC.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 9. The characters of learned abundances of three estimated abundances (Pavia University as example). (a) histograms, (b) the heatmaps of errors for the sum of abundance.</figcaption></figure>

## E. Number of Endmembers

In general, the reconstruction accuracy improves as the number of endmembers, *p*, increases. Fig. 7 shows the changes in the PSNR and SAM with the number of endmembers for the four experimental datasets. To fully explore the effect of the hyperparameter *p*, the experiment was repeated three times for each value of *p* using the same environment. The average values of the PSNR and SAM are shown in Fig. 7. In our model, the number of endmembers *p* represents the number of

feature size of abundance and also represents the kernel size of the shared convolution layer. Therefore, a larger number of endmembers allows the model to be more representative. Although the number of endmembers is assumed to be equal to the number of pure spectral bases in the linear unmixing, the number of endmembers can also be larger than the actual number of pure bases because the convolution weight matrix *E* can contain mixed material [32]. In addition, the convergence accuracy depends on the image complexity. In the experiments using the Indian Pines and Houston University data, the reconstruction accuracy began to converge at *p* = 30, and fast convergence also occurred with the Washington DC data. For the Pavia University data, although the convergence was slow, the results were acceptable for smaller values of *p*. For convenience, we set *p* = 100 in all cases when exploring the performance of the proposed model.

## F. Learned PSF Kernel

The LRHSI was simulated using a PSF with a filter corresponding to the Gaussian function equal to the ratio of the GSD resolutions. Using the Pavia University and Washington DC data, we evaluated the learned kernels under different conditions. For each dataset, the standard deviations of the Gaussian kernel were set to 0.5, 1, and 2.

The original PSF kernels and the estimated ones are shown in Fig. 8. Since the GSD ratios for the Pavia University and Washington DC data were 4 and 8, respectively, the kernel sizes for the two datasets were also 4 and 8, respectively. For Fig. 8, it can be seen that the estimated kernels are similar to the original ones. A large standard deviation produces a smoother kernel and a large weighted combination of contributions from the local pixel. From the above experimental results, it can be seen that the proposed method is highly suitable for estimating arbitrary PSFs for different datasets.

## G. Estimated Abundances

In our network, the abundance is constrained by the sparse and sum-to-one characteristics of the loss function. Therefore, we next used the Pavia University data to explore the three
<figure><img src="image_11.png" /><figcaption>Fig. 10. The visualization of Pavia University fusion results. The first column is the color-composite of fusion results; the second column is RMSE error of color-composite image; the third column is the MARE error of HSI cube; the fourth column is the SAM error of HSI cube.</figcaption></figure>

estimated abundances corresponding to $A_h^a$, $A_h^b$ and $A$. The histograms of these abundances are shown in Fig. 9(a). The figure clearly shows that the estimated abundances are sparse. Fig. 9(b) shows heatmaps of the errors relative to the sum of the abundances. The heatmap for $A_h^a$ shows that some of the areas at the edges did not completely satisfy the sum-to-one

<figure><img src="image_12.png" /><figcaption>Fig. 11. The visualization of Indian Pines fusion results. The first column is the color-composite of fusion results; the second column is RMSE error of color-composite image; the third column is the MARE error of HSI cube; the fourth column is the SAM error of HSI cube.</figcaption></figure>

constraint. The reason for this is that although the sum-to-one character is constrained by the loss function, it is possible that errors remain for some pixels. In contrast, the results for $A_h^b$ and $A$ are much better because they are constrained simultaneously by the LrHSI and HrMSI autoencoders. The abundances in our network, therefore, mostly do have sparse
<figure><img src="image_13.png" /><figcaption>Fig. 12. The visualization of Washington DC fusion results. The first column is the color-composite of fusion results; the second column is RMSE error of color-composite image; the third column is the MARE error of HSI cube; the fourth column is the SAM error of HSI cube.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 13. The visualization of Houston University fusion results. The first column is the color-composite of fusion results; the second column is RMSE error of color-composite image; the third column is the MARE error of HSI cube; the fourth column is the SAM error of HSI cube.</figcaption></figure>

and sum-to-one characteristics.

## H. Comparison with the State of the Art

1) Visual Performance: Following the work of Yokoya [56], in this study, a set of baseline methods were used for comparison. These included CNMF [32], GSOMP [35],

FUSE [64], GLPHS [65], GSA [25], HySure [29], Lanarass method (for convenience, we called it ICCV¹⁵ because it was published in proceedings of the 2015 International Convention on Computer Vision) [33], MAPSMM [27], SFIM-HS [26], and uSDN [49]. Because of the proposed fusion model is unsupervised algorithm and there is shortage of training sam-
<table><caption>TABLE II
QUANTITATIVE PERFORMANCE COMPARISON WITH THE DIFFERENT ALGORITHMS ON THE PAVIA UNIVERSITY DATA. THE BEST ONE IS SHOWN IN BOLD.</caption><thead><tr><th></th><th></th><th>CNMF</th><th>GSOMP</th><th>FUSE</th><th>GLPHS</th><th>GSA</th><th>HySure</th><th>ICCV15</th><th>MAPSMM</th><th>SFIMHS</th><th>uSDN</th><th>HyCoNet</th></tr></thead><tbody><tr><td rowspan="3">σ = 0.5</td><td>mSAM</td><td>4.8715</td><td>10.5319</td><td>4.7626</td><td>5.4793</td><td>4.1585</td><td>4.0756</td><td>4.4325</td><td>5.1705</td><td>5.4588</td><td>5.5614</td><td><strong>3.4107</strong></td></tr><tr><td>mPSNR</td><td>31.8724</td><td>29.0635</td><td>27.7470</td><td>30.4492</td><td>34.3943</td><td>34.5460</td><td>34.7154</td><td>30.8779</td><td>23.2267</td><td>33.6010</td><td><strong>38.7647</strong></td></tr><tr><td>ERGAS</td><td>4.3905</td><td>8.7599</td><td>6.4135</td><td>4.6189</td><td>3.1379</td><td>3.3403</td><td>3.2346</td><td>4.4652</td><td>25.2178</td><td>4.2374</td><td><strong>2.8285</strong></td></tr><tr><td rowspan="3">σ = 1</td><td>mSAM</td><td>4.4029</td><td>9.0941</td><td>4.4667</td><td>4.7169</td><td>3.4810</td><td>3.5015</td><td>4.0694</td><td>4.0270</td><td>4.8534</td><td>5.3427</td><td><strong>3.4002</strong></td></tr><tr><td>mPSNR</td><td>34.8032</td><td>31.0354</td><td>28.1708</td><td>32.9595</td><td>37.5222</td><td>36.9652</td><td>36.2953</td><td>34.9411</td><td>24.9070</td><td>33.1369</td><td><strong>38.9132</strong></td></tr><tr><td>ERGAS</td><td>3.5658</td><td>9.2775</td><td>6.1141</td><td>3.7639</td><td><strong>2.6955</strong></td><td>2.8068</td><td>2.8824</td><td>3.2969</td><td>20.8983</td><td>4.0075</td><td>2.7435</td></tr><tr><td rowspan="3">σ = 2</td><td>mSAM</td><td>3.7239</td><td>8.4744</td><td>4.4371</td><td>4.3608</td><td>3.4896</td><td>3.6169</td><td>4.0755</td><td>4.0317</td><td>4.5768</td><td>5.8754</td><td><strong>3.4286</strong></td></tr><tr><td>mPSNR</td><td>36.0128</td><td>32.0907</td><td>28.2257</td><td>34.6668</td><td>38.4449</td><td>37.5251</td><td>36.7325</td><td>37.4392</td><td>27.3538</td><td>32.9005</td><td><strong>38.5481</strong></td></tr><tr><td>ERGAS</td><td>3.1866</td><td>8.4769</td><td>6.0774</td><td>3.4117</td><td><strong>2.6507</strong></td><td>2.9766</td><td>2.8770</td><td>3.1937</td><td>16.5452</td><td>4.1957</td><td>2.8722</td></tr></tbody></table>

<table><caption>TABLE III
QUANTITATIVE PERFORMANCE COMPARISON WITH THE DIFFERENT ALGORITHMS ON THE INDIAN PINES DATA. THE BEST ONE IS SHOWN IN BOLD.</caption><thead><tr><th></th><th></th><th>CNMF</th><th>GSOMP</th><th>FUSE</th><th>GLPHS</th><th>GSA</th><th>HySure</th><th>ICCV15</th><th>MAPSMM</th><th>SFIMHS</th><th>uSDN</th><th>HyCoNet</th></tr></thead><tbody><tr><td rowspan="3">σ = 0.5</td><td>mSAM</td><td>2.4152</td><td>2.9762</td><td>3.4716</td><td>2.7551</td><td>2.4707</td><td>2.4194</td><td>2.4718</td><td>2.5734</td><td>2.9597</td><td>3.0325</td><td><strong>2.3211</strong></td></tr><tr><td>mPSNR</td><td>32.4564</td><td>32.2554</td><td>26.7902</td><td>30.3093</td><td>33.6348</td><td>32.9597</td><td>31.3817</td><td>31.1324</td><td>28.9370</td><td>32.9082</td><td><strong>34.0320</strong></td></tr><tr><td>ERGAS</td><td>1.4630</td><td>1.5615</td><td>2.6297</td><td>1.7735</td><td>1.3428</td><td>1.3927</td><td>1.6685</td><td>1.6691</td><td>2.0516</td><td>1.5496</td><td><strong>1.3236</strong></td></tr><tr><td rowspan="3">σ = 1</td><td>mSAM</td><td>2.2617</td><td>2.8014</td><td>3.3703</td><td>2.3955</td><td>2.2692</td><td>2.2854</td><td>2.3812</td><td>2.3034</td><td>2.5387</td><td>2.9288</td><td><strong>2.2447</strong></td></tr><tr><td>mPSNR</td><td>33.4572</td><td>32.4425</td><td>27.2244</td><td>32.9145</td><td>34.2168</td><td>33.6857</td><td>32.0225</td><td>33.5899</td><td>31.5298</td><td>31.1590</td><td><strong>34.3232</strong></td></tr><tr><td>ERGAS</td><td>1.3687</td><td>1.5447</td><td>2.5078</td><td>1.3870</td><td>1.2040</td><td>1.3392</td><td>1.5949</td><td>1.3654</td><td>1.5670</td><td>1.6710</td><td><strong>1.1946</strong></td></tr><tr><td rowspan="3">σ = 2</td><td>mSAM</td><td>2.2378</td><td>2.6955</td><td>3.3703</td><td>2.2745</td><td>2.2363</td><td>2.2939</td><td>2.3766</td><td>2.2395</td><td>2.3768</td><td>2.6423</td><td><strong>2.2022</strong></td></tr><tr><td>mPSNR</td><td>33.7191</td><td>32.8198</td><td>27.2985</td><td>34.4669</td><td>34.6838</td><td>33.5168</td><td>32.0668</td><td>34.5392</td><td>33.1865</td><td>31.9779</td><td><strong>34.7950</strong></td></tr><tr><td>ERGAS</td><td>1.3535</td><td>1.4978</td><td>2.4895</td><td>1.2376</td><td>1.1560</td><td>1.3768</td><td>1.5928</td><td>1.2865</td><td>1.3548</td><td>1.6187</td><td><strong>1.1203</strong></td></tr></tbody></table>

<table><caption>TABLE IV
QUANTITATIVE PERFORMANCE COMPARISON WITH THE DIFFERENT ALGORITHMS ON THE WASHINGTON DC DATA. THE BEST ONE IS SHOWN IN BOLD.</caption><thead><tr><th></th><th></th><th>CNMF</th><th>GSOMP</th><th>FUSE</th><th>GLPHS</th><th>GSA</th><th>HySure</th><th>ICCV15</th><th>MAPSMM</th><th>SFIMHS</th><th>uSDN</th><th>HyCoNet</th></tr></thead><tbody><tr><td rowspan="3">σ = 0.5</td><td>mSAM</td><td>8.8273</td><td>13.1330</td><td>8.6505</td><td>6.9509</td><td>7.1772</td><td>9.9828</td><td>8.8196</td><td>7.3361</td><td>7.2102</td><td>7.0720</td><td><strong>3.1984</strong></td></tr><tr><td>mPSNR</td><td>25.8781</td><td>25.3023</td><td>24.9366</td><td>27.5713</td><td>27.3163</td><td>26.4899</td><td>27.5002</td><td>26.9298</td><td>25.2003</td><td>33.0413</td><td><strong>34.8730</strong></td></tr><tr><td>ERGAS</td><td>2.8031</td><td>3.4564</td><td>3.0064</td><td>2.2598</td><td>2.2676</td><td>2.7298</td><td>2.6514</td><td>2.4260</td><td>3.4973</td><td>1.6220</td><td><strong>1.4334</strong></td></tr><tr><td rowspan="3">σ = 1</td><td>mSAM</td><td>7.2022</td><td>10.9944</td><td>8.2811</td><td>5.8102</td><td>5.8450</td><td>7.7087</td><td>8.0205</td><td>5.9814</td><td>6.0686</td><td>7.5995</td><td><strong>3.2828</strong></td></tr><tr><td>mPSNR</td><td>27.8536</td><td>27.9343</td><td>25.7046</td><td>29.3968</td><td>30.0065</td><td>28.7836</td><td>29.9389</td><td>28.8842</td><td>26.9397</td><td>31.1618</td><td><strong>34.9646</strong></td></tr><tr><td>ERGAS</td><td>2.2999</td><td>3.1083</td><td>2.8820</td><td>1.8748</td><td>1.7830</td><td>2.3315</td><td>2.2540</td><td>1.9589</td><td>36.0251</td><td>2.1821</td><td><strong>1.3959</strong></td></tr><tr><td rowspan="3">σ = 2</td><td>mSAM</td><td>5.8034</td><td>9.6210</td><td>7.9871</td><td>4.1380</td><td>4.8727</td><td>6.5962</td><td>7.3153</td><td>3.7289</td><td>4.4360</td><td>4.8462</td><td><strong>3.4581</strong></td></tr><tr><td>mPSNR</td><td>32.6577</td><td>28.9048</td><td>26.2474</td><td>33.4368</td><td>35.6293</td><td>29.8371</td><td>31.1686</td><td>33.7772</td><td>30.8051</td><td>32.8157</td><td><strong>35.8561</strong></td></tr><tr><td>ERGAS</td><td>1.6080</td><td>2.6798</td><td>2.7159</td><td>1.3183</td><td>1.3323</td><td>2.0079</td><td>1.9705</td><td><strong>1.2007</strong></td><td>2.5060</td><td>1.4754</td><td>1.4778</td></tr></tbody></table>

<table><caption>TABLE V
QUANTITATIVE PERFORMANCE COMPARISON WITH THE DIFFERENT ALGORITHMS ON THE UNIVERSITY OF HOUSTON DATA. THE BEST ONE IS SHOWN IN BOLD.</caption><thead><tr><th></th><th></th><th>CNMF</th><th>GSOMP</th><th>FUSE</th><th>GLPHS</th><th>GSA</th><th>HySure</th><th>ICCV15</th><th>MAPSMM</th><th>SFIMHS</th><th>uSDN</th><th>HyCoNet</th></tr></thead><tbody><tr><td rowspan="3">σ = 0.5</td><td>mSAM</td><td>4.1054</td><td>7.8904</td><td>4.5623</td><td>5.3563</td><td>5.0434</td><td>3.2919</td><td>4.7705</td><td>5.5041</td><td>5.3316</td><td>5.7001</td><td><strong>2.6070</strong></td></tr><tr><td>mPSNR</td><td>27.0969</td><td>28.5365</td><td>24.0123</td><td>25.7426</td><td>27.7226</td><td>33.0650</td><td>31.3685</td><td>25.2709</td><td>22.5418</td><td>29.1871</td><td><strong>35.2361</strong></td></tr><tr><td>ERGAS</td><td>1.9525</td><td>3.0026</td><td>2.5708</td><td>2.0911</td><td>1.7382</td><td>1.3990</td><td>1.5260</td><td>2.2232</td><td>4.5349</td><td>2.1537</td><td><strong>1.0251</strong></td></tr><tr><td rowspan="3">σ = 1</td><td>mSAM</td><td>3.2652</td><td>7.9230</td><td>4.3228</td><td>4.7508</td><td>4.3542</td><td>3.2576</td><td>4.7106</td><td>4.3736</td><td>4.8702</td><td>4.9899</td><td><strong>2.6670</strong></td></tr><tr><td>mPSNR</td><td>29.5577</td><td>28.9063</td><td>24.4412</td><td>26.9179</td><td>29.1212</td><td>34.0103</td><td>31.8105</td><td>26.9350</td><td>23.7143</td><td>29.0370</td><td><strong>35.1123</strong></td></tr><tr><td>ERGAS</td><td>1.4203</td><td>2.9949</td><td>2.4456</td><td>1.8281</td><td>1.5293</td><td>1.3078</td><td>1.4552</td><td>1.8481</td><td>3.3914</td><td>2.0462</td><td><strong>1.0252</strong></td></tr><tr><td rowspan="3">σ = 2</td><td>mSAM</td><td>3.2239</td><td>7.6917</td><td>4.1386</td><td>3.6522</td><td>3.5930</td><td>3.0419</td><td>4.1594</td><td>3.0458</td><td>4.0967</td><td>5.0680</td><td><strong>2.6976</strong></td></tr><tr><td>mPSNR</td><td>31.0872</td><td>30.0214</td><td>24.9240</td><td>29.6952</td><td>32.8512</td><td>35.4746</td><td>32.7251</td><td>31.2549</td><td>26.2060</td><td>29.4494</td><td><strong>35.5182</strong></td></tr><tr><td>ERGAS</td><td>1.2696</td><td>2.8725</td><td>2.3121</td><td>1.3591</td><td>1.2708</td><td>1.1233</td><td>1.3302</td><td>1.2350</td><td>2.3154</td><td>1.9161</td><td><strong>1.0379</strong></td></tr></tbody></table>

ples required for supervised learning, only unsupervised fusion methods were used for fair comparison. This is also in line with the application requirements of real scenes. In addition,

due to the fact that the proposed network is a fusion model where the SRF and PSF are unknown and the only assumption is that prior information about the spectral coverage of the MSI
<figure><img src="image_15.png" /><figcaption>Fig. 14. The PSNR of different bands for comparison methods. (a) Pavia University; (b) Indian Pines; (c) Washington DC; (d) University of Houston.</figcaption></figure>

is known, to conduct a fair comparison, the estimated SRFs obtained by the HySure SRF estimation method [29] were used for the compared methods that require a SRF as input. These included CNMF, GSOMP, ICCV15, FUSE, MAPSMM and uSDN. The HySure SRF estimation method also only requires to know the spectral coverage of the HrMSI bands. Based on the GSD ratio of the simulated input image, a Gaussian kernel with a kernel size equal to the GSD ratio was used for all the methods which required a PSF kernel [56]; otherwise, the default source code settings were used.

Firstly, we used color-composite and heatmap images to visually evaluate the performance of the fusion results when $\sigma = 0.5$. In Figs. 10-13, the first column is the color-composite image (RGB image), the second column is the RMSE heatmap of the color-composite image, the third column is the MRAE heatmap, which can be considered to show the pixel-wise error for the reconstructed image cube, and the fourth column is the SAM error, which represents the spectral consistency of each pixel in the reconstructed image.

For most of the methods, the results for the color-composite images in the first column are good. However, for the heatmap in the second column, there are big differences between these methods. The RGB images for GSOMP, HySure, ICCV15 and uSDN produce good results, indicating that these methods fully utilize the input HrMSI and the results retain more information about the HrMSI. Although the GSOMP method produces a good result for the RGB image, the heatmaps for MRAE and SAM are the poorest of those shown in Fig. 10, 12 and 13. A similar result was also reported by Yokoya et al [56], who explained that the reason for this is that, in GSOMP, the high-resolution abundance is only estimated by the HrMSI and the sparsity prior. The errors for MAPSMM and FUSE show a block error distribution, which indicates that estimation of the results using patch-by-patch processing is unstable. For GSA and GLPHS, the errors have an inhomogeneous plaque block distribution. This is because the high-resolution image is obtained by sharpening the low-resolution image by adding spatial detail information and the cumulative error in the sharpened image can induce local irregular errors.

In the MRAE and SAM images, edge errors for objects are unavoidable for all the methods and this situation is particularly obvious in Figs. 10, 12 and 13. This can be explained by the fact that the mixing effect of low-resolution images makes it difficult to achieve better results in heterogeneous regions. Nonetheless, the results of our proposed method achieve the best visual results.

2) *Quantitative Performance:* We investigated the quantitative performance of all the compared methods, and the quality measures obtained are shown in Tables II-V. The mSAM is the mean of the SAM for all pixels and is used to evaluate the spectral consistency of the reconstructed HSI. The mPSNR is the mean PSNR of all the bands and is a measure of the spatial quality. The ERGAS is a global statistical measure used to evaluate the dimensionless global error for fused data.

From Tables II-V, it can be seen that our proposed method produces stable results for different datasets and different PSF deviations. However, the performance is unstable for all of the compared methods, especially for the Washington DC dataset. Due to the complexity of the objects in this imagery, most methods cannot handle the local relationship between the LrHSI and HSI when $\sigma = 0.5$, this happened even to the methods that do not need prior knowledge of the PSF, e.g. GSA, SFIM-HS, GSOMP, HySure and MAP-SMM. For most methods, the performance improves as the standard deviation of the PSF increases. The results for the Indian Pines data are stable because the GSD for these data is the largest and the land objects are the simplest. The results for HySure, ICCV15 and uSDN are stable except for the Washington DC data. Although an adaptive PSF is also implemented in the HySure model, the performance varies depending on the dataset used. Compared with the other methods, CNMF, FUSE, GSA, HySure and ICCV15 have a better spectral consistency.

Fig. 14 shows the PSNR for the different bands of the HSI when $\sigma = 0.5$, showing the reconstructed spatial quality for each band. It is clear that our proposed method significantly outperforms the other methods that were tested. For the Pavia University and University of Houston data, the results for HySure, ICCV15 and GSOMP are good; however, again, all of the compared methods produce poor results for the Washington DC data. The results for the proposed method are not greatly affected by which dataset is used.

V. CONCLUSION

In this paper, we proposed a novel unsupervised deep learning method called HyCoNet to solve the HSI and MSI fusion problem for arbitrary PSFs and SRFs. Three coupled autoencoder nets were designed to extract spectral information from the LrHSI and spatialcontextual information from the HrMSI. Based on these autoencoder nets, the PSF was learned adaptively according to the correlation between the high- and low-resolution abundances, and the SRF was also learned by reconstruction of the autoencoder. Using the joint loss
function, the proposed method can easily be implemented in an end-to-end training manner and provide a straightforward training strategy. The experiments that were performed indicated that the proposed method solved the HSI and MSI fusion problem without knowing the prior information of the PSF and SRF, and produced stable and robust fusion results for arbitrary PSFs and SRFs.

## ACKNOWLEDGMENTS

The authors would like to thank Prof. David Landgrebe for providing the Indian Pines and Washington DC Mall data, Prof. Paolo Gamba for providing the Pavia University data, and the Hyperspectral Image Analysis group and the NSF Funded Center for Airborne Laser Mapping (NCALM) at the University of Houston for providing the CASI University of Houston data. The authors would like to express their appreciation to Dr. Naoto Yokoya for providing MATLAB codes for hyperspectral and multispectral data fusion toolbox and Dr. Ying Qu for providing the uSDN source code.

## REFERENCES

[1] B. Rasti, D. Hong, R. Hang, P. Ghamisi, X. Kang, J. Chanussot, and J. A. Benediktsson, "Feature extraction for hyperspectral imagery: The evolution from shallow to deep (overview and toolbox)," *IEEE Geosci. Remote Sens. Mag.*, 2020, dOI: 10.1109/MGRS.2020.2979764.

[2] L. Gao, J. Li, M. Khodadadzadeh, A. Plaza, B. Zhang, Z. He, and H. Yan, "Subspace-based support vector machines for hyperspectral image classification," *IEEE Geosci. Remote Sens. Lett.*, vol. 12, no. 2, pp. 349–353, 2014.

[3] D. Hong, X. Wu, P. Ghamisi, J. Chanussot, N. Yokoya, and X. X. Zhu, "Invariant attribute profiles: A spatial-frequency joint feature extractor for hyperspectral image classification," *IEEE Trans. Geosci. Remote Sens.*, vol. 58, no. 6, pp. 3791–3808, 2020.

[4] X. Cao, J. Yao, Z. Xu, and D. Meng, "Hyperspectral image classification with convolutional neural network and active learning," *IEEE Trans. Geosci. Remote Sens.*, 2020, dOI:10.1109/TGRS.2020.2964627.

[5] X. Cao, J. Yao, X. Fu, H. Bi, and D. Hong, "An enhanced 3-dimensional discrete wavelet transform for hyperspectral image classification," *IEEE Geosci. and Remote Sens. Lett.*, 2020, 10.1109/LGRS.2020.2990407.

[6] Q. Guo, B. Zhang, Q. Ran, L. Gao, J. Li, and A. Plaza, "Weighted-rxd and linear filter-based rxd: Improving background statistics estimation for anomaly detection in hyperspectral imagery," *IEEE J. Sel. Topics Appl. Earth Observ. Remote Sens.*, vol. 7, no. 6, pp. 2351–2366, 2014.

[7] C. Li, L. Gao, Y. Wu, B. Zhang, J. Plaza, and A. Plaza, "A real-time unsupervised background extraction-based target detection method for hyperspectral imagery," *J. Real-Time Image Process.*, vol. 15, no. 3, pp. 597–615, 2018.

[8] X. Wu, D. Hong, J. Tian, J. Chanussot, W. Li, and R. Tao, "Orsim detector: A novel object detection framework in optical remote sensing imagery using spatial-frequency channel features," *IEEE Trans. Geosci. Remote Sens.*, vol. 57, no. 7, pp. 5146–5158, 2019.

[9] X. Wu, D. Hong, J. Chanussot, Y. Xu, R. Tao, and Y. Wang, "Fourier-based rotation-invariant feature boosting: An efficient framework for geospatial object detection," *IEEE Geosci. Remote Sens. Lett.*, vol. 17, no. 2, pp. 302–306, 2020.

[10] W. He, H. Zhang, L. Zhang, W. Philips, and W. Liao, "Weighted sparse graph based dimensionality reduction for hyperspectral images," *IEEE Geosci. Remote Sens. Lett.*, vol. 13, no. 5, pp. 686–690, 2016.

[11] D. Hong, N. Yokoya, and X. X. Zhu, "Learning a robust local manifold representation for hyperspectral dimensionality reduction," *IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens.*, vol. 10, no. 6, pp. 2960–2975, 2017.

[12] H. Xu, H. Zhang, W. He, and L. Zhang, "Superpixel-based spatial-spectral dimension reduction for hyperspectral imagery classification," *Neurocomputing*, vol. 360, pp. 138–150, 2019.

[13] D. Hong, N. Yokoya, J. Chanussot, J. Xu, and X. X. Zhu, "Learning to propagate labels on graphs: An iterative multitask regression framework for semi-supervised hyperspectral dimensionality reduction," *ISPRS J. Photogramm. Remote Sens.*, vol. 158, pp. 35–49, 2019.

[14] Y. Xu, Z. Wu, J. Chanussot, P. Comon, and Z. Wei, "Nonlocal coupled tensor cp decomposition for hyperspectral and multispectral image fusion," *IEEE Trans. Geosci. Remote Sens.*, vol. 58, no. 1, pp. 348–362, 2019.

[15] J. Hu, D. Hong, and X. X. Zhu, "MIMA: Mapper-induced manifold alignment for semi-supervised fusion of optical image and polarimetric sar data," *IEEE Trans. Geosci. Remote Sens.*, vol. 57, no. 11, pp. 9025–9040, 2019.

[16] R. Hang, Z. Li, P. Ghamisi, D. Hong, G. Xia, and Q. Liu, "Classification of hyperspectral and lidar data using coupled cnns," *IEEE Trans. Geosci. Remote Sens.*, vol. 58, no. 7, pp. 4939–4950, 2020.

[17] M. Tang, L. Gao, A. Marinoni, P. Gamba, and B. Zhang, "Integrating spatial information in the normalized p-linear algorithm for nonlinear hyperspectral unmixing," *IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens.*, vol. 11, no. 4, pp. 1179–1190, 2017.

[18] D. Hong and X. X. Zhu, "SULoRA: Subspace unmixing with low-rank attribute embedding for hyperspectral data analysis," *IEEE J. Sel. Topics Signal Process.*, vol. 12, no. 6, pp. 1351–1363, 2018.

[19] J. Yao, D. Meng, Q. Zhao, W. Cao, and Z. Xu, "Nonconvex-sparsity and nonlocal-smoothness-based blind hyperspectral unmixing," *IEEE Trans. Image Process.*, vol. 28, no. 6, pp. 2991–3006, 2019.

[20] L. Loncan, L. B. De Almeida, J. M. Bioucas-Dias, X. Briottet, J. Chanussot, N. Dobigeon, S. Fabre, W. Liao, G. A. Licciardi, M. Simoes et al., "Hyperspectral pansharpening: A review," *IEEE Geosci. Remote Sens. Mag.*, vol. 3, no. 3, pp. 27–46, 2015.

[21] D. Hong, J. Chanussot, N. Yokoya, J. Kang, and X. X. Zhu, "Learning shared cross-modality representation using multispectral-lidar and hyperspectral data," *IEEE Geosci. Remote Sens. Lett.*, 2020, dOI:10.1109/LGRS.2019.2944599.

[22] R. B. Gomez, A. Jazaeri, and M. Kafatos, "Wavelet-based hyperspectral and multispectral image fusion," in *Geo-Spatial Image and Data Exploitation II*, vol. 4383. International Society for Optics and Photonics, 2001, pp. 36–42.

[23] Y. Zhang and M. He, "Multi-spectral and hyperspectral image fusion using 3-d wavelet transform," *J. Electron.*, vol. 24, no. 2, pp. 218–224, 2007.

[24] Z. Chen, H. Pu, B. Wang, and G.-M. Jiang, "Fusion of hyperspectral and multispectral images: A novel framework based on generalization of pan-sharpening methods," *IEEE Geosci. and Remote Sens. Lett.*, vol. 11, no. 8, pp. 1418–1422, 2014.

[25] B. Aiazzi, S. Baronti, and M. Selva, "Improving component substitution pansharpening through multivariate regression of $ms + pan$ data," *IEEE Trans. Geosci. Remote Sens.*, vol. 45, no. 10, pp. 3230–3239, 2007.

[26] J. Liu, "Smoothing filter-based intensity modulation: A spectral preserve image fusion technique for improving spatial details," *Int. J. Remote Sens.*, vol. 21, no. 18, pp. 3461–3472, 2000.

[27] M. T. Eismann and R. C. Hardie, "Hyperspectral resolution enhancement using high-resolution multispectral imagery with arbitrary response functions," *IEEE Trans. Geosci. Remote Sens.*, vol. 43, no. 3, pp. 455–465, 2005.

[28] Q. Wei, J. Bioucas-Dias, N. Dobigeon, and J.-Y. Tourneret, "Hyperspectral and multispectral image fusion based on a sparse representation," *IEEE Trans. Geosci. Remote Sens.*, vol. 53, no. 7, pp. 3658–3668, 2015.

[29] M. Simoes, J. Bioucas-Dias, L. B. Almeida, and J. Chanussot, "A convex formulation for hyperspectral image superresolution via subspace-based regularization," *IEEE Trans. Geosci. Remote Sens.*, vol. 53, no. 6, pp. 3373–3388, 2014.

[30] N. Akhtar, F. Shafait, and A. Mian, "Bayesian sparse representation for hyperspectral image super resolution," in *Proc. CVPR*. IEEE, 2015, pp. 3631–3640.

[31] R. Kawakami, Y. Matsushita, J. Wright, M. Ben-Ezra, Y.-W. Tai, and K. Ikeuchi, "High-resolution hyperspectral imaging via matrix factorization," in *Proc. CVPR*. IEEE, 2011, pp. 2329–2336.

[32] N. Yokoya, T. Yairi, and A. Iwasaki, "Coupled nonnegative matrix factorization unmixing for hyperspectral and multispectral data fusion," *IEEE Trans. Geosci. Remote Sens.*, vol. 50, no. 2, pp. 528–537, 2011.

[33] C. Lanaras, E. Baltsavias, and K. Schindler, "Hyperspectral super-resolution by coupled spectral unmixing," in *Proc. CVPR*. IEEE, 2015, pp. 3586–3594.

[34] E. Wycoff, T.-H. Chan, K. Jia, W.-K. Ma, and Y. Ma, "A non-negative sparse promoting algorithm for high resolution hyperspectral imaging," in *Proc. ICASSP*. IEEE, 2013, pp. 1409–1413.

[35] N. Akhtar, F. Shafait, and A. Mian, "Sparse spatio-spectral representation for hyperspectral image super-resolution," in *Proc. ECCV*. Springer, 2014, pp. 63–78.
[36] C. Yi, Y.-Q. Zhao, and J. C.-W. Chan, "Hyperspectral image super-resolution based on spatial and spectral correlation fusion," *IEEE Trans. Geosci. Remote Sens.*, vol. 56, no. 7, pp. 4165-4177, 2018.

[37] R. Dian, L. Fang, and S. Li, "Hyperspectral image super-resolution via non-local sparse tensor factorization," in *Proc. CVPR*. Springer, 2017, pp. 5344-5353.

[38] Y. Long, Y. Gong, Z. Xiao, and Q. Liu, "Accurate object localization in remote sensing images based on convolutional neural networks," *IEEE Trans. Geosci. Remote Sens.*, vol. 55, no. 5, pp. 2486-2498, 2017.

[39] X. Wu, D. Hong, P. Ghamisi, W. Li, and R. Tao, "MsRi-CCF: Multi-scale and rotation-insensitive convolutional channel features for geospatial object detection," *Remote Sens.*, vol. 10, no. 12, p. 1990, 2018.

[40] S. Jia, X. Deng, J. Zhu, M. Xu, J. Zhou, and X. Jia, "Collaborative representation-based multiscale superpixel fusion for hyperspectral image classification," *IEEE Trans. Geosci. Remote Sens.*, vol. 57, no. 10, pp. 7770-7784, 2019.

[41] J. Liang, J. Zhou, L. Tong, X. Bai, and B. Wang, "Material based salient object detection from hyperspectral images," *Pattern Recognit.*, vol. 76, pp. 476-490, 2018.

[42] C. Wang, Y. Liu, X. Bai, W. Tang, P. Lei, and J. Zhou, "Deep residual convolutional neural network for hyperspectral image super-resolution," in *Proc. ICIG*. Springer, 2017, pp. 370-380.

[43] X.-H. Han, B. Shi, and Y. Zheng, "Ssf-cnn: Spatial and spectral fusion with cnn for hyperspectral image super-resolution," in *Proc. ICIP*. IEEE, 2018, pp. 2506-2510.

[44] F. Palsson, J. R. Sveinsson, and M. O. Ulfarsson, "Multispectral and hyperspectral image fusion using a 3-d-convolutional neural network," *IEEE Geosci. and Remote Sens. Lett.*, vol. 14, no. 5, pp. 639-643, 2017.

[45] R. Dian, S. Li, A. Guo, and L. Fang, "Deep hyperspectral image sharpening," *IEEE Trans. Neural Netw. Learn. Syst.*, no. 99, pp. 1-11, 2018.

[46] Q. Xie, M. Zhou, Q. Zhao, D. Meng, W. Zuo, and Z. Xu, "Multispectral and hyperspectral image fusion by ms/hs fusion net," in *Proc. CVPR*. IEEE, 2019, pp. 1585-1594.

[47] W. Wang, W. Zeng, Y. Huang, X. Ding, and J. Paisley, "Deep blind hyperspectral image fusion," in *Proc. ICCV*, 2019, pp. 4150-4159.

[48] X.-H. Han, Y. Zheng, and Y.-W. Chen, "Multi-level and multi-scale spatial and spectral fusion cnn for hyperspectral image super-resolution," in *Proc. ICCV workshops*, 2019, pp. 0-0.

[49] Y. Qu, H. Qi, and C. Kwan, "Unsupervised sparse dirichlet-net for hyperspectral image super-resolution," in *Proc. CVPR*. IEEE, 2018, pp. 2511-2520.

[50] Y. Zhou, A. Rangarajan, and P. D. Gader, "An integrated approach to registration and fusion of hyperspectral and multispectral images," *IEEE Trans. Geosci. Remote Sens.*, 2019.

[51] Y. Fu, T. Zhang, Y. Zheng, D. Zhang, and H. Huang, "Hyperspectral image super-resolution with optimized rgb guidance," in *Proc. CVPR*, 2019, pp. 11661-11670.

[52] D. Hong, N. Yokoya, J. Chanussot, and X. X. Zhu, "An augmented linear mixing model to address spectral variability for hyperspectral unmixing," *IEEE Trans. Image Process.*, vol. 28, no. 4, pp. 1923-1938, 2019.

[53] J. Kang, D. Hong, J. Liu, G. Baier, N. Yokoya, and B. Demir, "Learning convolutional sparse coding on complex domain for interferometric phase restoration," *IEEE Trans. Neural Netw. Learn. Syst.*, 2020, doi: 10.1109/TNNLS.2020.2979546.

[54] Q. Wang and P. M. Atkinson, "The effect of the point spread function on sub-pixel mapping," *Remote Sens Environ.*, vol. 193, pp. 127-137, 2017.

[55] D. Hong, W. Liu, J. Su, Z. Pan, and G. Wang, "A novel hierarchical approach for multispectral palmprint recognition," *Neurocomputing*, vol. 151, pp. 511-521, 2015.

[56] N. Yokoya, C. Grohnfeldt, and J. Chanussot, "Hyperspectral and multispectral data fusion: A comparative review of the recent literature," *IEEE Geosci. Remote Sens. Mag.*, vol. 5, no. 2, pp. 29-56, 2017.

[57] B. Le Saux, N. Yokoya, R. Hansch, and S. Prasad, "2018 ieee grss data fusion contest: Multimodal land use classification [technical committees]," *IEEE Geosci. Remote Sens. Mag.*, vol. 6, no. 1, pp. 52-54, 2018.

[58] D. Hong, N. Yokoya, N. Ge, J. Chanussot, and X. X. Zhu, "Learnable manifold alignment (LeMA): A semi-supervised cross-modality learning framework for land cover and land use classification," *ISPRS J. Photogramm. Remote Sens.*, vol. 147, pp. 193-205, 2019.

[59] J. A. Barsi, K. Lee, G. Kvaran, B. L. Markham, and J. A. Pedelty, "The spectral response of the landsat-8 operational land imager," *Remote Sens.*, vol. 6, no. 10, pp. 10232-10251, 2014.

[60] B. Arad, O. Ben-Shahar, and R. Timofte, "Ntire 2018 challenge on spectral reconstruction from rgb images," in *Proc. CVPR*. IEEE, 2018, pp. 929-938.

[61] D. P. Kingma and J. Ba, "Adam: A method for stochastic optimization," in *Proc. ICLR*, 2014.

[62] A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, Z. Lin, N. Gimelshein, L. Antiga et al., "Pytorch: An imperative style, high-performance deep learning library," in *Proc. NIPS*, 2019, pp. 8024-8035.

[63] D. Hong, J. Chanussot, N. Yokoya, U. Heiden, W. Heldens, and X. X. Zhu, "Wu-net: A weakly-supervised unmixing network for remotely sensed hyperspectral imagery," in *Proc. IGARSS*. IEEE, 2019, pp. 373-376.

[64] Q. Wei, N. Dobigeon, and J.-Y. Tourneret, "Fast fusion of multi-band images based on solving a sylvester equation," *IEEE Trans. Image Process.*, vol. 24, no. 11, pp. 4109-4121, 2015.

[65] M. Selva, B. Aiazzi, F. Butera, L. Chiarantini, and S. Baronti, "Hyper-sharpening: A first approach on sim-ga data," *IEEE J. Sel. Topics Appl. Earth Observ. Remote Sens.*, vol. 8, no. 6, pp. 3008-3024, 2015.

<figure><img src="image_16.png" /></figure>

**Ke Zheng** received B.S. degree in geographic information system from Shandong Agricultural University, Taian, China, in 2008, and the M.S. and Ph.D. degrees in remote sensing from College of Geosciences and Surveying Engineering, China University of Mining and Technology (Beijing), Beijing, China, in 2016 and 2020, respectively.

Since 2020, he has been a Post-Doctoral Associate with the Key Laboratory of Digital Earth Science, Aerospace Information Research Institute, Chinese Academy of Science, Beijing, China. His research interests include image processing, machine learning, deep learning and their application in Earth Vision.

<figure><img src="image_17.png" /></figure>

**Lianru Gao** (M'12-SM'18) received the B.S. degree in civil engineering from Tsinghua University, Beijing, China, in 2002, and the Ph.D. degree in cartography and geographic information system from Institute of Remote Sensing Applications, Chinese Academy of Sciences (CAS), Beijing, China, in 2007.

He is currently a Professor with the Key Laboratory of Digital Earth Science, Aerospace Information Research Institute, CAS. He also has been a visiting scholar at the University of Extremadura, Cceres, Spain, in 2014, and at the Mississippi State University (MSU), Starkville, USA, in 2016. His research focuses on hyperspectral image processing and information extraction. In last ten years, he was the PI of 10 scientific research projects at national and ministerial levels, including projects by the National Natural Science Foundation of China (2010-2012, 2016-2019, 2018-2020), and by the Key Research Program of the CAS (2013-2015) et al. He has published more than 160 peer-reviewed papers, and there are more than 80 journal papers included by Science Citation Index (SCI). He was coauthor of an academic book entitled "Hyperspectral Image Classification And Target Detection". He obtained 28 National Invention Patents in China. He was awarded the Outstanding Science and Technology Achievement Prize of the CAS in 2016, and was supported by the China National Science Fund for Excellent Young Scholars in 2017, and won the Second Prize of The State Scientific and Technological Progress Award in 2018. He received the recognition of the Best Reviewers of the IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing in 2015, and the Best Reviewers of the IEEE Transactions on Geoscience and Remote Sensing in 2017.
<figure><img src="image_18.png" /></figure>

**Wenzhi Liao** (Senior Member, IEEE) received the B.Sc. degree in mathematics from Hainan Normal University, Haikou, China, in 2006, the Ph.D. degree in engineering from the South China University of Technology, Guangzhou, China, in 2012, and the Ph.D. degree in computer science engineering from Ghent University, Ghent, Belgium, in 2012.

From 2012 to 2019, he has worked as a Post-Doctoral Research Fellow first with Ghent University and then with the Research Foundation Flanders (FWO), Vlaanderen, Belgium. Since 2020, he has

been with the Sustainable Materials Management, Flemish Institute for Technological Research (VITO), Mol, Belgium. He also works as the Guest Professor of Ghent University. His research interests include image processing, pattern recognition, and remote sensing. In particular, his interests include mathematical morphology, multitask feature learning, multisensor data fusion, and hyperspectral image (HSI) restoration.

Dr. Liao was a recipient of the Best Paper Challenge Awards in both the 2013 IEEE GRSS Data Fusion Contest and the 2014 IEEE GRSS Data Fusion Contest. He serves as an Associate Editor for the IEEE JOURNAL OF SELECTED TOPICS IN APPLIED EARTH OBSERVATIONS AND REMOTE SENSING (JSTARS).

<figure><img src="image_19.png" /></figure>

**Danfeng Hong** (S'16-M'19) received the M.Sc. degree (summa cum laude) in computer vision, College of Information Engineering, Qingdao University, Qingdao, China, in 2015, the Dr. -Ing degree (summa cum laude) in Signal Processing in Earth Observation (SiPEO), Technical University of Munich (TUM), Munich, Germany, in 2019.

Since 2015, he also worked as a Research Associate at the Remote Sensing Technology Institute (IMF), German Aerospace Center (DLR), Oberpfaffenhofen, Germany. Currently, he is research scientist

and leads a Spectral Vision group at IMF, DLR, and also an adjunct scientist in GIPSA-lab, Grenoble INP, CNRS, Univ. Grenoble Alpes, Grenoble, France.

His research interests include signal / image processing and analysis, hyperspectral remote sensing, machine / deep learning, artificial intelligence and their applications in Earth Vision.

<figure><img src="image_20.png" /></figure>

**Bing Zhang** (M'11SM'12-F'19) received the B.S. degree in geography from Peking University, Beijing, China, in 1991, and the M.S. and Ph.D. degrees in remote sensing from the Institute of Remote Sensing Applications, Chinese Academy of Sciences (CAS), Beijing, China, in 1994 and 2003, respectively.

Currently, he is a Full Professor and the Deputy Director of the Aerospace Information Research Institute, CAS, where he has been leading lots of key scientific projects in the area of hyperspectral remote

sensing for more than 25 years. His research interests include the development of Mathematical and Physical models and image processing software for the analysis of hyperspectral remote sensing data in many different areas. He has developed 5 software systems in the image processing and applications. His creative achievements were rewarded 10 important prizes from Chinese government, and special government allowances of the Chinese State Council. He was awarded the National Science Foundation for Distinguished Young Scholars of China in 2013, and was awarded the 2016 Outstanding Science and Technology Achievement Prize of the Chinese Academy of Sciences, the highest level of Awards for the CAS scholars.

Dr. Zhang has authored more than 300 publications, including more than 170 journal papers. He has edited 6 books/contributed book chapters on hyperspectral image processing and subsequent applications. He is the IEEE fellow and currently serving as the Associate Editor for IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing. He has been serving as Technical Committee Member of IEEE Workshop on Hyperspectral Image and Signal Processing since 2011, and as the president of hyperspectral remote sensing committee of China National Committee of International Society for Digital Earth since 2012, and as the Standing Director of Chinese Society of Space Research (CSSR) since 2016. He is the Student Paper Competition Committee member in IGARSS from 2015-2019.

<figure><img src="image_21.png" /></figure>

**Ximin Cui** received the B.S. degree in mine surveying from China University of Mining and Technology, Xuzhou, China, in 1990, the M.S. degree in engineering surveying from China University of Mining and Technology, Beijing, China, in 1993, and the Ph.D. degree in mine engineering mechanics from China University of Mining and Technology, Beijing, China, in 1996.

He is currently a Professor with the College of Geosciences and Surveying Engineering, China University of Mining and Technology (Beijing). He

was the PI of 9 scientific research projects at national and ministerial levels, including projects by the National Natural Science Foundation of China (2000-2002, 2011-2013, 2015-2018) et al. He has published more than 170 peer-reviewed papers, and there are more than 70 papers included by Science Citation Index (SCI) and Engineering Village (EI).

His research interests include mining subsidence, deformation monitoring and industrial surveying.

<figure><img src="image_22.png" /></figure>

**Jocelyn Chanussot** (M'04SM'04-F'12) received the M.Sc. degree in electrical engineering from the Grenoble Institute of Technology (Grenoble INP), Grenoble, France, in 1995, and the Ph.D. degree from the Université de Savoie, Annecy, France, in 1998. Since 1999, he has been with Grenoble INP, where he is currently a Professor of signal and image processing. His research interests include image analysis, hyperspectral remote sensing, data fusion, machine learning and artificial intelligence. He has been a visiting scholar at Stanford University (USA),

KTH (Sweden) and NUS (Singapore). Since 2013, he is an Adjunct Professor of the University of Iceland. In 2015-2017, he was a visiting professor at the University of California, Los Angeles (UCLA). He holds the AXA chair in remote sensing and is an Adjunct professor at the Chinese Academy of Sciences, Aerospace Information research Institute, Beijing.

Dr. Chanussot is the founding President of IEEE Geoscience and Remote Sensing French chapter (2007-2010) which received the 2010 IEEE GRS-S Chapter Excellence Award. He has received multiple outstanding paper awards. He was the Vice-President of the IEEE Geoscience and Remote Sensing Society, in charge of meetings and symposia (2017-2019). He was the General Chair of the first IEEE GRSS Workshop on Hyperspectral Image and Signal Processing, Evolution in Remote sensing (WHISPERS). He was the Chair (2009-2011) and Cochair of the GRS Data Fusion Technical Committee (2005-2008). He was a member of the Machine Learning for Signal Processing Technical Committee of the IEEE Signal Processing Society (2006-2008) and the Program Chair of the IEEE International Workshop on Machine Learning for Signal Processing (2009). He is an Associate Editor for the IEEE Transactions on Geoscience and Remote Sensing, the IEEE Transactions on Image Processing and the Proceedings of the IEEE. He was the Editor-in-Chief of the IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing (2011-2015). In 2014 he served as a Guest Editor for the IEEE Signal Processing Magazine. He is a Fellow of the IEEE, a member of the Institut Universitaire de France (2012-2017) and a Highly Cited Researcher (Clarivate Analytics/Thomson Reuters, 2018-2019).