# Unsupervised and Unregistered Hyperspectral Image Super-Resolution With Mutual Dirichlet-Net

Ying Qu, Hairong Qi, Chiman Kwan, Naoto Yokoya, Jocelyn Chanussot

► To cite this version:

Ying Qu, Hairong Qi, Chiman Kwan, Naoto Yokoya, Jocelyn Chanussot. Unsupervised and Un-registered Hyperspectral Image Super-Resolution With Mutual Dirichlet-Net. IEEE Transactions on Geoscience and Remote Sensing, 2022, 60, pp.5507018. 10.1109/TGRS.2021.3079518. hal-03932319

HAL Id: hal-03932319

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

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
# Unsupervised and Unregistered Hyperspectral Image Super-Resolution with Mutual Dirichlet-Net

Ying Qu, Member, IEEE, Hairong Qi, Fellow, IEEE, Chiman Kwan Senior Member, IEEE, Naoto Yokoya, Member, IEEE, and Jocelyn Chanussot, Fellow, IEEE

**Abstract**—(Please find the final version from IEEE Transactions on Geoscience and Remote Sensing on IEEE Xplore. The code has been released on GitHub at https://github.com/yingutk/u2MDN.) Hyperspectral images (HSI) provide rich spectral information that has contributed to the successful performance improvement of numerous computer vision and remote sensing tasks. However, it can only be achieved at the expense of images' spatial resolution. Hyperspectral image super-resolution (HSI-SR) thus addresses this problem by fusing low resolution (LR) HSI with multispectral image (MSI) carrying much higher spatial resolution (HR). Existing HSI-SR approaches require the LR HSI and HR MSI to be well registered and the reconstruction accuracy of the HR HSI relies heavily on the registration accuracy of different modalities. In this paper, we propose an unregistered and unsupervised mutual Dirichlet-Net ($u^2$-MDN) to exploit the uncharted problem domain of HSI-SR *without the requirement of multi-modality registration*. The success of this endeavor would largely facilitate the deployment of HSI-SR since registration requirement is difficult to satisfy in real-world sensing devices. The novelty of this work is three-fold. First, to stabilize the fusion procedure of two unregistered modalities, the network is designed to extract spatial and spectral information of two modalities with different dimensions through a shared encoder-decoder structure. Second, the mutual information (MI) is further adopted to capture the non-linear statistical dependencies between the representations from two modalities (carrying spatial information) and their raw inputs. By maximizing the MI, spatial correlations between different modalities can be well characterized to further reduce the spectral distortion. We assume the representations follow a similar Dirichlet distribution for its inherent sum-to-one and non-negative properties. Third, a collaborative $l_{2,1}$ norm is employed as the reconstruction error instead of the more common $l_2$ norm to better preserve the spectral information. Extensive experimental results demonstrate the superior performance of $u^2$-MDN as compared to the state-of-the-art.

**Index Terms**—Hyperspectral image, unregistered, super-resolution, mutual information, unsupervised deep learning

## I. INTRODUCTION

Hyperspectral image (HSI) collects hundreds of contiguous spectral representations of objects, which demonstrates advantages over the conventional multispectral image (MSI) or RGB image with much less spectral information [1], [2]. Compared

Ying Qu, and Hairong Qi are with the Advanced Imaging and Collaborative Information Processing Group, Department of Electrical Engineering and Computer Science, University of Tennessee, Knoxville, TN 37996 USA (e-mail: yqu3@vols.utk.edu; hqi@utk.edu).

Chiman Kwan is with Applied Research LLC, Rockville, MD, 20850 USA (e-mail: chiman.kwan@arllc.net)

Naoto Yokoya is with RIKEN Center for Advanced Intelligence Project (AIP) Tokyo, 103-0027, Japan. (e-mail: naoto.yokoya@riken.jp)

Jocelyn Chanussot is with the Univ. Grenoble Alpes, Inria, CNRS, Grenoble INP, LJK, Grenoble, 38000, France. (e-mail: jocelyn@hi.is).

<figure><img src="image_2.png" /><figcaption>Fig. 1. Unregistered hyperspectral image super-resolution. (a) First band of the 20 degree rotated and cropped LR HSI with 38% information missing. (b) First band of the HR MSI. (c) First band of the reconstructed HR HSI by the proposed methods. (d) First band of the reference HR HSI.</figcaption></figure>

to conventional images, the rich spectral information of HSI can effectively distinguish visually similar objects that actually consist of different materials. Thus, HSI has been shown to enhance the performance of a wide range of computer vision and remote sensing tasks, such as, object recognition and classification [3]–[5], segmentation [6], tracking [7], environmental monitoring [8], and change detection [9].

During the HSI acquisition process, the finer the spectral resolution, the smaller the radiation energy that can reach the sensor for a particular spectral band within a narrow wavelength range. Thus, the high spectral resolution of HSI can only be achieved at the cost of its spatial resolution due to the hardware limitations [10], [11]. On the contrary, we can obtain conventional MSI or RGB with a much higher spatial resolution by integrating the radiation energy over broad spectral bands which inevitably reduces their spectral resolution significantly [12]. To improve the spatial resolution of HSI for better application performance, a natural way is to fuse the high spectral information extracted from HSI with the high-resolution spatial information extracted from conventional images to yield high resolution images in both spatial and spectral domains [4], [13]. This procedure is referred to as hyperspectral image super-resolution (HSI-SR) [11], [12].

HSI-SR can be broadly divided into three categories, traditional component substitution (CS) [14], [15] and multi-resolution analysis (MRA) based methods [16], matrix factorization based, and Bayesian-based approaches [4], [17]. Although HSI-SR has been intensively studied, spectral distortion can be easily introduced during the optimization procedure of methods from these categories. Recently, there have been several attempts to address the HSI-SR problem with deep learning where the mapping function between the LR HSI and HR HSI is learned using different frameworks [18], [19]. However, the deep learning-based approaches are generally limited to handle image pairs with large spatial-scale differ-
ences and the learned mapping function may not be readily
adapted to reconstruct HR HSI possessing different spectral
characteristics or acquired from different sensors.

Despite a plethora of works on HSI-SR, all current ap-
proaches have at least one pre-requisite to solving the problem
of HSI-SR, i.e., the two input modalities (HSI and MSI) must
be well registered, and the quality of the reconstructed HR HSI
relies heavily on the registration accuracy [2], [4], [20]–[22].
According to previous works, there are a few methods that
introduce registration as a pre-step before data fusion [17],
[23], [24]. However, these pre-steps can only handle small-
scale differences, e.g., two pixels/eight pixels offset in LR
HSI/HR MSI [20]. Moreover, even in the registration commu-
nity, HSI and MSI registration is a challenging problem itself
as one pixel in LR HSI may cover hundreds of pixels in the
corresponding HR MSI. The spectral difference is also large
that both the spectral response function (SRF) and multi-band
images have to be taken into consideration during registration
[20], [25]–[28].

In this paper, an unsupervised network structure is proposed,
aiming to solve the HSI-SR problem directly without multi-
modality registration. An example is shown in Fig. 1. We
address the problem based on the assumption that, the pixels
in the overlapped region of HR HSI and HR MSI can be
approximated by a linear combination of the same spectral in-
formation (spectral bases) with the same corresponding spatial
information (representations), which indicates how the spectral
basis is constructed for each pixel. Since LR HSI is the down-
sampled version of the HR HSI, ideally, its representations
should be correlated with that of the HR MSI and HR HSI, i.e.,
they should follow similar patterns and distributions although
possessing different resolutions, as shown in Fig. 2. Therefore,
to reconstruct HR HSI with minimum spectral distortion, the
network is designed to decouple both the LR HSI and HR MSI
into spectral bases and representations, such that their spectral
bases are shared and their representations are correlated with
each other.

<figure><img src="image_3.png" /><figcaption>Fig. 2. Learned hidden representations from unregistered (a) low resolution HSI and (b) high resolution MSI, respectively, as shown in Fig. 1.</figcaption></figure>

The novelty of this work is three-fold. First, to stabilize the
fusion procedure for two unregistered modalities, the network
extracts both the spectral and spatial information of the multi-
modalities through the same encoder-decoder structure, by
projecting the LR HSI onto the same statistical space as
HR MSI, as illustrated in Fig. 3. The representations of the
network are encouraged to follow a Dirichlet distribution
to naturally meet the non-negative and sum-to-one physical
constraints. Second, to prevent spectral distortion, we further

adopt mutual information (MI) to extract optimal and cor-
related representations from multi-modalities. Since the two-
modalities are unregistered, the correlated representations are
learned by maximizing the MI between the representations
and their own inputs during the network optimization. Third, a
collaborative l₂,₁ norm is employed as the reconstruction error
instead of commonly used l₂ loss, so that the network is able
to reconstruct individual pixels as accurately as possible. In
this way, the network preserves the spectral information better.
With the above design, the proposed network is able to work
directly on unregistered images and the spectral distortion of
the reconstructed HR HSI can be largely reduced. The pro-
posed method is referred to as unregistered and unsupervised
mutual Dirichlet Net, or u²-MDN for short.

u²-MDN is an extension of our previous work uSDN [21]. However, uSDN is only effective on general HSI-SR problem with well-registered LR HSI and HR MSI. Here, we have made substantial extensions to address the challenges of HSI-SR with unregistered multi-modalities. To the best of our knowledge, this is the first effort to solving the HSI-SR problem directly on unregistered image pairs with unsupervised deep learning. The major improvements can be summarized from three perspectives. First, the network structure is different from that of the uSDN. Instead of adopting two deep learning networks as in uSDN, the proposed u²-MDN is specifically designed to extract the representations of multi-modalities with only one encoder-decoder structure, which largely stabilizes the information extraction and fusion procedure given the unregistered multi-modalities. Second, uSDN minimizes spectral distortion of the reconstructed HR HSI by reducing the angular difference of the representations from multiple modalities, which fails to deal with unregistered cases, while the proposed u²-MDN is able to handle both well-registered and unregistered cases by extracting correlated representations with mutual information through the mutual discriminative network. Third, instead of commonly used l₂ loss adopted by uSDN, the collaborative l₂,₁ norm is introduced in the proposed u²-MDN to better preserve the spectral information.

## II. RELATED WORK

### A. Hyperspectral Image Super-Resolution

The problem of HSI-SR originates from multispectral im-
age super-resolution (MSI-SR) in the remote sensing field,
where the spatial resolution of MSI is further improved by a
high-resolution panchromatic image (PAN). Traditional widely
utilized MSI-SR methods can be roughly categorized into
two groups: the component substitution (CS) and the multi-
resolution analysis (MRA) based approaches. Generally, CS-
based approaches [14] project the given data onto a predefined
space where the spectral information and spatial information
are separated. Subsequently, the spatial component is substi-
tuted with the one extracted from PAN [15]. Several methods
based on CS have been proposed to address the problem of
hyper-sharpening and achieved promising results with different
criteria [29]–[31]. MRA-based approaches achieve the spatial
details by first applying a spatial filter to the HR images. Then
the spatial details are injected into the LR HSI [16], [17], [32]–
[34]. Although these traditional pan-sharpening approaches
can be extended to solve the HSI-SR problem, they usually
suffer from severe spectral distortions [11], [17], [35].

Recent approaches consist of Bayesian-based and matrix factorization-based methods [4], [17]. Bayesian approaches estimate the posterior distribution of the HR HSI given LR HSI and HR MSI. The unique framework of Bayesian offers a convenient way to regularize the solution space of HR HSI by employing a proper prior distribution such as Gaussian. Different methods vary according to the different prior distributions adopted. Wei et al. proposed a Bayesian Naive method [36] based on the assumption that the representation coefficients of HR HSI follow a Gaussian distribution. However, this assumption does not always hold especially when the ground truth HR HSI contains complex textures. Instead of using Gaussian prior, dictionary-based approaches solve the problem under the assumption that HR HSI is a linear combination of properly chosen over-complete dictionary and sparse coefficients [37]. Simoes et al. proposed HySure [38], which takes into account both the spatial and spectral characteristics of the given data. This approach solves the problem through vector-based total variation regularization. Akhtar et al. [11] introduced a non-parametric Bayesian strategy to solve the HSI-SR problem. The method first learns a spectral dictionary from LR HSI under the Bayesian framework. Then it estimates the spatial coefficients of the HR MSI by Bayesian sparse coding. Eventually, the HR HSI is generated by combining the spatial dictionary with the spatial coefficients. However, the spectral information extracted from LR HSI may not be the optimal spectral bases for MSI, since MSI is not utilized during the optimization procedure.

Matrix factorization-based approaches have been actively studied recently [10], [12], [39], [40], with Kawakami et al. [10] being the first that introduced matrix factorization to solve the HSI-SR problem. The method learns a spectral basis from LR HSI and then uses this basis to extract sparse coefficients from HR MSI with non-negative constraints. Similar to Bayesian-based approaches, the HR HSI is generated by linearly combining the estimated bases with the coefficients. Yokoya et al. [39] decomposed both the LR HSI and HR MSI alternatively to achieve the optimal non-negative bases and coefficients that are used to generate HR HSI. Wycoff et al. [41] solved the problem with alternating direction method of multipliers (ADMM). Lanaras et al. [12] further improved the fusion results by introducing a sparse constraint. However, most methods [12], [39], [41] are based on the same assumption that the down-sampling function between the spatial coefficients of HR HSI and LR HSI is known beforehand. In practice, this assumption is not always true due to the complex environmental conditions.

Most of these approaches focus on the spectral
characteristics of the HSI, where the spectral information of
the HSI is extracted while the spatial relationship between
pixels is untouched. Recently, there have been a few
approaches proposed to address the HSI-SR problem based
on tensor decomposition [42]–[45], which explored both the
spectral and spatial correlations of the HSI by learning a core
tensor and the dictionaries along three dimensions, i.e., the
spectral dimension, and two spatial dimensions. In this way,

the information of each dimension can be represented with
its own dictionary, while the core-tensor is shared among
multi-modalities. Although this formulation works well on
well-registered images, it is problematic on unregistered
image pairs, since the core-tensor cannot be shared between
HSI and MSI with large displacements. In addition, it might
limit the reconstruction ability of the method on remote
sensing images which have only a few redundant structures
on the spatial domain of HSI.

Chen et al. [46] proposed to simultaneously register im-
ages during the fusing process. However, it only works on
panchromatic and MSI. Zhou et al. [47] proposed an integrated
approach for registration and fusion, which addressed the
problem of HSI-SR on unregistered image pairs. However,
the registration is still a required step, and the fusion and reg-
istration are performed independently, which would introduce
additional errors during optimization.

B. Deep learning based Super-Resolution

Deep learning attracts increasing attention for natural image
super-resolution since 2014 when Dong et al. first introduced
convolution neural network (CNN) to solve the problem of
natural image super-resolution and demonstrated state-of-the-
art restoration quality [48]. Ledig et al. proposed a method
based on generative adversarial network and skipped residual
network. The method employed perceptual loss through the
VGG network which can recover photo-realistic textures from
heavily down-sampled images [49]. Usually, natural image
SR methods only work up to 8 times upscaling. There have
been several attempts to address the MSI-SR or HSI-SR with
deep learning in a supervised fashion. In 2015, a modified
sparse tied-weights denoising autoencoder was proposed by
Huang et al. [50] to enhance the resolution of MSI. The
method assumes that the mapping function between LR and
HR PAN is the same as the one between LR and HR MSI.
Masi et al. proposed a supervised three-layer SRCNN [51] to
learn the mapping function between LR MSI and HR MSI.
Similar to [51], Wei et al. [52] learned the mapping function
with deep residual network. Li et al. [53] solved the HSI-SR
problem by learning a mapping function with spatial constraint
strategy and convolutional neural network (CNN). Dian et
al. [54] initialized the HR HSI from the fusion framework
via the Sylvester equation. Then, the mapping function is
trained between the initialized HR-HSI and the reference HR
HSI through deep residual learning. Xie et al. [55] reduced
spectral distortions of the reconstructed HR HSI by exploiting
the approximate low-rankness prior along the spectral domain
of the HSI. However, these supervised deep learning-based
methods can not be readily adopted on HSI-SR for real
applications due to three reasons. First, the scale differences
between LR HSI and HR MSI can reach as large as 10, i.e., one
pixel in HSI covers 100 pixels in MSI. In some applications,
the scale difference can even be 25 [56] and 30 [57]. But
most existing super-resolution methods only work on up to 8
times upscaling. Second, they are designed to find an end-to-
end mapping function between the LR images and HR images
under the assumption that the mapping function is the same for different images. However, the mapping function may not remain the same for images acquired with different sensors. Even for the data collected from the same sensor, the mapping function for different spectral bands may not be the same. Thus the assumption may cause severe spectral distortion. Third, training a mapping function is a supervised problem that requires a large dataset, the down-sampling function, and the availability of the HR HSI, making supervised learning unrealistic for HSI.

Recently, we proposed an unsupervised uSDN [21], which addressed the problem of HSI-SR with deep network models. Specifically, it extracts the spectral and spatial information through two encoder-decoder networks from the two modalities. The angular difference between the LR HSI and HR MSI representations is minimized to reduce the spectral distortion for every ten iterations. Fu et al. [58] proposed an unsupervised CNN-based method for HSI super-resolution, which learns a mapping function between the RGB space and the spectral space with spatial constraint for the HR HSI. Zheng et al. [59] proposed an unsupervised method with learnable downsampling function based on the theory of linear unmixing. These methods can achieve promising results for different HSI datasets. However, they are specifically designed for well-registered image pairs.

### III. PROBLEM FORMULATION

Given the LR HSI, $\bar{Y}_h \in \mathbb{R}^{m \times n \times L}$, where $m, n$ and $L$ denote its width, height and number of spectral bands, respectively, and the unregistered HR MSI with overlapped region, $\bar{Y}_m \in \mathbb{R}^{M \times N \times l}$, where $M, N$ and $l$ denote its width, height and number of spectral bands, respectively, the goal is to reconstruct the HR HSI $\bar{X} \in \mathbb{R}^{M \times N \times L}$ based on the content of HR MSI. In general, MSI has much higher spatial resolution than HSI, and HSI has much higher spectral resolution than MSI, i.e., $M \gg m, N \gg n$ and $L \gg l$.

To facilitate the subsequent processing, we unfold the 3D images into 2D matrices, $Y_h \in \mathbb{R}^{mn \times L}, Y_m \in \mathbb{R}^{MN \times l}$ and $X \in \mathbb{R}^{MN \times L}$, such that each row represents the spectral reflectance of a single pixel. Since each pixel in both LR HSI and HR MSI can be approximated by a linear combination of c spectral bases $D$ [11], [12], [21], the matrices can be further decomposed as

$$Y_h = S_h D_h \quad (1)$$

$$Y_m = S_m D_m \quad (2)$$

$$X = S_m D_h \quad (3)$$

where $D_h \in \mathbb{R}^{c \times L}, D_m \in \mathbb{R}^{c \times l}$ denote the spectral bases of LR HSI and HR MSI, respectively. $S_h \in \mathbb{R}^{mn \times c}, S_m \in \mathbb{R}^{MN \times c}$ denote the coefficients of LR HSI and HR MSI, respectively. Since $S_h$ or $S_m$ indicate how the spectral bases are combined for individual pixels at specific locations, they preserve the spatial structure of HSI. Note that the benefit of unfolding the data into 2D matrices is that, the extraction procedure can decouple each pixel without changing the relationship of the pixel and its neighborhood pixels, thus the reconstructed image has less artifacts [11], [12], [21].

In real applications, although the areas captured by LR HSI and HR MSI might not be registered well, they always have overlapping regions, and the LR HSI includes all the spectral basis of HR MSI i.e., they share the same type of materials carrying specific spectral signatures. The relationship between LR HSI and HR MSI can be expressed as

$$C_h \neq C_m, C_h \cap C_m \neq \emptyset, D_m = D_h \mathcal{R}, \quad (4)$$

where $C_h$ and $C_m$ denote the contents of LR HSI and HR MSI, respectively. $\mathcal{R} \in \mathbb{R}^{L \times l}$ is the prior transformation matrix of sensor [10], [12], [13], [17], [21], [35], [37]–[39], which describes the relationship between HSI and MSI bases.

With $D_h \in \mathbb{R}^{c \times L}$ carrying the high-resolution spectral information and $S_m \in \mathbb{R}^{MN \times c}$ carrying the high-resolution spatial information, the desired HR HSI, $X$, is generated by Eq. (3).

The challenges to solve this problem are that 1) the ground truth $X$ is not available, and 2) the LR HSI and HR MSI do not cover the same region. To solve this unsupervised and unregistered HR-HSI problem, the key is to take advantage of the shared spectral information $D_h$ among different modalities. In addition, the representations of both modalities specifying the spatial information of scene should meet the non-negative and sum-to-one physical constraints. Moreover, in the ideal case, for the pixels in the overlapped region between LR HSI and HR MSI, their spatial information should follow similar patterns, because they carry the information of how the reflectance of shared materials (spectral basis) are mixed in each location. Therefore, the network should have the ability to learn correlated spatial and spectral information from unregistered multi-modality images to maximize its ability to prevent spectral distortion.

### IV. PROPOSED APPROACH

We propose an unsupervised architecture for unregistered LR HSI and HR MSI as shown in Fig. 3. Here, we highlight the structural uniqueness of the network. To extract correlated spectral and spatial information of unregistered multi-modalities, the network projects the LR HSI into the same statistical space as HR MSI, so that the two modalities can share the same encoder and decoder. The encoder enforces the representations (carrying spatial information) of both modalities to follow a Dirichlet distribution, to naturally meet the non-negative and sum-to-one physical properties. In order to prevent spectral distortion, mutual information is introduced during optimization to maximize the correlation between the representations of LR HSI and HR MSI. And the collaborative $l_{2,1}$ loss is adopted to encourage the network to extract accurate spectral and spatial information from both modalities.

#### A. Network Architecture

As shown in Fig. 3, the network reconstructs both the LR HSI $Y_h$ and HR MSI $Y_m$ by sharing the same encoder and decoder network structure. Since the number of the spectral band $L$ of the HSI $Y_h$ is much larger than that of the spectral band $l$ of MSI $Y_m$, we project $Y_h$ into an $l$ dimensional space by $\bar{Y}_h = Y_h \mathcal{R}$, such that $\bar{Y}_h$ represents the LR MSI lying in
<figure><img src="image_4.png" /><figcaption>Fig. 3. Simplified architecture of $u^2$-MDN.</figcaption></figure>

the same space as HR MSI. In this way, both modalities are linked to share the same encoder structure without additional parameters.

On the other hand, the spectral information $D_m$ of MSI is highly compressed from that of HSI, i.e., $D_m = D_h \mathcal{R}$. Thus, it is very unstable and difficult to directly extract $D_h$, carrying high spectral resolution from, MSI with low-spectral resolution. But the spectral basis of HR MSI can be transformed from those of LR HSI which possesses more spectral information, i.e., $\hat{\mathbf{Y}}_m = \mathbf{S}_m D_m = \mathbf{S}_m D_h \mathcal{R} = \mathbf{X} \mathcal{R}$. Therefore, in the network design, both modalities share the same decoder structure $D_h$. The transformation matrix $\mathcal{R}$ is added as fixed weights to reconstruct the HR MSI $\hat{\mathbf{Y}}_m$. Then the output of the layer before the fixed weights is actually $\mathbf{X}$, according to Eq. (3).

Let us define the input domain as $\mathcal{Y} = \{\hat{\mathbf{Y}}_h, \mathbf{Y}_m\}$, output domain as $\hat{\mathcal{Y}} = \{\hat{\mathbf{Y}}_h, \mathbf{X}\}$, and the representation domain as $\mathcal{S} = \{\mathbf{S}_h, \mathbf{S}_m\}$, the encoder of the network $E_\phi: \mathcal{Y} \to \mathcal{S}$, maps the input data to low-dimensional representations (latent variables on the Bottleneck hidden layer), i.e., $p_\phi(\mathcal{S}|\mathcal{Y})$ and the decoder $D_\psi: \mathcal{S} \to \hat{\mathcal{Y}}$ reconstructs the data from the representations, i.e., $p_\psi(\hat{\mathcal{Y}}|\mathcal{S})$. Note that the bottleneck hidden layer $\mathcal{S}$ behaves as the representation layer that reflects the spatial information, and the weights $\psi$ of the decoder $D_\psi$ serve as $D_h$ in Eq. (1), respectively. This correspondence is further elaborated below.

Taking the procedure of training LR HSI as an example. The LR HSI is reconstructed by $\hat{\mathbf{Y}}_h = \mathbf{D}_\psi(\mathbf{S}_h)$, where $\mathbf{S}_h = \mathbf{E}_\phi(\mathbf{Y}_h)$. Since $\mathbf{Y}_h$ carries the high-resolution spectral information, to better extract the spectral basis, part of the network should simulate the prior relationship described in Eq. (1). That is, the representation layer $\mathbf{S}_h$ acts as the proportional coefficients and the weights $\psi$ of the decoder correspond to the spectral basis $D_h$ in Eq. (1). Therefore, in the network structure, we define $\psi = \mathbf{W}_1\mathbf{W}_2...\mathbf{W}_k = D_h$ with identity activation function without bias, where $\mathbf{W}_k$ denotes the weights in the $k$th layer. In this way, $D_h$ preserves the spectral information of LR HSI, and the latent variables $\mathbf{S}_h$ preserves the spatial information effectively. More implementation details will be described in Sec. IV-B.

Eventually, the desired HR HSI is generated directly by $\mathbf{X} = \mathbf{S}_m D_h$. Note that the dashed lines in Fig. 3 show the path of back-propagation which will be elaborated in Sec. IV-C.

<figure><img src="image_5.png" /><figcaption>Fig. 4. Details of the encoder-decoder structure.</figcaption></figure>

## B. Mutual Dirichlet Network with Collaborative Constraint

To extract better spectral information and naturally incorporate the physical requirements of spatial information, i.e., non-negative and sum-to-one, the representations $\mathcal{S}$ are encouraged to follow a Dirichlet distribution. In addition, the network should have the ability to learn the correlated and optimized representations generated from the encoder $E_\phi$ for both modalities. Thus, in the network design, we maximize the mutual information (MI) between the representations of LR HSI, $\mathbf{S}_h$, and HR MSI, $\mathbf{S}_m$, by maximizing the MI between the input images and their own representations. To further reduce the spectral distortion, the collaborative $l_{2,1}$ loss is incorporated into the network instead of the traditional $l_2$ reconstruction loss. The detailed encoder-decoder structure and the MI structure are shown in Fig. 4 and Fig. 5, respectively.

1) **Dirichlet Structure:** To generate representations with Dirichlet distribution, we incorporate the stick-breaking structure between the encoder and representation layers. The stick-breaking process was first proposed by Sethuraman [60] back in 1994. It can be illustrated as breaking a unit-length stick into $c$ pieces, the length of which follows a Dirichlet distribution. Nalisnick and Smyth, and Qu et al. successfully coupled the expressiveness of networks with the Bayesian nonparametric model through a stick-breaking process [21], [61]. Here, we follow the work of [21], [61], which draw the samples of $\mathcal{S}$ from Kumaraswamy distribution [62].

The stick-breaking process is integrated into the network between the encoder $E_\phi$ and the decoder $D_\psi$, as shown in Fig. 3. Assuming that the generated representation row vector is denoted as $\mathbf{s}_i = \{s_{ij}\}_{1\le j\le c}$, we have $0 \le s_{ij} \le 1$, and $\sum_{j=1}^c s_{ij} = 1$. Each variable $s_{ij}$ can be defined as

$$ s_{ij} = \begin{cases} v_{i1} & \text{for } j = 1 \\ v_{ij} \prod_{k<j} (1 - v_{ik}) & \text{for } j > 1, \end{cases} \quad (5) $$

where $v_{ik} \sim \text{Beta}(u, \alpha, \beta)$. Since it is difficult to draw samples directly from the Beta distribution, we draw samples from the inverse transform of Kumaraswamy distribution. The benefit of the Kumaraswamy distribution is that it has a closed-form CDF, and it is equivalent to the Beta distribution when $\alpha = 1$ or $\beta = 1$. Let $\alpha = 1$, we have

$$ v_{ik} \sim 1 - (1 - u_{ik}^{\frac{1}{\beta_i}}). \quad (6) $$

Both parameters $u_{ik}$ and $\beta_i$ are learned through the network for each row vector as illustrated in Fig. 3. Because $\beta > 0$, a softplus is adopted as the activation function [63] at the $\beta$ layer. Similarly, a sigmoid [64] is used to map $u$ into $(0, 1)$ range at the $u$ layer. Due to the fact that the spectral
signatures of data are different for each image pair, the network only trains one group of data, i.e., LR HSI $\mathbf{Y}_h$ and HR MSI $\mathbf{Y}_m$, to reconstruct its own HR HSI $\mathbf{X}$. Therefore, to increase the representation power of the network, the encoder of the network is densely connected, i.e., each layer is fully connected with all its subsequent layers [65].

2) Mutual Dirichlet Network: Before further describing the details of the network, we first explain the reason that motivates this design. Given unregistered multi-modalities LR HSI, $\mathbf{Y}_h$ and HR MSI, $\mathbf{Y}_m$, and the desired HR HSI, $\mathbf{X}$, each pixel of which indicates the mixed spectral reflection of the captured area. The overlapped region of the three modalities is defined by $\mathcal{C}$. Ideally, each pixel in the overlapped region of these three modalities should possess the same spectral signatures. In addition, the corresponding proportional coefficients of $\mathbf{X}$ and $\mathbf{Y}_m$ should be the same for a given pixel within $\mathcal{C}$. Since $\mathbf{Y}_h$ is a down-sampling and transformed version of $\mathbf{X}$, its proportional coefficients (representations) should follow the same pattern as that of $\mathbf{X}$ and $\mathbf{Y}_m$, i.e., $\mathbf{S}_h$ and $\mathbf{S}_m$ should be highly correlated although with different resolution. One example is shown in Fig. 1. Therefore, to generate HR HSI with low spectral distortion, it is necessary to encourage the representations $\mathbf{S}_h$ and $\mathbf{S}_m$ to follow similar patterns. However, traditional constraints like correlation may not work properly, because the input LR HSI and HR MSI are not registered with each other and the mapping function $\mathbf{E}_\phi$, between the input $\mathcal{Y}$ and the representations $\mathbf{S}$, holds the non-linear property. Therefore, we introduce MI, which captures the non-linear statistical dependencies between variables [66], to reinforce the representations of LR HSI and HR MSI to follow similar patterns with statistics.

Mutual information has been widely used for multi-modality registrations [67], [68]. It is a Shannon-entropy-based measurement of mutual independence between two random variables, e.g., $\mathbf{S}_h$ and $\mathbf{S}_m$. The mutual information $\mathcal{I}(\mathbf{S}_h; \mathbf{S}_m)$ measures how much uncertainty of one variable ($\mathbf{S}_h$ or $\mathbf{S}_m$) is reduced given the other variable ($\mathbf{S}_m$ or $\mathbf{S}_h$). Mathematically, it is defined as

$$ \begin{aligned} \mathcal{I}(\mathbf{S}_h; \mathbf{S}_m) &= H(\mathbf{S}_h) - H(\mathbf{S}_h|\mathbf{S}_m) \\ &= \int_{\mathbf{S}_h \times \mathbf{S}_m} \log \frac{d\mathbb{P}_{\mathbf{S}_h|\mathbf{S}_m}}{d\mathbb{P}_{\mathbf{S}_h} \otimes d\mathbb{P}_{\mathbf{S}_m}} d\mathbb{P}_{\mathbf{S}_h|\mathbf{S}_m} \end{aligned} \quad (7) $$

where $H$ indicates the Shannon entropy, $H(\mathbf{S}_h|\mathbf{S}_m)$ is the conditional entropy of $\mathbf{S}_h$ given $\mathbf{S}_m$. $d\mathbb{P}_{\mathbf{S}_h|\mathbf{S}_m}$ is the joint probability distribution, and $\mathbb{P}_{\mathbf{S}_h}$, $\mathbb{P}_{\mathbf{S}_m}$ denote the marginals. Belghazi et al. [69] introduced an MI estimator, which allows neural network to estimate MI through back-propagation, by adopting the concept of Donsker-Varadhan representation [70].

In order to maximally preserve the spectral information of the reconstructed HR HSI, our goal is to encourage the two representations $\mathbf{S}_h$ and $\mathbf{S}_m$ to follow similar patterns by maximizing their MI, $\mathcal{I}(\mathbf{S}_h; \mathbf{S}_m)$, during the optimization procedure. Since $\mathbf{S}_h = \mathbf{E}_\phi(\mathbf{Y}_h)$ and $\mathbf{S}_m = \mathbf{E}_\phi(\mathbf{Y}_m)$, the MI can also be expressed as $\mathcal{I}(\mathbf{E}_\phi(\mathbf{Y}_h); \mathbf{E}_\phi(\mathbf{Y}_m))$. However, it is difficult to maximize such MI directly with neural networks, because the two modalities do not match with each other in our scenario. Therefore, we maximize the average MI between the representations and their own inputs, i.e., $\mathcal{I}(\mathbf{Y}_h, \mathbf{E}_\phi(\mathbf{Y}_h))$

<figure><img src="image_6.png" /><figcaption>Fig. 5. Details of the MI structure.</figcaption></figure>

and $\mathcal{I}(\mathbf{Y}_m, \mathbf{E}_\phi(\mathbf{Y}_m))$. The benefit of doing this is two-fold. First, by optimizing the encoder weights $\mathbf{E}_\phi$, it is able to greatly improve the quality of individual representations [71]. Thus it helps the network to preserve the spectral and spatial information better. Second, since the multi-modalities, i.e., $\mathbf{Y}_h$ and $\mathbf{Y}_m$, are correlated, and the dependencies (MI) between the representations and multi-modalities are maximized, it also maximizes the MI, $\mathcal{I}(\mathbf{S}_h; \mathbf{S}_m)$, between different modalities, such that $\mathbf{S}_h$ and $\mathbf{S}_m$ are encouraged to follow similar patterns. Let's explain it with a toy example. We assume that both $\mathbf{Y}_h$ and $\mathbf{Y}_m$ cover the same material 'brick', the spectral pixel of which in the image pairs are denoted by $\mathbf{y}_h$ and $\mathbf{y}_m$, respectively, and $\tilde{\mathbf{y}}_h = \mathbf{y}_h \mathcal{R}$. $\tilde{\mathbf{y}}_h$, and $\mathbf{y}_m$ may not be identical to each other in real applications, but they are correlated and should possess similar spectral information. By maximizing the MI between the image and their representations, we are able to find a better representation $\mathbf{s}_h$ which reduces the uncertainty of $\tilde{\mathbf{y}}_h$ to a large extent, and also a better representation $\mathbf{s}_m$, which reduces the uncertainty of $\tilde{\mathbf{y}}_m$ to a large extent. Since $\tilde{\mathbf{y}}_h$ and $\mathbf{y}_m$ are similar, $\mathbf{s}_m$ and $\mathbf{s}_h$ should also be similar. In this way, the MI can regularize the solution space, such that $\mathbf{S}_h$ and $\mathbf{S}_m$ have similar patterns.

Taking $\mathcal{I}(\mathbf{Y}_h, \mathbf{E}_\phi(\mathbf{Y}_h))$ as an example. It is equivalent to Kullback-Leibler (KL) divergence [69] between the joint distribution $\mathbb{P}_{\mathbf{Y}_h|\mathbf{E}_\phi(\mathbf{Y}_h)}$ and the product of the marginals $\mathbb{P}_{\mathbf{Y}_h} \otimes \mathbb{P}_{\mathbf{E}_\phi(\mathbf{Y}_h)}$. Let $\mathbb{P} = \mathbb{P}_{\mathbf{Y}_h|\mathbf{E}_\phi(\mathbf{Y}_h)}$ and $\mathbb{Q} = \mathbb{P}_{\mathbf{Y}_h} \otimes \mathbb{P}_{\mathbf{E}_\phi(\mathbf{Y}_h)}$, we can further express MI as

$$ \mathcal{I}(\mathbf{Y}_h, \mathbf{E}_\phi(\mathbf{Y}_h)) = \mathbb{E}_{\mathbb{P}}[\log \frac{d\mathbb{P}}{d\mathbb{Q}}] = D_{KL}(\mathbb{P}||\mathbb{Q}) \quad (8) $$

Such MI can be maximized by maximizing the KL-divergence's lower bound based on Donsker-Varadhan (DV) representation [70]. Since we do not need to calculate the exact MI, we introduce an alternative lower bound based on Jensen-Shannon which works better than the DV-based objective function [71].

In the network design, an additional network $\mathcal{T}_w: \mathcal{Y} \times \mathcal{S} \to \mathbb{R}$ is built with two fully-connected layers, whose weights are denoted as $w$. During the training procedure, the raw image and the extracted representations are stacked and fed into the network as shown in Fig. 5. Then the estimator can be defined as

$$ \mathcal{I}_{\phi,w}(\mathbf{Y}_h, \mathbf{E}_\phi(\mathbf{Y}_h)) := \mathbb{E}_{\mathbb{P}}[-sp(-\mathcal{T}_{w,\phi}(\mathbf{Y}_h, \mathbf{E}_\phi(\mathbf{Y}_h))] \quad (9) $$

where $sp(x) = \log(1 + e^x)$. Note that we ignore the negative samples in DV-based objective function [71], which are usually generated by shuffling the input data. Because it is unstable to train the network with random shifting input data given only two input data pairs. Since both $\mathbf{E}_\phi$ and $\mathcal{T}_w$ are used to find the
optimal representations, they are updated together. Combined with the MSI MI, the objective function is defined as

$$ \mathcal{L}_{\mathcal{I}}(\phi, w) = \mathcal{I}_{\phi,w}(\mathbf{Y}_h, \mathbf{E}_{\phi}(\mathbf{Y}_h)) + \mathcal{I}_{\phi,w}(\mathbf{Y}_m, \mathbf{E}_{\phi}(\mathbf{Y}_m)) \quad (10) $$

Since the encoder $\mathbf{E}_{\phi}$ and the estimation network of MI $\mathcal{T}_w$ for both LR HSI and HR MSI share the same weights $\phi$ and $w$, their optimized representations follow similar patterns. More optimization details are described in Sec. IV-C.

In order to extract better spectral information, we adopt the collaborative reconstruction loss with $l_{2,1}$ norm [72] instead of traditional $l_2$ norm for both LR HSI and HR MSI. The objective function for $l_{2,1}$ loss is defined as

$$ \mathcal{L}_{2,1}(\phi, \psi) = \|D_{\psi}(E_{\phi}(\mathbf{Y}_h)) - \mathbf{Y}_h\|_{2,1} + \|D_{\psi}(E_{\phi}(\mathbf{Y}_m)) - \mathbf{Y}_m\|_{2,1} \quad (11) $$

where $\|X\|_{2,1} = \sum_{i=1}^{m} \sqrt{\sum_{j=1}^{n} X_{i,j}^2}$. $l_{2,1}$ norm can be treated as the sequential application of the $l_2$ norm on each pixel vector, followed by the $l_1$ norm on the image to enforce the reconstruction errors of the entire image to be sparse, that is, most of the reconstruction errors of individual pixels to be zero, such that the individual pixels would be reconstructed as accurately as possible. In this way, it extracts better spectral information and further reduces the spectral distortion.

### C. Optimization and Implementation Details

The objective functions of the proposed network architecture can then be expressed as:

$$ \mathcal{L}(\phi, \psi, w) = \mathcal{L}_{2,1}(\phi, \psi) - \lambda \mathcal{L}_{\mathcal{I}}(\phi, w) + \mu \|\psi\|_{F}^{2} \quad (12) $$

where $l_2$ norm is applied on the decoder weights $\psi$ to prevent over-fitting. $\lambda$ and $\mu$ are the parameters that balance the trade-off between reconstruction error, negative of mutual information and weight loss, respectively.

Before feeding into the network, the spectral vectors in LR HSI and HR MSI are transformed to zero-mean vectors by reducing the vector mean of their own image. Since the spectral information of MSI has been compressed too much (e.g., HSI has 31 bands, but MSI has 3 bands), the decoder of the network is only updated by LR HSI data to stabilize the network. The number of the input nodes is equal to the band number of HR MSI $l$. LR HSI $\mathbf{Y}_h$ is projected into a $l$ dimensional space by $\tilde{\mathbf{Y}}_h = \mathbf{Y}_h \mathcal{R}$ before feeding into the network, while HR MSI is directly fed into the network. The number of the output nodes is chosen based on the band number of LR HSI $L$. When the input of the network is $\mathbf{Y}_h$, the output of the decoder is $\tilde{\mathbf{Y}}_h$. When the input of the network is $\mathbf{Y}_m$, the reconstructed $\tilde{\mathbf{Y}}_m$ is generated by multiplying the output of the decoder with fixed weights $\mathcal{R}$.

The encoder-decoder is constructed with fully-connected layers and the detailed structure is shown in Fig. 4. The input of the encoder has $l$ neurons carrying each pixel of the image, which is densely connected by stacking with all its subsequent layers. Let's take $l = 8$ as an example, the input layer has 8 neurons, and we assume that the second and the third layers have 3 neurons, respectively. The input layer is passed to the second layer by stacking the first layer on top of the second

layer. Then the stacked layer is passed to the third layer by stacking 11 neurons on top of the third layer. In this way, the encoder is densely connected. The layer $v$ is drawn with Eq. (6) given layer $u$ and layer $\beta$, which are learned by back-propagation. $\beta$ has only one node, which is learned by a two-layer densely-connected fully-connected neural network. It denotes the distribution parameter of each pixel. $u$ has 15 nodes, which are learned by a four-layer densely-connected neural-network. The representation layer $S$ with 15 nodes is constructed with $v$ and $\beta$, according to Eq. (5). The decoder has two fully-connected layers. The number of nodes and the activation functions for different layers are shown in Table I.

<table><caption>TABLE I
THE NUMBER OF LAYERS AND NODES IN THE PROPOSED NETWORK.</caption><thead><tr><th></th><th>u/β encoder</th><th>u/β/v</th><th>T<sub>w</sub></th><th>decoder</th></tr></thead><tbody><tr><td>#layers</td><td>4/2</td><td>1/1/1</td><td>2</td><td>2</td></tr><tr><td>#nodes</td><td>[3,3,3,3]/[3,3]</td><td>15/1/15</td><td>[18,1]</td><td>[15,15]</td></tr><tr><td>activation</td><td>linear</td><td>sigmoid/softplus/linear</td><td>sigmoid linear</td><td></td></tr></tbody></table>

The training is done in an unsupervised fashion without ground truth HR HSI. Given multi-modalities LR HSI and HR MSI, the network is optimized with back-propagation to extract their correlated spectral bases and representations, as illustrated in Fig. 3 with red-dashed lines. The training process stops when the reconstruction error of the network does not decrease anymore. Then we can feed the HR MSI into the trained network, and obtain the reconstructed HR HSI, $\mathbf{X}$, from the output of the decoder.

## V. EXPERIMENTS AND RESULTS

### A. Datesets

The proposed $u^2$-MDN has been extensively evaluated with two widely used benchmark datasets, CAVE [73] and Harvard [1], and five remote sensing datasets, Hyperspec Chikusei, CASI University of Houston, ROSIS-3 University of Pavia, HYDICE Washington DC Mall [4] and real data without simulation, as summarized in Table II.

1) **Cave dataset:** The CAVE dataset consists of 32 HR HSI images and each of which has a dimension of 512 × 512 with 31 spectral bands taken within the wavelength range 400–700 nm at an interval of 10 nm.

2) **Harvard dataset:** The Harvard dataset includes 50 HR HSI images with both indoor and outdoor scenes. The images are cropped to 1024 × 1024, with 31 bands taken at an interval of 10 nm within the wavelength range of 420–720 nm.

3) **Hyperspec Chikusei dataset:** The dataset was taken by Headwall's Hyperspec-VNIR-C sensor over Chikusei, Ibaraki, Japan. The image has a ground sampling distance (GSD) of 2.5 m and was cropped to 540 × 420 with 128 bands, covering the wavelength range from 363 to 1018 nm. Please refer to [74] for more details.

4) **University of Houston dataset:** This dataset was acquired by ITRES CASI-1500 sensor over the University of Houston campus with a GSD of 2.5 m [75]. It was cropped to 320 × 540 with 144 bands taken within the wavelength range 364–1046 nm.
5) University of Pavia dataset: The dataset was taken by the reflective optics spectrographic imaging system (ROSIS-3) sensor over the University of Pavia, Italy, with a GSD of 1.3 m. It was cropped to 560 × 320 with 103 spectral bands taken within the wavelength range 430–830 nm.

6) Washington DC Mall dataset: The dataset was acquired by the hyperspectral digital imagery collection experiment (HYDICE) sensor over the Mall in Washington DC, USA at a GSD of 2.5 m. The image was cropped to 420 × 300 with 191 bands covering the wavelength range from 400 to 2500 nm.

7) Real dataset without simulation: The LR HSI over the Cuprite mining district, Nevada, US, was acquired by Hyperion with a GSD of 30 m, the image size of which is 100×153 with 167 bands taken within the wavelength range from 426 to 2355 nm. The HR MSI is the SWIR data of WorldView3 with a GSD of 7.5 m, the image size of which is 460×670 with 8 bands covering the wavelength range from 1209 to 2329 nm. Both rigid and nonrigid deformation exist as shown in Figs. 14a and 14b.

## B. Experimental Setup

For real applications, the mis-registration of two modalities is crucial for HSI-SR [20], [22], [47]. To demonstrate how misregistration would influence the performance of HSI-SR, we conduct two groups of experiments to evaluate the various approaches, i.e., the experiments on well-registered image pairs, and on unregistered image pairs. By conducting experiments in these two scenarios, we intend to show that misregistration would influence the performance of HSI-SR significantly. Therefore, it is very important to develop algorithms that can directly work on unregistered image pairs.

The well-registered image pairs are generated in two different ways following the widely-used protocols for benchmark datasets [12], [42], [76] and the Walds protocol [4], [77] for remote sensing datasets.

* For benchmark HSI datasets, CAVE [73] and Harvard [1], the image pairs are generated with the extreme Super-Resolution (SR) ratio of 32, where the LR HSI $Y_h$ is obtained by averaging the HR HSI over 32 × 32 disjoint blocks. The HR MSI with 3 bands are generated by multiplying the HR HSI with the given spectral response matrix $R$ of Nikon D700 [11], [12], [21]. Note that we adopt this setting because it is the same protocol used by state-of-the-art methods [11], [12], [42], [76] on general hyperspectral images. In addition, for remote sensing applications, the scale difference can even be 25 [56] and 30 [57]. With such settings, we are able to evaluate the proposed method in extreme scenarios.

* For remote sensing datasets, the image pairs are simulated with the Walds protocol [77], where the LR HSI is generated by applying a Gaussian filter with its full width at half maximum (FWHM) equal to the SR ratio, to match a plausible system modulation transfer function (MTF) [4], [17], [31]. The MSI is generated by degrading the HR HSI in the spectral domain using MSI spectral reflection functions (SRFs) from different sensors as filters. The

datasets are listed in Table II. Please refer to [4] for more details. Note that since the scales are different between the real LR HSI and HR MSI for different sensors [4], [56], [57], the SR ratio is set to 4, 5, 6 and 8, to evaluate the robustness of the proposed method. The noise is added to the image with a signal-to-noise-ratio (SNR) of 30 dB in all bands.

The unregistered image pairs are generated in the same way as that of the well-registered image pairs, except that the LR HSI images are further distorted with rigid or nonrigid deformations.

* For benchmark HSI datasets, CAVE [73] and Harvard [1], it is easier to introduce rigid deformation. Thus, the LR HSI is further rotated with 5° and cropped by 15% of its surrounding pixels, e.g., for images in the CAVE dataset, 39,322 pixels of the MSI are not covered in the LR HSI; and for images in the Harvard dataset, 157,290 pixels of the MSI are not covered in the LR HSI.

* For remote sensing datasets, it is usually unavoidable to introduce nonrigid deformation [22]. Thus, following the protocol in [47], [78], the nonrigid distortion is emulated by introducing random shifts in pixels.

* For real data, the LR HSI is directly captured from Hyperion and the HR MSI is captured from WorldView3. Both rigid and nonrigid deformations exist as shown in Figs. 14a and 14b.

<table><caption>TABLE II
DATASET PAIRS FROM DIFFERENT SENSORS USED IN THE EXPERIMENTS.</caption><thead><tr><td>Dataset</td><td>HSI sensor</td><td>MSI sensor</td><td>SR Ratio</td></tr></thead><tbody><tr><td>CAVE</td><td>Apogee Alta U260</td><td>Nikon</td><td>32</td></tr><tr><td>Harvard</td><td>Nuance FX</td><td>Nikon</td><td>32</td></tr><tr><td>Chikusei</td><td>Hyperspec</td><td>WorldView2</td><td>6</td></tr><tr><td>Houston</td><td>CASI</td><td>Sentinel-2</td><td>5</td></tr><tr><td>Pavia</td><td>ROSIS-3</td><td>QuickBird</td><td>8</td></tr><tr><td>Washington</td><td>HYDICE</td><td>QuickBird</td><td>4</td></tr><tr><td>Real data</td><td>Hyperion</td><td>WorldView3</td><td>4</td></tr></tbody></table>

The results of the proposed method on individual images in Fig. 6 are compared with nine state-of-the-art methods, including traditional methods such as CS-based GSA [15] and MRA-based SFIM [33], matrix factorization based methods such as CNMF [39] and Lanaras' CSU [12], Bayesian-based methods such as HySure [38], sparse-coding based methods such as NSSR [76], tensor-based method [42], the integrated registration and fusion method [47], and the uSDN method [21] that belong to different categories of HSI-SR. These methods also reported the best performance [12], [17], [21], with the original code made available by the authors. Note that the proposed $u^2$-MDN is unsupervised, i.e., the HR HSI is not available during the training procedure. Thus, for a fair comparison, only unsupervised methods are included in the experiments. The average results on the datasets are also reported to evaluate the robustness of the proposed method.

For rigid deformation, since the resolution of HSI does not match that of the degraded MSI, i.e., there exists large displacement between two modalities, only five methods may reconstruct HR HSI from unregistered images without large errors. Thus, the proposed method is compared with these five state-of-the-art methods, i.e., GSA [15], SFIM [33],
CNMF [39], NSSR [76], and the integrated registration and fusion method [47] on unregistered image pairs. Note that, as discussed in Sec. III, in order to work on unregistered image pairs, the LR HSI should include all the spectral bases of HR MSI. For the CAVE and Harvard datasets, not all the image pairs meet this requirement after rotation and cropping. Thus, we choose seven commonly used image pairs from the benchmark dataset, where the LR HSI includes all the spectral bases of HR MSI even after rotation and cropping. The chosen image pairs are shown in Fig. 6. The remote sensing images are shown in Fig. 7.

<figure><img src="image_7.png" /><figcaption>Fig. 6. The HR MSI of individual test images from the CAVE [73] (top row) and Harvard [1] (bottom row) datasets.</figcaption></figure>

### C. Evaluation Metrics

For quantitative comparison, the error relative globale adimensionnelle de synthèse (ERGAS), the peak signal-to-noise ratio (PSNR), and the spectral angle mapper (SAM) are applied to evaluate the quality of the reconstructed HSI.

ERGAS provides a measurement of the band-wise normalized root of mean square error (RMSE) between the reference HSI, $\mathbf{X}$, and the reconstructed HSI, $\hat{\mathbf{X}}$, with the best value at 0 [77]. It is defined as

$$ \text{ERGAS}(\mathbf{X}, \hat{\mathbf{X}}) = \frac{100}{sr} \sqrt{\frac{1}{L} \sum_{i=1}^{L} \frac{\text{mean} \|\mathbf{X}_i - \hat{\mathbf{X}}_i\|_2^2}{(\text{mean}\mathbf{X}_i)^2}}, \quad (13) $$

where sr denotes the sr factor between the HR MSI and LR HSI, L denotes the number of spectral bands of the reconstructed $\hat{\mathbf{X}}$.

PSNR is the average ratio between the maximum power of the image and the power of the residual errors in all the spectral bands. A larger PSNR indicates a higher spatial quality of the reconstructed HSI. For each image band of HSI, the PSNR is defined as

$$ \text{PSNR}(\mathbf{X}_i, \hat{\mathbf{X}}_i) = 10 \cdot \log_{10} \left( \frac{\max(\mathbf{X}_i)^2}{\text{mean} \|\mathbf{X}_i - \hat{\mathbf{X}}_i\|_2^2} \right) \quad (14) $$

SAM [79] is commonly used to quantify the spectral distortion of the reconstructed HSI. The larger the SAM, the worse

<figure><img src="image_8.png" /><figcaption>Fig. 7. Color composite of the remote sensing datasets from [4]. The reference HR HSI of the (a) Chikusei, (b) Houston, (c) Pavia and (d) Washington datasets.</figcaption></figure>

the spectral distortion of the reconstructed HSI. For each HSI pixel $\hat{\mathbf{X}}_j$, the SAM is defined as

$$ \text{SAM}(\mathbf{X}_j, \hat{\mathbf{X}}_j) = \arccos \left( \frac{\mathbf{X}_j^T \hat{\mathbf{X}}_j}{\|\mathbf{X}_j\|_2 \|\hat{\mathbf{X}}_j\|_2} \right) \quad (15) $$

The global SAM is estimated by averaging the SAM over all the pixels in the entire image.

### D. Experimental Results on Registered Image Pairs

For a fair comparison, we first perform experiments on the general case when LR HSI and HR MSI are well registered. Table III show the experimental results of 7 groups of commonly benchmarked images from the CAVE and Harvard datasets [11], [12], [21], [76]. Table IV show the experimental results of the remote sensing images. The average results of the datasets are shown in Table V. Note that, in order to show how the method works in different scenarios, the data are not normalized for evaluation. Since the intensities of the Harvard dataset are quite small, the ERGAS of the reconstructed images is generally smaller than those of the CAVE dataset and remote sensing dataset.

We observe that CS-based GSA [15] is stable on both the benchmarked and remote sensing datasets. However, it could not preserve the spectral information well especially on the benchmarked datasets. Matrix-factorization-based CSU [12] works better than CNMF [39] on the benchmarked CAVE and Harvard datasets. However, its performance is worse than that
<table><caption>TABLE III
BENCHMARKED RESULTS IN TERMS OF ERGAS (E), PSNR (P) AND SAM (S) ON WELL-REGISTERED IMAGE PAIRS.</caption><thead><tr><th rowspan="3">Methods</th><th colspan="9">CAVE</th><th colspan="9">Harvard</th></tr><tr><th colspan="3">balloon</th><th colspan="3">cloth</th><th colspan="3">pompoms</th><th colspan="3">spool</th><th colspan="3">img1</th><th colspan="3">imgb5</th><th colspan="3">imgc5</th></tr><tr><th>E</th><th>P</th><th>S</th><th>E</th><th>P</th><th>S</th><th>E</th><th>P</th><th>S</th><th>E</th><th>P</th><th>S</th><th>E</th><th>P</th><th>S</th><th>E</th><th>P</th><th>S</th><th>E</th><th>P</th><th>S</th></tr></thead><tbody><tr><td>GSA</td><td>0.19</td><td>41.89</td><td>4.07</td><td>0.40</td><td>32.51</td><td>5.95</td><td>0.37</td><td>34.78</td><td>7.39</td><td>0.41</td><td>39.61</td><td>9.53</td><td>0.12</td><td>40.41</td><td>2.19</td><td>0.16</td><td>39.07</td><td>2.19</td><td>0.12</td><td>38.82</td><td><b>1.67</b></td></tr><tr><td>SFIM</td><td>0.59</td><td>33.52</td><td>8.45</td><td>0.54</td><td>30.59</td><td>5.25</td><td>3.76</td><td>25.39</td><td>11.89</td><td>2.93</td><td>28.63</td><td>19.71</td><td>0.23</td><td>32.62</td><td>2.10</td><td>0.29</td><td>33.15</td><td>3.52</td><td>0.23</td><td>35.62</td><td>2.84</td></tr><tr><td>CNMF</td><td>0.26</td><td>39.27</td><td>9.71</td><td>0.54</td><td>30.52</td><td>6.55</td><td>0.31</td><td>35.45</td><td>6.32</td><td>0.54</td><td>37.28</td><td>16.77</td><td>0.15</td><td>37.25</td><td>2.86</td><td>0.17</td><td>39.06</td><td>2.14</td><td>0.13</td><td>38.49</td><td>2.64</td></tr><tr><td>CSU</td><td>0.19</td><td>41.52</td><td>4.68</td><td>0.40</td><td>33.47</td><td>5.52</td><td>0.28</td><td>36.81</td><td>6.01</td><td>0.45</td><td>39.64</td><td>6.84</td><td>0.12</td><td>39.12</td><td>2.30</td><td>0.18</td><td>39.01</td><td>2.37</td><td>0.12</td><td>39.05</td><td>2.38</td></tr><tr><td>HySure</td><td>0.34</td><td>37.08</td><td>9.92</td><td>0.53</td><td>30.22</td><td>7.13</td><td>0.52</td><td>31.68</td><td>10.97</td><td>0.55</td><td>37.47</td><td>15.54</td><td>0.18</td><td>35.82</td><td>4.27</td><td>0.34</td><td>35.52</td><td>3.45</td><td>0.19</td><td>36.75</td><td>2.34</td></tr><tr><td>NSSR</td><td>0.16</td><td>43.2</td><td>3.35</td><td>0.31</td><td>33.3</td><td>4.58</td><td>0.26</td><td>37.71</td><td>5.31</td><td>0.45</td><td>39.41</td><td>6.91</td><td>0.14</td><td>39.91</td><td>2.24</td><td>0.17</td><td>39.12</td><td>2.17</td><td>0.12</td><td>38.87</td><td>1.87</td></tr><tr><td>CSTF</td><td><b>0.14</b></td><td><b>44.71</b></td><td>3.97</td><td>0.39</td><td>32.51</td><td>5.25</td><td>0.27</td><td>36.72</td><td>6.09</td><td>0.38</td><td><b>42.06</b></td><td>8.61</td><td>0.21</td><td>33.73</td><td>2.77</td><td>0.25</td><td>34.98</td><td>2.46</td><td>0.22</td><td>32.48</td><td>1.96</td></tr><tr><td>Integrated uSDN</td><td>0.28</td><td>37.75</td><td>2.64</td><td>1.47</td><td>21.55</td><td>8.73</td><td>0.52</td><td>30.29</td><td>5.99</td><td>1.03</td><td>30.94</td><td>6.77</td><td>0.32</td><td>29.81</td><td>2.68</td><td>0.63</td><td>26.29</td><td>2.31</td><td>0.27</td><td>30.47</td><td>1.79</td></tr><tr><td>u<sup>2</sup>-MDN</td><td>0.20</td><td>41.54</td><td>4.56</td><td>0.35</td><td>33.48</td><td><b>4.16</b></td><td>0.25</td><td>37.84</td><td>5.43</td><td>0.40</td><td>38.49</td><td>13.01</td><td>0.12</td><td>39.30</td><td>2.27</td><td>0.16</td><td>39.72</td><td>2.10</td><td><b>0.11</b></td><td>39.12</td><td>2.58</td></tr><tr><td></td><td>0.16</td><td>43.59</td><td><b>1.93</b></td><td><b>0.30</b></td><td><b>34.85</b></td><td>4.31</td><td><b>0.19</b></td><td><b>39.12</b></td><td><b>3.46</b></td><td><b>0.37</b></td><td>40.08</td><td><b>4.47</b></td><td><b>0.11</b></td><td><b>40.97</b></td><td><b>2.06</b></td><td><b>0.15</b></td><td><b>39.76</b></td><td><b>2.08</b></td><td><b>0.11</b></td><td><b>39.19</b></td><td>1.77</td></tr></tbody></table>

<table><caption>TABLE IV
REMOTE SENSING RESULTS IN TERMS OF ERGAS, PSNR AND SAM ON WELL-REGISTERED IMAGE PAIRS.</caption><thead><tr><th rowspan="2">Methods</th><th colspan="3">Chikusei</th><th colspan="3">Houston</th><th colspan="3">Pavia</th><th colspan="3">Washington</th></tr><tr><th>ERGAS</th><th>PSNR</th><th>SAM</th><th>ERGAS</th><th>PSNR</th><th>SAM</th><th>ERGAS</th><th>PSNR</th><th>SAM</th><th>ERGAS</th><th>PSNR</th><th>SAM</th></tr></thead><tbody><tr><td>GSA</td><td>1.432</td><td>42.1264</td><td>1.4478</td><td><b>2.7859</b></td><td>34.133</td><td>1.8443</td><td>1.0661</td><td>38.7949</td><td>3.5647</td><td>3.518</td><td>37.2308</td><td>2.187</td></tr><tr><td>SFIM</td><td>1.2284</td><td>47.4358</td><td><b>0.9379</b></td><td>2.9415</td><td>33.9958</td><td>0.9938</td><td>0.7274</td><td>42.9283</td><td>2.312</td><td>3.0356</td><td>39.2045</td><td>1.2382</td></tr><tr><td>CNMF</td><td>1.479</td><td>47.8427</td><td>1.1602</td><td>2.9896</td><td>33.1454</td><td>1.3882</td><td>0.7712</td><td>43.2417</td><td>2.3623</td><td><b>3.0341</b></td><td>39.1491</td><td>1.388</td></tr><tr><td>CSU</td><td>2.4705</td><td>35.8506</td><td>1.9208</td><td>3.2773</td><td>32.3793</td><td>2.0193</td><td>1.7283</td><td>33.9385</td><td>3.5754</td><td>4.2854</td><td>34.1841</td><td>1.9706</td></tr><tr><td>HySure</td><td><b>1.2216</b></td><td>48.7601</td><td>1.0934</td><td>2.9619</td><td><b>34.5328</b></td><td>1.7281</td><td>0.7767</td><td>43.2719</td><td>2.6094</td><td>3.3232</td><td>39.0</td><td>1.6808</td></tr><tr><td>NSSR</td><td>2.6427</td><td>33.5161</td><td>2.5263</td><td>4.7663</td><td>29.2931</td><td>5.3182</td><td>3.7068</td><td>28.8702</td><td>5.7786</td><td>9.1737</td><td>29.9297</td><td>4.0385</td></tr><tr><td>CSTF</td><td>1.9024</td><td>38.1548</td><td>1.7884</td><td>4.0207</td><td>29.5598</td><td>6.4</td><td>1.1877</td><td>37.3</td><td>4.0719</td><td>22.4659</td><td>20.4012</td><td>20.1433</td></tr><tr><td>Integrated uSDN</td><td>1.3854</td><td>43.4116</td><td>1.4104</td><td>4.0627</td><td>28.9168</td><td>3.9936</td><td>1.1773</td><td>37.9724</td><td>3.5066</td><td>5.8183</td><td>29.5106</td><td>4.3237</td></tr><tr><td>Proposed</td><td>1.7861</td><td>42.8702</td><td>1.3035</td><td>3.6198</td><td>32.5059</td><td>5.698</td><td>1.0221</td><td>39.4535</td><td>3.1874</td><td>6.7819</td><td>30.1769</td><td>5.3259</td></tr><tr><td></td><td>1.4717</td><td><b>50.2839</b></td><td>1.0578</td><td>2.8659</td><td>34.0584</td><td><b>0.8865</b></td><td><b>0.715</b></td><td><b>43.8022</b></td><td><b>2.3053</b></td><td>3.8988</td><td><b>39.2144</b></td><td><b>1.2298</b></td></tr></tbody></table>

<table><caption>TABLE V
THE AVERAGE ERGAS, PSNR AND SAM SCORES OVER WELL-REGISTERED BENCHMARKED AND REMOTE SENSING DATASETS.</caption><thead><tr><th rowspan="2">Methods</th><th colspan="3">CAVE</th><th colspan="3">Harvard</th><th colspan="3">Remote Sensing</th></tr><tr><th>ERGAS</th><th>PSNR</th><th>SAM</th><th>ERGAS</th><th>PSNR</th><th>SAM</th><th>ERGAS</th><th>PSNR</th><th>SAM</th></tr></thead><tbody><tr><td>GSA</td><td>0.34</td><td>37.2</td><td>6.74</td><td>0.13</td><td>39.43</td><td>2.02</td><td>2.2005</td><td>38.0713</td><td>2.261</td></tr><tr><td>SFIM</td><td>1.96</td><td>29.53</td><td>11.33</td><td>0.25</td><td>33.8</td><td>2.82</td><td><b>1.9832</b></td><td>40.8911</td><td>1.3705</td></tr><tr><td>CNMF</td><td>0.41</td><td>35.63</td><td>9.84</td><td>0.15</td><td>38.27</td><td>2.55</td><td>2.0685</td><td>40.8447</td><td>1.5747</td></tr><tr><td>CSU</td><td>0.33</td><td>37.86</td><td>5.76</td><td>0.14</td><td>39.06</td><td>2.35</td><td>2.9404</td><td>34.0881</td><td>2.3715</td></tr><tr><td>HySure</td><td>0.49</td><td>34.11</td><td>10.89</td><td>0.24</td><td>36.03</td><td>3.35</td><td>2.0709</td><td>41.3912</td><td>1.7779</td></tr><tr><td>NSSR</td><td>0.30</td><td>38.41</td><td>5.04</td><td>0.14</td><td>39.3</td><td>2.09</td><td>5.0724</td><td>30.4023</td><td>4.4154</td></tr><tr><td>CSTF</td><td>0.30</td><td>39.00</td><td>5.98</td><td>0.41</td><td>33.73</td><td>2.4</td><td>7.3942</td><td>31.3540</td><td>8.1009</td></tr><tr><td>Integrated uSDN</td><td>0.83</td><td>30.13</td><td>6.03</td><td>1.09</td><td>28.86</td><td>2.26</td><td>2.8672</td><td>33.3008</td><td>3.9413</td></tr><tr><td>u<sup>2</sup>-MDN</td><td>0.30</td><td>37.84</td><td>6.79</td><td>0.13</td><td>39.38</td><td>2.32</td><td>3.3025</td><td>36.2516</td><td>3.8787</td></tr><tr><td></td><td><b>0.26</b></td><td><b>39.41</b></td><td><b>3.54</b></td><td><b>0.12</b></td><td><b>39.97</b></td><td><b>1.97</b></td><td>2.2379</td><td><b>41.8397</b></td><td><b>1.3699</b></td></tr></tbody></table>

of CNMF on the remote sensing dataset, whose number of spectral bands is higher than that of the benchmarked dataset. MRA-based SFIM [33], Bayesian-based HySure [38] and the integrated fusion approach [47] could achieve relatively good performance on the remote sensing datasets, but their performance drops significantly on the benchmarked CAVE and Harvard datasets. On the contrary, sparse-coding-based NSSR [76] and tensor-based CSTF [42] could achieve much more competitive performance on the benchmarked datasets than on the remote sensing datasets. Note that for NSSR, the most effective step on the CAVE dataset is a post-processing step from [41], which actually degrades the performance on remote sensing datasets with more numbers of spectral bands. Thus, the post-processing step is disabled on the remote sensing datasets to improve the reconstruction accuracy. The tensor-based CSTF could achieve competitive results on the CAVE dataset, which has a redundant spatial structure. However, its performance drops on the remote sensing datasets with

less redundant spatial structure.

The deep-learning-based uSDN [21] preserves spectral information well on both the benchmarked and remote sensing datasets. However, it can only work on well-registered images due to its network design with angular difference regularization. Based on the average results shown in Table V, the proposed u²-MDN network powered by the mutual information and collaborative l₂₁ loss shows comparable, if not better, performance as compared to the state-of-the-art approaches in terms of ERGAS, PSNR, and SAM, and quite stable for different types of input images regardless of the number of spectral bands and SR ratios. In addition, it is very effective in preserving the spectral signature of the reconstructed HR HSI, showing much-improved performance, especially measured by SAM on the CAVE data. This further demonstrates the robustness of the proposed u²-MDN.
<table><caption>TABLE VI
RESULTS ON UNREGISTERED (RIGID DISTORTED) BENCHMARKED IMAGES IN TERMS OF ERGAS, PSNR AND SAM.</caption><thead><tr><th rowspan="3">Methods</th><th colspan="9">CAVE</th><th colspan="9">Harvard</th></tr><tr><th colspan="3">balloon</th><th colspan="3">cloth</th><th colspan="3">pompoms</th><th colspan="3">spool</th><th colspan="3">img1</th><th colspan="3">imgb5</th><th colspan="3">imgc5</th></tr><tr><th>E</th><th>P</th><th>S</th><th>E</th><th>P</th><th>S</th><th>E</th><th>P</th><th>S</th><th>E</th><th>P</th><th>S</th><th>E</th><th>P</th><th>S</th><th>E</th><th>P</th><th>S</th><th>E</th><th>P</th><th>S</th></tr></thead><tbody><tr><td>GSA</td><td>0.82</td><td>27.71</td><td>14.35</td><td>0.76</td><td>27.59</td><td>9.75</td><td>1.18</td><td>23.54</td><td>22.13</td><td>1.07</td><td>29.92</td><td>17.16</td><td>1.65</td><td>23.60</td><td>10.66</td><td>0.40</td><td>20.07</td><td>4.90</td><td>0.69</td><td>22.18</td><td>5.32</td></tr><tr><td>SFIM</td><td>1.51</td><td>22.47</td><td>12.69</td><td>1.01</td><td>24.74</td><td>9.65</td><td>1.82</td><td>19.50</td><td>14.89</td><td>1.88</td><td>25.30</td><td>21.02</td><td>1.33</td><td>17.41</td><td>3.28</td><td>0.68</td><td>25.38</td><td>4.44</td><td>0.89</td><td>19.93</td><td>3.96</td></tr><tr><td>CNMF</td><td>0.71</td><td>29.18</td><td>10.63</td><td>0.69</td><td>27.84</td><td>8.12</td><td>0.83</td><td>26.67</td><td>11.88</td><td>0.63</td><td>34.62</td><td>17.03</td><td>0.74</td><td>22.29</td><td>3.85</td><td>0.34</td><td>31.39</td><td>3.97</td><td>0.48</td><td>25.34</td><td>3.16</td></tr><tr><td>NSSR</td><td>0.52</td><td>32.59</td><td>8.07</td><td>0.72</td><td>27.16</td><td>8.05</td><td>0.76</td><td>27.45</td><td>10.22</td><td>1.03</td><td>32.80</td><td>15.94</td><td>0.61</td><td>25.83</td><td>5.29</td><td>0.50</td><td>29.72</td><td>6.80</td><td>0.35</td><td>28.66</td><td>2.64</td></tr><tr><td>Integrated<br>u<sup>2</sup>-MDN</td><td>1.14</td><td>24.68</td><td>9.34</td><td>1.68</td><td>19.84</td><td>11.25</td><td>1.82</td><td>19.38</td><td>17.60</td><td>1.65</td><td>29.81</td><td>13.60</td><td>1.00</td><td>19.80</td><td>4.05</td><td>0.68</td><td>25.40</td><td>3.24</td><td>0.87</td><td>20.25</td><td>2.40</td></tr><tr><td></td><td><strong>0.30</strong></td><td><strong>38.61</strong></td><td><strong>3.48</strong></td><td><strong>0.40</strong></td><td><strong>32.89</strong></td><td><strong>6.08</strong></td><td><strong>0.37</strong></td><td><strong>33.64</strong></td><td><strong>4.87</strong></td><td><strong>0.56</strong></td><td><strong>36.25</strong></td><td><strong>6.78</strong></td><td><strong>0.13</strong></td><td><strong>39.42</strong></td><td><strong>2.32</strong></td><td><strong>0.25</strong></td><td><strong>36.90</strong></td><td><strong>2.73</strong></td><td><strong>0.14</strong></td><td><strong>36.29</strong></td><td><strong>2.26</strong></td></tr></tbody></table>

<table><caption>TABLE VII
RESULTS ON UNREGISTERED (NONRIGID DISTORTED) REMOTE SENSING IMAGES IN TERMS OF ERGAS, PSNR AND SAM.</caption><thead><tr><th rowspan="2">Methods</th><th colspan="3">Chikusei</th><th colspan="3">Houston</th><th colspan="3">Pavia</th><th colspan="3">Washington</th></tr><tr><th>ERGAS</th><th>PSNR</th><th>SAM</th><th>ERGAS</th><th>PSNR</th><th>SAM</th><th>ERGAS</th><th>PSNR</th><th>SAM</th><th>ERGAS</th><th>PSNR</th><th>SAM</th></tr></thead><tbody><tr><td>GSA</td><td>4.4617</td><td>27.3028</td><td>7.6554</td><td>5.2167</td><td>27.6816</td><td>6.7581</td><td>3.5149</td><td>27.161</td><td>12.6803</td><td>7.1542</td><td>28.0469</td><td>7.8453</td></tr><tr><td>SFIM</td><td>6.8057</td><td>23.8192</td><td>5.8405</td><td>11.1495</td><td>22.5494</td><td>8.3767</td><td>5.2757</td><td>23.8248</td><td>9.2951</td><td>10.7509</td><td>23.9122</td><td>8.7243</td></tr><tr><td>CNMF</td><td>6.5318</td><td>25.5878</td><td>4.7582</td><td>5.6255</td><td>28.0016</td><td>5.5382</td><td>4.4232</td><td>25.2803</td><td>7.9551</td><td>25.8519</td><td>29.4765</td><td>5.9344</td></tr><tr><td>NSSR</td><td>5.8158</td><td>25.8475</td><td>5.4245</td><td>7.9435</td><td>23.7692</td><td>8.4696</td><td>5.0728</td><td>25.5778</td><td>8.2277</td><td>18.4162</td><td>22.3786</td><td>9.6114</td></tr><tr><td>Integrated<br>u<sup>2</sup>-MDN</td><td>8.9107</td><td>21.4644</td><td>7.5126</td><td>14.4247</td><td>18.7631</td><td>10.2734</td><td>7.4938</td><td>20.7214</td><td>11.9717</td><td>15.0896</td><td>20.9486</td><td>11.0361</td></tr><tr><td></td><td><strong>1.5843</strong></td><td><strong>45.4919</strong></td><td><strong>1.0899</strong></td><td><strong>2.9458</strong></td><td><strong>33.6519</strong></td><td><strong>1.0497</strong></td><td><strong>0.8898</strong></td><td><strong>40.4446</strong></td><td><strong>2.4371</strong></td><td><strong>4.0777</strong></td><td><strong>38.6361</strong></td><td><strong>1.2757</strong></td></tr></tbody></table>

## E. Experimental Results on Unregistered Image Pairs

In this section, two unregistered scenarios are studied, i.e., rigid distorted benchmarked datasets, and nonrigid distorted remote sensing datasets, as described in Sec. V-B. Note that, since the pixels in the HSI and MSI do not match with each other, the reconstruction errors are expected to be increased.

<figure><img src="image_9.png" /><figcaption>Fig. 8. The average PSNR of different wavelengths for the reconstructed HSI from the unregistered rigid distorted (a) CAVE dataset and (b) Harvard dataset, respectively.</figcaption></figure>

The performance of different methods on unregistered image pairs are reported in Tables VI and VII. Note that, only the methods that are able to work with unregistered image pairs are chosen in this group of experiments. Thus five state-of-the-art methods are compared in the tables. The traditional CS-based GSA [15] and MTF-based SFIM [33] fail in this scenario. This is because when the given two modalities are unregistered, the spatial details could not be directly added to improve the spatial resolution of LR HSI. The matrix-factorization-based CNMF and sparse-coding based NSSR are more robust than the traditional methods. However, their performance also drops for both benchmarked and remote sensing datasets. The reason is that the adopted predefined down-sampling function will introduce significant spectral distortion when the LR HSI and HR MSI are unregistered. The integrated fusion method could achieve good performance

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

<figure><img src="image_11.png" /><figcaption>Fig. 9. The PSNR of different wavelengths for the reconstructed HSI from unregistered nonrigid distorted from (a) Chikusei, (b) Houston, (c) Pavia and (d) Washington datasets, respectively.</figcaption></figure>

on remote sensing images with small distortion. However, its performance drops on images with large distortion. This is because the integrated fusion method performs registration before fusion, which may introduce additional distortion during optimization. The proposed u²-MDN is able to handle challenging scenarios much better than the state-of-the-art. The main reason that contributes to the success of the proposed approach is that, the network is able to extract the optimal and correlated spatial representations from two modalities through mutual information and collaborative loss. In this way, both the spatial and especially the spectral information are effectively
<table><caption>TABLE VIII
THE AVERAGE ERGAS, PSNR AND SAM SCORES OVER UNREGISTERED BENCHMARKED AND REMOTE SENSING DATASETS.</caption><thead><tr><th rowspan="2">Methods</th><th colspan="3">CAVE</th><th colspan="3">Harvard</th><th colspan="3">Remote Sensing</th></tr><tr><th>ERGAS</th><th>PSNR</th><th>SAM</th><th>ERGAS</th><th>PSNR</th><th>SAM</th><th>ERGAS</th><th>PSNR</th><th>SAM</th></tr></thead><tbody><tr><td>GSA</td><td>0.96</td><td>27.19</td><td>15.85</td><td>0.91</td><td>21.95</td><td>6.96</td><td>5.09</td><td>27.5481</td><td>8.7348</td></tr><tr><td>SFIM</td><td>1.56</td><td>23.00</td><td>14.56</td><td>0.97</td><td>20.91</td><td>3.89</td><td>8.50</td><td>23.5264</td><td>8.0592</td></tr><tr><td>CNMF</td><td>0.72</td><td>29.58</td><td>11.92</td><td>0.52</td><td>26.34</td><td>3.66</td><td>10.61</td><td>27.0866</td><td>6.0465</td></tr><tr><td>NSSR</td><td>0.76</td><td>30</td><td>10.57</td><td>0.49</td><td>28.07</td><td>4.91</td><td>9.31</td><td>24.3933</td><td>7.9333</td></tr><tr><td>Integrated</td><td>1.57</td><td>23.43</td><td>12.95</td><td>0.85</td><td>21.82</td><td>3.23</td><td>11.48</td><td>20.4744</td><td>10.1985</td></tr><tr><td>u<sup>2</sup>-MDN</td><td><b>0.41</b></td><td><b>35.35</b></td><td><b>5.30</b></td><td><b>0.17</b></td><td><b>37.54</b></td><td><b>2.44</b></td><td><b>2.3744</b></td><td><b>39.5561</b></td><td><b>1.4631</b></td></tr></tbody></table>

preserved. This demonstrates the representation capacity of the proposed structure.

To demonstrate the reconstruction performance in different spectral bands, the average PSNR of the benchmarked datasets on each wavelength is shown in Fig. 8. Since the numbers of spectral bands of the remote sensing datasets are different, we show their individual PSNR on each band in Fig. 9. We can observe that, regardless of the type of the datasets, the proposed method consistently outperforms the other methods for all the spectral bands on unregistered image pairs.

To visualize the reconstructed results for unregistered image pairs, we show the color composition of the reconstructed HR HSI in Figs. 10-13, among which Figs. 10, 11 demonstrate the results of the rigid distorted image pairs, while Figs. 12, 13 demonstrate the results of the nonrigid distorted image pairs. The first column of each figure presents the reference HR HSI, the distorted LR HSI, and the original LR HSI in (a), (h), and (o), respectively. The first through third rows show the reconstructed images, the absolute difference, and the spectral map of the results from different methods. We can observe that most approaches could not handle unregistered images pairs with large displacement well. The reconstructed results from SFIM have some blocking artifacts in most scenarios. The integrated fusion method has some smear effects on the reconstructed images due to the large displacement as shown in Fig. 10f and 13f. NSSR fails on the remote sensing datasets as shown in Figs. 12e and 13e, but it suffers relatively smaller spatial distortion on the benchmarked datasets. GSA could produce clear reconstructed images in most cases even though the images are unregistered, as shown in Figs. 10b, 12b and 13b. This observation is consistent with the conclusions drawn in [22]. However, we observe from the SAM maps that, it suffers from spectral distortion. The CNMF method handles unregistered image pairs better than the other approaches as shown in Figs. 10d, 11d, 12d and 13d. But its performance is limited by the predefined down-sampling function. The effectiveness of the proposed method can be readily observed from the reconstructed results of difference images shown in Figs. 10g, 11g, 12g and 13g, where the proposed approach has much less spectral and spatial distortion as compared to the state-of-the-art, regardless of the type of input images.

## F. Experimental Results on Unregistered Real Image Pairs

We further evaluate the proposed method on the real unregistered image pairs with both rigid and nonrigid distortions. Since there is no ground truth HR HSI in real applications,

we provide a visual inspection of the reconstructed results in Fig. 14. We can observe that, as long as the LR HSI includes all the spectral bases of HR MSI, the proposed method powered with mutual information is able to increase the spatial resolution of the LR HSI while preserving its spectral resolution well, even when the LR HSI and HR MSI have large pixel displacement.

## G. Ablation and Parameter Study

Taking the challenging rotated 'pompom' image from the CAVE dataset as an example, we further evaluate 1) the necessity of maximizing the mutual information between representations and input images and 2) the usage of collaborative $l_{2,1}$ loss. Since they are all designed to reduce the spectral distortion of the reconstructed image, we use SAM as the evaluation metric.

Fig. 15 illustrates the SAM of the reconstructed HR HSI when increasing the parameters of mutual information $\lambda$ in Eq. 12. We can observe that, if there is no mutual information maximization, i.e. $\lambda = 0$, the spectral information would not be preserved well. When we gradually increase $\lambda$, the reconstructed HR HSI preserves better spectral information, i.e., the SAM is largely reduced. The reason for that is, when we maximize the MI between the representations and their own inputs, it actually maximizes the mutual information of the representations of two modalities. Therefore, the network is able to correlate the extracted spectral and spatial information from unregistered HR MSI and LR MSI in an effective way, to largely reduce the spectral distortion. However, when the parameters are too large, it may hinder the reconstruction procedure of the image pairs. Therefore, we need to choose the proper parameters for the network. In our experiments, we keep $\mu = 1 \times 10^{-4}$ during the experiments to reduce overfitting. We set $\lambda = 1 \times 10^{-5}$ for general HSI dataset with less spectral bands and $\lambda = 1 \times 10^{-1}$ for remote sensing HSI with more spectral bands.

The effectiveness evaluation of the collaborative $l_{2,1}$ norm is demonstrated in Fig. 16. We can observe that with $l_1$ norm, the network converges much slower as compared to those using the $l_2$ norm and $l_{21}$ norm, and the $l_{21}$ norm converges to smaller spectral distortions than using the $l_2$ norm or the $l_1$ norm. Thus, $l_{2,1}$ norm can preserve the spectral information better and significantly reduce the spectral distortion of the restored HR HSI.
<figure><img src="image_12.png" /><figcaption>Fig. 10. Reconstructed results given unregistered rigid distorted image pairs from the CAVE dataset. (a) Color composite of the reference HR HSI. (b) Color composite of the distorted LR HSI. (c) Color composite of the LR HSI. (b)-(g): reconstructed results. (i)-(n): average absolute difference between the reconstructed HSI and reference HSI over different spectral bands, from different methods. (p)-(u) SAM of each pixel between the reconstructed HSI and reference HSI from different methods.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 11. Reconstructed results given unregistered rigid distorted image pairs from the Harvard dataset. (a) Color composite of the reference HR HSI. (b) Color composite of the distorted LR HSI. (c) Color composite of the LR HSI. (b)-(g): reconstructed results. (i)-(n): average absolute difference between the reconstructed HSI and reference HSI over different spectral bands, from different methods. (p)-(u) SAM of each pixel between the reconstructed HSI and reference HSI from different methods.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>Fig. 12. Reconstructed results given unregistered nonrigid distorted image pairs from the Chikusei dataset. (a) Color composite of the reference HR HSI. (h) Color composite of the distorted LR HSI. (o) Color composite of the LR HSI. (b)-(g): reconstructed results. (i)-(n): average absolute difference between the reconstructed HSI and reference HSI over different spectral bands, from different methods. (p)-(u) SAM of each pixel between the reconstructed HSI and reference HSI from different methods.</figcaption></figure>

## H. Tolerance Study

At last, we would like to examine how much spectral information can be preserved when the network deals with unregistered images. To preserve spectral information, the input LR HSI should cover all the spectral signatures of HR MSI. Thus, we choose the image in Fig. 1 from the Harvard dataset which has most of the spectral signatures centered in the image. The results are shown in Fig. 17. The image is rotated from 5 degrees to 30 degrees with 15% to 48% percent of information missing. We can observe that as long as the spectral bases are included in the LR HSI, no matter how small the overlapped region is between the LR HSI and HR MSI, we could always achieve the reconstructed image with small spectral distortion even for unregistered input images.

## VI. CONCLUSION

We proposed an unsupervised encoder-decoder network $u^2$-MDN to solve the problem of hyperspectral image super-resolution without multi-modality registration. The unique structure stabilizes the network training by projecting both modalities into the same space and extracting the spectral basis from LR HSI with rich spectral information as well

as spatial representations from HR MSI with high-resolution spatial information simultaneously. The network learns cor-related spatial information from two unregistered modalities by maximizing the mutual information between the represen-tations and their own raw inputs. In this way, it maximizes the MI between the two representations that largely reduces the spectral distortion. In addition, the collaborative $l_{2,1}$ norm is adopted to encourage the network to further preserve spectral information. Extensive experiments on two benchmark datasets demonstrated the superiority of the proposed approach over the state-of-the-art.

## ACKNOWLEDGMENT

The authors would like to thank all the developers of the evaluated methods who kindly offered their codes, and Dr. Danfeng Hong and Dr. Ke Zhang who provided suggestions on synthetic data generation. This publication was made possible by NASA grant NNX12CB05C and NNX16CP38P.

## REFERENCES

[1] A. Chakrabarti and T. Zickler, "Statistics of real-world hyperspectral images," *The IEEE Conference on Computer Vision and Pattern Recognition (CVPR)*, pp. 193–200, 2011.
<figure><img src="image_15.png" /><figcaption>Fig. 13. Reconstructed results given unregistered nonrigid distorted image pairs from the Pavia dataset. (a) Color composite of the reference HR HSI. (b) Color composite of the distorted LR HSI. (c) Color composite of the LR HSI. (d)-(g): reconstructed results. (i)-(n): average absolute difference between the reconstructed HSI and reference HSI over different spectral bands, from different methods. (p)-(u) SAM of each pixel between the reconstructed HSI and reference HSI from different methods.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Fig. 14. Color composite of (a) the LR HSI of the real data from Hyperion, (b) the HR MSI of the real data from WorldView3 (images courtesy Maxar), and (c) the reconstructed HR HSI from the proposed method.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>Fig. 15. Influence of MI</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Fig. 16. The effect of $l_{2,1}$</figcaption></figure>

<caption>Fig. 17. Tolerance study</caption>

[2] J. M. Bioucas-Dias, A. Plaza, N. Dobigeon, M. Parente, Q. Du, P. Gader, and J. Chanussot, "Hyperspectral unmixing overview: Geometrical, statistical, and sparse regression-based approaches," Selected Topics in Applied Earth Observations and Remote Sensing, IEEE Journal of, vol. 5, no. 2, 2012.

[3] C. Kwan, B. Ayhan, G. Chen, J. Wang, B. Ji, and C.-I. Chang, "A novel approach for spectral unmixing, classification, and concentration estimation of chemical and biological agents," IEEE Transactions on Geoscience and Remote Sensing, vol. 44, no. 2, pp. 409-419, 2006.

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

[5] J. M. Haut, M. E. Paoletti, J. Plaza, J. Li, and A. Plaza, "Active learning with convolutional neural networks for hyperspectral image classification using a new bayesian approach," IEEE Transactions on Geoscience and Remote Sensing, no. 99, pp. 1-22, 2018.

[6] P. S. S. Aydav and S. Minz, "Classification of hyperspectral images using self-training and a pseudo validation set," Remote Sensing Letters, vol. 9, no. 11, pp. 1109-1117, 2018.

[7] B. Uzkent, A. Rangnekar, and M. Hoffman, "Aerial vehicle tracking by adaptive fusion of hyperspectral likelihood maps," The IEEE Conference on Computer Vision and Pattern Recognition Workshops (CVPRW), July 2017.

[8] A. Plaza, Q. Du, J. M. Bioucas-Dias, X. Jia, and F. A. Kruse, "Foreword to the special issue on spectral unmixing of remotely sensed data," IEEE Transactions on Geoscience and Remote Sensing, vol. 49, no. 11, pp. 4103-4110, 2011.

[9] M. Borengasser, W. S. Hungate, and R. Watkins, Hyperspectral remote sensing: principles and applications, 2007.

[10] R. Kawakami, Y. Matsushita, J. Wright, M. Ben-Ezra, Y.-W. Tai, and K. Ikeuchi, "High-resolution hyperspectral imaging via matrix factorization," The IEEE Conference on Computer Vision and Pattern Recognition (CVPR), pp. 2329-2336, 2011.

[11] N. Akhtar, F. Shafait, and A. Mian, "Bayesian sparse representation for hyperspectral image super resolution," The IEEE Conference on Computer Vision and Pattern Recognition (CVPR), pp. 3631-3640, 2015.

[12] C. Lanaras, E. Baltsavias, and K. Schindler, "Hyperspectral super-resolution by coupled spectral unmixing," The IEEE Conference on Computer Vision and Pattern Recognition (CVPR), pp. 3586-3594, 2015.

[13] G. Vivone, L. Alparone, J. Chanussot, M. Dalla Mura, A. Garzelli, G. A. Licciardi, R. Restaino, and L. Wald, "A critical comparison among pansharpening algorithms," IEEE Transactions on Geoscience and Remote Sensing, vol. 53, no. 5, 2015.

[14] C. Thomas, T. Ranchin, L. Wald, and J. Chanussot, "Synthesis of multispectral images to high spatial resolution: A critical review of fusion methods based on remote sensing physics," IEEE Transactions on Geoscience and Remote Sensing, vol. 46, no. 5, pp. 1301-1312, 2008.

[15] B. Aiazzi, S. Baronti, and M. Selva, "Improving component substitution pansharpening through multivariate regression of ms+ pan data," IEEE Transactions on Geoscience and Remote Sensing, vol. 45, no. 10, 2007.

[16] B. Aiazzi, L. Alparone, S. Baronti, A. Garzelli, and M. Selva, "Mtf-tailored multiscale fusion of high-resolution ms and pan imagery," Photogrammetric Engineering & Remote Sensing, vol. 72, no. 5, pp. 591-596, 2006.

[17] 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 Geoscience and Remote Sensing Magazine, vol. 3, no. 3, 2015.

[18] Y. Chang, L. Yan, H. Fang, S. Zhong, and W. Liao, "Hsi-denet: Hyperspectral image restoration via convolutional neural network," IEEE Transactions on Geoscience and Remote Sensing, no. 99, pp. 1-16, 2018.

[19] P. Arun, K. M. Buddhiraju, A. Porwal, and J. Chanussot, "Cnn-based super-resolution of hyperspectral images," IEEE Transactions on Geoscience and Remote Sensing, 2020.

[20] Y. Zhou, A. Rangarajan, and P. D. Gader, "Nonrigid registration of hyperspectral and color images with vastly different spatial and spectral resolutions for spectral unmixing and pansharpening," 2017 IEEE Conference on Computer Vision and Pattern Recognition Workshops (CVPRW), pp. 1571-1579, 2017.

[21] Y. Qu, H. Qi, and C. Kwan, "Unsupervised sparse dirichlet-net for hyperspectral image super-resolution," The IEEE Conference on Computer Vision and Pattern Recognition (CVPR), pp. 2511-2520, 2018.

[22] S. Baronti, B. Aiazzi, M. Selva, A. Garzelli, and L. Alparone, "A theoretical analysis of the effects of aliasing and misregistration on pansharpened imagery," IEEE Journal of Selected Topics in Signal Processing, vol. 5, no. 3, pp. 446-453, 2011.

[23] F. D. Van der Meer, H. M. Van der Werff, F. J. Van Ruitenbeek, C. A. Hecker, W. H. Bakker, M. F. Noomen, M. Van Der Meijde, E. J. M. Carranza, J. B. De Smeth, and T. Woldai, "Multi-and hyperspectral geologic remote sensing: A review," International Journal of Applied Earth Observation and Geoinformation, vol. 14, no. 1, pp. 112-128, 2012.

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

[25] H. Chui and A. Rangarajan, "A new point matching algorithm for non-rigid registration," Computer Vision and Image Understanding, vol. 89, no. 2-3, pp. 114-141, 2003.

[26] X. Fan, H. Rhody, and E. Saber, "A spatial-feature-enhanced mmi algorithm for multimodal airborne image registration," IEEE Transactions on Geoscience and Remote Sensing, vol. 48, no. 6, pp. 2580-2589, 2010.

[27] A. Myronenko and X. Song, "Point set registration: Coherent point drift," IEEE Transactions on pattern analysis and machine intelligence, vol. 32, no. 12, pp. 2262-2275, 2010.

[28] J. Ma, H. Zhou, J. Zhao, Y. Gao, J. Jiang, and J. Tian, "Robust feature matching for remote sensing image registration via locally linear transforming," IEEE Transactions on Geoscience and Remote Sensing, vol. 53, no. 12, pp. 6469-6481, 2015.

[29] M. E. Schaepman, M. Jehle, A. Hueni, P. D'Odorico, A. Damm, J. Weyermann, F. D. Schneider, V. Laurent, C. Popp, F. C. Seidel et al., "Advanced radiometry measurements and earth science applications with the airborne prism experiment (apex)," Remote Sensing of Environment, vol. 158, pp. 207-219, 2015.

[30] M. Selva, B. Aiazzi, F. Butera, L. Chiarantini, and S. Baronti, "Hyper-sharpening: A first approach on sim-ga data," IEEE Journal of Selected Topics in Applied Earth Observations and Remote sensing, vol. 8, no. 6, pp. 3008-3024, 2015.

[31] M. Selva, L. Santurri, and S. Baronti, "Improving hypersharpening for worldview-3 data," IEEE Geoscience and Remote Sensing Letters, vol. 16, no. 6, pp. 987-991, 2019.

[32] S. G. Mallat, "A theory for multiresolution signal decomposition: the wavelet representation," IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 11, no. 7, 1989.

[33] J. Liu, "Smoothing filter-based intensity modulation: A spectral preserve image fusion technique for improving spatial details," International Journal of Remote Sensing, vol. 21, no. 18, pp. 3461-3472, 2000.

[34] P. Burt and E. Adelson, "The laplacian pyramid as a compact image code," IEEE Transactions on Communications, vol. 31, no. 4, pp. 532-540, 1983.
[35] R. Dian, L. Fang, and S. Li, "Hyperspectral image super-resolution via non-local sparse tensor factorization," *The IEEE Conference on Computer Vision and Pattern Recognition (CVPR)*, pp. 5344–5353, 2017.

[36] Q. Wei, N. Dobigeon, and J.-Y. Tourneret, "Bayesian fusion of hyperspectral and multispectral images," *2014 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)*, 2014.

[37] Q. Wei, J. Bioucas-Dias, N. Dobigeon, and J.-Y. Tourneret, "Hyperspectral and multispectral image fusion based on a sparse representation," *IEEE Transactions on Geoscience and Remote Sensing*, vol. 53, no. 7, 2015.

[38] M. Simões, J. Bioucas-Dias, L. B. Almeida, and J. Chanussot, "A convex formulation for hyperspectral image superresolution via subspace-based regularization," *IEEE Transactions on Geoscience and Remote Sensing*, vol. 53, no. 6, pp. 3373–3388, 2015.

[39] N. Yokoya, T. Yairi, and A. Iwasaki, "Coupled nonnegative matrix factorization unmixing for hyperspectral and multispectral data fusion," *IEEE Transactions on Geoscience and Remote Sensing*, vol. 50, no. 2, pp. 528–537, 2012.

[40] M. A. Veganzones, M. Simoes, G. Licciardi, N. Yokoya, J. M. Bioucas-Dias, and J. Chanussot, "Hyperspectral super-resolution of locally low rank images from complementary multisource data," *IEEE Transactions on Image Processing*, vol. 25, no. 1, pp. 274–288, 2016.

[41] E. Wycoff, T.-H. Chan, K. Jia, W.-K. Ma, and Y. Ma, "A non-negative sparse promoting algorithm for high resolution hyperspectral imaging," *Acoustics, Speech and Signal Processing (ICASSP), 2013 IEEE International Conference on*, pp. 1409–1413, 2013.

[42] S. Li, R. Dian, L. Fang, and J. M. Bioucas-Dias, "Fusing hyperspectral and multispectral images via coupled sparse tensor factorization," *IEEE Transactions on Image Processing*, vol. 27, no. 8, pp. 4118–4130, 2018.

[43] Y. Chang, L. Yan, X.-L. Zhao, H. Fang, Z. Zhang, and S. Zhong, "Weighted low-rank tensor recovery for hyperspectral image restoration," *IEEE Transactions on Cybernetics*, 2020.

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

[45] Y. Xu, Z. Wu, J. Chanussot, and Z. Wei, "Hyperspectral images super-resolution via learning high-order coupled tensor ring representation," *IEEE Transactions on Neural Networks and Learning Systems*, 2020.

[46] Chen, Yeqing, Li, Wei, Liu, Junzhou, and Huang, "Sirf: Simultaneous satellite image registration and fusion in a unified framework." *IEEE Transactions on Image Processing A Publication of the IEEE Signal Processing Society*, 2015.

[47] Y. Zhou, A. Rangarajan, and P. D. Gader, "An integrated approach to registration and fusion of hyperspectral and multispectral images," *IEEE Transactions on Geoscience and Remote Sensing*, vol. 58, no. 5, pp. 3020–3033, 2020.

[48] C. Dong, C. C. Loy, K. He, and X. Tang, "Image super-resolution using deep convolutional networks," *IEEE Transactions on Pattern Analysis and Machine Intelligence*, vol. 38, no. 2, pp. 295–307, 2016.

[49] C. Ledig, L. Theis, F. Huszár, J. Caballero, A. Cunningham, A. Acosta, A. Aitken, A. Tejani, J. Totz, Z. Wang et al., "Photo-realistic single image super-resolution using a generative adversarial network," *arXiv preprint arXiv:1609.04802*, 2016.

[50] W. Huang, L. Xiao, Z. Wei, H. Liu, and S. Tang, "A new pan-sharpening method with deep neural networks," *IEEE Geoscience and Remote Sensing Letters*, vol. 12, no. 5, pp. 1037–1041, 2015.

[51] G. Masi, D. Cozzolino, L. Verdoliva, and G. Scarpa, "Pansharpening by convolutional neural networks," *Remote Sensing*, vol. 8, no. 7, p. 594, 2016.

[52] Y. Wei, Q. Yuan, H. Shen, and L. Zhang, "Boosting the accuracy of multi-spectral image pan-sharpening by learning a deep residual network," *arXiv preprint arXiv:1705.07556*, 2017.

[53] Y. Li, J. Hu, X. Zhao, W. Xie, and J. Li, "Hyperspectral image super-resolution using deep convolutional neural network," *Neurocomputing*, vol. 266, pp. 29–41, 2017.

[54] R. Dian, S. Li, A. Guo, and L. Fang, "Deep hyperspectral image sharpening," *IEEE Transactions on Neural Networks and Learning Systems*, 2018.

[55] Q. Xie, M. Zhou, Q. Zhao, D. Meng, W. Zuo, and Z. Xu, "Multispectral and hyperspectral image fusion by ms/hss fusion net," *2019 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR)*, 2019.

[56] C. Kwan, B. Budavari, A. C. Bovik, and G. Marchisio, "Blind quality assessment of fused worldview-3 images by using the combinations of pansharpening and hypersharpening paradigms," *IEEE Geoscience and Remote Sensing Letters*, vol. 14, no. 10, pp. 1835–1839, 2017.

[57] C. Kwan, C. Haberle, A. Echavarren, B. Ayhan, B. Chou, B. Budavari, and S. Dickenshied, "Mars surface mineral abundance estimation using themis and tes images," *IEEE Ubiquitous Computing, Electronics and Mobile Communication Conference, New York City*, November 2018.

[58] Y. Fu, T. Zhang, Y. Zheng, D. Zhang, and H. Huang, "Hyperspectral image super-resolution with optimized rgb guidance," *2019 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR)*, pp. 11661–11670, 2019.

[59] K. Zheng, L. Gao, W. Liao, D. Hong, B. Zhang, X. Cui, and J. Chanussot, "Coupled convolutional neural network with adaptive response function learning for unsupervised hyperspectral super resolution," *IEEE Transactions on Geoscience and Remote Sensing*, 2020.

[60] J. Sethuraman, "A constructive definition of dirichlet priors," *Statistica Sinica*, pp. 639–650, 1994.

[61] E. Nalisnick and P. Smyth, "Deep generative models with stick-breaking priors," *International Conference on Machine Learning (ICML)*, 2017.

[62] P. Kumaraswamy, "A generalized probability density function for double-bounded random processes," *Journal of Hydrology*, vol. 46, no. 1-2, pp. 79–88, 1980.

[63] C. Dugas, Y. Bengio, F. Bélisle, C. Nadeau, and R. Garcia, "Incorporating second-order functional knowledge for better option pricing," *Advances in neural information processing systems*, pp. 472–478, 2001.

[64] J. Han and C. Moraga, "The influence of the sigmoid function parameters on the speed of backpropagation learning," *From Natural to Artificial Neural Computation*, pp. 195–201, 1995.

[65] G. Huang, Z. Liu, K. Q. Weinberger, and L. van der Maaten, "Densely connected convolutional networks," *arXiv preprint arXiv:1608.06993*, 2016.

[66] J. B. Kinney and G. S. Atwal, "Equitability, mutual information, and the maximal information coefficient," *Proceedings of the National Academy of Sciences*, p. 201309933, 2014.

[67] B. Zitova and J. Flusser, "Image registration methods: a survey," *Image and vision computing*, vol. 21, no. 11, pp. 977–1000, 2003.

[68] J. Woo, M. Stone, and J. L. Prince, "Multimodal registration via mutual information incorporating geometric and spatial context," *IEEE Transactions on Image Processing*, vol. 24, no. 2, pp. 757–769, 2015.

[69] I. Belghazi, S. Rajeswar, A. Baratin, R. D. Hjelm, and A. Courville, "Mine: mutual information neural estimation," *arXiv preprint arXiv:1801.04062*, 2018.

[70] M. D. Donsker and S. S. Varadhan, "Asymptotic evaluation of certain markov process expectations for large time. iv," *Communications on Pure and Applied Mathematics*, vol. 36, no. 2, pp. 183–212, 1983.

[71] R. D. Hjelm, A. Fedorov, S. Lavoie-Marchildon, K. Grewal, A. Trischler, and Y. Bengio, "Learning deep representations by mutual information estimation and maximization," *arXiv preprint arXiv:1808.06670*, 2018.

[72] F. Nie, H. Huang, X. Cai, and C. H. Ding, "Efficient and robust feature selection via joint $L_2$, 1-norms minimization," *Advances in neural information processing systems*, pp. 1813–1821, 2010.

[73] F. Yasuma, T. Mitsunaga, D. Iso, and S. K. Nayar, "Generalized assorted pixel camera: postcapture control of resolution, dynamic range, and spectrum," *IEEE Transactions on Image Processing*, vol. 19, no. 9, pp. 2241–2253, 2010.

[74] N. Yokoya and A. Iwasaki, "Airborne hyperspectral data over chikusei," *Space Appl. Lab., Univ. Tokyo, Tokyo, Japan, Tech. Rep. SAL-2016-05-27*, 2016.

[75] C. Debes, A. Merentitis, R. Heremans, J. Hahn, N. Frangiadakis, T. van Kasteren, W. Liao, R. Bellens, A. Piżurica, S. Gautama et al., "Hyperspectral and lidar data fusion: Outcome of the 2013 grss data fusion contest," *IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing*, vol. 7, no. 6, pp. 2405–2418, 2014.

[76] W. Dong, F. Fu, G. Shi, X. Cao, J. Wu, G. Li, and X. Li, "Hyperspectral image super-resolution via non-negative structured sparse representation," *IEEE Transactions on Image Processing*, vol. 25, no. 5, pp. 2337–2352, 2016.

[77] L. Wald, T. Ranchin, and M. Mangolini, "Fusion of satellite images of different spatial resolutions: Assessing the quality of resulting images," *Photogrammetric Engineering & Remote Sensing*, vol. 63, no. 6, pp. 691–699, 1997.

[78] R. C. Hardie, M. A. Rucci, A. J. Dapore, and B. K. Karch, "Block matching and wiener filtering approach to optical turbulence mitigation and its application to simulated and real imagery with quantitative error analysis," *Optical Engineering*, vol. 56, no. 7, p. 071503, 2017.

[79] F. A. Kruse, A. Lefkoff, J. Boardman, K. Heidebrecht, A. Shapiro, P. Barloon, and A. Goetz, "The spectral image processing system (sips)—interactive visualization and analysis of imaging spectrometer data," *Remote sensing of environment*, vol. 44, no. 2-3, pp. 145–163, 1993.
<figure><img src="image_19.png" /></figure>

**Ying Qu** (S'16–M'18) received the B.S. degree in automatics and M.S. degree in pattern recognition & artificial intelligence from Northeastern University, Shenyang, China in 2008 and 2010, respectively, and the Ph.D. degree in computer engineering from the University of Tennessee, Knoxville, in 2017. She is currently a research associate with the Department of Electrical Engineering and Computer Science at the University of Tennessee, Knoxville. Her current research interests are remote sensing, artificial intelligence and computer vision. She was the recipient

of the Best Student Paper Awards at The International Geoscience and Remote Sensing Symposium (IGARSS) in 2016.

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

**Hairong Qi** (IEEE Fellow since 2017) received the B.S. and M.S. degrees in computer science from Northern JiaoTong University, Beijing, China in 1992 and 1995, respectively, and the Ph.D. degree in computer engineering from North Carolina State University, Raleigh, in 1999. She is currently the Gonzalez Family Professor with the Department of Electrical Engineering and Computer Science at the University of Tennessee, Knoxville. Her cur-

urrent research interests are in advanced imaging and collaborative processing in resource-constrained distributed environment, hyperspectral image analysis, and automatic target recognition. Dr. Qi's research is supported by National Science Foundation (NSF), DARPA, Office of Naval Research (ONR), Department of Homeland Security (DHS), U.S. Army Space and Missile Defense Command, and U.S. Army Medical Research and Materiel Command. Dr. Qi is the recipient of the NSF CAREER Award. She also received the Best Paper Awards at the 18th International Conference on Pattern Recognition (ICPR) in 2006, the 3rd ACM/IEEE International Conference on Distributed Smart Cameras (ICDSC) in 2009, and IEEE Workshop on Hyperspectral Image and Signal Processing: Evolution in Remote Sensor (WHISPERS) in 2015. She is awarded the Highest Impact Paper from the IEEE Geoscience and Remote Sensing Society in 2012.

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

**Chiman Kwan** (S'85–M'93–SM'98) received his BS with honors in Electronics from the Chinese University of Hong Kong in 1988, and MS and Ph.D. degrees in electrical engineering from the University of Texas at Arlington in 1989 and 1993, respectively. Currently, he is the Chief Technology Officer of Signal Processing, Inc. and Applied Research LLC, leading research and development efforts in chemical agent detection, biometrics, speech processing, im- age fusion, mission planning, and fault diagnostics and prognostics. His primary research areas include

robust and adaptive control methods, signal and image processing, commu- nications, neural networks, and pattern recognition applications. From April 1991 to February 1994, he worked in the Beam Instrumentation Department of the SSC (Superconducting Super Collider Laboratory) in Dallas, Texas, where he was heavily involved in the modeling, simulation and design of modern digital controllers and signal processing algorithms for the beam control and synchronization system. He later joined the Automation and Robotics Research Institute in Fort Worth, where he applied intelligent control methods such as neural networks and fuzzy logic to the control of power systems, robots, and motors. Between July 1995 and April 2006, he was with Intelligent Automation, Inc. in Rockville, Maryland. He has served as Principal Investigator/Program Manager for more than 120 different projects. Dr. Kwan has 15 patents, 65 invention disclosures, more than 120 papers in archival journals and more than 250 additional papers published in major conference proceedings. He is listed in the New Millennium edition of Who's Who in Science and Engineering and is a member of Tau Beta Pi. He also received several awards from IEEE related to fault diagnostics and prognostics, and a certificate of recognition from NASA for the health monitoring of Auxiliary Power Units in the Space Shuttle.

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

**Naoto Yokoya** (S'10–M'13) received the M.Eng. and Ph.D. degrees from the Department of Aero- nautics and Astronautics, the University of Tokyo, Tokyo, Japan, in 2010 and 2013, respectively. He is currently a Lecturer at the University of Tokyo and a Unit Leader at the RIKEN Center for Ad- vanced Intelligence Project, Tokyo, Japan, where he leads the Geoinformatics Unit. He was an Assistant Professor at the University of Tokyo from 2013 to 2017. In 2015-2017, he was an Alexander von Humboldt Fellow, working at the German Aerospace

Center (DLR), Oberpfaffenhofen, and Technical University of Munich (TUM), Munich, Germany. His research is focused on the development of image processing, data fusion, and machine learning algorithms for understanding remote sensing images, with applications to disaster management. Dr. Yokoya won the first place in the 2017 IEEE Geoscience and Remote Sensing Society (GRSS) Data Fusion Contest organized by the Image Analysis and Data Fusion Technical Committee (IADF TC). He is the Chair (2019-2021) and was a Co-Chair (2017-2019) of IEEE GRSS IADF TC and also the secretary of the IEEE GRSS All Japan Joint Chapter since 2018. He is an Associate Editor for the IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing (JSTARS) since 2018. He is/was a Guest Editor for the IEEE JSTARS in 2015-2021, for Remote Sensing in 2016-2021, and for the IEEE Geoscience and Remote Sensing Letters (GRSL) in 2018-2019.

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

**Jocelyn Chanussot** (M'04–SM'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 GRSS 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, since 2018).