# I_0-I_1 Hybrid Total Variation Regularization and its Applications on Hyperspectral Image Mixed Noise Removal and Compressed Sensing

Minghua Wang, Qiang Wang, Jocelyn Chanussot, Danfeng Hong

► To cite this version:

Minghua Wang, Qiang Wang, Jocelyn Chanussot, Danfeng Hong. I_0-I_1 Hybrid Total Variation Regularization and its Applications on Hyperspectral Image Mixed Noise Removal and Compressed Sensing. IEEE Transactions on Geoscience and Remote Sensing, 2021, 59 (9), pp.7695-7710. 10.1109/TGRS.2021.3055516. hal-03429673

HAL Id: hal-03429673

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

Submitted on 25 Sep 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# $l_0-l_1$ Hybrid Total Variation Regularization and Its Applications on Hyperspectral Image Mixed Noise Removal and Compressed Sensing

Minghua Wang<sup>1,†</sup>, Graduate Student Member, IEEE, Qiang Wang<sup>1,†</sup>, Member, IEEE, Jocelyn Chanussot<sup>1,†</sup>, Fellow, IEEE, and Danfeng Hong<sup>1,†</sup>, Member, IEEE

**Abstract**—The total variation (TV) regularization has been widely used in various applications related to hyperspectral (HS) signal and image processing due to its potential in modeling the underlying smoothness of HS data. However, most existing TV norms usually tend to generate spatial oversmoothing or artifacts. To this end, we propose a novel $l_0-l_1$ hybrid TV ($l_0-l_1$HTV) regularization with the applications to HS mixed noise removal and compressed sensing (CS). More specifically, $l_0-l_1$HTV can be regarded as a globally and locally integrated TV regularizer, where the $l_0$ gradient constraint is incorporated into the $l_1$ spatial-spectral TV ($l_1$SSTV). $l_1$SSTV is capable of exploiting the local structure information across both spatial and spectral domains, while the $l_0$ gradient can promote a globally spectral-spatial smoothness by directly controlling the number of nonzero gradients of HS images. This efficient combination considers more comprehensive prior knowledge of HS images, yielding sharper edge preservation and resolving the above drawbacks of existing pure TV norms. More significantly, $l_0-l_1$HTV can be easily injected into HS-related processing models, and an effective algorithm based on the alternating direction method of multipliers (ADMM) is developed to solve the optimization problems. Extensive experiments conducted on several HS data sets substantiate the superiority and effectiveness of the proposed method in comparison with many state-of-the-art methods.

**Index Terms**—$l_0-l_1$ hybrid total variation ($l_0-l_1$HTV), $l_0$ gradient, $l_1$ spatial-spectral TV ($l_1$SSTV), alternating direction method of multipliers (ADMM), compressed sensing (CS), hyperspectral image (HSI) denoising.

Manuscript received June 12, 2020; revised September 17, 2020 and December 10, 2020; accepted January 11, 2021. This work was supported in part by the National Natural Science Foundation of China under Grant 61876054, in part by the Délégation Générale de l'Armement (Project ANR-DGA APHYPIS) under Grant ANR-16 ASTR-0027-01, and in part by the China Scholarship Council (CSC). (Corresponding author: Minghua Wang.)

Minghua Wang and Qiang Wang are with the Department of Control Science and Engineering, Harbin Institute of Technology, Harbin 150001, China (e-mail: minghuawang1993@163.com; wangqiang@hit.edu.cn).

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

Danfeng Hong is with the Remote Sensing Technology Institute (IMF), German Aerospace Center (DLR), 82234 Weßling, Germany, and also with the GIPSA-Lab, Grenoble INP, CNRS, Université Grenoble Alpes, 38000 Grenoble, France (e-mail: danfeng.hong@dlr.de).

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

Digital Object Identifier 10.1109/TGRS.2021.3055516

## I. INTRODUCTION

HYPERSPECTRAL (HS) imagery acquires hundreds of spectral bands ranging from ultraviolet to infrared wavelength for the same scene on the Earth surface. Compared with other traditional imagery techniques, hyperspectral image (HSI) provides more abundant spatial and spectral information and, thus, has its wide range of applications in agriculture, military, environment monitoring, and food industry [1]–[3]. However, due to various factors of the imaging system and environment, the observed HSIs are unavoidably degraded by noises in real-world scenes. Moreover, in some high-resolution or long-distance transmission HS imaging systems, the increasing amount of 3-D HS data causes the high cost of imaging, storage, or transmission. HSI compressed sensing (CS) is designed to reduce the computation of imaging, and the final recovered HSI images are estimated from a few measurements. As essential preprocessing steps, these HSI restorations have direct influence on the performances of further HSI image processing including classification [4]–[6], unmixing [7], [8], super-resolution [9], [10], feature learning [11], [12], and target detection [13]–[15]. Therefore, it is critical to design effective techniques for HSI restoration.

In general, the optimization problem of HSI restoration methods consists of data fidelity and regularization terms. The former term is derived from a forward model of the object observation. The regularization term can be interpreted as the prior information from underlying properties on recovered HS images. One of the prevalent regularization techniques for image restoration is based on total variation (TV). The TV regularization was introduced to remove noise and preserve image edges for natural image denoising [16]–[18] and CS reconstruction [19]–[21]. Yuan et al. [22] explored local spatial prior and developed a spectral-spatial adaptive TV (SSAHTV) model for HSI denoising. However, SSAHTV just reduces the Gaussian noise and loses useful image detail information. Li et al. [23] minimized the TV of the abundance fractions subject for HSI CS and unmixing. Eason and Andrews [24] proposed an iterative TV method for HSI CS. However, these TV norms strongly enhance spatial piecewise smoothness, which ignores HSI spectral correlation and leads to spatial oversmoothing.

To overcome the drawbacks of the above TV norms, the spectral correlation of 3-D HSIs is taken into account, and
many new TV norms are proposed. Ono et al. [25] proposed an
arranged structure tensor total variation (ASTV) for multichan-
nel image restoration. However, this regularization based on
the singular value decomposition of the low-rank (LR) matrix
spends much computational time. Chang et al. [26] proposed
an anisotropic spectral-spatial total variation (ASSTV) regu-
larization to promote the smoothness across both the spectral
and spatial differences. Unfortunately, ASSTV causes spectral
smoothing due to its strong suppression of direct spectral
difference. Aggarwal and Majumdar [27] combined 2-D TV
along the spatial dimension and 1-D TV along the spec-
tral dimension and proposed a spatiospectral total variation
(SSTV) model. Although SSTV can remove several types of
noises, it produces obvious artifacts in the final recovered HSI
images. Takeyama et al. [28] proposed a hybrid spatiospectral
total variation (HSSTV), which is designed to fuse two kinds
of local differences with a balancing weight. Specifically,
HSSTV1 and HSSTV2 utilize local l₁ spatial norm and
l₂ spatial norm, respectively. However, the authors need to
analyze and determine an extra parameter and, unfortunately,
lose the consideration of HSI global spectral information.

By taking advantage of another HSI prior, i.e., LR, TV
regularization is merged with LR models. LR representation
reveals that the spectral bands of HSI lie in the same LR space
and the HSI data can be represented by a linear combination
of finite endmembers. Zhang et al. [29] rearranged a 3-D
HSI data as a 2-D matrix and proposed an LR matrix
approximation (LRMR) method. Subsequently, He et al.
[30] proposed a band-by-band TV regularized LR matrix
factorization (LRTV) method, and Wu et al. [31] fused the
band-by-band TV regularization into the weighted nuclear
norm minimization (WNMM). Golbabaee and Vandergheynst
[32] proposed a convex function that penalizes both the
trace norm and the band-by-band TV norm of the LR
HSI matrix for HSI CS. However, this band-by-band TV
regularization handles each HSI band separately and, thus,
yields unsatisfying spectral distortions. Wu et al. [33]
incorporated the structure tensor total variation (STV) into
the above WNNM to proposed an STWNNM model. To
exploit spatiospectral total variation (SSTV), Wang et al. [34]
proposed a novel LR constraint and SSTV regularization
model, and Wang et al. [35] integrated this kind of SSTV into
multidirectional LR frameworks. Chen et al. [36] proposed a
weighted group sparsity regularized l₂,₁ norm to exploit the
shared sparse pattern of the difference image for different
bands and combined this regularization into the well-known
LR Tucker decomposition. By using a 3-D anisotropic total
variation (3DATV), He et al. [37] combined it with the local
LR matrix and proposed an LLRSSTV method for HSI
mixed noise removal. Wang et al. [38] merged this 3DATV
into an LR tensor Tucker decomposition model. Peng et al.
[39] calculated the sparsity on the subspace of gradient maps
and proposed an enhanced 3-D total variation (E-3DTV)
regularization for HSI denoising and CS. Wang et al. [40]
propose a novel tensor-based approach named JTenRe3-DTV
for HSI CS, by jointing the tensor Tucker decomposition and
3DTV. Nevertheless, all the above TVs lean upon l₁-norms

and penalize large image gradient magnitudes. In real cases,
HSI edges can be adversely affected.

To effectively sharpen image edges, Li et al. [41] first proposed a new l₀ gradient minimization for natural image smoothing. This minimization globally controlling how many nonzero gradients of results to recover prominent structure have been used in image denoising [41], image deblurring [42], and image inpainting [43]. Xu et al. [42] proposed an effective method based on a generalized l₀ sparse expression for motion deblurring. Xue et al. [43] extended the above l₀ gradient into a low gradient and converted this regularization into an LR framework. Wang et al. [44] developed the l₀ gradient to explore smoothness across both the spatial and spectral domains and combined it with the LR tensor Tucker decomposition model. Xiong et al. [45] embedded the l₀ gradient regularization into an LR tensor factorization. However, these LR-based methods are computationally expensive since they require singular value decomposition. Utilizing either spatial smoothness or spatial-spectral smoothness will inevitably cause oversmoothing or artifacts, respectively.

To consider more prior HSI structure information and intro-
duce less computational burden, we combine spatial–spectral
TV and l₀ gradient regularization. The aforementioned l₀
gradient has a fatal flaw: a user-give parameter penalizing
the degree of smoothness has no physical meaning. Therefore,
Ono [46] proposed a novel l₀ gradient projection employing
a parameter with a clear meaning: the l₀ gradient value of
the output image to control the degree of smoothness. This
l₀ gradient projection variant, l₀ HS TV (l₀hTV), is extended
for HSI restoration and is effective for extracting image edges
information. However, it is difficult for l₀ HS TV to distinguish
outliers, such as sparse noise, when these outliers locate on
HSI image edges. Based on the abovementioned discussion,
we proposed a l₀-l₁ Hybrid Total Variation (l₀-l₁HTV) tech-
nique for HSI mixed noise removal and compressed sensing.
Compared with existing techniques, the main contributions of
this study are summarized as follows.

1) We extend the $l_0$ gradient minimization into HSI restoration, named $l_0hTV$, for directly controlling the number of nonzero image gradients and establishing a physical relationship with the $l_0$ gradient value of the recovered result. $l_0hTV$ characterizes the global spectral-spatial smoothness to complement image details, such as sharp edges and valuable texture.

2) To take full consideration of the local spatial-spectral structure and the global spectral-spatial information of HSI, we design a globally and locally integrated TV regularizer, where $l_0hTV$ is added on SSTV. SSTV captures the local information across spatial and spectral domains to recover the overall structure of HSIs. This TV regularizer compensates mutually to conquer the shortcomings of artifacts or oversmoothing.

3) We integrate SSTV and $l_0hTV$ of HSI into a consolidated framework of $l_0-l_1HTV$, which can be formulated for two HSI tasks: HSI noise removal and compressed sensing. The Alternating Direction Method of Multipliers (ADMM) is adapted to solve the optimization problems.
Experimental results with both simulated and real data demonstrate the superiority of $l_0-l_1$HTV. In particular, the classification accuracy of the results by $l_0-l_1$HTV is improved significantly.

The remainder of this article is organized as follows. Some notations and preliminaries of tensors are introduced in Section II. Section III presents the study of $l_0$ gradient minimization and the new $l_0-l_1$HTV. Then, Sections IV and V provide HSI denoising and compressed sensing models by involving the $l_0-l_1$HTV regularization, respectively. The experimental results and analysis are reported in Section VI. Finally, the conclusions are drawn in Section VII.

## II. NOTATIONS AND PRELIMINARIES

In this article, scalars, vectors, matrices, and tensors are denoted by lowercase letters, e.g., $x$, boldface lowercase letters, e.g., $\mathbf{x}$, boldface capital letters, e.g., $\mathbf{X}$, and Euler script letters, e.g., $\mathcal{X}$, respectively. For a three-way HSI tensor $\mathcal{X} \in \mathbb{R}^{h \times v \times z}$, we denote its $(h_i, v_i, z_i)$-element as $\mathcal{X}_{h_i, v_i, z_i}$ and use notation $\mathcal{X}(i, \cdot; \cdot)$, $\mathcal{X}(\cdot, i; \cdot)$ and $\mathcal{X}(\cdot; \cdot, i)$ to represent the $i$th horizontal, lateral, and frontal slices, respectively. The $l_1$-norm is calculated as $||\mathcal{X}||_1 = \sum_{h_i, v_i, z_i} |\mathcal{X}_{h_i, v_i, z_i}|$, and the Frobenius norm is written as $||\mathcal{X}||_F = (\sum_{h_i, v_i, z_i} |\mathcal{X}_{h_i, v_i, z_i}|^2)^{1/2}$. The $l_1$-norm of a vector denotes $||\mathbf{x}||_1 = \sum_i |\mathbf{x}_i|$, and the $l_2$-norm of a vector is $||\mathbf{x}||_2 = (\sum_i |\mathbf{x}_i|^2)^{1/2}$. A matrix $\mathbf{X} \in \mathbb{R}^{m \times n}$ can be vectorized as $\mathbf{x} = \text{vec}(\mathbf{X}) \in \mathbb{R}^t$, where $t = m * n$ and $\text{vec}(\mathbf{X})$ stacks the columns of $\mathbf{X}$.

**Definition 1 (Mixed $l_{1,0}$ Pseudonorm [46]):** Let $\mathbf{y}$ be a vector of $\mathbb{R}^m$, and index sets $\theta_1, \dots, \theta_i, \dots, \theta_n (1 \le n \le m)$ satisfy as follows.

1) Each $\theta_i$ is a subset of $1, \dots, m$.

2) $\theta_i \cap \theta_j = \emptyset$ for any $i \neq l$.

3) $\bigcup_{i=1}^n \theta_i = 1, \dots, m$.

The mixed $l_{1,0}$ pseudonorm of $\mathbf{y}$ is defined as follows:

$$ ||\mathbf{y}||_{1,0}^{\theta} = ||(|\mathbf{y}_{\theta_1}||_1, \dots, |\mathbf{y}_{\theta_i}||_1, \dots, |\mathbf{y}_{\theta_n}||_1)||_0 \quad (1) $$

where $\mathbf{y}_{\theta_i}$ denotes a subvector of $\mathbf{y}$ with its entries specified by $\theta_i$.

**Definition 2 (Indicator Function I(y) of $l_{1,0}$ Mixed Pseudonorm [46]):**

$$ I_{||\mathcal{B}||_{1,0}^{\theta}}(\mathbf{y}) = \begin{cases} 0, & ||\mathcal{B}\mathbf{y}||_{1,0}^{\theta} \le \gamma \\ \infty, & \text{otherwise} \end{cases} \quad (2) $$

where $\mathcal{B}$ is a given operator. In this article, $\mathcal{B}$ is a diagonal tensor with 0 or 1 entries.

**Proposition 1 (Projection Onto $l_{1,0}$ Mixed Pseudonorm Ball With Binary Mask [46]):** Set $\mathbf{y} \in \mathbb{R}^m$ as a known vector and set $\gamma$ as a nonnegative integer. Let $\mathbf{W}$ be a known diagonal binary matrix, and let $\theta_1, \dots, \theta_i, \dots, \theta_n (1 \le n \le m)$ be index sets satisfying the conditions from **definition 1**. Without loss of generality, $\mathbf{W}\mathbf{y} = (\mathbf{y}_{\theta_1}^T, \dots, \mathbf{y}_{\theta_n}^T)^T$ are assumed. $\mathbf{y}_{\theta_1}, \dots, \mathbf{y}_{\theta_n}$ are sorted in descending order according to $||\mathbf{y}_{\theta_{-1}}||_2 > ||\mathbf{y}_{\theta_{-2}}||_2 > \dots > ||\mathbf{y}_{\theta_{-n}}||_2$, where $\mathbf{y}_{\theta_{-1}}, \dots, \mathbf{y}_{\theta_{-n}}$ are obtained with the new order, and the original index sets $i$ have $\vec{i}$ that corresponds to it one by one.

For the following problem:

$$ \min_{\mathbf{y}^*} ||\mathbf{y}^* - \mathbf{y}||^2 \quad \text{s.t.} \quad ||\mathbf{W}\mathbf{y}^*||_{1,0}^{\theta} \le \gamma \quad (3) $$

one of the optimal solutions is given by

$$ \mathbf{y}^* = \begin{cases} \mathbf{y}, & ||\mathbf{W}\mathbf{y}||_{1,0}^{\theta} \le \gamma \\ (\tilde{\mathbf{y}}_{\theta_{-1}}^{\mathrm{T}}, \dots, \tilde{\mathbf{y}}_{\theta_{-n}}^{\mathrm{T}})^{\mathrm{T}} + (\mathbf{I} - \mathbf{W})\mathbf{y}, & ||\mathbf{W}\mathbf{y}||_{1,0}^{\theta} > \gamma \end{cases} \quad (4) $$

where

$$ \tilde{\mathbf{y}}_{\theta_{-1}}^{\mathrm{T}} = \begin{cases} \mathbf{y}_{\theta_1}^{\mathrm{T}}, & \text{if } i \in \{1, \dots, \gamma\} \\ 0, & \text{if } i \in \{\gamma + 1, \dots, n\}. \end{cases} \quad (5) $$

**Definition 3 (Mode-$k_1, k_2$ Unfolding of a Tensor [47]):** For a $m$-way tensor $\mathcal{X} \in \mathbb{R}^{n_1 \times n_2 \times \dots \times n_m}$, this operator is denoted by $\mathcal{X}(k_1, k_2) \in \mathbb{R}^{n_{k_1} \times n_{k_2} \times \prod_{i \neq k_1, k_2} n_i}$, whose frontal slices are the lexicographic ordering of the mode-$k_1, k_2$ slice of $\mathcal{X}$.

## III. $l_0-l_1$ HYBRID TOTAL VARIATION

### A. $l_0$ Hyperspectral TV

Xu et al. [42] first proposed the $l_0$ gradient model to preserve natural image edges. For a natural 2-D image $\mathbf{X} \in \mathbb{R}^{m \times n}$, the $l_0$ Gradient is defined as follows [42]:

$$ l_0 \text{TV}(\mathbf{X}) = \sum_{i=1}^{m*n} C(||\text{vec}(D_h \mathbf{X})||_1 + ||\text{vec}(D_v \mathbf{X})||_1) \\ = \sum_{i=1}^{m} \sum_{j=1}^{n} C(|\mathbf{X}_{i+1,j} - \mathbf{X}_{i,j}| + |\mathbf{X}_{i,j+1} - \mathbf{X}_{i,j}|) \quad (6) $$

where $D_h$ and $D_v$ are linear difference operators along image horizontal and vertical directions. $C(\mathbf{X})$ is a binary function counting the number of nonzero image gradients and $C(\mathbf{X}_{i,j}) = 1$, if $\mathbf{X}_{i,j} \neq 0$; $C(\mathbf{X}_{i,j}) = 0$, otherwise.

For an HSI $\mathcal{X} \in \mathbb{R}^{h \times d \times z}$, $\sum_{k=1}^z l_0 \text{TV}(\mathcal{X}^k)$ can be a band-by-band model by the simplest way to extend. Although the spatial smoothness of each HSI band is utilized separately, HSI spectral correlation is neglected. Simultaneously, different HSI bands are acquired from the same scene, and the spatial information of all spectral bands should be fused and formulated. $l_0$hTV is defined as follows:

$$ l_0 \text{hTV}(\mathcal{X}) = \sum_{i}^{h} \sum_{j}^{v} C \left( \sum_{k}^{z} \left( |\mathcal{X}_{i+1,j,k} - \mathcal{X}_{i,j,k}| + |\mathcal{X}_{i,j+1,k} - \mathcal{X}_{i,j,k}| \right) \right) \quad (7) $$

where boundary values of gradients are written as $\mathcal{X}_{i,j+1,k} - \mathcal{X}_{i,j,k} = 0$, if $j = v$, and $\mathcal{X}_{i+1,j,k} - \mathcal{X}_{i,j,k} = 0$, if $i = h$. $l_0$hTV($\mathcal{X}$), actually counting the nonzero gradients in the spatial dimension with the assistance of all spectral bands, greatly benefits for capturing the spatial edges of the HSI.
<figure><img src="image_2.png" /><figcaption>Fig. 1. Illustration of the proposed $l_0-l_1$ hybrid total variation regularization.</figcaption></figure>

For the mixed $l_{1,0}$ pseudonorm defined by Definition 1, another formulation is adopted to express $l_0hTV$ as follows:

$$l_0hTV(\mathcal{X}) = ||BD\mathcal{X}||_{1,0}^{\theta} \quad (8)$$

where operator $D$ is an operator to fuse the differences along the horizontal and vertical directions into one tensor. Operator $B$ is an operator that forces boundary values of gradients to be zero when $i = h$ and $j = v$. A simple example is given as follows. $D_h\mathcal{X}$ and $D_v\mathcal{X}$ denote the difference along image horizontal and vertical directions, respectively. $D_h\mathcal{X}$ and $D_v\mathcal{X}$ are stored into one tensor $D\mathcal{X}$, i.e., $D\mathcal{X}(:, : , 1, : ) = D_h\mathcal{X}$ and $D\mathcal{X}(:, : , 2, : ) = D_v\mathcal{X}$. The boundary values of $BD\mathcal{X}(:, : , 1, 1)$, $BD\mathcal{X}(:, : , 1, 2)$, and $BD\mathcal{X}(:, : , 2, 1)$, $BDX(:, : , 2, 2)$ are set as 0 when $i = h = 3$ and $j = v = 3$, respectively.

**Example 1:** A tensor $\mathcal{X} \in \mathbb{R}^{3 \times 3 \times 2}$ and its $BD\mathcal{X} \in \mathbb{R}^{3 \times 3 \times 2 \times 2}$ can be expressed as follows:

$$\mathcal{X}(:, : , 1) = \begin{bmatrix} \mathcal{X}_{111} & \mathcal{X}_{121} & \mathcal{X}_{131} \\ \mathcal{X}_{211} & \mathcal{X}_{221} & \mathcal{X}_{231} \\ \mathcal{X}_{311} & \mathcal{X}_{321} & \mathcal{X}_{331} \end{bmatrix} \quad (9)$$

$$\mathcal{X}(:, : , 2) = \begin{bmatrix} \mathcal{X}_{112} & \mathcal{X}_{122} & \mathcal{X}_{132} \\ \mathcal{X}_{212} & \mathcal{X}_{222} & \mathcal{X}_{232} \\ \mathcal{X}_{312} & \mathcal{X}_{322} & \mathcal{X}_{332} \end{bmatrix} \quad (10)$$

$$(BD\mathcal{X})(:, : , 1, 1) = \begin{bmatrix} \mathcal{X}_{121} - \mathcal{X}_{111} & \mathcal{X}_{131} - \mathcal{X}_{121} & 0 \\ \mathcal{X}_{221} - \mathcal{X}_{211} & \mathcal{X}_{231} - \mathcal{X}_{221} & 0 \\ \mathcal{X}_{321} - \mathcal{X}_{311} & \mathcal{X}_{331} - \mathcal{X}_{321} & 0 \end{bmatrix} \quad (11)$$

$$(BD\mathcal{X})(:, : , 1, 2) = \begin{bmatrix} \mathcal{X}_{122} - \mathcal{X}_{112} & \mathcal{X}_{132} - \mathcal{X}_{122} & 0 \\ \mathcal{X}_{222} - \mathcal{X}_{212} & \mathcal{X}_{232} - \mathcal{X}_{222} & 0 \\ \mathcal{X}_{322} - \mathcal{X}_{312} & \mathcal{X}_{332} - \mathcal{X}_{322} & 0 \end{bmatrix} \quad (12)$$

$$(BD\mathcal{X})(:, : , 2, 1) = \begin{bmatrix} \mathcal{X}_{211} - \mathcal{X}_{111} & \mathcal{X}_{221} - \mathcal{X}_{121} & \mathcal{X}_{231} - \mathcal{X}_{131} & 0 \\ \mathcal{X}_{311} - \mathcal{X}_{211} & \mathcal{X}_{321} - \mathcal{X}_{221} & \mathcal{X}_{331} - \mathcal{X}_{231} & 0 \\ 0 & 0 & 0 & 0 \end{bmatrix} \quad (13)$$

$$(BD\mathcal{X})(:, : , 2, 2) = \begin{bmatrix} \mathcal{X}_{212} - \mathcal{X}_{112} & \mathcal{X}_{222} - \mathcal{X}_{122} & \mathcal{X}_{232} - \mathcal{X}_{132} \\ \mathcal{X}_{312} - \mathcal{X}_{212} & \mathcal{X}_{322} - \mathcal{X}_{222} & \mathcal{X}_{332} - \mathcal{X}_{232} \\ 0 & 0 & 0 \end{bmatrix}. \quad (14)$$

With Definition 1, here, $n = 9$, and each $y_{\theta_i}$ is determined by mode-1,2 unfolding $(BD\mathcal{X})_{(1,2)} \in \mathbb{R}^{9 \times 2 \times 2}$, for example, $y_{\theta_1} = [\mathcal{X}_{121} - \mathcal{X}_{111}, \mathcal{X}_{122} - \mathcal{X}_{112}, \mathcal{X}_{211} - \mathcal{X}_{111}, \mathcal{X}_{212} - \mathcal{X}_{112}]^T$. It is easily proven that $l_0hTV(\mathcal{X})$ and $||BD\mathcal{X}||_{1,0}^{\theta}$ are mathematically equivalent.

## B. $l_0-l_1$ Hybrid Total Variation Regularization

Fig. 1 represents the visual illustration of the $l_0-l_1$HTV regularization. As shown at the top of Fig. 1, the global spatial smoothness of all spectral bands is restrained by $l_0hTV(\mathcal{X}) = ||BD\mathcal{X}||_{1,0}^{\theta}$, where $\gamma$ denotes the $l_0$ gradient value of the output image. As the number of iterations increases, the image edges of HSI gradually become clearer. Meanwhile, if $l_0$-norm is excessively penalized, the final results will appear oversmoothing phenomenon. At the bottom of Fig. 1, the SSTV handling local spatial-spectral differences combine 2-D spatial TV (defined by $2DTV = ||D_h\mathcal{X}||_1 + ||D_v\mathcal{X}||_1$) with 1-D spectral TV. Although most of structure information is restored, SSTV unavoidably caused artifacts. For more comprehensively exploiting the global spatial and the local spatial-spectral information, $l_0hTV(\mathcal{X})$ and SSTV are incorporated as a globally and locally integrated regularizer, and with **Definition 2**, the proposed $l_0-l_1$ Hybrid Total Variation ($l_0-l_1$HTV) is defined as follows:

$$\begin{aligned} l_0-l_1HTV(\mathcal{X}) &= l_0hTV(\mathcal{X}) + SSTV(\mathcal{X}) \\ &= I_{||B \cdot ||_{1,0}^{\theta}}(BD\mathcal{X}) + (||D_h\mathcal{X}D_z||_1 + ||D_v\mathcal{X}D_z||_1) \end{aligned} \quad (15)$$
where $D_h$ and $D_v$ are two operators to calculate differences along the horizontal and vertical directions. $D_z$ is a 1-D finite differencing operator to extract the spectral signature of each pixel. In practical applications, it is hard to find a proper parameter $\lambda$ for minimizing $\lambda l_0 h$TV($\mathcal{X}$) since this user-given parameter $\lambda$ has no relationship with the degree of the restored image $\mathcal{X}$. Therefore, $||BD\mathcal{X}||_{1,0} \le \gamma$ is used for HSI restoration, and $\gamma$ directly represents the $l_0$ gradient value of the output image.

#### IV. $l_0-l_1$ HYBRID TOTAL VARIATION FOR HSI DENOISING

A clear HSI cube $\mathcal{X} \in \mathbb{R}^{h \times v \times z}$ is corrupted by mixed noise, including the Gaussian noise, salt and pepper noise, deadlines, and strips. The Gaussian noise can be denoted by $\mathcal{N}$, and the others are unified as sparse noise $\mathcal{S}$. Thus, an observed 3D HSI $\mathcal{T} \in \mathbb{R}^{h \times v \times z}$ is represented as follows:

$$ \mathcal{T} = \mathcal{X} + \mathcal{S} + \mathcal{N}. \quad (16) $$

With (15) in mind, the proposed $l_0-l_1$HTV for removing HSI mixed noise is formulated as

$$ \min_{\mathcal{X}, \mathcal{S}} \|\mathcal{T} - \mathcal{X} - \mathcal{S}\|_{F}^{2} + \lambda \|\mathcal{S}\|_{1} + \lambda_{tv} (\|\mathcal{P}\|_{1} + \|\mathcal{Q}\|_{1}) \\ \times (\|D_h \mathcal{X} D_z\|_{1} + \|D_v \mathcal{X} D_z\|_{1}) + I_{\|B\|_{1,0}^{\theta}} (BD\mathcal{X}) \quad (17) $$

where $\lambda$ and $\lambda_{tv}$ are two positive constants to balance the contributions of each term. The first term $\|\mathcal{T} - \mathcal{X} - \mathcal{S}\|_{F}^{2}$ is a data fidelity term. $\|\mathcal{S}\|_{1}$ is the $l_1$-norm of the sparse noise $\mathcal{S}$ to remove salt and pepper noise, deadlines, and strips. $\|D_h \mathcal{X} D_z\|_{1} + \|D_v \mathcal{X} D_z\|_{1}$ is used to characterize the local smoothness in the spatial and spectral directions. $I_{\|B\|_{1,0}^{\theta}} (BD\mathcal{X})$ is the mixed $l_{1,0}$ pseudonorm to enhance image edges and remove the Gaussian noise.

To solve the above HSI denoising problem, three auxiliary variable $\mathcal{V}$, $\mathcal{P}$, and $\mathcal{Q}$ are introduced, and (17) is reformulated as follows:

$$ \min_{\mathcal{X}, \mathcal{S}, \mathcal{V}, \mathcal{P}, \mathcal{Q}} \|\mathcal{T} - \mathcal{X} - \mathcal{S}\|_{F}^{2} + \lambda \|\mathcal{S}\|_{1} + \lambda_{tv} (\|\mathcal{P}\|_{1} + \|\mathcal{Q}\|_{1}) \\ + I_{\|B\|_{1,0}^{\theta}} (\mathcal{V}) \\ \text{s.t. } \mathcal{P} = D_h \mathcal{X} D_z, \mathcal{Q} = D_v \mathcal{X} D_z, \mathcal{V} = D\mathcal{X}. \quad (18) $$

The ADMM method [48], [49] is applied for the aforementioned problem, and the unconstrained function of (18) is rewritten as:

$$ \min_{\mathcal{X}, \mathcal{S}, \mathcal{V}, \mathcal{P}, \mathcal{Q}} \|\mathcal{T} - \mathcal{X} - \mathcal{S}\|_{F}^{2} + \lambda \|\mathcal{S}\|_{1} + \lambda_{tv} (\|\mathcal{P}\|_{1} + \|\mathcal{Q}\|_{1}) \\ + I_{\|B\|_{1,0}^{\theta}} (\mathcal{V}) + \beta \|\mathcal{P} - D_h \mathcal{X} D_z - \mathcal{Y}_1\|_{F}^{2} \\ + \beta \|\mathcal{Q} - D_v \mathcal{X} D_z - \mathcal{Y}_2\|_{F}^{2} + \beta \|\mathcal{V} - D\mathcal{X} - \mathcal{Y}_3\|_{F}^{2}. \quad (19) $$

Equation (19) is divided into five subproblems:

$$ \begin{aligned} & 1) \min_{\mathcal{X}} \|\mathcal{T} - \mathcal{X} - \mathcal{S}\|_{F}^{2} + \beta \|\mathcal{P} - D_h \mathcal{X} D_z - \mathcal{Y}_1\|_{F}^{2} \\ & \qquad + \beta \|\mathcal{Q} - D_v \mathcal{X} D_z - \mathcal{Y}_2\|_{F}^{2} + \beta \|\mathcal{V} - D\mathcal{X} - \mathcal{Y}_3\|_{F}^{2} \quad (20) \end{aligned} $$

$$ \begin{aligned} & 2) \min_{\mathcal{S}} \lambda \|\mathcal{S}\|_{1} + \|\mathcal{T} - \mathcal{X} - \mathcal{S}\|_{F}^{2} \quad (21) \\ & 3) \min_{\mathcal{P}} \lambda_{tv} \|\mathcal{P}\|_{1} + \beta \|\mathcal{P} - D_h \mathcal{X} D_z - \mathcal{Y}_1\|_{F}^{2} \quad (22) \end{aligned} $$

$$ \begin{aligned} & 4) \min_{\mathcal{Q}} \lambda_{tv} \|\mathcal{Q}\|_{1} + \beta \|\mathcal{Q} - D_v \mathcal{X} D_z - \mathcal{Y}_2\|_{F}^{2} \quad (23) \\ & 5) \min_{\mathcal{V}} I_{\|B\|_{1,0}^{\theta}} (\mathcal{V}) + \beta \|\mathcal{V} - D\mathcal{X} - \mathcal{Y}_3\|_{F}^{2}. \quad (24) \end{aligned} $$

$$ \begin{aligned} & 4) \min_{\mathcal{Q}} \lambda_{tv} \|\mathcal{Q}\|_{1} + \beta \|\mathcal{Q} - D_v \mathcal{X} D_z - \mathcal{Y}_2\|_{F}^{2} \quad (23) \\ & 5) \min_{\mathcal{V}} I_{\|B\|_{1,0}^{\theta}} (\mathcal{V}) + \beta \|\mathcal{V} - D\mathcal{X} - \mathcal{Y}_3\|_{F}^{2}. \quad (24) \end{aligned} $$

### Algorithm 1 $l_0-l_1$HTV for HSI Mixed Noise Removal

**Input:** $\mathcal{T}$, $\lambda$, $\lambda_{tv}$, $\gamma$

1: **Initialization:** $\mathcal{X} = \mathcal{S} = 0, \mathcal{V} = \mathcal{P} = \mathcal{Q} = \mathcal{Y}_1 = \mathcal{Y}_2 = \mathcal{Y}_3 = 0, i = 0, \epsilon = 10^{-4} * h\nu$

2: **while** $0 \le i \le \text{maxiter}$ or $||BD\mathcal{X}||_{1,0}^{\theta} - \gamma| > \epsilon$

3: $i = i + 1$;

4: update $\mathbf{x}$ by solving (26):

$$ (1 + \beta D^{*} D + \beta D_{1}^{*} D_{1} + \beta D_{2}^{*} D_{2}) \mathbf{x} = \mathbf{t} - \mathbf{s} $$

5: $+ \beta(\mathbf{p} - \mathbf{y}_1 + \mathbf{q} - \mathbf{y}_2 + \mathbf{v} - \mathbf{y}_3)$

6: update $\mathcal{X}$: $\mathcal{X}$=rearrange($\mathbf{x}$)

7: update $\mathcal{S}$, $\mathcal{P}$ and $\mathcal{Q}$ by solving (30), (31) and (32):

8: $\mathcal{S} = S_{\lambda/2}[\mathcal{T} - \mathcal{X}]$

9: $\mathcal{P} = S_{\lambda_{tv}/2}[D_h \mathcal{X} D_z + \mathcal{Y}_1]$

10: $\mathcal{Q} = S_{\lambda_{tv}/2}[D_v \mathcal{X} D_z + \mathcal{Y}_2]$

11: update $\mathcal{V}$ by solving (33):

12: $\mathcal{V}' = B(D\mathcal{X} + \mathcal{Y}_3)$,

13: Sort $\mathcal{V}'_{\theta_1}, \dots, \mathcal{V}'_{\theta_n}$ in descending order in terms of $l_2$ norms,

14: Set $\mathcal{V}'_{\theta_{j+1}} = 0, \dots, \mathcal{V}'_{\theta_n} = 0$ and then obtain a updated $\mathcal{V}'$,

15: $\mathcal{V} = \mathcal{V}' + (I - B)(D\mathcal{X} + \mathcal{Y}_3)$

16: update the multipliers

$$ \mathcal{Y}_1 = \mathcal{Y}_1 + (D_h \mathcal{X} D_z - \mathcal{P}) $$

$$ \mathcal{Y}_2 = \mathcal{Y}_2 + (D_h \mathcal{X} D_z - \mathcal{Q}) $$

$$ \mathcal{Y}_3 = \mathcal{Y}_3 + (D\mathcal{X} - \mathcal{V}) $$

17: **end while**

**Output:** $\mathcal{X}$.

The $\mathcal{X}$-related subproblem is rewritten into its vector

$$ \begin{aligned} \mathbf{x} &= \arg \min_{\mathbf{x}} f(\mathbf{x}) \\ &= \arg \min_{\mathbf{x}} \|\mathbf{t} - \mathbf{x} - \mathbf{s}\|_{F}^{2} + \beta \|\mathbf{p} - D_{1} \mathbf{x} - \mathbf{y}_{1}\|_{F}^{2} \\ &\quad + \beta \|\mathbf{q} - D_{2} \mathbf{x} - \mathbf{y}_{2}\|_{F}^{2} + \beta \|\mathbf{v} - D \mathbf{x} - \mathbf{y}_{3}\|_{F}^{2} \quad (25) \end{aligned} $$

where $D_1 = D_h \otimes D_z$, $D_2 = D_v \otimes D_z$, and $\otimes$ denotes the Kronecker product. Set the partial derivative of $f(\mathbf{x})$ as zero, and yield the following linear equation:

$$ (1 + \beta D^{*} D + \beta D_{1}^{*} D_{1} + \beta D_{2}^{*} D_{2}) \mathbf{x} = \mathbf{t} - \mathbf{s} + \beta (\mathbf{p} - \mathbf{y}_{1} + \mathbf{q} - \mathbf{y}_{2} + \mathbf{v} - \mathbf{y}_{3}). \quad (26) $$

Equation (26) can be solved by the preconditioned conjugate gradient (PCG) [50] or the algorithm for sparse linear equations and least squares (LSQR) [51], [52].

Equations (21), (22), and (23) about $\mathcal{S}$, $\mathcal{P}$, and $\mathcal{Q}$ are the same form and, respectively, rewritten as

$$ \min_{\mathcal{S}} \lambda \|\mathcal{S}\|_{1} + \|\mathcal{S} - (\mathcal{T} - \mathcal{X})\|_{F}^{2} \quad (27) $$

$$ \min_{\mathcal{P}} \lambda_{tv} \|\mathcal{P}\|_{1} + \beta \|\mathcal{P} - (D_h \mathcal{X} D_z + \mathcal{Y}_1)\|_{F}^{2} \quad (28) $$

$$ \min_{\mathcal{Q}} \lambda_{tv} \|\mathcal{Q}\|_{1} + \beta \|\mathcal{Q} - (D_v \mathcal{X} D_z + \mathcal{Y}_2)\|_{F}^{2} \quad (29) $$

where the sign operator $\operatorname{sgn}(x)$ is defined as follows:

$$ \operatorname{sgn}(x) = \begin{cases} 1, & x > 0 \\ 0, & x = 0 \\ -1, & x < 0. \end{cases} $$
By using the operator $S_{\varepsilon}[x] = \operatorname{sgn}(x) \cdot \max(|x| - \varepsilon, 0)$, the solutions are obtained as

$$S = S_{\lambda/2}[T - \mathcal{X}] \quad (30)$$

$$\mathcal{P} = S_{\lambda_{tv}/2}[D_h \mathcal{X}D_z + \mathcal{Y}_1] \quad (31)$$

$$\mathcal{Q} = S_{\lambda_{tv}/2}[D_v \mathcal{X}D_z + \mathcal{Y}_2]. \quad (32)$$

The subproblem of $\mathcal{V}$ is equivalent to the following constrained minimization problem:

$$\begin{array}{ll} \min_{\mathcal{V}} & ||\mathcal{V} - D\mathcal{X} - \mathcal{Y}_3||_F^2 \\ \text{s.t.} & ||B\mathcal{V}||_{1,0}^\theta \le \gamma. \end{array} \quad (33)$$

The parameter $\beta$ is omitted because it has no impact on the optimization. Proposition 1 is adopted to solve (33). First, $\mathcal{V}'$ is initialized to $\mathcal{V}' = B(D\mathcal{X} + \mathcal{Y}_3)$, and at the $i$th iteration, $\mathcal{V}^{(i)}$ is sorted in descending order. $\mathcal{V}'_{\gamma+1} = 0, \dots, \mathcal{V}'_{\theta_n}$ are set as zero vectors, which effectively removes most of noise while preserving sharper image edges. Simultaneously, artifacts of the denoised result are avoided since the information from all spectral bands is used to promote spatial smoothness. The detail of the $l_0-l_1$HTV method for HSI mixed noise removal is summarized in Algorithm 1.

## V. $l_0-l_1$ HYBRID TOTAL VARIATION FOR HSI COMPRESSED SENSING

HSI compressed sensing aims at reconstructing a 3-D HSI data $\mathcal{X}$ from compressive measurements $\mathbf{y} \in \mathbb{R}^m$. In real acquisition, disturbed by the Gaussian noise $\mathbf{n}$, the compressive measurements $\mathbf{y}$ can be obtained by

$$\mathbf{y} = \Psi(\mathcal{X}) + \mathbf{n}. \quad (34)$$

$\Psi$ is instantiated as $\Psi = \mathbf{D} \cdot \mathbf{H} \cdot \mathbf{P}$, where $\mathbf{D}$ is a random down sampling operator, $\mathbf{H}$ is a random permutation matrix, $\mathbf{P}$ is the Walsh-Hadamard transform, and the mapping of $\Psi$ is $\mathbb{R}^{h \times d \times z} \rightarrow \mathbb{R}^m$. This compressed operator has been successfully applied for various CS cases [32], [53], [54]. However, it is an ill-posed inverse problem to recover $\mathcal{X}$ from (34). As the prior information is efficient to find the optimal solution, $l_0-l_1$HTV regularization is adopted, and the proposed HSI CS model is formulated as

$$\min_{\mathcal{X}} \|\mathbf{y} - \Psi(\mathcal{X})\|_F^2 + \lambda_{tv} (\|D_h \mathcal{X} D_z\|_1 + \|D_v \mathcal{X} D_z\|_1) + I_{||\mathbf{B} \cdot \mathbf{1}_{1,0}^{\theta}(\mathcal{V})} (\mathbf{B} \mathcal{D} \mathcal{X}). \quad (35)$$

By introducing three auxiliary variables $\mathcal{P}$, $\mathcal{Q}$ and $\mathcal{V}$, we rewrite the (35) as follows:

$$\begin{array}{ll} \min_{\mathcal{X}, \mathcal{V}, \mathcal{P}, \mathcal{Q}} & \|\mathbf{y} - \Psi(\mathcal{X})\|_F^2 + \lambda_{tv} (\|\mathcal{P}\|_1 + \|\mathcal{Q}\|_1) + I_{||\mathbf{B} \cdot \mathbf{1}_{1,0}^{\theta}(\mathcal{V})} \\ \text{s.t.} & \mathcal{P} = D_h \mathcal{X} D_z, \mathcal{Q} = D_v \mathcal{X} D_z, \mathcal{V} = D \mathcal{X}. \end{array} \quad (36)$$

Based on the ADMM method [48], the unconstrained function of (36) is given as follows:

$$\begin{array}{ll} \min_{\mathcal{X}, \mathcal{V}, \mathcal{P}, \mathcal{Q}} & \|\mathbf{y} - \Psi(\mathcal{X})\|_F^2 + \lambda_{tv} (\|\mathcal{P}\|_1 + \|\mathcal{Q}\|_1) + I_{||\mathbf{B} \cdot \mathbf{1}_{1,0}^{\theta}(\mathcal{V})} \\ & + \beta \|\mathcal{P} - D_h \mathcal{X} D_z - \mathcal{Y}_1\|_F^2 + \beta \|\mathcal{Q} - D_v \mathcal{X} D_z - \mathcal{Y}_2\|_F^2 \\ & + \beta \|\mathcal{V} - D \mathcal{X} - \mathcal{Y}_3\|_F^2. \end{array} \quad (37)$$

Equation (37) is divided into four subproblems:

$$\begin{array}{ll} 1) & \min_{\mathcal{X}} \|\mathbf{y} - \Psi(\mathcal{X})\|_F^2 + \beta \|\mathcal{P} - D_h \mathcal{X} D_z - \mathcal{Y}_1\|_F^2 \\ & + \beta \|\mathcal{Q} - D_v \mathcal{X} D_z - \mathcal{Y}_2\|_F^2 + \beta \|\mathcal{V} - D \mathcal{X} - \mathcal{Y}_3\|_F^2. \end{array} \quad (38)$$

$$2) \min_{\mathcal{P}} \lambda_{tv} \|\mathcal{P}\|_1 + \beta \|\mathcal{P} - D_h \mathcal{X} D_z - \mathcal{Y}_1\|_F^2 \quad (39)$$

$$3) \min_{\mathcal{Q}} \lambda_{tv} \|\mathcal{Q}\|_1 + \beta \|\mathcal{Q} - D_v \mathcal{X} D_z - \mathcal{Y}_2\|_F^2. \quad (40)$$

$$4) \min_{\mathcal{V}} I_{||\mathbf{B} \cdot \mathbf{1}_{1,0}^{\theta}(\mathcal{V})} (\mathcal{V}) + \beta \|\mathcal{V} - D \mathcal{X} - \mathcal{Y}_3\|_F^2. \quad (41)$$

## Algorithm 2 $l_0-l_1$HTV for HSI Compressed Sensing

**Input:** $\mathbf{y}$, $\lambda_{tv}$, $\gamma$

1. **Initialization:** $\mathcal{X} = 0$, $\mathcal{V} = \mathcal{P} = \mathcal{Q} = \mathcal{Y}_1 = \mathcal{Y}_2 = \mathcal{Y}_3 = 0$, $i = 0, \epsilon = 10^{-4} * h\nu$

2. **while** $0 \le i \le \maxiter$ or $||\mathcal{B} \mathcal{D} \mathcal{X}||_{1,0}^{\theta} - \gamma| > \epsilon$

3. $i = i + 1$;

4. update $\mathbf{x}$ by solving (43):

$$(1 + \beta D^* D + \beta D_1^* D_1 + \beta D_2^* D_2)\mathbf{x} = \mathbf{y} - \Psi(\mathbf{x})$$

5. $+ \beta(\mathbf{p} - \mathbf{y}_1 + \mathbf{q} - \mathbf{y}_2 + \mathbf{v} - \mathbf{y}_3)$

6. update $\mathcal{X}$: $\mathcal{X}$=rearrange($\mathbf{x}$)

7. update $\mathcal{P}$ and $\mathcal{Q}$ by solving (44) and (45):

8. $\mathcal{P} = S_{\lambda_{tv}/2}[D_h \mathcal{X} D_z + \mathcal{Y}_1]$

9. $\mathcal{Q} = S_{\lambda_{tv}/2}[D_v \mathcal{X} D_z + \mathcal{Y}_2]$

10. update $\mathcal{V}$ by solving (33):

11. $\mathcal{V}' = B(D\mathcal{X} + \mathcal{Y}_3)$,

12. Sort $\mathcal{V}'_{\theta_1}, \dots, \mathcal{V}'_{\theta_n}$ in descending order in terms of $l_2$ norms,

13. Set $\mathcal{V}'_{\theta_{\gamma+1}} = 0, \dots, \mathcal{V}'_{\theta_n} = 0$ and then obtain a updated $\mathcal{V}'$,

14. $\mathcal{V}' = \mathcal{V}' + (I - B)(D\mathcal{X} + \mathcal{Y}_3)$

15. update the multipliers

$$\mathcal{Y}_1 = \mathcal{Y}_1 + (D_h \mathcal{X} D_z - \mathcal{P})$$

$$\mathcal{Y}_2 = \mathcal{Y}_2 + (D_h \mathcal{X} D_z - \mathcal{Q})$$

$$\mathcal{Y}_3 = \mathcal{Y}_3 + (D \mathcal{X} - \mathcal{V})$$

16. **end while**

**Output:** $\mathcal{X}$.

The $\mathcal{X}$-related subproblem is rewritten into its vector

$$\begin{align} \mathbf{x} &= \operatorname*{arg\,min}_{\mathbf{x}} f(\mathbf{x}) \nonumber \\ &= \operatorname*{arg\,min}_{\mathbf{x}} \|\mathbf{y} - \Psi(\mathbf{x})\|_F^2 + \beta \|\mathbf{p} - D_1 \mathbf{x} - \mathbf{y}_1\|_F^2 \nonumber \\ &\quad + \beta \|\mathbf{q} - D_2 \mathbf{x} - \mathbf{y}_2\|_F^2 + \beta \|\mathbf{v} - D \mathbf{x} - \mathbf{y}_3\|_F^2. \tag{42} \end{align}$$

The partial derivative of $f(\mathbf{x})$ is set as zero, and optimizing $\mathbf{x}$ is equivalent to solve the following linear equation:

$$(1 + \beta D^* D + \beta D_1^* D_1 + \beta D_2^* D_2)\mathbf{x} = \mathbf{y} - \Psi(\mathbf{x}) + \beta(\mathbf{p} - \mathbf{y}_1 + \mathbf{q} - \mathbf{y}_2 + \mathbf{v} - \mathbf{y}_3). \quad (43)$$

Equation (43) is also solved by the PCG [50] or the algorithm for sparse linear equations and least squares (LSQR) [52]. Based on the operator $S_{\varepsilon}[x] = \operatorname{sgn}(x) \cdot \max(|x| - \varepsilon, 0)$, the solutions of (39) and (40) are calculated as

$$\mathcal{P} = S_{\lambda_{tv}/2}[D_h \mathcal{X} D_z + \mathcal{Y}_1] \quad (44)$$

$$\mathcal{Q} = S_{\lambda_{tv}/2}[D_v \mathcal{X} D_z + \mathcal{Y}_2]. \quad (45)$$
<table><caption>TABLE I
ADDED MIXED NOISES IN NINE CASES</caption><thead><tr><th>Noise Case</th><th>Mixed Noise</th></tr></thead><tbody><tr><td>Case 1</td><td>Gaussian (SNR = 10 ~ 20 dB)</td></tr><tr><td>Case 2</td><td>Gaussian (SNR = 10 ~ 20 dB)+salt and pepper(10%) in all bands</td></tr><tr><td>Case 3</td><td>Gaussian (SNR = 10 ~ 20 dB)+salt and pepper(20%) in all bands + four vertical and horizontal deadlines in each of four random bands</td></tr><tr><td>Case 4</td><td>Gaussian (mean SNR = 12.12 dB, and noise variance is calculated by the Eq.(50))</td></tr><tr><td>Case 5</td><td>Gaussian (mean SNR = 12.12 dB, and noise variance is calculated by the Eq.(50)) +salt and pepper(10%) in all bands</td></tr><tr><td>Case 6</td><td>Gaussian (mean SNR = 12.12 dB, and noise variance is calculated by the Eq.(50)) +salt and pepper(20%) in all bands + eight vertical and horizontal deadlines in each of eight random bands</td></tr><tr><td>Case 7</td><td>Gaussian (zero-mean, and noise variance is different for different bands and is selected from 0.1~0.2 randomly)</td></tr><tr><td>Case 8</td><td>Gaussian (zero-mean, and noise variance is different for different bands and is selected from 0.1~0.2 randomly) +salt and pepper(selected from 0 to 20% randomly) in all bands</td></tr><tr><td>Case 9</td><td>Gaussian (zero-mean, and noise variance is different for different bands and is selected from 0.1~0.2 randomly) +salt and pepper(selected from 0 to 20% randomly) in all bands + dead-lines in 61~70 bands (with the number of stripes randomly selected from 1 to 10 and the width of the dead-lines randomly generated from 1 to 3)</td></tr></tbody></table>

Extracting all items containing $V$ from (37), similar to the process of solving $V$ in the HSI denoising model, the detail of the $l_0-l_1$ HTV method for HSI CS is summarized in Algorithm 2.

## VI. EXPERIMENTAL RESULTS AND DISCUSSION

In this section, extensive experiments on HSI denoising and HSI compressed sensing are conducted to evaluate the performances of the proposed method. In the denoising experiments, we select four types of methods: TV based method, i.e., 2DTV, ASSTV [26], SSTV [27], HSSTV [28], E-3DTV [39], and LR-based method, i.e., LRMR [29]. In the CS experiments, $l_0-l_1$HTV is compared with including 2DTV, ASSTV [26], SSTV [27], and LR-based method, i.e., LR and joint-sparse recovery (SparCS) [55]. Three TV norms are replaced in (35) and also solved by ADMM.

In simulated experiments, we employ the mean peak signal-to-noise ratio (MPSNR), the mean structural similarity index (MSSIM), the Errur Relative Globale Adimensionnelle de Synthese (ERGAS), and the mean spectral angle distance (MSAD) to evaluate HSI restoration qualities as follows:

$$ \text{MPSNR} = \frac{1}{z} \sum_{i=1}^{z} \text{PSNR}(\mathcal{X}(:, :; i), \bar{\mathcal{X}}(:, :; i)) \quad (46) $$

$$ \text{MSSIM} = \frac{1}{z} \sum_{i=1}^{z} \text{SSIM}(\mathcal{X}(:, :; i), \bar{\mathcal{X}}(:, :; i)) \quad (47) $$

$$ \text{ERGAS} = \sqrt{\frac{1}{z} \sum_{i=1}^{z} \frac{\text{mse}(\mathcal{X}(:, :; i), \bar{\mathcal{X}}(:, :; i))}{\text{Mean}_2(\mathcal{X}(:, :; i))}} \quad (48) $$

$$ \text{MSAD} = \frac{1}{hv} \sum_{i=1}^{hv} \frac{180}{\pi} \times \arccos \frac{(\mathbf{X}_{(3)}(i, :))^T \cdot (\bar{\mathbf{X}}_{(3)}(i, :))}{\|\mathbf{X}_{(3)}(i, :))\| \cdot \|\bar{\mathbf{X}}_{(3)}(i, :))\|} \quad (49) $$

where $\mathcal{X}(:, :; i)$ and $\bar{\mathcal{X}}(:, :; i)$ denote the $i$th band of the reference and recovered images, and $\mathbf{X}_{(3)}(i, :)$ and $\bar{\mathbf{X}}_{(3)}(i, :)$ are the $i$th spectral signatures of two images.

## A. HSI Mixed Noise Removal Experiments

Two HSI data sets are adopted including the Washington dc Mall (WDC) data set and the Indian synthetic (Indian-s) data. The WDC data set was captured by the HS digital imagery collection experiment. A subregion image in the size of 200 × 200 × 160 holds complex structure and texture information. The Indian synthetic (Indian-s) data has been considered in [29], [33], [38], which is generated by USGS splib06 using the ground truth of the Indian Pines data set. The Indian-s data contain 145 × 145 pixels and 175 bands, with evident local smoothness information.

The noise of HSI data generally manifests as a mixture of Gaussian noise, salt and pepper noise, deadline noise, and strips. The details of noise cases are given in Table I. Three different kinds of the Gaussian noise are taken into account. In cases 1~3, the Gaussian noise is randomly added to all bands of HSIs. The noisy intensity of different bands is different, and the SNR values of the Gaussian noise vary from 10~20 dB randomly. For cases 4~6, the variance $\sigma_z^2$ of the Gaussian noise is defined as follows [46]:

$$ \sigma_z^2 = \delta^2 \frac{\exp\{-(z - Z/2)^2/(2\eta^2)\}}{\sum_{z=1}^{Z} \exp\{-(z - Z/2)^2/(2\eta^2)\}} \quad (50) $$

where $\delta$ and $\eta$ are two adjustable parameters, and the mean noise SNR is set as 12.12 dB.

1) WDC Data Set: In terms of visual quality, Fig. 2 shows band 70 of WDC data denoised by eight different methods under Case 5. ASSTV fails to remove all noise completely. Although noise is largely suppressed by 2DTV, its result yields some blurry areas and loses much detailed information. HSSTV and E-3DTV perform comparatively better than SSTV and LRMR in mixed noise removal, but their results exist more or less undesired artifacts in the green rectangle area. As observed, the proposed $l_0-l_1$HTV produces superior denoised results without any oversmoothing or artifacts and, meanwhile, preserves clearer structure and texture information, such as four vertical boundaries in the middle of the recovered image.
<figure><img src="image_3.png" /><figcaption>Fig. 2. Denoised results by all the compared methods for WDC data set. (a) Original band 70. (b) Simulated noise band of case 5. (c) 2DTV. (d) ASSTV. (e) SSTV. (f) LRMR. (g) HSSTV1. (h) HSSTV2. (i) E-3DTV. (j) $l_0-l_1$HTV.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 3. (a) PSNR and (b) SSIM values of each band in Case 5 for the WDC data set.</figcaption></figure>

To quantitatively compare the performances of these methods, we list the assessment results in Table II, with the best ones marked in bold and the second ones underlined. Consistent with visual comparison, the proposed method achieves the highest evaluation scores among seven denoised approaches. SSTV, HSSTV, and E-3DTV outperform LRMR under cases 1~6 since these state-of-the-art TV-based methods ponder more local smoothness information. E-3DTV and LRMR have similar denoised results under cases 7~9. Nevertheless, LRMR explores the LR property of HSI and enables to handle complex Gaussian noise, leading to higher MSSIM values than HSSTV and E-3DTV. As ASSTV utilizes local information along with three directions of HSI, ASSTV gains higher scores than 2DTV. Fig. 3 presents the PSNR and SSIM values of each band in Case 5 for the WDC data set. Due to the same curves of $l_0-l_1$HTV, LRMR, HSSTV, and E-3DTV, the curves in the black dashed frames are enlarged. It is obvious that the proposed method shown in red provides the best performances.

2) *Indian-s Data Set:* Fig. 4 shows the denoised performances of compared methods under Case 9. As shown in Fig. 4(a), Indian-s data has much sharper edges than the WDC data set. Similarly, ASSTV still remains the amount of noise, and 2DTV generates serious oversmoothing. ASSTV

<table><caption>TABLE II
QUANTITATIVE RESULTS OF ALL THE METHODS UNDER DIFFERENT MIXED NOISE CASES FOR WDC DATA SET</caption><thead><tr><th rowspan="2">Noise Case</th><th rowspan="2">Index</th><th rowspan="2">Noisy</th><th colspan="6">WDC Dataset</th></tr><tr><th>2DTV</th><th>ASSTV</th><th>SSTV</th><th>LRMR</th><th>HSSTV1</th><th>HSSTV2</th><th>E-3DTV</th><th>$l_0-l_1$HTV</th></tr></thead><tbody><tr><td rowspan="3">Case 1</td><td>MPSNR</td><td>27.22</td><td>29.96</td><td>32.95</td><td>37.84</td><td>36.85</td><td><strong>38.30</strong></td><td>38.09</td><td>38.18</td><td><strong>38.89</strong></td></tr><tr><td>MSSIM</td><td>0.7973</td><td>0.8779</td><td>0.9476</td><td>0.9776</td><td>0.9797</td><td>0.9809</td><td>0.9794</td><td>0.9819</td><td><strong>0.9834</strong></td></tr><tr><td>EGRAS</td><td>11.9726</td><td>7.6773</td><td>5.5183</td><td>3.3347</td><td>3.2881</td><td>2.9490</td><td>3.0544</td><td>3.0633</td><td><strong>2.7484</strong></td></tr><tr><td rowspan="3">Case 2</td><td>MSAD</td><td>10.6152</td><td>7.5971</td><td>4.7666</td><td>3.1015</td><td>3.6719</td><td><strong>2.8716</strong></td><td>2.9465</td><td>3.0112</td><td><strong>2.6935</strong></td></tr><tr><td>MPSNR</td><td>14.11</td><td>26.14</td><td>22.76</td><td>37.15</td><td>36.00</td><td>37.17</td><td>36.71</td><td><strong>38.13</strong></td><td><strong>38.27</strong></td></tr><tr><td>MSSIM</td><td>0.3140</td><td>0.7781</td><td>0.7005</td><td>0.9767</td><td>0.9752</td><td>0.9776</td><td>0.9763</td><td><strong>0.9804</strong></td><td><strong>0.9813</strong></td></tr><tr><td rowspan="3">Case 3</td><td>EGRAS</td><td>46.9504</td><td>11.1352</td><td>17.8865</td><td>3.5510</td><td>3.6063</td><td>3.2904</td><td>3.4705</td><td><strong>3.1541</strong></td><td><strong>2.9667</strong></td></tr><tr><td>MSAD</td><td>36.9857</td><td>11.6210</td><td>14.4817</td><td>3.3381</td><td>4.0624</td><td>3.2414</td><td>3.4029</td><td><strong>3.0299</strong></td><td><strong>2.9030</strong></td></tr><tr><td>MPSNR</td><td>11.13</td><td>25.04</td><td>19.85</td><td>36.04</td><td>35.15</td><td>33.16</td><td>33.53</td><td><strong>36.06</strong></td><td><strong>36.82</strong></td></tr><tr><td rowspan="3">Case 4</td><td>MSSIM</td><td>0.1756</td><td>0.6995</td><td>0.5508</td><td>0.9689</td><td>0.9679</td><td>0.9575</td><td>0.9564</td><td>0.9702</td><td><strong>0.9750</strong></td></tr><tr><td>EGRAS</td><td>46.9504</td><td>12.7272</td><td>25.3070</td><td>3.7629</td><td>4.1288</td><td>4.8508</td><td>4.7852</td><td><strong>3.9624</strong></td><td><strong>3.4041</strong></td></tr><tr><td>MSAD</td><td>43.6584</td><td>13.3058</td><td>19.1071</td><td>3.7519</td><td>4.6473</td><td>4.7344</td><td>4.7852</td><td><strong>3.8913</strong></td><td><strong>3.4023</strong></td></tr><tr><td rowspan="3">Case 5</td><td>MPSNR</td><td>24.34</td><td>27.40</td><td>31.26</td><td>36.40</td><td>35.13</td><td>36.35</td><td><strong>36.86</strong></td><td>36.19</td><td><strong>37.00</strong></td></tr><tr><td>MSSIM</td><td>0.7008</td><td>0.8029</td><td>0.9099</td><td>0.9674</td><td>0.9689</td><td>0.9713</td><td><strong>0.9747</strong></td><td>0.9698</td><td><strong>0.9734</strong></td></tr><tr><td>EGRAS</td><td>15.2816</td><td>10.3531</td><td>6.9993</td><td>3.8969</td><td>3.6063</td><td>3.6474</td><td><strong>3.6205</strong></td><td>3.8483</td><td><strong>3.5670</strong></td></tr><tr><td rowspan="3">Case 6</td><td>MSAD</td><td>14.1859</td><td>10.0657</td><td>6.1004</td><td>3.6533</td><td>4.0624</td><td>3.5669</td><td><strong>3.4223</strong></td><td>3.7338</td><td><strong>3.3771</strong></td></tr><tr><td>MPSNR</td><td>11.16</td><td>24.56</td><td>19.73</td><td><strong>34.52</strong></td><td>33.28</td><td>32.50</td><td>32.08</td><td><strong>33.67</strong></td><td><strong>35.03</strong></td></tr><tr><td>MSSIM</td><td>0.1662</td><td>0.6796</td><td>0.5383</td><td>0.9541</td><td><strong>0.9539</strong></td><td>0.9428</td><td>0.9376</td><td>0.9481</td><td><strong>0.9605</strong></td></tr><tr><td rowspan="3">Case 7</td><td>MPSNR</td><td>16.48</td><td>20.43</td><td>25.67</td><td>30.10</td><td>31.17</td><td><strong>31.27</strong></td><td>29.43</td><td>31.14</td><td><strong>31.38</strong></td></tr><tr><td>MSSIM</td><td>0.3543</td><td>0.5043</td><td>0.8107</td><td>0.8853</td><td><strong>0.9178</strong></td><td>0.9170</td><td>0.8741</td><td>0.9049</td><td><strong>0.9198</strong></td></tr><tr><td>EGRAS</td><td>34.9656</td><td>23.1351</td><td>12.3470</td><td>7.3824</td><td>6.2509</td><td><strong>6.2464</strong></td><td>7.7847</td><td>6.6029</td><td><strong>6.2090</strong></td></tr><tr><td rowspan="3">Case 8</td><td>MSAD</td><td>32.4990</td><td>22.9750</td><td>10.5510</td><td>7.9050</td><td>7.2501</td><td><strong>6.6695</strong></td><td>8.3656</td><td>6.9146</td><td><strong>6.6289</strong></td></tr><tr><td>MPSNR</td><td>12.58</td><td>15.39</td><td>22.20</td><td>29.41</td><td>29.68</td><td>29.00</td><td>28.92</td><td><strong>29.71</strong></td><td><strong>30.07</strong></td></tr><tr><td>MSSIM</td><td>0.2220</td><td>0.3104</td><td>0.6337</td><td>0.8708</td><td><strong>0.9054</strong></td><td>0.8825</td><td>0.8785</td><td>0.8726</td><td><strong>0.8951</strong></td></tr><tr><td rowspan="3">Case 9</td><td>MPSNR</td><td>41.8212</td><td>34.0091</td><td>15.9943</td><td>8.4385</td><td><strong>8.2500</strong></td><td>8.2950</td><td>8.4206</td><td>8.3783</td><td><strong>7.5561</strong></td></tr><tr><td>MSSIM</td><td>0.2199</td><td>0.3763</td><td>0.6315</td><td>0.8160</td><td><strong>0.9033</strong></td><td>0.8600</td><td>0.8554</td><td>0.8632</td><td><strong>0.8949</strong></td></tr><tr><td>EGRAS</td><td>58.3482</td><td>42.5358</td><td>18.8190</td><td>7.9522</td><td><strong>7.5299</strong></td><td>8.9066</td><td>9.0515</td><td>9.5305</td><td><strong>7.2896</strong></td></tr><tr><td></td><td>MSAD</td><td>42.1183</td><td>33.6726</td><td>16.0386</td><td>8.4101</td><td><strong>8.3486</strong></td><td>8.8264</td><td>8.9889</td><td>8.6956</td><td><strong>7.5511</strong></td></tr></tbody></table>

and SSTV fail to remove strips completely. Unlike SSTV, LRMR restores a better image without lots of artifacts since LRMR exploits LR prior information. HSSTV compensates for artifacts of a pure SSTV algorithm to a certain extent. However, HSSTV just exploits the local differences of HSI, which causes wrong boundary judgment. Although E-3DTV outperforms the abovementioned algorithms, there exist fewer artifacts in the smoothing areas since E-3DTV only utilizes the local spatial information and the sparsity insight of 3DTV. The white area on the left reveals that $l_0-l_1$HTV obtains the best smoothness and the sharpest edges. $l_0-l_1$HTV is superior to the other approaches.

Table II gives the quantitative results of eight denoised algorithms and further proves the best performances of our
<figure><img src="image_5.png" /><figcaption>Fig. 4. Denoised results by all the compared methods for Indian-s data set. (a) Original band 65. (b) Simulated noise band of case 9. (c) 2DTV. (d) ASSTV. (e) SSTV. (f) LRMR. (g) HSSTV1. (h) HSSTV2. (i) E-3DTV. (j) $l_0-l_1$ HTV.</figcaption></figure>

<caption>TABLE III</caption>

QUANTITATIVE RESULTS OF ALL THE METHODS UNDER DIFFERENT MIXED NOISE CASES FOR INDIAN-S DATA SET

<table><thead><tr><th rowspan="2">Noise Case</th><th colspan="10">Indian-s Dataset</th></tr><tr><th>Index</th><th>Noisy</th><th>2DTV</th><th>ASSTV</th><th>SSTV</th><th>LRMR</th><th>HSSTV1</th><th>HSSTV2</th><th>E-3DTV</th><th>$l_0-l_1$HTV</th></tr></thead><tbody><tr><td rowspan="3">Case 1</td><td>MPSNR</td><td>26.34</td><td>31.80</td><td>34.04</td><td>42.11</td><td>41.50</td><td>42.20</td><td>40.84</td><td>42.09</td><td>48.36</td></tr><tr><td>MSSIM</td><td>0.5125</td><td>0.9092</td><td>0.7987</td><td>0.9706</td><td>0.9554</td><td>0.9853</td><td>0.9723</td><td>0.9864</td><td>0.9976</td></tr><tr><td>ERGAS</td><td>11.4981</td><td>5.3534</td><td>5.1290</td><td>1.6352</td><td>1.7970</td><td>1.6658</td><td>1.9309</td><td>2.3059</td><td>0.8380</td></tr><tr><td rowspan="3">Case 2</td><td>MSAD</td><td>1.8373</td><td>5.6082</td><td>4.0230</td><td>1.7533</td><td>1.8373</td><td>1.7766</td><td>2.0184</td><td>1.6364</td><td>0.8329</td></tr><tr><td>MPSNR</td><td>14.10</td><td>30.55</td><td>25.89</td><td>40.76</td><td>40.18</td><td>40.92</td><td>39.17</td><td>41.73</td><td>46.86</td></tr><tr><td>MSSIM</td><td>0.1416</td><td>0.8869</td><td>0.6865</td><td>0.9457</td><td>0.9264</td><td>0.9775</td><td>0.9506</td><td>0.9882</td><td>0.9967</td></tr><tr><td rowspan="3">Case 3</td><td>ERGAS</td><td>45.2641</td><td>6.2247</td><td>12.7398</td><td>1.8873</td><td>2.0617</td><td>1.8495</td><td>2.3058</td><td>1.9456</td><td>0.9965</td></tr><tr><td>MSAD</td><td>34.4532</td><td>6.5157</td><td>8.9093</td><td>2.0229</td><td>2.1259</td><td>1.8881</td><td>2.2816</td><td>1.7073</td><td>0.9866</td></tr><tr><td>MPSNR</td><td>11.38</td><td>29.36</td><td>22.10</td><td>39.33</td><td>38.89</td><td>38.36</td><td>36.12</td><td>40.63</td><td>45.50</td></tr><tr><td rowspan="3">Case 4</td><td>MSSIM</td><td>0.0734</td><td>0.8613</td><td>0.5166</td><td>0.9457</td><td>0.9264</td><td>0.9676</td><td>0.9142</td><td>0.9843</td><td>0.9958</td></tr><tr><td>ERGAS</td><td>63.1225</td><td>7.1611</td><td>20.3016</td><td>2.2232</td><td>2.4029</td><td>2.4883</td><td>3.2396</td><td>2.2755</td><td>1.0183</td></tr><tr><td>MSAD</td><td>39.9088</td><td>7.4905</td><td>12.3116</td><td>2.3942</td><td>2.4565</td><td>2.4307</td><td>3.1494</td><td>1.9189</td><td>1.0323</td></tr><tr><td rowspan="3">Case 5</td><td>MPSNR</td><td>23.51</td><td>30.79</td><td>30.83</td><td>39.17</td><td>38.84</td><td>40.78</td><td>39.97</td><td>40.07</td><td>45.89</td></tr><tr><td>MSSIM</td><td>0.4147</td><td>0.8796</td><td>0.6741</td><td>0.9682</td><td>0.9081</td><td>0.9726</td><td>0.9317</td><td>0.9742</td><td>0.9974</td></tr><tr><td>ERGAS</td><td>16.0971</td><td>6.0026</td><td>8.1258</td><td>2.3333</td><td>2.5499</td><td>1.9414</td><td>2.5508</td><td>2.6449</td><td>1.1403</td></tr><tr><td rowspan="3">Case 6</td><td>MSAD</td><td>14.6163</td><td>6.1425</td><td>6.0690</td><td>2.3942</td><td>2.4827</td><td>1.9183</td><td>2.3623</td><td>2.0342</td><td>1.0600</td></tr><tr><td>MPSNR</td><td>13.62</td><td>29.78</td><td>25.25</td><td>38.14</td><td>37.91</td><td>38.95</td><td>36.76</td><td>39.59</td><td>44.69</td></tr><tr><td>MSSIM</td><td>0.1264</td><td>0.8606</td><td>0.6348</td><td>0.9075</td><td>0.8965</td><td>0.9632</td><td>0.8987</td><td>0.9739</td><td>0.9940</td></tr><tr><td rowspan="3">Case 7</td><td>ERGAS</td><td>16.0887</td><td>6.7552</td><td>13.1283</td><td>2.6097</td><td>2.8033</td><td>2.3932</td><td>3.2468</td><td>2.9820</td><td>1.2872</td></tr><tr><td>MSAD</td><td>36.5602</td><td>6.9482</td><td>9.3958</td><td>2.6529</td><td>2.7355</td><td>2.3290</td><td>3.0067</td><td>2.1509</td><td>1.2322</td></tr><tr><td>MPSNR</td><td>11.13</td><td>28.40</td><td>21.67</td><td>36.81</td><td>36.56</td><td>36.54</td><td>34.09</td><td>37.93</td><td>42.85</td></tr><tr><td rowspan="3">Case 8</td><td>MSSIM</td><td>0.0681</td><td>0.8275</td><td>0.4806</td><td>0.8953</td><td>0.8715</td><td>0.9458</td><td>0.8492</td><td>0.9695</td><td>0.9902</td></tr><tr><td>ERGAS</td><td>46.3208</td><td>9.7260</td><td>20.6494</td><td>3.0334</td><td>3.4358</td><td>3.1618</td><td>4.3455</td><td>3.0903</td><td>1.6262</td></tr><tr><td>MSAD</td><td>35.5602</td><td>8.1290</td><td>12.7695</td><td>3.0798</td><td>3.1478</td><td>2.9339</td><td>3.9652</td><td>2.5740</td><td>1.5174</td></tr><tr><td rowspan="3">Case 9</td><td>MPSNR</td><td>16.07</td><td>28.40</td><td>27.64</td><td>32.04</td><td>34.66</td><td>36.09</td><td>32.53</td><td>39.64</td><td>40.25</td></tr><tr><td>MSSIM</td><td>0.1499</td><td>0.7921</td><td>0.6297</td><td>0.7656</td><td>0.8819</td><td>0.9468</td><td>0.8385</td><td>0.9804</td><td>0.9880</td></tr><tr><td>ERGAS</td><td>36.0327</td><td>8.3164</td><td>9.9599</td><td>5.9355</td><td>3.7659</td><td>3.4308</td><td>5.1724</td><td>2.5135</td><td>2.1761</td></tr><tr><td rowspan="3">Case 10</td><td>MSAD</td><td>32.9956</td><td>9.2680</td><td>9.2830</td><td>6.6732</td><td>4.2082</td><td>3.8704</td><td>5.8814</td><td>2.3930</td><td>2.4115</td></tr><tr><td>MPSNR</td><td>12.93</td><td>27.46</td><td>23.17</td><td>31.32</td><td>33.86</td><td>33.19</td><td>30.08</td><td>36.92</td><td>39.37</td></tr><tr><td>MSSIM</td><td>0.0956</td><td>0.7982</td><td>0.5024</td><td>0.7403</td><td>0.8634</td><td>0.9231</td><td>0.7656</td><td>0.9724</td><td>0.9838</td></tr><tr><td rowspan="3">Case 11</td><td>ERGAS</td><td>54.4352</td><td>9.0489</td><td>18.7659</td><td>6.3556</td><td>4.1185</td><td>5.3218</td><td>7.2719</td><td>3.2275</td><td>2.4493</td></tr><tr><td>MSAD</td><td>36.5381</td><td>10.0114</td><td>12.2540</td><td>7.0768</td><td>4.5191</td><td>5.1050</td><td>7.2911</td><td>3.1632</td><td>2.6785</td></tr><tr><td>MPSNR</td><td>12.87</td><td>27.19</td><td>23.10</td><td>31.22</td><td>33.78</td><td>32.96</td><td>29.87</td><td>36.80</td><td>39.28</td></tr><tr><td rowspan="3">Case 12</td><td>MSSIM</td><td>0.0951</td><td>0.7942</td><td>0.5014</td><td>0.7403</td><td>0.8615</td><td>0.9208</td><td>0.7603</td><td>0.9720</td><td>0.9852</td></tr><tr><td>ERGAS</td><td>54.5384</td><td>9.2997</td><td>18.8065</td><td>6.3824</td><td>4.1608</td><td>5.3789</td><td>7.3462</td><td>3.2776</td><td>2.4148</td></tr><tr><td>MSAD</td><td>36.5962</td><td>10.1866</td><td>12.2709</td><td>7.1077</td><td>4.5275</td><td>5.1436</td><td>7.3547</td><td>3.1931</td><td>2.6458</td></tr></tbody></table>

proposed $l_0-l_1$ HTV. For this kind of data set with sharper edges, E-3DTV and HSSTV bring better results than 2DTV, ASSTV, and LRMR in most of the cases. E-3DTV obtains the second-highest metric values among seven comparing methods. HSSTV deals with complex Gaussian noise cases, such as cases 8 and 9, a little worse than LRMR. Due to the good use of local HSI information, HSSTV obtains better MSSIM scores than LRMR. Meanwhile, as we can see in Fig. 5, the PSNR and SSIM values of each band obtained by the proposed $l_0-l_1$ HTV are much higher than the other state-of-the-art algorithms.

3) Real-World Data Set: The Indian Pines data set was acquired by the NASA AVIRIS instrument in Northwestern Indiana, containing 145 × 145 and 224 spectral reflectance bands. Fig. 6 shows that the original band 220 of this scene is corrupted by mixed noise, and target objects can barely be recognized, which directly hinders further HSI analysis.

<figure><img src="image_6.png" /><figcaption>Fig. 5. (a) PSNR and (b) SSIM values of each band in case 9 for Indian-s data set.</figcaption></figure>

The above eight methods are implemented for HSI denoising. 2DTV and ASSTV lose their utility when they encounter heavy noise cases. SSTV and LRMR remove lots of noise, but they produce artifacts. The result denoised by E-3DTV appears blurred image edges and fewer artifacts. In contrast, HSSTV and $l_0-l_1$ HTV restore the clearest image without artifacts.

In real data experiments, the above image assessments are useless since the reference images are unknown. Meanwhile, preserving sharp image edges benefits further HSI tasks, such as classification. Therefore, the supervised classification experiment is designed, and the well-known support vector machines (SVM) [56]–[58] is applied to quantitatively compare the performances of different methods. The first and second columns of Table IV show 16 different classes and the number of training (random 10%) and test samples. The accuracy for each class, kappa coefficient, overall accuracy (OA), and average accuracy (AA) are depicted in Table IV, and eight denoising techniques can be quantitatively evaluated. Fig. 7 presents the visual comparison of classification results before and after denoising. It is conspicuously seen that classification accuracies are significantly fortified after the denoising process. Due to the oversmoothing phenomenon, the results by 2DTV and ASSTV appear large misclassified areas. Therefore, their kappa, OA, and AA scores have risen less. E-3DTV presents higher classification accuracies than LRMR since E-3DTV explores the sparsity of subspace on the gradient maps and is one of the advanced 3DTV variant algorithms. SSTV provides higher metric values, but artifacts cause incorrect classification results for each class. As kappa, OA, and AA show, two HSSTV methods have similar classification results, which are a little lower than the proposed method. $l_0-l_1$ HTV produces the highest evaluation scores, with kappa, OA, and AA of 0.9594, 0.9644, and 0.9666, respectively.
<figure><img src="image_7.png" /><figcaption>Fig. 6. Denoised results by all the compared methods on the Indian Pine data. (a) Original band 220. (b) 2DTV. (c) ASSTV. (d) SSTV. (e) LRMR. (f) HSSTV1. (g) HSSTV2. (h) E-3DTV. (i) $l_0-l_1$HTV.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 7. Classification results by all the compared methods in the Indian Pine data. (a) Original. (b) 2DTV. (c) ASSTV. (d) SSTV. (e) LRMR. (f) HSSTV1. (g) HSSTV2. (h) E-3DTV. (i) $l_0-l_1$HTV.</figcaption></figure>

<table><caption>TABLE IV
CLASSIFICATION RESULTS ON INDIAN PINES DATA BY DIFFERENT METHODS</caption><thead><tr><th>Class</th><th>Train/Test</th><th>Noisy</th><th>2DTV</th><th>ASSTV</th><th>SSTV</th><th>LRMR</th><th>HSSTV1</th><th>HSSTV2</th><th>E-3DTV</th><th>$l_0 - l_1$ HTV</th></tr></thead><tbody><tr><td>Alfalfa</td><td>5/41</td><td>0.6341</td><td>0.6585</td><td>0.7073</td><td>0.7804</td><td>0.7317</td><td><u>0.8780</u></td><td>0.8537</td><td>0.6585</td><td><strong>0.9512</strong></td></tr><tr><td>Corn-notill</td><td>143/1285</td><td>0.8202</td><td>0.9253</td><td>0.8677</td><td>0.9167</td><td>0.9105</td><td><u>0.9183</u></td><td><u>0.9198</u></td><td>0.8996</td><td><strong>0.9447</strong></td></tr><tr><td>Corn-mintill</td><td>83/747</td><td>0.7161</td><td>0.9076</td><td>0.7778</td><td>0.9518</td><td>0.8955</td><td><u>0.9572</u></td><td><strong><u>0.9625</u></strong></td><td>0.8983</td><td>0.9394</td></tr><tr><td>Corn</td><td>24/213</td><td>0.8028</td><td><strong>0.9848</strong></td><td>0.8967</td><td>0.9154</td><td>0.9248</td><td><u>0.9765</u></td><td>0.9390</td><td>0.9577</td><td>0.9437</td></tr><tr><td>Grass-pasture</td><td>49/434</td><td>0.9377</td><td><strong>0.9931</strong></td><td>0.9470</td><td>0.9723</td><td>0.9354</td><td><u>0.9700</u></td><td>0.9677</td><td>0.9493</td><td><u>0.9816</u></td></tr><tr><td>Grass-trees</td><td>73/657</td><td>0.9497</td><td>0.9833</td><td>0.9696</td><td>0.9863</td><td>0.9619</td><td><u>0.9897</u></td><td>0.9848</td><td>0.9833</td><td><strong>0.9985</strong></td></tr><tr><td>Grass-pasture-mowed</td><td>3/25</td><td>0.6400</td><td>0.7200</td><td>0.7600</td><td>0.8400</td><td>0.6400</td><td><u>0.8800</u></td><td>0.800</td><td>0.5600</td><td><strong>0.9600</strong></td></tr><tr><td>Hay-windrowed</td><td>48/430</td><td>0.9906</td><td>0.9930</td><td><strong>1.0000</strong></td><td>0.9953</td><td>0.9976</td><td><u>0.9977</u></td><td><u>0.9977</u></td><td><u>0.9977</u></td><td><strong>1.0000</strong></td></tr><tr><td>Oats</td><td>2/18</td><td>0.7222</td><td><u>0.9444</u></td><td><strong>1.0000</strong></td><td><u>0.9444</u></td><td>0.8333</td><td>0.9359</td><td>0.9444</td><td><strong>1.0000</strong></td><td><strong>1.0000</strong></td></tr><tr><td>Soybean-notill</td><td>98/874</td><td>0.7059</td><td>0.8272</td><td>0.7780</td><td><strong>0.9405</strong></td><td>0.8821</td><td><u>0.9359</u></td><td>0.9176</td><td>0.8833</td><td>0.9245</td></tr><tr><td>Soybean-mintill</td><td>246/2209</td><td>0.8415</td><td>0.9140</td><td>0.8873</td><td><u>0.9601</u></td><td>0.9388</td><td>0.9520</td><td>0.9579</td><td>0.9511</td><td><strong>0.9765</strong></td></tr><tr><td>Soybean-clean</td><td>60/533</td><td>0.8480</td><td>0.9587</td><td>0.9043</td><td><u>0.9699</u></td><td>0.9324</td><td><strong>0.9737</strong></td><td>0.9625</td><td>0.9043</td><td>0.9437</td></tr><tr><td>Wheat</td><td>21/184</td><td><u>0.9945</u></td><td><strong>0.9946</strong></td><td><strong>0.9946</strong></td><td><u>0.9945</u></td><td>0.9891</td><td>0.9736</td><td>0.9891</td><td><strong>0.9946</strong></td><td>0.9783</td></tr><tr><td>Woods</td><td>127/1138</td><td>0.9516</td><td><u>0.9824</u></td><td>0.9754</td><td>0.9674</td><td>0.9455</td><td>0.9736</td><td><u>0.9824</u></td><td>0.9807</td><td><strong>0.9851</strong></td></tr><tr><td>Buildings-Grass-Trees-Drives</td><td>39/347</td><td>0.6311</td><td>0.8559</td><td>0.6513</td><td>0.8587</td><td>0.8040</td><td><strong>0.9510</strong></td><td><u>0.9424</u></td><td>0.9049</td><td>0.9308</td></tr><tr><td>Stone-Steel-Towers</td><td>10/83</td><td>0.9518</td><td><strong>1.0000</strong></td><td>0.9518</td><td>0.9638</td><td>0.9759</td><td><strong>1.0000</strong></td><td><strong>1.0000</strong></td><td><strong>1.0000</strong></td><td><strong>1.0000</strong></td></tr><tr><td>kappa</td><td>-</td><td>0.8197</td><td>0.9197</td><td>0.8677</td><td>0.9451</td><td>0.9131</td><td><u>0.9510</u></td><td>0.9496</td><td>0.9271</td><td><strong>0.9594</strong></td></tr><tr><td>OA</td><td>-</td><td>0.8423</td><td>0.9293</td><td>0.8842</td><td>0.9518</td><td>0.9237</td><td><u>0.9569</u></td><td>0.9557</td><td>0.9360</td><td><strong>0.9644</strong></td></tr><tr><td>AA</td><td>-</td><td>0.8212</td><td>0.9192</td><td>0.8793</td><td>0.9349</td><td>0.8911</td><td><u>0.9588</u></td><td>0.9451</td><td>0.9077</td><td><strong>0.9666</strong></td></tr></tbody></table>

## B. HSI Compressed Sensing Experiments

Three HSI data experiments are employed to validate the effectiveness of $l_0-l_1$HTV on HSI CS with seven different sample ratio, i.e., 0.05, 0.1, 0.2, 0.4, and 0.6. Three data sets include Pavia University, CAVE Toy, and Reno, and their cubes used for HSI CS experiments are of sizes 250 × 121 × 130, 300 × 300 × 31, and 150 × 150 × 100.

1) Pavia University: In terms of visual quality, the representative bands of eight HSI CS results in the sampling ratio $r = 0.4$ are presented in Fig. 8. As shown in the green square of the enlarged areas, reconstructed images by 2DTV, ASSTV, and HSSTV are smoothed excessively, which causes the loss of street detailed information. As observed from the white areas, SSTV and SparCS lead to more or fewer

artifacts. Meanwhile, HSSTV overcomes the artifacts of SSTV, and E-3DTV compensates for the defect of oversmoothing. $l_0-l_1$HTV provides the best result without any artifact, and the comparison of PSNR and SSIM values of all the bands in Fig. 9 also illustrates the superiority of our $l_0-l_1$HTV method.

Table V presents the quantitative indices of eight different sampling ratios on the Pavia university data set. SSTV and SparCS can reconstruct great images in high sampling ratios, but poor results are obtained in low sampling ratios. Nevertheless, the results of 2DTV are just the opposite. As two modified TV-based methods, HSSTV and E-3DTV always obtain better results than SSTV, but they still have quality decays with sample ratios decreasing. These indicate the instabilities of the above three methods. $l_0-l_1$HTV provides much
<figure><img src="image_9.png" /><figcaption>Fig. 8. HSI CS results on the Pavia University data set, and the sample ratio is 0.4. (a) Original band 1. (b) 2DTV. (c) ASSTV. (d) SSTV. (e) SparCS. (f) HSSTV1. (g) HSSTV2. (h) E-3DTV. (i) $l_0-l_1$ HTV.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 9. (a) PSNR and (b) SSIM values of each band for Pavia University under $r = 0.4$.</figcaption></figure>

<table><caption>TABLE V
QUANTITATIVE RESULTS OF ALL THE METHODS UNDER DIFFERENT CS CASES FOR THE PAVIA UNIVERSITY DATA SET</caption><thead><tr><th rowspan="2">Ratio</th><th rowspan="2">Index</th><th colspan="7">Pavia university Dataset</th></tr><tr><th>2DTV</th><th>ASSTV</th><th>SSTV</th><th>SparCS</th><th>HSSTV1</th><th>HSSTV2</th><th>E-3DTV</th><th>$l_0-l_1$HTV</th></tr></thead><tbody><tr><td rowspan="5">0.05</td><td>MPSNR</td><td>21.48</td><td>20.37</td><td>21.17</td><td>19.71</td><td>25.15</td><td>22.20</td><td>26.89</td><td><strong>29.63</strong></td></tr><tr><td>MSSIM</td><td>0.4587</td><td>0.3174</td><td>0.4446</td><td>0.3002</td><td>0.7009</td><td>0.6681</td><td>0.7508</td><td><strong>0.8492</strong></td></tr><tr><td>ERGAS</td><td>20.0014</td><td>22.8979</td><td>20.9138</td><td>25.0856</td><td>13.3918</td><td>18.5314</td><td>10.8490</td><td><strong>7.8779</strong></td></tr><tr><td>MSAD</td><td>23.8498</td><td>27.7240</td><td>25.0375</td><td>29.2799</td><td>14.6353</td><td>19.8596</td><td>11.7016</td><td><strong>8.8695</strong></td></tr><tr><td>MPSNR</td><td>23.43</td><td>22.04</td><td>23.01</td><td>20.83</td><td>28.57</td><td>27.93</td><td>30.39</td><td><strong>36.33</strong></td></tr><tr><td rowspan="5">0.1</td><td>MSSIM</td><td>0.5374</td><td>0.4216</td><td>0.5270</td><td>0.3597</td><td>0.9036</td><td>0.8964</td><td>0.8613</td><td><strong>0.9411</strong></td></tr><tr><td>ERGAS</td><td>15.8906</td><td>18.8693</td><td>16.9183</td><td>23.2884</td><td>9.2553</td><td>9.8463</td><td>7.2841</td><td><strong>3.7945</strong></td></tr><tr><td>MSAD</td><td>18.7013</td><td>22.4843</td><td>19.9202</td><td>20.8280</td><td>10.2367</td><td>10.9056</td><td>7.6556</td><td><strong>4.2857</strong></td></tr><tr><td>MPSNR</td><td>25.69</td><td>25.57</td><td>27.79</td><td>26.19</td><td>36.49</td><td>32.09</td><td>34.24</td><td><strong>41.11</strong></td></tr><tr><td>MSSIM</td><td>0.6531</td><td>0.5768</td><td>0.7175</td><td>0.6363</td><td>0.9666</td><td>0.9560</td><td>0.9390</td><td><strong>0.9717</strong></td></tr><tr><td rowspan="5">0.2</td><td>ERGAS</td><td>12.2199</td><td>12.5174</td><td>9.7927</td><td>12.5922</td><td>4.0403</td><td>6.4183</td><td>4.7032</td><td><strong>2.4319</strong></td></tr><tr><td>MSAD</td><td>14.2670</td><td>14.6661</td><td>11.3162</td><td>13.7806</td><td>4.3824</td><td>7.0270</td><td>5.3836</td><td><strong>2.5706</strong></td></tr><tr><td>MPSNR</td><td>28.80</td><td>29.13</td><td>37.60</td><td>39.23</td><td>43.11</td><td>42.14</td><td>42.82</td><td><strong>46.15</strong></td></tr><tr><td>MSSIM</td><td>0.7841</td><td>0.7160</td><td>0.9389</td><td>0.9173</td><td>0.9803</td><td>0.9801</td><td>0.9826</td><td><strong>0.9844</strong></td></tr><tr><td>ERGAS</td><td>8.5431</td><td>8.3255</td><td>3.3359</td><td>4.5526</td><td>2.1228</td><td>2.3508</td><td>1.9911</td><td><strong>1.6831</strong></td></tr><tr><td rowspan="5">0.4</td><td>MSAD</td><td>9.8675</td><td>9.7278</td><td>3.7420</td><td>3.9402</td><td>2.1740</td><td>2.4210</td><td>2.1164</td><td><strong>1.5689</strong></td></tr><tr><td>MPSNR</td><td>31.15</td><td>34.69</td><td>45.58</td><td>47.04</td><td>47.25</td><td>46.32</td><td>46.81</td><td><strong>47.92</strong></td></tr><tr><td>MSSIM</td><td>0.8533</td><td>0.8875</td><td>0.9836</td><td>0.9875</td><td>0.9858</td><td>0.9852</td><td>0.9909</td><td><strong>0.9881</strong></td></tr><tr><td>ERGAS</td><td>6.5299</td><td>4.4228</td><td>1.7491</td><td>1.6358</td><td><strong>1.4127</strong></td><td>1.5741</td><td>1.3848</td><td><strong>1.4497</strong></td></tr><tr><td>MSAD</td><td>7.5053</td><td>5.1355</td><td>1.6679</td><td>1.4617</td><td>1.3862</td><td>1.5693</td><td>1.3894</td><td><strong>1.3116</strong></td></tr></tbody></table>

higher metrics than the other approaches, with lower sampling ratios.

2) *Toy*: Compared with the above Toy HSIs, more regular clothes texture information exists in the magnified areas of Fig. 10(a). Unfortunately, ASSTV and SSTV fail to recover

the final HSI under a low sampling ratio $r = 0.2$. Although a few texture details can be found in the results by 2DTV and SparCS, it is hard to distinguish the colors of texture information. As the zoomed areas of the yellow frame show, much white background information is lost by two HSSTV algorithms. E-3DTV produces a better background result, but there still exist few artifacts. The $l_0-l_1$HTV method provides the best result with reconstructing most information of the original image. The PSNR and SSIM values of each band in Case 9 for the Indian-s data set are shown in Fig. 9. We can observe that, in almost all bands, the metric values of our proposed method (the red curves) are higher than those of the other approaches.

Table VI lists the four metric values of different methods. Under almost all cases, the metric scores obtained by using $l_0-l_1$HTV are the highest among all the techniques. 2DTV, ASSTV, and SSTV outperform SparCS under most sampling ratios, which demonstrates that TV methods are more suitable for recovering HSIs with more texture information than pure LR methods. The MSSIM, ERGAS, and MSAD values of HSSTV are the highest when high sample ratio $r = 0.6$, but the metrics of HSSTV decrease more than E-3DTV and $l_0-l_1$HTV with decreasing sample ratios.

3) *Reno Data Set*: Fig. 12 depicts the reconstructed results of band 2 of the Reno data set by different methods, with the random sample ratio $r = 0.1$. 2DTV causes spatial oversmoothing. SparCS exploits the LR property of HSI, but SparCS ignores spatial smoothness, which leads to the lack of local structure information. Due to the existence of artifacts, blurry image edges can hardly be distinguished in the result recovered by SSTV. Based on local spatial and spatiospectral information, HSSTV recovers some clear image edges but loses lots of pixels in the white areas. It is evident that E-3DTV and $l_0-l_1$HTV get the best reconstruction, but $l_0-l_1$HTV recovers the complete white areas and the sharpest edges. The PSNR and SSIM of each band restored by different approaches are shown in Fig. 13, and the results demonstrate that $l_0-l_1$HTV marked in red achieved the best results in all the bands.
<figure><img src="image_11.png" /><figcaption>Fig. 10. HSI CS results on the Toy data set and sample ratio is 0.2. (a) Original band 2. (b) 2DTV. (c) ASSTV. (d) SSTV. (e) SparCS. (f) HSSTV1. (g) HSSTV2. (h)E-3DTV. (i) $l_0-l_1$HTV.</figcaption></figure>

<table><caption>TABLE VI
QUANTITATIVE RESULTS OF ALL THE METHODS UNDER DIFFERENT CS CASES FOR THE TOY DATA SET</caption><thead><tr><th rowspan="2">Ratio</th><th rowspan="2">Index</th><th colspan="8">Toy Dataset</th></tr><tr><th>2DTV</th><th>ASSTV</th><th>SSTV</th><th>SparCS</th><th>HSSTV1</th><th>HSSTV2</th><th>E-3DTV</th><th>$l_0-l_1$HTV</th></tr></thead><tbody><tr><td rowspan="5">0.05</td><td>MPSNR</td><td>13.18</td><td>12.39</td><td>12.91</td><td>10.89</td><td>14.62</td><td>13.36</td><td>20.42</td><td><b>23.74</b></td></tr><tr><td>MSSIM</td><td>0.2033</td><td>0.1415</td><td>0.1299</td><td>0.0757</td><td>0.3027</td><td>0.6318</td><td>0.3896</td><td><b>0.6234</b></td></tr><tr><td>ERGAS</td><td>47.0874</td><td>51.5439</td><td>48.4195</td><td>61.4419</td><td>39.7529</td><td>45.9325</td><td>20.5011</td><td><b>13.8817</b></td></tr><tr><td>MSAD</td><td>39.2882</td><td>44.2176</td><td>40.7690</td><td>53.3116</td><td>27.8701</td><td>33.6779</td><td>10.7912</td><td><b>16.6667</b></td></tr><tr><td>MPSNR</td><td>15.83</td><td>13.78</td><td>14.51</td><td>12.82</td><td>18.55</td><td>16.10</td><td>27.66</td><td><b>32.23</b></td></tr><tr><td rowspan="5">0.1</td><td>MSSIM</td><td>0.2746</td><td>0.1716</td><td>0.1572</td><td>0.1091</td><td>0.8110</td><td>0.7610</td><td>0.7519</td><td><b>0.8800</b></td></tr><tr><td>ERGAS</td><td>34.8259</td><td>44.0011</td><td>40.2711</td><td>49.3815</td><td>25.2887</td><td>33.5346</td><td>8.8715</td><td><b>5.2783</b></td></tr><tr><td>MSAD</td><td>27.0173</td><td>36.0451</td><td>32.6395</td><td>41.5952</td><td>17.8904</td><td>24.3612</td><td>6.47281</td><td><b>4.0486</b></td></tr><tr><td>MPSNR</td><td>26.30</td><td>19.65</td><td>18.62</td><td>19.11</td><td>29.60</td><td>25.42</td><td>33.85</td><td><b>36.55</b></td></tr><tr><td>MSSIM</td><td>0.7171</td><td>0.3640</td><td>0.2635</td><td>0.2877</td><td>0.8326</td><td>0.7993</td><td>0.8957</td><td><b>0.9462</b></td></tr><tr><td rowspan="5">0.2</td><td>ERGAS</td><td>10.6142</td><td>22.4768</td><td>25.2059</td><td>24.5570</td><td>7.0920</td><td>11.4444</td><td>4.4208</td><td><b>3.3489</b></td></tr><tr><td>MSAD</td><td>7.5504</td><td>16.1117</td><td>19.5294</td><td>41.5952</td><td>5.4126</td><td>8.6127</td><td>3.2960</td><td><b>2.5473</b></td></tr><tr><td>MPSNR</td><td>37.23</td><td>31.63</td><td>31.15</td><td>30.6924</td><td>39.84</td><td>38.26</td><td>40.52</td><td><b>41.23</b></td></tr><tr><td>MSSIM</td><td>0.9480</td><td>0.7557</td><td>0.7128</td><td>0.7133</td><td>0.9842</td><td>0.9803</td><td>0.9673</td><td><b>0.9811</b></td></tr><tr><td>ERGAS</td><td>3.0369</td><td>5.6096</td><td>6.0629</td><td>89464</td><td>2.2174</td><td>2.6428</td><td>2.1044</td><td><b>1.6689</b></td></tr><tr><td rowspan="5">0.4</td><td>MSAD</td><td>2.2502</td><td>4.2858</td><td>4.6485</td><td>5.4921</td><td>1.7069</td><td>2.0394</td><td>1.5899</td><td><b>1.2115</b></td></tr><tr><td>MPSNR</td><td>43.85</td><td>33.97</td><td>44.60</td><td>35.90</td><td><b>45.27</b></td><td>45.01</td><td>44.18</td><td><b>47.92</b></td></tr><tr><td>MSSIM</td><td>0.9884</td><td>0.8182</td><td>0.9816</td><td>0.8457</td><td><b>0.9894</b></td><td>0.9891</td><td>0.9822</td><td>0.9881</td></tr><tr><td>ERGAS</td><td>1.3753</td><td>4.2852</td><td>1.3168</td><td>5.0260</td><td><b>1.1959</b></td><td>1.1866</td><td>1.4131</td><td>1.4497</td></tr><tr><td>MPSNR</td><td>1.0637</td><td>3.3165</td><td>1.0021</td><td>3.1134</td><td>0.9190</td><td><b>0.9015</b></td><td>1.0589</td><td>1.3116</td></tr></tbody></table>

<figure><img src="image_12.png" /><figcaption>Fig. 11. (a) PSNR and (b) SSIM values of each band for Toy under $r = 0.2$.</figcaption></figure>

Table VII provides the quantitative results of eight different sampling ratios on the Reno data set in terms of four evaluation indices. It can be easily seen that the proposed $l_0-l_1$HTV gives a significantly fortified performance under five sampling ratios compared with other competing methods. Since important structure details of the final HSIs are barely reconstructed, 2DTV and ASSTV get similar assessment scores. By using more local differences, E-3DTV and HSSTV recover the

second-best and third-best results, respectively. Especially, E-3DTV is more beneficial for the CS recovery with lower sampling ratios. However, $l_0-l_1$HTV achieves the best reconstructions under all sampling ratios.

In summary, the structure information of the three data sets is different, while our proposed method always obtains the best-recovered images. $l_0-l_1$HTV can be a useful regularizer to ameliorate the performances for HSI denoising and CS tasks.

### C. Discussion

1) **Convergence of the Proposed Methods:** The convergence is verified on HSI denoising and CS experiment data sets. Fig. 14(a) displays the RalCha of two denoised results for WDC and Indian-s data sets under noisy cases 5 and 9, respectively. RalCha is the relative change of two restored results, i.e., $||\mathcal{X}(i+1) - \mathcal{X}(i)||_F^2/||\mathcal{X}(i)||_F^2$, where $\mathcal{X}(i)$ is the $i$th recovered HSI image. Fig. 14(b) shows the RalCha of two recovered results for PaviaUA and Toy data sets under sampling ratio $r = 0.1$ and $r = 0.2$, respectively. It can be spotted that, as the iterations are conducted, the RalCha values of $l_0-l_1$HTV for HSI denoising and CS reconstruction monotonically decrease, finally approaching zero.

2) **Parameter Analysis:** In the $l_0-l_1$HTV model, four parameters should be carefully identified for the HSI denoising and CS. We analyze the impact of each parameter on the restoration results and discuss how to choose such parameters in our experiments. The denoised results are based on the simulated experiments in WDC and Indian-s data sets. The effect of each parameter on the CS is discussed under sampling ratio $r = 0.1$.

Parameter $\gamma$ plays an important role to control the number of $l_0$ gradients. Due to the relationship with an HSI size, $\gamma$ is replaced by the parameter $\gamma'$ that is set to various percentages of total HSI size in Figs. 15(a) and 16(a). $\gamma'$ is chosen from the set [0.1, 1.0] with 0.01 as step size. The Indian-s and Toy data sets have simpler texture information than the WDC and
<figure><img src="image_13.png" /><figcaption>Fig. 12. HSI CS results on the Reno data set, and the sample ratio is 0.1. (a) Original band 2. (b) 2DTV. (c) ASSTV. (d) SSTV. (e) SparCS. (f) HSSTV1. (g) HSSTV2. (h) E-3DTV. (i) $l_0-l_1$ HTV.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 13. (a) PSNR and (b) SSIM values of each band for the Reno data set under $r = 0.1$.</figcaption></figure>

<caption>TABLE VII</caption>

<table><caption>QUANTITATIVE RESULTS OF ALL THE METHODS UNDER DIFFERENT CS CASES FOR THE RENO DATA SET</caption><thead><tr><th rowspan="2">Ratio</th><th rowspan="2">Index</th><th colspan="9">Reno Dataset</th></tr><tr><th>2DTV</th><th>ASSTV</th><th>SSTV</th><th>SparCS</th><th>HSSTV1</th><th>HSSTV2</th><th>E-3DTV</th><th>$l_0-l_1$ HTV</th></tr></thead><tbody><tr><td rowspan="3">0.05</td><td>MPSNR</td><td>15.76</td><td>15.96</td><td>16.21</td><td>13.56</td><td>18.55</td><td>17.04</td><td>22.80</td><td><b>25.11</b></td></tr><tr><td>MSSIM</td><td>0.2574</td><td>0.2743</td><td>0.2891</td><td>0.1488</td><td><u>0.7336</u></td><td>0.6472</td><td>0.6073</td><td><b>0.7724</b></td></tr><tr><td>ERGAS</td><td>33.2416</td><td>32.5009</td><td>31.5179</td><td>45.2039</td><td>24.0110</td><td>28.4047</td><td><u>14.5842</u></td><td><b>11.1754</b></td></tr><tr><td rowspan="3">0.1</td><td>MSAD</td><td>32.1112</td><td>31.2797</td><td>30.2264</td><td>41.6306</td><td>20.7043</td><td>24.7868</td><td><u>13.3769</u></td><td><b>9.8166</b></td></tr><tr><td>MPSNR</td><td>16.87</td><td>17.62</td><td>17.77</td><td>14.48</td><td>24.67</td><td>23.08</td><td><u>27.98</u></td><td><b>32.05</b></td></tr><tr><td>MSSIM</td><td>0.3240</td><td>0.3867</td><td>0.3866</td><td>0.2189</td><td><u>0.9070</u></td><td>0.8743</td><td>0.8876</td><td><b>0.9230</b></td></tr><tr><td rowspan="3">0.2</td><td>ERGAS</td><td>29.2555</td><td>26.8520</td><td>26.3792</td><td>38.1143</td><td>12.1361</td><td>14.4872</td><td><u>8.0589</u></td><td><b>5.0768</b></td></tr><tr><td>MSAD</td><td>27.7898</td><td>25.1334</td><td>24.8107</td><td>14.4849</td><td>10.5626</td><td>12.6363</td><td><u>6.5857</u></td><td><b>4.6492</b></td></tr><tr><td>MPSNR</td><td>20.10</td><td>22.98</td><td>21.84</td><td>22.18</td><td>35.54</td><td>32.76</td><td><u>35.72</u></td><td><b>38.94</b></td></tr><tr><td rowspan="3">0.4</td><td>MSSIM</td><td>0.5054</td><td>0.6348</td><td>0.5921</td><td>0.6234</td><td><u>0.9782</u></td><td>0.9744</td><td>0.9733</td><td><b>0.9792</b></td></tr><tr><td>ERGAS</td><td>20.1403</td><td>14.3473</td><td>16.5678</td><td>16.5344</td><td>3.5288</td><td>5.0153</td><td><u>3.3380</u></td><td><b>2.3849</b></td></tr><tr><td>MSAD</td><td>18.4749</td><td>12.8383</td><td>15.1398</td><td>14.8583</td><td>3.1360</td><td>4.3568</td><td><u>2.9699</u></td><td><b>2.1355</b></td></tr><tr><td rowspan="3">0.6</td><td>MPSNR</td><td>29.54</td><td>30.77</td><td>36.19</td><td>36.84</td><td>38.11</td><td>36.60</td><td><u>42.82</u></td><td><b>45.81</b></td></tr><tr><td>MSSIM</td><td>0.8679</td><td>0.8669</td><td>0.9537</td><td>0.9306</td><td>0.9887</td><td>0.9858</td><td><u>0.9914</u></td><td><b>0.9961</b></td></tr><tr><td>ERGAS</td><td>6.7873</td><td>5.8024</td><td>3.2823</td><td>5.1925</td><td>2.8389</td><td>3.3429</td><td><u>1.5316</u></td><td><b>1.1066</b></td></tr><tr><td rowspan="3">0.8</td><td>MSAD</td><td>6.1679</td><td>5.3277</td><td>2.9512</td><td>3.6990</td><td>2.4193</td><td>2.8556</td><td><u>1.3672</u></td><td><b>0.9465</b></td></tr><tr><td>MPSNR</td><td>33.83</td><td>33.69</td><td>47.64</td><td>47.92</td><td><u>48.34</u></td><td>47.41</td><td>47.78</td><td><b>48.81</b></td></tr><tr><td>MSSIM</td><td>0.9413</td><td>0.9179</td><td>0.9957</td><td>0.9890</td><td><u>0.9974</u></td><td>0.9950</td><td>0.9964</td><td><b>0.9969</b></td></tr><tr><td rowspan="3">1.0</td><td>ERGAS</td><td>4.1420</td><td>4.1573</td><td>0.9421</td><td>1.7419</td><td><u>0.9277</u></td><td>1.0385</td><td>0.9366</td><td><b>0.8553</b></td></tr><tr><td>MSAD</td><td>3.7852</td><td>3.8306</td><td>0.8052</td><td>0.9994</td><td><u>0.7756</u></td><td>0.8616</td><td>0.8017</td><td><b>0.7165</b></td></tr></tbody></table>

Pavia data sets; thus, $\gamma'$ values for the Indian-s and Toy data sets are lower than the others.

Parameter $\lambda$ is related to the strength of sparse noise norm. $l_0-l_1$ HTV is tested using different values of $\lambda$ varied among [0.01,0.1]. It is obvious that $l_0-l_1$ HTV produces an appealing result when $\lambda$ lies in the range of [0.05, 0.06].

Parameter $\lambda_{td}$ restricts the influence of the SSTV regularization on denoising or CS effectiveness. Figs. 15(c) and 16(b) present the robustness of the proposed method to the changes of $\lambda_{td}$. When $\lambda_{td}$ is given within [0.07, 0.08] for denoising and [8, 10] for CS, respectively, PSNR stabilizes at a high value.

Parameter $\beta$ is used for the ADMM multipliers. First, $\beta$ is set as 0.1 for denoising and 0.01 for CS, respectively, and then

<figure><img src="image_15.png" /><figcaption>Fig. 14. RelCha values with respect to the number of iterations. (a) HSI denoising. (b) HSI CS.</figcaption></figure>

tuned gradually. $\beta$ is fixed as 0.3 and 0.05 for denoising and CS, respectively, which achieves the best performances.

3) Computational Complexity: The computational complexities of the proposed algorithm for two HSI restorations are discussed in this section, and the first three parts for denoising and CS processes are universal as follows.

1) Solving the subproblem (22), (23), or (39), (40) includes two thresholds and two SSTV operators, which takes about $O(2hvz + 2hvz)$ iterations.

2) For (24) or (41), the sorting of the L2-norms of $hv$ sub-vectors spends the most time and requires $O(hv(\log hv))$ time.

3) Most of the computation time for updating the multipliers lies on the SSTV and $l_0$HTV operators, which requires about $O(3hvz)$ iterations.

For HSI denoising, the threshold for solving S-related problem requires $O(hvz)$, and the computation of LSQR algorithm used to update $x$ costs $O(4hv + 5z)$.

For HSI CS, the LSQR algorithm is also utilized to update $x$, and its computation costs $O(4hv + hvz)$.

Therefore, the computational complexities of the HSI denoising process and CS the process are in total $O(8hvz + 5z + hv(\log hv))$ and $O(8hvz + 4hv + hv(\log hv))$, respectively.

The running times of different methods on HSI denoising and CS reconstruction experiments are reported in Table VIII. All the experiments are conducted in MATLAB 2018b on the PC with an Intel i7 CPU at 2.60 GHz and 32 GB of memory. For HSI denoising, ASSTV is the fastest one among all the comparing methods, but most of the mixed
<figure><img src="image_16.png" /><figcaption>Fig. 15. MPSNR as a function of the parameters (i.e., (a) $\gamma'$, (b) $\lambda'$, (c) $\lambda_{TV}$, and (d) $\beta'$) in the proposed method for HSI denoising on the WDC and Indian-s data sets.</figcaption></figure>

<table><caption>TABLE VIII
RUNNING TIME COMPARISON OF DIFFERENT METHODS</caption><thead><tr><th>HSI data sets</th><th>size</th><th>2DTV</th><th>ASSTV</th><th>SSTV</th><th>LRMR</th><th>HSSTV1</th><th>HSSTV2</th><th>E-3DTV</th><th>l<sub>0</sub>-l<sub>1</sub>HTV</th></tr></thead><tbody><tr><td>WDC</td><td>200 × 200 × 160</td><td>314.5s</td><td>109.9s</td><td>564.3s</td><td>131.6s</td><td>1710.7s</td><td>1805.6s</td><td>390.5s</td><td>571.9s</td></tr><tr><td>Indian-s</td><td>145 × 145 × 175</td><td>189.9s</td><td>57.4s</td><td>432.3s</td><td>99.2s</td><td>1054.2s</td><td>1061.8s</td><td>193.10s</td><td>349.6s</td></tr><tr><td>Indian Pines</td><td>145 × 145 × 220</td><td>235.5s</td><td>72.5s</td><td>466.5s</td><td>115.1s</td><td>1124.6s</td><td>1286.3s</td><td>211.8s</td><td>400.6s</td></tr><tr><th>HSI data sets</th><th>size</th><th>2DTV</th><th>ASSTV</th><th>SSTV</th><th>SparCS</th><th>HSSTV1</th><th>HSSTV2</th><th>E-3DTV</th><th>l<sub>0</sub>-l<sub>1</sub>HTV</th></tr><tr><td>Pavia University</td><td>250 × 121 × 103</td><td>832.1s</td><td>1247.7s</td><td>941.3s</td><td>1516.9s</td><td>1126.1s</td><td>1274.9s</td><td>288.1s</td><td>778.3s</td></tr><tr><td>Toy</td><td>300 × 300 × 31</td><td>1025.7s</td><td>1303.7s</td><td>1066.2s</td><td>1355.1s</td><td>1334.9s</td><td>1465.4s</td><td>307.4s</td><td>819.7s</td></tr><tr><td>Reno</td><td>150 × 150 × 100</td><td>708.9s</td><td>1012.1s</td><td>734.7s</td><td>1027.5s</td><td>1028.4s</td><td>1135.9s</td><td>267.7s</td><td>679.4s</td></tr></tbody></table>

<figure><img src="image_17.png" /><figcaption>Fig. 16. MPSNR as a function of the parameters (i.e., (a) $\gamma'$, (b) $\lambda'$, (c) $\lambda_{TV}$, and (d) $\beta'$) in the proposed method for HSI CS on the PaviaUA and Toy data sets.</figcaption></figure>

noise remains. LRMR is the second-best and hardly handles heavy Gaussian cases. 2DTV and E-3DTV spend similar computational time because the update processing of 2DTV is simple, and E-3DTV converges quickly. HSSTV uses two local differences and has a slow convergence, which leads to an expensive computation cost. In contrast, although $l_0-l_1$HTV contains two TV norms, $l_0$hTV and SSTV, the running times of $l_0-l_1$HTV and SSTV are similar since the pure SSTV model converges slowly. For HSI CS, ASSTV, SparCS, and HSSTV are more computationally expensive. Due to its fast convergence speed, E-3DTV is the fastest one among all the CS algorithms. Our proposed method is slightly slower than SSTV. However, $l_0-l_1$HTV obtains the best performances in the HSI restoration experiments.

## VII. CONCLUSION

In this article, we have proposed a novel $l_0-l_1$HTV regularization for HSI denoising and CS restoration. $l_0$ HS TV

restrains the spatial smoothness from all spectral bands, and SSTV exploits local structure information along with spatial and spectral directions. The hybrid TV is designed as a globally and locally integrated regularizer for HSI restorations. Experimental results on HSI mixed noise removal and CS recovery demonstrate the superior performances of our proposed methods.

For further work, two extensions of our regularization technique can be explored: 1) $l_0-l_1$HTV can be applied to other HSI processing tasks, such as classification and 2) in this study, the simple use of $l_0-l_1$HTV has been focused. In fact, this kind of HTV can be injected into LR matrix or tensor models.

## REFERENCES

[1] J. M. Bioccas-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, Jun. 2013.

[2] B. Rasti, P. Scheunders, P. Ghamisi, G. Licciardi, and J. Chanussot, “Noise reduction in hyperspectral imagery: Overview and application,” *Remote Sens.*, vol. 10, no. 3, p. 482, Mar. 2018.

[3] B. Rasti et al., “Feature extraction for hyperspectral imagery: The evolution from shallow to deep: Overview and toolbox,” *IEEE Geosci. Remote Sens. Mag.*, vol. 8, no. 4, pp. 60–88, Dec. 2020.

[4] Y. Gu, J. Chanussot, X. Jia, and J. A. Benediktsson, “Multiple Kernel learning for hyperspectral image classification: A review,” *IEEE Trans. Geosci. Remote Sens.*, vol. 55, no. 11, pp. 6547–6565, Nov. 2017.

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

[6] F. Liu and Q. Wang, “A sparse tensor-based classification method of hyperspectral image,” *Signal Process.*, vol. 168, Mar. 2020, Art. no. 107361.

[7] J. M. Bioccas-Dias et al., “Hyperspectral unmixing overview: Geometrical, statistical, and sparse regression-based approaches,” *IEEE J. Sel. Topics Appl. Earth Observ. Remote Sens.*, vol. 5, no. 2, pp. 354–379, Apr. 2012.

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

[9] H. Zhang, L. Zhang, and H. Shen, “A super-resolution reconstruction algorithm for hyperspectral images,” *Signal Process.*, vol. 92, no. 9, pp. 2082–2096, Sep. 2012.
[10] 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.

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

[12] D. Hong, N. Yokoya, J. Chanussot, and X. X. Zhu, "CoSpace: Common subspace learning from hyperspectral-multispectral correspondences," *IEEE Trans. Geosci. Remote Sens.*, vol. 57, no. 7, pp. 4349-4359, Jul. 2019.

[13] Y. Xu, Z. Wu, J. Chanussot, and Z. Wei, "Joint reconstruction and anomaly detection from compressive hyperspectral images using Mahalanobis distance-regularized tensor RPCA," *IEEE Trans. Geosci. Remote Sens.*, vol. 56, no. 5, pp. 2919-2930, May 2018.

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

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

[16] L. I. Rudin, S. Osher, and E. Fatemi, "Nonlinear total variation based noise removal algorithms," *Phys. D, Nonlinear Phenomena*, vol. 60, nos. 1-4, pp. 259-268, Nov. 1992.

[17] D. M. Strong, P. Blomgren, and T. F. Chan, "Spatially adaptive local-feature-driven total variation minimizing image restoration," *Proc. SPIE*, vol. 3167, pp. 222-233, Oct. 1997.

[18] A. Chopra and H. Lian, "Total variation, adaptive total variation and nonconvex smoothly clipped absolute deviation penalty for denoising blocky images," *Pattern Recognit.*, vol. 43, no. 8, pp. 2609-2619, Aug. 2010.

[19] D. Needell and R. Ward, "Stable image reconstruction using total variation minimization," *SIAM J. Imag. Sci.*, vol. 6, no. 2, pp. 1035-1058, Jan. 2013.

[20] C. Poon, "On the role of total variation in compressed sensing," *SIAM J. Imag. Sci.*, vol. 8, no. 1, pp. 682-720, Jan. 2015.

[21] D. Needell and R. Ward, "Near-optimal compressed sensing guarantees for total variation minimization," *IEEE Trans. Image Process.*, vol. 22, no. 10, pp. 3941-3949, Oct. 2013.

[22] 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, Oct. 2012.

[23] C. Li, T. Sun, K. F. Kelly, and Y. Zhang, "A compressive sensing and unmixing scheme for hyperspectral data processing," *IEEE Trans. Image Process.*, vol. 21, no. 3, pp. 1200-1210, Mar. 2012.

[24] D. T. Eason and M. Andrews, "Total variation regularization via continuation to recover compressed hyperspectral images," *IEEE Trans. Image Process.*, vol. 24, no. 1, pp. 284-293, Jan. 2015.

[25] S. Ono, K. Shirai, and M. Okuda, "Vectorial total variation based on arranged structure tensor for multichannel image restoration," in *Proc. IEEE Int. Conf. Acoust., Speech Signal Process. (ICASSP)*, Mar. 2016, pp. 4528-4532.

[26] Y. Chang, L. Yan, H. Fang, and H. Liu, "Simultaneous destriping and denoising for remote sensing images with unidirectional total variation and sparse representation," *IEEE Geosci. Remote Sens. Lett.*, vol. 11, no. 6, pp. 1051-1055, Jun. 2014.

[27] H. K. Aggarwal and A. Majumdar, "Hyperspectral image denoising using spatio-spectral total variation," *IEEE Geosci. Remote Sens. Lett.*, vol. 13, no. 3, pp. 442-446, Mar. 2016.

[28] S. Takeyama, S. Ono, and I. Kumazawa, "A constrained convex optimization approach to hyperspectral image restoration with hybrid spatio-spectral regularization," 2019, *arXiv:1907.13357*. [Online]. Available: http://arxiv.org/abs/1907.13357

[29] 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, Aug. 2014.

[30] 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, Jan. 2016.

[31] Z. Wu, Q. Wang, Z. Wu, and Y. Shen, "Total variation-regularized weighted nuclear norm minimization for hyperspectral image mixed denoising," *J. Electron. Imag.*, vol. 25, no. 1, Feb. 2016, Art. no. 013037.

[32] M. Golbabaee and P. Vanderghynst, "Joint trace/TV norm minimization: A new efficient approach for spectral compressive imaging," in *Proc. 19th IEEE Int. Conf. Image Process.*, Sep. 2012, pp. 933-936.

[33] 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, Feb. 2017.

[34] Q. Wang, Z. Wu, J. Jin, T. Wang, and Y. Shen, "Low rank constraint and spatial spectral total variation for hyperspectral image mixed denoising," *Signal Process.*, vol. 142, pp. 11-26, Jan. 2018.

[35] M. Wang, Q. Wang, J. Chanussot, and D. Li, "Hyperspectral image mixed noise removal based on multidirectional low-rank modeling and spatial-spectral total variation," *IEEE Trans. Geosci. Remote Sens.*, vol. 59, no. 1, pp. 488-507, Jan. 2021, doi: 10.1109/TGRS. 2020.2993631.

[36] 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.

[37] W. He, H. Zhang, H. Shen, and L. Zhang, "Hyperspectral image denoising using local low-rank matrix recovery and global spatial-spectral total variation," *IEEE J. Sel. Topics Appl. Earth Observ. Remote Sens.*, vol. 11, no. 3, pp. 713-729, Mar. 2018.

[38] 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.

[39] J. Peng, Q. Xie, Q. Zhao, Y. Wang, L. Yee, and D. Meng, "Enhanced 3DTV regularization and its applications on HSI denoising and compressed sensing," *IEEE Trans. Image Process.*, vol. 29, pp. 7889-7903, 2020.

[40] Y. Wang, L. Lin, Q. Zhao, T. Yue, D. Meng, and Y. Leung, "Compressive sensing of hyperspectral images via joint tensor tucker decomposition and weighted total variation regularization," *IEEE Geosci. Remote Sens. Lett.*, vol. 14, no. 12, pp. 2457-2461, Dec. 2017.

[41] L. Xu, C. Lu, Y. Xu, and J. Jia, "Image smoothing via $L_0$ gradient minimization," *ACM Trans. Graph.*, vol. 30, no. 6, pp. 174:1-174:12, 2011.

[42] L. Xu, S. Zheng, and J. Jia, "Unnatural $L_0$ sparse representation for natural image blurring," in *Proc. IEEE Conf. Comput. Vis. Pattern Recognit.*, Jun. 2013, pp. 1107-1114.

[43] H. Xue, S. Zhang, and D. Cai, "Depth image inpainting: Improving low rank matrix completion with low gradient regularization," *IEEE Trans. Image Process.*, vol. 26, no. 9, pp. 4311-4320, Sep. 2017.

[44] M. Wang, Q. Wang, and J. Chanussot, "$L_0$ gradient regularized low-rank tensor model for hyperspectral image denoising," in *Proc. 10th Workshop Hyperspectral Imag. Signal Process., Evol. Remote Sens. (WHISPERS)*, Sep. 2019, pp. 1-6.

[45] F. Xiong, J. Zhou, and Y. Qian, "Hyperspectral restoration via $L_0$ gradient regularized low-rank tensor factorization," *IEEE Trans. Geosci. Remote Sens.*, vol. 57, no. 12, pp. 10410-10425, Dec. 2019.

[46] S. Ono, "$L_0$ gradient projection," *IEEE Trans. Image Process.*, vol. 26, no. 4, pp. 1554-1564, Apr. 2017.

[47] Y.-B. Zheng, T.-Z. Huang, X.-L. Zhao, T.-X. Jiang, T.-Y. Ji, and T.-H. Ma, "Tensor N-tubal rank and its convex relaxation for low-rank tensor recovery," 2018, *arXiv:1812.00688*. [Online]. Available: http://arxiv.org/abs/1812.00688

[48] S. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein, "Distributed optimization and statistical learning via the alternating direction method of multipliers," *Found. Trends Mach. Learn.*, vol. 3, no. 1, pp. 1-122, 2010.

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

[50] C. C. Paige and M. A. Saunders, "LSQR: An algorithm for sparse linear equations and sparse least squares," *ACM Trans. Math. Softw.*, vol. 8, no. 1, pp. 43-71, Mar. 1982.

[51] R. Barrett et al., *Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods*, vol. 43. Philadelphia, PA, USA: SIAM, 1994.

[52] D. Hong, N. Yokoya, J. Chanussot, J. Xu, and X. X. Zhu, "Joint and progressive subspace analysis (JPSA) with spatial-spectral manifold alignment for semisupervised hyperspectral dimensionality reduction," *IEEE Trans. Cybern.*, early access, Nov. 11, 2020, doi: 10.1109/TCYB. 2020.3028931.
[53] 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.

[54] A. Gogna, A. Shukla, H. K. Agarwal, and A. Majumdar, "Split Bregman algorithms for sparse/joint-sparse and low-rank signal recovery: Application in compressive hyperspectral imaging," in *Proc. IEEE Int. Conf. Image Process.*, Oct. 2014, pp. 1302-1306.

[55] M. Golbabaee and P. Vandergheynst, "Hyperspectral image compressed sensing via low-rank and joint-sparse matrix recovery," in *Proc. IEEE Int. Conf. Acoust., Speech Signal Process.*, Mar. 2012, pp. 2741-2744.

[56] D. Hong et al., "More diverse means better: Multimodal deep learning meets remote-sensing imagery classification," *IEEE Trans. Geosci. Remote Sens.*, early access, Aug. 24, 2020, doi: 10.1109/TGRS.2020.3016820.

[57] F. Melgani and L. Bruzzone, "Classification of hyperspectral remote sensing images with support vector machines," *IEEE Trans. Geosci. Remote Sens.*, vol. 42, no. 8, pp. 1778-1790, Aug. 2004.

[58] D. Hong, L. Gao, J. Yao, B. Zhang, A. Plaza, and J. Chanussot, "Graph convolutional networks for hyperspectral image classification," *IEEE Trans. Geosci. Remote Sens.*, early access, Aug. 18, 2020, doi: 10.1109/TGRS.2020.3015157.

<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. His research interests include image analysis, hyperspectral remote sensing, data fusion, machine learning, and artificial intelligence. He has been a Visiting Scholar with Stanford University, Stanford, CA, USA; KTH Royal Institute of Technology, Stockholm, Sweden; and National University of Singapore (NUS), Singapore. Since 2013, he has 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 is also an Adjunct Professor with the Aerospace Information Research Institute, Chinese Academy of Sciences, Beijing, China.

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

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

**Minghua Wang** (Graduate Student Member, IEEE) received the B.S. degree from the School of Automation, Harbin Institute of Technology (HIT), Harbin, China, in 2016, where she is pursuing the Ph.D. degree in control science and engineering.

She is also a Visiting Ph.D. Student with the GIPSA-Lab, Grenoble INP, CNRS, Université Grenoble Alpes, Grenoble, France, from 2019 to 2020. Her research interests include hyperspectral image denoising, tensor completion, and low-rank representation.

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

**Danfeng Hong** (Member, IEEE) received the M.Sc. degree (*summa cum laude*) in computer vision from the College of Information Engineering, Qingdao University, Qingdao, China, in 2015, and the Dr.-Ing. degree (*summa cum laude*) in signal processing in Earth observation (SiPEO) from the Technical University of Munich (TUM), Munich, Germany, in 2019.

Since 2015, he worked as a Research Associate with Remote Sensing Technology Institute (IMF), German Aerospace Center (DLR), Weßling, Germany. He is a Research Scientist with IMF, DLR, where leads the Spectral Vision working group. He is also an Adjunct Scientist with the GIPSA-Lab, Grenoble INP, CNRS, Université Grenoble Alpes, Grenoble, France. His research interests include signal/image processing and analysis, hyperspectral remote sensing, machine/deep learning, artificial intelligence, and their applications in Earth vision.

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

**Qiang Wang** (Member, IEEE) received the B.S., M.S., and Ph.D. degrees in control science and engineering from the Harbin Institute of Technology (HIT), Harbin, China, in 1998, 2000, and 2004, respectively.

Since 2008, he has been a Professor with the Department of Control Science and Engineering, HIT. His research interests include hyperspectral image denoising, signal/image processing, multisensor data fusion, wireless sensor networks, and intelligent detection technology.