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

HAL
open science

# Hyperspectral Restoration and Fusion With Multispectral Imagery via Low-Rank Tensor-Approximation

Na Liu, Lu Li, Wei Li, Ran Tao, James Fowler, Jocelyn Chanussot

► To cite this version:

Na Liu, Lu Li, Wei Li, Ran Tao, James Fowler, et al. Hyperspectral Restoration and Fusion With Multispectral Imagery via Low-Rank Tensor-Approximation. IEEE Transactions on Geoscience and Remote Sensing, 2021, 59 (9), pp.7817-7830. 10.1109/TGRS.2020.3049014. hal-03429690

HAL Id: hal-03429690

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

Submitted on 25 Sep 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Hyperspectral Restoration and Fusion With Multispectral Imagery via Low-Rank Tensor-Approximation

Na Liu<sup>D</sup>, Student Member, IEEE, Lu Li<sup>D</sup>, Wei Li<sup>D</sup>, Senior Member, IEEE, Ran Tao<sup>D</sup>, Senior Member, IEEE, James E. Fowler<sup>D</sup>, Fellow, IEEE, and Jocelyn Chanussot<sup>D</sup>, Fellow, IEEE

**Abstract**—Tensor-based fusion that couples the high spatial resolution of a multispectral image (MSI) to the high spectral resolution of a hyperspectral image (HSI) is considered. The fusion problem is first formulated mathematically as a convex optimization of a tensor trace norm imposing low-rank spatially as well as spectrally, with an alternating-directions optimization featuring linearization providing the solution. Although prior tensor-based fusion approaches typically resort to tensor decomposition, the proposed algorithm exploits ideas from the field of tensor completion to directly impose a low-rank property spatially and spectrally while avoiding the computationally complex patch clustering and dictionary learning common to competing fusion techniques. Additionally, small modifications to the basic optimization permit a fusion process robust to missing hyperspectral values such as those that can result from dead stripes in real hyperspectral sensors. The experimental evaluations on both synthetic imagery as well as real imagery demonstrate that the resulting low-rank tensor-approximation (LRTA) fusion algorithm preserves both spatial details and texture, yielding significantly improved image quality when compared to other state-of-the-art fusion methods as well as effective restoration under conditions of missing stripes within the HSI.

**Index Terms**—Data fusion, hyperspectral imagery (HSI), low-rank tensor.

Manuscript received October 15, 2020; revised December 20, 2020; accepted December 29, 2020. This paper is supported by the Beijing Natural Science Foundation (Grant no. JQ20021), National Natural Science Foundation of China (Grant no. 61922013). (Na Liu and Lu Li contributed equally to this work.) (Corresponding author: Wei Li.)

Na Liu is with the College of Information Science and Technology, Beijing University of Chemical Technology, Beijing 100029, China, and also with the Université Grenoble Alpes, 38000 Grenoble, France (e-mail: liuna@ieee.org).

Lu Li is with the School of Automation, Beijing Information Science and Technology University, Beijing 100192, China (e-mail: 1141017448@qq.com).

Wei Li and Ran Tao are with the School of Information and Electronics, Beijing Institute of Technology, Beijing 100081, China, and also with the Beijing Key Laboratory of Fractional Signals and Systems, Beijing 100081, China (e-mail: liwei089@ieee.org; rantao@bit.edu.cn).

James E. Fowler is with the Department of Electrical and Computer Engineering, Mississippi State University, Starkville, MS 39762 USA (e-mail: fowler@ece.msstate.edu).

Jocelyn Chanussot is with the LJK, CNRS, INRIA, Grenoble INP, University of Grenoble Alpes, 38000 Grenoble, France, and also with the Aerospace Information Research Institute, Chinese Academy of Sciences, Beijing 100094, China (e-mail: jocelyn.chanussot@gipsa-lab.grenoble-inp.fr).

Color versions of one or more figures in this article are available at https://doi.org/10.1109/TGRS.2020.3049014.

Digital Object Identifier 10.1109/TGRS.2020.3049014

## I. INTRODUCTION

DUE to an inevitable tradeoff between spatial and spectral resolution in optical remote-sensing systems, hyperspectral imagery (HSI) is acquired with very high spectral resolution, but with relatively poor spatial resolution, while multispectral imagery (MSI) typically has higher spatial resolution at the expense of poorer spectral resolution. Consequently, using MSI to enhance the spatial resolution of HSI is of great interest in numerous applications [1]–[4]. For example, Yokaya et al. [3] surveyed a number of recent efforts toward fusing HSI with MSI in order to generate an image with both high spectral and spatial resolution, which we denote as HS$^2$I.

Many such hyperspectral and multispectral (HS-MS) fusion methods entail casting the HSI and MSI image cubes as 2-D matrices by raster-scanning the image spatially, that is, an image cube of spatial size of $M \times N$ pixels, wherein each pixel has $D$ spectral bands, is represented as $D \times MN$ matrix $X = [x_1 \dots x_{MN}]$, where pixel vector $x_i \in \mathbb{R}^D$, $1 \le i \le MN$. However, such a matricization of the 3-D image cube tends to obfuscate both the spatial texture structure as well as the spatial-spectral dependence within the image cube [5], [6]. This necessitates the addition of explicit strategies to exploit joint spatial-spectral characteristics—such as various regularizers (e.g., [7])—which no doubt increases the complexity of algorithms.

As an alternative to matricization, there has been increasing interest in treating an HSI or MSI image cube as being a multidimensional array representation of an order-3 tensor, such that tensor theory and processes can be brought to bear on the HS-MS fusion task. A popular strategy is to deploy a tensor decomposition, such as the well-known Tucker decomposition, e.g. [5], [6], [8]. Specifically, the target HS$^2$I is formed from a core tensor multiplied by its three factor matrices, where two factor matrices are learned from the MSI in the two spatial modes, while the third is learned from the HSI in the spectral mode. While such tensor-decomposition approaches to HS-MS fusion are conceptually straightforward and can yield effective results, they tend to be fraught with several difficulties, including the well-known nonuniqueness of the Tucker decomposition (making it less than straightforward how to formulate the decomposition), the high computational
complexity entailed in learning the factor matrices (typically via patch-based dictionary-learning schemes), and the fact that finding a factor matrix which is optimal with respect to some criterion does not necessarily yield a correspondingly optimal HS²I after multiplication with the core tensor (due to the lack of an equivalent to Parseval's theorem for general tensor decompositions<fnref n="1" />). Accordingly, additional spatial-spectral regularizers have been coupled to tensor decompositions to further incorporate spatial and spectral information into the fused HS²I (e.g., [10], [11]. Alternatively, patch clustering has also been used to drive tensor-based fusion directly without resorting to tensor decompositions [12]. Finally, both coupled-tensor-ring representation and factorization models have also recently been introduced for HSI-fusion problems [13], [14].

However, it has been established [15]–[17] that HSI possesses not only global correlation across the spectrum but also nonlocal self-similarity over space, guaranteeing low rank in both the spatial and spectral directions. Specifically, exploiting low-rank aspects of HSI has shown significant potential in a myriad of HSI applications, including hyperspectral dimensionality reduction, anomaly detection, and noise reduction (e.g., [16]–[20]). Given that computing the rank of a tensor is NP-hard [21], low-rank methods commonly employ a proxy for tensor rank, for example, Liu et al. [22] adopted a tensor trace norm which is demonstrated to be effective at capturing geometric structure and texture for image inpainting via tensor completion.

Apart from the issue of low spatial resolution mentioned above, another challenge for HSI recovery is pixel degradation caused by optical aberrations and misalignments in imaging systems, especially stripe noise [23], [24]. Such striping artifacts severely degrade the quality of HSI and ultimately influence the HS-MS fusion result (i.e., significant noise can remain in the fused image). As described in [3], a straightforward approach to reduce noise in the fused image is through postprocessing with denoising. Yet, although denoising after fusion can improve image quality to a certain extent, it can also introduce new errors or fail to remove all the noise. Consequently, noise-robust HS-MS fusion is considered here in order to improve spatial resolution while simultaneously restoring HSI degradations. While we focus specifically on stripe noise, restoration of other types of artifacts (e.g., missing or dead patches, isolated pixels, cloud cover, etc.) would also be feasible.

To this end, in this article, we recast the tensor-trace-norm formulation used in [22] for tensor completion to drive HS-MS fusion, producing an HS²I cube that couples the spatial geometry and texture of the constituent MSI with the spectral detail from the HSI. Our resulting low-rank tensor-approximation (LRTA) algorithm imposes a low-rank tensor trace norm directly on the target HS²I cube, thereby circumventing the difficulties outlined above surrounding tensor decompositions. We also avoid the computationally complex patch clustering and dictionary learning common to existing techniques by

devising an augmented Lagrangian to optimize the tensor trace norm directly, effectively adapting ideas from the field of tensor completion to the HS-MS fusion problem. This has the added benefit that the proposed LRTA fusion algorithm can be readily modified to be able to complete the fusion task while simultaneously recovering missing HSI values due to stripe noise. We believe that the resulting LRTA algorithm represents the first approach to appear in the literature that applies such tensor-completion concepts to the HS-MS fusion problem.

In the remainder of the article, we describe our proposed LRTA approach to HS-MS fusion and experimentally evaluate its effectiveness. After introducing some basic mathematical preliminaries in Section II, we provide the details on the LRTA HS-MS fusion algorithm in Section III. Section IV follows with experimental results that compare the performance of the proposed algorithm against that of several other state-of-the-art fusion approaches. Finally, Section V closes the article with some concluding remarks.

## II. MATHEMATICAL PRELIMINARIES

A real-valued order-n tensor is denoted as $\mathcal{X} \in \mathbb{R}^{I_1 \times \cdots \times I_k \times \cdots \times I_n}$, with each element being $x_{i_1, \dots, i_k, \dots, i_n}$, where $I_k$ is the dimension of mode k, $1 \le k \le n$. A tensor may be unfolded into a matrix along a single mode; the “unfold” operation along mode k is defined as $\text{unfold}_k(\mathcal{X}) = \mathcal{X}_{(k)} \in \mathbb{R}^{I_k \times (I_1 \times \cdots I_{k-1} I_{k+1} \cdots I_n)}$. Correspondingly, the “fold” operation reverses the unfolding, $\text{fold}_k(\mathcal{X}_{(k)}) = \mathcal{X}$.

As is commonly done with matrices, we use a trace norm as a convex proxy for rank; for tensors, we adopt the trace norm from [22],

$$ \|\mathcal{X}\|_* = \sum_{k=1}^{n} \alpha_k \|\mathcal{X}_{(k)}\|_* \quad (1) $$

where $\|\cdot\|$ is the matrix spectral norm and $\alpha_k \ge 0$ are weights such that $\sum_{k=1}^n \alpha_k = 1$. We also define the Frobenius norm of a tensor as $\|\mathcal{X}\|_F = (\sum_{i_1, \dots, i_n} |x_{i_1, \dots, i_n}|^2)^{1/2}$ such that $\|\mathcal{X}\|_F = \|\mathcal{X}_{(k)}\|_F$, $1 \le k \le n$.

The multiplication of a tensor $\mathcal{X}$ with a matrix $\mathbf{A} \in \mathbb{R}^{J_k \times I_k}$ on mode k is defined by $\mathcal{Y} = \mathcal{X} \times_k \mathbf{A}$, where $\mathcal{Y} \in \mathbb{R}^{I_1 \times \cdots \times I_{k-1} \times I_k \times I_{k+1} \cdots \times I_n}$ and $y_{i_1, \dots, i_{k-1}, j_k, i_{k+1}, \dots, i_n} = \sum_{i_k} I_k = x_{i_1, \dots, i_k, \dots, i_n} a_{j_k, i_k}$. Alternatively, after an unfold operation along mode k, the multiplication can be defined as $\mathcal{Y}_{(k)} = \mathbf{A}\mathcal{X}_{(k)}$. The element-wise product operator of tensor $\mathcal{X}$ with tensor $\mathcal{W}$ of the same size is denoted as $\mathcal{X} \odot \mathcal{W}$. Finally, the interpolation operation of tensor $\mathcal{X}$ by ratio R is defined by $\mathcal{Z} = \text{interpolation}(\mathcal{X}, R, \text{"method"})$, which returns a tensor in which the first two orders are resized, that is, $\mathcal{Z} \in \mathbb{R}^{RI_1 \times RI_2 \times \cdots \times I_k \times \cdots \times I_n}$. Here, "method" indicates the interpolation process used, such as nearest-neighbor interpolation, bilinear interpolation, bicubic interpolation, and so on.

## III. PROPOSED LRTA FUSION ALGORITHM

The core principle in estimating an HS²I is to extract spectral and spatial information from the HSI and MSI, respectively. To do so, the spatial and spectral degradation from the HS²I to the observed HSI and MSI, respectively, are modeled. The HS-MS fusion task is then formulated as

<footnote n="1">¹This is effectively a consequence of the fact that tensor decompositions into Parseval frames are known for only certain symmetric tensors (e.g., [9]); the HS²I will be, in general, nonsymmetric.</footnote>
<figure><img src="image_3.png" /><figcaption>Fig. 1. Proposed LRTA HS-MS fusion framework. The point-spread function (PSF) and spectral response function (SRF) are simplified here for visualization purposes. On the left is the spatial and spectral degradation of the HSI and MSI; on the right is the proposed fusion algorithm wherein spectral correlation and spatial nonlocal similarity are integrated into the desired HS²I via LRTA.</figcaption></figure>

an LRTA problem which is solved using the recent linearized alternating direction method with adaptive penalty (LADMAP) [25]–[27], a variant of alternating direction method of multipliers (ADMM) useful for when constituent subproblems do not have closed-form solutions.

Below, we introduce the specifics of the proposed LRTA HS-MS fusion algorithm. We start in Section III-A with a presentation of the mathematical model of the HSI/MSI sensor observation process along with a mathematical formulation of the HS²I fusion problem in terms of a tensor trace norm. The proposed LRTA algorithm is described in Section III-B to solve this fusion problem. Finally, the LRTA algorithm is extended in Section III-C to a more realistic scenario wherein the HSI image cube is degraded by missing values, such as when the sensor induces striping artifacts. In this latter case, LRTA learns the fused HS²I while effectively recovering the missing HSI values.

### A. Problem Formulation

As is common in HS-MS fusion, we model both the observed HSI and MSI as degradations of the real ground target, wherein the HSI is the result of a spatial-degradation process applied to the HS²I, while the MSI results from spectral degradation; both degradation processes are a consequence of the inevitable tradeoff between spatial and spectral resolution in real optical remote-sensing systems. As illustrated in Fig. 1, the low-spatial-resolution HSI (denoted as $\mathcal{L} \in \mathbb{R}^{M' \times N' \times D}$) and the high-spatial-resolution MSI (denoted as $\mathcal{H} \in \mathbb{R}^{M \times N \times D'}$) are regarded as spatially and spectrally degraded versions, respectively, of the unknown HS²I, $\mathcal{X} \in \mathbb{R}^{M \times N \times D}$. Hence, the HSI and MSI are modeled as

$$ \begin{aligned} \mathcal{X}_{(3)}\mathbf{D} &= \mathcal{L}_{(3)} \\ \mathbf{S}\mathcal{X}_{(3)} &= \mathcal{H}_{(3)} \end{aligned} \quad (2) $$

where $\mathbf{S} \in \mathbb{R}^{D' \times D}$ denotes a matrix spanned by a spectral response function (SRF) that constitutes a spectral-downsampling process, while $\mathbf{D} \in \mathbb{R}^{M' \times M' \times N'}$ is a matrix generated from an $R \times R$ point-spread kernel (PSK) $\Psi$, $R = M/M' = N/N'$. Specifically, $\mathbf{D}$ encapsulates the hyperspectral sensor's spatial point-spread function (PSF) describing the sensor's spatial degradation process (represented discretely by $\Psi$) coupled with the subsampling entailed by the imaging process.

The underlying low-rank nature of images has been widely exploited for various image- and video-processing applications [28]. We expect that the real ground-target HS²I is characterized by strong correlations in spectrally adjacent bands as well as high similarity between spatially neighboring pixels; both of these manifest mathematically as low-rankness both spatially and spectrally. Thus, we formulate the HS-MS fusion problem as an estimation of the desired $\mathcal{X}$ with low rank under constraints of the two sensor observation models. Specifically, we cast the HS-MS fusion problem as follows:

$$ \begin{aligned} & \min_{\mathcal{X}'} \sum_{k=1}^{3} \alpha_k \|\mathcal{X}_{(k)}\|_{*} \\ & \text{s.t.} \quad \mathcal{X}_{(3)}\mathbf{D} = \mathcal{L}_{(3)} \text{ and } \mathbf{S}\mathcal{X}_{(3)} = \mathcal{H}_{(3)} \end{aligned} \quad (3) $$

where parameters $\alpha_k$ controls the intensity of low-rank constraints on each mode of the expected fusion result $\mathcal{X}$ (setting these parameters will be considered later in Section IV-C). As the trace norm is the tightest convex envelope for the rank, (3) is a convex optimization which we solve below. We note that the HS-MS fusion problem (3) can be regarded, in essence, as the completion problem considered in [22]. As discussed in [22], tensor completion builds up the relationship between the known and unknown elements to estimate the latter. The construction of the relationship facilitates the estimation method's ability to capture the information from both neighbor entries as well as global structure. Inspired by
this core principle, (3) is expected to estimate the unknown $\mathcal{X}$ by building up its relationship to the known $\mathcal{L}$ and $\mathcal{H}$ through observation models.

## B. Proposed LRTA HS-MS Fusion Algorithm

The problem in (3) is difficult to solve directly via differential calculus because elements in the tensor are interdependent. Consequently, we follow the LADMAP paradigm in which we introduce, for $k \in \{1, 2, 3\}$, auxiliary variables $\mathbf{M}_k$ for $\mathcal{X}_{(k)}$. Additionally, we have Lagrange multipliers $\mathbf{Y}'_1, \mathbf{Y}'_2$, and $\mathbf{Y}_k$, where, again, $k \in \{1, 2, 3\}$, as well as penalty parameters $\mu, \gamma$, and $\beta$. This renders (3) separable such that it can be solved via LADMAP. Specifically, we form the augmented Lagrangian as follows:

$$ 
\begin{aligned}
L_{\mu, \beta, \gamma}(\mathcal{X}, \mathbf{M}_1, \mathbf{M}_2, \mathbf{M}_3, \mathbf{Y}_1, \mathbf{Y}_2, \mathbf{Y}_3, \mathbf{Y}'_1, \mathbf{Y}'_2) \\
= \sum_{k=1}^{3} \alpha_k \|\mathbf{M}_k\|_{*} + \sum_{k=1}^{3} \langle \mathbf{Y}_k, \mathcal{X}_{(k)} - \mathbf{M}_k \rangle \\
+ \langle \mathbf{Y}'_1, \mathbf{M}_3 \mathbf{D} - \mathcal{L}_{(3)} \rangle + \langle \mathbf{Y}'_2, \mathbf{SM}_3 - \mathcal{H}_{(3)} \rangle \\
+ \frac{\mu}{2} \sum_{k=1}^{3} \|\mathcal{X}_{(k)} - \mathbf{M}_k\|_{\mathrm{F}}^2 + \frac{\beta}{2} \|\mathbf{M}_3 \mathbf{D} - \mathcal{L}_{(3)}\|_{\mathrm{F}}^2 \\
+ \frac{\gamma}{2} \|\mathbf{SM}_3 - \mathcal{H}_{(3)}\|_{\mathrm{F}}^2
\end{aligned}
\quad (4) $$

which we can rewrite as

$$ 
\begin{aligned}
L_{\mu, \beta, \gamma}(\mathcal{X}, \mathbf{M}_1, \mathbf{M}_2, \mathbf{M}_3, \mathbf{Y}_1, \mathbf{Y}_2, \mathbf{Y}_3, \mathbf{Y}'_1, \mathbf{Y}'_2) \\
= \sum_{k=1}^{3} \alpha_k \|\mathbf{M}_k\|_{*} - \frac{1}{2\mu} \sum_{k=1}^{3} \|\mathbf{Y}_k\|_{\mathrm{F}}^2 - \frac{1}{2\beta} \|\mathbf{Y}'_1\|_{\mathrm{F}}^2 - \frac{1}{2\gamma} \|\mathbf{Y}'_2\|_{\mathrm{F}}^2 \\
+ \frac{\mu}{2} \sum_{k=1}^{3} \|\mathcal{X}_{(k)} - \mathbf{M}_k + \frac{\mathbf{Y}_k}{\mu}\|_{\mathrm{F}}^2 + \frac{\beta}{2} \|\mathbf{M}_3 \mathbf{D} - \mathcal{L}_{(3)} + \frac{\mathbf{Y}'_1}{\beta}\|_{\mathrm{F}}^2 \\
+ \frac{\gamma}{2} \|\mathbf{SM}_3 - \mathcal{H}_{(3)} + \frac{\mathbf{Y}'_2}{\gamma}\|_{\mathrm{F}}^2.
\end{aligned}
\quad (5) $$

Optimization problem (3) is then solved iteratively via (5) with the following subproblems:

1) Update $\mathbf{M}_k$ for $k = 1$ and $2$ via

$$ \mathbf{M}_{k,t+1} = \arg \min_{\mathbf{M}_k} \alpha_k \|\mathbf{M}_k\|_{*} + \frac{\mu}{2} \left\| \mathcal{X}_{(k,t)} - \mathbf{M}_k + \frac{\mathbf{Y}_{k,t}}{\mu} \right\|_{\mathrm{F}}^2. \quad (6) $$

Problems of the form of (6) have a closed-form solution via the singular-value thresholding operator [29]

$$ D_{\tau}(\mathbf{X}) = U S_{\tau} [\mathbf{\Sigma}] \mathbf{V}^{H} \quad (7) $$

where $\mathbf{X} = \mathbf{U}\mathbf{\Sigma}\mathbf{V}^{H}$ is a singular-value decomposition of matrix $\mathbf{X}$, and the matrix-valued shrinkage operator $S_{\tau}[\mathbf{X}]$ is the scalar shrinkage operator $S_{\tau}[x] = \operatorname{sgn}(x) \max\{|x| - \tau, 0\}$ applied to each element of $\mathbf{X}$. Therefore, $\mathbf{M}_k$ is updated as follows:

$$ \mathbf{M}_{k,t+1} = D_{\frac{\alpha_k}{\mu}} \left( \mathcal{X}_{(k,t)} + \frac{\mathbf{Y}_{k,t}}{\mu} \right), \quad k = 1, 2. \quad (8) $$

2) Update $\mathbf{M}_3$:

$$ 
\begin{aligned}
\mathbf{M}_{3,t+1} = & \arg \min_{\mathbf{M}_3} \alpha_3 \|\mathbf{M}_3\|_{*} \\
& + \frac{\mu}{2} \left\| \mathcal{X}_{(3,t)} - \mathbf{M}_3 + \frac{\mathbf{Y}_{3,t}}{\mu} \right\|_{\mathrm{F}}^2 + \frac{\beta}{2} \left\| \mathbf{M}_3 \mathbf{D} - \mathcal{L}_{(3)} + \frac{\mathbf{Y}'_{1,t}}{\beta} \right\|_{\mathrm{F}}^2 \\
& + \frac{\gamma}{2} \left\| \mathbf{SM}_3 - \mathcal{H}_{(3)} + \frac{\mathbf{Y}'_{2,t}}{\gamma} \right\|_{\mathrm{F}}^2.
\end{aligned}
\quad (9) $$

Unlike (6), (9) usually does not have a closed-form solution because $\mathbf{D}$ and $\mathbf{S}$ are not identities, which precludes direct application of singular-value thresholding. To overcome this difficulty, we propose to linearize the quadratic terms, thereby solving the problem while guaranteeing convergence. This linearization, coupled with an adaptive rule for updating the penalty parameter, has been widely deployed to derive approximate solutions with closed form for subproblems considered in ADMM [25]. Specifically, let $h(\mathbf{M}_3)$ equal the sum of the quadratic terms in (9):

$$ h(\mathbf{M}_3) = \frac{\mu}{2} \left\| \mathcal{X}_{(3,t)} - \mathbf{M}_3 + \frac{\mathbf{Y}_{3,t}}{\mu} \right\|_{\mathrm{F}}^2 \\
+ \frac{\beta}{2} \left\| \mathbf{M}_3 \mathbf{D} - \mathcal{L}_{(3)} + \frac{\mathbf{Y}'_{1,t}}{\beta} \right\|_{\mathrm{F}}^2 \\
+ \frac{\gamma}{2} \left\| \mathbf{SM}_3 - \mathcal{H}_{(3)} + \frac{\mathbf{Y}'_{2,t}}{\gamma} \right\|_{\mathrm{F}}^2. \quad (10) $$

Then, by linearizing, we have

$$ h(\mathbf{M}_3) \approx h(\mathbf{M}_{3,t}) + \langle \nabla h(\mathbf{M}_{3,t}), \mathbf{M}_3 - \mathbf{M}_{3,t} \rangle + \frac{\tau}{2} \|\mathbf{M}_3 - \mathbf{M}_{3,t}\|_{F}^{2} \quad (11) $$

where $\tau > 0$ is a proximal parameter, and $\nabla h(\mathbf{M}_{3,t})$ is the gradient of $h(\mathbf{M}_3)$ at $\mathbf{M}_{3,t}$ which is

$$ 
\begin{aligned}
\nabla h(\mathbf{M}_{3,t}) = & -\mu \left( \mathcal{X}_{(3,t)} - \mathbf{M}_{3,t} + \frac{\mathbf{Y}_{3,t}}{\mu} \right) \\
& + \beta \left( \mathbf{M}_{3,t} \mathbf{D} - \mathcal{L}_{(3)} + \frac{\mathbf{Y}'_{1,t}}{\beta} \right) \mathbf{D}^{T} \\
& + \gamma \mathbf{S}^{T} \left( \mathbf{SM}_{3,t} - \mathcal{H}_{(3)} + \frac{\mathbf{Y}'_{2,t}}{\gamma} \right).
\end{aligned}
\quad (12) $$

Particularly, $\tau$ is computed via imposing a Lipschitz continuity constraint on $h(\mathbf{M}_3)$, yielding

$$ \tau = \mu + \beta \|\mathbf{DD}^T\|_2 + \gamma \|\mathbf{S}^T \mathbf{S}\|_2. \quad (13) $$

Replacing the quadratic term in (9) with (11) gives

$$ 
\begin{aligned}
\mathbf{M}_{3,t+1} = & \arg \min_{\mathbf{M}_3} \alpha_3 \|\mathbf{M}_3\|_{*} + h(\mathbf{M}_{3,t}) \\
& + \langle \nabla h(\mathbf{M}_{3,t}), \mathbf{M}_3 - \mathbf{M}_{3,t} \rangle + \frac{\tau}{2} \|\mathbf{M}_3 - \mathbf{M}_{3,t}\|_{F}^{2}
\end{aligned}
\quad (14) $$
which can then be rewritten as

$$ \mathbf{M}_{3,t+1} = \arg\min_{\mathbf{M}_3} \alpha_3 \|\mathbf{M}_3\|_* + h(\mathbf{M}_{3,t}) + \frac{\tau}{2} \left\| \mathbf{M}_3 - \mathbf{M}_{3,t} + \frac{\nabla h(\mathbf{M}_{3,t})}{\tau} \right\|_F^2 - \frac{1}{2\tau} \left\| \nabla h(\mathbf{M}_{3,t}) \right\|_F^2. \quad (15) $$

Introducing singular-value thresholding, the update rule for $\mathbf{M}_3$ becomes

$$ \mathbf{M}_{3,t+1} = \mathcal{D}_{\frac{\alpha_3}{\tau}} \left( \mathbf{M}_{3,t} - \frac{\nabla h(\mathbf{M}_{3,t})}{\tau} \right) \quad (16) $$

or

$$ \mathbf{M}_{3,t+1} = \mathcal{D}_{\frac{\alpha_3}{\tau}} \left( \mathbf{M}_{3,t} + \frac{\mu}{\tau} \left( \mathcal{X}_{(3),t} - \mathbf{M}_{3,t} + \frac{\mathbf{Y}_{3,t}}{\mu} \right) - \frac{\beta}{\tau} \left( \mathbf{M}_{3,t}\mathbf{D} - \mathcal{L}_{(3)} + \frac{\mathbf{Y}'_{1,t}}{\beta} \right) \mathbf{D}^T - \frac{\gamma}{\tau} \mathbf{S}^T \left( \mathbf{S}\mathbf{M}_{3,t} - \mathcal{H}_{(3)} + \frac{\mathbf{Y}'_{2,t}}{\gamma} \right) \right). \quad (17) $$

3) Update $\mathcal{X}$: We have

$$ \mathcal{X}_{(k),t+1} = \arg\min_{\mathcal{X}_{(k)}} \left\langle \mathbf{Y}_{k,t}, \mathcal{X}_{(k)} - \mathbf{M}_{k,t+1} \right\rangle + \frac{\mu}{2} \|\mathcal{X}_{(k)} - \mathbf{M}_{k,t+1}\|_F^2 \quad (18) $$

for $k \in \{1, 2, 3\}$, which can be rewritten as

$$ \mathcal{X}_{(k),t+1} = \arg\min_{\mathcal{X}_{(k)}} \frac{\mu}{2} \left\| \mathcal{X}_{(k)} - \mathbf{M}_{k,t+1} + \frac{\mathbf{Y}_{k,t}}{\mu} \right\|_{\mathrm{F}}^{2}. \quad (19) $$

The closed-form solution of (19) is

$$ \mathcal{X}_{(k),t+1} = \mathbf{M}_{k,t+1} - \frac{\mathbf{Y}_{k,t}}{\mu}. \quad (20) $$

Hence, we have ultimately

$$ \begin{aligned} \mathcal{X}_{t+1} &= \frac{1}{3} \sum_{k=1}^{3} \mathrm{fold}_k(\mathcal{X}_{(k),t+1}) \\ &= \frac{1}{3} \sum_{k=1}^{3} \mathrm{fold}_k \left( \mathbf{M}_{k,t+1} - \frac{\mathbf{Y}_{k,t}}{\mu} \right). \end{aligned} \quad (21) $$

4) Update $\mathbf{Y}_k$ for $k \in \{1, 2, 3\}$ as well as $\mathbf{Y}'_1$, and $\mathbf{Y}'_2$:

$$ \mathbf{Y}_{k,t+1} = \mathbf{Y}_{k,t} + \mu (\mathcal{X}_{(k),t+1} - \mathbf{M}_{k,t+1}) \quad (22) $$

$$ \mathbf{Y}'_{1,t+1} = \mathbf{Y}'_{1,t} + \beta (\mathbf{M}_{3,t+1}\mathbf{D} - \mathcal{L}_{(3)}) \quad (23) $$

$$ \mathbf{Y}'_{2,t+1} = \mathbf{Y}'_{2,t} + \gamma (\mathbf{S}\mathbf{M}_{3,t+1} - \mathcal{H}_{(3)}). \quad (24) $$

The final optimization procedure for the proposed LRTA fusion framework is presented in Algorithm 1. We note that we initialize the algorithm with an estimation of $\mathcal{X}$ created by a bicubic interpolation of HSI $\mathcal{L}$ with the ratio $R = M/M'$ applied in each spectral band; this initial $\mathcal{X}$ is then used to initialize each $\mathbf{M}_k$ via unfolding.

We now consider the convergence of the proposed LRTA. As the trace norm defined in (1) is a convex proxy for rank, the objective function considered in the fusion problem (3)

Algorithm 1 The LRTA HS-MS Fusion Algorithm for (3)

1: **Input:** HSI $\mathcal{L}$, MSI $\mathcal{H}$, PSF $\mathbf{D}$, PSK $\Psi$, and SRF $\mathbf{S}$; parameters $\alpha_k$, $\mu$, $\beta$, $\gamma$, $T_{\max}$;

2: **Initialization:** $\mathbf{Y}'_{1,0} = \mathbf{Y}'_{2,0} = \mathbf{Y}_{k,0} = \mathbf{0}$, $k \in \{1, 2, 3\}$; $\eta_1 = 10^{-4}$, $\eta_2 = 10^{-5}$, $t = 0$;

3: $\mathcal{X}_t$ = interpolation($\mathcal{L}$, $\mathcal{R}$, “bicubic”);

4: $\mathbf{M}_{k,t} = \mathrm{unfold}_k(\mathcal{X}_t)$ for $k \in \{1, 2, 3\}$;

5: **repeat**

6: update $\mathbf{M}_{k,t+1}$ via (8) for $k = 1, 2$;

7: update $\mathbf{M}_{3,t+1}$ via (17);

8: update $\mathcal{X}_{t+1}$ via (21);

9: update Lagrange multipliers $\mathbf{Y}_{k,t+1}$, $\mathbf{Y}'_{1,t+1}$ and $\mathbf{Y}'_{2,t+1}$ via (22), (23), and (24);

10: $t \leftarrow t + 1$;

11: **until** $\max(\|\mathbf{M}_{3,t}\mathbf{D} - \mathcal{L}_{(3)}\|_{\mathrm{F}}, \|\mathbf{S}\mathbf{M}_{3,t} - \mathcal{H}_{(3)}\|_{\mathrm{F}}) < \eta_1$ **and** $\|\mathcal{X}_t - \mathcal{X}_{t-1}\|_{\mathrm{F}} < \eta_2$ **or** $t > T_{\max}$;

12: **return** $\mathcal{X} = \mathcal{X}_t$.

is convex, but not differentiable. However, based on block coordinate descent theory, the introduction of variables $\mathbf{M}_k$ in the augmented Lagrangian (4) leads to closed solutions for both $\mathcal{X}$ [i.e., (21)] and $\mathbf{M}_k$ [i.e., (8) and (17)] [22], [29].

C. Extension to Stripe Noise or Other Missing Values

The estimation of missing values within degraded visual data sets has arisen in many applications (e.g., [22], [29], [30]). Specifically in the remote-sensing area, the HSI tends to be degraded with missing stripes due to dead sensor elements being swept across the scene in pushbroom acquisition. Other malfunctions during the acquisition process may give rise to other missing values due to dead pixels in a line or across a patch; additionally, various outliers may be identified within the image as well as [31], or portions of the image may be obscured by cloud cover. We now generalize the LRTA HS-MS fusion framework discussed above in order to fuse HSI with MSI while recovering corrupted regions within the HSI.

Specifically, let us suppose that the values of HSI $\mathcal{L}$ in the set $\Omega$ are known while the remaining elements are missing, which is denoted as $\mathcal{L}_{\Omega}$. The operation that maps $\mathcal{L}$ into $\mathcal{L}_{\Omega}$ is denoted as $\mathcal{L}_{\Omega} = \mathcal{W} \odot \mathcal{L}$, where $\mathcal{W}$ is a binary indicator tensor, where $\mathcal{W}$ is 0 denotes that the corresponding HSI value is missing, while 1 means the HSI value is present and known. We note that, in real applications, such as for MODIS and Landsat products, poor-quality data can often be identified via provided quality-assurance (QA) layer documents, which, in turn, can be used to determine $\mathcal{W}$ [32].

Under this missing-value scenario, we generalize (3) as

$$ \begin{aligned} & \min_{\mathcal{X}} \sum_{k=1}^{3} \alpha_k \|\mathcal{X}_{(k)}\|_{*} \\ & \text{s.t. } \mathcal{W}_{(3)} \odot (\mathcal{X}_{(3)}\mathbf{D}) = (\mathcal{W} \odot \mathcal{L})_{(3)} \text{ and } \mathbf{S}\mathcal{X}_{(3)} = \mathcal{H}_{(3)}. \end{aligned} \quad (25) $$

The solution to (25) is very similar to that outlined above for the original problem (3). Specifically, the update of $\mathbf{M}_3$ in (9)
is revised to

$$ \mathbf{M}_{3,t+1} = \arg \min_{\mathbf{M}_3} \alpha_3 \|\mathbf{M}_3\|_* + \frac{\mu}{2} \left\| \chi_{(3),t} - \mathbf{M}_3 + \frac{\mathbf{Y}_{3,t}}{\mu} \right\|_{\mathrm{F}}^2 \\ + \frac{\beta}{2} \left\| W_{(3)} \odot (\mathbf{M}_3 \mathbf{D} - \mathcal{L}_{(3)}) + \frac{\mathbf{Y}'_{1,t}}{\beta} \right\|_{\mathrm{F}}^2 \\ + \frac{\gamma}{2} \left\| \mathbf{S}\mathbf{M}_3 - \mathcal{H}_{(3)} + \frac{\mathbf{Y}'_{2,t}}{\gamma} \right\|_{\mathrm{F}}^2. \quad (26) $$

Following a linearization strategy similar to that which produced (17) yields

$$ \mathbf{M}_{3,t+1} = \mathcal{D}_{\frac{\alpha_3}{\tau}} \left( \mathbf{M}_{3,t} - \frac{\nabla g(\mathbf{M}_{3,t})}{\tau} \right) \quad (27) $$

where the gradient $\nabla g(\mathbf{M}_{3,t})$ is

$$ \nabla g(\mathbf{M}_{3,t}) = -\mu \left( \chi_{(3),t} - \mathbf{M}_{3,t} + \frac{\mathbf{Y}_{3,t}}{\mu} \right) \\ + \beta \left[ W_{(3)} \odot (\mathbf{M}_{3,t} \mathbf{D} - \mathcal{L}_{(3)}) + \frac{\mathbf{Y}'_{1,t}}{\beta} \right] \mathbf{D}^T \\ + \gamma \mathbf{S}^T \left( \mathbf{S}\mathbf{M}_{3,t} - \mathcal{H}_{(3)} + \frac{\mathbf{Y}'_{2,t}}{\gamma} \right). \quad (28) $$

Additionally, the update of $\mathbf{Y}'_{1,t+1}$ in (23) is revised to

$$ \mathbf{Y}'_{1,t+1} = \mathbf{Y}'_{1,t} + \beta W_{(3)} \odot (\mathbf{M}_{3,t+1} \mathbf{D} - \mathcal{L}_{(3)}). \quad (29) $$

The other steps of Algorithm 1 remain unchanged.

## IV. EXPERIMENTAL STUDY

We now compare the LRTA approach proposed in Section III to several classic and state-of-the-art HS-MS fusion algorithms. To maximize the fairness and transparency of the comparison, we first analyze the resolution-enhancement performance of all methods both quantitatively as well as visually. Then, we further verify the restoration ability of the proposed LRTA for the case in which the HSI is degraded with stripe noise. Specifically, we compare with coupled nonnegative matrix factorization (CNMF) [33], HySure [7], coupled sparse tensor factorization (CSTF) [6], and low tensor multi-rank regularization (LTMR) [12]. CNMF and HySure are representatives of more classical techniques that rely on matricization of the HSI and MSI image cubes. In the former, nonnegative matrix factorization is applied to matrix versions of both image cubes, and the factors interchanged to complete the fusion, while, in the latter, linear unmixing factors the matricized HSI within a convex optimization regularized by total variation. We also compare with CSTF and LTMR, both of which are tensor-based methods. CSTF is representative of the common tensor-decomposition approach to fusion which involves the estimation of three-factor matrices and a corresponding sparse core tensor with dictionary learning applied to produce the factor matrices. Finally, LTMR is a recent tensor-based technique not driven by tensor decomposition; instead, a low tensor multi-rank regularization is applied to a subspace learned from the HSI along with patch clustering applied to the MSI to produce the fused image.

### A. Experimental Data Sets

To gauge the performance of the proposed LRTA HS-MS fusion algorithm, two synthetic data sets and one real data set from different categories of scenes (e.g., vegetation, urban, suburban) are considered, featuring different geographical and sensor characteristics. The well-known Pavia University and Indian Pines HSI data sets are used to create synthetic data for fusion, while the Zhangye data set consists of a real, co-located HSI/MSI image pair. Each of these data sets are described below.

1) Synthetic Data Sets: The Pavia University data set<fnref n="2" /> was acquired by a Reflective Optics System Imaging Spectrometer (ROSIS) sensor over the University of Pavia in northern Italy. The original data has 610 × 610 pixels with 103 spectral bands. However, some of the pixels were discarded before use as they contain no information. The commonly used version covers 610 × 340 pixels with spectral wavelength range of 430–860 nm and 1.3-m spatial resolution. Ninety three spectral bands are preserved after removing the water-absorption bands, and the top-left 256 × 256 × 93 subimage of data set is adopted as the reference imagery.

The Indian Pines data set [34] was gathered by an Airborne Visible/Infrared Imaging Spectrometer (AVIRIS) sensor in northwestern Indiana. The image comprises 145 × 145 pixels and covers 224 spectral bands with a wavelength range of 400–2450 nm and a 20-m spatial resolution. Two-thirds of agriculture region and one-third of forest region or regions of natural perennial vegetation constitute the majority of the scene. The available 220 spectral bands were reduced to 200 bands by removing water-absorption bands (bands 104–108, 150–163, and 220).

In prior literature, the Pavia University and Indian Pines data sets are frequently used as reference images for both qualitative and quantitative assessments of HS-MS fusion (e.g., [3]). A typical approach applies sensor simulations both spatially and spectrally to the original HSI to generate artificial HSI and MSI cubes, respectively; the original HSI image is then treated as the ground-truth reference HS²I for quantitative comparisons.

Specifically, spatial simulation is performed to generate the low-spatial-resolution HSI using the $R \times R$ PSK $\Psi$; here, we use both a spatial Gaussian as well as a simple block average as $\Psi$; downsampling ratios $R$ of both 8 and 4 are evaluated to correspond to different levels of spatial resolution. Spectral simulation is performed to generate the MSI by degrading the reference image in the spectral domain using a filter. For Pavia University, an IKONOS-like spectral-response filter is used to generate a 4-band MSI as was done in [12]. For Indian Pines, a uniform spectral-response filter corresponding to Landsat TM bands 1–5 and 7 generates the MSI, which covers the spectral regions 450–520, 520–600, 630–690, 760–900, 1550–1750, and 2080–2350 nm, respectively, as was done in [33].

2) Real Data Set: The Zhangye data set consists of a real co-located pair of HSI and MSI image cubes captured by an

<footnote n="2">²http://www.ehu.eus/ccwintco/index.php?title=Hyperspectral_Remote_Sensing_Scenes</footnote>
airborne hyperspectral imaging system over a suburban area of Zhangye, Gansu Province, China. The acquired HSI and MSI cover the visible and near-infrared spectrum, ranging over 383–1055 nm. The HSI has 48 spectral bands, whereas the MSI has 4 bands. The spatial size of the HSI is 1201 × 200 pixels, whereas the MSI has size 6001 × 1001 pixels, corresponding to a spatial resolution of 5 m and 1 m, respectively. In our experiments, we crop a 400 × 400-subimage from the HSI and a corresponding 80 × 80-subimage from the MSI, choosing an area which contains substantial detail and texture content.

## B. Image-Quality Assessments

Since the original Pavia University and Indian Pines data sets provide a reference HS²I, we can employ a variety of full-reference image-quality measures to quantitatively gauge performance of the fusion techniques under consideration. Specifically, we use peak signal-to-noise ratio (PSNR), root mean square error (RMSE), spectral angle mapper (SAM), the dimensionless global relative error of synthesis (ERGAS) [35], the universal image quality index (UIQI) [36], and the structural similarity index measure (SSIM) [37]. We note that PSNR, ERGAS, UIQI, and SSIM are 2-D image-distortion measures; here, we calculate these measures for each spectral band and then average the results across all bands. On the other hand, SAM is defined between two hyperspectral pixel vectors; we thus average SAM across all pixel vectors in the image. Finally, we note that higher PSNR, UIQI, and SSIM values imply better-quality fusion, whereas lower RMSE, SAM, and ERGAS values correspond to the better results.

## C. Parameter Tuning

We now analyze the effect of various parameters on the performance of LRTA before turning our attention to comparison with the other techniques. In the proposed LRTA, two groups of parameters must be considered: trace-norm parameters $\alpha_k$ and penalty parameters $\mu$, $\beta$, and $\gamma$.

The $\alpha_k$ parameters in (3) control the relative intensity of the low-rank constraints on each mode of the fusion result $\mathcal{X}$. Hence, the values of $\alpha_k$ can influence the fusion and restoration results significantly. In [22], which considers traditional RGB imagery, the authors set $\alpha_1 = \alpha_2 = \alpha_3 = 1/3$ and do not discuss them further. In contrast, for HSI, the $\alpha_k$ values relate to the rank in both spatial and spectral directions, and there have been several works [10], [38], [39] which focus on estimating these ranks. As is common, we set $\alpha_1$ equals to $\alpha_2$ as they both correspond to spatial directions and note that, typically, for HSI, the spatial-mode ranks are much larger than that of the spectral mode [10], [38]. We follow [10] in setting the weights as

$$ \alpha_k = \omega_k \sqrt{\frac{I_{\max}}{I_k}}, \quad k \in \{1, 2, 3\} \qquad (30) $$

where $I_k$ is the size of $\mathcal{X}$ on each mode, that is, $(I_1, I_2, I_3) = (M, N, D)$, $I_{\max} = \max_k(I_k)$, and $\omega_k$ is a ratio parameter. The tuning of $\omega_k$ is illustrated in Fig. 2, wherein we fix $\omega_3$ to

<figure><img src="image_4.png" /><figcaption>Fig. 2. Tuning of $\omega_1 = \omega_2$ for LRTA. (a, b) Indian Pines; (c, d) Pavia University. Red: 4 × 4 average $\Psi$; Blue: 4 × 4 Gaussian $\Psi$.</figcaption></figure>

be 100, and tune such that $\omega_1 = \omega_2$. Four cases are considered during the tuning for each data set: Fig. 2(a) and (c) is pure HS-MS fusion result wherein the HSI is degraded spatially by only the average or Gaussian kernel without missing values (see Section IV-D for details), while Fig. 2(b) and (d) combines fusion and restoration, wherein the HSI is degraded the spatial kernels as well as missing values (see Section IV-E). It can be observed that in all situations, performance for values $\omega_1, \omega_2 \le 0.1$ shows little variation, whereas performance falls off dramatically in all cases for $\omega_1, \omega_2 \ge 2.0$. With a goal to set $\omega_1$ and $\omega_2$ such that universally good performance is achieved, we would want to choose a value in the range $\omega_1, \omega_2 \in [0.1, 2.0]$. Accordingly, we set $\omega_1 = \omega_2 = 1.0$ for the remainder of the experimental results, which is the same choice made in [10]. The final $\alpha_k$ values are then obtained via (30) and then normalizing such that $\sum_{k=1}^3 \alpha_k = 1$.

As for the LRTA penalty parameters, $\mu$ controls the weighting of the fidelity term in (4), while parameters $\beta$ and $\gamma$ control the tradeoff between spectral and spatial information in (4). As such, penalty parameters typically relate more to the speed of convergence than to the final optimization result, we seek appropriate values that yield convergence without excessive iteration. To this end, we arbitrarily set $\mu = 0.01$ and then tune $\beta$ and $\gamma$ over 100 iterations, choosing values for both from the set {0.01, 0.05, 0.1, 0.5, 1.0, 5.0}. The resulting relative mean error is depicted in Fig. 3, indicating performance is near optimal for both data sets between 0.5 and 1.0, so we set $\beta = \gamma = 0.5$ for the remainder of the experiments. We note that these nonzero values of $\beta$ and $\gamma$ demonstrate that the spatial and spectral terms in (4) impact the performance of LRTA as intended.

Finally, Fig. 4 illustrates convergence of LRTA in terms of the RMSE. We observe that LRTA largely reaches convergence after about 30 iterations for the two synthetic data sets. Consequently, we double this for a margin of error, setting $T_{\max} = 60$ in Algorithm 1 for the remainder of the experiments.
<table><caption>TABLE I
FUSION PERFORMANCE FOR THE SYNTHETIC DATA SETS
CREATED USING THE VARIOUS PSK OPERATORS Ψ</caption><thead><tr><th colspan="12">Pavia University</th></tr><tr><th rowspan="2"></th><th colspan="6">8 × 8 Gaussian Ψ</th><th colspan="6">8 × 8 Average Ψ</th></tr><tr><th>PSNR</th><th>RMSE</th><th>SAM</th><th>ERGAS</th><th>UIQI</th><th>SSIM</th><th>PSNR</th><th>RMSE</th><th>SAM</th><th>ERGAS</th><th>UIQI</th><th>SSIM</th></tr></thead><tbody><tr><td>Bicubic</td><td>22.41</td><td>19.492</td><td>8.5019</td><td>5.7228</td><td>0.4527</td><td>0.5108</td><td>22.59</td><td>19.079</td><td>8.5228</td><td>5.5981</td><td>0.447</td><td>0.5121</td></tr><tr><td>CNMF</td><td>39.04</td><td>2.8846</td><td>2.2358</td><td>0.8610</td><td>0.9909</td><td>0.9839</td><td>41.44</td><td>2.2319</td><td>1.9232</td><td>0.6397</td><td>0.9936</td><td>0.9885</td></tr><tr><td>HySure</td><td>41.69</td><td>2.2725</td><td>2.2788</td><td>0.6389</td><td>0.9923</td><td>0.9868</td><td>42.01</td><td>2.1885</td><td>2.2537</td><td>0.6186</td><td>0.9929</td><td>0.9885</td></tr><tr><td>CSTF</td><td>43.17</td><td>1.9608</td><td>2.0480</td><td>0.5515</td><td>0.9932</td><td>0.9871</td><td>43.05</td><td>2.0109</td><td>2.0861</td><td>0.5623</td><td>0.9930</td><td>0.9872</td></tr><tr><td>LTMR</td><td>42.38</td><td>2.5125</td><td>2.3400</td><td>0.6378</td><td>0.9916</td><td>0.9865</td><td>43.81</td><td>2.0772</td><td>2.0947</td><td>0.5250</td><td>0.9933</td><td>0.9860</td></tr><tr><td>LRTA</td><td><b>43.85</b></td><td><b>1.9122</b></td><td><b>1.9083</b></td><td><b>0.5143</b></td><td><b>0.9938</b></td><td><b>0.9893</b></td><td><b>44.38</b></td><td><b>1.7778</b></td><td><b>1.7942</b></td><td><b>0.4805</b></td><td><b>0.9946</b></td><td><b>0.9906</b></td></tr></tbody></table><table><thead><tr><th colspan="12">Indian Pines</th></tr><tr><th rowspan="2"></th><th colspan="6">4 × 4 Gaussian Ψ</th><th colspan="6">4 × 4 Average Ψ</th></tr><tr><th>PSNR</th><th>RMSE</th><th>SAM</th><th>ERGAS</th><th>UIQI</th><th>SSIM</th><th>PSNR</th><th>RMSE</th><th>SAM</th><th>ERGAS</th><th>UIQI</th><th>SSIM</th></tr></thead><tbody><tr><td>Bicubic</td><td>40.76</td><td>4.8746</td><td>2.4072</td><td>1.2517</td><td>0.7951</td><td>0.8845</td><td>40.74</td><td>4.8883</td><td>2.4205</td><td>1.2555</td><td>0.7928</td><td>0.8835</td></tr><tr><td>CNMF</td><td>45.10</td><td>2.5363</td><td>1.6965</td><td>0.7768</td><td>0.8718</td><td>0.9645</td><td>45.23</td><td>2.5214</td><td>1.6922</td><td>0.7725</td><td>0.8720</td><td>0.9646</td></tr><tr><td>HySure</td><td>46.12</td><td>2.6203</td><td>1.7748</td><td>0.7683</td><td>0.8773</td><td>0.9640</td><td>46.09</td><td>2.6195</td><td>1.7735</td><td>0.7700</td><td>0.8750</td><td>0.9636</td></tr><tr><td>CSTF</td><td>43.16</td><td>3.1561</td><td>2.0158</td><td>0.9171</td><td>0.8443</td><td>0.9491</td><td>42.83</td><td>3.2945</td><td>2.1046</td><td>0.9537</td><td>0.8467</td><td>0.9491</td></tr><tr><td>LTMR</td><td>43.86</td><td>3.0929</td><td>1.9466</td><td>0.8747</td><td>0.8692</td><td>0.9569</td><td>43.88</td><td>3.0837</td><td>1.9476</td><td>0.8727</td><td>0.8696</td><td>0.9571</td></tr><tr><td>LRTA</td><td><b>47.66</b></td><td><b>2.4176</b></td><td><b>1.6148</b></td><td><b>0.6941</b></td><td><b>0.8962</b></td><td><b>0.9663</b></td><td><b>47.68</b></td><td><b>2.4161</b></td><td><b>1.6136</b></td><td><b>0.6937</b></td><td><b>0.8962</b></td><td><b>0.9664</b></td></tr></tbody></table>

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

<figure><img src="image_6.png" /><figcaption>Fig. 3. Tuning of parameters $\gamma$ and $\beta$ for LRTA over 100 iterations. (a) Indian Pines, (b) Pavia University. Note: for visualization purposes, the figures plot the reciprocal of relative mean error (i.e., larger is better).</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 4. Convergence in terms of relative mean error for two synthetic data sets. Indian Pines: orange, Pavia University: blue.</figcaption></figure>

## D. Results and Analysis on Spatial Resolution Enhancement

Table I tabulates quantitative performance of the various algorithms on the two synthetic data sets. It can be observed that the proposed LRTA achieves the best performance across all the cases. Furthermore, the fusion results are depicted visually in Figs. 5 and 6 wherein we observe that the proposed LRTA achieves better representation of spatial edge and texture details while more successfully achieving an HS²I with a satisfactory high spatial resolution and spectrum preservation. Table I as well as Figs. 5 and 6 also include, as a simple benchmark, straightforward spatial bicubic interpolation of the HSI to produce a high-resolution HS²I (denoted as “Bicubic”), which, unsurprising, yields uniformly poor performance.

We now investigate performance of the various fusion algorithms on a real HSI/MSI pair using the Zhangye data set. In the preceding synthetic experiments, we had known spatial and spectral sensor observation models; however, in real applications, these are rarely known exactly. Consequently, we must estimate both the PSF as well as the SRF from the observed MSI and HSI data pair. Here, we use the estimation method proposed in [7] to estimate the SRF. For the PSF,

As for the fusion methods to which we compare, the parameters of CNMF and HySure are adjusted as reported in Table III of [3]. For the tensor-based CSTF, we use the parameters determined in [6] (i.e., the number of dictionary atoms are set as $n_s = 12$ and $n_w = n_h = 240$). Finally, for LTMR, the parameters are set as discussed in [12] (i.e., patch size of $7 \times 7$; $K = 201$ groups).
<figure><img src="image_8.png" /><figcaption>Fig. 5. Pseudocolor visualization of fusion results for the Pavia University data set; pseudocolor composed of bands 45, 15, and 20. Colorbar gives the scale for the error maps. “Bicubic” indicates simple spatial bicubic interpolation applied to the HSI to generate the HS²I. Detailed subimages are enlarged in red squares.</figcaption></figure>

<table><caption>TABLE II
EXECUTION TIME (IN SECONDS) FOR THE TWO SYNTHETIC DATA SETS</caption><thead><tr><th colspan="3">Pavia University</th></tr><tr><th></th><th>Gaussian kernel (8 × 8)</th><th>Average kernel (8 × 8)</th></tr></thead><tbody><tr><th>CNMF</th><td>36.6</td><td>37.6</td></tr><tr><th>HySure</th><td>40.6</td><td>42.4</td></tr><tr><th>CSTF</th><td>104.8</td><td>108.4</td></tr><tr><th>LTMR</th><td>86.0</td><td>85.9</td></tr><tr><th>LRTA</th><td>49.3</td><td>50.9</td></tr></tbody></table>

<table><thead><tr><th colspan="3">Indian Pines</th></tr><tr><th></th><th>Gaussian kernel (4 × 4)</th><th>Average kernel (4 × 4)</th></tr></thead><tbody><tr><th>CNMF</th><td>16.6</td><td>17.2</td></tr><tr><th>HySure</th><td>12.9</td><td>12.5</td></tr><tr><th>CSTF</th><td>103.4</td><td>102.5</td></tr><tr><th>LTMR</th><td>40.6</td><td>41.1</td></tr><tr><th>LRTA</th><td>26.7</td><td>27.0</td></tr></tbody></table>

we simply adopt a 5 × 5 average kernel. Due to the lack of a reference image, quantitative results are not possible; we thus illustrate the fusion results visually in Fig. 7. We observe that all methods achieve satisfactory fusion results.

Table II summarizes the computational complexity of the fusion methods under consideration. All experiments were carried out using MATLAB on an Intel® Core i7-4970 computer with 16GB of memory. Note that the computational

cost of CSTF is much higher than that of other techniques due to the dictionary learning on each mode. The proposed LRTA achieves a satisfactory fusion result with significantly less computation than LTMR or CSTF, the other tensor-based methods under consideration.

## E. Experimental Results and Analysis on Restoration Ability

Under the situation wherein the HSI is degraded by stripes (or random noise, or missing data across a spatial patch), fusion methods should be noise resilient. Here, we demonstrate that the proposed LRTA fusion is capable of recovering corrupted data during the HS-MS fusion process. We remark that matricization-based approaches—such as CNMF and HySure—cannot remove noise; hence, independent denoising must be applied as preprocessing to the HSI before fusion for these methods. Similarly, although CSTF and LTMR can remove additive Gaussian noise, they cannot complete missing values due to striping or other sensor failures which must also be corrected via preprocessing.

For the experimental results here, we generate artificial stripes within the HSI. Specifically, stripes are added according to the process outlined in [24]. Given an HSI $L$ of size $M' \times N' \times D$, synthetic stripes are generated by randomly selecting $dN'$ columns in a given band to be striped. In each column to be striped, we choose a random value $v$ and create
<figure><img src="image_9.png" /><figcaption>Fig. 6. Pseudocolor visualization of fusion results for the Indian Pines data set; pseudocolor composed of bands 45, 15, and 20. Colorbar gives the scale for the error maps. “Bicubic” indicates simple spatial bicubic interpolation applied to the HSI to generate the HS²I. Detailed subimages are enlarged in red squares.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 7. Pseudocolor visualization of fusion results for the Zhangye data set; pseudocolor composed of bands 25, 15, and 10. “Bicubic” indicates simple spatial bicubic interpolation applied to the HSI to generate the HS²I. “MSI” denotes a pseudocolor visualization of the original multispectral bands 4, 3, and 2. Detailed subimages are enlarged in red squares.</figcaption></figure>

the stripe by adding $v$ to every pixel value in that column. Here, $d$ is the parameter that controls the spatial density of the striping, $0 < d < 1$, while $v$ is the relative intensity

of a stripe. For each stripe, $v$ is randomly chosen using a uniform distribution over the range $[-r, +r]$. Thus, parameter $r$, $0 < r < 1$, controls the overall intensity of the striping.
<figure><img src="image_11.png" /><figcaption>Fig. 8. Restoration and fusion results of band 50 for the Pavia University data set under striping degradation of the HSI. Colorbar gives the scale for the error maps. “Bicubic” indicates simple spatial bicubic interpolation applied to the HSI; “ASSTV+Bicubic” indicates description of the HSI via ASSTV prior to bicubic interpolation. Detailed subimages are enlarged in red squares.</figcaption></figure>

<table><caption>TABLE III
RESTORATION & FUSION PERFORMANCE FOR THE SYNTHETIC DATA SETS CREATED USING THE VARIOUS PSK OPERATORS Ψ</caption><thead><tr><th rowspan="3"></th><th colspan="10">Pavia University</th></tr><tr><th colspan="5">8 × 8 Gaussian Ψ</th><th colspan="5">8 × 8 Average Ψ</th></tr><tr><th>PSNR</th><th>RMSE</th><th>SAM</th><th>ERGAS</th><th>UIQI</th><th>SSIM</th><th>PSNR</th><th>RMSE</th><th>SAM</th><th>ERGAS</th><th>UIQI</th><th>SSIM</th></tr></thead><tbody><tr><th>Bicubic</th><td>20.34</td><td>24.7511</td><td>23.3626</td><td>7.1560</td><td>0.3789</td><td>0.4228</td><td>20.33</td><td>24.8072</td><td>23.1910</td><td>7.1646</td><td>0.3635</td><td>0.4169</td></tr><tr><th>ASSTV+Bicubic</th><td>24.24</td><td>15.8259</td><td>9.3132</td><td>4.5977</td><td>0.5875</td><td>0.5679</td><td>24.12</td><td>16.0520</td><td>9.4719</td><td>4.6699</td><td>0.5692</td><td>0.5589</td></tr><tr><th>ASSTV+CNMF</th><td>32.92</td><td>6.2253</td><td>5.9327</td><td>1.9027</td><td>0.9706</td><td>0.9618</td><td>33.51</td><td>6.0510</td><td>6.3288</td><td>1.8469</td><td>0.9747</td><td>0.9652</td></tr><tr><th>ASSTV+HySure</th><td>37.20</td><td>3.9491</td><td>3.8889</td><td>1.1719</td><td>0.9812</td><td>0.9732</td><td>38.66</td><td>3.4270</td><td>3.6785</td><td>1.0454</td><td>0.9873</td><td>0.9815</td></tr><tr><th>ASSTV+CSTF</th><td>27.55</td><td>12.2319</td><td>9.9529</td><td>3.4986</td><td>0.8222</td><td>0.7771</td><td>28.09</td><td>11.4745</td><td>9.3862</td><td>3.3050</td><td>0.8466</td><td>0.8056</td></tr><tr><th>ASSTV+LTMR</th><td>32.15</td><td>7.3854</td><td>6.2375</td><td>2.1804</td><td>0.9398</td><td>0.9346</td><td>32.96</td><td>7.2979</td><td>6.6571</td><td>2.1379</td><td>0.9429</td><td>0.9381</td></tr><tr><th>LRTA</th><td><strong>43.34</strong></td><td><strong>2.0296</strong></td><td><strong>2.0402</strong></td><td><strong>0.5485</strong></td><td><strong>0.9933</strong></td><td><strong>0.9889</strong></td><td><strong>43.63</strong></td><td><strong>1.9295</strong></td><td><strong>1.9855</strong></td><td><strong>0.5236</strong></td><td><strong>0.9941</strong></td><td><strong>0.9901</strong></td></tr></tbody></table>

<table><thead><tr><th rowspan="3"></th><th colspan="10">Indian Pines</th></tr><tr><th colspan="5">4 × 4 Gaussian Ψ</th><th colspan="5">4 × 4 Average Ψ</th></tr><tr><th>PSNR</th><th>RMSE</th><th>SAM</th><th>ERGAS</th><th>UIQI</th><th>SSIM</th><th>PSNR</th><th>RMSE</th><th>SAM</th><th>ERGAS</th><th>UIQI</th><th>SSIM</th></tr></thead><tbody><tr><th>Bicubic</th><td>22.08</td><td>20.4823</td><td>13.8273</td><td>10.8219</td><td>0.1851</td><td>0.4792</td><td>22.23</td><td>20.1057</td><td>13.5508</td><td>10.6397</td><td>0.1737</td><td>0.4727</td></tr><tr><th>ASSTV+Bicubic</th><td>29.37</td><td>12.5173</td><td>8.4614</td><td>5.4936</td><td>0.7246</td><td>0.8672</td><td>30.19</td><td>12.3817</td><td>8.3414</td><td>5.3580</td><td>0.7214</td><td>0.8640</td></tr><tr><th>ASSTV+CNMF</th><td>30.41</td><td>11.7047</td><td>8.1657</td><td>5.3852</td><td>0.8173</td><td>0.9405</td><td>31.39</td><td>11.5541</td><td>8.0394</td><td>5.2424</td><td>0.8181</td><td>0.9430</td></tr><tr><th>ASSTV+HySure</th><td>31.05</td><td>11.0532</td><td>7.7482</td><td>5.6558</td><td>0.8187</td><td>0.9135</td><td>31.92</td><td>11.1057</td><td>7.7833</td><td>5.6201</td><td>0.8322</td><td>0.9436</td></tr><tr><th>ASSTV+CSTF</th><td>32.47</td><td>10.6920</td><td>7.4553</td><td>4.8720</td><td>0.7723</td><td>0.9123</td><td>32.61</td><td>10.9289</td><td>7.6316</td><td>4.9586</td><td>0.7876</td><td>0.9276</td></tr><tr><th>ASSTV+LTMR</th><td>30.43</td><td>11.3700</td><td>7.9384</td><td>5.4317</td><td>0.7537</td><td>0.8764</td><td>32.05</td><td>11.0760</td><td>7.7361</td><td>5.1626</td><td>0.8184</td><td>0.9345</td></tr><tr><th>LRTA</th><td><strong>47.40</strong></td><td><strong>2.5616</strong></td><td><strong>1.7326</strong></td><td><strong>0.7337</strong></td><td><strong>0.8884</strong></td><td><strong>0.9644</strong></td><td><strong>46.87</strong></td><td><strong>2.6150</strong></td><td><strong>1.7749</strong></td><td><strong>0.7560</strong></td><td><strong>0.8843</strong></td><td><strong>0.9639</strong></td></tr></tbody></table>

This process is repeated independently in each of the *D* bands of the HSI. Thus, the synthetic stripes are aperiodic spatially and uncorrelated across the spectral bands.

For the techniques that must employ denoising pre-processing, we use anisotropic spatio-spectral total variation (ASSTV) [40] to remove stripes in the HSI prior to fusion.

ASSTV treats the HSI as an entire spatial-spectral cube and imposes total-variation optimization across both spatial and spectral modes simultaneously; meanwhile, the total-variation measure is anisotropic spatially in the sense that it exploits the known directional nature of the stripes to remove them while preserving image edges. We note that we use ASSTV since
<figure><img src="image_12.png" /><figcaption>Fig. 9. Restoration and fusion results of band 50 for the Indian Pines data set under striping degradation of the HSI. Colorbar gives the scale for the error maps. “Bicubic” indicates simple spatial bicubic interpolation applied to the HSI; “ASSTV+Bicubic” indicates description of the HSI via ASSTV prior to bicubic interpolation. Detailed subimages are enlarged in red squares.</figcaption></figure>

it is a relatively prominent description method that has been extensively verified to be effective at removing comprehensive stripes as well as random noise while preserving the edges and detailed information; however, it would be feasible to employ most other description or restoration methods as well.

For LRTA, we directly fuse the HSI and MSI without any preprocessing, while for CNMF, HySure, CSTF, and LTMR, we apply ASSTV to the striped HSI before fusion. As additional, non-fusion-based benchmarks, we also consider simple bicubic interpolation applied directly to the striped HSI (denoted “Bicubic”) as well as bicubic interpolation applied after ASSTV description of the HSI (“ASSTV+Bicubic”). Table III tabulates quantitative recovery and fusion performance of the various algorithms on the two data sets using $r = 0.2$ and $d = 0.6$ to model heavy dense stripes. From Table III, it can be observed that LRTA is capable of successfully fusing the HSI and MSI while simultaneously recovering from the heavy striping that degrades the HSI, as is evident from the fact that the PSNR values, for example, for LRTA reported in Table III are only slightly lower than the corresponding fusion-only results in Table I. On the other hand, the PSNR values for the other fusion methods reported in Table III are on the order of 5–15 dB lower than the corresponding values in Table I. This conclusion is reinforced by Figs. 8 and 9 which present the restoration plus fusion results visually for a single spectral band. It can be observed that while ASSTV+CNMF achieves a smooth result for the

Pavia University data set, it incurs a relatively high degree of error for the Indian Pines data set. Meanwhile, the other fusion-based methods exhibit significant blocking artifacts for both data sets; we believe that this is due to patch- and dictionary-based processing (LTMR and CSTF) as well as nonorthogonal subspace bases (CSTF and HySure) which permit smooth results in the optimization subspace to produce nonsmooth results in the target $HS^2I$ space. On the other hand, LRTA achieves low reconstructed error with fewer visual artifacts.

V. CONCLUSION

In this article, we proposed the LRTA algorithm for simultaneous HSI restoration and spatial resolution enhancement through HS-MS fusion, leveraging a low-rank tensor model to couple the high spatial information from an MSI image to the high spectral correlation emanating from an HSI image. A tensor trace norm was adopted to impose low-rankness spatially as well as spectrally, and a LADMAP-driven solution was devised. While tensor-based fusion has appeared in prior literature, typical approaches resorted to a tensor decomposition that often led to a complicated patch-based dictionary-learning fusion hindered by the lack of a tensor equivalent to Parseval’s theorem. In contrast, the proposed LRTA exploits concepts from the field of tensor completion to directly impose a low-rank property spatially and spectrally on the target $HS^2I$ while avoiding the computationally complex
patch clustering and dictionary learning common to competing techniques. Additionally, small modifications to the basic LRTA optimization resulted in a fusion process robust to missing HSI values such as those that can result from heavy striping in real hyperspectral sensors. In experimental results on both synthetic imagery as well as real imagery from an actual hyperspectral sensor, the proposed LRTA algorithm demonstrated effective fusion as well as preservation of both spatial details and texture, yielding significantly improved image quality when compared to other state-of-the-art fusion methods, along with effective restoration under conditions of HSI striping.

## REFERENCES

[1] Y. Zhang, S. De Backer, and P. Scheunders, "Noise-resistant wavelet-based Bayesian fusion of multispectral and hyperspectral images," *IEEE Trans. Geosci. Remote Sens.*, vol. 47, no. 11, pp. 3834-3843, Nov. 2009.

[2] L. Loncan et al., "Hyperspectral pansharpening: A review," *IEEE Trans. Geosci. Remote Sens.*, vol. 3, no. 3, pp. 27-46, Sep. 2015.

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

[4] Y. Xu, Z. Wu, J. Chanussot, and Z. Wei, "Nonlocal patch tensor sparse representation for hyperspectral image super-resolution," *IEEE Trans. Image Process.*, vol. 28, no. 6, pp. 3034-3047, Jun. 2019.

[5] R. Dian, L. Fang, and S. Li, "Hyperspectral image super-resolution via non-local sparse tensor factorization," in *Proc. IEEE Conf. Comput. Vis. Pattern Recognit. (CVPR)*, Honolulu, HI, USA, Jul. 2017, pp. 3862-3871.

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

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

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

[9] L. Oeding, E. Robeva, and B. Sturmfels, "Decomposing tensors into frames," *Adv. Appl. Math.*, vol. 73, pp. 125-153, Feb. 2016.

[10] K. Zhang, M. Wang, S. Yang, and L. Jiao, "Spatial-spectral-graph-regularized low-rank tensor decomposition for multispectral and hyperspectral image fusion," *IEEE J. Sel. Topics Appl. Earth Observ. Remote Sens.*, vol. 11, no. 4, pp. 1030-1040, Apr. 2018.

[11] T. Xu, T.-Z. Huang, L.-J. Deng, X.-L. Zhao, and J. Huang, "Hyperspectral image superresolution using unidirectional total variation with tucker decomposition," *IEEE J. Sel. Topics Appl. Earth Observ. Remote Sens.*, vol. 13, pp. 4381-4398, Jul. 2020.

[12] R. Dian and S. Li, "Hyperspectral image super-resolution via subspace-based low tensor multi-rank regularization," *IEEE Trans. Image Process.*, vol. 28, no. 10, pp. 5135-5146, Oct. 2019.

[13] Y. Xu, Z. Wu, J. Chanussot, and Z. Wei, "Hyperspectral images super-resolution via learning high-order coupled tensor ring representation," *IEEE Trans. Neural Netw. Learn. Syst.*, vol. 31, no. 11, pp. 4747-4760, Nov. 2020.

[14] W. He, Y. Chen, N. Yokoya, C. Li, and Q. Zhao, "Hyperspectral super-resolution via coupled tensor ring factorization," 2020, arXiv:2001.01547. [Online]. Available: http://arxiv.org/abs/2001.01547

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

[16] J. Xue, Y. Zhao, W. Liao, and J. C.-W. Chan, "Nonlocal low-rank regularized tensor decomposition for hyperspectral image denoising," *IEEE Trans. Geosci. Remote Sens.*, vol. 57, no. 7, pp. 5174-5189, Jul. 2019.

[17] J. Xue, Y. Zhao, W. Liao, and J. Chan, "Nonlocal tensor sparse representation and low-rank regularization for hyperspectral image compressive sensing reconstruction," *Remote Sens.*, vol. 11, no. 2, p. 193, Jan. 2019.

[18] Y. Wang, J. Peng, Q. Zhao, Y. Leung, X.-L. Zhao, and D. Meng, "Hyperspectral image restoration via total variation regularized low-rank tensor decomposition," *IEEE J. Sel. Topics Appl. Earth Observ. Remote Sens.*, vol. 11, no. 4, pp. 1227-1243, Apr. 2018.

[19] Z. Huang, S. Li, L. Fang, H. Li, and J. A. Benediktsson, "Hyperspectral image denoising with group sparse and low-rank tensor decomposition," *IEEE Access*, vol. 6, pp. 1380-1390, Dec. 2018.

[20] Y. Chen, W. He, N. Yokoya, and T.-Z. Huang, "Hyperspectral image restoration using weighted group sparsity-regularized low-rank tensor decomposition," *IEEE Trans. Cybern.*, vol. 50, no. 8, pp. 3556-3570, Aug. 2020.

[21] C. J. Hillar and L.-H. Lim, "Most tensor problems are NP-hard," *J. ACM*, vol. 60, no. 6, pp. 1-39, Nov. 2013.

[22] J. Liu, P. Musialski, P. Wonka, and J. Ye, "Tensor completion for estimating missing values in visual data," *IEEE Trans. Pattern Anal. Mach. Intell.*, vol. 35, no. 1, pp. 208-220, Jan. 2013.

[23] P. Rakwatin, W. Takeuchi, and Y. Yasuoka, "Stripe noise reduction in MODIS data by combining histogram matching with facet filter," *IEEE Trans. Geosci. Remote Sens.*, vol. 45, no. 6, pp. 1844-1856, Jun. 2007.

[24] N. Liu, W. Li, R. Tao, and J. E. Fowler, "Wavelet-domain low-rank/group-sparse destriping for hyperspectral imagery," *IEEE Trans. Geosci. Remote Sens.*, vol. 57, no. 12, pp. 10310-10321, Dec. 2019.

[25] J. Yang and X. Yuan, "Linearized augmented lagrangian and alternating direction methods for nuclear norm minimization," *Math. Comput.*, vol. 82, no. 281, pp. 301-329, Jan. 2013.

[26] Z. Lin, R. Liu, and Z. Su, "Linearized alternating direction method with adaptive penalty for low-rank representation," in *Advances in Neural Information Processing Systems 24*, J. Shawe-Taylor, R. S. Zemel, P. L. Bartlett, F. Pereira, and K. Q. Weinberger, Eds. Red Hook, NY, USA: Curran Associates, 2011, pp. 612-620.

[27] W. Li, J. Liu, and Q. Du, "Sparse and low-rank graph for discriminant analysis of hyperspectral imagery," *IEEE Trans. Geosci. Remote Sens.*, vol. 54, no. 7, pp. 4094-4105, Jul. 2016.

[28] B. Madathil et al., "Tensor low rank modeling and its applications in signal processing," 2019, arXiv:1912.03435. [Online]. Available: http://arxiv.org/abs/1912.03435

[29] J.-F. Cai, E. J. Candès, and Z. Shen, "A singular value thresholding algorithm for matrix completion," *SIAM J. Optim.*, vol. 20, no. 4, pp. 1956-1982, Jan. 2010.

[30] M. Signoretto, R. Van de Plas, B. De Moor, and J. A. K. Suykens, "Tensor versus matrix completion: A comparison with application to spectral data," *IEEE Signal Process. Lett.*, vol. 18, no. 7, pp. 403-406, Jul. 2011.

[31] H. Shen and L. Zhang, "A MAP-based algorithm for destriping and inpainting of remotely sensed images," *IEEE Trans. Geosci. Remote Sens.*, vol. 47, no. 5, pp. 1492-1502, May 2009.

[32] F. Gao, J. Masek, M. Schwaller, and F. Hall, "On the blending of the landsat and MODIS surface reflectance: Predicting daily landsat surface reflectance," *IEEE Trans. Geosci. Remote Sens.*, vol. 44, no. 8, pp. 2207-2218, Aug. 2006.

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

[34] P. Gamba, "A collection of data for urban area characterization," in *Proc. IEEE Int. Geosci. Remote Sens. Symp. (IGARSS)*, vol. 1, Anchorage, AK, USA, Sep. 2004, pp. 69-72.

[35] L. Wald, *Data Fusion: Definitions and Architectures*. Paris, France: Presses de l'École des Mines de Paris, 2002.

[36] Z. Wang and A. C. Bovik, "A universal image quality index," *IEEE Signal Process. Lett.*, vol. 9, no. 3, pp. 81-84, Mar. 2002.

[37] Z. Wang, A. C. Bovik, H. R. Sheikh, and E. P. Simoncelli, "Image quality assessment: From error visibility to structural similarity," *IEEE Trans. Image Process.*, vol. 13, no. 4, pp. 600-612, Apr. 2004.

[38] W. Cao et al., "Total variation regularized tensor RPCA for background subtraction from compressive measurements," *IEEE Trans. Image Process.*, vol. 25, no. 9, pp. 4075-4090, Sep. 2016.

[39] T. Yokota, N. Lee, and A. Cichocki, "Robust multilinear tensor rank estimation using higher order singular value decomposition and information criteria," *IEEE Trans. Signal Process.*, vol. 65, no. 5, pp. 1196-1206, Mar. 2017.

[40] Y. Chang, L. Yan, H. Fang, and C. Luo, "Anisotropic spectral-spatial total variation model for multispectral remote sensing image destriping," *IEEE Trans. Image Process.*, vol. 24, no. 6, pp. 1852-1866, Jun. 2015.
<figure><img src="image_13.png" /></figure>

**Na Liu** (Student Member, IEEE) received the B.E. degree from the Beijing University of Chemical Technology, Beijing, China, in 2016, where she is pursuing the Ph.D. degree.

She is also a Visiting Ph.D. Student with Université Grenoble Alpes, hosted by Dr. W. Thuiller in the Laboratoire d'écologie alpine, Grenoble, France, funded by the Chinese Scholarship Council. Her research interests include remote sensing data analysis, hyperspectral data enhancement, and feature extraction. Her supervisors are Dr. W. Li and Dr. J. Chanussot.

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

**Lu Li** received the B.Sc. and M.Sc. degrees from Wuhan University, Wuhan, China, in 2003 and 2007, respectively, and the Ph.D. degree from the Institute of Remote Sensing and Digital Earth, Chinese Academy of Sciences, Beijing, China, in 2017.

From 2017 to 2019, he was a Post-Doctoral Researcher with the College of Information Science and Technology, Beijing University of Chemical Technology, Beijing. He is with the School of Automation, Beijing Information Science and Technology University, Beijing. His research interests include hyperspectral image classification, image processing, and photogrammetry.

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

**Wei Li** (Senior Member, IEEE) received the B.E. degree in telecommunications engineering from Xidian University, Xi'an, China, in 2007, the M.S. degree in information science and technology from Sun Yat-sen University, Guangzhou, China, in 2009, and the Ph.D. degree in electrical and computer engineering from Mississippi State University, Starkville, MS, USA, in 2012.

Subsequently, he spent one year as a Post-Doctoral Researcher at the University of California—Davis, Davis, CA, USA. He was a Professor with the College of Information Science and Technology, Beijing University of Chemical Technology, Beijing, China, from 2013 to 2019. He is a Professor with the School of Information and Electronics, Beijing Institute of Technology, Beijing. His research interests include hyperspectral image analysis, pattern recognition, and data compression.

Dr. Li received the 2015 Best Reviewer Award from the IEEE Geoscience and Remote Sensing Society (GRSS) for his service for IEEE JOURNAL OF SELECTED TOPICS IN APPLIED EARTH OBSERVATIONS AND REMOTE SENSING (JSTARS) and the Outstanding Paper Award from the IEEE International Workshop on Hyperspectral Image and Signal Processing: Evolution in Remote Sensing (Whispers) 2019. He is also serving as an Associate Editor for the IEEE SIGNAL PROCESSING LETTERS and the IEEE JOURNAL OF SELECTED TOPICS IN APPLIED EARTH OBSERVATIONS AND REMOTE SENSING (JSTARS). He has served as Guest Editor for special issue of Journal of Real-Time Image Processing, Remote Sensing, and IEEE JOURNAL OF SELECTED TOPICS IN APPLIED EARTH OBSERVATIONS AND REMOTE SENSING (JSTARS).

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

**Ran Tao** (Senior Member, IEEE) received the B.S. degree from the Electronic Engineering Institute of the PLA, Hefei, China, in 1985, and the M.S. and Ph.D. degrees from the Harbin Institute of Technology, Harbin, China, in 1990 and 1993, respectively.

He has been a Senior Visiting Scholar with the University of Michigan, Ann Arbor, MI, USA, and the University of Delaware, Newark, DE, USA, in 2001 and 2016, respectively. He is a Professor with the School of Information and Electronics, Beijing Institute of Technology, Beijing, China. His research interests include fractional Fourier transform and its applications, theory, and technology for radar and communication systems.

Dr. Tao is also a Fellow of the Institute of Engineering and Technology and the Chinese Institute of Electronics. He is also a member of Wireless Communication and Signal Processing Commission, International Union of Radio Science (URSI). He was a recipient of the National Science Foundation of China for Distinguished Young Scholars in 2006 and a Distinguished Professor of Changjiang Scholars Program in 2009. He was also a recipient of the First Prize of Science and Technology Progress in 2006 and 2007 and the First Prize of Natural Science in 2013, both awarded by the Ministry of Education. He has been a Chief-Professor of the Creative Research Groups with the National Natural Science Foundation of China since 2014. He was a Chief Professor of the Program for Changjiang Scholars and Innovative Research Team in University from 2010 to 2012. He is also the Vice-Chair of the IEEE China Council and the URSI China Council.

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

**James E. Fowler** (Fellow, IEEE) received the B.S. degree in computer and information science engineering and the M.S. and Ph.D. degrees in electrical engineering from The Ohio State University, Columbus, OH, USA, in 1990, 1992, and 1996, respectively.

In 1997, he held a Post-Doctoral Assignment at the Université de Nice-Sophia Antipolis, Nice, France. In 2004, he was a Visiting Professor with the Télécom ParisTech, Paris, France. He is a Billie J. Ball Professor with the Department of Electrical & Computer Engineering, Mississippi State University, Starkville, MS, USA.

Dr. Fowler is also the Vice Chair of the Computational Imaging Technical Committee of the IEEE Signal Processing Society. He is also a Former Chair of the Image, Video, and Multidimensional Signal Processing Technical Committee of the IEEE Signal Processing Society. He was/has been a General Co-Chair of the 2014 IEEE International Conference on Image Processing, Paris, and the Speech, Image, and Video Processing Track Chair of the Asilomar Conference on Signals, Systems, and Computers since 2013. He is also the Publicity Chair of the Data Compression Conference. He is the Editor-in-Chief of IEEE SIGNAL PROCESSING LETTERS. He was previously a Senior Area Editor of IEEE TRANSACTIONS ON IMAGE PROCESSING and Associate Editor of IEEE TRANSACTIONS ON COMPUTATIONAL IMAGING, IEEE TRANSACTIONS ON IMAGE PROCESSING, IEEE TRANSACTIONS ON MULTIMEDIA, and IEEE SIGNAL PROCESSING LETTERS. He is also an Associate Editor for the EURASIP Journal on Image and Video Processing.

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

**Jocelyn Chanussot** (Fellow, IEEE) received the M.Sc. degree in electrical engineering from the Grenoble Institute of Technology (Grenoble INP), Grenoble, France, in 1995, and the Ph.D. degree from the Université de Savoie, Annecy, France, in 1998.

Since 1999, he has been with Grenoble INP, where he is a Professor of signal and image processing. He has also been a Visiting Scholar with Stanford University, Stanford, CA, USA, KTH, Stockholm, Sweden, and NUS, Singapore. Since 2013, he has also been an Adjunct Professor with the University of Iceland, Reykjavik, Iceland. From 2015 to 2017, he was a Visiting Professor with the University of California at Los Angeles (UCLA), Los Angeles, CA, USA. He holds the AXA Chair in Remote Sensing and an Adjunct Professor with the Chinese Academy of Sciences, Aerospace Information Research Institute, Beijing, China. His research interests include image analysis, hyperspectral remote sensing, data fusion, machine learning, and artificial intelligence.

Dr. Chanussot was a member of the Institut Universitaire de France from 2012 to 2017 and a Highly Cited Researcher (Clarivate Analytics/Thomson Reuters). He was also a member of the Machine Learning for Signal Processing Technical Committee of the IEEE Signal Processing Society from 2006 to 2008 and the Program Chair of the IEEE International Workshop on Machine Learning for Signal Processing in 2009. He is also the Founding President of the IEEE Geoscience and Remote Sensing French Chapter, which received the 2010 IEEE GRS-S Chapter Excellence Award. He received multiple outstanding paper awards. He was the Vice President of the IEEE Geoscience and Remote Sensing Society (GRSS), in charge of meetings and symposia from 2017 to 2019. He is also the General Chair of the first IEEE GRSS Workshop on Hyperspectral Image and Signal Processing, Evolution in Remote Sensing (WHISPERS). He was the Co-Chair from 2005 to 2008 and the Chair of the GRS Data Fusion Technical Committee from 2009 to 2011. He is also an Associate Editor of the IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, the IEEE TRANSACTIONS ON IMAGE PROCESSING, and the PROCEEDINGS OF THE IEEE.