A Novel Tensor-Based Hyperspectral Image Restoration Method With Low-Rank Modeling in Gradient Domains

The hyperspectral image (HSI) is easily contaminated by various kinds of mixed noise (such as Gaussian noise, impulse noise, stripes, and deadlines) during the process of data acquisition and conversion, which significantly affect the quality and applications of HSI. As an important and effective scheme for the quality improvement of HSI, the HSI restoration problem aims to recover a clean HSI from the noisy HSI with mixed noise. Thus, based on the tensor modeling of HSI, we propose a novel tensor-based HSI restoration model with low-rank modeling in gradient domains in a unified tensor representation framework in this article. First, for the spectral low-rank modeling of HSI in spectral gradient domain, we particularly exploit the low-rank property of spectral gradient, and propose the spectral gradient-based weighted nuclear norm low-rank prior term. Second, for the spatial-mode low-rank modeling of HSI in spatial gradient domain, we particularly exploit the low-rank property of spatial gradient tensors via the discrete Fourier transform, and propose the spatial gradient-based tensor nuclear norm low-rank prior term. Then, we use the alternative direction method of multipliers to solve the proposed model. Finally, the restoration results on both the simulated and real HSI datasets demonstrate that the proposed method is superior to many state-of-the-art methods in the aspects of visual and quantitative comparisons.

, [6]. However, due to the influence of acquisition equipment and external environment, HSIs are usually contaminated by various kinds of mixed noise (such as, Gaussian noise, impulse noise, deadlines, and stripes), which not only affects the image quality, but also limits the later processing and applications of HSI. Therefore, HSI restoration, which is recently an important scheme to remove the mixed noise so as to improve the image quality of HSI, has very important research and application significance. Specifically, as an effective way to remove the noise, the core of image priors-based HSI restoration model is to make full use of the effective spectral and spatial priors of HSI, such as, nonlocal self-similarity, spatial-spectral smoothness, and low-rank priors [7], [8], [9].
As we know, the nonlocal self-similarity was earlier used for grayscale image denoising, which assumed that a grayscale image contained lots of similar 2-D blocks (i.e., patches) at different positions of the image. Then, many traditional methods using spatial nonlocal self-similarity were generalized to denoise HSI band by band, such as block matching and 3-D filtering (BM3D) [10], dictionary learning [11], and nonlocal means filter [12], [13]. However, these methods treated each band of HSI as a grayscale image, which did not consider the important correlation between spectral bands, so the restoration performance was not very ideal and could be further improved. To solve this problem, scholars investigated the spatial and spectral nonlocal self-similarity-based restoration methods using similar 3-D cubes to replace the common 2-D blocks, such as MSPCA-BM3D [14], BM4D [15], 3-D nonlocal sparse (3DNLS) [16], and nonlocal tensor dictionary learning [17].
More clearly, the methods based on low-rank priors have been extensively studied and achieved great success in the field of HSI restoration. Specifically, inspired by the robust principal component analysis (RPCA) model [18], which modeled the observed data as the sum of a low-rank matrix and a sparse matrix, Zhang et al. [19] proposed the low-rank matrix restoration (LRMR) model. Moreover, Chen et al. [20] proposed a nonconvex low-rank matrix approximation (NonLRMA) model, which used a nonconvex regularizer to replace the traditional nuclear norm, so the rank function can be better approximated. Although LRMR can effectively remove the mixed noise in HSI, it does not impose any spatial constraints on the neighboring pixels of HSI. To solve this problem, Zhu et al. [21] considered both spectral and spatial constraints to preserve the fine structure of This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ HSI, and proposed the spectral nonlocal low-rank model. Based on the LRMR model, He et al. [22] proposed the noise-adjusted iterative low-rank matrix approximation (NAILRMA) method, which considered the different noise in different bands and can well remove the noise from HSI. However, when the local noise is larger, the denoising performance will become worse. To address this shortcoming, Wang et al. [23] proposed the group low-rank representation model, where the local similarity within patches and the nonlocal similarity between groups of patches were both used to preserve the HSI image features and remove the noise. What is more, Fan et al. [24] proposed the superpixel segmentation and low-rank representation (SSLRR) model to remove different types of noise.
Moreover, since the total variation (TV) regularization method [25] can effectively preserve shape edges and enhance piecewise smoothness, then various TV-based spatial-spectral smoothness prior models are also proposed for HSI restoration, such as cubic TV (CTV) [26], spectral-spatial adaptive hyperspectral TV (SSAHTV) [27], spectrally adaptive multidimensional nonlocal total variation (SAMNLTV) [28], and spatio-spectral TV (SSTV) [29]. Although these methods fully use the local spatial and spectral information of HSI, they ignore the low-rank characteristics of HSI. Based on this, He et al. [30] proposed the TV-regularized low-rank matrix factorization (LRTV) model, which introduced the band-by-band TV regularization term and the low-rank prior together for HSI restoration. However, LRTV considered the TV model separately for each spectral band, ignoring the spectral smoothness of HSI. To overcome this drawback, Wang et al. [31] used the SSTV [29] to model the spatial-spectral smoothness of HSI, and proposed the low-rank constraint and spatial spectral total variation (LSSTV) model, which can effectively remove the mixed noise and preserve the edge structure. Moreover, He et al. [32] further used the anisotropic spatial-spectral TV to model the global spatial and spectral smoothness and consistency of HSI, and proposed the local low-rank matrix recovery and global spatial-spectral TV (LLRSSTV) model, which can effectively separate the sparse noise.
However, these above low rank-based methods usually need to rearrange the HSI cube into the matrix for low-rank modeling, it is easy to destroy the correlation between spectral bands. To overcome these drawbacks, treating HSI as a third-order tensor had been recently and widely studied to maintain the cube structure of HSI. Thus, many tensor low-rank-based HSI restoration methods had been proposed, such as the low-rank tensor approximation (LRTA) model [33], TV regularized low-rank tensor decomposition (LRTDTV) [34], double-factor-regularized low-rank tensor factorization (LRTFDFR) [35], low-rank tensor dictionary learning denoising method [36], weighted tensor lowrank restoration (WLRTR) model [37], global spatial-spectral TV regularized nonconvex local low-rank tensor approximation (LLxRGTV) [38], etc.
Among the tensor decomposition-based HSI denoising methods, the Tucker decomposition [39] and CANDE-COMP/PARAFAC (CP) decomposition [40] schemes are commonly used to model the low-rank property of tensors, such as rank-1 tensor decomposition (R1TD) [41], nonlocal low-rank regularized R1TD (NLR-R1TD) [42], nonlocal similarity-based nonnegative Tucker decomposition [43], PARAFAC model [44], and nonlocal low-rank regularized tensor CP decomposition [45]. However, these two decompositions cannot give the best approximation to the tensor, Chen et al. [46] further proposed the nonlocal tensor-ring (TR) approximation, and exploited the nonlocal self-similarity and global spectral correlation of HSI at the same time to improve the restoration performance. More specifically, different from CP decomposition and Tucker decomposition, Kilmer et al. [47] defined tube rank and multirank of tensor based on tensor singular value decomposition (t-SVD), which can better describe the low-rank property of tensor. On this basis, by applying the discrete Fourier transform (DFT) along the mode-3 of the tensor, Zhang et al. [48] then defined the tensor nuclear norm (TNN) as the sum of singular values of all the frontal slices as the convex approximation of tensor rank function. However, TNN only describes the low rankness of spectral mode, but ignores the spatial modes. To model the low rankness of the spectral and spatial modes of HSI at the same time, Zheng et al. [49] extended t-SVD to mode-k t-SVD, then exploited the tensor low-rank property of HSI via the DFT along spectral and spatial modes, and particularly proposed the tensor fiber rank model and its convex relaxation model, namely, the 3-D tensor nuclear norm (3DTNN) model. Furthermore, some HSI restoration methods based on low-rank modeling in the transformed domain are also proposed, such as the framelet-based three-modal tensor nuclear norm (F-3MTNN) model [50] and the factor group sparsity-regularized nonconvex low-rank approximation (FGSLR) model [51]. What is more, some new low-rank tensor representation models are recently proposed for HSI restoration especially for HSI completion and cloud removal, such as the multilayer sparsity-based tensor decomposition (MLSTD) model [52], the parametric tensor sparsity model based on Laplacian scale mixture modeling via three-layer transform (LSM-TLT) [53], and the multimodal core tensor factorization (MCTF) model [54].
More recently, the deep learning-based HSI restoration approaches have attracted more attention and became more popular, such as the convolutional neural network (CNN)-based HSI denoising method (HSI-DeNet) [55], which modeled HSIs as tensors and well treated various noise simultaneously, and the HSI denoising with spatial-spectral deep residual CNN method (HSID-CNN) [56], which used spatial-spectral multiscale feature extraction and deep CNN-based residual learning to obtain a nonlinear mapping between the clean and noisy HSIs.
More clearly, the tensor construction and its tensor low-rank property modeling of HSI are the core of the tensor low-rank prior-based HSI restoration methods. Thus, we also focus on the tensor low-rank-based prior methods for HSI restoration, and particularly aim to investigate more effective scheme for the tensor construction and tensor low-rank modeling of spectral and spatial modes in this article.
To the best of our knowledge, there are some methods based on the gradient-based low-rank modeling for image processing. For example, the enhanced 3DTV (E-3DTV) model [57] has exploited the low-rankness of unfolded spatial gradients, i.e., the matrix-based low-rank property, not in the tensor form, and the gradient-based low-rank approximation (Grad-LR) model [58] just exploits the nonlocal 2-D patched-based low rankness of the spatial gradients for 2-D image inpainting, not in the tensor form either. More clearly, there are no relative works to investigate the spatial gradient tensors based low-rank properties in the tensor form for the HSI restoration recently. Thus, the low-rank modeling of HSI in the spatial gradient domain in the tensor form is highly valuable to be exploited for spatial-spectral structure preservation.
Moreover, by comparing with the strong tensor fiber lowrank property of spatial mode, the tensor fiber low-rank property of spectral mode are relatively unsuitable for spectral low-rank modeling of HSI in the 3DTNN method [49], owing to the HSI clearly shows strong spectral low-rank property of spectral mode in the original domain [19], [30]. Thus, the spectral low-rank property of spectral mode of HSI needs to be suitably modeled.
Based on the above analysis, instead of directly treating the HSI as a third-order tensor for tensor low-rank modeling in the original domain, we thus construct the spatial gradient tensors and the spectral gradient tensor for HSI, and particularly investigate the low-rank properties of HSI in gradient domains in this article. Specifically, the main contributions of this article are described as follows: 1) For the HSI, we first obtain its spatial horizontal gradient tensor and spatial vertical gradient tensor in spatial gradient domain, and the spectral gradient tensor in spectral gradient domain. Then, for the spectral low-rank prior modeling of HSI in spectral gradient domain, we particularly investigate the low-rank property of the spectral gradient of HSI to establish the spectral gradient-based weighted nuclear norm low-rank prior term (see Fig. 1), which can effectively remove the structural noise, such as stripes and deadlines.
2) For the spatial low-rank prior modeling of HSI in spatial gradient domain, we consider to use the DFT along the spatial mode-1 of spatial horizontal gradient tensor and the spatial mode-2 of spatial vertical gradient tensor of HSI, and obtain the spatial mode-1 gradient frequency tensor and the spatial mode-2 gradient frequency tensor, respectively. Furthermore, the lowrank property of the spatial mode-1 gradient frequency tensor and the spatial mode-2 gradient frequency tensor is modeled by the unified spatial gradient-based TNN low-rank prior term (see Fig. 1).
3) Consequently, based on the spectral gradient-based weighted nuclear norm low-rank prior term and the spatial gradient-based TNN low-rank prior, we propose a novel unified tensor-based HSI restoration model with low-rank modeling in gradient domains in this article.
The rest part of the article are listed as follows. We proposed the unified tensor-based HSI restoration model with low-rank modeling in gradient domains in Section II, and solved the proposed model in Section III. Section IV showed the overall experimental results and analysis. Section V provided the discussion in detail. Finally, the conclusion was given in Section VI.

II. PROPOSED MODEL
In this section, a unified tensor-based HSI restoration model with low-rank modeling in gradient domains will be proposed. For better modeling and analyzing, we first give some important symbols in Table I.

A. Formulation of HSI Restoration Problem
In fact, HSI is a cube data, which can be modeled as a third-order tensor. More clearly, let X ∈ n 1 ×n 2 ×n 3 to be the original clean HSI, and Y ∈ n 1 ×n 2 ×n 3 to be the noisy HSI with Gaussian noise and sparse noise (such as impulse noise, stripes, and deadlines), where n 1 × n 2 is the size of spatial dimension, and n 3 is the number of spectral bands. Therefore, the HSI with Gaussian noise and sparse noise contaminated can be expressed as where N ∈ n 1 ×n 2 ×n 3 and S ∈ n 1 ×n 2 ×n 3 are the Gaussian noise and sparse noise, respectively. The goal of HSI restoration is to restore the clean HSI X from the noisy HSI Y. By imposing the effective prior knowledge of the clean HSI X and the sparse noise S, thus, the image prior-based HSI restoration model can be expressed as where 1 2 Y − X − S 2 F is the fidelity term modeling the Gaussian noise N , R (S) = S 1 is the sparse prior regularization term modeling the sparse noise S, J(X ) is the prior regularization term modeling the spatial and spectral properties of the clean HSI X , λ 1 and β are the tradeoff parameters used to balance the two prior regularization terms.

1) Motivation:
Clearly, we know that the common idea for spectral low-rank modeling used in the previous methods is that the spectra of HSI X ∈ n 1 ×n 2 ×n 3 is low rank, namely, the matrix unfold 3 (X ) ∈ n 3 ×n 1 n 2 is low rank [19], [30], and then the nuclear norm-based spectral low-rank prior is usually enforced as Moreover, we also considered that the previous spatialspectral smoothness-based methods usually used the spectral gradient-based term ∇ 3 X 1 to enforce the sparsity of spectral gradient for spectral information preserving, where ∇ 3 is the spectral gradient operator. Thus, we mainly investigate the lowrank prior modeling of HSI in spectral gradient domain in this article.
Inspired by them, we mainly aim to model the spectral smoothness and spectral low-rankness of HSI in the spectral gradient domain at the same time in this article. To this end, we particularly investigate the low-rank prior modeling of the spectral gradient of HSI and establish the spectral gradient-based weighted nuclear norm low-rank prior model.
2) Spectral Gradient-Based Weighted Nuclear Norm Low-Rank Prior Term: Specifically, ∇ 3 represents the spectral gradient operator along the mode-3 spectral dimension of HSI X , then the spectral gradient of HSI X is expressed as As clearly shown in Fig. 1, for HSI X ∈ n 1 ×n 2 ×n 3 , we can obtain its spectral gradient tensor ∇ 3 X ∈ n 1 ×n 2 ×n 3 . Then, by arranging ∇ 3 X into a matrix unfold 3 (∇ 3 X ) ∈ n 3 ×n 1 n 2 via the mode-3 unfolding operation, and performing singular value decomposition (SVD) on unfold 3 (∇ 3 X ), we can clearly obtain the low-rank property of unfold 3 (∇ 3 X ) in Fig. 1. Therefore, the low-rank property of unfold 3 (∇ 3 X ) is particularly modeled by the following spectral gradient-based weighted nuclear norm low-rank prior term, which is formulated as ε is a very small positive scalar, and λ 2 ≥ 0.
C. Low-Rank Modeling in Spatial Gradient Domain 1) Motivation: First, we considered that the 3DTNN method [49] exploited the tensor low-rank property of HSI in original domain by applying the DFT along each mode of HSI X to form the DFT-based HSIs X 1 = fft(X , [ ], 1), X 2 = fft(X , [ ], 2), and 3), and obtained that the slices of each mode of the DFT-based HSIs are low rank. Clearly, the 3DTNN prior term was proposed as Moreover, we also considered that the previous spatialspectral smoothness-based methods usually used the spatial gradients-based prior term 2 i ρ i ∇ i X 1 to enforce the sparsity of spatial gradients, where ∇ 1 and ∇ 2 are the spatial horizontal and vertical gradient operators. Specifically, as Section I mentioned, the enhanced 3DTV (E-3DTV) model [57] has exploited the low rankness of unfolded spatial gradients, i.e., the matrix-based low-rank property, not in the tensor form, and the gradient-based low-rank approximation (Grad-LR) model [58] just exploits the nonlocal 2-D patched-based low rankness of the spatial gradients for 2-D image inpainting, not in the tensor form either. More clearly, there are no relative works to investigate the spatial gradient tensors-based low-rank properties in the tensor form for HSI restoration recently. Thus, the low-rank modeling of HSI in the spatial gradient domain in the tensor form is highly valuable to be exploited for spatial-spectral structure preservation.
Thus, for the spatial-mode low-rank modeling, instead of directly treating the HSI as a third-order tensor for tensor lowrank modeling of 3DTNN method [49] in the original domain, we mainly investigate the spatial gradient tensors and tensor low-rank priors of HSI in the spatial gradient domain in this article. Specifically, inspired by the 3DTNN method, we still apply the DFT but to the spatial gradients of HSI, and obtain that the spatial gradient tensors in the DFT domain are also low rank, and particularly propose the following spatial gradient-based TNN low-rank prior term.
2) Spatial Gradient-Based TNN Low-Rank Prior Term: Thus, for the tensor construction, we first apply the spatial horizontal gradient operator ∇ 1 and the spatial vertical gradient operator ∇ 2 to the HSI X ∈ n 1 ×n 2 ×n 3 , and obtain its spatial horizontal gradient tensor ∇ 1 X ∈ n 1 ×n 2 ×n 3 and spatial vertical gradient tensor ∇ 2 X ∈ n 1 ×n 2 ×n 3 , which are defined as Then, we use the DFT along the spatial mode-1 of spatial horizontal gradient tensor ∇ 1 X and the spatial mode-2 of spatial vertical gradient tensor ∇ 2 X , respectively. Thus, we obtain the spatial mode-1 horizontal gradient frequency tensor ∇ 1 X 1 = fft(∇ 1 X , [ ], 1) ∈ n 1 ×n 2 ×n 3 and the spatial mode-2 vertical gradient frequency tensor ∇ 2 X 2 = fft(∇ 2 X , [ ], 2) ∈ n 1 ×n 2 ×n 3 , respectively. Finally, we particularly exploit the tensor low-rank property of ∇ 1 X 1 and ∇ 2 X 2 in Fig. 1. As clearly shown in Fig. 1, the slices of ∇ 1 X 1 along mode-1 and the slices of ∇ 2 X 2 along mode-2 are low rank, respectively. Thus, the low rankness of the spatial mode-1 horizontal gradient frequency tensor ∇ 1 X 1 and the spatial mode-2 vertical gradient frequency tensor ∇ 2 X 2 is modeled by the mode-1 and mode-2 TNNs, respectively. Furthermore, we propose the unified spatial gradient-based TNN low-rank (SGTNNLR) prior term, which is formulated as where ∇ 1 X 1 (i, :, :) and ∇ 2 X 2 (:, i, :) are the ith mode-1 slice of ∇ 1 X 1 and ith mode-2 slice of ∇ 2 X 2 , respectively.

D. Proposed Model
Based on above model (2), we combine the spectral gradientbased weighted nuclear norm low-rank prior term (5) and the spatial gradient-based TNN low-rank prior term (9), and propose the unified tensor-based HSI restoration model with low-rank modeling in gradient domains as III. ALGORITHM Clearly, we here apply the ADMM method to solve the proposed model (11). By setting three auxiliary variables L, W 1 , and W 2 , the model (11) is equivalently transformed into the following form: Then, its augmented Lagrange function is expressed as 13) where η, μ 1 , and μ 2 are the penalty parameters and U , H 1 , and H 2 are the Lagrange multipliers.
Thus, by initializing k , (k = 1, 2) and using the iterative procedure, the optimization of (13) is decomposed into the optimization of the following four subproblems under the ADMM framework.
1) Optimization of S: The subproblem of S (i+1) is given as which can be efficiently solved by the soft-thresholding method as where shrink(Q, δ)= sign(Q) max(0, |Q| − δ) is the softthresholding operator with threshold value δ.

A. Experiment Setting
In order to verify the effectiveness of our proposed tensorbased HSI restoration method with low-rank modeling in gradient domains, we have carried out many denoising experiments on various simulated and real HSI datasets by comparing with some state-of-the-art methods. In this article, two simulated datasets and three real datasets are used for the HSI denoising experiments.
All the experiments are carried out on Windows 11 and MATLAB(R2016a), which uses Intel Core i5-10210 U 1.6 GHz and 16 GB RAM.

B. Experiments With Simulated Datasets
In this section, the visually and quantitatively experimental results of simulated HSI datasets are clearly given.  Then, three quantitative image quality indexes, namely, the mean peak signal-to-noise ratio (MPSNR), mean structural similarity (MSSIM), and spectral angle mapper (SAM) are used to evaluate the quality of denoised HSIs.
Specifically, we add six cases of different noise to the clean HSI to generate the simulated PaviaC and WDC datasets, which are the following six cases.
Case 1: (Gaussian Noise) Zero-mean Gaussian noise with different signal-to-noise ratio randomly selected between 10 and 20 dB is added to each band of the PaviaC and WDC dataset.
Case 2: (Gaussian Noise + Salt and Pepper) Gaussian noise is added in the same way as Case 1, and then the salt and pepper noise with a density of 20 is randomly added to 20 bands.
Case 3: (Gaussian Noise + Deadlines) Add Gaussian noise in the same way as Case 1, and then randomly select 10 bands to add the deadlines. Specifically, Fig. 2 shows the denoising results of all methods on the given band 41, band 49, band 1, and band 30 of the simulated PaviaC dataset under Case 1, Case 2, Case 4, and Case 6, respectively.
As shown in the top-row of Fig. 2, i.e., band 41 of the simulated PaviaC dataset under Case 1 with Gaussian noise, all the methods can effectively remove the Gaussian noise, but LRMR, LRTV, LRTDTV, 3DTNN, LRTFDFR, ITSR, and NG-meet methods more or less blur the details of the image. Moreover, FGSLR 1/2 and the proposed method restore the original clean image to the greatest extent while better preserving the image details and structures.
As shown in the second-row of Fig. 2, i.e., band 49 of the simulated PaviaC dataset under Case 2 with Gaussian and Salt and Pepper Noise, we can clearly see that ITSR can not effectively remove the noise, and LRMR still has a small amount of noise. LRTV, LRTDTV, 3DTNN, LRTFDFR, NG-meet, and FGSLR 1/2 can well remove the noise, but still have some shortcomings in preserving the image details, which slightly blur few details of the image at different degrees. By contrast, the proposed method shows better denoising performance, which can better remove the noise and at the same time restore the image texture structure and details finer.
As shown in the third-row of Fig. 2, i.e., band 1 of the simulated PaviaC dataset under Case 4 with Gaussian Noise and Stripes, LRMR, 3DTNN, ITSR, and NG-meet can not effectively remove the stripe noise, whereas LRTV, LRTDTV, and LRTFDFR can well remove the stripe noise, but still blur some image details at different degrees. The denoising result of FGSLR 1/2 is satisfactory, but a few local details are not well preserved in the enlarged region. What is more, the proposed method performs much better results in removing stripe noise and preserving image details.
As shown in the bottom-row of Fig. 2, i.e., band 30 of the simulated PaviaC dataset under Case 6 with Gaussian noise, salt and pepper, deadlines, and stripes, four kinds of mixed noise are added to the original image to generate the simulated noisy image. However, the mixed noise can not be effectively removed by LRMR, 3DTNN, ITSR, and NG-meet, and their restored results are not satisfactory. Specifically, LRMR, 3DTNN, ITSR, and NG-meet can not well remove the stripe and deadline noise. LRTV, LRTDTV, and LRTFDFR can well remove the stripe and deadline noise, and they also have some information loss and blur artifacts at different degrees. FGSLR 1/2 can obtain very satisfactory restored result. Compared with other methods, the proposed method achieves the best visual performance, which can best remove the stripe and deadline noise, and also preserve the image details finest.
Furthermore, Fig. 3 shows the denoising results on the given band 43, band 31, band 22, and band 37 of the simulated WDC dataset under Case 1, Case 2, Case 4, and Case 6, respectively.
As shown in the top-row of Fig. 3, i.e., band 43 of the simulated WDC dataset under Case 1 with Gaussian noise, all the methods can well remove the Gaussian noise. Meanwhile, the results of LRTV, LRTDTV, and ITSR are too smooth. LRMR, 3DTNN, LRTFDFR, NG-meet, and FGSLR 1/2 all have the problem of loss of image details in different degrees. The proposed method can better preserve the global and local details of the image, and has the best visual effect.
As shown in the second-row of Fig. 3, i.e., band 31 of the simulated WDC dataset under Case 2 with Gaussian and Salt and Pepper Noise, the image restored by ITSR still contains a small amount of noise. LRTV and LRTDTV blur more details of image, LRMR, 3DTNN, LRTFDFR, and NG-meet blur less details of image. The denoising result of FGSLR 1/2 is satisfactory, but it can be seen from the enlarged region that the details of the image are also slightly blurred. Moreover, the proposed method shows better results in removing noise and preserving details.
As shown in the third-row of Fig. 3, i.e., band 22 of the simulated WDC dataset under Case 4 with Gaussian Noise and Stripes, 3DTNN and ITSR can not well remove the stripes, and the denoising results of LRTV and LRTDTV are not very ideal, thus losing some image details. LRMR, LRTFDFR, NG-meet, FGSLR 1/2 and the proposed method show good performance in removing the stripes, and the image recovered by the proposed method is the closest to the original HSI.
As shown in the bottom-row of Fig. 3, i.e., band 37 of the simulated WDC dataset under Case 6 with Gaussian noise, salt and pepper, deadlines, and stripes, ITSR only removes a small portion of the noise. Moreover, ITSR and 3DTNN can not well remove the deadlines, and stripes, and LRMR and NG-meet can remove some deadlines and stripes but also remain few stripes in the denoised images with a much closer look. LRTV, LRTDTV, LRTFDFR, and FGSLR 1/2 can well remove the deadlines and stripes, thus showing satisfactory results. Clearly, the visual performance of the proposed method is the most satisfactory, which can best remove the stripe and deadline noise, and preserve image details at the same time.
In summary, for above visual comparisons of different cases of two simulated PaviaC and WDC datasets, the denoising results of our proposed method are still excellent. In all cases, it can well remove the mixed noise, restore the clean image, and preserve the original information of the image. With the variety of mixed noise added, the superiority of the proposed method is more obvious and robust in removing the deadlines and stripes.  Tables II and III, the best results of each quality index are shown in bold, and the second best results are shown in underline.
As shown in Table II, the MPSNR, MSSIM, and SAM values of the proposed method are the best for all these six cases, which thus show that the proposed method achieves the best denoising results on the simulated PaviaC dataset.
As shown in Table III, the proposed method achieves the best results for most of the cases. Clearly, the proposed method gives the best MPSNR and MSSIM results for Case 1, Case 2, and Case 4, and the second best SAM result for Case 4, where 3DTNN gives the best SAM results for Case 1 and Case 2, and FGSLR 1/2 gives the best SAM result for Case 4. What is more, the proposed method consistently gives the best MPSNR, MSSIM, and SAM  results for the rest Case 3, Case 5, and Case 6. Moreover, in most cases of the simulated WDC dataset, the proposed method is much better than the other compared methods.
More specifically, Figs. 4 and 5 show the PSNR and SSIM values of each band of different methods on the simulated PaviaC and WDC datasets under Case 1, Case 2, Case 4, and Case 6, respectively. It can be seen from Figs. 4 and 5 that the proposed method achieves the best PSNR and SSIM values in most bands of the simulated PaviaC and WDC datasets. It is worth noting that the curves of PSNR and SSIM of the proposed method are smooth and stable, while the other methods cause the large fluctuations owing to the deadlines or stripes noise at these bands are not well removed. Therefore, the proposed method has better denoising ability to remove deadlines and stripes in mixed noise.   In general, by comparing with the other methods on the simulated PaviaC and WDC datasets, the proposed method has advantages in both visual evaluation and quantitative evaluation. Thus, it further proves that the proposed method is more effective and feasible for removing mixed noise of HSI.  , by comparing with the sparse noise S of noisy HSI, (i.e., the original sparse noise S as reference), the proposed method can best restore the sparse noise S among these compared methods.

5) Comparison With Deep Learning Method:
Moreover, to further illustrate the effectiveness of the proposed method against the deep learning method, we also compare our proposed method with the HSID-CNN method [56]. As we clearly know, the HSID-CNN method is just proposed for the additive Gaussian noise removal from HSIs [56]. Thus, to make the comparison more fair, we particularly use the aforementioned simulated PaviaC dataset of size 200 × 200 × 80 with the additive Gaussian noise of five different variances σ 2 n , i.e., σ n = 5, σ n = 25, σ n = 50, σ n = 75, and σ n = 100, which are the same as the simulated experiments conducted by the HSID-CNN method. In this sense, we give the correspondingly quantitative results in Table IV. Clearly, our proposed method gives better MPSNR, MSSIM, and SAM results than those of HSID-CNN for most of the cases, which shows that our proposed method can give better denoising performance than HSID-CNN generally.

C. Experiments With Real Datasets
In this section, the experimental results of real HSI datasets are also given.
Moreover, the denoising results of different methods on the real Urban, EO-1, and PaviaU datasets are shown in Figs. 10-12, respectively.
As shown in Fig. 10, i.e., band 108 of the real Urban dataset, it is obvious that the original band 108 of Urban is very blurred after being polluted by the noise. From the whole parts and the red enlarged parts of the denoised images, we can clearly see that LRMR and LRTFDFR have only removed a little amount of noise, and their denoised images are still blurred. LRTV, LRTDTV, 3DTNN, ITSR, NG-meet, and FGSLR 1/2 remove more of the noise, and the objects in the red box become also visible. However, compared with the proposed method, the restored results of the LRTV, LRTDTV, 3DTNN, ITSR, NG-meet, and FGSLR 1/2 methods are still not satisfied, which also lose many image details. More clearly, the proposed method can best remove the noise and preserve image details. Therefore, the proposed method shows the best visual denoising performance on the real Urban dataset.
Then, Fig. 11 shows the denoising results of the band 95 of the real EO-1 dataset. From the whole parts and the enlarged red box of the images, we can clearly see that the denoising performance of LRMR, LRTV, and LRTFDFR on the real EO-1 dataset is not ideal, and their restored images are still blurred and noisy. LRTDTV, 3DTNN, ITSR, NG-meet, and FGSLR 1/2 can remove more of the noise. However, by comparing with the proposed method, LRTDTV, 3DTNN, ITSR, NG-meet, and FGSLR 1/2 can not well preserve the details of the image, which are also some blurred and are not as fine as the restored image of the proposed method. Clearly, the proposed method can best remove the noise and preserve image details finest. Therefore, the denoising performance of the proposed method on the real EO-1 dataset is the best.
Moreover, Fig. 12 shows the denoising results of the band 3 of the real PaviaU dataset, where the denoising performance of each method is similar with the cases of real Urban and EO-1 datasets.
Meanwhile, we use the no-reference HSI quality assessment (NHSIQA) index [62] for quantitatively assessing the real data denoising experiments. Then, Table V gives the NHSIQA results of different approaches on the real Urban, EO-1, and PaviaU datasets, where the smaller NHSIQA stands for the higher quality of HSI. As given in Table V, the proposed approach consistently gives the smallest NHSIQA results, i.e., the highest quality of restored HSI.    In conclusion, for these results on the real Urban, EO-1, and PaviaU datasets, the proposed method shows the best visual and quantitative denoising performance.

A. Parameter Analysis
As we know, the results of the prior-based denoising models are usually influenced by the model parameters, and choosing the appropriate parameters also has very important influence to the final denoising results.
Thus, Fig. 13 shows the sensitivity analysis of the proposed method by studying the MPSNR and MSSIM results versus parameters α 1 , α 2 , λ 1 , and λ 2 particularly on the simulated PaviaC dataset under Case 6. It can be clearly seen from Fig. 13 that the different setting of parameters α 1 , α 2 , λ 1 , and λ 2 can affect the MPSNR and MSSIM results of the proposed method. Based on this, for the setting of parameters α 1 , α 2 , λ 1 , and λ 2 in our experiments, it is suggested that the values of α 1 , α 2 , λ 1 , and λ 2 should be selected in the intervals of Moreover, we also discuss the influence of the proposed spatial gradient-based TNN low-rank (SGTNNLR) prior term (9) with different setting of α 1 and α 2 on the simulated PaviaC dataset under Case 6, which aims to analyze how the spatial horizontal gradient-based tensor low-rank term ∇ 1 X TNN 1  and the spatial vertical gradient-based tensor low-rank term ∇ 2 X TNN 2 affect the denoising performance separately.
As we can see from the results of the above three methods on the simulated PaviaC dataset under Case 6 given in Table VI, the SGTNNLR term with both spatial horizontal and vertical gradients-based tensor low-rank terms (i.e., SGTNNLR with  α 1 = 0 and α 2 = 0) gives the best results, and the spatial horizontal gradient-based tensor low-rank term (i.e., SGTNNLR with α 2 = 0) and the spatial vertical gradient-based tensor lowrank term (i.e., SGTNNLR with α 1 = 0) can both affect the denoising performance at different degrees, which both play important roles in the proposed SGTNNLR prior term (9) for HSI restoration.

C. Influence of the Spectral and Spatial Gradient Transforms of the Proposed Model
Clearly, the low-rank modeling of proposed model is investigated in the spectral and spatial gradient domains. Thus, we further discuss the influence of spectral and spatial gradient transforms of the proposed model (11) on the simulated PaviaC dataset under Case 6.
As clearly given in Table VII, Model 3, i.e., the proposed model (11) with both the spectral and spatial gradient transforms, gives the best results, and the spectral and spatial gradient transforms can both affect the performance at different degrees, which shows than the spectral and spatial gradient transforms both play important roles in the proposed model (11) for the HSI restoration.

D. Influence of the DFT of Proposed SGTNNLR Prior Term (9)
Moreover, we discuss the influence of DFT of the proposed SGTNNLR prior term (9) on the simulated PaviaC dataset under Case 6, where the following two methods are considered: i) SGTNNLR without DFT, and ii) SGTNNLR with DFT.
As clearly provided in Table VIII, SGTNNLR with DFT give better results than those of SGTNNLR without DFT, which thus shows that the DFT also plays an important role in the proposed SGTNNLR prior term (9) for HSI restoration.

E. Convergence Analysis
Moreover, the convergence of the proposed method is also analyzed. To this end, the MPSNR and MSSIM results versus iterations on the simulated PaviaC and WDC datasets under Case 6 are representatively provided in Fig. 14 for the convergence analysis.
As obviously shown in Fig. 14, the values of MPSNR and MSSIM first increase with the increase of iterations, and then reach a relatively stable state, which is convergent finally. Therefore, the convergence of the proposed method can be further validated.

F. Analysis of Computational Complexity
Furthermore, the computational complexity of our proposed Algorithm 1 is also analyzed, which can be separately analyzed by the following five subproblems: 1) For the solution of S (i+1) subproblem via (15), it has the complexity of O(n 1 n 2 n 3 ) for each iteration.

G. Comparisons of CPU Time
Finally, the detailed comparisons of CPU time of different approaches on different HSI datasets are given in Table IX. As given in Table IX, LRMR always spends the shortest CPU time, whereas ISTR always spends the longest CPU time. Clearly, since the proposed model is much more complex than LRMR, then the proposed method spends much more CPU time than LRMR, but much less than ISTR, which shows some competitiveness with the other approaches.

VI. CONCLUSION
In this article, we have proposed a novel and effective tensorbased HSI restoration model with low-rank modeling in gradient domains by using the tensor modeling theory. Clearly, we propose the spectral gradient-based weighted nuclear norm low-rank prior term for the spectral low-rank modeling of HSI in spectral gradient domain. More specifically, we particularly propose the spatial gradient-based TNN low-rank prior term for the spatial low-rank modeling of HSI in spatial gradient domain. Moreover, the proposed model is effectively optimized by the ADMM method. Finally, a lot of denoising results on both simulated and real HSI datasets clearly demonstrate that the proposed method is superior to many state-of-the-art denoising methods in the aspects of visual and quantitative comparisons. More importantly, the proposed method is much more effective and robust for the removal of the stripes and deadlines.