Spectral-Smoothness and Non-Local Self-Similarity Regularized Subspace Low-Rank Learning Method for Hyperspectral Mixed Denoising

: During the imaging process, hyperspectral image (HSI) is inevitably affected by various noises, such as Gaussian noise, impulse noise, stripes or deadlines. As one of the pre-processing steps, the removal of mixed noise for HSI has a vital impact on subsequent applications, and it is also one of the most challenging tasks. In this paper, a novel spectral-smoothness and non-local self-similarity regularized subspace low-rank learning (termed SNSSLrL) method was proposed for the mixed noise removal of HSI. First, under the subspace decomposition framework, the original HSI is decomposed into the linear representation of two low-dimensional matrices, namely the subspace basis matrix and the coefﬁcient matrix. To further exploit the essential characteristics of HSI, on the one hand, the basis matrix is modeled as spectral smoothing, which constrains each column vector of the basis matrix to be a locally continuous spectrum, so that the subspace formed by its column vectors has continuous properties. On the other hand, the coefﬁcient matrix is divided into several non-local block matrices according to the pixel coordinates of the original HSI data, and block-matching and 4D ﬁltering (BM4D) is employed to reconstruct these self-similar non-local block matrices. Finally, the formulated model with all convex items is solved efﬁciently by the alternating direction method of multipliers (ADMM). Extensive experiments on two simulated datasets and one real dataset verify that the proposed SNSSLrL method has greater advantages than the latest state-of-the-art methods.


Introduction
Hyperspectral image (HSI) has played an important role in many modern scenarios like urban planning, agricultural exploration, criminal investigation, military surveillance, etc. [1]. However, the existing mixed noise generated during the process of digital imaging, e.g., Gaussian noise, salt and pepper, stripes and deadlines, seriously diminishes the accuracy of the above applications [2]. In addition, a series of HSI processing applications (e.g., unmixing [3,4], classification [5][6][7], super-resolution [8], target detection [9,10]) are greatly dependent on the imaging quality of HSI. Therefore, as a key preprocessing step of the above-mentioned applications, the studies on HSI mixed denoising methods have become a hot research topic and aroused widespread interest in the field. Now, many methods of HSI denoising have been proposed. As each band of HSI can be regarded as a two-dimensional matrix, the simplest way is to directly use the traditional and classic two-dimensional image denoising method, such as Total Variation (TV) [11], K-SVD [12], NLM [13], BM3D [14], WNNM [15]. Part of the methods are based on the domain transform theory, that is, data are transformed from the original spatial domain to a transform domain for denoising, and then the denoising results are inversely transformed to the spatial domain, for example, based on curvelet transform [16], wavelet transform [17], principal component analysis (PCA) [18]. However, the disadvantage of these methods is that they ignore the strong correlation between the spatial dimension and spectral dimension of HSI, which leads to a poor denoising effect [19]. Therefore, to make full use of spectral and spatial information, a series of methods are proposed. On the basis of BM3D [14], Maggioni et al. [20] proposed BM4D method by using nonlocal similarity and Wiener filtering. Compared with using TV [21] method in each band, the adaptive TV method is applied to the spatial-spectral dimension in [22,23], which improves the denoising effect. Similarly, Qian et al. [24] proposed a local spatial structure sparse representation method to denoise HSI in groups. However, although the above-mentioned traditional denoising methods have good performance on HSI denoising, they work insufficiently well when the noise intensity is large and the noise category is specific, such as deadlines and stripes.
Recently, deep learning techniques and convolutional neural networks (CNN) have made great achievements in visual tasks, e.g., image classification [25][26][27][28], object detection [29][30][31], image forensics [32][33][34][35] and semantic segmentation [36][37][38]. Moreover, CNN is also used for various image restoration tasks [39][40][41][42], especially for denoising HSI [43,44]. In these methods, strict encoding and decoding structures [43] are adopted to construct a deep convolutional neural network through end-to-end training mode. Due to the ability of CNN to learn the nonlinear mapping function from a large amount of training data [45], the methods based on CNN also obtain fine denoising results. CNN has been extended for HSI denoising in [44]. However, each spectral band was denoised separately, and the spectral correlation between HSI is not fully utilized. These approaches model HSI by learning the multichannel features with two-dimensional convolution, which loses the rich information of spectral dimension. Different from the previous deep learning method that always takes two-dimensional convolution as the basic operation of the network, Wei et al. [46] introduces a three-dimensional convolution and a quasi-recurrent pooling function, which makes full use of the rich spectral information in HSI, thus significantly improving the accuracy and efficiency of the model.
In recent years, low-rank matrix restoration (LRMR) has become an important method of HSI restoration. This method involves learning low dimensional subspace from high dimensional data [47]. The most representative is the HSI denoising method based on LRMR proposed by Zhang et al. [48]. Xu [49] introduces a low-rank regularizer, which uses matrix factorization and weighted nuclear norm to better approximate the original low-rank hypothesis. However, in this framework, singular value decomposition is needed to solve these weighted kernel norm optimization problems. Not only is the computational complexity of the algorithm relatively high, but it is also not robust to the removal of mixed noises [50]. To solve these problems, Cao et al. [51] adaptively fit the noise through the exponential power distribution of the mixed noise, and proposed a more flexible method to apply in HSI mixed denoising. However, HSI is three-dimensional data and the spectral relationship between different bands is ignored based on matrix decomposition. The algorithm based on low-rank tensor decomposition maintains the global structure of hyperspectral data. A series of HSI restoration methods based on tensor decomposition have been proposed. For example, using CANDECOMP/PARAFAC (CP) decomposition [52], Tucker decomposition [53] and tensor-SVD decomposition [54] to establish the HSI restoration model, the noise removal results have been greatly improved. On this basis, people use various regularization methods to explore its prior. For example, Wang et al. [55] introduced Total Variation (TV) regularization into the low-rank tensor restoration framework to preserve the spatial piecewise smoothness. Sun et al. [56] propose a non-local enabled tensor subspace low-rank learning strategy to achieve mixed denoising. Lin et al. [57] introduce a deep prior to the low-rank decomposition framework to explore the nonlinear property of noise. Zheng et al. [58] propose a tensor fiber rank model for HSI mixed noise removal to explore the high-order low-rankness of HSI.
Although the above-mentioned methods have achieved a certain effect on the problem of mixed denoising, directly processing the original HSI often leads to inefficient execution of the algorithm, which is extremely fatal in real applications. In this paper, a spectralsmoothness regularized subspace low-rank learning method with parameter-free BM4D denoiser is proposed for HSI mixed denoising, and the detailed flowchart of the proposed model is shown in Figure 1.
The main contributions of this paper are summarized as follows.

1.
The sampled HSI is first projected into a low-dimensional subspace, which greatly alleviates the complexity of the denoising algorithm. Then, spectral smoothing regularization is enforced to the basis matrix of the subspace, which constrains the reconstructed HSI to maintain spectral smoothness and continuity.

2.
The plug-and-play non-local BM4D denoiser is enforced to the coefficient matrix of the subspace to fully utilize the self-similarity of the spatial dimension of HSI. Furthermore, the L 1 norm is used to separate the sparse noise. In the process of alternate optimization, the latent clean HSI is gradually learned from the degraded HSI. The remainder of this paper is organized as follows. Section 2 formulates the preliminary problem of HSI mixed denoising and elaborates the proposed SNSSLrL model. In Sections 3 and 4, experimental evaluations compared to the latest denoising methods are conducted to demonstrate the performance of our proposed SNSSLrL method. Eventually, the conclusion and future work are summarized in Section 5.

Degradation Model
In the complex process of signal acquisition and digital imaging, HSI is often contaminated by Gaussian noise and sparse noise. Gaussian noise refers to noise whose probability density obeys Gaussian distribution, e.g., dark noise and read noise. Sparse noise refers to a small proportion of noise in the sampled HSI, e.g., stripes, salt & pepper, deadlines, etc.
According to the additive principle of noise, the denoising problem of an HSI with the size of n 1 × n 2 × n 3 can be formulated as: where Y ∈ R n 3 ×m is the matrix formed by the vectorization of the third-order tensor Y ∈ R n 1 ×n 2 ×n 3 , which is used to represent the degraded HSI (m = n 1 × n 2 ). L ∈ R n 3 ×m is a latent clean image, N ∈ R n 3 ×m is the Gaussian noise exist in HSI, and S ∈ R n 3 ×m is the sparse noise, i.e., salt and pepper, stripes, and deadlines. Given the sampled image Y, clean L can be restored by separating N and S. Evidently, in Equation (1), X cannot be solved directly from Y. Therefore, the priors of each component are required to compress the solution domain of the uncertain equation. Through the introduction of regularization constraints, the original additive model is transformed into the following regularization model: where R 1 (L) and R 2 (S) are the regularization terms of clean HSI and sparse noise. The idea of regularization refers to introducing certain constraints by using the prior information of the image and noise to compressing the solution space of the uncertainty problem. R 1 (L) describes the prior information of the latent noise-free image L, such as low-rank property. R 2 (S) describes the prior information of sparse noise to ensure that the observation image Y is decomposed into low-rank component L and sparse noise S without excessive deviation, such as the L 1 norm. The low-rank regularization and non-local regularization of the latent clean HSI L have been proven to achieve fine denoising results. Imposing the rank constraints on the potential clean image, low-rank regularization can separate the latent clean low-rank component from the full-rank noise component. As there are similar patches that exist in HSI, many denoising studies use non-local regularization to explore this self-similarity. For sparse noise S, L 1 norm is able to fully use its sparse property. γ is a regularization parameter which is employed to control the contribution of different regularization terms. Fidelity term Y − L − S 2 F ≤ ε is used to constrain the Gaussian noise, where ε is the standard deviation of noise.

Subspace Low-Rank Regularization
Even working favorably on the HSI denoising problem, the following disadvantages still exist in above regularization model:

•
If heavy Gaussian noise exists in the sampled image, it is difficult to obtain a fine denoising result by low-rank regularization. Meanwhile, since some sparse noises will also show a certain extent of low-rankness, low-rank regularization is also powerless in such a case. • For the non-local regularization, the denoising result depends heavily on the choice of the search window and neighborhood window, and the cost of fine denoising results is often exchanged for greater execution time. In real scenarios, too much calculation time is often undesirable.
To address the above problems, we have developed the subspace low-rank learning methods in our previous work [2,19,[59][60][61]. The principle of linear mixing in HSI instructs us that high-dimensional HSI exists in a subspace spanned by endmembers. Projecting the original high-dimensional image data into a low-dimensional subspace spanned by an orthogonal basis makes full use of the spectral low-rankness. Then imposing the non-local regularization constraints on the coefficient component of the subspace to realize the utility of the spatial self-similarity of HSI. Therefore, the subspace low-rank regularization model is formulated as arg min where E ∈ R n 3 ×k is the basis of subspace and E T E = I is the orthogonal constraint designed to drive E spans more of the original space. Z * is the low-rank regularization term imposed on coefficient Z ∈ R k×m of subspace instead of the original image. The function σ NL (Z) represents that a non-local transformation is performed on the coefficient Z to utilize the non-local self-similarity of the image. γ and λ is regularization parameters which are used to control the contribution of nuclear norm and sparse norm. The fidelity term Y − EZ − S 2 F ≤ ε indicates the subspace low-rank decomposition, i.e., L = EZ. In this way, the regularization operations imposed on the original high-dimensional image are transformed into the operations executed in the subspace. Not only enhances the exploration of low-rankness but also greatly reduces the computational complexity of the non-local algorithm.

Proposed Model
To further improve the performance of the subspace representation model, in this section, we proposed the spectral-smoothness regularized subspace low-rank learning model with non-local denoiser. Firstly, in our design, the original high-dimensional image data are projected into a low-rank subspace to obtain the basis matrix and coefficient matrix of the subspace. Then, spectral smoothing regularization is applied to the obtained basis matrix instead of orthogonal constraints. Meanwhile, the parameter-free non-local BM4D denoiser is embedded into the subspace model in a plug-and-play manner. Since the spectral continuous basis can promote the reconstruction of the continuous clean image data, and the non-local algorithm can remove the heavy Gaussian noise, the proposed model can achieve better restoration of the latent clean data. The proposed model is formulated as arg min where DE 2 F is the spectral smoothing term, D is the first-order difference operator along the spectral dimension. For σ NL (Z), the well-known plug-and-play parameter-free BM4D denoiser is employed. The idea of the BM4D denoising algorithm is to use the non-local self-similar property to achieve denoising. Because there are many similar non-local regions that exist in the natural image, BM4D extracts three-dimensional similar cubes from image tensor and stacks them to form a four-dimensional non-local tensor. Then the noise in the four-dimensional tensor is removed through the least square filtering [20]. In the process of solving, a spectral smoothing basis and a clean coefficient are learned. Our designed model not only inherited well the high computational efficiency of the subspace representation model but also further improved the accuracy of denoising. The embedding of parameter-free non-local BM4D denoiser also greatly reduces the parameter complexity of the model, which is very critical in practical applications.

Optimization
It is extremely difficult to solve Equation (4) directly. Here we design an optimization approach to solve the proposed SNSSLrL model. Through introducing a Lagrangian multiplier Λ, the original model is transformed into the following Lagrangian equation: where µ is the penalty parameter, Λ is a Lagrangian multiplier. The above Lagrangian equation can be solved by alternate iteration [62]. In each iteration, fixing the remaining variables and only updating the current variable, the optimization of Equation (5) can be transformed into the optimization of the following three subproblems. Through alternate updating, the basis and coefficient of the subspace are learned iteratively, thus realizing the reconstruction of potentially latent data. The main steps of solving Equation (5) are summarized in Algorithm 1.

Algorithm 1
Optimization procedure for solving Model 4.

Ensure:
The latent noise-free data X.

5:
Update S (t+1) by Equation (12). 6: Update Λ (t+1) by Equation (13). Update iteration by t = t + 1 8: end while The subproblem of updating Z is given by: which can be solved by the parameter-free plug-and-play BM4D denoiser: • Update E (Line 4): The subproblem of updating E is given by: which can be transformed into solving the following Sylvester equation: where D T D is a circulant matrix, Z (t+1) Z (t+1)T is a symmetric matrix. Therefore, by using first-order Fourier transform to diagonalize D T D and first-order SVD to diagonalize Z (t+1) Z (t+1)T , that is D T D = F T ΨF and Z (t+1) Z (t+1)T = U T ΣU, a closedform solution of Equation (7) is obtained by where • Update S (Line 5): The subproblem of updating S is given by: Here we use the soft-thresholding function softTH(A, τ) = sign(A) max(0, |A| − τ) to solve this subproblem efficiently: where µλ is the thresholding parameter. |A| is the absolute operation of matrix A. • Update Λ (Line 6): Updating Lagrangian multiplier Λ is given by:

Simulation Configurations
As shown in Figure 2, two datasets were used in the simulated numeric experiments, i.e., HYDICE Washington DC Mall (WDC), ROSIS Pavia University (PU) and HYDICE Urban. They have 191 and 103 high-quality spectral bands, respectively. Sub-images with the spatial size of 256 × 256 are cropped for simulation. The well-known HYDICE Urban dataset whose size is 307 × 307× 210 is used for the real experiment. Many bands of Urban dataset are severely degraded by watervapour absorption and mixed noise. To verify the denoising performance of the proposed algorithm, six representative state-ofthe-art HSI denoising methods are used as competitors, i.e., SLRL4D [61], 3DlogTNN [58], LRTDGS [63], LRTDTV [55], LLRSSTV [64], FSSLRL [60]. Among them, SLRL4D and FSSLRL are subspace-based methods, 3DlogTNN, LRTDTV, and LRTDGS are the latest tensor representation-based methods. All experiments are conducted on a workstation with dual Intel Silver 4210 and 128 G of RAM. To quantitatively assess the performance of denoisers, indexes such as peak-signal-tonoise-ratio (PSNR), structural similarity (SSIM), feature similarity index mersure (FSIM), Erreur relative global adimensionnelle de synthèse (ERGAS), and mean spectral angle (MSA) are used in the subsequent four cases of simulation.
• case 1: zero-mean Gaussian noise with the randomly selected variance in the range from 0.2 to 0.25 is first added to all bands. Meanwhile, in 20 continuous bands, 10% of pixels are contaminated by salt and pepper noise. • case 2: zero-mean Gaussian noise is added as the same condition in case 1, and deadlines with the randomly selected number in the range [3,10] and widths in [1,3] are added to the continuous 20 bands.
• case 3: zero-mean Gaussian noise is added as the same condition in case 1, and stripes with the randomly selected number in the range [2,8]  In order to reproduce the results of algorithms in this article for other interested readers, Table 1 shows the parameter settings of all the comparison methods on the simulated datasets. The parameters of the comparison methods are first set to the recommended parameter settings in the corresponding references and then those values associated to the optimal results for each dataset are adopted after fine-tuning.

Visual Evaluation Results
Among the four different cases, we chose representative case 1 and case 4 to evaluate the visual effects of seven different methods. Figures 3 and 4 show the recovery of the WDC data set. Figures 5 and 6 Figure 4b shows that severe mixed noise (impulse noise, stripes, and deadlines) destroys the high-frequency information of band 105 of the WDC under case 4, and the ground information is almost in an invisible state. It is clear that the FSSLRL method and the LLRSSTV method do not completely deal with the impulse noise, the image restored by the LRTDGS method is blurred, that is, a lot of detailed image information is lost. The 3DlogTNN method cannot complete the task of removing noise at edges. The denoising effect of the LRTDTV method and the SLRL4D method is superior to other methods, but compared with our SNSSLrL method, both in the smoothness of image restoration and detail preservation, there are obvious shortcomings. In Figure 3, we can get similar conclusions, where the 3DlogTNN method does not have the processing task of fringe noise in case 1, but its processing effect is still not ideal, mainly reflected in the loss of details. Although 3DlogTNN considers the third-order tensor of HSI, that is, considering both spatial and spectral information, it is only designed to remove Gaussian noise. Therefore, it appears powerless for removing the complex mixed noise in HSI. Due to the lack of spatial constraints, the FSSLRL achieves the lowest PSNR value. Unsurprisingly, among all the comparison methods, the proposed SNSSLrL denoiser achieves the best denoising effect. Figure 3 shows the PSNR and SSIM values of all methods on band 28 of the WDC dataset. Figure 5 shows that in the PU dataset, without stripes or bad lines, the recovery effect of the FSSLRL, LRTDGS, and LLRSSTV methods is obviously not as good as the other four methods. Meanwhile, 3DLogTNN, LRTDGS, and SLRL4D are obviously not as good as our SNSSLrL method. Compared with the second-placed SLRL4D method, the proposed SNSSLrL method achieves 1.16 dB and 0.0768 higher on PSNR and SSIM, respectively. The performance of each algorithm in the PU dataset in removing mixed noise is similar to their performance in the WDC dataset. It is clear that Figure 6b shows that band 82 of the PU dataset has been added with strong impulse noise, stripes, and deadlines. Both 3DLogTNN and LLRSSTV methods failed to effectively remove the stripes and deadlines. Although FSSLRL, LRTDGS, and LRTDTV can remove all the stripes and deadlines, they have poor performance in maintaining the details of the image. The PSNR values obtained by them are 26.51, 29.07, and 29.83 dB, respectively. The SLRL4D method and the proposed SNSSLrL method have achieved better mixed noise removal effects, but SNSSLrL is better in maintaining image details, especially the edge details of geometric shapes in the red box.
To further compare the performance of each algorithm, Figure 7 shows the spectral curve reconstructed by each algorithm at pixel (100,100) in case 4 of the PU dataset, where the red curve represents the true spectrum, and the blue curve represents the reconstructed one. It is clear that the LRTDTV method fails to effectively reconstruct the corresponding spectrum. The spectrum reconstructed by 3DLogTNN contains a lot of violent fluctuations, indicating that its performance is unstable on more bands, and there are still some noises that have not been removed. For LRTDGS, LLRSSTV, and FSSLRL methods, their reconstruction accuracy on the 1-10 band and 80-100 band is much lower than the proposed SNSSLrL method. For the SLRL4D method, its reconstruction accuracy in the 60-70 band is lower than the proposed method. Therefore, benefiting from the spectral smoothing constraint on the subspace basis, the proposed SNSSLrL method achieves the closest reconstruction to the true spectrum.

Quantitative Evaluation Results
To quantitatively evaluate the performance of each algorithm and give the reader a quantitative reference, Figure 8 shows the PSNR values as a function of bands for each algorithm under the four cases of WDC and PU datasets. In the first row of Figure 8, it is evident that the proposed SNSSLrL method is relatively stable in removing Gaussian noise and impulse noise. Its performance on the band 1-60 and 120-188 exceeds all comparison algorithms, but on the 60-120 band, the performance is relatively ordinary. Meanwhile, for mixed noise, especially stripes and deadlines (as shown in Figure 8b,d), although the proposed method does not achieve the best results on all bands, its stability and average performance is the best. The reason is that the subspace decomposition can effectively suppress the influence of noise, furthermore, the spectral smoothing constraint on the subspace basis enables the reconstructed image to maintain a certain continuity between the upper and lower bands, thereby further maintaining image details. In the second row of Figure 8, it is clear that the proposed method has achieved the best PSNR values on almost all bands, which shows that for PU data, whether it is Gaussian noise, impulse noise, or mixed noises, the proposed SNSSLrL method can achieve a better denoising performance than the other six competitors.  Table 2 lists the quantitative denoising results of the above seven algorithms in four cases on the WDC and PU datasets, with a total of five indicators. For case 1 of WDC, the PSNR value of SNSSLrL is 0.03 dB higher than the second-place 3DLogTNN method, while the ERGAS value is 0.69 higher (the lower the ERGAS value, the better the algorithm performance). For the other cases of the two datasets, the proposed SNSSLrL method has achieved the best quantitative results on the five indicators. For example, for case 4 of PU data, the proposed SNSSLrL method achieves a PSNR value of 0.78 dB higher than the second place LRTDTV method.
In conclusion, whether it is simple Gaussian noise and impulse noise or complex mixed noise, the proposed SNSSLrL method achieves the best results in all terms of visual effects and quantitative evaluation of HSI recovery compared to all six competitors.  Figure 9a) contains severe Gaussian noise, impulse noise, and vertical stripe noise; band 144 (as shown in Figure 10a) contains severe Gaussian noise, impulse noise, and horizontal stripe noise; band 210 (as shown in Figure 11a) is one of the most severely contaminated bands by noise, and it is almost impossible to see that it contains any useful information. From the reconstruction results in the above figure, it can be seen that the 3DLogTNN, LRTDGS, and LLRSSTV methods failed to restore the 105, 144, and 210 bands and most of the information of the reconstructed image has been lost. Thanks to the subspace decomposition framework and the spatial-spectral total variation constraint, the SLRL4D, FSSLRL, and LRTDTV methods can recover part of the image information, but most of the details are lost. For the above three severely polluted bands, only the proposed method reconstructed satisfactory results, which contain many fine details in the image.  Figure 12 plots the horizontal mean profiles of band 210 of the urban dataset. The horizontal mean profile can effectively reflect the mixed noise in the image, especially the impulse noise and band noise. As shown in Figure 12a, due to the pollution of heavy mixed noise in band 210, the horizontal mean profile of the original band 210 fluctuates rapidly. From Figure 12, the following conclusions can be drawn: The horizontal mean profiles of the LRTDGS, LLRSSTV, and FSSLRL methods contain higher-frequency fluctuations, indicating that the serious stripe noise has not been removed. Compared with SLRL4D, LRTDTV, and the proposed SNSSLrL method, the horizontal mean profile of 3DLogTNN is significantly different. The reason may be that due to the serious mixed noise pollution on the continuous bands, the 3DlogTNN method cannot effectively use the inter-band information to reconstruct the image. Compared with these methods, our proposed SNSSLrL achieves a flatter curve, which indicates that the cross stripes and various other noises of the urban data set are more effectively eliminated.

Experimental Results on Real Datasets
In summary, the competitive recovery results of our proposed SNSSLrL denoiser have been verified.

Quantitative Evaluation Results
To further verify the effectiveness of the SNSSLrL in real scenarios, Q-metric was used to quantitatively evaluate the denoising effect of each algorithm on urban data. Q-metric is a blind image spatial information index [65]. The greater value of Q-metric represents a better preservation of spatial details in the image. Table 3 lists the average Q-metrics of all bands in the urban dataset, the optimal is bolded. It is evident that the proposed SNSSLrL denoiser obtains the highest Q value, which indicates that it delivers the best restoration results in preserving the spatial fine details compared to the other six denoisers.
The running time of an algorithm is one of the important indicators used to measure the possibility of the real application of the algorithm. Table 4 lists the execution time of all competing methods on the real urban dataset, the optimal is bolded. The results in Table 4 are the average of ten running times, and all algorithms are performed in the same experimental environment. As shown in the table, the subspace decomposition techniques, such as FSSLRL, SLRL4D, and the proposed SNSSLrL method, have outstanding performance in computational efficiency and take less time than the other four methods. Thanks to the spectral smoothing constraint on subspace basis, SNSSLrL achieves better denoising effects with fewer iterations, which further saves running time.

Discussion
In the proposed SNSSLrL model, there are four important parameters, that is, dimension of subspace p, regularization parameters γ and λ, and the Lagrangian parameter µ. Since BM4D denoiser is a parameterless method, we do not need to make strict parameter adjustments for the non-local item in the model.

The Regularization Parameters µ, γ and λ
As shown in Formula (5), µ is a Lagrangian parameter. In theory, when 1 2µ tends to positive infinity, Formula (5) has the same solution as model (4). However, in actual applications, we usually set µ = 10/ max(m, n 3 )( √ m + √ n 3 ) to get a good result.
The regularization parameters γ and λ indeed have a great impact on the experimental results. Figure 13 shows the sensitivity analysis of the parameters γ and λ in case 4 of the simulated WDC dataset. It can be seen that when γ is located in the interval of [15][16][17][18][19][20][21][22][23][24][25] and λ is located in the interval of [0.1-0.2], the proposed SNSSLrL algorithm can achieve a relatively stable and satisfactory result on PSNR, SSIM, and ERGAS indicators, respectively.

The Dimension of Subspace p
The parameter p is an important parameter related to the subspace dimension. According to the linear mixing theory, its ideal value should be equal to the number of end-members existing in the scene. However, in actual scenarios, due to the complex nonlinear mixture, the ideal value of parameter p is difficult to guarantee. Therefore, in our experiments, the value of p is fixed to 4. we suggest choosing the optimal value of p in the range of [3][4][5][6][7][8]. Figure 14 plots the MPSNR and MSSIM values as a function of iterations under case 2 of the simulated PU dataset. It is clear that both the MPSNR value and the MSSIM value quickly enter a stable state after the first few iterations, which demonstrates that the proposed SNSSLrL algorithm has good convergence.

Conclusions and Future Work
In the framework of subspace decomposition, in the paper, a mixed noise removal algorithm is proposed for HSI that uses both spectral smoothing prior and non-local self-similar prior. The proposed SNSSLrL model was solved quickly and efficiently by the ADMM algorithm. Experiments on two simulated datasets and a real dataset show that the proposed spectral smoothing constraint on the subspace basis can effectively characterize the intrinsic properties of the subspace, and further improve the effect of mixed denoising. All indicators surpass several state-of-the-art HSI denoising algorithms with high time efficiency.
Future work will focus on utilizing advanced technology to exploit non-local similarities under the framework of subspace decomposition, and perform subspace decomposition in the sense of tensors, to achieve more accurate spatial-spectral information mining, thereby further improving the efficiency and accuracy of mixed noise removal.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: