# Interpretable Hyperspectral AI: When Non-Convex Modeling meets Hyperspectral
Remote Sensing

Danfeng Hong, Wei He, Naoto Yokoya, Jing Yao, Lianru Gao, Liangpei Zhang,
Jocelyn Chanussot, and Xiao Xiang Zhu D. Hong is with the Remote Sensing
Technology Institute (IMF), German Aerospace Center (DLR), 82234 Wessling,
Germany, and also with the Univ. Grenoble Alpes, CNRS, Grenoble INP, GIPSA-
lab, 38000 Grenoble, France. (e-mail: danfeng.hong@dlr.de)W. He is with the
Geoinformatics Unit, RIKEN Center for Advanced Intelligence Project (AIP),
RIKEN, 103-0027 Tokyo, Japan. (e-mail: wei.he@riken.jp)N. Yokoya is with
Graduate School of Frontier Sciences, the University of Tokyo, 277-8561 Chiba,
Japan, and also with the Geoinformatics Unit, RIKEN Center for Advanced
Intelligence Project (AIP), RIKEN, 103-0027 Tokyo, Japan. (e-mail:
naoto.yokoya@riken.jp)L. Gao and J. Yao are with the Key Laboratory of Digital
Earth Science, Aerospace Information Research Institute, Chinese Academy of
Sciences, Beijing 100094, China. (e-mail: gaolr@aircas.ac.cn)L. Zhang is with
the State Key Laboratory of Information Engineering in Surveying, Mapping and
Remote Sensing, Wuhan University, Wuhan 430072, China.
(e-mail:zlp62@whu.edu.cn)J. Chanussot is with the Univ. Grenoble Alpes, INRIA,
CNRS, Grenoble INP, LJK, 38000 Grenoble, France, also with the Aerospace
Information Research Institute, Chinese Academy of Sciences, 100094 Beijing,
China. (e-mail: jocelyn@hi.is)X. Zhu is with the Remote Sensing Technology
Institute (IMF), German Aerospace Center (DLR), 82234 Wessling, Germany, and
Data Science in Earth Observation (SiPEO), Technical University of Munich
(TUM), 80333 Munich, Germany. (e-mail<EMAIL_ADDRESS>
###### Abstract

This is the pre-acceptance version, to read the final version please go to
IEEE Geoscience and Remote Sensing Magazine on IEEE Xplore. Hyperspectral
imaging, also known as image spectrometry, is a landmark technique in
geoscience and remote sensing (RS). In the past decade, enormous efforts have
been made to process and analyze these hyperspectral (HS) products mainly by
means of seasoned experts. However, with the ever-growing volume of data, the
bulk of costs in manpower and material resources poses new challenges on
reducing the burden of manual labor and improving efficiency. For this reason,
it is, therefore, urgent to develop more intelligent and automatic approaches
for various HS RS applications. Machine learning (ML) tools with convex
optimization have successfully undertaken the tasks of numerous artificial
intelligence (AI)-related applications. However, their ability in handling
complex practical problems remains limited, particularly for HS data, due to
the effects of various spectral variabilities in the process of HS imaging and
the complexity and redundancy of higher dimensional HS signals. Compared to
the convex models, non-convex modeling, which is capable of characterizing
more complex real scenes and providing the model interpretability technically
and theoretically, has been proven to be a feasible solution to reduce the gap
between challenging HS vision tasks and currently advanced intelligent data
processing models.

This article mainly presents an advanced and cutting-edge technical survey for
non-convex modeling towards interpretable AI models covering a board scope in
the following topics of HS RS:

  * •

HS image restoration,

  * •

dimensionality reduction,

  * •

data fusion and enhancement,

  * •

spectral unmixing,

  * •

cross-modality learning for large-scale land cover mapping.

Around these topics, we will showcase the significance of non-convex
techniques to bridge the gap between HS RS and interpretable AI models with a
brief introduction on the research background and motivation, an emphasis on
the resulting methodological foundations and solution, and an intuitive
clarification of illustrative examples. At the end of each topic, we also pose
the remaining challenges on how to completely model the issues of complex
spectral vision from the perspective of intelligent ML combined with physical
priors and numerical non-convex modeling, and accordingly point out future
research directions.

This paper aims to create a good entry point to the advanced literature for
experienced researchers, Ph.D. students, and engineers who already have some
background knowledge in HS RS, ML, and optimization. This can further help
them launch new investigations on the basis of the above topics and
interpretable AI techniques for their focused fields.

###### Index Terms:

Artificial intelligence, image processing, interpretability, hyperspectral,
machine learning, modeling, non-convex, remote sensing, signal processing.

##  I Introduction

###  I-A Background and Significance of Hyperspectral Remote Sensing

Imaging spectroscopy, which was first-ever to be conceptualized by Goetz1 et
al. in 1980’s [1], is a seminal hyperspectral (HS) imaging technique of truly
achieving the integration of the 1-D spectrum and the 2-D image for earth
remote sensing (RS). Imaging spectroscopy is a typical “passive” RS technique,
which assembles spectroscopy and digital photography into a unified system.
Fig. 1 shows the data acquisition process of two different imaging patterns:
“active” RS and “passive” RS [2]. The resulting HS image collects hundreds of
2-D images finely sampled from the approximately contiguous wavelength across
the whole electromagnetic spectrum [3] (see Fig. 2). This enables the
recognition and identification of the materials, particularly for those that
have extremely similar spectral signatures in visual cues (e.g., RGB) [4], at
a more accurate and finer level. As a result, HS RS has been significantly
advanced and widely applied in many challenging tasks of earth observation
[5], such as fine-grained land cover classification, mineral mapping, water
quality assessment, precious farming, urban planning and monitoring, disaster
management and prediction, and concealed target detection.

Figure 1: An illustration to clarify the data acquisition process of two
different imaging patterns, i.e., “active” RS and “passive” RS.  Figure 2: A
showcase to clarify the electromagnetic spectrum: the order from low to high
frequency is Long-waves, Radio-waves, Micro-waves, Infrared, Visible,
Ultraviolet, X-rays, and Gamma-rays, where four widely-concerned intervals,
e.g., Radio-waves, Infrared, Visible, and X-rays, are finely partitioned.

More specifically, characterized by the distinctive 3-D signaling structure,
the advantages of the HS image over conventional 1-D or 2-D signal products
can be summarized as

  * •

compared to the common 1-D signal system, the HS 2-D spatial pattern provides
us more structured information, enabling the discrimination of underlying
objects of interest at a more semantically meaningful level;

  * •

Beyond the 2-D natural images, the rich spectral information in HS images is
capable of detecting materials through tiny discrepancies in the spectral
domain, since the HS imaging instruments exploit the sensors to collect
hundreds or thousands of wavelength channels with an approximately continuous
spectral sampling at a subtle interval (e.g., 10nm).

Furthermore, the significance of HS images compared to other RS imaging
techniques with a lower spectral resolution, e.g., multispectral (MS) or RGB
imaging, mainly embodies in the following three aspects:

  * 1)

HS images are capable of finely discriminating the different classes that
belong to the same category, such as Bitumen and Asphalt, Stressed Grass and
Synthetic Grass, Alunite and Kaolin. For those optical broadband imaging
products (e.g., MS imagery), they can only identify certain materials with the
observable differences in the spectral signatures, e.g., Water, Trees, Soil.

  * 2)

The higher spectral resolution creates the possibilities to some challenging
applications that can be hardly achieved by only depending on formerly imaging
techniques, e.g., parameter extraction of biophysics and biochemistry,
biodiversity conservation, monitoring and management of the ecosystem,
automatic detection of food safety, which provides new insight into RS and
geoscience fields.

  * 3)

Due to the limitations of image resolution either in spectral or spatial
domains, physically and chemically atmospheric effects, and environmental
conditions (e.g., the interference of soil background, illumination, the
uncontrolled shadow caused by clouds or building occlusion, topography change,
complex noises), those traditional RS imaging techniques were to a great
extent dominated by qualitative analysis. As HS RS arises, quantitative or
semi-quantitative analysis becomes increasingly plausible in many practical
cases.

###  I-B An Ever-Growing Relation between Non-convex Modeling and
Interpretable AI in Hyperspectral Remote Sensing

In recent years, a vast number of HS RS missions (e.g., MODIS, HypSEO, DESIS,
Gaofen-5, EnMap, HyspIRI, etc.) have been launched to enhance our
understanding and capabilities to the Earth and environment, contributing to
the rapid and better development in a wide range of relevant applications,
such as land cover land use classification, spectral unmixing, data fusion,
image restoration, and multimodal data analysis. With the ever-growing
availability of RS data sources from both satellite and airborne sensors on a
large scale and even global scale, expert system-centric data processing and
analysis mode has run into bottlenecks and can not meet the demand of the big
data era. For this reason, data-driven signal and image processing, machine
learning (ML), AI models have been garnering growing interest and attention
from the researchers in the RS community.

Supported by well-established theory and numerical optimization, convex models
have been proven to be effective to model a variety of HS tasks under highly
idealized assumptions. However, there exist unknown, uncertain, and
unpredictable factors in the complex real scenes. Due to these factors that
lead to the lack of sound understanding and modeling capability to the scene,
convex models fail to work properly. The specific reasons could be two-fold.
On the one hand, integrating the benefits of 1-D and 2-D signals, the 3-D
structurized HS images offer greater potential and better solutions (compared
to natural images) to deal with the varying situation, but simultaneously
increase the model’s complexity and uncertainty to some extent. On the other
hand, due to unprecedented spatial, spectral, and temporal resolutions of HS
images in remotely sensed HS imaging, the difficulties and challenges in the
sophisticated HS vision approaches are mainly associated with the volume of
the HS data, complex material (spectral) mixing behavior, and uncontrolled
degradation mechanisms in data acquisition caused by illumination, noise, and
atmospheric effects.

The aforementioned factors, to a great extent, limit convex models to be
intelligent approaches for fully understanding and interpreting the real-life
scenario. Therefore, this naturally motivates us to investigate the
possibility of processing and analyzing the HS data in a non-convex modeling
fashion. In the following, we briefly make a qualitative comparison between
convex and non-convex models to clarify that non-convex modeling might be an
optimally feasible solution towards interpretable AI models in HS RS.

  * •

Convex models are theoretically guaranteed to converge to the global optimal
solution, yet most tasks related to HS RS are complex in reality and hardly
simplified to an equivalent and perfect convex formulation. This to some
extent makes convex models inapplicable to practical tasks, due to the lack of
interpretability and completeness for problem modeling.

  * •

Rather, non-convex models are capable of characterizing the complex studied
scene in HS RS more finely and completely, thereby more tending to achieve
automatization and intelligentization in the real world. Moreover, by
excavating the intrinsic properties from the HS data effectively to yield
physically meaningful priors, the solution space of non-convex models can be
shrunk to a “good” region bit by bit.

  * •

Although non-convex models are complex by considering more complicated prior
knowledge, possibly leading to the lack of stable generalization ability, they
hold the higher potential that convex models do not have, particularly in
explaining models, understanding scenes, and achieving intelligent HS image
processing and analysis. Furthermore, this might be able to provide
researchers with a broader range of HS vision related topics, making it
applicable for more real cases in a variety of HS RS-related tasks.

###  I-C Contributions

With the advent of the big data era, an ever-increasing data bulk and
diversity brings rare opportunities and challenges for the development of HS
RS in earth observation. Data-driven AI approaches, e.g., ML-based, deep
learning (DL)-based, have occupied a prominent place in manifold HS RS
applications. Nevertheless, how to open the “model” and give them
interpretability remains unknown yet. In this article, we raise a bold and
understandable standpoint, that is, non-convex modeling might be an effective
means to bridge the gap between interpretable AI models and HS RS. To support
the opinion, this article provides a detailed and systematic overview by
reviewing advanced and latest literature with an emphasis on non-convex
modeling in terms of five classic and burgeoning topics related to HS RS. More
specifically,

  * •

We present a comprehensive discussion and analysis related to non-convex
modeling in five well-noticed and promising HS RS-related applications, such
as HS image restoration, dimensionality reduction and classification, spectral
unmixing, data fusion and enhancement, and cross-modality learning.

  * •

For each topic, a few representative works are emphatically and detailedly
introduced by the attempts to make a connection between non-convex modeling
and intelligent/interpretable models. Moreover, the example experiments
(qualitatively or quantitatively) are subsequently performed after the
detailed method description. Those selected methods that engage in the
comparative experiments are accompanied by available code and data links for
the sake of reproducibility. Finally, the remaining challenges are highlighted
to further clarify the gap between interpretable ML/AI models and practical HS
RS applications.

  * •

Regarding the three aspects of non-convex modeling, interpretable AI, and HS
RS, the end of this article concludes with some remarks, makes summary
analysis, and hints at plausible future research work.

We need to point out and emphasize, however, that this paper features food for
thoughts for advanced readers, and it is not an introduction for beginners
entering the field. The goal of this paper is to provide a cutting-edge survey
rather than a real tutorial. As a result, readers are expected to have some
prior knowledge across multidisciplinary, such as HS RS, convex optimization,
non-convex modeling, ML, and AI, where some basic principle, definitions, and
deductions need to be mastered. For beginners who are willing to start with
new researches on non-convex modeling for HS RS applications, we recommend
reading and learning following materials and references to get to know, for
example,

  * •

what is the HS imaging or HS RS (e.g., principles, superiority,
practicability) and its relevant applications (topics) [5];

  * •

convex optimization and its solutions (including a detailed description of
physical meaningful priors, e.g., low-rank, sparsity, graph regularization,
non-negativity, sum-to-one, etc.) as well as its relationship with non-convex
modeling [6];

  * •

a general guideline on why and how to build non-convex models and how to solve
non-convex minimization problems [7];

  * •

classic and advanced ML algorithms, including essential ideas, designing
thought, implementation process [8];

  * •

a big picture about what is the AI and how to build the basic AI models [9].

Moreover, we hope that this paper can be also regarded as a good starting
point and evolve many novels, interesting, and noteworthy research issues
around the fields of non-convex modeling, interpretable ML/AI, and HS RS,
serving to more application cases in reality.

Figure 3: Illustration of five promising topics in HS RS, including image
restoration, dimensionality reduction and classification, data fusion and
enhancement, spectral unmixing, and cross-modality learning.

##  II Outline and Basic Preparation

This paper starts with a brief introduction to a general non-convex
minimization problem in signal or image processing, and then specifies the
non-convex modeling with each topic in HS RS from an ML perspective. According
to the characteristics of HS imaging, these different issues may bring fresh
challenges to researches on non-convex modeling and optimization, which
contributes to boosting the development of both HS RS and ML for intelligent
data processing and analysis. Very recently, some successful showcases have
garnered growing attention from researchers who engage in HS data processing
and analysis, ML, or statistical optimization, leading to many newly-developed
non-convex models and their solutions. They can be roughly categorized into
the following several groups in a wide range of real applications related to
HS RS. Fig. 3 gives the illustration for each topic.

###  II-A A Brief Introduction from Convex to Non-convex Models

As the name suggests, in the convex model, the shape of the area represented
by the function (or model) is convex. Accordingly, let the function
$f:\mathbb{R}_{n}\rightarrow\mathbb{R}$ be convex, if the domain of $f$ is a
convex set, and for the variable $\theta\in[0,1]$, any two points (e.g., $x$
and $y$) meet the following condition:

$\displaystyle f(y)\geq\theta f(x)+(1-\theta)f(y),$ (1)

then the function $f$ is convex, whose necessary and sufficient condition is
$f(y)\geq f(x)+\bigtriangledown f(x)(x-y)$.

Within the domain, the globally optimal solution of the convex model can be
obtained by common convex optimization techniques, such as linear programming
[10], quadratic programming [11], second order cone programming [12].

Informally, although convex methods have been widely used to model various
tasks, owing to the existence and uniqueness of solutions and the relatively
low model complexity, the real scenes are complex and changeable, inevitably
leading to many uncertainties and difficulties in the modeling process. For
this reason, non-convex modeling is capable of providing stronger modeling
power to the algorithm designer and can fully meet the demand of
characterizing the real complex scenes well. This naturally motivates us to
shift our emphases on some key issues related to non-convex modeling.

A general non-convex model usually consists of a smooth objective function
(e.g., Euclidean loss, negative log-likelihood) with Lipschitz gradient
$f(\mathbf{X})$ and the non-convex constraints $\mathcal{C}$, which can be
generalized by optimizing the following minimization problem:

$\displaystyle\mathop{\min}_{\mathbf{X}}\frac{1}{2}f(\mathbf{X})\;\;{\rm
s.t.}\;\mathbf{X}\in\mathcal{C},$ (2)

where $\mathbf{X}$ is the to-be-estimated variable, which can be defined as a
vector (1-D signal), a matrix (2-D image), or an unfolded matrix (3-D HS
image). The constraints $C$ could be sparsity-promoting variants, low-rank,
TV, and others, which need to be determined by the specific tasks.

Unlike convex models, there exist many local minimums in non-convex models.
This poses a big challenge on finding globally optimal solutions. For possible
solutions of non-convex methods, one strategy is to relax the non-convex
problem to an approximately convex model [13]. Another can break the non-
convex problems down into several convex subproblems and solve them in
parallel by the means of convex ways [14].

In the light of the different research goals, the general model in Eq. (2) can
be extended to the task-driven variants covering a broad scope within the HS
RS, including image restoration, dimensionality reduction, data fusion and
enhancement, spectral unmixing, and cross-modality learning. It should be
noted that in the following sections, some representative methods will be
introduced with a focus on non-convex modeling, while the alternating
direction method of multipliers (ADMM) optimization framework is recommended
for use as a general solver to solve these non-convex models. For specific
solutions of these models in each topic, please refer to the cited references.

###  II-B Main Abbreviation

AI: | artificial intelligence.
---|---
ANN: | artificial neural network.
CML: | cross-modality learning.
DL: | deep learning.
DR: | dimensionality reduction.
GSD: | ground sampling distance.
HS: | hyperspectral.
KNN: |  $k$ nearest neighbors.
LDA: | linear discriminant analysis.
LMM: | linear mixing model.
MA: | manifold alignment.
ML: | machine learning.
MML: | multimodality learning.
MS: | multispectral.
NMF: | non-negative matrix factorization.
RS: | remote sensing.
SAR: | synthetic aperture radar.
SOTA: | state-of-the-art.
SSL: | shared subspace learning.
SU: | spectral unmixing.
1-D: | one-dimensional.
2-D: | two-dimensional.
3-D: | three-dimensional.

###  II-C Nomenclature

$\boldsymbol{\mathcal{X}}$: | to-be-estimated 3-D HS image.
---|---
$\boldsymbol{\mathcal{Y}}$: | observed 3-D HS image.
$\boldsymbol{\mathcal{N}}_{G}$: | 3-D Gaussian noise.
$\boldsymbol{\mathcal{N}}_{S}$: | 3-D sparse noise.
$\boldsymbol{\mathcal{O}}$: | core tensor.
$r$: | rank of matrix.
$\mathbf{X}$: | unfolded 2-D matrix of $\boldsymbol{\mathcal{X}}$.
$\mathbf{x}_{i}$: | the $i$-th pixel (1-D vector) of $\mathbf{X}$.
$\mathbf{Y}$: | unfolded 2-D matrix of $\boldsymbol{\mathcal{Y}}$.
$\mathbf{N}_{S}$: | unfolded 2-D matrix of $\boldsymbol{\mathcal{N}}_{S}$.
$\mathbf{H}$: | first-order difference matrix.
$\mathcal{C}$: | model constraint set.
$\Phi$: | to-be-estimated variable set.
$f$: | transformation functions.
$\mathbf{Q}$: | combination coefficients of NMF.
$\mathbf{W}$: | graph or manifold structure.
$\mathbf{L}$: | Laplacian matrix.
$\mathbf{D}$: | degree matrix.
$\mathbf{U}$: | subspace projection matrix.
$\mathbf{I}$: | identity matrix.
$\mathbf{d}$: | distance or similarity matrix.
$\sigma$: | standard derivation.
$C_{k}$: | sample set of the $k$-th class.
$\mathbf{M}$: | one-hot encoded matrix.
$\mathbf{P}$: | regression matrix.
$\mathbf{E}$: | endmember matrix.
$\mathbf{E}_{0}$: | reference endmember matrix.
$\mathbf{A}$: | abundance matrix.
$\mathbf{S}$: | scaling factors (matrix).
$\mathbf{V}$: | spectral variability dictionary (matrix).
$\mathbf{J}$: | coefficients corresponding to $\mathbf{V}$.
$\mathbf{R}$: | spatial degradation function.
$\mathbf{G}$: | spectral response function.
$\mathbf{N}_{H}$: | HS noise.
$\mathbf{N}_{H}$: | MS noise.
$\mathbf{Z}$: | high spatial resolution MS image.
$m$: | the number of the considered modality.
$c$: | scaling constant.

###  II-D Notation

$\lVert\mathbf{X}\rVert_{\operatorname{F}}$ | Forbenius norm of $\mathbf{X}$, obtained by $\sqrt{\sum_{i,j}\mathbf{X}_{i,j}^{2}}$
---|---
$\lVert\mathbf{X}\rVert_{1,1}$ |  $\ell_{1}$ norm of $\mathbf{X}$, obtained by $\sum_{i,j}|\mathbf{X}_{i,j}|$
$\lVert\mathbf{X}\rVert_{2,1}$ |  $\ell_{2,1}$ norm of $\mathbf{X}$, obtained by $\sum_{i}|\sqrt{\sum_{j}\mathbf{X}_{i,j}^{2}}|$
$\lVert\mathbf{X}\rVert_{1/2}$ |  $\ell_{1/2}$ norm of $\mathbf{X}$, obtained by $\sum_{i,j}\mathbf{x}_{j}(i)^{1/2}$
$\lVert\mathbf{X}\rVert_{q}$ |  $\ell_{q}$ norm of $\mathbf{X}$, obtained by $\sum_{i,j}\mathbf{x}_{j}(i)^{q}$
$\lVert\mathbf{X}\rVert_{{\rm TV}}$ |  ${\rm TV}$ norm of $\mathbf{X}$, obtained by $\lVert\mathbf{H}_{h}\mathbf{X}+\mathbf{H}_{v}\mathbf{X}\rVert_{2,1}$
$\lVert\mathbf{X}\rVert_{0}$ |  $\ell_{0}$ norm of $\mathbf{X}$, obtained by $\lim_{p\rightarrow 0}\sum_{i,j}|\mathbf{X}_{i,j}|^{p}$
$\operatorname{tr}(\mathbf{X})$ | trace of $\mathbf{X}$, obtained by $\sum_{i}\mathbf{X}_{i,i}$
$\lVert\mathbf{X}\rVert_{*}$ | nuclear norm of $\mathbf{X}$, obtained by $\operatorname{tr}(\sqrt{\mathbf{X}^{\top}\mathbf{X}})$
$\lVert\mathbf{x}\rVert_{2}$ |  $\ell_{2}$ norm of $\mathbf{x}$, obtained by $\sqrt{\sum_{j}\mathbf{x}_{j}^{2}}$
$\odot$ | the element-wise multiplication operator
$\phi_{i}$ | the neighbouring pixels of the target pixel $i$

##  III Hyperspectral Image Restoration

Owing to the wealthy spectral information, the HS image has been widely used
in different kinds of applications, including urban planning, agriculture,
forestry, target detection, and so on. However, due to the limitation of
hyperspectral imaging systems and the weather conditions, HS images are always
suffering from the pollution of various noise. Fig. (4) illustrates different
noise types of HS images observed from airborne and spaceborne sensors.
Therefore, HS image denoising and restoration is a necessary pre-processing
for the subsequent applications to assist the noise.

Figure 4: Examples for different noise types of HS Images observed from
airborne and spaceborne sensors.

The statistical distribution of hyperspectral noise is complicated. For
instance, the readout noise, which is assumed to obey the Gaussian
distribution, is produced by the imaging device (Charge-coupled Device) during
the conversation from electrons to the final image [15]. The stripes are
generated in the hyperspectral data collected by the pushroom hyperspectral
sensors [16]. Due to the weather environment, the obtained HS data are always
suffering from the cloud and cloud shadow [17]. Besides, the HS images are
also suffering from the signal-dependent noise [18], multiplicative noise,
impulsive noise, Laplace noise and so on. Furthermore, In this paper, we
follow the mainstream and focus on the Gaussian noise removal [19, 20] and
mixed noise removal problem [21, 22].

Since 1980s, researchers have paid attention to the HS noise analysis. For
example, maximum noise fraction (MNF) transformation [23] noise-adjusted
principal component analysis [19] are utilized to extract high-quality
components for the subsequent classification and reject the low-quality
components. Following the mainstream of gray/color image denoising in computer
vision society, various state-of-the-art technologies, such as wavelets [24],
sparse representation [25, 26], TV [20, 27], non-local means processing [28,
29], low-rank matrix representation [21, 30, 31, 22], tensor representation
[32, 33, 34], DL [35, 36] and so on.

The non-convex regularized methods have also been developed for HS image
restoration. From [21] the HS images are assumed to be corrupted by the
additive noise, including Gaussian noise, stripes, deadlines, pixel missing,
impulse noise and so on. The observation model is formulated as:

$\boldsymbol{\mathcal{Y}}=\boldsymbol{\mathcal{X}}+\boldsymbol{\mathcal{N}}_{G}+\boldsymbol{\mathcal{N}}_{S},$
(3)

where $\boldsymbol{\mathcal{Y}}$ represents the observed noisy image,
$\boldsymbol{\mathcal{X}}$ stands for the latent clean image,
$\boldsymbol{\mathcal{N}}_{G}$ is the Gaussian noise, and
$\boldsymbol{\mathcal{N}}_{S}$ is the sparse noise, including stripes,
deadlines, pixel missing and impulse noise. Typically, when sparse noise
$\boldsymbol{\mathcal{N}}_{S}$ is omitted, the model (3) is degraded to the
Gaussian noise removal problem. From [21], the low-rank property exploration
of clean image $\boldsymbol{\mathcal{X}}$ has attracted much attention and
achieved state-of-the-art HS image denoising performance [29, 34]. Generally
speaking, the mainstreams are from two points. Firstly, how to explore the
low-rank property of clean image $\boldsymbol{\mathcal{X}}$. Until now,
spectral low-rank property [21, 30], spatial low-rank property [37, 33] and
non-local low-rank property [32, 29] have been well studied. Furthermore, how
to balance the contribution of different low-rank properties is also a key
problem [38, 29]. Secondly, the low-rank constraint of
$\boldsymbol{\mathcal{X}}$ formulates a non-convex optimization problem.
Therein, how to efficiently solve the rank constraint non-convex problem is
another key problem. In the next subsection, we review several outstanding
works and illustrate how these works formulate the low-rank modeling, and how
to solve the non-convex problem.

TABLE I: Prior properties of the selected methods. One, two, and three $\bullet$ denote the low, medium, and high intensity of prior information, respectively. Methods | Low-rankness | Local smoothness
---|---|---
spatial | spectral | non-local | spatial | spectral
LRTA | $\bullet\bullet$ | $\bullet\bullet$ |  |  |
NAILRMA |  | $\bullet\bullet$ |  |  |
TDL |  | $\bullet$ | $\bullet\bullet\bullet$ |  |
FastHyDe |  | $\bullet\bullet$ | $\bullet$ |  |
NGmeet |  | $\bullet\bullet$ | $\bullet\bullet$ |  |
LRMR |  | $\bullet\bullet$ |  |  |
LRTV |  | $\bullet\bullet$ |  | $\bullet\bullet$ |
LRTDTV | $\bullet$ | $\bullet\bullet$ |  | $\bullet\bullet$ | $\bullet\bullet$
LRTDGS | $\bullet$ | $\bullet\bullet\bullet$ |  | $\bullet\bullet\bullet$ |
LRTF-FS |  | $\bullet\bullet\bullet$ |  | $\bullet\bullet\bullet$ | $\bullet\bullet$
TABLE II: The restoration results of different selected methods on Gaussian noise and mixed noised, respectively. Index | Gaussian noise removal | Mixed noise removal
---|---|---
LRTA | NAILRMA | TDL | FastHyDe | NGmeet | LRMR | LRTV | LRTDTV | LRTDGS | LRTF-FS
PSNR | 25.99 | 32.81 | 32.11 | 33.51 | 33.82 | 32.22 | 33.05 | 32.34 | 33.33 | 33.26
SSIM | 0.7095 | 0.9519 | 0.9443 | 0.9601 | 0.9607 | 0.9401 | 0.9459 | 0.9335 | 0.9596 | 0.9549
MSA | 11.35 | 4.75 | 4.72 | 4.42 | 4.38 | 4.95 | 5.06 | 4.09 | 4.24 | 4.44

###  III-A Gaussian Noise Removal

In this section, five methods are selected to represent the state-of-the-art
HS image Gaussian noise removal approaches. These methods utilize different
low-rank matrix/tensor decomposition models to exploit the spatial, spectral,
or non-local low-rank properties of the original clean HS image. The
properties of these five methods are summarized in Table I. The five methods
are briefly described in the following.

####  III-A1 LRTA

111https://www.sandia.gov/tgkolda/TensorToolbox/

On the basis of the observation model (3) with ignoring sparse noise $S$, low-
rank tensor approximation (LRMA) [39] tries to restore the HS image from the
following objective function

$\min_{\boldsymbol{\mathcal{X}}}\|\boldsymbol{\mathcal{Y}}-\boldsymbol{\mathcal{X}}\|_{\operatorname{F}}^{2}\;\;{\rm
s.t.}\;\boldsymbol{\mathcal{X}}=\boldsymbol{\mathcal{O}}\times_{1}\mathbf{A}\times_{2}\mathbf{B}\times_{3}\mathbf{C},$
(4)

where
$\boldsymbol{\mathcal{X}}=\boldsymbol{\mathcal{O}}\times_{1}\mathbf{A}\times_{2}\mathbf{B}\times_{3}\mathbf{C},$
is the Tucker decomposition,
$\boldsymbol{\mathcal{O}}\in\mathbb{R}^{r_{1}\times r_{2}\times r_{3}}$ stands
for the core tensor, and $\mathbf{A}\in\mathbb{R}^{M\times
r_{1}},\mathbf{B}\in\mathbb{R}^{N\times
r_{2}},\mathbf{C}\in\mathbb{R}^{B\times r_{3}}$ are the factors related to
different dimensions. With the rank $(r_{1},r_{2},r_{3})$ of Tucker
decomposition set in advance, the LRTA model (4) can simultaneously capture
the global spatial and spectral low-rank properties. (4) provides a simple
general model for different kinds of low-rank matrix/tensor decomposition
based HS image denoising methods, that’s to say, we can change the Tucker
decomposition constraint of $\boldsymbol{\mathcal{X}}$ to different kinds of
matrix/tensor decomposition, such as canonical polyadic (CP) decomposition,
tensor train decomposition, tensor ring decomposition and so on.

####  III-A2 NAILRMA

222https://sites.google.com/site/rshewei/home

Noise-adjusted iterative LRMA (NAILRMA) [30] method assumes that the spectral
low-rank property is more important than that of spatial ones, and therefore
simply spectral low-rank regularizer is utilized to restrict the original
spectral image $\boldsymbol{\mathcal{X}}$. From other works [40], it also
indicates that the spatial TV regularizer is more important than that of the
spatial low-rank regularizer. In HS images, the similar signatures
representing the same class also appear in the nearby spatial location. To
enhance the spectral low-rank property, NAILRMA segmented the HS image into
spatial overlapping patches and process each patch individually. The noise
intensity in different bands of the HS image is different, which is a big
challenge appearing in the HS image Gaussian noise removal, is mitigated by
the noise-adjusted iterative strategy [30]. At last, the randomized singular
value decomposition (RSVD) is utilized to solve the non-convex low-rank
approximation problem.

####  III-A3 TDL

333http://gr.xjtu.edu.cn/web/dymeng/

Tensor dictionary learning (TDL) combines the non-local regularization and
low-rank tensor approximation. The noisy HS image is firstly segmented into
spatial overlapping patches, and the similar patches are clustered together to
formulate a higher-order tensor. In this way, the non-local spatial
information is collected. Then the higher-order tensors are denoised as the
same of (4), and finally the denoised higher-order tensors are utilized to
formulate the final denoised HS image. TDL represents the first method to
exploit non-local low-rank property, and the subsequent methods LLRT [38], KBR
[32], and NLTR [34] also achieve remarkable HS image Gaussian noise removal
results.

####  III-A4 FastHyDe

444https://github.com/LinaZhuang/FastHyDe_FastHyIn

The main difference between HS images and color/multispectral images is the
number of spectral bands. To eliminate this difference and utilize well
developed color/multispectral denoising methods for HS image denoising, Zhuang
et.al proposed the fast HS denoising (FastHyDe) [41] method by translating the
HS image to low-dimensional reduced image via SVD. By this translation,
various state-of-the-art color/multispectral image denoising methods, such as
wavelets [42] and BM3D [41] are used to denoise the reduced image. Finally,
the denoised reduced image is translated back to the denoised HS image via
inverse SVD. Generally speaking, under the framework of FastHyDe, the HS image
noise removal task is linked to the development of color/multispectral image
noise removal tasks.

####  III-A5 NGmeet

555https://github.com/quanmingyao/NGMeet

Spatial non-local low-rank regularizer can produce state-of-the-art HS image
noise removal performance. However, as the increase of spectral number, the
time cost of non-local related methods also increase incredibly [38, 32]. Non-
local meets global (NGmeet) method also tries to translate the HS image to the
reduced image, and utilizes a non-local low-rank method to denoise the reduced
image. Different from FastHyDe, NGmeet tries to perfect the framework by
iteratively eliminating the error caused by SVD on the noisy HS image, and
automatically estimating the spectral rank of the reduced image.

###  III-B Mixed Noise Removal

In this section, we select five representative methods for the HS image mixed
noise removal. These methods are on the basis of the observation model (3). We
focus on the non-convex low-rank regularizer of original image
$\boldsymbol{\mathcal{X}}$. The properties of these five methods are
summarized in Table I.

####  III-B1 LRMR

Zhang et al. firstly introduced the observation model (3) to analysis of
complex HS noise [21]. LRMR tries to restore the original clean image
$\boldsymbol{\mathcal{X}}$ from the noisy image via low-rank and sparse
decomposition model as follows:

$\min_{\mathbf{X}}\lVert\mathbf{Y}-\mathbf{X}-\mathbf{N}_{S}\rVert_{\operatorname{F}}^{2}+\lambda_{1}rank(\mathbf{X})+\lambda_{2}card(\mathbf{N}_{S}),$
(5)

where $\mathbf{Y},\mathbf{X},\mathbf{N}_{S}$ are the reshaped matrices of
$\boldsymbol{\mathcal{Y}},\boldsymbol{\mathcal{X}},\boldsymbol{\mathcal{N}}_{S}$
along the spectral dimension, respectively, $\lambda_{1}$ and $\lambda_{2}$
are the parameters to trade-off the contributions of $rank(\mathbf{X})$ and
non-zero elements $card(\mathbf{N}_{S})$. LRMR utilizes “GoDec” algorithm [43]
to alternatively update non-convex constraint $\mathbf{X}$ and
$\mathbf{N}_{S}$, and finally obtains the restored image. To improve the
efficiency of the optimization to (5), several non-convex substitutions, such
as reweighted nuclear norm [31], $\gamma$-norm [44], smooth rank approximation
[45] and normalized $\epsilon$-Penalty [46], are further developed to exploit
the spectral low-rank property of $\mathbf{X}$.

####  III-B2 LRTV

Low-rank total variation (LRTV) claimed that spectral low-rank property is not
enough to describe the property and HS images, and therein introduced TV to
explore the spatial smoothness via TV. Generally, low-rank regularization and
TV are the two most studied regularizers, and the combination of them to
produce the state-of-the-art HS image mixed noise removal is becoming popular.
Most of the following works try to either improve the low-rank modeling [47,
48] or the smoothness modeling [49, 33, 50] of the HS image to further improve
the restoration accuracy. To further combine low-rank and TV, the low-rank
exploration of the HS difference image is also developed [51, 52].

####  III-B3 LRTDTV

Low-rank tensor decomposition with TV (LRTDTV) [48] tries to improve LRTV by
utilizing low-rank tensor decomposition to exploit the low-rank property of HS
images, meanwhile spatial-spectral TV (SSTV) to explore the spatial and
spectral smoothness simultaneously. Although LRTDTV achieved better mixed
noise removal results as reported in [48], the spatial rank utilized in LRTDTV
is much larger than that of spectral rank. This is mainly because the spatial
low-rank property of HS images is not so important compared to the spectral
low-rank property. From another side, the spatial non-local low-rank
regularization is proved to be more efficient [40] than spatial low-rank
property for HS restoration problem.

####  III-B4 LRTDGS

666https://chenyong1993.github.io/yongchen.github.io/

Low-rank tensor decomposition with group sparse regularization (LRTDGS) [33]
also utilizes low-rank rank tensor decomposition to exploit the low-rank
property of HS images. Differently, LRTDGS explores the group sparsity of the
difference image instead of SSTV in LRTDTV. From the mathematical modeling,
LRTDGS utilizes weighted $ell_{2,1}$ norm regularization to fulfill the row-
group sparsity of the difference image.

####  III-B5 LRTF-FR

777https://yubangzheng.github.io/homepage/

Following the idea of NGmeet [29], factor-regularized low-rank tensor
factorization (LRTF-FR) [53] also utilizes matrix decomposition to decouple
the spatial and spectral priors. From one side, the spectral signatures of the
HS image are assumed to be of smooth structure. From another side, the reduced
image is assumed to have a group sparse structure in the difference domain.
The optimization model of LRTF-FR is

$\displaystyle\min_{\boldsymbol{\mathcal{X}},\boldsymbol{\mathcal{N}}_{S}}$
$\displaystyle\lVert\boldsymbol{\mathcal{Y}}-\boldsymbol{\mathcal{X}}-\boldsymbol{\mathcal{N}}_{S}\rVert_{\operatorname{F}}^{2}+\lambda_{1}\lVert\boldsymbol{\mathcal{X}}\times_{3}\mathbf{H}_{3}\rVert_{2,1}$
(6)
$\displaystyle+\lambda_{2}\sum_{k=1}^{2}\lVert\boldsymbol{\mathcal{X}}\times_{k}\mathbf{H}_{k}\rVert_{\operatorname{F}}^{2}+\lambda_{3}\lVert\boldsymbol{\mathcal{N}}_{S}\rVert_{1,1},$

where $\mathbf{H}_{k},k=1,2,3$ are the first-order difference matrices.
Furthermore, in the optimization to (6), the reweighted strategy is utilized
to update $\ell_{2,1}$ norm and $L\ell_{1}$ norm to further improve the
restored results.

###  III-C Experimental Study

We choose the HS image from DLR Earth Sensing Imaging Spectrometer (DESIS)
installed on the International Space Station (ISS)[54] for the experimental
study to compare different methods on the Gaussian and mixed noise removal
tasks. We remove the noisy bands and select a sub-image of size $400\times
400\times 200$ as the clean reference image, which is normalized to $[0,1]$.
Firstly, we add Gaussian noise of noise variance $0.1569$ to simulate the
Gaussian noisy image, and apply different Gaussian noise removal methods to
remove the Gaussian noise. Furthermore, we add salt & pepper noise and stripes
to simulate the mixed noisy image, and apply mixed noise removal methods to
remove the mixed noise. As similar in [33], we choose the mean of peak signal-
to-noise rate (PSNR) over all bands, the mean of structural similarity (SSIM)
over all bands, and the mean of spectral angle mapping (MSA) overall spectral
vectors to evaluate the restored results.

Figure 5: The illustration of different methods on the noise removal results.
The first row shows the results of different methods on the Guassian noise
remove experiments (R:70, G:100, B:36). The second row displays the results of
different methods on the mixed noise remove experiments (R:31, G:80, B:7).

Table II presents the evaluation values of different methods on Gaussian noise
and mixed noise removal results, respectively. For the Gaussian noise removal
task, NGmeet achieves the best values of three evaluation indices. However,
the gap between NGmeet and FastHyDe is limited. For the mixed noise removal
task, LRTDGS achieves the best accuracy in PSNR and SSIM values, meanwhile
LRTDTV achieved the best MSA value. Combining Tables I and II, we can conclude
that, firstly, the spectral low-rank prior information is important for HS
restoration. Secondly, the contribution of spatial low-rank prior information
for HS restoration is limited. Thirdly, on the basis of spectral low-rank
regularization, spatial and spectral smoothness prior can further improve the
final HS restoration results.

###  III-D Remaining Challenges

Up to date, many non-convex regularized methods have been proposed to develop
the low-rank priors and local smoothness priors, and achieve remarkable HS
restoration results for Gaussian and mixed noise removal. However, these
methods still face several challenges for further work. We summary these
challenges as the following.

  * •

Efficiency. Although low-rank related methods have achieved state-of-the-art
restoration results, they are time-consuming. For instance, NGmeet and LRTVGS
speed more than 10 minutes to process the HS image of size $400\times
400\times 200$. Furthermore, the related state-of-the-art restoration methods
always exploit multiple priors of the HS image, resulting in the confusion of
the parameter chosen. Therein, how to reduce the model complexity and improve
the optimization efficiency of the HS image restoration is the key challenge.

  * •

Scalability. Previous non-convex related methods always focus on the small HS
image processing. However, HS images are used to observe the earth and the
spatial size of one scene is usually very large. How to improve the
scalability of the restoration approaches is the key challenge. DL provides
the possibility for fast and large scale processing of HS images. Whereas DL
approaches always rely on the quality of training samples, and the applicable
scope of the test data is always limited. To improve the scalability, how to
embed well studied non-convex regularizers to the DL architecture should also
be further analyzed.

  * •

Real Application. Until now, most HS image restoration methods are evaluated
on the simulated experiments. However, in most cases, the evaluation indices
fail to predict the accuracy of the real HI image restoration results. From
another side, the noise distribution in the real noisy HS images is complex.
How to testify the related methods on the real HS images should be also
further analyzed. From another side, the training samples in the real
application are always limited. The blind and unsupervised approaches will
become the mainstream of future real HS image restoration.

##  IV Dimensionality Reduction

The HS dimensionality reduction (DR) and feature extraction have long been a
fundamental but challenging research topic in HS RS [55, 56]. The main reasons
mainly lie in the following aspects. Due to the highly-correlated
characteristic between spectral bands, the HS images are subjected to
information redundancy, which could hurt the ability to discriminate the
materials under the certain extremely-conditioned cases (curse of
dimensionality). Plus, as the HS dimension gradually increases along with the
spectral domain, large storage capability and high-performance computing are
needed. Furthermore, these dimension-reduced features are usually applied for
the high-level classification or detection task [57, 58]. Recently, many works
based on non-convex modeling have shown to be effective for automatically
extracting dimension-reduced features of HS images. Linking with Eq. (2), the
DR task can be generalized to the following optimization problem:

$\displaystyle\mathop{\min}_{f_{\Phi},\mathbf{X}}\frac{1}{2}\lVert
f_{\Phi}(\mathbf{Y})-\mathbf{X}\rVert_{\operatorname{F}}^{2}\;\;{\rm
s.t.}\;\mathbf{X},f_{\Phi}\in\mathcal{C},$ (7)

where $f_{\Phi}(\bullet)$ denotes the transformation from the original HS
space to dimension-reduced subspaces with the respect to the variable set
$\Phi$, and $\mathbf{X}$ is the low-dimensional representations of
$\mathbf{Y}$. Revolving around the general form in Eq. (7), we review
currently advanced DR methods from three different aspects: unsupervised,
supervised, and semi-supervised models.

Figure 6: An illustration for supervised DR models in HS images with two
different groups: discriminant analysis based DR and regression-induced DR.

###  IV-A Unsupervised Model

Non-negative matrix factorization (NMF) [59] is a common unsupervised learning
tool, which has been widely applied in HS DR. These works can be well
explained by Eq. (7), the NMF-based DR problem can be then formulated as

$\displaystyle\mathop{\min}_{\mathbf{Q}\geq\mathbf{0},\mathbf{X}\geq\mathbf{0}}\frac{1}{2}\lVert\mathbf{Y}-\mathbf{X}\mathbf{Q}\rVert_{\operatorname{F}}^{2}+\Psi(\mathbf{X})+\Omega(\mathbf{Q}),$
(8)

where $\mathbf{Q}$ denotes the combination coefficients, $\Phi(\mathbf{X})$
and $\Omega(\mathbf{Q})$ are the potential regularization terms for the
variables $\mathbf{X}$ and $\mathbf{Q}$, respectively. Until the current,
there have been some advanced NMF-based works in HS DR. Gillis et al. [60]
used sparse NMF under approximations for HS data analysis. Yan et al. [61]
proposed a graph-regularized orthogonal NMF (GONMF) model with the application
to spatial-spectral DR of HS images. Wen et al. [62] further extended the
GONMF with combining multiple features for HS DR. Rasti et al. [63] designed
an orthogonal total variation component analysis (OTVCA) approach for HS
feature extraction. Moreover, the HS data are directly regarded as a high-
dimensional tensor structure in [64], where the low-rank attribute is fully
considered in the process of low-dimensional embedding. In detail, we
summarize the regularization and constraints of the above methods as follows:

  * •

Sparsity [60]: $\Omega(\mathbf{Q})=\lVert\mathbf{Q}\rVert_{0}$;

  * •

Graph Regularization [61]:
$\Psi(\mathbf{X})=\operatorname{tr}(\mathbf{X}\mathbf{L}\mathbf{X}^{\top}),\;\rm{s.t.}\;\mathbf{X}\mathbf{X}^{\top}=\mathbf{I}$;

  * •

Multi-graph Regularization [62]:
$\Psi(\mathbf{X})=\sum_{i=1}^{s}\operatorname{tr}(\mathbf{X}\mathbf{L}^{i}\mathbf{X}^{\top}),\;\rm{s.t.}\;\mathbf{X}\mathbf{X}^{\top}=\mathbf{I}$;

  * •

Total Variation [63]:
$\Psi(\mathbf{X})=\sum_{i=1}^{r}\lVert\sqrt{(\mathbf{H}_{h}\mathbf{X}_{i})^{2}+(\mathbf{H}_{v}\mathbf{X}_{i})^{2}}\rVert_{1},\\\
\rm{s.t.}\;\mathbf{Q}\mathbf{Q}^{\top}=\mathbf{I}$;

  * •

Low-rank Graph [64]:
$\Psi(\mathbf{X})=\lVert\mathbf{X}\rVert_{*}+\operatorname{tr}(\mathbf{X}\mathbf{L}\mathbf{X}^{\top})$.

$\mathbf{L}=\mathbf{D}-\mathbf{W}$ is the Laplacian matrix, where
$\mathbf{D}_{i,i}=\sum_{j}\mathbf{W}_{i,j}$ is the degree matrix and
$\mathbf{W}$ is the graph (or manifold) structure of $\mathbf{X}$ [65].
$\lVert\bullet\rVert_{0}$, $\lVert\bullet\rVert_{2,1}$, and
$\lVert\bullet\rVert_{*}$ denote the $\ell_{0}$-norm [66], $\ell_{2,1}$-norm
[67], and nuclear norm [68], respectively.

Another type of unsupervised DR approaches is graph embedding, also known as
manifold learning, which can be also grouped into Eq. (7) well (according to
[69]):

$\displaystyle\mathop{\min}_{\mathbf{U},\mathbf{X}}\frac{1}{2}\lVert\mathbf{X}-\mathbf{U}\mathbf{Y}\rVert_{\operatorname{F}}^{2}+\Psi(\mathbf{X})+\Omega(\mathbf{U})\;{\rm
s.t.}\;\mathbf{X}\mathbf{X}^{\top}=\mathbf{I},$ (9)

where $\mathbf{U}$ denotes the to-be-estimated projection matrix that bridges
the high-dimensional data $\mathbf{Y}$ with the low-dimensional embedding
$\mathbf{X}$. The regularization term for the variable $\mathbf{U}$ can be
usually expressed as

$\displaystyle\Omega(\mathbf{U})=\operatorname{tr}(\mathbf{U}\mathbf{Y}\mathbf{L}\mathbf{Y}^{\top}\mathbf{U}^{\top})+\lVert\mathbf{U}\rVert_{\operatorname{F}}^{2},$
(10)

while the regularizer with respect to $\mathbf{X}$ can be given by

$\displaystyle\Psi(\mathbf{X})=\operatorname{tr}(\mathbf{X}\mathbf{L}\mathbf{X}^{\top}).$
(11)

The main difference between these manifold learning-based DR approaches lies
in the graph construction, i.e., $\mathbf{W}$. Ma et al. [70] integrated the
KNN classifier with several representative manifold learning algorithms, i.e.,
locally linear embedding [71], Laplacian eigenmaps [65], and local tangent
space alignment [72], for HS image classification. Huang et al. [73] embedded
the sparse graph structure, which is performed by solving a $\ell_{1}$-norm
optimization problem, for the DR of HS images. He et al. [74] extended the
work of [73] by generating a weighted sparse graph. Hong et al. [75] developed
a new spatial-spectral graph for the DR of HS images, called RLMR, by jointly
taking neighbouring pixels of a target pixel in spatial and spectral domains
into account. An et al. [76] attempted to learn the low-dimensional tensorized
HS representations on a sparse and low-rank graph. To sum up, core graphs of
the aforementioned methods can be obtained by

  * •

Sparse Graph [73]:
$\min_{\mathbf{W}}\lVert\mathbf{W}\rVert_{1,1},\;\rm{s.t.}\;\mathbf{Y}=\mathbf{Y}\mathbf{W}$;

  * •

Weighted Sparse Graph [74]:
$\min_{\mathbf{W}}\lVert\mathbf{d}\odot\mathbf{W}\rVert_{1,1},\;\rm{s.t.}\;\mathbf{Y}=\mathbf{Y}\mathbf{W},$
where $\mathbf{d}$ denotes a weighted matrix on $\mathbf{W}$ and $\odot$ is
the element-wise multiplication operator;

  * •

Spatial-spectral Graph [75]:
$\mathop{\min}_{\mathbf{w}_{i,0}}\sum_{j\in\phi_{i}^{spa}}\lVert\mathbf{y}_{i,j}-\sum_{k\in\phi_{i}^{spe}}\mathbf{y}_{i,k}w_{i,k,j}\rVert_{2}^{2}\\\
{\rm
s.t.}\;\lVert\sum_{k\in\phi_{i}^{spe}}\mathbf{y}_{i,k}(4w_{i,k,0}-\sum_{k=1}^{4}w_{i,k,j})\rVert_{2}^{2}\leq\eta,\\\
\qquad\mathbf{w}_{i,j}^{\operatorname{T}}\mathbf{w}_{i,j}=1,$
where $\phi_{i}^{spa}$ and $\phi_{i}^{spe}$ denote the neighbouring pixels in
spatial and spectral spaces, respectively;

  * •

Sparse and Low-rank Graph [76]:
$\min_{\mathbf{W}}\lVert\mathbf{W}\rVert_{1,1}+\lVert\mathbf{W}\rVert_{*},\;\rm{s.t.}\;\mathbf{Y}=\mathbf{Y}\mathbf{W}$.

###  IV-B Supervised Model

Unlike unsupervised DR that relies on embedding various priors to reduce the
dimension of HS data, supervised models are capable of better learning class-
separable low-dimensional representations via the use of label information.
The supervised DR models can be described from two different categories in
this subsection, as shown in Fig. 6. A typical group is the discriminant
analysis [55] closely related to graph embedding and manifold learning.
Intuitively speaking, these methods belong to a special case of unsupervised
graph embedding, which means they can be well explained by Eq. (9). The main
difference lies in that the labels are used for constructing the graph
structure, i.e., $\mathbf{W}$, thereby yielding a more discriminative
subspace.

In the supervised DR, a direct graph structure is written as

${\bf W}_{ij}=\begin{cases}\begin{aligned} 1,\;\;&\text{if ${\bf y}_{i}$ and
${\bf y}_{j}\in C_{k}$;}\\\ 0,\;\;&\text{otherwise,}\end{aligned}\end{cases}$
(12)

where $C_{k}$ means the sample set of the $k$-th class. Furthermore, more
advanced supervised graphs have been developed to better represent the HS data
in a low-dimensional subspace, such as sparse graph discriminant analysis
[77], collaborative graph discriminant analysis [78], feature space
discriminant analysis (FSDA) [79], spatial-spectral local discriminant
embedding [80]. These approaches sought to construct a soft graph instead of a
hard graph in Eq. (12). That is, the graph is built by using radial basis
function (RBF) to measure the similarity between samples belonging to the same
class [81]:

${\bf W}_{ij}=\begin{cases}\begin{aligned}
\exp\frac{-\lVert\mathbf{y}_{i}-\mathbf{y}_{j}\rVert_{2}^{2}}{2\sigma^{2}},\;\;&\text{if
${\bf y}_{i}$ and ${\bf y}_{j}\in C_{k}$;}\\\
0,\;\;&\text{otherwise,}\end{aligned}\end{cases}$ (13)

or by solving $\ell_{1}$-norm or $\ell_{2}$-norm optimization functions in the
same class set, e.g., [77], [78].

The DR behavior can be also modeled from a regression perspective by directly
connecting samples and labels [82], which provides a new insight into the
research of the supervised HS DR. A general form for the regression-based
supervised DR model can be formulated as

$\displaystyle\mathop{\min}_{\mathbf{P},\mathbf{U}}$
$\displaystyle\frac{1}{2}\lVert\mathbf{M}-\mathbf{P}\mathbf{X}\rVert_{\operatorname{F}}^{2}+\Psi(\mathbf{P})+\Omega(\mathbf{U})$
(14) $\displaystyle{\rm
s.t.}\;\mathbf{X}=\mathbf{U}\mathbf{Y},\;\mathbf{U}\mathbf{U}^{\top}=\mathbf{I},$

where the variable $\mathbf{P}$ denotes the regression coefficients or basis
signals, and $\mathbf{M}$ is the one-hot encoded matrix obtained by labels.
Eq. (14) can be, to some extent, regarded as an interpretable linearized
artificial neural network (ANN) mode (shallow network). Ji et al. [83] jointly
performed DR and classification, which is a good fit for Eq. (14) with
$\Psi(\mathbf{P})=\lVert\mathbf{P}\rVert_{\operatorname{F}}^{2}$. To enhance
the spectrally discriminative ability, Hong et al. [84] employed a LDA-like
graph on the basis of [83] to regularize the low-dimensional representations
in a Laplacian matrix form, i.e.,
$\Omega(\mathbf{U})=\operatorname{tr}(\mathbf{U}\mathbf{Y}\mathbf{L}\mathbf{Y}^{\top}\mathbf{U}^{\top})$.
In the same work [84], Hong et al. further extended their model to a deep
version, called JPlay, with a $k$-layered linear regression:

$\displaystyle\mathop{\min}_{\mathbf{P},\\{\mathbf{U}_{i}\\}_{i=1}^{k}}\frac{1}{2}\lVert\mathbf{M}-\mathbf{P}\mathbf{X}_{i}\rVert_{\operatorname{F}}^{2}+\Psi(\mathbf{P})+\Omega(\\{\mathbf{U}_{i}\\}_{i=1}^{k})$
(15) $\displaystyle{\rm
s.t.}\;\mathbf{X}_{i}=\mathbf{U}_{i}\mathbf{X}_{i-1},\;\mathbf{X}_{1}=\mathbf{U}_{1}\mathbf{Y},\;\mathbf{X}_{i}\geq\mathbf{0},\;\lVert\mathbf{x}_{i}\rVert_{2}\leq
1,$

with $\Psi(\mathbf{P})=\lVert\mathbf{P}\rVert_{\operatorname{F}}^{2}$ and

$\displaystyle\Omega(\\{\mathbf{U}_{i}\\}_{i=1}^{k})$
$\displaystyle=\sum_{i=1}^{k}\operatorname{tr}(\mathbf{U}_{i}\mathbf{X}_{i-1}\mathbf{L}\mathbf{X}_{i-1}^{\top}\mathbf{U}_{i}^{\top})$
$\displaystyle+\sum_{i=1}^{k}\lVert\mathbf{X}_{i-1}-\mathbf{U}_{i}^{\top}\mathbf{U}_{i}\mathbf{X}_{i-1}\rVert_{\operatorname{F}}^{2}.$

The J-Play attempts to open the “black box” of deep networks in an explainable
way by multi-layered linearized modeling. With explicit mappings and
physically meaningful priors, the non-convex J-Play takes a big step towards
the interpretable AI model.

###  IV-C Semi-supervised Model

Due to the fact that labeling samples is extremely expensive, particularly for
RS images covering a large geographic region, the joint use of labeled and
unlabeled information then becomes crucial in DR and classification.

A simple and feasible strategy for semi-supervised learning is to integrate
supervised and unsupervised techniques, e.g., LDA and locality preserving
projections [85]. By simultaneously constructing graphs of labeled and
unlabeled samples (e.g., using Eqs. (12) and (13), respectively), Eq. (9) can
be easily extended to a semi-supervised version, leading to semi-supervised
discriminant analysis (SSDA) [86]. Zhao et al. [87] further improved the SSDA
performance by using “soft” (or “pseudo”) labels predicted by label
propagation instead of directly using unsupervised similarities between
unlabeled samples. Similarly, Wu et al. [88] generated pseudo-labels using the
Dirichlet process mixing model and achieved a novel SSDA approach to learn the
low-dimensional HS embedding. These above-mentioned methods are performed
surrounding various hand-crafted graph structures ($\mathbf{W}$).

Figure 7: An example to clarify the graph structure of JPSA method, where
$\mathbf{W}^{p}$ and $\mathbf{W}^{sp}$ denote the pixel-wise and superpixel-
wise subgraphs as well as $\mathbf{W}^{a}$ is the aligned graph between pixels
and superpixels.

A different idea is to simulate the brain-like or human-like behaviors in the
semi-supervised DR task. It is well known that the feedback reward is a key
component that forms the intelligent information processing system. Inspired
by it, [89] developed an iterative multitask learning (IMR) framework by
adaptively learning the label propagation (LP) on graphs to simulate the
feedback mechanism, thereby achieving the HS DR process more effectively and
efficiently. The IMR is a semi-supervised extension of Eq. (14) with graph
learning, which can be generally modeled as

$\displaystyle\mathop{\min}_{\mathbf{P},\mathbf{U},\mathbf{L}}$
$\displaystyle\sum_{j=1}^{2}\lVert\mathbf{M}_{j}-\mathbf{P}\mathbf{U}\mathbf{Y}_{j}\rVert_{\operatorname{F}}^{2}+\Psi(\mathbf{P})+\Omega(\mathbf{U})$
(16)
$\displaystyle\mathrm{s.t.}\;\mathbf{U}\mathbf{U}^{\top}=\mathbf{I},\;\mathbf{L}\in\mathcal{C},$

where $\mathbf{Y}_{1}$ and $\mathbf{Y}_{2}$ denote the labeled and unlabeled
samples from $\mathbf{Y}$, respectively.
$\Psi(\mathbf{P})=\lVert\mathbf{P}\rVert_{\operatorname{F}}^{2}$ and
$\Omega(\mathbf{U})=\operatorname{tr}(\mathbf{U}\mathbf{Y}\mathbf{L}\mathbf{Y}^{\top}\mathbf{U}^{\top})$.
The non-convex constraint $\mathcal{C}$ with the respect to the variable
$\mathbf{L}$ can be summarized as

$\displaystyle\mathcal{C}:=\\{\mathbf{L}=\mathbf{L}^{\top},\;\mathbf{L}_{p,q,p\neq
q}\preceq 0,\;\mathbf{L}_{p,q,p=q}\succeq
0,\;\operatorname{tr}(\mathbf{L})=c\\},$

where $c>0$ is a scaling constant. Eq. (16) is a typical data-driven graph
learning model, which is capable of automatically learning the graph structure
from the data without any hand-crafted priors. By using the iterative strategy
to simulate the feedback system, $\mathbf{M}_{2}^{(t+1)}$ in the
$t$$+$$1$-step can be updated by the graph-based LP on the learned graph of
the $t$-step $\mathbf{W}^{(t)}$:

$\displaystyle\cdots\cdots\mathbf{M}_{2}^{(t+1)}\leftarrow\mathbf{W}^{(t)}\leftarrow\mathbf{M}_{2}^{(t)}\cdots\cdots.$
(17)

Besides, another intelligent feature extraction algorithm, named JPSA, which
is extended from [84], was presented in [90] by the attempts to align pixels
and superpixels for spatial-spectral semi-supervised HS DR. JPSA basically
follows the JPlay framework and the major difference is the graph structure
$\mathbf{W}$. The graph in JPSA consists of not only pixel-wise and
superpixel-wise similarities but also aligned components between pixels and
superpixels. Fig. 7 gives an example to clarify the graph structure of JPSA.
Note that the JPSA’s graph can be seen as a full data-driven structure, which
can, to a great extent, self-express the intrinsic properties of HS data and
further achieves intelligent information extraction and DR.

TABLE III: Quantitative comparison of different DR algorithms in terms of OA,
AA, and $\kappa$ using the NN classifier on the Indian Pines dataset. The best
one is shown in bold.

Methods | dimension | OA (%) | AA (%) | $\kappa$
---|---|---|---|---
OSF | 220 | 65.89 | 75.71 | 0.6148
OTVCA [63] | 16 | 74.18 | 77.61 | 0.7228
RLMR [75] | 20 | 83.75 | 86.90 | 0.8147
FSDA [79] | 15 | 64.14 | 74.52 | 0.5964
JPlay [84] | 20 | 83.92 | 89.35 | 0.8169
IMR [89] | 20 | 82.80 | 86.27 | 0.8033
JPSA [90] | 20 | 92.98 | 95.40 | 0.9197

###  IV-D Experimental Study

Classification is explored as a potential application to evaluate the
performance of state-of-the-art (SOTA) DR algorithms, including original
spectral features (OSF), OTVCA888https://github.com/danfenghong/HyFTech [63],
RLMR999https://github.com/danfenghong/IEEE_JSTARS_RLML [75], FSDA [79],
JPlay101010https://github.com/danfenghong/ECCV2018_J-Play [84], IMR [89], and
JPSA [90]. Experiments are performed on the Indian Pine data using the nearest
neighbor (NN) classifier in terms of three indices: Overall Accuracy (OA),
Average Accuracy (AA), and Kappa Coefficient ($\kappa$). The scene consists of
$145\times 145$ pixels and $220$ spectral bands ranging from $0.4\mu m$ to
$2.5\mu m$. More details regarding training and test samples can be found in
[91].

Table III lists the quantitative results of different DR methods. Overall, OSF
without feature extraction or DR yields the worst classification performance,
compared to those SOTA DR methods. This, to a great extend, demonstrates the
effectiveness and necessity of DR in the HS image classification task. It is
worth noting that the approaches with spatial-spectral modeling, e.g., OTVCA,
RLMR, JPSA, tend to obtain better classification results. The performance of
RLMR is superior to that of OTVCA, owing to the full consideration of the
neighboring information in a graph form rather than the smoothing operation
only modeled by the TV regularization. As a linearized deep model, supervised
JPlay obviously performs better than others, especially FSDA that is also a
supervised DR model. More importantly, the JPSA with a semi-supervised
learning strategy dramatically outperforms other competitors, since it can
jointly learn richer representations from both pixels and superpixels by means
of spatial-spectral manifold alignment and deep regressive regression.

###  IV-E Remaining Challenges

Although extensive SOTA methods have recently shown the effectiveness and
superiority in the HS DR and classification, there is still a long way to go
towards the AI-guided intelligent information processing. We herein summarize
the potential remaining challenges briefly.

  * •

Optimal Subspace Dimension. Subspace dimension is a crucial parameter in DR,
which is determined experimentally and empirically in most of existing
methods. Despite some parameter estimation algorithms, e.g., intrinsic
dimension [92], subspace identification [93], they fail to avoid the pre-
survey of various prior knowledge and human intervention in the dimension
estimation process.

  * •

Effects of Noises. HS images usually suffer from noise degradation in remotely
sensed imaging. These noises are complex and closely associated with spectral
signatures. Therefore, separating noises from HS data effectively and reducing
the noise sensitivity (or preserving spectral discrimination) in the DR
process remains challenging.

  * •

Robustness and Generalization. Robust estimation and advanced non-convex
regularizers have been widely applied to model the DR behavior, yet the
complex noise type, the limited training samples, and noisy labels hinder the
robustness and generalization ability to be further improved. For this reason,
more robust and intelligent models should be developed in either theory or
practice emphatically in the next generation DR technique.

Figure 8: A showcase in a real HS scene (Pavia City Centre) to have a quick
look at the 3-D HS cube, spectral signals, and material mixture as well as
pure pixels (i.e., endmember) and mixed pixels. In the studied scene, the pure
pixels correspond to two spectral reflectance curves of vegetation and water,
respectively, while the examples of mixed pixels explain the case of spectral
mixing, e.g., the two mixed pixels comprise of three pure components
(endmembers) in varying proportions. In addition, the figure located in the
right upper gives two toy examples to illustrate the material miscibility.

##  V Spectral Unmixing

Spectral unmixing (SU) can be usually seen as a special case of blind source
separation (BSS) problem in ML, referring to a procedure that decomposes the
observed pixel spectrum of the HS image into a series of constituent spectral
signals (or endmembers) of pure materials and a set of corresponding abundance
fractions (or abundance maps) [94]. Due to the meter-level ground sampling
distance (GSD) of HS imaging, the spectral signatures for most pixels in HS
images are acquired in the form of a complex mixture that consists of at least
two types of materials. Fig. 8 gives a showcase to visualize the HS cube,
spectral signatures, and material mixing process as well as pure and mixed
pixels. Different from general signals, e.g., digital signals, speech signals,
there are specific absorption properties in the spectrum signals of different
materials. Plus, HS images suffer from miscellaneous unknown degradation,
either physically or chemically, in the remotely sensed imaging, inevitably
bringing many uncertainties in SU. Therefore, SU plays a unique role in HS RS,
yielding many challenging researchable tasks compared to BSS in ML.

Ideally, a linear mixing model (LMM) can be used to accurately describe the SU
process [95], which is modeled as the following constrained optimization
problem:

$\displaystyle\mathop{\min}_{\mathbf{E},\mathbf{A}}\frac{1}{2}\lVert\mathbf{Y}-\mathbf{E}\mathbf{A}\rVert_{\operatorname{F}}^{2}\;\;\rm{s.t.}\;\mathbf{E},\mathbf{A}\in\mathcal{C}.$
(18)

The variables $\mathbf{E}$ and $\mathbf{A}$ in Eq. (18) stand for the
endmembers and abundance maps in the SU issue, respectively. According to the
endmembers ($\mathbf{E}$) that are available (given) or not in the process of
SU, existing SU models can be loosely divided into blind SU and endmember-
guided SU.

###  V-A Blind Spectral Unmixing

NMF is a baseline model in a wide range of applications, and the same is true
in SU. Up to the present, NMF-based interpretable models have been developed
extensively for pursing the intelligent SU with the consideration of
physically meaningful priors with respect to $\mathbf{E}$ and $\mathbf{A}$,
e.g., the abundance non-negative constraint (ANC), the abundance sum-to-one
constraint (ASC). The resulting basic blind SU model can be written as

$\displaystyle\mathop{\min}_{\mathbf{E},\mathbf{A}}\frac{1}{2}\lVert\mathbf{Y}-\mathbf{E}\mathbf{A}\rVert_{\operatorname{F}}^{2}+\Phi(\mathbf{E})+\Omega(\mathbf{A})\;\;\rm{s.t.}\;\mathbf{E},\mathbf{A}\in\mathcal{C},$
(19)

where the constraint $\mathcal{C}$ is

$\displaystyle\mathcal{C}:=\\{\mathbf{E}\geq\mathbf{0},\;\mathbf{A}\geq\mathbf{0},\;\mathbf{1}^{\top}\mathbf{A}=\mathbf{1}\\}.$

On the basis of the model (19), Yang et al. [96] proposed sparse NMF for SU
with a well-designed S-measure sparseness. Qian et al. [97] imposed the
sparsity constraint on abundances and used $\ell_{1/2}$-regularized NMF for
blind SU, which has shown to be more effective than $\ell_{0}$\- and
$\ell_{1}$-norm terms. In [98], Sigurdsson et al. relaxed $\ell_{1/2}$-norm to
$\ell_{q}$-norm ($0\leq q\leq 1$) for a better estimation of abundances.
Thouvenin et al. [99] developed an improved LMM, called perturbed LMM (PLMM),
by the attempts to model spectral variabilities as perturbed information that
simply meets the Gaussian distribution. A similar work is presented in [100],
where the scaling factor, as a major spectral variability (SV), is modeled
into LMM to yield an extended LMM (ELMM) for the blind SU task. He et al.
[101] employed total variation (TV) and weighted $\ell_{1}$-norm terms to
further enhance the smoothness and sparseness of abundances. Yao et al. [102]
sought to explain the NMF-based SU model by simulating human observations on
HS images, e.g., sparsity, non-local, smooth properties, in a non-convex
modeling fashion. Another type of interesting SU strategy is to embed the
graph or topological structure of the HS data. The local neighboring relation
is introduced into the NMF model, showing robust SU results [103]. Similarly,
Lu et al. [104] enforced the abundances to follow the manifold structure of
spectral signatures in the form of Laplacian regularization form for HS
unmixing. Wang et al. [105] used a structuralized hypergraph regularization in
sparse NMF to better depict the underlying manifolds of the HS data. Very
recently, Qin et al. [106] proposed a novel graph TV regularization to
estimate endmembers and abundances more effectively. There are still other
variants that directly unmix the 3-D HS tensor by preserving spatial structure
information as much as possible. For that, Qian et al. [107] proposed a
matrix-vector non-negative tensor factorization framework for blind SU.
Imbiriba et al. [108] modeled the low-rank properties in the HS tensor to
address the SV for robust SU. A further modified work based on [108] is
proposed via weighted non-local low-rank tensor decomposition for sparse HS
unmixing.

Broadly, these key non-convex priors of the above models can be briefly
summarized as follows:

  * •

$\ell_{1/2}$-NMF [97]:
$\Omega(\mathbf{A})=\lVert\mathbf{A}\rVert_{1/2}=\sum_{k,n=1}^{K,N}\mathbf{a}_{n}(k)^{1/2}$;

  * •

$\ell_{q}$-NMF [98]:
$\Omega(\mathbf{A})=\lVert\mathbf{A}\rVert_{q}=\sum_{k,n=1}^{K,N}\mathbf{a}_{n}(k)^{q}$;

  * •

PLMM [99]:
$\Phi(\mathbf{E})=\frac{1}{2}\lVert\mathbf{E}-\mathbf{E}_{0}\rVert_{\operatorname{F}}^{2}$,
$\Omega(\mathbf{A})=\frac{1}{2}\lVert\mathbf{A}\mathbf{H}\rVert_{\operatorname{F}}^{2}$,
$\Psi(\Delta)=\frac{1}{2}\sum_{n=1}^{N}\lVert\Delta_{n}\rVert_{\operatorname{F}}^{2}$,
where $\mathbf{E}_{0}$, $\mathbf{H}$, and $\Delta$ denote the reference
endmembers, the matrix differences in spatial four nearest neighbors, and
pixel-wise perturbed information, respectively.

  * •

ELMM [100]:
$\Phi(\mathbf{E})=\sum_{n=1}^{N}\lVert\mathbf{E}_{n}-\mathbf{E}_{0}\mathbf{S}_{n}\rVert_{\operatorname{F}}^{2}$,
$\Omega(\mathbf{A})=\lVert\mathbf{H}_{h}(\mathbf{A})\rVert_{2,1}+\lVert\mathbf{H}_{v}(\mathbf{A})\rVert_{2,1}$,
$\Psi(\mathbf{S})=\lVert\mathbf{H}_{h}(\mathbf{S})\rVert_{\operatorname{F}}^{2}+\lVert\mathbf{H}_{v}(\mathbf{S})\rVert_{\operatorname{F}}^{2}$,
where $\mathbf{H}_{h}$ and $\mathbf{H}_{v}$ are the horizontal and vertical
gradients;

  * •

TV-RSNMF [101]:
$\Omega(\mathbf{A})=\lVert\mathbf{d}\odot\mathbf{A}\rVert_{1,1}+\lVert\mathbf{A}\rVert_{\rm
TV}$;

  * •

NLHTV [102]:
$\Omega(\mathbf{A})=\sum_{n=1}^{N}\lVert
J_{w}\mathbf{a}_{n}\rVert_{\mathcal{S}_{1}}+\sum_{i,j}\log(|x_{i,j}|+\epsilon)$,
where $J_{w}$ and $\lVert\bullet\rVert_{\mathcal{S}_{1}}$ are defined as non-
local Jacobian operator and the Schatten-1 norm, respectively.

  * •

Graph $\ell_{1/2}$-NMF [104]:
$\Omega(\mathbf{A})=\lVert\mathbf{A}\rVert_{1/2}+\operatorname{tr}(\mathbf{A}\mathbf{L}\mathbf{A}^{\top})$;

  * •

Graph TV [106]: $\Omega(\mathbf{A})=\lVert\mathbf{A}\rVert_{\rm
TV}+\operatorname{tr}(\mathbf{A}\mathbf{L}\mathbf{A}^{\top})$.

Owing to the powerful data fitting ability, DL-based SU approaches have
recently been paid increasing attention and achieved better unmixing results
[109, 110, 111, 112]. Although these methods still suffer from the effects of
“black box”, i.e., the lack of model interpretability, yet their performances
have preliminary shown the effectiveness and feasibility in unmixing the HS
data more accurately.

###  V-B Endmember-Guided Spectral Unmixing

A mass of blind SU methods has been developed and shown to be effective to
simultaneously obtain endmembers and abundance maps. However, these blind
methods tend to extract physical meaningless endmembers, e.g., noisy signals,
spectral signatures corresponding to non-existent materials, due to the lack
of certain interpretable model guidance or prior knowledge. A straightforward
solution is to provide nearly real endmembers extracted from HS images. This
naturally leads to the researches on endmember-guided SU. As the name
suggests, the SU process is performed with given reference endmembers or the
guidance of extracted endmembers from the HS image. That is, the endmembers
$\mathbf{E}$ in Eq. (19) are known. Accordingly, the endmember-guided SU can
be implemented in a three-stage way.

  * •

Firstly, the number of endmembers can be estimated by subspace estimation
algorithms, e.g., HySime [93];

  * •

Secondly, the endmembers can be extracted based on geometric observations of
HS data structure. Several well-known methods are vertex component analysis
(VCA) [113], pixel purity index (PPI) [114], and fast autonomous endmember
extraction (N-FINDER) [115].

  * •

Lastly, the abundances of materials are estimated using regression-based
methods, which can generally written as

$\displaystyle\mathop{\min}_{\mathbf{A}}\frac{1}{2}\lVert\mathbf{Y}-\mathbf{E}\mathbf{A}\rVert_{\operatorname{F}}^{2}+\Omega(\mathbf{A})\;\;{\rm
s.t.}\;\mathbf{A}\in\mathcal{C}.$ (20)

Following the three steps, many well-working non-convex models have been
successfully developed to estimate the abundance maps of different materials
at a more accurate level. Heinz et al. [116] thoroughly analyzed the spectral
mixture in the SU issue, yielding a fully constrained least-squares unmixing
(FCLSU) algorithm. Due to the hard ASC, the abundances can not be fully
represented in a simplex. For this reason, a partial constraint least-squares
unmixing (PCLSU) [117] model emerges as required without ASC. Bioucas-Dias et
al. [118] relaxed the strong $\ell_{0}$-norm to the solvable $\ell_{1}$-norm
in the sparse HS unmixing model and designed a fast and generic optimization
algorithm based on the ADMM framework [14], called sparse unmixing by variable
splitting and augmented Lagrangian (SUnSAL). In [119], a TV spatial
regularization is considered to further enhance the unmixing performance.
Iordache et al. [120] extended the sparse regression model to the
collaborative version regularized by $\ell_{2,1}$-norm for SU. Fu et al. [121]
proposed a semi-blind HS unmixing model by correcting the mismatches between
estimated endmembers and pure spectral signatures from the library. Huang et
al. [122] jointly imposed sparsity and low-rank properties on the abundances
for better estimating abundance maps. Hong et al. [123] devised an interesting
and effective subspace-based abundance estimation model. The model neatly
sidesteps to directly decompose the HS data in the complex high-dimensional
space instead of projecting the HS data into a more robust subspace, where the
SV tends to be removed in a more generalized way with low-rank attribute
embedding. Beyond the current framework, Hong et al. [124] further augmented
the basic LMM by fully modeling SVs, e.g., principal scaling factors and other
SVs that should be incoherent or low-coherent with endmembers, in order to
yield an interpretable and more intelligent SU model, called augmented LMM
(ALMM).

Figure 9: A visual example to clarify SVs in a real HS scene (Pavia City
Centre). An image patch cropped from the scene is select to show the spectral
bundles involving spectral variations of trees in (a). (b) shows a pure
spectral signature (i.e., endmember) of trees acquired from the laboratory.
(c) represents the differences between (a) and (b), which is seen as SVs.
Figure 10: Visualization of abundance maps estimated by different SOTA SU
algorithms on the Urban data, where SAM is computed to generate the
classification-like maps regarded as the GT to measure the shape similarity of
abundance maps obtained by different SU methods.

The non-convexity of these methods on priors, constraints, or modeling can be
summarized as follows:

  * •

FCLSU [116]: $\mathbf{A}\geq\mathbf{0}$,
$\mathbf{1}^{\top}\mathbf{A}=\mathbf{1}$;

  * •

PCLSU [117]: $\mathbf{A}\geq\mathbf{0}$;

  * •

SUnSAL [118]: $\Omega(\mathbf{A})=\lVert\mathbf{A}\rVert_{1,1}$,
$\mathbf{A}\geq\mathbf{0}$, $\mathbf{1}^{\top}\mathbf{A}=\mathbf{1}$;

  * •

SUnSAL-TV [119]:
$\Omega(\mathbf{A})=\lVert\mathbf{A}\rVert_{1,1}+\lVert\mathbf{A}\rVert_{\rm
TV}$, $\mathbf{A}\geq\mathbf{0}$;

  * •

CSR [120]:
$\Omega(\mathbf{A})=\lVert\mathbf{a}\rVert_{2,1}=\sum_{n=1}^{N}\lVert\mathbf{a}_{n}\rVert_{2}$,
$\mathbf{A}\geq\mathbf{0}$;

  * •

DANSER [121]:
$\Omega(\mathbf{A})=\sum_{n=1}^{N}(\lVert\mathbf{a}_{n}\rVert_{2}^{2}+\tau)^{p/2}$,
$\mathbf{A}\geq\mathbf{0}$,
$\Phi(\mathbf{E})=\lVert\mathbf{E}-\mathbf{E}_{0}\rVert_{\operatorname{F}}^{2}$;

  * •

SULoRA [123]:
$\Psi(\mathbf{U})=\lVert\mathbf{Y}-\mathbf{U}\mathbf{Y}\rVert_{\operatorname{F}}^{2}+\lVert\mathbf{U}\rVert_{*}$,
$\Omega(\mathbf{A})=\lVert\mathbf{A}\rVert_{1,1}$, $\mathbf{A}\geq\mathbf{0}$,
where $\mathbf{U}$ denotes the subspace projection and
$\lVert\bullet\rVert_{*}$ is the nuclear norm that approximates the rank
property of the matrix $\bullet$;

  * •

ALMM [124]: $\Phi(\mathbf{A})=\lVert\mathbf{A}\rVert_{1,1}$,
$\mathbf{A}\geq\mathbf{0}$,
$\Gamma(\mathbf{J})=\lVert\mathbf{J}\rVert_{\operatorname{F}}^{2}$,
$\Psi(\mathbf{V})=\lVert\mathbf{A}^{\top}\mathbf{V}\rVert_{\operatorname{F}}^{2}+\lVert\mathbf{V}^{\top}\mathbf{V}-\mathbf{I}\rVert_{\operatorname{F}}^{2}$,
where $\mathbf{V}$ and $\mathbf{J}$ denote the SV dictionary and corresponding
coefficients, respectively.

###  V-C Experimental Study

A real urban HS data acquired by the HYDICE over the urban area, Texas, USA,
in 2015 (the latest version111111http://www.tec.army.mil/Hypercube) is used to
evaluate the performance of several selected SOTA unmixing methods
qualitatively, including $\ell_{1/2}$-NMF [97],
PLMM121212https://pthouvenin.github.io/unmixing-plmm/ [99],
ELMM131313https://openremotesensing.net/knowledgebase/spectral-variability-
and-extended-linear-mixing-model/ [100], NLHTV [102], FCLSU [116],
SUnSAL141414http://www.lx.it.pt/~bioucas/ [118],
SULoRA151515https://github.com/danfenghong/IEEE_JSTSP_SULoRA [123], and
ALMM161616https://github.com/danfenghong/ALMM_TIP [124]. The HS image consists
of $307\times 307$ pixels and 162 spectral bands after removing noisy bands in
the wavelength range of $0.4\mu m$ to $2.5\mu m$ at a $2m$ GSD. Moreover, four
main materials (or endmembers) are investigated in the studied scene, i.e.,
asphalt, grass, trees, and roof. Furthermore, HySime [93] and VCA [113]
algorithms are adopted to determine the number of endmembers and extract
endmembers from the HS image (as the initialization for blind SU methods) for
all compared algorithms, respectively.

Fig. 10 shows the visual comparison between different SOTA unmixing algorithms
in terms of abundance maps. Owing to the consideration of real endmembers
extracted from the HS scene, the last four endmember-guided SU methods perform
evidently better than the blind SU ones. ELMM models the scaling factors,
tending to better capture the distributions of different materials. The
embedding of non-local spatial information makes the NLHTV method obtain a
more similar shape of abundance maps to the GT, yielding comparable unmixing
performance with ELMM. Remarkably, the unmixing results with regard to
abundance maps of SULoRA and ALMM algorithms are superior to those of other
methods, since the SVs can be fully considered by robustly embedding low-rank
attributes in a latent subspace using SULoRA and characterizing complex real
scenes more finely using ALMM.

###  V-D Remaining Challenges

SU has long been a challenging and widely concerned topic in HS RS. Over the
past decades, tons of SU works have been proposed by the attempts to unmix
these mixed spectral pixels more effectively. Yet, some key and essential
issues and challenges still remain to be solved.

  * •

Benchmark Data. Unlike classification, recognition, and detection tasks, the
ground truth of material abundances is able to be hardly collected, due to the
immeasurability of abundance values in reality. On the other hand, spectral
signatures (i.e., endmembers) of pure materials are often acquired in the lab.
This usually leads to uncertain mismatches between real endmembers and lab
ones. It turns to be urgent to establish benchmark datasets for SU by drawing
support from more advanced imaging techniques or developing interpretable
ground truth generation models or processing chain.

  * •

Evaluation Criteria. Reconstruction errors (RE) or spectral angle mapper (SAM)
are the two commonly used evaluation indices in SU. It should be noted,
however, that the results of RE or SAM are not equivalent to those of
unmixing. Linking to the issue of benchmark data, the measurement between real
results and estimated ones is the optimal choice, if we have the ground truth
for abundances and endmembers. If not, developing meaningful and reasonable
evaluation indices (e.g., classification accuracy) should give the top
priority in future work.

  * •

Spectral Variability. Spectral signatures inevitably suffer from various SVs
caused by illumination and topography change, noise effects from external
conditions and internal equipment, atmospheric interference, and complex
mixing of materials in the process of imaging. Fig. 9 shows a visual example
to specify the SVs (e.g., trees) in a real HS scene. Considerable
uncertainties brought by these factors have a big negative impact on accurate
estimation of abundances and endmembers in SU.

  * •

Nonlinearity. The complex interactions (e.g., intimate mixing, multilayered
mixing [94]) between multiple materials, also known as nonlinearity,
inevitably occur in the process of HS imaging. The nonlinearity in SU is a
longstanding and pending challenge. Most of existing nonlinear unmixing models
only attempt to consider certain special cases [125], e.g., bilinear mixing,
intimate mixtures, etc. Consequently, there is still lack of a general and
powerful model that can robustly address various nonlinearities in SU.

  * •

Model Explainability. The non-negativity and the sum-to-one constraint
considered in LMM are the basic priors for spectral signals in HS images.
However, only the two constraints fail to model the complex unmixing process
in an explainable fashion. To further enhance the explainability, new spectral
mixture models should be developed beyond the classic LMM by fully excavating
the intrinsic attribute knowledge that lies in the HS image.

##  VI Data Fusion and Enhancement

The high spectral resolution of HS images enables the identification and
discrimination of materials, meanwhile the high spatial resolution can provide
the possibility of the derivation of surface parameters [126]. However, due to
the equipment limitation, there is usually a trade-off between the spatial and
spectral resolutions, and the HS images obtained by the spaceborne imaging
spectrometers are usually with a moderate ground sampling distance [126]. To
enhance the spatial resolution, one popular way is to fuse the HS images with
high spatial MS images to generate new high spatial-spectral HS (HrHS) images.
In particular, enormous effects have been recently made to enhance the spatial
or spectral resolutions of HS images by means of ML techniques. Fig. 11
illustrates the fusion process of HS-MS images to generate the HrHS image.
Suppose we have the low-spatial resolution HS image
$\boldsymbol{\mathcal{Y}}\in\mathbb{R}^{m\times n\times B}$ and high-spatial
resolution MS image $\boldsymbol{\mathcal{Z}}\in\mathbb{R}^{M\times N\times
b}$ with $M\gg m$, $N\gg n$ and $B\gg b$, the fusion purpose is to generate
the high-spatial resolution HS image
$\boldsymbol{\mathcal{X}}\in\mathbb{R}^{M\times N\times B}$. The degradation
models from $\boldsymbol{\mathcal{X}}$ to $\boldsymbol{\mathcal{Y}}$ and
$\boldsymbol{\mathcal{Z}}$ are formulated as

$\displaystyle\mathbf{Y}=\mathbf{X}\mathbf{R}+\mathbf{N}_{H}$ (21)
$\displaystyle\mathbf{Z}=\mathbf{G}\mathbf{X}+\mathbf{N}_{M}$ (22)

where $\mathbf{X},\mathbf{Y},\mathbf{Z}$ are the reshaped matrices along the
spectral dimension of
$\boldsymbol{\mathcal{X}},\boldsymbol{\mathcal{Y}},\boldsymbol{\mathcal{Z}}$,
respectively, $\mathbf{R}$ is the mixed cyclic convolution and downsampling
operator, $\mathbf{G}$ is the spectral response function (SRF) of the MS image
sensor, $\mathbf{N}_{H}$ and $\mathbf{N}_{M}$ are the corresponding MS-HS
noise. To unify different observation models [127, 128, 129, 130, 126, 131],
$\mathbf{N}_{H}$ and $\mathbf{N}_{M}$ are assumed to be the independent
identically distributed Gaussian noise. Via the maximum a posteriori (MAP)
estimation method and Bayes rule [127, 129, 130], the following non-convex
optimization model is obtained

$\displaystyle\min_{\mathbf{X}}\|\mathbf{Y}-\mathbf{X}\mathbf{R}\|_{\operatorname{F}}^{2}+\|\mathbf{Z}-\mathbf{G}\mathbf{X}\|_{\operatorname{F}}^{2},$
(23)

where $\mathbf{R}$ and $\mathbf{G}$ are assumed to be known (in [129, 126],
$\mathbf{R}$ and $\mathbf{G}$ are estimated in advance of the optimization).
and As mentioned in [129, 130], the optimization of (23) is a NP-hard problem,
and over-estimation of $\mathbf{Z}$ will result in the unstable fusion
results. Therein, additional property of $\mathbf{X}$ and prior regularizers
should be exploited in the optimization model (23). It should be noted,
however, that the two functions $\mathbf{R}$ and $\mathbf{G}$ can be given
according to known sensors and also can be learned or automatically estimated
from the data itself.

Figure 11: Illustration of MS-HS image fusion to generate the HrHS image.

HS pansharpening is a heuristic way to perform the HS-MS fusion [132], which
has been widely applied in the HS image enhancement task. Component
substitution (CS) and multiresolution analysis (MRA) are the two main types of
pansharpening techniques. The former one aims to inject detailed information
of MS images into the low-resolution HS image, thereby generating the high-
resolution HS product. The latter one is to pansharpen the HS image by
linearly combining MS bands to synthesize a high-resolution HS band using
regression techniques. Another group for the HS-MS fusion task is the
subspace-based model, which roughly consists of Bayesian and unmixing based
methods (see [126]). Different from pansharpening, the subspace-based
approaches project the to-be-fused MS and HS images to a new space where the
dimension is generally smaller than that of the unknown high-resolution HS
image, by the means of the probability-driven Bayesian estimation (Bayesian-
based methods) or SU-guided matrix joint factorization (unmixing-based
methods).

In the following, we focus on the subspace methods, and review the related HS-
MS image fusion methods from the non-convex modeling perspective. A more
detailed review can be referred to [126, 132].

Figure 12: The fusion results of different methods with Chikusei image. The
color is composed of bands 70, 100, 36. An enlarged region is framed in green
and the corresponding residual image between the fused image and MS-GT is
framed in red.

###  VI-A Unmixing based methods

Hyperspectral unmixing (HU) [5, 101] assumes that the mixed class of an HS
image can be decomposed to the collection of constitute spectra (endmembers),
and their corresponding proportions (abundances). With LMM assumption, the
different endmembers do not interfere with each other [5]. By embedding the
LMM model into (23), we can obtain the following general unmixing based
approaches

$\displaystyle\min_{\mathbf{X},\mathbf{E},\mathbf{A}}\|\mathbf{Y}-\mathbf{X}\mathbf{R}\|_{\operatorname{F}}^{2}+\|\mathbf{Z}-\mathbf{G}\mathbf{X}\|_{\operatorname{F}}^{2},$
(24) $\displaystyle{\rm
s.t.}\;\;\mathbf{X}=\mathbf{E}\mathbf{A},\;\mathbf{E},\mathbf{A}\geq
0,\;\mathbf{1}_{R}^{\top}\mathbf{A}=\mathbf{1}_{MN},$

where $\mathbf{E}$, $\mathbf{A}$ are the endmember matrix and abundance
matrix, which are assumed to obey the non-negative and abundance sum-to-one
constraints. Generally, nonlinear unmixing models [5] can be also utilized for
the fusion task of HS-MS images. However, due to the generality of the LMM
model, we focus on the review of LMM based fusion approaches.

Eismann et al. proposed a maximum a posteriori estimation method to deduce the
cost objective function, and introduced a stochastic mixing model (MAP-SMM) to
embed LMM into the cost function [127]. MAP-SMM method tries to estimate the
prior probabilities for all the mixture classes, including the mean vectors
and covariance matrices of the endmember classes. The learned prior
probabilities are passed to the cost function to help the reconstruction of
the final HrHS image $\mathbf{X}$.

Yokoya et al., regarded (24) as the coupled NMF (CNMF) problem [128] and
introduced the multiplicative update rules to optimize (24). Firstly, CNMF
utilizes
$\|\mathbf{Y}-\mathbf{E}\mathbf{A}\mathbf{R}\|_{\operatorname{F}}^{2}$ as the
cost function to update $\mathbf{E}$ and $\mathbf{A}_{h}$ with the endmember
matrix $\mathbf{E}$ initialized by vertex component analysis (VCA). Here,
$\mathbf{A}_{h}=\mathbf{A}\mathbf{R}$ is the abundance matrix from HS images.
Secondly, by initializing $\mathbf{E}_{m}=\mathbf{G}\mathbf{E}$ which is the
endmember matrix from the MS image, CNMF again utilizes the multiplicative
update rules to update $\mathbf{E}$ from the cost function
$\|\mathbf{Z}-\mathbf{G}\mathbf{E}\mathbf{A}\|_{\operatorname{F}}^{2}$.
Finally, the HrHS image $\mathbf{X}$ is reconstructed from
$\mathbf{E}\mathbf{A}$. The following works [133, 134, 135] also utilize the
CNMF framework to fuse HS-MS images. Differently, [133, 135] introduced a non-
negative dictionary learning strategy, while [134] proposed the proximal
alternating linearized minimisation algorithm to update $\mathbf{E}$ and
$\mathbf{A}$.

On the basis of (24), Wang et al., further regularized $\mathbf{X}$ with non-
local low-rank Tucker decomposition [136]. The improved non-local Tucker
decomposition regularized CNMF model [136] was solved by the multi-block ADMM,
and achieved remarkable fusion results. It indicates that additional
regularizers on $\mathbf{X}$ can further improve the fusion accuracy. From
another side, it is necessary to make a trade-off between the complex models
with higher accuracy and the computation efficiency for real large scale HS-MS
image fusion task.

###  VI-B Orthogonal subspace based methods

Another common assumption in HS-MS fusion is that the spectral information of
$\mathbf{X}$ underlies a orthogonal subspace, whose dimension is much smaller
than the number of bands $B$ [129, 101], i.e.,
$\mathbf{X}=\mathbf{E}\mathbf{A}$ with $\mathbf{E}\in\mathbb{R}^{B\times k}$,
$\mathbf{A}\in\mathbb{R}^{k\times MN}$, and $k\ll B$. $\mathbf{E}^{\top}$ is
an orthogonal matrix with $\mathbf{E}^{\top}\mathbf{E}=\mathbf{I}_{k}$.
Therefore, the subspace based model is formulated as

$\displaystyle\min_{\mathbf{X},\mathbf{E},\mathbf{A}}\|\mathbf{Y}-\mathbf{X}\mathbf{R}\|_{\operatorname{F}}^{2}+\|\mathbf{Z}-\mathbf{G}\mathbf{X}\|_{\operatorname{F}}^{2},$
(25) $\displaystyle{\rm
s.t.}\;\;\mathbf{X}=\mathbf{E}\mathbf{A},\;\mathbf{E}^{\top}\mathbf{E}=\mathbf{I}_{k}.$

Although additional spectral subspace prior is exploited, the optimization of
(25) still faces several challenges. Firstly, if $k\gg b$, that’s to say, the
dimension number of the subspace is larger than the bands’ number of the MS
image, the optimization of (25) is the under-estimate problem. Therefore, to
ensure a reasonable solution, prior information of coefficient $\mathbf{A}$
need to be exploited. [129] pre-trains a dictionary to represent $\mathbf{A}$,
and updates $\mathbf{A}$ via sparse representation. Hyperspectral super-
resolution (HySure) in [130] assumes that $\mathbf{A}$ appears the spatial
smoothness structure and regularize $\mathbf{A}$ with band-by-band TV. [137]
translates the optimization of $\mathbf{A}$ to a Sylvester equation and
proposes a fast fusion method for (25) (FUSE).

Secondly, the optimization of orthogonal matrix $\mathbf{E}^{\top}$ is another
challenge due to the non-convex of (25). One appearing approach [129, 130,
137] is to pre-estimate $\mathbf{E}$ from $\mathbf{Y}$ in advance, and fix the
variable $\mathbf{E}$ during the optimization of (25). Specially, FUSE [137]
adopted principal component analysis (PCA), meanwhile, HySure utilized VCA to
extract $\mathbf{E}$ from $\mathbf{Y}$. Another strategy is to regularize the
update of $\mathbf{E}$ and $\mathbf{A}$ as a coupled matrix factorization
problem, and blind dictionary learning strategy is utilized to update
$\mathbf{E}$ [138]. A hybrid inexact block coordinate descent [139] is
introduced to exactly estimate $\mathbf{E}$.

###  VI-C Tensor based methods

The above subspace based methods utilize low-rank matrix decomposition to
exploit the low-rank property of the reshaped high-spatial resolution HS image
$\mathbf{X}$. However, the original HS image is a 3-D tensor, and therein, the
researchers introduce the tensor decomposition to simultaneously capture the
spatial and spectral low-rank property. The coupled sparse tensor
factorization (CSTF) approach [140] utilized Tucker decomposition, presented
as follows:

$\displaystyle\boldsymbol{\mathcal{X}}=\boldsymbol{\mathcal{O}}\times_{1}\mathbf{E}_{1}\times_{2}\mathbf{E}_{2}\times_{3}\mathbf{E}_{3},$
(26) $\displaystyle{\rm
subject\;to}\;\;\mathbf{E}_{i}^{\top}\mathbf{E}_{i}=\mathbf{I},\;\|\boldsymbol{\mathcal{O}}\|_{0}\leq\mathcal{C},$

to regularize the high-spatial resolution HS image $\boldsymbol{\mathcal{X}}$.
In (26), the core tensor $\boldsymbol{\mathcal{O}}$ is assumed to obey the
sparse property, and $\mathbf{E}_{i}$ is the orthogonal matrix of the $i$-th
dimension. Subsequently, CP decomposition [141], tensor train decomposition
[142], tensor ring decomposition [143, 131], and so on, are utilized to
regularize $\boldsymbol{\mathcal{X}}$. Furthermore, non-local LRTD is also
investigated for the fusion task [144, 145, 146].

It is worth noting that unmixing, orthogonal subspace, and tensor based
methods share the common idea that spectral space of $\mathbf{X}$ should lie
in the low-dimensional space. Unmixing based approaches interpret the low-rank
property as the endmembers and abundances, which are assumed to be non-
negative, meanwhile Orthogonal subspace and tensor based methods ignore the
non-negative restrict. Unmixing based approaches are interpretable from the
physical meaning, but suffering from the unstable convergence in the
optimization. Orthogonal subspace and tensor based methods lose physical
meaning, but can be optimized more elegantly.

Very recently, there are some preliminary works to perform the fusion task by
means of DL-based methods [147, 148, 149, 150, 151, 152, 153] and show the
effective and competitive fusion performance. A similar problem existed in
these methods is the model interpretability and rationality. Clearly
explaining the intrinsic meaning in each layer of deep networks would
contribute to better modeling the fusion task and further obtaining higher-
quality products.

TABLE IV: Quantitative comparison of different algorithms on the HS-MS image
fusion experiments. The best one is shown in bold.

Methods | RMSE | ERGAS | SA | SSIM
---|---|---|---|---
CNMF [128] | 6.404 | 0.715 | 4.89 | 0.8857
ICCV’15 [154] | 5.203 | 0.589 | 4.64 | 0.9139
HySure [130] | 8.537 | 0.812 | 9.45 | 0.8527
FUSE [137] | 8.652 | 0.869 | 9.51 | 0.8401
CSTF [140] | 8.32 | 0.841 | 8.34 | 0.8419
STEREO [141] | 9.4425 | 0.891 | 9.78 | 0.8231
NLSTF [144] | 8.254 | 0.819 | 8.36 | 0.8424

(a) Training for MML and CML (b) Testing for MML (c) Testing for CML

Figure 13: An illustration for the model’s training and testing in MML- and
CML-based classification tasks (take the bi-modality as an example). (a) They
share the same training process, i.e., two modalities are used for model
training. The main difference lies in the testing phase, (b) MML still needs
the input of two modalities, (c) while one modality is absent in CML.

###  VI-D Experimental Study

In this section, we select unmixing based methods:
CNMF171717http://naotoyokoya.com/Download.html [128],
ICCV’15181818https://github.com/lanha/SupResPALM [154]; subspace based
methods: HySure191919https://github.com/alfaiate [130],
FUSE202020https://github.com/qw245/BlindFuse [137]; tensor decomposition
regularized methods: STEREO212121https://sites.google.com/site/harikanats/
[141], CSTF [140]; finally non-local tensor decomposition regularized method:
NLSTF222222https://sites.google.com/view/renweidian/ [144] for the comparison
and analysis. We use evaluation indices, including the root mean square error
(RMSE), relative dimensional global error in synthesis (ERGAS) [155], MSA, and
SSIM [156] as evaluation criteria for the fusion results of different methods.

The selected dataset for the experiment is the Chikusei dataset obtained at
Chikusei, Ibaraki, Japan, on 29 July 2014 [126]. The selected high-spatial
resolution HS image is of size $448\times 448\times 128$, and the simulated
HS-MS images are of size $448\times 448\times 3$ and $14\times 14\times 128$,
respectively. Tab. (IV) presents the quantitative comparison results of
different algorithms on the HS-MS image fusion, meanwhile Fig. (12) presents
the visual illustration. From the results, it can be observed that even the HS
image is spatially degraded by 32 times, the fusion methods can efficiently
reconstruct the spatial details with the help of a 3 band MS image. On this
tested toy dataset, ICCV’15 performed the best. However, different datasets
need different kinds of regularizers. The fusion of HS-MS images for efficient
and large scale applications is still a challenge for further research.

###  VI-E Remaining Challenges

Subspace based non-convex methods for the fusion of HS-MS images have been
well developed. However, most remarkable results are achieved on the simulated
experiments. For the real applications with HS-MS images from two different
satellite sensors, there still remain several challenges.

  * •

Blind. Most fusion methods assume the linear spatial and spectral downsampling
from HR-HS image to HS-MS images. However, in real applications, the
degradation is complex and unknown in advance. how to blindly reconstruct the
HR-HS image is a challenge in future research.

  * •

Regularizer. We reviewed the subspace based fusion methods from unmixing,
orthogonal subspace, and tensor decomposition perspectives. Different
assumptions are suitable for the exploited for the different structure of the
HS image. How to mine the essence of HS images and develop efficient
regularizers for large scale processing still remains a challenge.

  * •

Evaluation. In the real cases, the enhanced HrHS images from HS and MS images
are not existed as the reference images in the real scenario. How to evaluate
the final enhanced HrHS images is also a key problem for the future fusion
approach development of HS-MS.

##  VII Cross-modality Learning for Large-scale Land Cover Mapping

With the ever-growing availability of diverse RS data sources from both
satellite and airborne sensors, multimodal data processing and analysis in RS
[157, 158] can provide potential possibilities to break the performance
bottleneck in many high-level applications, e.g., land cover classification.

HS data are featured by rich spectral information, which enables the high
discrimination ability for material recognition at a more accurate and fine
level. It should be noted, however, that the HS image coverage from space is
much narrow compared to MS imaging due to the limitations of imaging
principles and devices. That means the HS-dominated multimodal learning (MML)
fails to identify the materials on a large geographic coverage and even global
scale [159]. But fortunately, large-scale MS or synthetic aperture radar (SAR)
images are openly available from e.g., Sentinel-1, Sentinel-2, Landsat-8.
This, therefore, drives us to ponder over a problem: can HS images acquired
only in a limited area improve the land cover mapping performance using a
larger area covered by the MS or SAR images? This is a typical issue of cross-
modality learning (CML) from a ML’s point of view.

Take the bi-modality as an example, CML for simplicity refers to that training
a model using two modalities and one modality is absent in the testing phase,
or _vice versa_ (only one modality is available for training and bi-modality
for testing) [160]. Such a CML problem that exists widely in a variety of RS
tasks is more applicable to real-world cases. Fig. 13 illustrates the
differences between MML and CML in terms of training and testing process. The
core idea of CML is to find a new data space, where the information can be
exchanged effectively across different modalities. Thereupon, we formulate
this process in a general way as follows:

$\displaystyle\mathop{\min}_{\mathbf{X},\\{\mathbf{U}_{s}\\}_{s=1}^{m}}\sum_{s=1}^{m}\frac{1}{2}\lVert\mathbf{X}-\mathbf{U}_{s}\mathbf{Y}_{s}\rVert_{\operatorname{F}}^{2}\;\;{\rm
s.t.}\;\mathbf{X},\\{\mathbf{U}_{s}\\}_{s=1}^{m}\in\mathcal{C},$ (27)

where $m$ is the number of input modality. For simplicity, we only consider
the bi-modality case in this topic, i.e., $m=2$. According to different
learning strategies on modalities, CML can be roughly categorized into two
groups: manifold alignment (MA) and shared subspace learning (SSL). The
differences between the two types of approaches mainly lie in

  * •

MA learns the low-dimensional embedding by preserving the aligned manifold (or
graph) structure between different modalities. In the process of graph
construction, the similarities between samples (unsupervised MA) and indirect
label information (supervised or semi-supervised MA) are used. Despite the
competitive performance obtained by MA-based approaches for the CML task, the
discrimination ability of learned features remains limited due to the lack of
directly bridging low-dimensional features with label information.

  * •

SSL, as the name suggests, aims to find a latent shared subspace, where the
features of different modalities are linked via a manifold alignment
regularizer. Also, the learned features are further connected with label
information. The two steps are jointly optimized in a SSL model, tending to
yield more discriminative feature representations.

More specifically, we will briefly review and detail some representative
approaches belonging to the aforementioned two groups as follows.

###  VII-A Manifold Alignment based Approach

As the name suggests, MA is capable of aligning multiple modalities on
manifolds into a latent subspace, achieving a highly effective knowledge
transfer [161]. Due to the interactive learning ability, MA has a good fit for
large-scale RS image classification. In [162], the domain adaptation was
investigated to reduce the gap between the source and target domains of HS
data for land cover classification. By simultaneously considering labeled and
unlabeled samples, Tuia et al. [163] used semi-supervised MA (SSMA) techniques
[164] to align the multi-view RS images onto the manifold space by the
attempts to eliminate the effects of image variants caused by different views.
Matasci et al. [165] modified the classic transfer component analysis [166],
making it applicable to land cover classification of RS images. Moreover, a
kernelized MA approach presented in [167] projected the multimodal RS data to
a higher dimensional space and aligned them in a nonlinear way. Hu et al.
[168] deeply reviewed the semi-supervised MA methods with respect to the
fusion classification of HS and polarimetric SAR images. Based on the work in
[168], the same investigators made full use of topological data analysis and
designed a new graph structure for optical (e.g., HS) and SAR data fusion
[169].

Mathematically, the MA idea can be implemented by solving the following non-
convex model:

$\displaystyle\mathop{\min}_{\\{\mathbf{U}\\}_{s=1}^{m}}\frac{A+C}{B},$ (28)

where $A$, $B$, and $C$ are

$\displaystyle
A=\frac{1}{2}\sum_{p=1}^{m}\sum_{q=1}^{m}\sum_{i=1}^{n}\sum_{j=1}^{n}\lVert\mathbf{U}_{p}\mathbf{y}_{p}^{i}-\mathbf{U}_{q}\mathbf{y}_{q}^{j}\rVert_{2}^{2}\mathbf{W}_{sim}^{i,j},$
$\displaystyle
B=\frac{1}{2}\sum\limits_{p=1}^{m}\sum\limits_{q=1}^{m}\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n}\lVert\mathbf{U}_{p}\mathbf{y}_{p}^{i}-\mathbf{U}_{q}\mathbf{y}_{q}^{j}\rVert_{2}^{2}\mathbf{W}_{dis}^{i,j},$
$\displaystyle
C=\frac{1}{2}\sum_{t=1}^{m}\sum_{i=1}^{n}\sum_{j=1}^{n}\lVert\mathbf{U}_{t}\mathbf{y}_{t}^{i}-\mathbf{U}_{t}\mathbf{y}_{t}^{j}\rVert_{2}^{2}\mathbf{W}_{t}^{i,j}.$

By minimizing the problem (28), the $\\{\mathbf{U}\\}_{s=1}^{m}$ can be
estimated via generalized eigenvalues decomposition. We then have
$\mathbf{X}=\mathbf{U}_{s}\mathbf{Y}_{s}$. Three different graphs need to be
pre-computed in Eq. (28), including the similarity graph, i.e.,
$\mathbf{W}_{sim}$:

$\displaystyle\mathbf{W}_{sim}=\left[\begin{matrix}\mathbf{W}_{sim}^{1,1}&\mathbf{W}_{sim}^{1,2}&\cdots&\mathbf{W}_{sim}^{1,m}\\\
\mathbf{W}_{sim}^{2,1}&\mathbf{W}_{sim}^{2,2}&\cdots&\mathbf{W}_{sim}^{2,m}\\\
\vdots&\vdots&\ddots&\vdots\\\
\mathbf{W}_{sim}^{m,1}&\mathbf{W}_{sim}^{m,2}&\cdots&\mathbf{W}_{sim}^{m,m}\\\
\end{matrix}\right],$ (29)

the dissimilarity matrix, i.e., $\mathbf{W}_{dis}$:

$\displaystyle\mathbf{W}_{dis}=\left[\begin{matrix}\mathbf{W}_{dis}^{1,1}&\mathbf{W}_{dis}^{1,2}&\cdots&\mathbf{W}_{dis}^{1,m}\\\
\mathbf{W}_{dis}^{2,1}&\mathbf{W}_{dis}^{2,2}&\cdots&\mathbf{W}_{dis}^{2,m}\\\
\vdots&\vdots&\ddots&\vdots\\\
\mathbf{W}_{dis}^{m,1}&\mathbf{W}_{dis}^{m,2}&\cdots&\mathbf{W}_{dis}^{m,m}\\\
\end{matrix}\right],$ (30)

and the topology structure for each single modality obtained by $knn$ graph,
i.e., $\mathbf{W}_{t}$:

$\displaystyle\mathbf{W}_{t}=\left[\begin{matrix}\mathbf{W}_{t}^{1,1}&\mathbf{0}&\cdots&\mathbf{0}\\\
\mathbf{0}&\mathbf{W}_{d}^{2,2}&\cdots&\mathbf{0}\\\
\vdots&\vdots&\ddots&\vdots\\\
\mathbf{0}&\mathbf{0}&\cdots&\mathbf{W}_{t}^{m,m}\\\ \end{matrix}\right].$
(31)

In Eqs. (29)$-$(31), $\mathbf{W}_{sim}^{i,j}$, $\mathbf{W}_{dis}^{i,j}$, and
$\mathbf{W}_{t}^{i,j}$ are given, respectively, by

$\mathbf{W}_{sim}^{i,j}=\begin{cases}\begin{aligned} 1,\;\;&\text{if
$\mathbf{y}_{p}^{i}$ and $\mathbf{y}_{q}^{j}\in C_{k}$}\\\
0,\;\;&\text{otherwise,}\end{aligned}\end{cases}$
$\mathbf{W}_{dis}^{i,j}=\begin{cases}\begin{aligned} 1,\;\;&\text{if
$\mathbf{y}_{p}^{i}$ and $\mathbf{y}_{q}^{j}\notin C_{k}$ }\\\
0,\;\;&\text{otherwise,}\end{aligned}\end{cases}$
$\mathbf{W}_{t}^{i,j}=\begin{cases}\begin{aligned}
\exp\frac{\lVert\mathbf{y}^{i}-\mathbf{y}^{j}\rVert_{2}^{2}}{2\sigma^{2}},\;\;&\text{if
$\mathbf{y}_{p}^{i}\in\phi_{k}(\mathbf{y}_{q}^{j})$;}\\\
0,\;\;&\text{otherwise,}\end{aligned}\end{cases}$

where $\phi_{k}(\bullet)$ denotes the $k$ nearest neighbors of $\bullet$.

###  VII-B Shared Subspace Learning based Approach

Due to the lack of the direct relational modeling between the learned features
and label information, MA-based approaches fail to activate the connections
across modalities effectively [160], thereby yielding the relatively weak
transferability between different modalities, particularly heterogeneous data.
There have been some tentative works in recent years, providing potential
solutions to overcome the aforementioned challenges. For example, Hong et al.
[170] for the first time proposed a supervised CoSpace model to learn a latent
discriminative subspace from HS-MS correspondences for the CML-related
classification problem. Beyond it, the same authors [171] fully tapped the
potential of the CoSpace by learning the data-driven graph structure from both
labeled and unlabeled samples, yielding a learnable manifold alignment (LeMA)
approach. Moreover, [172] deeply investigated and analyzed different
regression techniques, i.e., $\ell_{2}$-norm ridge regression, $\ell_{1}$-norm
sparse regression, in CoSpace. In [173], a semi-supervised graph-induced
aligned learning (GiAL) was developed by jointly regressing labels and pseudo-
labels.

Accordingly, these methods can be generalized to be a unified model [170] to
address the CML’s problem in a regression-based fashion:

$\displaystyle\mathop{\min}_{\mathbf{P},\\{\mathbf{U}_{s}\\}_{s=1}^{m}}$
$\displaystyle\frac{1}{2}\lVert\mathbf{M}-\mathbf{P}\mathbf{U}_{s}\mathbf{Y}_{s}\rVert_{\operatorname{F}}^{2}+\Psi(\mathbf{P})+\Omega(\\{\mathbf{U}_{s}\\}_{s=1}^{m})$
(32) $\displaystyle{\rm
s.t.}\;\;\mathbf{U}_{s}\mathbf{U}_{s}^{\top}=\mathbf{I},\;s=1,\cdots,m,$

where $\\{\mathbf{U}_{s}\\}_{s=1}^{m}$ denote the projections linking to the
shared features for different modalities. To avoid the over-fitting of the
model and stabilize the learning process, $\mathbf{P}$ can be regularized by
the Frobenius-norm [170] or $\ell_{1,1}$-norm [172]:

$\displaystyle\Psi(\mathbf{P})=\lVert\mathbf{P}\rVert_{\operatorname{F}}^{2},\;\text{or}\;\lVert\mathbf{P}\rVert_{1,1},$
(33)

and $\Omega(\\{\mathbf{U}_{s}\\}_{s=1}^{m})$ is specified as a manifold
alignment term on the multimodal data, which is written as

$\displaystyle\Omega(\\{\mathbf{U}_{s}\\}_{s=1}^{m})=\operatorname{tr}(\mathbf{U}\mathbf{Y}\mathbf{L}\mathbf{Y}^{\top}\mathbf{U}^{\top}),$
(34)

where $\mathbf{U}=[\mathbf{U}_{1},\mathbf{U}_{2},\cdots,\mathbf{U}_{m}]$ and

$\displaystyle\mathbf{Y}=\left[\begin{matrix}\mathbf{Y}_{1}&\mathbf{0}&\cdots&\mathbf{0}\\\
\mathbf{0}&\mathbf{Y}_{2}&\cdots&\mathbf{0}\\\ \vdots&\vdots&\ddots&\vdots\\\
\mathbf{0}&\mathbf{0}&\cdots&\mathbf{Y}_{m}\\\ \end{matrix}\right].$

Similar to Fig. 7, $\mathbf{L}$ is a joint Laplacian matrix.

Using the general model in Eq. (32),

  * •

[170] considers the HS-MS correspondences that exist in an overlapped region
as the model input. The learned shared representations (e.g.,
$\mathbf{X}=\mathbf{U}_{s}\mathbf{Y}_{s}$) can be then used for classification
on a larger area, even though only MS data are available in the inference
phase;

  * •

Differently, [171] inputs not only the labeled HS-MS pairs but also unlabeled
MS data in large quantity. With the graph learning, i.e., the variable
$\mathbf{L}$ is to be learned from the data rather than fixed by a given RBF,
the unlabeled information can be made use of to find a better decision
boundary. According to the equivalent form of Eq. (34), we then have

$\displaystyle\operatorname{tr}(\mathbf{U}\mathbf{Y}\mathbf{L}\mathbf{Y}^{\top}\mathbf{U}^{\top})=\frac{1}{2}\operatorname{tr}(\mathbf{W}\mathbf{d})=\frac{1}{2}\lVert\mathbf{W}\odot\mathbf{d}\rVert_{1,1},$
(35)

where $\mathbf{d}_{i,j}=\lVert\mathbf{x}_{i}-\mathbf{x}_{j}\rVert_{2}^{2}$
denotes the pair-wise distance in Euclidean space. Using Eq. (35), the
resulting optimization problem with respect to the variable $\mathbf{W}$ is

$\displaystyle\frac{1}{2}$
$\displaystyle\lVert\mathbf{W}\odot\mathbf{d}\rVert_{1,1}$ (36)
$\displaystyle{\rm s.t.}\;\mathbf{W}=\mathbf{W}^{\top},\;\mathbf{W}_{i,j}\geq
0,\;\lVert\mathbf{W}\rVert_{1,1}=c.$

  * •

Inspired by the brain-like feedback mechanism presented in [89], a more
intelligent CML model was proposed [173]. With the joint use of labels and
pseudo-labels updated by the graph feedback in each iteration, more
representative features can be also learned (even if a certain modality is
absent, i.e., the CML case).

TABLE V: Quantitative comparison of SOTA algorithms related to the CML’s issue
in terms of OA, AA, and $\kappa$ using the NN classifier on the Houston2013
datasets. The best one is shown in bold.

Methods | OA (%) | AA (%) | $\kappa$
---|---|---|---
O-Baseline | 62.12 | 65.97 | 0.5889
USMA [85] | 65.54 | 68.81 | 0.6251
SMA [161] | 68.01 | 70.50 | 0.6520
SSMA [164] | 69.29 | 72.00 | 0.6659
CoSpace [170] | 69.38 | 71.69 | 0.6672
LeMA [171] | 73.42 | 74.76 | 0.7110
GiAL [173] | 80.66 | 81.31 | 0.7896

###  VII-C Experimental Study

We evaluate the performance of several SOTA algorithms related to the CML’s
issue both quantitatively and qualitatively. They are O-Baseline (i.e., using
original image features), unsupervised MA (USMA) [85], supervised MA
(SMA)232323https://sites.google.com/site/changwangnk/home/ma-html [161], SSMA
[164], CoSpace242424https://github.com/danfenghong/IEEE_TGRS_CoSpace [170],
LeMA252525https://github.com/danfenghong/ISPRS_LeMA [171], and GiAL [173].
Three common indices, e.g., OA, AA, and $\kappa$, are adopted to quantify the
classification performance using the SVM classifier on the Houston2013 HS-MS
datasets that have been widely used in many researches [170, 171, 172, 173].

Table V gives the quantitative comparison between the above-mentioned methods
for the CML-related classification, while Fig. 14 visualizes a region of
interest (ROI) of classification maps. By and large, the classification
accuracy of O-Baseline, i.e., only using MS data, is much lower than other
methods. By aligning multimodal data on manifolds, MA-based approaches perform
better than O-Baseline with the approximated increase of $3\%$ OA in USMA,
$6\%$ OA in SMA, and $7\%$ OA in SSMA. As expected, the classification
performance of SSL-based models, e.g., CoSpace, LeMA, and GiAL, is obviously
superior to that of MA-based ones. In particular, GiAL dramatically
outperforms other competitors, owing to the use of the brain-like feedback
mechanism and graph-driven pseudo-label learning. Visually, shared learning
methods tend to capture more robust spectral properties and achieve more
realistic classification results. As can be seen from Fig. 14, the shadow
region covered by clouds can be finely classified by CoSpace, LeMA, and GiAL,
while MA-based models fail to identify the materials well in the region.

Figure 14: ROI visualization of classification maps using different SOTA
methods related to the CML’s issue.

###  VII-D Remaining Challenges

CML has drawn growing interest from researchers in computer vision and ML, yet
it is rarely investigated in the RS community. In other words, CML is a
relatively emerging topic in RS, which means there are lots of difficulties
(or challenges) to be overcome. In detail,

  * •

Data Preparation. Since multimodal data are acquired by different contexts,
sensors, resolutions, etc., this inevitably poses a great challenge to data
collection and processing. For example, the errors caused by interpolation of
different resolutions, registration methods of geographical coordinates,
pixel-wise biases of different sensors, and uncertainties of image degradation
in the imaging process easily generate unregistered multimodal data to a great
extent.

  * •

Model Transferability. Due to different imaging mechanisms and principles, the
viscosity between pixels from the same modality is much stronger than from the
different modalities. This might lead to difficulties in fusing multimodal
information at a deep level, particularly heterogeneous data (e.g., HS and SAR
data), further limiting the model’s transferability.

  * •

Labeling. Unlike natural images or street view images that are relatively easy
and accurate to be labeled manually, labeling RS scenes (field trips needed)
is extremely expensive and time-consuming. Consequently, a limited number of
labeled samples are available for training and even worse, there are many
noisy labels in these samples. These problems will be noteworthy and to-be-
solved key points of the next generation interpretable AI models in the RS-
related CML task.

##  VIII Conclusion and Future Prospect

Characterized by the nearly continuous spectral profile that is capable of
sampling and representing the whole electromagnetic spectrum, HS images play
an important role in both promoting developments of new techniques and
accelerating the practical applications, not only limiting to the fields of
RS, geoscience, signal and image processing, numerical optimization and
modeling, ML, and AI. However, there still exist severe difficulties and
challenges that need to be carefully considered in the development and
application of HS RS techniques. One sign reveals that HS data analysis
methods dominated by expert systems have been unable to meet the demand of an
ever-growing volume of HS data whether in performance gain or in processing
efficiency. Another sign is that despite the currently unprecedented progress
made on computer vision, ML, and AI techniques, the model compatibility and
interpretability for HS RS applications remain limited.

Due to the SVs of HS data caused by various degradation mechanisms (e.g.,
environmental condition, atmospheric effects, spectral nonlinear mixing,
etc.), the redundancy of high-dimensional HS signals, and the complex
practical cases underlying in the HS products (e.g., low spatial resolution,
narrow imaging range, instrumental noises), convex models under ideal
circumstances usually fail to extract useful and diagnostic information from
HS images (especially those products that are corrupted seriously) and thereby
understand our environment. Considering that non-convex modeling is capable of
characterizing more complex real scenes and better providing the model
interpretability, in this article we present a comprehensive and technical
survey over five promising and representative research topics related to HS RS
with a focus on non-convex modeling, such as HS image restoration,
dimensionality reduction and classification, data fusion and enhancement,
spectral unmixing, and cross-modality learning. Among these topics, we review
the current state-of-the-art methods with illustrations, show the significance
and superiority of non-convex modeling to bridge the gap between HS RS and
interpretable AI, and point out remaining challenges and future research
directions.

It is well-known that the HS image processing and analysis chain is wide-
ranging. Apart from the five topics covered in this paper, we are not able to
detailedly report all the important and promising applications related to HS
RS missions. There are several very noteworthy and active fields that should
be paid more attention in future work, including target/change detection,
time-series analysis, multitemporal fusion/classification, physical parameter
inversion, image quality assessment, and various practical applications (e.g.,
precious farming, disaster management and response). Moreover, some crucial
steps for the HS image pre-processing algorithms are also missing, such as
atmospheric and geometric corrections, geographic coordinate registration,
etc. Furthermore, the methodologies summarized and reported in this article
mainly focus on the survey of shallow non-convex models. Undeniably, deep
models, e.g., DL-based methods, are capable of excavating deeper and intrinsic
properties of HS data. There is, therefore, room for improvement in the
development of more intelligent DL-related non-convex modeling with
application to HS RS. For example, embedding more physically meaningful priors
and devising advanced and novel deep unfolding [174] or unrolling [175]
strategies to closely integrate data-driven DL and theoretically-guaranteed
optimization technique is to open and interpret the so-called “black box” in
DL models.

Finally, we have to admit that non-convex modeling and optimization is a
powerful tool across multidisciplinary and the relevant studies along the
direction have been made tremendous progress theoretically and technically.
This provides the possibility of creating new methodologies and implementing
interpretable AI for various HS RS applications. In this paper, we attempt to
“intellectualize” these models by introducing more interpretable and
physically meaningful knowledge to meet the actual needs in a non-convex
modeling fashion. In other words, we hope that non-convex modeling can play
the role as a bridge to connect interpretable AI models and various research
topics in HS RS. Our efforts in this paper are made to foster curiosity and
create a good starting point for post-graduate, Ph.D. students, and senior
researchers working in the HS-related fields, thereby further looking for new
and advanced research directions in the interdisciplinary involving signal and
image processing, ML, AI, and RS.

## Acknowledgement

The authors would like to thank Prof. D. Landgrebe from Purdue University for
providing the AVIRIS Indian Pines data, Prof. P. Gamba from the University of
Pavia for providing the ROSIS-3 Pavia University and Centre data, the
Hyperspectral Image Analysis group at the University of Houston for providing
the CASI University of Houston dataset used in the IEEE GRSS DFC2013 and
DFC2018, and the Hyperspectral Digital Imagery Collection Experiment (HYDICE)
for sharing the urban dataset free of charge.

This work from the D. Hong and X. Zhu sides is jointly supported by the German
Research Foundation (DFG) under grant ZH 498/7-2, by the Helmholtz Association
through the framework of Helmholtz Artificial Intelligence (HAICU) - Local
Unit “Munich Unit @Aeronautics, Space and Transport (MASTr)” and Helmholtz
Excellent Professorship “Data Science in Earth Observation - Big Data Fusion
for Urban Research”, by the German Federal Ministry of Education and Research
(BMBF) in the framework of the international AI future lab “AI4EO – Artificial
Intelligence for Earth Observation: Reasoning, Uncertainties, Ethics and
Beyond”. This work from the L. Gao side is supported by the National Natural
Science Foundation of China under Grant 42030111 and Grant 41722108. This work
from the W. He and N. Yokoya sides is supported by the Japan Society for the
Promotion of Science under KAKENHI 19K20308 and KAKENHI 18K18067. This work
from the J. Chanussot side has been partially supported by MIAI@Grenoble
Alpes, (ANR-19-P3IA-0003) and by the AXA Research Fund.

The corresponding authors of this paper are Dr. Wei He and Prof. Lianru Gao.

## References

  * [1] A. Goetz, G. Vane, J. Solomon, and B. Rock, “Imaging spectrometry for earth remote sensing,” Science, vol. 228, no. 4704, pp. 1147–1153, 1985.
  * [2] L. Tsang, J. A. Kong, and R. T. Shin, “Theory of microwave remote sensing,” 1985\.
  * [3] W. Turner, S. Spector, N. Gardiner, M. Fladeland, E. Sterling, and M. Steininger, “Remote sensing for biodiversity science and conservation,” Trends Ecol. Evol., vol. 18, no. 6, pp. 306–314, 2003.
  * [4] D. Hong, W. Liu, J. Su, Z. Pan, and G. Wang, “A novel hierarchical approach for multispectral palmprint recognition,” Neurocomputing, vol. 151, pp. 511–521, 2015.
  * [5] J. Bioucas-Dias, A. Plaza, G. Camps-Valls, P. Scheunders, N. Nasrabadi, and J. Chanussot, “Hyperspectral remote sensing data analysis and future challenges,” IEEE Geosci. Remote Sens. Mag., vol. 1, no. 2, pp. 6–36, 2013\.
  * [6] S. Boyd, S. P. Boyd, and L. Vandenberghe, Convex optimization.  Cambridge university press, 2004.
  * [7] I. Ekeland, “Nonconvex minimization problems,” Bull. Amer. Math. Soc., vol. 1, no. 3, pp. 443–474, 1979.
  * [8] M. Mohri, A. Rostamizadeh, and A. Talwalkar, Foundations of machine learning.  MIT press, 2018.
  * [9] N. J. Nilsson, Principles of artificial intelligence.  Morgan Kaufmann, 2014.
  * [10] V. Chvatal, V. Chvatal, et al., Linear programming.  Macmillan, 1983.
  * [11] M. Frank, P. Wolfe, et al., “An algorithm for quadratic programming,” Nav. Res. Logist. Q., vol. 3, no. 1-2, pp. 95–110, 1956.
  * [12] M. S. Lobo, L. Vandenberghe, S. Boyd, and H. Lebret, “Applications of second-order cone programming,” Linear Algebra Appl., vol. 284, no. 1-3, pp. 193–228, 1998.
  * [13] P. Jain, P. Kar, et al., “Non-convex optimization for machine learning,” Found. Trends® Mach. Learn., vol. 10, no. 3-4, pp. 142–363, 2017.
  * [14] S. Boyd, N. Parikh, and E. Chu, Distributed optimization and statistical learning via the alternating direction method of multipliers.  Now Publishers Inc, 2011.
  * [15] J. Martin-Herrero, “Anisotropic diffusion in the hypercube,” IEEE Trans. Geosci. Remote Sens., vol. 45, no. 5, pp. 1386–1398, 2007.
  * [16] B. Datt, T. R. McVicar, T. G. Van Niel, D. L. Jupp, and J. S. Pearlman, “Preprocessing eo-1 hyperion hyperspectral data to support the application of agricultural indexes,” IEEE Trans. Geosci. Remote Sens., vol. 41, no. 6, pp. 1246–1259, 2003.
  * [17] L. Gómez-Chova, J. Amorós, G. Camps-Valls, J. D. Martin, J. Calpe, L. Alonso, L. Guanter, J. C. Fortea, and J. Moreno, “Cloud detection for chris/proba hyperspectral images,” in Remote Sensing of Clouds and the Atmosphere X, vol. 5979, p. 59791Q, International Society for Optics and Photonics, 2005.
  * [18] N. Acito, M. Diani, and G. Corsini, “Signal-dependent noise modeling and model parameter estimation in hyperspectral images,” IEEE Trans. Geosci. Remote Sens., vol. 49, no. 8, pp. 2957–2971, 2011.
  * [19] C.-I. Chang and Q. Du, “Interference and noise-adjusted principal components analysis,” IEEE Trans. Geosci. Remote Sens., vol. 37, pp. 2387–2396, Sep. 1999.
  * [20] Q. Yuan, L. Zhang, and H. Shen, “Hyperspectral image denoising employing a spectral–spatial adaptive total variation model,” IEEE Trans. Geosci. Remote Sens., vol. 50, no. 10, pp. 3660–3677, 2012.
  * [21] H. Zhang, W. He, L. Zhang, H. Shen, and Q. Yuan, “Hyperspectral image restoration using low-rank matrix recovery,” IEEE Trans. Geosci. Remote Sens., vol. 52, no. 8, pp. 4729–4743, 2013.
  * [22] W. He, H. Zhang, L. Zhang, and H. Shen, “Total-variation-regularized low-rank matrix factorization for hyperspectral image restoration,” IEEE Trans. Geosci. Remote Sens., vol. 54, no. 1, pp. 178–188, 2015.
  * [23] A. A. Green, M. Berman, P. Switzer, and M. D. Craig, “A transformation for ordering multispectral data in terms of image quality with implications for noise removal,” IEEE Trans. Geosci. Remote Sens., vol. 26, no. 1, pp. 65–74, 1988.
  * [24] B. Rasti, J. R. Sveinsson, and M. O. Ulfarsson, “Wavelet-based sparse reduced-rank regression for hyperspectral image restoration,” IEEE Trans. Geosci. Remote Sens., vol. 52, no. 10, pp. 6688–6698, 2014.
  * [25] Y. Qian and M. Ye, “Hyperspectral imagery restoration using nonlocal spectral-spatial structured sparse representation with noise estimation,” IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens., vol. 6, no. 2, pp. 499–515, 2012.
  * [26] J. Li, Q. Yuan, H. Shen, and L. Zhang, “Noise removal from hyperspectral image with joint spectral–spatial distributed sparse representation,” IEEE Trans. Geosci. Remote Sens., vol. 54, no. 9, pp. 5425–5439, 2016.
  * [27] Z. Wu, Q. Wang, J. Jin, and Y. Shen, “Structure tensor total variation-regularized weighted nuclear norm minimization for hyperspectral image mixed denoising,” Signal Process., vol. 131, pp. 202–219, 2017.
  * [28] Y. Peng, D. Meng, Z. Xu, C. Gao, Y. Yang, and B. Zhang, “Decomposable nonlocal tensor dictionary learning for multispectral image denoising,” in Proc. CVPR, pp. 2949–2956, 2014.
  * [29] W. He, Q. Yao, C. Li, N. Yokoya, Q. Zhao, H. Zhang, and L. Zhang, “Non-local meets global: An integrated paradigm for hyperspectral image restoration,” IEEE Trans. Pattern Anal. Mach. Intell., pp. 1–1, 10.1109/TPAMI.2020.3027563, 2020.
  * [30] W. He, H. Zhang, L. Zhang, and H. Shen, “Hyperspectral image denoising via noise-adjusted iterative low-rank matrix approximation,” IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens., vol. 8, no. 6, pp. 3050–3061, 2015.
  * [31] Y. Xie, Y. Qu, D. Tao, W. Wu, Q. Yuan, W. Zhang, et al., “Hyperspectral image restoration via iteratively regularized weighted schatten p-norm minimization.,” IEEE Trans. Geosci. Remote Sens., vol. 54, pp. 4642–4659, Aug. 2016.