# Hyperspectral Image Mixed Noise Removal Based on Multidirectional Low-Rank Modeling and Spatial-Spectral Total Variation

Minghua Wang, Qiang Wang, Jocelyn Chanussot, Dan Li

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

To cite this version:

Minghua Wang, Qiang Wang, Jocelyn Chanussot, Dan Li. Hyperspectral Image Mixed Noise Removal Based on Multidirectional Low-Rank Modeling and Spatial-Spectral Total Variation. IEEE Transactions on Geoscience and Remote Sensing, 2021, 59 (1), pp.488-507. 10.1109/TGRS.2020.2993631. hal-03429691

HAL Id: hal-03429691

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

Submitted on 26 Sep 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Hyperspectral Image Mixed Noise Removal Based on Multidirectional Low-Rank Modeling and Spatial-Spectral Total Variation

Minghua Wang<sup>D</sup>, Qiang Wang<sup>D</sup>, Member, IEEE, Jocelyn Chanussot<sup>D</sup>, Fellow, IEEE, and Dan Li

**Abstract**— Conventional low-rank (LR)-based hyperspectral image (HSI) denoising models generally convert high-dimensional data into 2-D matrices or just treat this type of data as 3-D tensors. However, these pure LR or tensor low-rank (TLR)-based methods lack flexibility for considering different correlation information from different HSI directions, which leads to the loss of comprehensive structure information and inherent spatial-spectral relationship. To overcome these shortcomings, we propose a novel multidirectional LR modeling and spatial-spectral total variation (MLR-SSTV) model for removing HSI mixed noise. By incorporating the weighted nuclear norm, we obtain the weighted sum of weighted nuclear norm minimization (WSWNNM) and the weighted sum of weighted tensor nuclear norm minimization (WSWTNNM) to estimate the more accurate LR tensor, especially, to remove the dead-line noise better. Gaussian noise is further denoised and the local spatial-spectral smoothness is preserved effectively by SSTV regularization. We develop an efficient algorithm for solving the derived optimization based on the alternating direction method of multipliers (ADMM). Extensive experiments on both synthetic data and real data demonstrate the superior performance of the proposed MLR-SSTV model for HSI mixed noise removal.

**Index Terms**— Hyperspectral image (HSI) denoising, multidirectional low-rank (MLR) modeling, spatial-spectral total variation (SSTV), weighted sum of weighted nuclear norm minimization (WSWNNM), weighted sum of weighted tensor nuclear norm minimization (WSWTNNM).

## I. INTRODUCTION

HYPERSPECTRAL images (HSIs) are composed of hundreds of contiguous narrow spectral bands for the same

Manuscript received September 13, 2019; revised February 13, 2020 and April 7, 2020; accepted April 26, 2020. This work was supported in part by the China Scholarship Council (CSC), in part by the Delegation Generale de l'Armement (Project ANR-DGA APHYPIS) under Grant ANR-16 ASTR-0027-01, and in part by the National Natural Science Foundation of China under Grant 61876054 and Grant 61801214. (Corresponding author: Minghua Wang.)

Minghua Wang is with the Department of Control Science and Engineering, Harbin Institute of Technology, Harbin 150001, China, and also with the Université Grenoble Alpes, CNRS, Grenoble INP, GIPSA-lab, Grenoble 38000, France (e-mail: minghuawang1993@163.com).

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

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

Dan Li is with the College of Astronautics, Nanjing University of Aeronautics and Astronautics, Nanjing 260001, China (e-mail: lidanhit@163.com).

Color versions of one or more of the figures in this article are available online at http://ieeexplore.ieee.org.

Digital Object Identifier 10.1109/TGRS.2020.2993631

area on the earth surface. These images have been widely used in various applications [1]–[3], such as agriculture, urban planning, environment monitoring, and so on. However, HSIs are unavoidably corrupted by several kinds of noise, for example, Gaussian noise, salt and pepper noise, stripes, or dead-line noise, which severely cause the degradation of HSIs and hinder the effectiveness of subsequent processing analysis including HSI dimensionality reduction [4], super resolution [5], [6], unmixing [7], [8], classification [9], [10], feature learning [11], [12], and target detection [13]–[16]. Therefore, HSI denoising has become an important and challenging research topic.

In the last decades, various HSI denoising models and methods have been proposed [17]–[21]. A 3-D HSI containing hundreds of bands is regarded as an extension of an RGB image and the traditional 2-D denoising algorithms are used into the HSI band by band [22]–[24]. Green et al. [22] used the maximum noise fraction (MNF) transformation to remove the noise of each band. The transformation of image denoising is gradually improved. Elad and Aharon [23] proposed the approach based on sparse and redundant representations over trained dictionaries. Dabov et al. [24] proposed the denoising strategy based on an enhanced sparse representation in the transform domain. Nevertheless, this kind of processing approach ignores the high correlations of spectral bands, which leads to low-quality recovered results.

Researchers have managed to exploit the low-rank (LR) property of HSI. For example, Hong et al. [25] utilized the LR attributes and proposed a milestone HSI reconstruction algorithm, achieving the best performance in spectral unmixing against spectral variabilities. The high-dimensional data $\mathcal{X} \in \mathbb{R}^{h \times d \times z}$ is reformatted as a Casorati matrix $\mathbf{X} \in \mathbb{R}^{h \times v \times z}$ (a matrix whose columns are vectorized bands from the HSI) [26] where each column represents a vectorized HSI spectral band. $\mathbf{X}$ is represented as $\mathbf{X} = \mathbf{A}\mathbf{D}$, where $\mathbf{A} \in \mathbb{R}^{h \times v \times z}$ is the abundance matrix and $\mathbf{D} \in \mathbb{R}^{r \times z}$ is the endmember matrix. In the general case, the number endmembers $r$ is far smaller than $hv$ or $z$, that is, $r \ll \min(hv, z)$ [27]–[29]. Based on the LR property of hyperspectral imagery, LR matrix approximation and LR matrix factorization based methods [26], [27], [30] have been successfully used in HSI denoising. In [29], LR matrix recovery (LRMR) method was proposed to separate the LR matrix $\mathbf{X}$ from the noisy image $\mathbf{T}$, where $\mathbf{T} = \mathbf{X} + \mathbf{S}$ and $\mathbf{S}$ denotes the unknown sparse noise part. To give a better approximation of the LR part, Xie et al. [30] introduced a weighted schatten $p$-norm into LRMA model and Chen et al. [31] designed a nonconvex norm
for a tighter LR approximation. Xu et al. [32] exploited an LR matrix factorization instead of an LR matrix approximation for avoiding the high computational cost problem caused by calculating the singular value decomposition of large matrices.

To simultaneously consider the spectral correlation and the spatial structure information, He et al. [26] proposed a TV-regularized LR (LRTV) matrix factorization algorithm, which utilized the band-by-band total variation (TV) with the spatial information. Gao et al. [33] proposed a LR representation based denoising method to reduce the influence of noise in mineral mapping. The proposed method named HyDRoS shows high application value to satellite HSI processing and can markedly improve the mineral mapping accuracy. Wu et al. [34] proposed an HSI denoising model named TWNNM, which combined the band-by-band TV regularization and weighted nuclear norm minimization. LRTV and TWNNM methods provided unsatisfactory results in some heavy mixed noisy situations [35] and the bandwise spatial TV regularization led to some undesired “oil painting” effects. For solving this problem, Wu et al. [36] proposed a modified HSI denoising model named STWNNM, whose structure tensor TV explored the gradient information from all spectral bands. However, these approaches only add the spatial smoothness prior and lack the spectral smoothness prior. Aggarwal and Majumdar [37] proposed the spatial-spectral total variation (SSTV) model, which extended the traditional TV to explore both the spatial and the spectral information. This model performed better than the other techniques in mixed noise cases and [38] gave the reasonable explanation for the superior performance of this regularization. But SSTV model ignores the essential global LR property of HSI. Wang et al. [38] combined the LR model with SSTV and proposed the LSSTV method. Although the spectral correlation has been considered, high-dimensional HSI data are converted into a 2-D matrix, which will cause the loss of useful multiorder information.

Recently, many tensor LR-based approaches have been presented by considering an HSI as a third-order tensor. Tensor decomposition (TD) and CANDECOMP/PARAFAC decomposition (CPD) are two common methods. Renard et al. [20] proposed a LR tensor approximation (LRTA) method based on Tucker decomposition to produce a denoised result. Wang et al. [35] used Tucker decomposition and TV to remove mixed noise. Another TD approach was based on the global correlation along spectrum (GCS) and nonlocal self-similarity across space [39]. However, it is hard to predefine three ranks along three modes before the algorithm starts. Liu et al. [40] first proposed CPD for HSI denoising, but the CP rank is difficult to estimate. Xue et al. [41] proposed a nonlocal LR regularized CP TD and a new manner of rank determination.

Based on tensor tubal rank, Zhang et al. [42] proposed the tensor nuclear norm (TNN)-based tensor robust principal component analysis (TRPCA) models and Lu et al. [43] further deduce the TNN for TRPCA model, which is effective to remove sparse noise, but, is intractable for Gaussian noise. Fan et al. [44] introduced a Frobenius norm term to model Gaussian noise and proposed an HSI denoising method based on low-rank tensor recovery (LRTR). To better use the spatial spectral smooth structure of the datas, Chen et al. [45] directly

<figure><img src="image_3.png" /><figcaption>Fig. 1. Illustration of three unfoldings of the third-order Pavia city HSI tensor and corresponding singular values of the unfolding matrices.</figcaption></figure>

combined TV regularizations with TNN-based LR approximation and proposed two tensor low-rank (TLR)-TV algorithms. They just have a good performance for mixed noise including Gaussian noise and salt and pepper noise. Chang et al. [46] proposed hyper-Laplacian regularized unidirectional LR tensor recovery method (LLRT). Although LLRT is an efficient way for HSI Gaussian noise removal, it fails to exploit the global spectral information and to remove the mixed noise. He et al. [47] proposed a unified HSI denoising paradigm to combine the spatial nonlocal similarity with global spectral LR property simultaneously. Fan et al. [48] integrated spatial spectral TV into LR tensor factorization to remove mixed noise.

Recent studies pay more attention to the LR property of a 2-D HSI matrix or 3-D tensor. However, there are still two shortcomings for these pure LR or TLR approaches: 1) due to dead-line noise located at the same position in some bands, pure LR or TLR methods treat this kind of noise as belonging to the LR part and hence retain it and 2) both LR and TLR methods are not capable of getting rid of Gaussian noise, since this kind of noise has an independent distribution. To alleviate the first limitation, the LR property of three unfolding matrices along three directions can be considered [35], [39], [46], as illustrated in Fig. 1. This can be directly expressed as the weighted sum of the nuclear norm (WSNN) of the mode-k matrix. The TLR property along three directions has also considered in Fig. 2 and Zheng et al. [49] proposed the weighted sum of TNN (WSTNN). These two methods are effective to remove the salt and pepper noise in simulated experiments [49]. Simultaneously, dead-line noise will not be regarded as LR features if LR or TLR approaches are applied in the other two modes. But, actually, some dead-line noise still exists since standard nuclear norm or TNN regularize each singular value equally. We introduce a weighted nuclear norm and a weighted tensor nuclear norm into WSNN and WSTNN models to promote the low rankness of the solution and completely remove dead-line noise. To mitigate the second flaw, additional structure information should be involved. TV is recognized as one of the best tool for the removal of Gaussian noise. SSTV regularization is hence considered to exploit both
<figure><img src="image_4.png" /><figcaption>Fig. 2. Illustration of the TLR property along three directions. (a) Third-order Pavia city $\mathcal{X} \in \mathbb{R}^{128 \times 128 \times 80}$. (b) $\hat{\mathcal{X}}^k$ generated by fft($\mathcal{X}^k$, 3). (c) Singular values of the first frontal slice of $\hat{\mathcal{X}}^k$. (d) Singular values from the second to the last frontal slice of $\hat{\mathcal{X}}^k$.</figcaption></figure>

the spatial and the spectral smoothness prior. In this article, a novel HSI denoising model, termed multidirectional LR modeling and spatial-spectral total variation (MLR-SSTV), is proposed. Compared with the existing techniques, our main contributions can be summarized as follows.

1) The multidirectional LR constraints are studied and injected into RPCA and TRPCA for exploiting the essential structure information of the 3-D HSI. Weighted nuclear norm is used into MLR models for treating each singular value adaptively, which can improve the capability and flexibility of the weighted sum of weighted nuclear norm minimization (WSWNNM) and the weighted sum of weighted tensor nuclear norm minimization (WSWTNNM) models. These models can use global LR property from three directions to get rid of dead-lines and estimate the uncorrupted HSI efficiently.

2) To capture both the local spatial and the spectral information, the SSTV regularization is incorporated into the multidirectional LR models. This regularization can simultaneously preserve the spatial-spectral detail smoothness information, and remove Gaussian noise.

3) The multidirectional LR constraints and local spatial-spectral regularization are formulated into a consistent framework in MLR-SSTV. We develop an efficient numerical framework based on the alternating direction method of multipliers (ADMM) to solve the MLR-SSTV. Both the theoretical analysis and experimental results demonstrate the superior denoising performances of our proposed method.

The remainder of this article is organized as follows. In Section II, we briefly introduce some notations and preliminaries of tensors. Then, the MLR-SSTV model is described in Section III. The experimental results and analysis are provided in Section IV. Finally, the conclusions are drawn in Section V.

## II. NOTATIONS AND PRELIMINARIES

In this article, scalars, vectors, and matrices are denoted by lowercase letters, $x$; boldface lowercase letters, $\mathbf{x}$; and boldface capital letters, $\mathbf{X}$, respectively. Tensors are treated as multi-index arrays [50], which are denoted by Euler script letters, $\mathcal{X}$. An $m$-way real-valued tensor is denoted by $\mathbb{R}^{n_1 \times n_2 \times \dots \times n_m}$ and its order $m$ is defined as the number of its dimensions or modes. 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. Some norms of vectors, matrices, and tensors are used. The $l_1$ norm is written as $||\mathcal{X}||_1 = \sum_{h_i, v_i, z_i} |\mathcal{X}_{h_i, v_i, z_i}|$ and the Frobenius norm is $||\mathcal{X}||_F = (\sum_{h_i, v_i, z_i} |\mathcal{X}_{h_i, v_i, z_i}|^2)^{1/2}$. The nuclear norm of a matrix is $||\mathbf{X}||_*$ = $\sum_i \sigma_i(\mathbf{X})$ where $\sigma_i(\mathbf{X})$ is the singular values of matrix $\mathbf{X}$ and the $l_2$ norm of a vector is $||\mathbf{x}||_2 = (\sum_i |\mathbf{x}_i|^2)^{1/2}$. $\hat{\mathcal{X}}$ is the result of discrete Fourier transformation of $\mathcal{X}$ along third dimension, that is $\hat{\mathcal{X}} = \text{fft}(\mathcal{X}, [], 3)$. Naturally, $\mathcal{X}$ can be computed via $\mathcal{X} = \text{ifft}(\hat{\mathcal{X}}, [], 3)$.

**Definition 1 (Two Tensor Product [51]):** Given two three-way tensors $\mathcal{A} \in \mathbb{R}^{n_1 \times n_2 \times n_3}$ and $\mathcal{B} \in \mathbb{R}^{n_2 \times n_4 \times n_3}$, the result $\mathcal{C} \in \mathbb{R}^{n_1 \times n_4 \times n_3}$ of the product of $\mathcal{A}$ and $\mathcal{B}$ is

$$ \mathcal{C}(i, k, \cdot) = \mathcal{A} * \mathcal{B} = \sum_{j=1}^{n_2} \mathcal{A}(i, j, \cdot) \mathcal{B}(j, k, \cdot). \quad (1) $$

**Definition 2 (Special Tensors [42]):** Conjugate Transpose: For a three-way tensor $\mathcal{X} \in \mathbb{R}^{h \times v \times z}$, its conjugate transpose is the tensor conj($\mathcal{X}$) = $\mathcal{X}^* \in \mathbb{R}^{v \times h \times z}$ obtained by conjugate transposing each of the front slices and reversing the order of transposed frontal 2 through $z$.

**Identity Tensor:** The identity tensor $\mathcal{I} \in \mathbb{R}^{h \times v \times z}$ is the tensor whose first frontal slice is the identity matrix and other frontal slices are all zeros.

**Orthogonal Tensor:** A three-way tensor $\mathcal{Q}$ is orthogonal if it satisfies $\mathcal{Q}^* * \mathcal{Q} = \mathcal{Q} * \mathcal{Q}^* = \mathcal{I}$.
**F-Diagonal Tensor:** A three-way tensor $\mathcal{S}$ is f-diagonal if its slices are all diagonal matrices.

**Definition 3 (Mode-k Matricization/Unfolding [52], [53]):** This operator noted unfold($\mathcal{X}$, k) reorders the elements of a tensor $\mathcal{X}$ into a matrix $\mathbf{X}_{(k)}$. As an example of HSI, Fig. 1 shows a 3-D Pavia city HSI data (128 × 128 × 80), together with its three unfoldings, $\mathbf{X}_{(1)} \in \mathbb{R}^{h \times v \times z}$, $\mathbf{X}_{(2)} \in \mathbb{R}^{v \times z \times h}$ and $\mathbf{X}_{(3)} \in \mathbb{R}^{z \times h \times v}$. Inversely, fold($\mathbf{X}_{(k)}$, k) denotes the folding of the matrix into a tensor.

**Definition 4 (Tensor Trace Norm [54]):** It is the sum of the nuclear norm (SNN) of the mode-k unfolding matrix for a three-way HSI tensor

$$||\mathcal{X}||_{\text{SNN}} := \sum_{k=1}^{3} \alpha_k ||\mathbf{X}_{(k)}||_{*} \quad (2)$$

where weighted constants $\alpha_k$ satisfy $\alpha_k \ge 0(k = 1, 2, 3)$ and $\sum_{k=1}^3 \alpha_k = 1$. Due to the existence of the weighted constants $\alpha_k$, we use the WSNN to represent the SNN or tensor trace norm.

**Definition 5 (TNN [43]):** Let $\mathcal{X} = \mathcal{U} * \mathcal{S} * \mathcal{V}^*$ be the t-SVD of $\mathcal{X} \in \mathbb{R}^{h \times v \times z}$, TNN is the sum of singular values of $\mathcal{X}$, that is

$$||\mathcal{X}||_{*} := \sum_{k=1}^{z} S(k, k, 1) \quad (3)$$

and also can be expressed as the SNN of all the frontal slices of $\hat{\mathcal{X}}$

$$||\mathcal{X}||_{*} := \sum_{k=1}^{z} ||\hat{\mathcal{X}}(\cdot, \cdot, k)||_{*} \quad (4)$$

**Theorem 1 (t-SVD [43]):** $\mathcal{X} \in \mathbb{R}^{h \times v \times z}$ can be factorized as

$$\mathcal{X} = \mathcal{U} * \mathcal{S} * \mathcal{V}^* \quad (5)$$

where $\mathcal{U} \in \mathbb{R}^{h \times v \times z}$, $\mathcal{V} \in \mathbb{R}^{v \times v \times z}$ are orthogonal tensors and $\mathcal{S} \in \mathbb{R}^{h \times v \times z}$ is an f-diagonal tensor. The details of t-SVD are described in Algorithm 1.

**Algorithm 1 t-SVD**

**Input:** $\mathcal{X} \in \mathbb{R}^{h \times v \times z}$

1: $\hat{\mathcal{X}} = \text{fft}(\mathcal{X}, \ll, 3)$;

2: **for** $i = 0, 1, \dots, \lfloor \frac{z+1}{2} \rfloor$

3: $[\hat{\mathcal{U}}(\cdot, \cdot, i), \hat{\mathcal{S}}(\cdot, \cdot, i), \hat{\mathcal{V}}(\cdot, \cdot, i)] = \text{SVD}(\hat{\mathcal{X}}(\cdot, \cdot, i))$;

4: **end for**

5: **for** $i = \lfloor \frac{z+1}{2} \rfloor + 1, \dots, z$

6: $\hat{\mathcal{U}}(\cdot, \cdot, i) = \text{conj}(\hat{\mathcal{U}}(\cdot, \cdot, z - i + 2))$;

7: $\hat{\mathcal{S}}(\cdot, \cdot, i) = (\hat{\mathcal{S}}(\cdot, \cdot, z - i + 2))$

8: $\hat{\mathcal{V}}(\cdot, \cdot, i) = \text{conj}(\hat{\mathcal{V}}(\cdot, \cdot, z - i + 2))$;

9: **end for**

10: $\mathcal{U} = \text{ifft}(\hat{\mathcal{U}}, \ll, 3)$, $\mathcal{S} = \text{ifft}(\hat{\mathcal{S}}, \ll, 3)$, $\mathcal{V} = \text{ifft}(\hat{\mathcal{V}}, \ll, 3)$;

**Output:** $\mathcal{U}, \mathcal{S}, \mathcal{V}$.

**Definition 6 (Mode-k Permutation of a Tensor [49]):** As Fig. 3 shows, for a three tensor $\mathcal{X} \in \mathbb{R}^{h \times v \times z}$, this operator noted by $\mathcal{X}^k = \text{permutation}(\mathcal{X}, k)$ changes its permutation order of $(h, v, z)$ with k times, i.e., $\mathcal{X}^1 \in \mathbb{R}^{v \times z \times h}$, $\mathcal{X}^2 \in$

<figure><img src="image_5.png" /><figcaption>Fig. 3. Illustration of mode-k permutation of a tensor $\mathcal{X} \in \mathbb{R}^{h \times v \times z}$.</figcaption></figure>

$\mathbb{R}^{z \times h \times v}$, $\mathcal{X}^3 \in \mathbb{R}^{h \times v \times z}$. The inverse operator is defined as $\mathcal{X} = \text{ipermutation}(\mathcal{X}^k, k)$.

**Definition 7 (WSTNN [49]):** It is the WSTNN of each mode-k unfolding tensor for a three-way HSI tensor

$$||\mathcal{X}||_{\text{WSTNN}} := \sum_{k=1}^{3} \alpha_k ||\mathcal{X}^k||_{*} \quad (6)$$

where weighted constants $\alpha_k$ satisfy $\alpha_k \ge 0(k = 1, 2, 3)$ and $\sum_{k=1}^3 \alpha_k = 1$. As is shown in Fig. 2, the global LR property along each direction can be observed visually.

### III. PROPOSED MODEL

#### A. MLR-SSTV Model

In general, HSIs are assumed to be contaminated by mixed noise and the degradation model can be formulated as $\mathcal{T} = \mathcal{X} + \mathcal{S} + \mathcal{N}$, where $\mathcal{T}$, $\mathcal{X}$, $\mathcal{S}$, and $\mathcal{N}$ have the same size $h \times v \times z$ and correspond to the observed noisy HSI, the latent clean HSI, the sparse noise, and the Gaussian noise, respectively.

By more comprehensively using the global LR property along three directions of HSI, the multidirectional LR (MLR) model is summarized as

$$ \begin{gathered} \min_{\mathcal{X}, \mathcal{S}, \mathcal{N}} \sum_{k=1}^{3} \alpha_k E_k(\mathcal{X}) + \lambda \|\mathcal{S}\|_1 + \mu \|\mathcal{N}\|_F^2 \\ \text{s.t. } \mathcal{T} = \mathcal{X} + \mathcal{S} + \mathcal{N} \end{gathered} \quad (7) $$

where the function $E_k(\mathcal{X})$ is set either as $||\mathbf{X}_{(k)}||_{*}$ when matrix form WSNN is used or as $||\mathcal{X}^k||_{*}$ when tensor form WSTNN is used.

However, it is not fair to treat each singular value equally during the minimization process since larger singular values usually contain most of the information, such as image edges and thus, these values should be shrank less than the smaller ones [55]. To improve the flexibility of the above standard forms, weighted nuclear norm is introduced into the MLR model

$$ \begin{gathered} \min_{\mathcal{X}, \mathcal{S}, \mathcal{N}} \sum_{k=1}^{3} \alpha_k E_k(\mathcal{X})_{\omega} + \lambda \|\mathcal{S}\|_1 + \mu \|\mathcal{N}\|_F^2 \\ \text{s.t. } \mathcal{T} = \mathcal{X} + \mathcal{S} + \mathcal{N} \end{gathered} \quad (8) $$

where the function $E_k(\mathcal{X})_{\omega} = ||\mathbf{X}_{(k)}||_{\omega,*}$ is used for the matrix form WSWNN and $E_k(\mathcal{X})_{\omega} = ||\mathcal{X}^k||_{\omega,*}$ is used for the tensor form WSWTNN. A test is used to point out that weighted nuclear norm can solve the overshrunk problem of WSNN and WSTNN to remove dead-line noise completely, which is illustrated in Section IV-A.

The weighted TNN of $\mathcal{X}$, denoted as $||\mathcal{X}^k||_{\omega,*}$, is the sum of the WNNM of all the frontal slices of $\mathcal{X}$, that is, $||\mathcal{X}^k||_{\omega,*} = \sum_{k=1}^z ||\hat{\mathcal{X}}(\cdot, \cdot, k)||_{\omega,*} = \sum_{k=1}^z \sum_{i=1}^{\min(h,v)} \omega(i, i, k)\sigma(i, i, k)$,
where $\sigma(i, i, k)$ is the singular value of each frontal slice of $\hat{X}$ and $\omega$ is the weighting parameter.

The MLR model for HSI denoising only considers the different LR properties from different directions, which may lead to a little residual Gaussian noise. To solve this problem, the TV regularization is incorporated into our model. Considering both the spatial and the spectral information, we introduce the SSTV into the MLR model. The SSTV is defined as follows [37]:

$$||\mathcal{X}||_{\text{SSTV}} = ||D_h \mathcal{X} D_z||_1 + ||D_v \mathcal{X} D_z||_1 \quad (9)$$

where $D_h$ and $D_v$ are two differential operators along the horizontal and vertical directions. $D_z$ is a 1-D finite differencing operator, which is used to extract the spectral signature of each pixel.

The MLR framework is useful to exploit the global LR property of HSIs from three directions and SSTV is a local model for spatial-spectral texture extraction and preservation. Combining these two complementary models, we propose the novel following MLR-SSTV model:

$$\begin{gathered} \min_{\mathcal{X}, \mathcal{S}, \mathcal{N}} \sum_{k=1}^{3} \alpha_k E_k(\mathcal{X})_{\omega} + \lambda \|\mathcal{S}\|_1 + \lambda_{\text{tv}} \|\mathcal{X}\|_{\text{SSTV}} + \mu \|\mathcal{N}\|_F^2 \\ \text{s.t. } \mathcal{T} = \mathcal{X} + \mathcal{S} + \mathcal{N} \end{gathered} \quad (10)$$

where $\alpha_k$, $\lambda$, $\mu$, and $\lambda_{tv}$ are four positive constants to control the tradeoff among four terms. The first term $E_k(\mathcal{X})_\omega$ corresponds to WSWNN or WSWTNN model based on the multidirection LR property to separate the clean HSI from the noisy one. The second term $\|\mathcal{S}\|_1$ is used to find the sparse noise, including salt and pepper noise and dead-line noise. The $||\mathcal{X}||_{\text{SSTV}}$ term is utilized to characterize the smooth structure of spatial-spectral domains for removing Gaussian noise and enhancing the local detail information. The last Frobenius norm term $||\mathcal{N}||_F^2$ is further used to remove heavy Gaussian noise. The block scheme of the proposed MLR-SSTV model is shown in Fig. 4.

## B. Optimization Procedure

By introducing three auxiliary variables $\mathcal{M}_k$ with $k = 1, 2, 3, \mathcal{P}$ and $\mathcal{Q}$, we first rewrite the model from (10) into the equivalent function

$$\begin{gathered} \min_{\mathcal{X}, \mathcal{M}_k, \mathcal{S}, \mathcal{N}} \sum_{k=1}^{3} \alpha_k E_k(\mathcal{M}_k)_{\omega} \\ + \lambda \|\mathcal{S}\|_1 + \lambda_{\text{tv}} (\|\mathcal{P}\|_1 + \|\mathcal{Q}\|_1) + \mu \|\mathcal{N}\|_F^2 \\ \text{s.t. } \mathcal{M}_1 = \mathcal{M}_2 = \mathcal{M}_3 = \mathcal{X}, \mathcal{P} = D_h \mathcal{X} D_z, \mathcal{Q} = D_v \mathcal{X} D_z \\ \mathcal{T} = \mathcal{X} + \mathcal{S} + \mathcal{N}. \end{gathered} \quad (11)$$

Based on the ADMM [56], [57], the unconstrained augmented Lagrangian function of (11) is formulated as

$$\begin{gathered} \min_{\mathcal{X}, \mathcal{M}_k, \mathcal{S}, \mathcal{N}} \sum_{k=1}^{m} \alpha_k E_k(\mathcal{M}_k)_{\omega} + \lambda \|\mathcal{S}\|_1 + \lambda_{\text{tv}} (\|\mathcal{P}\|_1 + \|\mathcal{Q}\|_1) \\ + \mu \|\mathcal{N}\|_F^2 + \sum_{k=1}^{3} \left( \langle \mathcal{Y}_{1k}, \mathcal{X}_k - \mathcal{M}_k \rangle + \frac{\beta}{2} \|\mathcal{X}_k - \mathcal{M}_k\|_F^2 \right) \end{gathered}$$

$$\begin{aligned} & + \langle \mathcal{Y}_2, D_h \mathcal{X} D_z - \mathcal{P} \rangle + \frac{\beta}{2} \|D_h \mathcal{X} D_z - \mathcal{P}\|_F^2 \\ & + \langle \mathcal{Y}_3, D_v \mathcal{X} D_z - \mathcal{Q} \rangle + \frac{\beta}{2} \|D_v \mathcal{X} D_z - \mathcal{Q}\|_F^2 \\ & + \langle \mathcal{Y}_4, \mathcal{X} + \mathcal{S} + \mathcal{N} - \mathcal{T} \rangle + \frac{\beta}{2} \|\mathcal{X} + \mathcal{S} + \mathcal{N} - \mathcal{T}\|_F^2 \end{aligned} \quad (12)$$

where $\mathcal{Y}_{11}, \mathcal{Y}_{12}, \mathcal{Y}_{13}, \mathcal{Y}_2, \mathcal{Y}_3$ and $\mathcal{Y}_4$ are six multipliers used in the ADMM iterations. Equation (12) can be alternatively optimized over one variable when fixing the others. The whole procedure is divided into five parts.

1) Update $\mathcal{M}_k, k = 1, 2, 3$ by minimizing

$$\min_{\mathcal{M}_k} \frac{\alpha_k}{\beta} E_k(\mathcal{M}_k)_{\omega} + \frac{1}{2} \|\mathcal{M}_k - \left(\mathcal{X} + \frac{\mathcal{Y}_{1k}}{\beta}\right)\|_F^2. \quad (13)$$

For WSWNN form, (13) can be rewritten as

$$\min_{\mathcal{M}_k} \frac{\alpha_k}{\beta} \|\mathbf{M}_{(k)}\|_{\omega, *} + \frac{1}{2} \|\mathcal{M}_k - \left(\mathcal{X} + \frac{\mathcal{Y}_{1k}}{\beta}\right)\|_F^2. \quad (14)$$

$\mathbf{M}_{(k)}$ is actually a WNNM model and has led to a closed form, which is proven in [58]. The weighted shrinking operator $\mathbf{R}_{(\alpha_k/\beta)}(\mathcal{X} + \mathcal{Y}_1^k) = \text{diag}((\Sigma(i, i) - (\alpha_k/\beta)\omega_i)_+) \text{ is applied to the singular value matrix } \Sigma \text{ of } \mathbf{SVD}(\text{unfold}(\mathcal{X} + (\mathcal{Y}_{1k}/\beta), k)) = [\mathbf{U}, \Sigma, \mathbf{V}], \text{ where } \omega_i \text{ is the weight of WNNM calculated by } \omega_i = 1/(\Sigma(i, i) + \epsilon), \epsilon > 0 \text{ is a small positive value. The estimated } \mathbf{M}_{(k)} \text{ and } \mathcal{M}_k \text{ are expressed as}$

$$\mathbf{M}_{(k)} = \mathbf{U} \mathbf{R}_{\frac{\alpha_k}{\beta}} \left( \mathcal{X} + \frac{\mathcal{Y}_{1k}}{\beta} \right) \mathbf{V}^T, \mathcal{M}_k = \text{fold}(\mathbf{M}_{(k)}, k). \quad (15)$$

## Algorithm 2 Weighted Tensor Singular Value Thresholding

Input: $\mathcal{L}_k \in \mathbb{R}^{h \times v \times z}$, $\frac{\alpha_k}{\beta}$

1: $\hat{\mathcal{L}}_k = \text{fft}(\mathcal{L}_k, [1, 3])$;

2: **for** $i = 0, 1, \dots, \lfloor \frac{z+1}{2} \rfloor$

3: $[\hat{\mathcal{U}}(:, : , i), \hat{\mathcal{S}}(:, : , i), \hat{\mathcal{V}}(:, : , i)] = \text{SVD}(\hat{\mathcal{L}}_k(:, : , i));$

4: $\omega = \frac{1}{\text{diag}(\hat{\mathcal{S}}(:, : , i)) + \epsilon}$;

5: $\hat{\mathcal{S}}_{\text{new}}(:, : , i) = \max(\hat{\mathcal{S}}(:, : , i) - \frac{\alpha_k}{\beta}\omega, 0)$;

6: **end for**

7: **for** $i = [\lfloor \frac{z+1}{2} \rfloor + 1], \dots, z$

8: $\hat{\mathcal{U}}(:, : , i) = \text{conj}(\hat{\mathcal{U}}(:, : , z - i + 2));$

9: $\hat{\mathcal{S}}_{\text{new}}(:, : , i) = \text{conj}(\hat{\mathcal{S}}_{\text{new}}(:, : , i));$

10: $\hat{\mathcal{V}}(:, : , i) = \text{conj}(\hat{\mathcal{V}}(:, : , z - i + 2));$

11: **end for**

12: $\mathcal{U} = \text{ifft}(\hat{\mathcal{U}}, [1, 3]), \mathcal{S}_{\text{new}} = \text{ifft}(\hat{\mathcal{S}}_{\text{new}}, [1, 3]), \mathcal{V} = \text{ifft}(\hat{\mathcal{V}}, [1, 3]),$

13: $\mathcal{M}_k^k = \mathcal{U} * \mathcal{S}_{\text{new}} * \mathcal{V}^*$;

14: $\mathcal{M}_k = \text{ipermutation}(\mathcal{M}_k^k, k)$;

Output: $\mathcal{M}_k$.

For WSWTNN form, (13) can be transformed as

$$\min_{\mathcal{M}_k} \frac{\alpha_k}{\beta} \|\mathcal{M}_k^k\|_{\omega, *} + \frac{1}{2} \|\mathcal{M}_k - \left(\mathcal{X} + \frac{\mathcal{Y}_{1k}}{\beta}\right)\|_F^2. \quad (16)$$

$\mathcal{M}_k^k$ is the mode-k permutation of tensor $\mathcal{M}_k$. Let $\mathcal{L}_k = \mathcal{X} + \frac{\mathcal{Y}_{1k}}{\beta}$ and the above subproblem can be solved by the following weighted tensor singular value thresholding
<figure><img src="image_6.png" /><figcaption>Fig. 4. HSI tensor denoised model based on MLR- SSTV.</figcaption></figure>

method. According to Algorithm 2, SVD is applied for each front slice matrix of DFT of $\mathcal{X}$.

2) Update $\mathcal{P}$, $\mathcal{Q}$ and $\mathcal{S}$ by minimization. Three subproblems are rewritten as

$$ 
\begin{aligned}
\min_{\mathcal{P}} & \lambda_{\text{tv}} \|\mathcal{P}\|_1 + \frac{\beta}{2} \left\| \mathcal{P} - \left(D_h \mathcal{X} D_z + \frac{\mathcal{Y}_2}{\beta}\right) \right\|_F^2 \\
\min_{\mathcal{Q}} & \lambda_{\text{tv}} \|\mathcal{Q}\|_1 + \frac{\beta}{2} \left\| \mathcal{Q} - \left(D_v \mathcal{X} D_z + \frac{\mathcal{Y}_3}{\beta}\right) \right\|_F^2 \\
\min_{\mathcal{S}} & \lambda_1 \|\mathcal{S}\|_1 + \frac{\beta}{2} \left\| \mathcal{S} - \left(\mathcal{T} - \mathcal{X} - \mathcal{N} - \frac{\mathcal{Y}_4}{\beta}\right) \right\|_F^2
\end{aligned}
\quad (17) $$

By defining the soft-thresholding operator $S_\varepsilon[x] = \text{sgn}(x) \ast \max(|x| - \varepsilon, 0)$ [59], the solutions of three subproblems are obtained as

$$ 
\begin{aligned}
\mathcal{P} &= \mathbf{S}_{\lambda_{\text{tv}}/(\beta)} \left[ D_h \mathcal{X} D_z + \frac{\mathcal{Y}_2}{\beta} \right] \\
\mathcal{Q} &= \mathbf{S}_{\lambda_{\text{tv}}/(\beta)} \left[ D_v \mathcal{X} D_z + \frac{\mathcal{Y}_3}{\beta} \right] \\
\mathcal{S} &= \mathbf{S}_{\lambda/(\beta)} \left[ \mathcal{T} - \mathcal{X} - \mathcal{N} - \frac{\mathcal{Y}_4}{\beta} \right]
\end{aligned}
\quad (18) $$

where sign operator $\text{sgn}(x)$ is defined as

$$ \text{sgn}(x) = \begin{cases} 1, & x > 0 \\ 0, & x = 0 \\ -1, & x < 0. \end{cases} $$

3) Update $\mathcal{N}$. Extract all $\mathcal{N}$-related items from (12)

$$ \min_{\mathcal{N}} \mu \|\mathcal{N}\|_F^2 + \frac{\beta}{2} \left\| \mathcal{T} - \mathcal{X} - \mathcal{S} - \mathcal{N} - \frac{\mathcal{Y}_4}{\beta} \right\|_F^2 \quad (19) $$

and through a simple reduction, the solution is shown as follows:

$$ \mathcal{N} = \frac{\beta(\mathcal{T} - \mathcal{X} - \mathcal{S}) - \mathcal{Y}_4}{\mu + \beta}. \quad (20) $$

4) Update $\mathcal{X}$. The $\mathcal{X}$-related subproblem is rewritten into its vector

$$ 
\begin{aligned}
\mathbf{x} = \operatorname*{arg\,min}_{\mathbf{x}} f(\mathbf{x}) = & \operatorname*{arg\,min}_{\mathbf{x}} \frac{\beta}{2} \left\| 3\mathbf{x} - \sum_{k=1}^{3} \left( m_k - \frac{\mathbf{y}_{1k}}{\beta} \right) \right\|_2^2 \\
& + \frac{\beta}{2} \left\| \mathbf{p} - D_1 \mathbf{x} - \frac{\mathbf{y}_2}{\beta} \right\|_2^2 + \frac{\beta}{2} \left\| \mathbf{q} - D_2 \mathbf{x} - \frac{\mathbf{y}_3}{\beta} \right\|_2^2 \\
& + \frac{\beta}{2} \left\| \mathbf{t} - \mathbf{x} - \mathbf{s} - \mathbf{n} - \frac{\mathbf{y}_4}{\beta} \right\|_2^2
\end{aligned}
\quad (21) $$

where $D_1 = D_h \otimes D_z$, $D_2 = D_v \otimes D_z$, and $\otimes$ denotes the Kronecker product. Then (21) is solved by setting the partial derivative of $f(\mathbf{x})$ as zero, that is

$$ 
\begin{aligned}
\frac{\partial f(\mathbf{x})}{\partial \mathbf{x}} &= \left[ \mathbf{x} - \left( \mathbf{t} - \mathbf{s} - \mathbf{n} - \frac{\mathbf{y}_4}{\beta} \right) \right] + 3\mathbf{x} - \sum_{k=1}^{3} \left( m_k - \frac{\mathbf{y}_{1k}}{\beta} \right) \\
&\quad + D_1^T \left( D_1 \mathbf{x} + \frac{\mathbf{y}_2}{\beta} - \mathbf{p} \right) + D_2^T \left( D_2 \mathbf{x} + \frac{\mathbf{y}_3}{\beta} - \mathbf{q} \right) = 0.
\end{aligned}
\quad (22) $$
This problem can be treated as the following linear system:

$$ (4 + D_1^T D_1 + D_2^T D_2)\mathbf{x} = \left(\mathbf{t} - \mathbf{s} - \mathbf{n} - \frac{\mathbf{y}_4}{\beta}\right) + \sum_{k=1}^3 \left(\mathbf{m}_k - \frac{\mathbf{y}_{1k}}{\beta}\right) + D_1^T \left(\mathbf{p} - \frac{\mathbf{y}_2}{\beta}\right) + D_2^T \left(\mathbf{q} - \frac{\mathbf{y}_3}{\beta}\right). \quad (23) $$

Obviously, (23) can be solved by the preconditioned conjugate gradient (PCG) [60] or the algorithm for sparse linear equations and least squares (LSQR) [61]. The solution $\mathbf{x}$ is rearranged into the recovered tensor $\mathcal{X}$.

Summarizing the aforementioned description, we arrive at the ADMM iteration for the proposed MLR-SSTV in Algorithm 3.

**Algorithm 3 MLR-SSTV for HSI Mixed Noise Removal**

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

1: Initialization: $\mathcal{X} = \mathcal{S} = 0$, $\mathcal{P} = \mathcal{Q} = 0$, $M_k = \mathcal{Y}_{1k} = 0$, $\mathcal{Y}_2 = \mathcal{Y}_3 = \mathcal{Y}_4 = 0$, $\epsilon = 10^{-6}$, maxiter = 100, $i = 0$

2: for $i = 0, 1, \dots, \text{maxiter}$

3: $i = i + 1$;

4: for WSWNN, update $M_k$ by solving Eq.(15):

$$ \mathbf{M}_{(k)} = \mathbf{UR}_{\frac{\mathcal{N}}{\beta}}(\mathcal{X} + \mathcal{Y}_1^k)\mathbf{V}^T, \quad M_k = \text{fold}(\mathbf{M}_{(k)}). $$

5: for WSWTNN, solve Eq.(16) via weighted tensor singular value thresholding method;

6: update $\mathcal{P}$, $\mathcal{Q}$ and $\mathcal{S}$ by solving Eq.(18)

$$ \mathcal{P} = \mathbf{S}_{\lambda_{tv}/(\beta)}[D_h \mathcal{X} D_z + \frac{\mathbf{y}_2}{\beta}] $$

$$ \mathcal{Q} = \mathbf{S}_{\lambda_{tv}/(\beta)}[D_v \mathcal{X} D_z + \frac{\mathbf{y}_3}{\beta}] $$

$$ \mathcal{S} = \mathbf{S}_{\lambda/(\beta)}[\mathcal{T} - \mathcal{X} - \mathcal{N} - \frac{\mathbf{y}_4}{\beta}] $$

7: update $\mathcal{N}$ by solving Eq.(20)

$$ \mathcal{N} = \frac{\beta(\mathcal{T} - \mathcal{X} - \mathcal{S} - \mathcal{Y}_4)}{\mu + \beta} $$

8: update $\mathcal{x}$ by solving Eq.(22) via PDG or LSQR algorithm:

$$ (4 + D_1^T D_1 + D_2^T D_2)\mathbf{x} = \left(\mathbf{t} - \mathbf{s} - \mathbf{n} - \frac{\mathbf{y}_4}{\beta}\right) + \sum_{k=1}^3 \left(\mathbf{m}_k - \frac{\mathbf{y}_{1k}}{\beta}\right) + D_1^T \left(\mathbf{p} - \frac{\mathbf{y}_2}{\beta}\right) + D_2^T \left(\mathbf{q} - \frac{\mathbf{y}_3}{\beta}\right) $$

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

10: update the multipliers

$$ \mathcal{Y}_{1k} = \mathcal{Y}_{1k} + \beta(\mathcal{X}_k - \mathcal{M}_k) $$

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

$$ \mathcal{Y}_3 = \mathcal{Y}_3 + \beta(D_v \mathcal{X} D_z - \mathcal{Q}) $$

$$ \mathcal{Y}_4 = \mathcal{Y}_4 + \beta(\mathcal{X} + \mathcal{S} + \mathcal{N} - \mathcal{T}) $$

11: end for

**Output:** $\mathcal{X}$, $\mathcal{S}$, $\mathcal{N}$

2) LRMR [29] unfolds a 3-D HSI image as a matrix via vectorizing each band as a column and uses the LR property of this 2-D matrix.

3) LRTV method [26] combines LR property and simple TV for mixed noise removal.

4) SSTV [37] extends the traditional TV model for the spatial and spectral correlation.

5) LR constraint and SSTV (LSSTV) [38] integrates LR constraint into the above SSTV model to remove mixed noise.

6) LRTA [20] employs Tucker decomposition and is a parameter-free method.

7) Parallel factor analysis (PARAFAC) [40] reduces the number of estimated rank to one and estimates the optimal rank using the covariance matrix of the n-mode unfolding denoised matrix.

8) TV regularized LR TD (LRTDTV) [35] exploits the global correlation based on Tucker decomposition and uses SSTV to describe the structure information from the spatial and spectral modes.

9) LR tensor recovery (LRTR) [44] formulates the clean HSI as a low-tubal-rank model and simultaneously removes Gaussian noise and sparse noise.

10) LLRT [46] considers unidirectional LR tensor recovery and introduces analysis-based hyper-Laplacian prior into a new spectral TV.

Their implementation codes are obtained from the authors' websites or paperworks.

## A. Simulated HSI Data Experiments

The degraded simulated HSIs are generated from two public HSI data sets including the Pavia City Center data set and the Washington DC Mall data set. The Pavia City Center data set is captured by the Reflective Optics System Imaging Spectrometer (ROSIS) sensor. Since the first 22 bands are heavily destroyed by the noise, these images hardly qualify as the clean reference ones. The bands are selected from 23 to 102 and the size of the selected subimage is 128 × 128 × 80. The Washington DC Mall data set was collected by the hyperspectral digital imagery collection experiment, containing 1209×307×191 pixels. A subregion image 256×256×191 is used in our experiment. In the experiments, the values of HSI data are normalized into [0,1]. To better evaluate the denoised results, we adopt four quality assessments including the mean peak signal-to-noise ratio (MPSNR), mean structural similarity index (MSSIM), the Erreur Relative Globale Adimensionnelle de Synthese (ERGAS), and the mean spectral angle distance (MSAD). MPSNR and MSSIM are denoted by the averages of the PSNR and SSIM values between noisy bands and denoised ones, which pays attention to the recovery precision of spatial pixels

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

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

# IV. EXPERIMENTAL RESULTS AND DISCUSSION

Both simulated and real data experiments are carried out to access the performances of the proposed MLR-SSTV for HSI mixed noise removal. Herein, different state-of-the-art HSI denoising methods selected for comparison include the following.

1) Block matching 4-D (BM4D) filtering [62], as an extension of block matching and 3-D (BM3D) filtering, exploits nonlocal similarity based in the transform domain and removes Gaussian noise efficiently.
<table><caption>TABLE I
ADDED MIXED NOISE IN SIX CASES</caption><thead><tr><th>Noise Case</th><th>PSNR(dB)</th></tr></thead><tbody><tr><td>Case 1</td><td>Gaussian (SNR = 10~20 dB)+salt and pepper(20%) in all bands<br>+ four vertical and horizontal deadlines in each of four random bands</td></tr><tr><td>Case 2</td><td>Gaussian (SNR = 10~20 dB)+salt and pepper(selected from 0 to 20% randomly) in all bands<br>+ dead-lines in 61~70 bands (with the number of stripes randomly selected from 1 to 10<br>and the width of the dead-lines randomly generated from 1 to 3)</td></tr><tr><td>Case 3</td><td>Gaussian (mean SNR = 12.12 dB, and noise variance is calculated by Eq.(28))<br>+salt and pepper(20%) in all bands<br>+ eight vertical and horizontal deadlines in each of eight random bands</td></tr><tr><td>Case 4</td><td>Gaussian (mean SNR = 12.12 dB, and noise variance is calculated by Eq.(28))<br>+salt and pepper(selected from 0 to 20% randomly) in all bands<br>+ dead-lines in 61~70 bands (with the number of stripes randomly selected from 1 to 10<br>and the width of the dead-lines randomly generated from 1 to 3)</td></tr><tr><td>Case 5</td><td>Gaussian (zero-mean, and noise variance is different for different bands and is selected from 0.1~0.2 randomly)<br>+salt and pepper(20%) in all bands+ eight vertical and horizontal deadlines in each of eight random bands</td></tr><tr><td>Case 6</td><td>Gaussian (zero-mean, and noise variance is different for different bands and is selected from 0.1~0.2 randomly)<br>+salt and pepper(selected from 0 to 20% randomly) in all bands<br>+ dead-lines in 61~70 bands (with the number of stripes randomly selected from 1 to 10<br>and the width of the dead-lines randomly generated from 1 to 3)</td></tr></tbody></table>

<table><caption>TABLE II
MPSNR RESULTS OF HSI IMAGE DENOISED BY DIFFERENT METHODS</caption><thead><tr><th rowspan="2">Data set</th><th rowspan="2">Noise case</th><th colspan="10">MPSNR(dB)</th></tr><tr><th>Noise</th><th>BM4D</th><th>LRMR</th><th>LRTV</th><th>SSTV</th><th>LSSTV</th><th>LRTA</th><th>PARAFAC</th><th>LRTDTV</th><th>LRTR</th><th>LLRT</th><th>WSWNN-SSTV</th><th>WSWTNN-SSTV</th></tr></thead><tbody><tr><td rowspan="6">Pavia</td><td>Case 1</td><td>11.29</td><td>23.13</td><td>37.53</td><td>38.82</td><td>39.06</td><td>39.92</td><td>21.38</td><td>20.81</td><td>39.28</td><td>31.96</td><td>23.22</td><td><b>40.20</b></td><td><u>40.10</u></td></tr><tr><td>Case 2</td><td>14.69</td><td>24.46</td><td>38.68</td><td>39.80</td><td>40.09</td><td>40.60</td><td>23.00</td><td>24.91</td><td>40.56</td><td>34.48</td><td>24.23</td><td><b>41.30</b></td><td><u>41.26</u></td></tr><tr><td>Case 3</td><td>11.27</td><td>23.18</td><td>33.93</td><td>36.41</td><td>36.69</td><td>37.25</td><td>21.17</td><td>20.65</td><td>36.47</td><td>29.16</td><td>23.29</td><td><b>38.13</b></td><td><u>37.79</u></td></tr><tr><td>Case 4</td><td>14.52</td><td>24.82</td><td>36.59</td><td>38.22</td><td>38.37</td><td>38.81</td><td>22.96</td><td>24.65</td><td>38.79</td><td>31.74</td><td>25.09</td><td><b>39.27</b></td><td><u>39.15</u></td></tr><tr><td>Case 5</td><td>11.41</td><td>23.38</td><td>31.05</td><td>41.32</td><td>41.58</td><td>42.12</td><td>21.42</td><td>20.91</td><td>40.53</td><td>35.22</td><td>26.09</td><td><u>43.56</u></td><td><b>44.12</b></td></tr><tr><td>Case 6</td><td>12.76</td><td>24.26</td><td>26.76</td><td>30.52</td><td>29.37</td><td>30.39</td><td>22.07</td><td>22.41</td><td>31.17</td><td>24.82</td><td>25.09</td><td><u>31.23</u></td><td><b>32.09</b></td></tr><tr><td rowspan="6">WDC</td><td>Case 1</td><td>10.70</td><td>23.15</td><td>39.13</td><td>45.78</td><td>42.27</td><td>40.43</td><td>20.83</td><td>20.46</td><td>47.07</td><td>36.44</td><td>22.55</td><td><u>49.36</u></td><td><b>50.90</b></td></tr><tr><td>Case 2</td><td>14.66</td><td>20.85</td><td>38.23</td><td>46.45</td><td>44.42</td><td>45.19</td><td>24.91</td><td>26.60</td><td>47.56</td><td>38.61</td><td>25.97</td><td><u>50.31</u></td><td><b>51.61</b></td></tr><tr><td>Case 3</td><td>10.68</td><td>23.46</td><td>39.16</td><td>43.30</td><td>41.11</td><td>39.09</td><td>20.73</td><td>20.40</td><td>45.54</td><td>35.97</td><td>23.51</td><td><u>47.72</u></td><td><b>48.78</b></td></tr><tr><td>Case 4</td><td>14.56</td><td>21.06</td><td>38.27</td><td>45.26</td><td>43.38</td><td>45.79</td><td>24.80</td><td>26.51</td><td>46.61</td><td>38.24</td><td>25.51</td><td><u>48.77</u></td><td><b>49.59</b></td></tr><tr><td>Case 5</td><td>10.10</td><td>20.44</td><td>29.66</td><td>35.41</td><td>29.68</td><td>32.41</td><td>18.62</td><td>18.29</td><td><u>35.46</u></td><td>24.83</td><td>21.14</td><td>32.39</td><td><b>37.05</b></td></tr><tr><td>Case 6</td><td>12.79</td><td>23.16</td><td>29.96</td><td>35.73</td><td>30.57</td><td>33.28</td><td>21.57</td><td>21.72</td><td><u>36.60</u></td><td>25.09</td><td>22.37</td><td>33.13</td><td><b>37.53</b></td></tr></tbody></table>

ERGAS calculates the relative error of all the bands and MSAD is a spectral-based metric

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

$$ \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 (27) $$

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

To simulate mixed noise, Gaussian noise, salt and pepper noise, dead-line noise, and stripes are added to the Pavia City Center and the Washington DC Mall HSI data. The specific

noise configuration is shown in Table I. For cases 1 and 2, Gaussian noise is randomly added to all bands and the range of the SNR is from 10 to 20 dB. For cases 3 and 4, the variance $\sigma_z^2$ of Gaussian noise is defined as [63]

$$ \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 (28) $$

where $\delta$ and $\eta$ are two adjustable parameters and the mean noise SNR is set as 12.12 dB. In other cases, the Gaussian noise with zero mean and different variances randomly sampled from 0 to 0.2 is added into each band. In cases 1, 3, and 5, each band is corrupted by the salt and pepper noise with proportion 20% and in the other cases, salt and pepper noise is randomly selected from 0 to 0.2. In cases 2, 4, and 6, 61~70 bands are heavily corrupted by dead-lines with the 3~10 stripes and the width generated from 1 to 3 while several
<figure><img src="image_7.png" /><figcaption>Fig. 5. Denoised results by all the compared methods. (a) Original band 65 of Pavia City Center. (b) Simulated noise band of case 6. (c) BM4D. (d) LRMR. (e) LRTV. (f) SSTV. (g) LSSTV. (h) LRTA. (i) PARAFAC. (j) LRTR. (k) LRTDTV. (l) LLRT. (m) WSWNN-SSTV. (n) WSWTNN-SSTV.</figcaption></figure>

<table><caption>TABLE III
MSSIM RESULTS OF HSI IMAGE DENOISED BY DIFFERENT METHODS</caption><thead><tr><th rowspan="2">Data set</th><th rowspan="2">Noise case</th><th colspan="11">MSSIM</th></tr><tr><th>Noise</th><th>BM4D</th><th>LRMR</th><th>LRTV</th><th>SSTV</th><th>LSSTV</th><th>LRTA</th><th>PARAFAC</th><th>LRTDTV</th><th>LRTR</th><th>LLRT</th><th>WSWNN-SSTV</th><th>WSWTNN-SSTV</th></tr></thead><tbody><tr><td rowspan="6">Pavia</td><td>Case 1</td><td>0.0969</td><td>0.6157</td><td>0.9563</td><td>0.9691</td><td>0.9762</td><td><u>0.9789</u></td><td>0.5704</td><td>0.4943</td><td>0.9786</td><td>0.8855</td><td>0.6703</td><td><b>0.9820</b></td><td><b>0.9820</b></td></tr><tr><td>Case 2</td><td>0.2558</td><td>0.6210</td><td>0.9494</td><td>0.9756</td><td>0.9811</td><td>0.9820</td><td>0.5373</td><td>0.6682</td><td><u>0.9841</u></td><td>0.9379</td><td>0.6374</td><td><b>0.9860</b></td><td><b>0.9860</b></td></tr><tr><td>Case 3</td><td>0.0949</td><td>0.6108</td><td>0.9121</td><td>0.9370</td><td>0.9602</td><td>0.9637</td><td>0.5414</td><td>0.4797</td><td>0.9404</td><td>0.7974</td><td>0.6415</td><td><b>0.9717</b></td><td><u>0.9702</u></td></tr><tr><td>Case 4</td><td>0.2378</td><td>0.6351</td><td>0.9491</td><td>0.9689</td><td>0.9724</td><td>0.9753</td><td>0.5356</td><td>0.6594</td><td>0.9760</td><td>0.8838</td><td>0.6365</td><td><u>0.9775</u></td><td><b>0.9777</b></td></tr><tr><td>Case 5</td><td>0.1029</td><td>0.6223</td><td>0.9106</td><td>0.9648</td><td>0.9872</td><td>0.9885</td><td>0.5688</td><td>0.4991</td><td>0.9895</td><td>0.9410</td><td>0.7689</td><td><u>0.9930</u></td><td><b>0.9952</b></td></tr><tr><td>Case 6</td><td>0.1192</td><td>0.6400</td><td>0.7123</td><td>0.8552</td><td>0.8097</td><td>0.8444</td><td>0.5185</td><td>0.5462</td><td>0.8746</td><td>0.5752</td><td>0.6984</td><td><u>0.8633</u></td><td><b>0.8910</b></td></tr><tr><td rowspan="6">WDC</td><td>Case 1</td><td>0.0462</td><td>0.3866</td><td>0.8627</td><td>0.9549</td><td>0.9554</td><td>0.9382</td><td>0.3705</td><td>0.3015</td><td>0.9790</td><td>0.8589</td><td>0.4107</td><td><u>0.9814</u></td><td><b>0.9883</b></td></tr><tr><td>Case 2</td><td>0.2028</td><td>0.3745</td><td>0.8424</td><td>0.9612</td><td>0.9695</td><td>0.9710</td><td>0.4067</td><td>0.5060</td><td>0.9800</td><td>0.9057</td><td>0.4480</td><td><u>0.9833</u></td><td><b>0.9897</b></td></tr><tr><td>Case 3</td><td>0.0449</td><td>0.3837</td><td>0.8616</td><td>0.9286</td><td>0.9447</td><td>0.9210</td><td>0.3538</td><td>0.2951</td><td>0.9612</td><td>0.8497</td><td>0.3803</td><td><u>0.9714</u></td><td><b>0.9776</b></td></tr><tr><td>Case 4</td><td>0.1916</td><td>0.3781</td><td>0.8413</td><td>0.9464</td><td>0.9620</td><td>0.9605</td><td>0.3979</td><td>0.5009</td><td>0.9750</td><td>0.9003</td><td>0.4221</td><td><u>0.9759</u></td><td><b>0.9821</b></td></tr><tr><td>Case 5</td><td>0.0300</td><td>0.3222</td><td>0.5845</td><td>0.7879</td><td>0.5568</td><td>0.6988</td><td>0.3309</td><td>0.2413</td><td><u>0.8063</u></td><td>0.3408</td><td>0.3204</td><td>0.7061</td><td><b>0.8574</b></td></tr><tr><td>Case 6</td><td>0.0577</td><td>0.3804</td><td>0.5954</td><td>0.7949</td><td>0.6056</td><td>0.7328</td><td>0.3592</td><td>0.3588</td><td><u>0.8354</u></td><td>0.3477</td><td>0.3421</td><td>0.7347</td><td><b>0.8700</b></td></tr></tbody></table>

vertical and horizontal deadlines are randomly added to four bands in the others cases.

The denoised results are compared from two aspects. In terms of visual quality, the restoration results for Pavia City Center data in case 6 and for Washington DC Mall data in case 4 with heavier noise are displayed in Figs. 5 and

6, respectively. As observed, the proposed methods perform better than the other compared methods and are capable of better preserving the image structure information. Although Gaussian noise is removed by the BM4D method, important spatial information is lost and sparse noise still exists in Figs. 5(c) and 6 (c). The results obtained by LRMR, LRTA, and
<figure><img src="image_8.png" /><figcaption>Fig. 6. Denoised results by all the compared methods. (a) Original band 20 of Washington DC Mall data. (b) Simulated noise band of Case 4. (c) BM4D. (d) LRMR. (e) LRTV.(f) SSTV. (g) LSSTV. (h) LRTA. (i) PARAFAC. (j)LRTR. (k) LRTDTV. (l) LLRT. (m) WSWNN-SSTV. (n) WSWTNN-SSTV.</figcaption></figure>

<table><caption>TABLE IV
ERGAS RESULTS OF HSI IMAGE DENOISED BY DIFFERENT METHODS</caption><thead><tr><th rowspan="2">Data set</th><th rowspan="2">Noise case</th><th colspan="12">ERGAS</th></tr><tr><th>Noise</th><th>BM4D</th><th>LRMR</th><th>LRTV</th><th>SSTV</th><th>LSSTV</th><th>LRTA</th><th>PARAFAC</th><th>LRTDTV</th><th>LRTR</th><th>LLRT</th><th>WSWNN-SSTV</th><th>WSWTNN-SSTV</th></tr></thead><tbody><tr><td rowspan="6">Pavia</td><td>Case 1</td><td>68.0015</td><td>17.3808</td><td>3.8193</td><td>3.5958</td><td>2.7593</td><td>2.5688</td><td>21.3255</td><td>22.9294</td><td>2.6853</td><td>6.3983</td><td>16.3244</td><td><strong>2.4456</strong></td><td><strong>2.4646</strong></td></tr><tr><td>Case 2</td><td>50.7069</td><td>16.5926</td><td>8.2915</td><td>3.7784</td><td>3.7837</td><td>2.3592</td><td>19.0469</td><td>16.3088</td><td>2.3230</td><td>5.0862</td><td>15.5835</td><td><strong>2.1693</strong></td><td><strong>2.1941</strong></td></tr><tr><td>Case 3</td><td>68.1115</td><td>17.2600</td><td>6.9121</td><td>5.9264</td><td>3.616</td><td>3.4060</td><td>21.8697</td><td>23.4157</td><td>5.6618</td><td>8.8285</td><td>16.1223</td><td><strong>2.6514</strong></td><td><strong>3.1912</strong></td></tr><tr><td>Case 4</td><td>51.0999</td><td>16.0824</td><td>6.0535</td><td>3.6810</td><td>2.977</td><td>2.8466</td><td>18.9157</td><td>16.5251</td><td>2.8309</td><td>6.5567</td><td>13.2526</td><td><strong>2.6801</strong></td><td><strong>2.7025</strong></td></tr><tr><td>Case 5</td><td>67.0918</td><td>16.8994</td><td>8.0590</td><td>5.0908</td><td>2.1953</td><td>2.0820</td><td>21.3482</td><td>6.4678</td><td>6.4678</td><td>2.4283</td><td>12.1074</td><td><strong>1.7457</strong></td><td><strong>1.6349</strong></td></tr><tr><td>Case 6</td><td>59.2396</td><td>15.1093</td><td>12.2334</td><td>7.5667</td><td>8.4823</td><td>7.5448</td><td>20.5060</td><td>20.5667</td><td>8.5630</td><td>14.1691</td><td>13.3445</td><td><strong>6.7959</strong></td><td><strong>6.1625</strong></td></tr><tr><td rowspan="6">WDC</td><td>Case 1</td><td>355.074</td><td>84.8500</td><td>7.7708</td><td>3.5463</td><td>6.2720</td><td>14.5107</td><td>113.843</td><td>123.291</td><td>2.9365</td><td>17.1452</td><td>102.0863</td><td><strong>2.9004</strong></td><td><strong>1.9805</strong></td></tr><tr><td>Case 2</td><td>243.747</td><td>185.633</td><td>28.1338</td><td>3.2507</td><td>4.7539</td><td>6.3406</td><td>76.7945</td><td>73.1653</td><td>6.2493</td><td>13.9843</td><td>62.1659</td><td><strong>2.6877</strong></td><td><strong>1.8168</strong></td></tr><tr><td>Case 3</td><td>355.772</td><td>84.7557</td><td>7.7842</td><td>4.7409</td><td>6.5240</td><td>16.2167</td><td>115.175</td><td>123.864</td><td>3.8942</td><td>16.2452</td><td>102.3613</td><td><strong>3.3615</strong></td><td><strong>2.6618</strong></td></tr><tr><td>Case 4</td><td>244.146</td><td>184.434</td><td>11.9877</td><td>3.8577</td><td>4.9710</td><td>4.7246</td><td>77.2755</td><td>73.4074</td><td>5.6274</td><td>13.8686</td><td>62.2557</td><td><strong>3.0940</strong></td><td><strong>2.4120</strong></td></tr><tr><td>Case 5</td><td>374.466</td><td>139.085</td><td>21.6976</td><td>10.2586</td><td>38.4029</td><td>23.9023</td><td>168.195</td><td>174.502</td><td><strong>11.9841</strong></td><td>74.0830</td><td>119.350</td><td>23.3147</td><td><strong>10.4755</strong></td></tr><tr><td>Case 6</td><td>273.712</td><td>107.268</td><td>20.9605</td><td>9.9922</td><td>33.3881</td><td>21.4427</td><td>122.052</td><td>124.234</td><td><strong>10.6099</strong></td><td>71.6210</td><td>116.255</td><td>21.2984</td><td><strong>10.2868</strong></td></tr></tbody></table>

PARAFAC methods remain an amount of noise in heavy Gaussian noise case. LRTV and LLRT methods remove most noise, but LRTV loses structure spatial information due to the oversmoothing of TV regularization and the results by LLRT contain unexpectedly spectral distortion. SSTV, LSSTV, and LRTDTV methods perform comparatively better, while their

results contain some of Gaussian noise. Our two methods have similar performances for the first four cases, but WSWTNN-SSTV can better get rid of the complex mixed noise under cases 5 and 6. In particular, the proposed WSWTNN-SSTV method also preserves most of the image details such as sharp image edges and texture information, which demonstrates that
<table><caption>TABLE V
MSAD RESULTS OF HSI IMAGE DENOISED BY DIFFERENT METHODS</caption><thead><tr><th rowspan="2">Data set</th><th rowspan="2">Noise case</th><th colspan="10">MSAD</th></tr><tr><th>Noise</th><th>BM4D</th><th>LRMR</th><th>LRTV</th><th>SSTV</th><th>LSSTV</th><th>LRTA</th><th>PARAFAC</th><th>LRTDTV</th><th>LRTR</th><th>LLRT</th><th>WSWNN-SSTV</th><th>WSWTNN-SSTV</th></tr></thead><tbody><tr><td rowspan="6">Pavia</td><td>Case 1</td><td>48.1487</td><td>15.7919</td><td>4.3515</td><td>3.7876</td><td>3.4403</td><td>3.1755</td><td>17.6794</td><td>18.4703</td><td>3.3633</td><td>7.8053</td><td>14.1644</td><td><strong>3.0364</strong></td><td><strong>3.0858</strong></td></tr><tr><td>Case 2</td><td>40.7009</td><td>16.7266</td><td>4.6830</td><td>3.4024</td><td>3.4201</td><td>2.9265</td><td>18.2259</td><td>14.3863</td><td>2.9158</td><td>5.9385</td><td>15.8364</td><td><strong>2.6852</strong></td><td><strong>2.7453</strong></td></tr><tr><td>Case 3</td><td>48.4834</td><td>16.0268</td><td>6.8753</td><td>5.6626</td><td>4.4344</td><td>4.1709</td><td>18.5498</td><td>16.0390</td><td>5.5187</td><td>10.9026</td><td>16.0021</td><td><strong>3.0100</strong></td><td><strong>3.1912</strong></td></tr><tr><td>Case 4</td><td>41.1833</td><td>16.1055</td><td>4.8646</td><td>3.8939</td><td>3.6842</td><td>3.5407</td><td>18.2546</td><td>14.8096</td><td>3.5326</td><td>7.9818</td><td>14.8535</td><td><strong>3.3226</strong></td><td><strong>3.3874</strong></td></tr><tr><td>Case 5</td><td>48.0419</td><td>15.6301</td><td>7.0040</td><td>3.9529</td><td>2.6094</td><td>2.5024</td><td>17.8114</td><td>18.6730</td><td>3.0280</td><td>5.8036</td><td>12.9438</td><td>2.1150</td><td><strong>2.0039</strong></td></tr><tr><td>Case 6</td><td>45.5189</td><td>15.3750</td><td>14.4909</td><td>8.1627</td><td>10.5755</td><td>9.0383</td><td>15.4537</td><td>17.8667</td><td>8.9017</td><td>17.1997</td><td>13.3898</td><td><strong>8.5456</strong></td><td><strong>7.9373</strong></td></tr><tr><td rowspan="6">WDC</td><td>Case 1</td><td>60.2336</td><td>24.6210</td><td>15.2407</td><td>7.5015</td><td>14.4071</td><td>18.6880</td><td>24.8818</td><td>26.9924</td><td>7.5870</td><td>18.2557</td><td>24.8532</td><td><strong>6.0128</strong></td><td><strong>4.3518</strong></td></tr><tr><td>Case 2</td><td>57.5608</td><td>43.5927</td><td>15.2278</td><td>6.9783</td><td>11.6810</td><td>9.7260</td><td>29.3104</td><td>23.7123</td><td>7.4912</td><td>16.9463</td><td>25.4127</td><td><strong>5.5935</strong></td><td><strong>3.9811</strong></td></tr><tr><td>Case 3</td><td>60.3342</td><td>24.6435</td><td>15.4149</td><td>10.1867</td><td>15.0217</td><td>17.8179</td><td>25.3337</td><td>27.2228</td><td>8.6638</td><td>18.5253</td><td>26.4487</td><td><strong>6.7609</strong></td><td><strong>5.3567</strong></td></tr><tr><td>Case 4</td><td>57.7820</td><td>43.0781</td><td>16.5679</td><td>7.9379</td><td>12.2217</td><td>8.5295</td><td>29.5569</td><td>23.8139</td><td>8.0230</td><td>17.3248</td><td>25.4018</td><td><strong>6.2661</strong></td><td><strong>4.9149</strong></td></tr><tr><td>Case 5</td><td>58.1542</td><td>25.1780</td><td>21.3647</td><td>18.9256</td><td>25.2838</td><td>22.3527</td><td>25.8341</td><td>27.2882</td><td><u>19.6061</u></td><td>29.6429</td><td>26.3824</td><td>22.3257</td><td><strong>18.6128</strong></td></tr><tr><td>Case 6</td><td>55.8042</td><td>25.8897</td><td>21.6142</td><td>18.4107</td><td>23.6010</td><td>21.4145</td><td>25.2073</td><td>24.9394</td><td><u>18.6209</u></td><td>29.9820</td><td>25.3526</td><td>21.5023</td><td><strong>18.1870</strong></td></tr></tbody></table>

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

<figure><img src="image_10.png" /><figcaption>Fig. 7. (a) PSNR and (b) SSIM values of each band in case 6 of Pavia data.</figcaption></figure>

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

<figure><img src="image_12.png" /><figcaption>Fig. 8. (a) PSNR and (b) SSIM values of each band in case 4 of Washington DC data.</figcaption></figure>

more structure information is retained by using the tensor form of MLR-SSTV model.

Tables II–V present the denoised results by all the compared methods, where the best ones are marked in bold and the second ones are underlined. It is easy to see that our proposed MLR-SSTV methods perform the best among all the methods in terms of MPSNR, MSSIM, ERGAS, and MSAD.

TV-based approaches including LRTV, SSTV, LSSTV, LRT-DTV, and MLR-SSTV get higher metrics scores than other pure LR or TLR ones. The degradation in denoised image quality with the heavier mixed noise case is intuitively observed. In a relative sense, however, two MLR-SSTV methods suffer slower decays. Especially, WSWTNN-SSTV has a superior performance over other methods in case 6 whose
<figure><img src="image_13.png" /><figcaption>Fig. 9. MPSNR values with respect to different values of parameter $\theta$. (a) WSWNN-SSTV. (b) WSWTNN-SSTV.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 11. MPSNR values with respect to different values of parameter $\mu$. (a) WSWNN-SSTV. (b) WSWTNN-SSTV.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Fig. 10. MPSNR values with respect to different values of parameter $\lambda$. (a) WSWNN-SSTV. (b) WSWTNN-SSTV.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Fig. 12. MPSNR values with respect to different values of parameter $\lambda_{tv}$. (a) WSWNN-SSTV. (b) WSWTNN-SSTV.</figcaption></figure>

mixed noise has different statistics in different parts of the image. This is attributed to the use of MLR constraints and spatial-spectral prior information. Figs. 7 and 8 give the PSNR and SSIM values of each band in case 6 of Pavia data and case 4 of Washington data. Clearly, the PSNR and SSIM values of BM4D shown in black are lower than ones of other LR-based methods. LRMR, LRTV, LRTDTV, and SSTV perform worse when they encounter heavy dead-line noise. Our proposed WSWNN-SSTV and WSWTNN-SSTV shown in yellow and red provide the best performances.

1) **Parameter Analysis:** In the proposed MLR-SSTV solvers, five parameters, i.e., $\alpha$, $\lambda$, $\nu$, $\lambda_{tv}$, and $\beta$, need to be carefully identified. All the experiments were done under case 2.

The weights $\alpha = [\alpha_1, \alpha_2, \alpha_3]$ are important parameters for MLR-SSTV and are chosen differently for different cases. Particularly, for HSIs, the correlation along the spectral mode should be much stronger than those along the other two modes. Herein, we denote $\alpha$ as $[1, 1, \theta]/(2 + \theta)$ and only $\theta$ needs to be predetermined. As Fig. 9 shows, the MPSNR changes with respect to $\theta$. When $\theta = 10000$, WSWNN-SSTV achieves the highest MPSNR values. In WSWTNN-SSTV, the best results are obtained at $\theta = 0.001$ for HSIs Pavia City Center data and $\theta = 0.0001$ Washington DC Mall data. When $\theta > 10000$ and $\theta < 0.001$, WSWNN-SSTV and WSWTNN-SSTV can produce appealing results, respectively.

The regularization parameter $\lambda$ controls the strength of the sparse noise $S$. As observed from Fig. 10, WSWNN-SSTV can produce appealing results when $\lambda$ lies in the range of $[0.03 \sim 0.05]$. In Fig. 10(b), MPSNR values show the best performances when $\lambda$ is set as 0.02 for Pavia data and 0.003 for

Washington data. With the larger image size, $\lambda$ needs to be smaller.

Parameter $\mu$ is used to adjust the Gaussian noise strength. As Fig. 11 shows, WSWNN-SSTV is sensitive to $\mu$ and achieves the highest values for two HSIs data when $\mu = 1$ and $\mu = 5$. The MPSNR value curves of WSWTNN-SSTV have less fluctuations and the best denoised results are obtained when $\mu$ are set as 0.2 and 1, respectively. Parameter $\mu$ can be set as 1 firstly and then be adjusted according to the real condition of HSI noise.

Parameter $\lambda_{tv}$ restricts the impact of SSTV regularization on denoising effectiveness. As presented in Fig. 12, when $\lambda_{tv}$ increases from 0.01 to 0.05 for two data sets, the MPSNR values show upward trends. After peaking at $\lambda_{tv} = 0.05$, the performances fall off. Therefore, $\lambda_{tv} = 0.05$ is used in our experiments. The MPSNR values of WSWTNN-SSTV remain almost stable and the best ones are obtained at $\lambda_{tv} = 0.018$ and $\lambda_{tv} = 0.004$. When $\lambda_{tv}$ lies in the range of $[0.004 \sim 0.02]$, WSWTNN-SSTV can produce an appealing result.

Penalty parameter $\beta$ is used in the ADMM iteration. In Fig. 13, the MRSNR values increase slowly and drop slightly. The best scores of WSWNN-SSTV are obtained when $\beta = 1.2$. WSWTNN-SSTV has the best performances on Pavia data and Washington data at $\beta = 0.008$ and $\beta = 0.02$, respectively.

2) **Convergence Analysis:** The convergence of MLR-SSTV is empirically analyzed on the testing HSIs Pavia City Center under case 6 and Washington DC Mall under Noise case 4. Fig. 14 displays, RalCha, the relative change of two denoised results, such as, $||X(iter + 1) - X(iter)||_F^2/||X(iter)||_F^2$ where
<figure><img src="image_17.png" /><figcaption>Fig. 13. MPSNR values with respect to different values of parameter $\beta$. (a) WSWNN-SSTV. (b) WSWTNN-SSTV.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Fig. 14. RelCha values with respect to the number of iterations. (a) WSWNN-SSTV. (b) WSWTNN-SSTV.</figcaption></figure>

$\mathcal{X}$(iter) is the ith denoised result. It is easy to see that the RelCha values rapidly tend to zero, which numerically proves the convergence of two methods.

3) *Effectiveness of Weighted Nuclear Norm Used in Multidirectional Framework:* We provide evidence for the importance of effectively accounting for the weighted nuclear norm used in a multidirectional framework via an illustrative example. To describe this problem intuitively, we use: 1) WNN; 2) WSNN; 3) WSWNN; 4) WSTNN; and 5) WSWTNN without any regularization to denoise the HSI image corrupted by the Gaussian noise (SNR = 10~20), the salt and pepper noise (5%), and four vertical and horizontal dead-lines randomly in four bands, respectively.

As shown in Fig. 15, the denoised results with these four methods are compared with the original clean image and we zoom on the region marked with the red frame. The WNN, WSNN, and WTNN models are able to reduce Gaussian noise and some salt and pepper noise, but horizontal and vertical dead-line noise are still more or less present. In other words, it is not suitable for mixed noise removal to just apply weighted nuclear norm or MLR models. WSTNN and WSNN have better performances for removing dead-lines than WNN. As we can see, all these kinds of noise are denoised in the fourth and sixth images recovered by WSWNN and WSWTNN. Consistent with the visual observations, the MPSNR and MSSIM scores of the denoised image by WSWNN and WSWTNN methods are 3.68 dB, 0.0273 and 2.5 dB, 0.0080 better than those by WSNN and WSTNN methods, respectively. Moreover, WSWTNN is superior to WSWNN, in the recovery of texture information. Therefore, integrating weighted nuclear norm into multidirectional framework indeed appears a very appropriate way for remove dead-line noise.

## B. Real-World HSI Data Experiments

In this section, two real-world HSI data sets are presented, such as, the Hyperspectral Digital Imagery Collection Experiment (HYDICE) Urban data set and the Airborne Visible/Infrared Imaging Spectrometer (AVIRIS) Indian Pines data set. The above-compared methods are still implemented.

1) *HYDICE Urban Data Set:* This size of the original data is 307 × 307 × 210. This data set is corrupted by Gaussian noise, salt and pepper noise, stripes, the atmosphere, and

water absorption in Figs. 16(a) and 17(a). In this experiment, all images are selected as the test data including the bands polluted by the atmosphere and water absorption.

In Fig. 16, the first row gives the noisy and denoised images of band 139. The second row presents the magnified part from the red frame of the first row. This band is heavily contaminated by mixed noise. The BM4D and LRTA produce relatively worse results due to lots of residual Gaussian noise and dead-lines. The LRMR, SSTV, PAFARAC, and LRTR obtain better results, but dead-line noise still exists in the denoised results. With the help of TV regularization, LRTV and LRTDTV get rid of more dead-line noise. Although LLRT removes all mixed noise, it loses lots of useful structure information. For the LSSTV and the proposed method WSWNN-SSTV, most of the noise is removed and the image details are preserved. In Fig. 16(m) recovered by WSWTNN-SSTV, heavy mixed noise is removed completely and the trees, house roofs, and streets can be also distinguished easily. Similar results can be observed in Fig. 17. In Table VI, the images before and after denoising are compared by a no-reference HSI quantitative evaluation based on quality-sensitive features extraction, which is introduced in [64]. The lower score represents higher quality. Therefore, our proposed method WSWTNN-SSTV has the best performances according to its lowest score.

2) *AVIRIS Indian Pines Data Set:* The size of this data set is 145 × 145 × 220. As Figs. 18(a) and 19(a) show, the HSI data are corrupted by mixed noise. Similar to the Urban data experiment, we use all images as the test data and the same algorithms to remove the noise.

Figs. 18 and 19 give the denoised results by different methods. As Figs. 18(g) and 19(g) show, some of the mixed noise is removed and the LRTA method obtains the worst results. PAFARAC produces artifacts and LRTR, but LLRT are not capable of handling different noise cases. For the BM4D, LRMR, and LRTV, some dead-line noise exists, as presented in Figs. 18(c) and (d) and 19(c) and (d). Although SSTV, LRTDTV, and LSSTV get rid of mixed noise better, they cause spatial over smoothing and lose some image details, as shown in the second row of Figs. 18 and 19. Our methods achieve the best performances, and WSWTNN-SSTV preserves most image details such as the clearest image edges, which is beneficial for the consequent HSI analysis work.
<figure><img src="image_19.png" /><figcaption>Fig. 15. Resulting HSIs with their (Left) MPSNR and (Right) MSSIM. (a) Ground truth. (b) Noisy HSI. (c) WNN. (d) WSNN. (e) WSWNN. (f) WSTNN. (g) WSWTNN.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Fig. 16. Denoised results by all the compared methods in the Urban data. (a) Original band 139. (b) BM4D. (c) LRMR. (d) LRTV. (e) SSTV. (f) LSSTV. (g) LRTA. (h) PAFARAC. (i) LRTDTV. (j) LRTR. (k) LLRT. (l) WSWNN-SSTV. (m) WSWTNN-SSTV.</figcaption></figure>

<table><caption>TABLE VI
BLIND QUANTITATIVE EVALUATION ON THE HYDICE URBAN DATA BY DIFFERENT METHODS</caption><thead><tr><th>Methods</th><th>BM4D</th><th>LRMR</th><th>LRTV</th><th>SSTV</th><th>LSSTV</th><th>LRTA</th><th>PAFARAC</th><th>LRTDTV</th><th>LRTR</th><th>LLRT</th><th>WSWNN-SSTV</th><th>WSWTNN-SSTV</th></tr></thead><tbody><tr><th>Metric Values</th><td>12.4500</td><td>12.3946</td><td>12.4789</td><td>12.4073</td><td>12.0369</td><td>12.4435</td><td>13.6753</td><td>12.2077</td><td>12.1133</td><td>14.7524</td><td><u>11.8505</u></td><td><b>11.4324</b></td></tr></tbody></table>

To further compare the performances of these methods, we apply the well-known support vector machine (SVM) [65] algorithm into the supervised classification experiment. As the first column of Table VII shows, there are 16 different classes

employed to test the classification accuracy. The random 10% of samples from each class generate the training sets and the others are used to be the test sets. The classification results before and after denoising can be observed in Fig. 20. The
<figure><img src="image_21.png" /><figcaption>Fig. 17. Denoised results by all the compared methods in the Urban data. (a) Original band 207. (b) BM4D. (c) LRMR. (d) LRTV. (e) SSTV. (f) LSSTV. (g) LRTA. (h) PAFARAC. (i) LRTDTV. (j) LRTR. (k) LLRT. (l) WSWNN-SSTV. (m) WSWTNN-SSTV.</figcaption></figure>

<table><caption>TABLE VII
CLASSIFICATION RESULTS ON THE AVIRIS INDIAN PINES DATA BY DIFFERENT METHODS</caption><thead><tr><th>Class</th><th>Train/Test</th><th>Origin</th><th>BM4D</th><th>LRMR</th><th>LRTV</th><th>SSTV</th><th>LSSTV</th><th>LRTA</th><th>PAFARAC</th><th>LRTDTV</th><th>LRTR</th><th>LLRT</th><th>WSWNN-SSTV</th><th>WSWTNN-SSTV</th></tr></thead><tbody><tr><td>Alfalfa</td><td>5/41</td><td>0.6341</td><td>0.7560</td><td>0.7317</td><td>0.6829</td><td>0.7804</td><td><u>0.8048</u></td><td>0.6341</td><td>0.6585</td><td>0.6829</td><td>0.6585</td><td>0.6341</td><td><b>0.9268</b></td><td>0.7808</td></tr><tr><td>Corn-notill</td><td>143/1285</td><td>0.8202</td><td>0.8708</td><td>0.9105</td><td>0.8287</td><td>0.9167</td><td>0.9229</td><td>0.8249</td><td>0.7992</td><td>0.8708</td><td>0.7642</td><td>0.8148</td><td><u>0.9361</u></td><td><b>0.9385</b></td></tr><tr><td>Corn-mintill</td><td>83/747</td><td>0.7161</td><td>0.8728</td><td>0.8955</td><td>0.8888</td><td>0.9518</td><td>0.9598</td><td>0.7175</td><td>0.6854</td><td>0.8969</td><td>0.8032</td><td>0.7604</td><td><u>0.9732</u></td><td><b>0.9772</b></td></tr><tr><td>Corn</td><td>24/213</td><td>0.8028</td><td>0.9014</td><td>0.9248</td><td>0.8309</td><td>0.9154</td><td>0.8873</td><td>0.8122</td><td>0.8028</td><td><u>0.9201</u></td><td>0.8403</td><td>0.9061</td><td>0.9061</td><td><b>0.9436</b></td></tr><tr><td>Grass-pasture</td><td>49/434</td><td>0.9377</td><td>0.9608</td><td>0.9354</td><td>0.9700</td><td>0.9723</td><td>0.9700</td><td>0.9378</td><td>0.9493</td><td>0.9746</td><td>0.9770</td><td>0.9401</td><td><u>0.9815</u></td><td><b>0.9862</b></td></tr><tr><td>Grass-trees</td><td>73/657</td><td>0.9497</td><td>0.9847</td><td>0.9619</td><td>0.9786</td><td>0.9863</td><td>0.9847</td><td>0.9482</td><td>0.9589</td><td>0.9878</td><td>0.9863</td><td>0.9513</td><td><b>0.9939</b></td><td><u>0.9924</u></td></tr><tr><td>Grass-pasture-mowed</td><td>3/25</td><td>0.6400</td><td>0.7200</td><td>0.6400</td><td>0.7600</td><td><u>0.8400</u></td><td><u>0.8400</u></td><td>0.6400</td><td>0.6800</td><td>0.7200</td><td>0.7200</td><td>0.6400</td><td><b>0.9200</b></td><td><u>0.8400</u></td></tr><tr><td>Hay-windrowed</td><td>48/430</td><td>0.9906</td><td>0.9930</td><td>0.9976</td><td><b>1.0000</b></td><td><u>0.9953</u></td><td><u>0.9953</u></td><td>0.9907</td><td><b>1.0000</b></td><td>0.9930</td><td><b>1.0000</b></td><td>0.9953</td><td><b>1.0000</b></td><td><b>1.0000</b></td></tr><tr><td>Oats</td><td>2/18</td><td>0.7222</td><td>0.8333</td><td>0.8333</td><td><b>1.0000</b></td><td><u>0.9444</u></td><td><u>0.9444</u></td><td>0.7778</td><td>0.4444</td><td><b>1.0000</b></td><td>0.7222</td><td>0.7222</td><td><b>1.0000</b></td><td><b>1.0000</b></td></tr><tr><td>Soybean-notill</td><td>98/874</td><td>0.7059</td><td>0.8478</td><td>0.8821</td><td>0.7276</td><td><u>0.9405</u></td><td><b>0.9530</b></td><td>0.7037</td><td>0.7700</td><td>0.8329</td><td>0.7483</td><td>0.7334</td><td>0.8981</td><td>0.9302</td></tr><tr><td>Soybean-mintill</td><td>246/2209</td><td>0.8415</td><td>0.9198</td><td>0.9388</td><td>0.8632</td><td>0.9601</td><td><u>0.9655</u></td><td>0.8483</td><td>0.8891</td><td>0.9157</td><td>0.8420</td><td>0.8520</td><td>0.9479</td><td><b>0.9720</b></td></tr><tr><td>Soybean-clean</td><td>60/533</td><td>0.8480</td><td>0.9530</td><td>0.9324</td><td>0.9643</td><td><u>0.9699</u></td><td>0.9756</td><td>0.7880</td><td>0.9062</td><td>0.9437</td><td>0.8480</td><td>0.8668</td><td>0.9662</td><td><b>0.9737</b></td></tr><tr><td>Wheat</td><td>21/184</td><td>0.9945</td><td>0.9945</td><td>0.9945</td><td>0.9945</td><td><u>0.9945</u></td><td><b>0.9946</b></td><td>0.9945</td><td>0.9945</td><td>0.9891</td><td><u>0.9945</u></td><td>0.9945</td><td><b>0.9946</b></td><td><b>0.9946</b></td></tr><tr><td>Woods</td><td>127/1138</td><td>0.9516</td><td>0.9666</td><td>0.9455</td><td>0.9472</td><td>0.9674</td><td>0.9648</td><td>0.9468</td><td>0.9745</td><td><b>0.9964</b></td><td>0.9517</td><td>0.9622</td><td><u>0.9912</u></td><td>0.9886</td></tr><tr><td>Buildings-Grass-Trees-Drives</td><td>39/347</td><td>0.6311</td><td>0.7608</td><td>0.7809</td><td>0.8040</td><td>0.8587</td><td>0.8443</td><td>0.8040</td><td>0.8991</td><td>0.8703</td><td>0.6311</td><td>0.6571</td><td><b>0.9654</b></td><td><u>0.9135</u></td></tr><tr><td>Stone-Steel-Towers</td><td>10/83</td><td>0.9518</td><td>0.9518</td><td>0.9518</td><td><b>1.0000</b></td><td>0.9638</td><td>0.9638</td><td>0.9759</td><td>0.9277</td><td><u>0.9879</u></td><td>0.9518</td><td>0.9518</td><td><b>1.0000</b></td><td>0.9638</td></tr><tr><td>kappa</td><td>-</td><td>0.8197</td><td>0.9019</td><td>0.9131</td><td>0.8675</td><td>0.9451</td><td>0.9482</td><td>0.8273</td><td>0.8618</td><td>0.9114</td><td>0.8618</td><td>0.8360</td><td><u>0.9526</u></td><td><b>0.9603</b></td></tr><tr><td>OA</td><td>-</td><td>0.8423</td><td>0.9140</td><td>0.9237</td><td>0.8835</td><td>0.9518</td><td>0.9545</td><td>0.8484</td><td>0.8790</td><td>0.9223</td><td>0.8790</td><td>0.8564</td><td><u>0.9583</u></td><td><b>0.9652</b></td></tr><tr><td>AA</td><td>-</td><td>0.8212</td><td>0.893</td><td>0.8911</td><td>0.8901</td><td>0.9349</td><td>0.9357</td><td>0.8191</td><td>0.8632</td><td>0.9114</td><td>0.8632</td><td>0.8364</td><td><u>0.9626</u></td><td><b>0.9688</b></td></tr></tbody></table>

accuracy for each class, kappa coefficient, overall accuracy (OA), and average accuracy (AA) are also given in Table VII. In Fig. 20(a), the classification results in the original noised case appear interleaved, and the kappa, OA, and AA are just

0.8197, 0.8423, and 0.8212. After denoising, the accuracy and three metrics are more or less improved. Consistence with the intuitive images, the metric values of the BM4D, LRTA, PAFARAC, LRTR, and LLRT rise less due to their incomplete
<figure><img src="image_22.png" /><figcaption>Fig. 18. Denoised results by all the compared methods in the Indian Pine data. (a) Original band 108. (b) BM4D. (c) LRMR. (d) LRTV. (e) SSTV. (f) LSSTV. (g) LRTA. (h) PAFARAC. (i) LRTDTV. (j) LRTR. (k) LLRT. (l) WSWNN-SSTV. (m) WSWTNN-SSTV.</figcaption></figure>

<table><caption>TABLE VIII
RUNNING TIME COMPARISON OF DIFFERENT METHODS</caption><thead><tr><th>HSI data sets</th><th>size</th><th>BM4D</th><th>LRMR</th><th>LRTV</th><th>SSTV</th><th>LSSTV</th><th>LRTA</th><th>PAFARAC</th><th>LRTDTV</th><th>LRTR</th><th>LLRT</th><th>WSWNN-SSTV</th><th>WSWTNN-SSTV</th></tr></thead><tbody><tr><td>Pavia Data</td><td>128 × 128 × 80</td><td>117.0s</td><td>9.8s</td><td>15.8s</td><td>95.9s</td><td>633.4s</td><td>2.8s</td><td>163.5s</td><td>42.1s</td><td>16.1s</td><td>281.9s</td><td>427.2s</td><td>495.9s</td></tr><tr><td>WashingtonDC data</td><td>256 × 256 × 191</td><td>1305.4s</td><td>142.9s</td><td>232.5s</td><td>1413.2s</td><td>5258.5s</td><td>16.5s</td><td>5687.5s</td><td>655.2s</td><td>163.0s</td><td>1954.8s</td><td>1250.5s</td><td>2065.1s</td></tr><tr><td>Urban</td><td>307 × 307 × 210</td><td>1817.1s</td><td>179.1s</td><td>317.3s</td><td>1515.2s</td><td>8654.4s</td><td>10.2s</td><td>16654.7s</td><td>646.5s</td><td>419.5s</td><td>3453.5s</td><td>4256.1s</td><td>5056.9s</td></tr><tr><td>Indian Pines</td><td>145 × 145 × 220</td><td>431.5s</td><td>115.1s</td><td>64.5s</td><td>466.5s</td><td>2559.5s</td><td>2.7s</td><td>2036.6s</td><td>52.7s</td><td>72.7s</td><td>946.9s</td><td>1476.9s</td><td>2220.6s</td></tr></tbody></table>

noise removal. LRMR, LRTV, and LRTDTV also cause similar evaluation results since they all exploit the global LR property of HSI. SSTV, LSSTV, and the proposed MLR-SSTV models produce better results, whereas WSWTNN-SSTV gets the highest metric values with kappa, OA, and AA of 0.9603, 0.9652, and 0.9688.

## C. Discussion

The computational complexity of the proposed algorithm mainly consists of the following parts.

1) The computation of solving $M_k$ depends on the SVD. For WSWNN-SSTV and WSWTNN-SSTV, the complexities of these two operators are $O(3h^2v^2z+3z^3)$ and $O(hvz(\log(hvz) + \min(h, v) + \min(v, z) + \min(z, h)))$, respectively.

2) Solving the subproblem (18) includes three threshold and two SSTV operators, which takes about $O(3hvz + 2hvz)$ iterations.

3) Updating $N$ requires about $O(hvz)$ iterations.

4) LSQR algorithm is used to update $x$ and its computation costs $O(3hv + 5z)$.

5) Most of the computation time for updating the multipliers lies on the SSTV operators, which requires about $O(2hvz)$ iterations.

Therefore, the computational complexity of the whole process are in total $O(3h^2v^2z+3z^3+8hvz+3hv+5z)$ for WSWNN-SSTV and $O(hvz(\log(hvz) + \min(h, v) + \min(v, z) + \min(z, h))) + 8hvz + 3hv + 5z)$ for WSWTNN-SSTV.

The running times of different methods are compared in Table VIII. All of the experiments are conducted in MATLAB 2018b on the same PC with an Intel i7 CPU at 2.60 GHz
<figure><img src="image_23.png" /><figcaption>Fig. 19. Denoised results by all the compared methods in the Indian Pine data. (a) Original band 220. (b) BM4D. (c) LRMR. (d) LRTV. (e) SSTV. (f) LSSTV. (g) LRTA. (h) PAFARAC. (i) LRTDTV. (j) LRTR. (k) LLRT. (l) WSWNN-SSTV. (m) WSWTNN-SSTV.</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>Fig. 20. Classification results by all the compared methods in the Indian Pine data. (a) Original. (b) BM4D. (c) LRMR. (d) LRTV. (e) SSTV. (f) LSSTV. (g) LRTA. (h) PAFARAC. (i) LRTDTV. (j) LRTR. (k) LLRT. (l) WSWNN-SSTV. (m) WSWTNN-SSTV.</figcaption></figure>

and 32 GB of memory. The LRTA method is the fastest one among the state-of-the-art methods, but it can hardly remove all the noise. LRMR and LRTR are the second and third faster, but they retain lots of dead-lines noise. LRTV and LRTDTV

also belong to the fasters. LRTV just uses the band-by-band TV while LRTDTV needs a number of **SVD** operators. The running times of LSSTV, WSWNN-SSTV and WSWTNN-SSTV are more than the SSTV due to the combinations
of the LR constraint. But WSWNN-SSTV and WSWTNN-SSTV converge more rapidly, utilize fewer iterations and cost less time than LSSTV since our methods consider the more comprehensive MLR framework.

## V. CONCLUSION

In this article, we have proposed a novel HSI mixed noise removal method MLR-SSTV by combining the MLR-SSTV. The MLR relationship of HSI data is explored from 2-D matrices and 3-D tensors. A weighted nuclear norm is integrated into this multidirectional framework to remove dead-line noise completely. SSTV regularization is employed to get rid of heavy Gaussian noise and enforce the smoothness in the spatial and spectral domains. Especially, WSWTNN-SSTV can preserve more image details such as sharp edges and texture information. The optimization problem is solved by the ADMM framework efficiently. In both simulated and real-world data experiments, the proposed MLR-SSTV provides the best performances among all the compared state-of-the-art methods. Meanwhile, the effects of HSI denoising as a pre-processing step for HSI classification have been verified and classification accuracy is improved most when the proposed MLR-SSTV methods are used.

For the future work, two extensions could be developed to improve our method: 1) we will investigate more efficient multidirectional model, such as more modifications for the LR tensor model and 2) the $l_1$ norm of the SSTV regularization has been widely applied for the HSI denoising. Another available norm, $l_0$ norm, has been applied in our other work [66] and it will be developed into the multidirectional model.

## REFERENCES

[1] J. M. Bioucas-Dias, A. Plaza, G. Camps-Valls, P. Scheunders, N. Nasrabadi, and J. Chanussot, "Hyperspectral remote sensing data analysis and future challenges," *IEEE Geosci. Remote Sens. Mag.*, vol. 1, no. 2, pp. 6-36, 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, 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.*, early access, Apr. 29, 2020, doi: 10.1109/NGRS.2020.2979764.

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

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

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

[7] J. M. Bioucas-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] 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.

[10] Z. He, L. Liu, R. Deng, and Y. Shen, "Low-rank group inspired dictionary learning for hyperspectral image classification," *Signal Process.*, vol. 120, pp. 209-221, Mar. 2016.

[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] B. Du, R. Zhao, L. Zhang, and L. Zhang, "A spectral-spatial based local summation anomaly detection method for hyperspectral images," *Signal Process.*, vol. 124, pp. 115-131, Jul. 2016.

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

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

[17] H. Othman and S.-E. Qian, "Noise reduction of hyperspectral imagery using hybrid spatial-spectral derivative-domain wavelet shrinkage," *IEEE Trans. Geosci. Remote Sens.*, vol. 44, no. 2, pp. 397-408, Feb. 2006.

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

[19] J. Li, Q. Yuan, H. Shen, and L. Zhang, "Hyperspectral image recovery employing a multidimensional nonlocal total variation model," *Signal Process.*, vol. 111, pp. 230-248, Jun. 2015.

[20] N. Renard, S. Bourennane, and J. Blanc-Talon, "Denoising and dimensionality reduction using multilinear tools for hyperspectral images," *IEEE Geosci. Remote Sens. Lett.*, vol. 5, no. 2, pp. 138-142, Apr. 2008.

[21] B. Priego, R. J. Duro, and J. Chanussot, "4DCAF: A temporal approach for denoising hyperspectral image sequences," *Pattern Recognit.*, vol. 72, pp. 433-445, Dec. 2017.

[22] A. A. Green, M. Berman, P. Switzer, and M. D. Craig, "A transformation for ordering multispectral data in terms of image quality with implications for noise removal," *IEEE Trans. Geosci. Remote Sens.*, vol. GRS-26, no. 1, pp. 65-74, Jan. 1988.

[23] M. Elad and M. Aharon, "Image denoising via sparse and redundant representations over learned dictionaries," *IEEE Trans. Image Process.*, vol. 15, no. 12, pp. 3736-3745, Dec. 2006.

[24] K. Dabov, A. Foi, V. Katkovnik, and K. Egiazarian, "Image denoising by sparse 3-D transform-domain collaborative filtering," *IEEE Trans. Image Process.*, vol. 16, no. 8, pp. 2080-2095, Aug. 2007.

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

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

[27] Y.-Q. Zhao and J. Yang, "Hyperspectral image denoising via sparse representation and low-rank constraint," *IEEE Trans. Geosci. Remote Sens.*, vol. 53, no. 1, pp. 296-308, Jan. 2015.

[28] J. M. Bioucas-Dias and J. M. P. Nascimento, "Hyperspectral subspace identification," *IEEE Trans. Geosci. Remote Sens.*, vol. 46, no. 8, pp. 2435-2445, Aug. 2008.

[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] Y. Xie, S. Gu, Y. Liu, W. Zuo, W. Zhang, and L. Zhang, "Weighted schatten p-norm minimization for image denoising and background subtraction," *IEEE Trans. Image Process.*, vol. 25, no. 10, pp. 4842-4857, Oct. 2016.

[31] Y. Chen, Y. Guo, Y. Wang, D. Wang, C. Peng, and G. He, "Denoising of hyperspectral images using nonconvex low rank matrix approximation," *IEEE Trans. Geosci. Remote Sens.*, vol. 55, no. 9, pp. 5366-5380, Sep. 2017.

[32] F. Xu, Y. Chen, C. Peng, Y. Wang, X. Liu, and G. He, "Denoising of hyperspectral image using low-rank matrix factorization," *IEEE Geosci. Remote Sens. Lett.*, vol. 14, no. 7, pp. 1141-1145, Jul. 2017.
[33] L. Gao, D. Yao, Q. Li, L. Zhuang, B. Zhang, and J. Bioucas-Dias, "A new low-rank representation based hyperspectral image denoising method for mineral mapping," *Remote Sens.*, vol. 9, no. 11, p. 1145, 2017.

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

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

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

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

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

[39] Q. Xie et al., "Multispectral images denoising by intrinsic tensor sparsity regularization," in *Proc. IEEE Conf. Comput. Vis. Pattern Recognit. (CVPR)*, Jun. 2016, pp. 1692-1700.

[40] X. Liu, S. Bourennane, and C. Fossati, "Denoising of hyperspectral images using the PARAFAC model and statistical performance analysis," *IEEE Trans. Geosci. Remote Sens.*, vol. 50, no. 10, pp. 3717-3724, Oct. 2012.

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

[42] Z. Zhang, G. Ely, S. Aeron, N. Hao, and M. Kilmer, "Novel methods for multilinear data completion and de-noising based on tensor-SVD," in *Proc. IEEE Conf. Comput. Vis. Pattern Recognit.*, Jun. 2014, vol. 44, no. 9, pp. 3842-3849.

[43] C. Lu, J. Feng, Y. Chen, W. Liu, Z. Lin, and S. Yan, "Tensor robust principal component analysis: Exact recovery of corrupted low-rank tensors via convex optimization," in *Proc. IEEE Conf. Comput. Vis. Pattern Recognit. (CVPR)*, Jun. 2016, pp. 5249-5257.

[44] H. Fan, Y. Chen, Y. Guo, H. Zhang, and G. Kuang, "Hyperspectral image restoration using low-rank tensor recovery," *IEEE J. Sel. Topics Appl. Earth Observ. Remote Sens.*, vol. 10, no. 10, pp. 4589-4604, Oct. 2017.

[45] Y. Chen, S. Wang, and Y. Zhou, "Tensor nuclear norm-based low-rank approximation with total variation regularization," *IEEE J. Sel. Topics Signal Process.*, vol. 12, no. 6, pp. 1364-1377, Dec. 2018.

[46] Y. Chang, L. Yan, and S. Zhong, "Hyper-laplacian regularized unidirectional low-rank tensor recovery for multispectral image denoising," in *Proc. IEEE Conf. Comput. Vis. Pattern Recognit. (CVPR)*, Jul. 2017, pp. 4260-4268.

[47] W. He, Q. Yao, C. Li, N. Yokoya, and Q. Zhao, "Non-local meets global: An integrated paradigm for hyperspectral denoising," in *Proc. IEEE/CVF Conf. Comput. Vis. Pattern Recognit. (CVPR)*, Jun. 2019, pp. 6868-6877.

[48] H. Fan, C. Li, Y. Guo, G. Kuang, and J. Ma, "Spatial-spectral total variation regularized low-rank tensor decomposition for hyperspectral image denoising," *IEEE Trans. Geosci. Remote Sens.*, vol. 56, no. 10, pp. 6196-6213, Oct. 2018.

[49] 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

[50] S. Gandy, B. Recht, and I. Yamada, "Tensor completion and low-n-rank tensor recovery via convex optimization," *Inverse Problems*, vol. 27, no. 2, Feb. 2011, Art. no. 025010.

[51] Z. Zhang, D. Liu, S. Aeron, and A. Vetro, "An online tensor robust PCA algorithm for sequential 2D data," in *Proc. IEEE Int. Conf. Acoust., Speech Signal Process. (ICASSP)*, Mar. 2016, pp. 2434-2438.

[52] T. G. Kolda and B. W. Bader, "Tensor decompositions and applications," *SIAM Rev.*, vol. 51, no. 3, pp. 455-500, Aug. 2009.

[53] R. Zhao, Q. Wang, J. Fu, and L. Ren, "Exploiting block-sparsity for hyperspectral kronecker compressive sensing: A tensor-based Bayesian method," *IEEE Trans. Image Process.*, vol. 29, pp. 1654-1668, Jul. 2020. [Online]. Available: https://ieeexplore.ieee.org/document/8861298

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

[55] S. Gu, L. Zhang, W. Zuo, and X. Feng, "Weighted nuclear norm minimization with application to image denoising," in *Proc. IEEE Conf. Comput. Vis. Pattern Recognit.*, Jun. 2014, pp. 2862-2869.

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

[57] S. Boyd, "Distributed optimization and statistical learning via the alternating direction method of multipliers," *Found. Trends Mach. Learn.*, vol. 3, no. 1, pp. 1-122, 2010.

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

[59] Z. Lin, M. Chen, and Y. Ma, "The augmented Lagrange multiplier method for exact recovery of corrupted low-rank matrices," 2010, arXiv:1009.5055. [Online]. Available: http://arxiv.org/abs/1009.5055

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

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

[62] M. Maggioni, V. Katkovnik, K. Egiazarian, and A. Foi, "Nonlocal transform-domain filter for volumetric data denoising and reconstruction," *IEEE Trans. Image Process.*, vol. 22, no. 1, pp. 119-133, Jan. 2013.

[63] J. Li, Q. Yuan, H. Shen, and L. Zhang, "Noise removal from hyperspectral image with joint spectral-spatial distributed sparse representation," *IEEE Trans. Geosci. Remote Sens.*, vol. 54, no. 9, pp. 5425-5439, Sep. 2016.

[64] J. Yang, Y. Zhao, C. Yi, and J. C.-W. Chan, "No-reference hyperspectral image quality assessment via quality-sensitive features learning," *Remote Sens.*, vol. 9, no. 4, p. 305, 2017.

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

[66] M. Wang, Q. Wang, and J. L. B. F. J. de Vries, "L0 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.

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

Minghua Wang 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 Univ. Grenoble Alpes, CNRS, Grenoble INP, GIPSA-lab, Grenoble, France, from 2019 to 2020. Her research interests include hyperspectral image denoising, tensor completion, and low-rank representation.

<figure><img src="image_26.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, multi-sensor data fusion, wireless sensor networks, and intelligent detection technology.
<figure><img src="image_27.png" /></figure>

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

Since 1999, he has been with Grenoble INP, where he is a Professor of signal and image processing. He has been a Visiting Scholar with Stanford University, Stanford, CA, USA, KTH, Stockholm, Sweden, and 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 an Adjunct professor with the Chinese Academy of Sciences, Aerospace Information research Institute, Beijing, China.

Dr. Chanussot was the Founding President of the IEEE Geoscience and Remote Sensing French Chapter from 2007 to 2010, which received the 2010 IEEE GRS-S Chapter Excellence Award. He was a member of the Machine Learning for Signal Processing Technical Committee of the IEEE Signal Processing Society from 2006 to 2008. He was the Vice-President of the IEEE Geoscience and Remote Sensing Society, in charge of meetings and symposia from 2017 to 2019. He has received multiple outstanding paper awards. 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 from 2009 to 2011 and the Co-Chair of the

GRS Data Fusion Technical Committee from 2005 to 2008. He was 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 from 2011 to 2015. In 2014, he served as a Guest Editor for the *IEEE Signal Processing Magazine*. He is a member of the Institut Universitaire de France from 2012 to 2017 and a Highly Cited Researcher (Clarivate Analytics/Thomson Reuters, 2018–2019).

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

**Dan Li** 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 2012, 2014, and 2018, respectively.

Since 2018, she has been a Lecturer with the College of Astronautics, Nanjing University of Aeronautics and Astronautics (NUAA), Nanjing, China. Her research interests include hyperspectral image classification, signal processing, sparse sampling, and reconstruction technology.