Abstract
Hyperspectral images (HSIs) can help deliver more reliable representations of real scenes than traditional images and enhance the performance of many computer vision tasks. However, in real cases, an HSI is often degraded by a mixture of various types of noise, including Gaussian noise and impulse noise. In this paper, we propose a logarithmic nonconvex regularization model for HSI mixed noise removal. The logarithmic penalty function can approximate the tensor fibered rank more accurately and treats singular values differently. An alternating direction method of multipliers (ADMM) is also presented to solve the optimization problem, and each subproblem within ADMM is proven to have a closedform solution. The experimental results demonstrate the effectiveness of the proposed method.
1. Introduction
A hyperspectral image (HSI) consists of multiple discrete bands at specific frequencies. It can deliver additional information that the human eye fails to capture for real scenes and has been attracting much interest from researchers in a wide range of application fields, such as land use analysis, environmental monitoring, and field surveillance [1–3]. However, HSIs always suffer from various degradations, such as Gaussian noise, impulse noise, and random noise, which can affect the subsequent image processing, such as unmixing, classification, and target detection [4, 5]. Improving the HSI quality merely through a hardware scheme is unsustainable and impractical. Therefore, it is natural to introduce image processingbased approaches to obtain a highquality HSI before subsequent applications.
The numerical optimization algorithm plays an important role in HSI denoising, such as Liu et al. [6] proposed a twostep waveletdomain estimation method to extract the noise map, and Lu et al. [7] presented some representative highorder variational models in the context of image denoising. From the perspective of the prior data format, we classify the existing HSI restoration methods into three categories: (1) 1D vectorbased sparse representation methods [8–13]; (2) 2D matrixbased lowrank matrix recovery (LRMR) methods [14–21]; and, (3) 3D tensorbased approximation methods [22–33]. Although the existing works have made significant progress in HSI restoration, there are still drawbacks that need to be improved, such as when the multidimensional HSI data is transformed into a vector or matrix, it usually breaks the spectralspatial structural correlation. The tensor lowrankness characterization of an HSI is expected to explore the global correlation and preserve the intrinsic structural information. The tensors’ recovery under a limited number of measurements is an important problem that arises in a variety of applications, such as computer vision [34–36]. Based on low tubalrank tensor recovery, Zhang et al. [37] proposed an HSI mixed noise removal model. However, the framework of the tensor singular value decomposition (tSVD) lacks flexibility for handling different correlations along the different modes of HSIs, leading to suboptimal denoising performance. Then, Zheng et al. [38] proposed an HSI mixed noise removal model with tensor fibered rank, which is based on the modek tSVD. Moreover, Zheng et al. [38] also proposed a threedirectional tensor nuclear norm (3DTNN) as its convex relaxation to provide an efficient numerical solution and a threedirectional logbased tensor nuclear norm (3DLogTNN) as its nonconvex relaxation to promote the low rankness of the solution. Compared to 3DTNN, 3DLogTNN has two advantages. First, it is a closer approximation to the fibered rank than 3DTNN. Second, by using the sum of the log function of singular values, 3DLogTNN can better approximate to the fibered rank than 3DTNN.
It is well known that suitable nonconvex penalty functions induce sparsity among the singular values more effectively. However, the use of nonconvex penalty functions will lead to nonconvex optimization problems. Then, it suffers from numerous issues such as spurious local minima in the subproblem, for example, 3DLogTNN in [38]. To avoid the intrinsic difficulties, we introduce a new nonconvex logarithmic regularization model, which allows the use of nonconvex penalty function while maintaining convexity of the subproblem within ADMM. Also, the new model can provide a good approximation for the fibered rank and preserve the major information well.
The rest of the paper is structured as follows. Section 2 presents some notations, explains tSVD and defines modek tSVD. Section 3 introduces the proposed ADMM based on the parametric penalty function for HSI denoising. The experimental results and analysis are reported in Section 4. Finally, the conclusion is given in Section 5.
2. Brief Overview of Tensor Singular Value Decomposition
In this section, we first describe the notations used throughout the paper and then introduce the tensor decomposition proposed in [39–41] and modek tSVD proposed in [38].
2.1. Notation and Indexing
For a thirdorder tensor , denote as the horizontal, lateral, and frontal slices, respectively. denote the mode1, mode2, and mode3 fibers, respectively; Figures 1 and 2 show the horizontal, lateral, and frontal slides, denoted by , respectively, of a thirdorder tensor . Suppose ; we adopt the definition of the Frobenius norm of a tensor and the definition of the norm of a tensor, .
(a)
(b)
(c)
(a)
(b)
(c)
Let denote , that is, the frontal slice of . Then,
Similarly, , and
It can be seen in Figure 3, for , let ; then, are the frontal slices of tensor , where is computed by applying the Fast Fourier Transform (FFT) along each tube of , i.e., and .
2.2. tSVD
In this section, we exploit the proposed tSVD. A tSVD interprets thirdorder tensors as linear operators on the space of oriented matrices [39]. Based on a tSVD, O. Semerci exploited the decomposition, completion, and recovery of multilinear data [42], and Zhang et al. [34] proposed novel methods for completion and denoising of multilinear data and, as an application, considered 3D and 4D (color) video data completion and denoising.
Definition 1. (tproduct). For and , the tproduct is a tensor of size . For and ,The tproduct is analogous to matrix multiplication except that circular convolution replaces the multiplication operations between the elements, which are represented by tubes.
Theorem 1. (tSVD). For , the tSVD of is given bywhere and are orthogonal tensors, is a rectangular diagonal tensor, and denotes the tproduct.
Figure 4 illustrates the decomposition for the 3D case. Additionally, one can obtain this decomposition by computing matrix SVDs in the Fourier domain, see Algorithm 1.
However, when handling different correlations along different modes of an HSI, the tSVD and the induced tubal rank lack flexibility. This inflexible HSI characterization usually does not have good denoising effects. Zheng et al. [38] proposed a novel tensor decomposition by generalizing the tSVD to the modek tSVD, which achieves a more flexible and accurate HSI lowrankness characterization.

2.3. Modek tSVD
In this section, we introduce the modek tSVD and the related definitions.
For a thirdorder tensor , the modek block circulation operation is denoted aswhere is the modek slice of .
The modek block diagonalization operation and its inverse operation are defined as
The modek block vectorization operation and its inverse operation are defined as
Definition 2. (Modek tproduct). For and , the mode1 tproduct is a tensor of size :For and , the mode2 tproduct is a tensor of size :For and , the mode3 tproduct is a tensor of size :
Definition 3. (Modek identity tensor). is the modek identity tensor, whose first modek slice is an identity matrix and other modek slices are all zeros.
Definition 4. (Modek conjugate transpose). For is the modek conjugate transpose of , which is obtained by transposing each of the modek slices and then reversing the order of transposed modek slices 2 through .
Definition 5. (Modek diagonal tensor). For , each of its modek slices is a diagonal matrix, and then, is a modek diagonal.
Definition 6. (Modek orthogonal tensor). Ifwhere the tensor is modek orthogonal.
Definition 7. (Tensor modek permutation). For , is the modek permutation of , whose mode3 slice is the modek slice of , and its inverse operation is .
Theorem 2. (Modek tSVD). The factorization of iswhere and are the modek orthogonal tensors and is the modek diagonal tensor.
Theorem 2 is proven in [38] (Th. 2).
Definition 8. (Tensor fibered rank). For , is the fibered rank of , whose element is the modek tensor fibered rank . , where the element of is the rank of the modek slice of .
The modek tSVD can be efficiently obtained by computing a series of matrix SVDs in the Fourier domain and and achieves more flexible and accurate HSI lowrankness characterization (see Algorithm 2).

3. HSI Denoising Model and Its ADMM
In this section, we show our new denoising model and ADMM for solving the proposed model is also derived.
3.1. The Logarithmic Penalty Function
This article mainly proposes logarithmic penalty function [43] that serves as the model for the penalty function developed in the HSI denoising model below and is designed to have less bias. The logarithmic penalty is given bywhere controls the degree of nonconvexity of the penalty function. This function satisfies the following conditions:(A1)(A2)(A3)(A4)
The proximity operator associated with the nonconvex function is
In [44–46], the authors prove that, for , the function in equation (14) is convex. Therefore, the proximity operator finds an optimal solution for convex minimization problem (14). The proximity operator associated with logarithmic penalty equation (13) is a continuous nonlinear threshold function with as the threshold value [44], namely,and is given by
3.2. HSI Denoising Model
Zhang et al. proposed an HSI mixed noise removal model based on low tubalrank tensor recovery (LRTR) [37]. It can address the mixed noise in HSIs and decompose a noisy HSI into three parts, i.e., a lowtubalrank part (the clean HSI), a Gaussian noise part, and a sparse noise part. Zheng [38] proposed an HSI denoising model based on a lowfiberedrank prior, and it is formulated aswhere is the observed noisy HSI, is Gaussian noise, is sparse noise, is the clean hyperspectral signal, which has the lowfiberedrank tensor property, and and are regularization parameters.
The Gaussian noise model (large degree of freedom) corresponds to a dense noise type [47, 48]. should not be low in rank. Otherwise, cannot be recovered from random noise. Due to the stochastic nature of Gaussian noise, it is assumed that there is no correlation (or a weak correlation) among noise components. Thus, the rank of is normally full and much larger than the rank of . Therefore, the lowfiberedrank component can be recovered from equation (17) by properly choosing the tuning parameters and .
Directly minimizing the tensor fibered rank is NPhard. Based on logarithmic penalty equation (13), we propose that the following 3DNLogATNN HSI denoising model is commonly formulated:where and is the weight of the fibered rank. are set as in the HIS denoising model:where is the number of singular values of , is the singular values of , and is the modek slice of .
3.3. ADMM for Solving Model (18)
We use the ADMM to solve equation (18). By casting the three auxiliary variables , equation (18) can be rewritten as follows:
The augmented Lagrangian function of equation (18) iswhere and are the Lagrange multipliers and and are positive scalars. Now we can solve the problem within the ADMM framework.
With the other parameters fixed, can be updated by solvingwhich is equivalent to the following subproblem for :
To solve equation (23), we can rewrite it aswhere , and . From equation (16),where . Therefore,
Though model (20) is based on a nonconvex penalty function, the parameters of which are constrained to ensure the convexity of the subproblem .
With other parameters fixed, can be updated by solvingwhich is equivalent to the following subproblem:
It has the following closedform solution:
With the other parameters fixed, can be updated by solvingwhich is equivalent to the following subproblem:and it has the following closedform solution:
With the other parameters fixed, can be updated by solvingwhich is equivalent to the following subproblem:
Using the tensor soft thresholding operator, the following solution can be obtained [38]:
The Lagrange multipliers and can be updated by solving
Hence, the proposed algorithm for HSI denoising is summarized in Algorithm 3.

4. Experiment Results
To validate the effectiveness of the proposed method for HSI denoising, we perform experiments on simulated data and compare the experimental results both quantitatively and visually. The Washington DC Mall data, Pavia City Center data, and the Indian Pines data are used. In our experiments, the Washington DC Mall data uses only a subimage (191 bands and size of each band is ). The Pavia City Center data uses only a subimage (80 bands and size of each band is ). And, the synthetic data by Zhang et al. [37] was generated using the ground truth of the Indian Pines dataset; the size of the synthetic HSI was . In addition, to facilitate the numerical calculation and visualization, all the bands of the HSI are normalized into [0, 1], and they will be stretched to the original level after restoration.
The three evaluation measures are the mean peak signaltonoise ratio (MPSNR), mean structure similarity (MSSIM), and spectral angle mapping (SAM). The three metrics are defined as follows to measure the quality of the denoised result:where represents the image size in the space, represents the original image, represents the distortion image, is a pixel that can achieve the maximum value, and is the number of PSNR:where is a constant for is the same as and represent the and standard deviations, respectively, represent the original image and the distorted image, respectively, and represent the jth local window contents, and is the number of local windows:where is the number of wavelengths, and represent the spectrum vectors, and
The PSNR and structural similarity index measure (SSIM) are two conventional perceptual quality indexes (PQIs) in image processing and computer vision. They evaluate the similarities between the target image and the reference image based on the mean square error (MSE) and structural consistency. The larger these two measures are, the closer the target HSI is to the reference HSI. The SAM is a physically based spectral classification method that uses an Ndimensional angle to match pixels to reference spectra. Different from the former two measures, the smaller the SAM is, the more similar the target HSI is to the reference HSI.
Real HSIs usually include several different types of noise. To simulate realnoise scenarios, we consider the variance of the Gaussian noise and the variance of the impulse noise . We use statistical structures to simulate different types of noise, including independent and identically distributed (i.i.d.) and noni.i.d. noise, which are elaborated as follows:(1)Case 1 (noni.i.d. Gaussian noise): all the bands of the test dataset are contaminated by zeromean Gaussian noise with different intensities. The variance in the Gaussian noise is randomly changed from U (0.1, 0.2) and U (0.55, 0.65).(2)Case 2 (noni.i.d. impulse noise): in this case, all bands are contaminated by impulse noise with different ratios. The ratios of impulse noise are randomly changed from U (0.35, 0.45), U (0.45, 0.55), and U (0.55, 0.65).(3)Case 3 (i.i.d. Gaussian + i.i.d. impulse noise): in this case, all bands are corrupted by zeromean Gaussian noise and impulse noise. The variance in the Gaussian noise is 0.3, and the ratio of the impulse noise is 0.1; the variance in the Gaussian noise is 0.1, and the ratio of the impulse noise is 0.4; the variance in the Gaussian noise is 0.3, and the ratio of the impulse noise is 0.5.(4)Case 4 (noni.i.d. Gaussian + i.i.d. impulse noise): in this case, all bands are corrupted by zeromean Gaussian noise and impulse noise with different intensities. The variance in the Gaussian noise is randomly changed from U (0.3, 0.4), and the ratio of the impulse noise is 0.2; the variance in the Gaussian noise is randomly changed from U (0.2, 0.3), and the ratio of impulse noise is 0.3; the variance in the Gaussian noise is randomly changed from U (0.4, 0.5), and the ratio of impulse noise is 0.1.(5)Case 5 (i.i.d. Gaussian + noni.i.d. impulse noise): in this case, all bands are corrupted by zeromean Gaussian noise and impulse noise with different intensities. The variance in the Gaussian noise is 0.1, and the ratios of impulse noise are randomly varied from U (0.5, 0.6); the variance in the Gaussian noise is 0.3, and the ratios of impulse noise are randomly varied from U (0.3, 0.4); the variance in the Gaussian noise is 0.2, and the ratios of impulse noise are randomly varied from U (0.4, 0.5).(6)Case 6 (noni.i.d. Gaussian + noni.i.d. impulse noise): in this case, all bands are corrupted by zeromean Gaussian noise and impulse noise with different intensities. The variance in the Gaussian noise is randomly changed from U (0.2, 0.3), U (0.1, 0.2), and U (0.4, 0.5), and the ratios of impulse noise are randomly varied from U (0.2, 0.3), U (0.3, 0.4), U (0.1, 0.2), and U (0.4, 0.5).
Parameter setting: we analyze the parameters involved in the proposed method on HSIs’ Washington DC Mall, Pavia City Center, and Indian Pines, i.e., the weight , the regularization parameters and , the threshold parameter , the penalty parameter , and a constant in 3DLogATNN. In all the following experiments, the parameters in these compared methods were manually adjusted according to their default strategies.
The regularization parameter for 3DLogATNN: it is easy to see that is the parameter used to restrict the sparsity of the Gaussian noise. We set , where is the standard deviation of Gaussian noise and is a tuning parameter. The results were based on the simulated data experiment in case 11. Figure 5 shows the restoration results as varied in the set {0.0001, 0.0005, 0.0009, 0.001, 0.002, 0.003, 0.004, 0.005, 0.007, 0.009, 0.11, 0.5, 0.15}. It can be clearly observed from this figure that the results of the 3DLogATNN solver are relatively stable in terms of both MPSNR and MSSIM values, with the value of changing from 0.002 to 0.003. Therefore, we suggest the use of in all the simulated data experiments.
(a)
(b)
The regularization parameter for 3DLogATNN: it is easy to see that is the parameter used to restrict the sparsity of the impulse noise. We set , whereand is a tuning parameter. The results were based on the simulated data experiment in case 21. Figure 6 shows the restoration results as varied in the set {0.05, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9}. It can be clearly observed from this figure that the results of the 3DLogATNN solver are relatively stable in terms of both MPSNR and MSSIM values, with the value of changing from 0.4 to 0.5. Therefore, we suggest the use of in all the simulated data experiments.
(a)
(b)
The constant for 3DLogATNN: the parameter controls the degree of nonconvexity of the penalty function. The results were based on the simulated data experiment in case 21. Figure 7 shows the restoration results as varied in the set {0.01, 0.04, 0.05, 0.06, 0.1, 0.5, 1, 5, 10}. It can be clearly observed from this figure that the results of the 3DLogATNN solver are relatively stable in terms of both MPSNR and MSSIM values, with the value of changing from 0.05 to 0.06. Therefore, we suggest the use of in all the simulated data experiments.
(a)
(b)
We adjust the parameters to achieve the best visual result, and the parameter setting is presented in Table 1.
Compared with the stateoftheart methods, including TRPCA + BM4D [36, 49], LRMR [37], LRTR [37], LRTDTV [50], and NMoG [51], on lowrank matrix/tensor approximation and noise modeling, the extensive experimental results demonstrate that the 3DTNN and 3DLogTNN [38] methods are better at removing the mixed noise. Therefore, the denoising results of the proposed method are quantitatively and visually compared with two stateoftheart HSI denoising methods, i.e., 3DTNN and 3DLogTNN. The denoising results of all the methods in six cases are shown in Table 2. Three typical bands of the denoised HSIs in the mixture noise case obtained with different methods are shown in Figures 8–10. Figure 8 shows the denoising results at band 71 of the Washington DC Mall HSI, Figure 9 shows the denoising results at band 52 of the Pavia City Center HSI, and Figure 10 shows the denoising results at band 28 of the Indian Pines HSI. It can be seen that the proposed 3DLogATNN can effectively remove the mixed noise and preserve the detailed information of the original image. The proposed method obtains the best visual quality by removing all the mixture noise and preserving the details well. Table 2 shows that the 3DLogATNN method converges faster than the 3DLogTNNbased method on all the Washington DC Mall, Pavia City Center data, and Indian Pines data, and our method outperforms the compared ones for the Pavia City Center data. Besides, 3DLogATNN is more stable than the other two algorithms, because it can get the best results in most cases.
(a)
(b)
(a)
(b)
(a)
(b)
5. Conclusion
In this paper, we present a new 3DLogATNN method for HSI denoising by modek tSVD. The logarithmic penalty function is introduced in 3DLogATNN, which enables it to extract lowrank and sparse components more accurately from a degraded HSI contaminated by several types of noise. In addition, the ADMMbased algorithm is applied to effectively solve the proposed HSI denoising model. The experiments have substantiated the superiority of the proposed method over stateoftheart methods.
Data Availability
The data used to support the findings of this study are included within the article.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This work was supported by the National Natural Science Foundation of China (61967004, 11901137, 11961011, and 72061007), Guangxi Natural Science Foundation (2018GXNSFBA281023), China Postdoctoral Science Foundation (2020M682959), Guangxi Key Laboratory of Cryptography and Information Security (GCIS201621 and GCIS201927), Guangxi Key Laboratory of Automatic Detecting Technology and Instruments (YQ20113 and YQ20114), and Promotion Project of Basic Ability of Young and MiddleAged Teachers in Universities of Guangxi under Grant (2019KY0253).