SLRL4D: Joint Restoration of Subspace Low-Rank Learning and Non-Local 4-D Transform Filtering for Hyperspectral Image

: During the process of signal sampling and digital imaging, hyperspectral images (HSI) inevitably suffer from the contamination of mixed noises. The ﬁdelity and efﬁciency of subsequent applications are considerably reduced along with this degradation. Recently, as a formidable implement for image processing, low-rank regularization has been widely extended to the restoration of HSI. Meanwhile, further exploration of the non-local self-similarity of low-rank images are proven useful in exploiting the spatial redundancy of HSI. Better preservation of spatial-spectral features is achieved under both low-rank and non-local regularizations. However, existing methods generally regularize the original space of HSI, the exploration of the intrinsic properties in subspace, which leads to better denoising performance, is relatively rare. To address these challenges, a joint method of subspace low-rank learning and non-local 4-d transform ﬁltering, named SLRL4D, is put forward for HSI restoration. Technically, the original HSI is projected into a low-dimensional subspace. Then, both spectral and spatial correlations are explored simultaneously by imposing low-rank learning and non-local 4-d transform ﬁltering on the subspace. The alternating direction method of multipliers-based algorithm is designed to solve the formulated convex signal-noise isolation problem. Finally, experiments on multiple datasets are conducted to illustrate the accuracy and efﬁciency of SLRL4D.


Introduction
It began with the forging of the great imaging spectrometer. Many merits were given to the generated hyperspectral image (HSI)-multi bands, high resolution and rich details [1]. The vision of mankind was no longer merely limited to the spatial dimension with the assistance of HSI. An increasing number of application fields now rely on HSI, for instance, urban planning, geological exploration, military surveillance and precision agriculture [2]. Among these application fields, many critical tasks of computer vision processing such as unmixing [3], target detection [4], Inspired by the robust principal component analysis model (RPCA) [20], He et al. conducted research to explore the low-rank property of HSI and modelled the HSI restoration as a low-rank matrix recovery (LRMR) problem, thus realizing the removal of various mixed noises [21]. This work proved that the probability is extremely high of restoring the low-rank latent clean HSI matrix from the observation HSI matrix, which is degraded by mixed noises. During the process of restoration, performance would be much improved if the global and local redundancies of HSI were both thoroughly utilized [22]. In Reference [23], a sparse method is employed to encode the global and local redundancy of the spectral dimension. Meanwhile, the global redundancy of the spectral dimension is explored by a low-rank constraint. Chen et al. [24] proposed a low-rank non-negative matrix factorization method to capture the spectral geometric features while maximizing the restoration of spatial detail. However, the spatial smoothing information of HSI is not fully utilized by the studies above. Due to the powerful ability to maintain the spatial smoothness of the image, the total variation (TV) regularized low-rank methods [25] are widely used in HSI restoration research. Aggarwal et al. [26] extended the traditional TV model to the three-dimensional space, and proposed a spatio-spectral total variation (SSTV) method for HSI restoration. The proposed SSTV method considers the high correlations in both the spatial and spectral dimensions, and performance enhancement of the TV-based restoration is achieved. However, in SSTV, only local spatial and spectral features have been explored. There is still no further research on the global low-rank property of HSI. Based on LRMR, Wang et al. continued their previous research [27] and proposed a combinatorial method with low-rank constraints named LSSTV [28]. The HSI restoration problem is modelled as a weighted nuclear norm minimization problem, and the jagged distortion left in SSTV is solved under certain circumstances. As an improvement of the LRMR, He et al. designed a TV regularized low-rank factorization method [29]. For further enhancement, they divided the HSI cube into many local patches, then imposed low-rank constraints on these patches [30]. This is about images in local fields have a high probability of containing the same features, thus, they usually share similar spectral features [31]. Tensor representation methods also shows great potential to explore the low-rank property [32]. Fan et al. extended the low-rank constraints to tensor represented HSI, and proposed a low-rank tensor recovery method that maintains correlations between three dimensions of HSI [33]. Wang et al. [34] combined tensor Tucker decomposition and anisotropic spatial-spectral TV, achieving exciting results in the removal of mixed noises.
However, the lack of spatial information often leads to poor restoration performance when strong Gaussian noise arrives. Most of the low-rank methods above only consider the piecewise smoothness or local similarity of the spatial dimension. Meanwhile, TV-based methods often cause distortion in the spectral domain, and some noises also appear smooth property (e.g., white Gaussian noise) [35]. The structural information of restored HSI is not guaranteed in these methods.

Restoration by Non-Local Similarity
Non-local self-similarity (NSS) theory holds that a patch of an image has many homogeneous patches with similar characteristics at other positions of the image. Since the structure of the clean image is redundant, similar patches of the target patch may exist at any position [36]. The NSS-based restoration method exploits this information redundancy of HSI to get rid of noises, where details of the spatial domain are kept to the greatest extent in course of denoising [37]. The priors of low-rankness and NSS are bonded in Reference [38]. By imposing non-local priors to the LRMR model, both pixels and spectrum in HSI are double corrected to remove mixed noises while preserving the fine spatial structure of HSI, thus improving the denoising performance. Bai et al. [39] applied a fixed window to extract non-local patches from the original image, and characterized the non-local similarity of HSI in the spatial domain and the global similarity in the spectral domain. In Reference [40], Xue et al. embedded the combined prior terms of low-rank and NSS into a sparse representation model. For further research, the tensor CANDECOMP/PARAFAC (CP) decomposition is applied to the joint denoising model [41]. Under this circumstance, the spatial information, for instance, image edges and textures of HSI, was further exploited, and the low-rank property was further utilized by the clustering of non-local patches. Both the global correlation along the spectrum and NSS were explored by Xie et al. [42]. In this work, all non-local patches were stacked to a third-order tensor and achieved excellent denoising performance. Similar research has been conducted in Reference [43]; the group sparsity and low-rank property of HSI are further utilized. Zhang et al. hypothesized that the piecewise smooth feature exists in the restored non-local patches of HSI, and thus conducted research on the combination of NSS and TV [44]. Considering the poor maintenance of CP decomposition on the intrinsic structure of HSI, Chen et al. applied a novel tensor decomposition approach named tensor ring decomposition to the non-local based HSI restoration model [45].
Unfortunately, although NSS-based low-rank methods have delivered gratifying restoration performance, the setting of certain parameters impacts the result of restoration excessively. For instance, various sizes of the search window and number of non-local patches often bring different results, some of which are unacceptable. Maggioni et al. designed a block-matching and 4D filtering (BM4D) restoration algorithm, which extended planar space denoising to three-dimensional space, and provided new ideas for HSI restoration [46,47]. This work simultaneously exploits the local correlation and non-local correlation between voxels instead of pixels, high-dimensional data is thus divided into many patches with similar features. Then, these patches are stacked together to form a four-dimensional hyperrectangle. Eventually, two diverse steps of estimation occur, the signal and noise are isolated by order of thresholding or Wiener filtering. The best legacy of this work is that actual parameter-free denoising is accomplished when applying BM4D to the HSI restoration problem. However, as with all NSS based restoration methods, the processing procedure of BM4D is relatively intricate, which ineluctably results in excessive time consumption when confronting the huge size of the original image space of HSI.

Restoration by Subspace Projection
Recently, the subspace low-rank (SLR) representation method has shown great potential to address the above challenges [48][49][50][51][52]. The theory of manifold learning suggests that the data in high-dimensional space is immoderately redundant, and principal information lives in the low-dimensional subspace [53]. Similar theories exist in the HSI compression field. Since the high-dimensional data is redundant, according to the low-rank property of HSI, the dimensions of the dictionary matrix and coefficient matrix generated from the low-rank decomposition are necessarily much lower than the dimension of the original image, thereby the purpose of data compression is realized. The reconstruction process is the inverse of compression, that is, the linear combination of dictionary and coefficient is the original image. Therefore, through this low-rank decomposition method, it is no longer necessary to store the original high-dimensional data, but only the compressed low-dimensional data. For low-rank image restoration, Figure 1 illustrates the distribution of singular value in the original space and low-dimensional subspace of the HYDICE Urban dataset. Since the spectrum of HSI lives in the low-dimensional subspace, as shown in Figure 1, the subspace has a lower singular value than the original space, which indicates that the subspace decomposition can yield a more accurate low-rank approximation. Part of this subspace low-rank property is given by the lower dimension of the projected matrix, and the other is given by the finiteness of the linear combination of end-members. For any pixel in HSI, it is impossible to be mixed by all end-members in the scene, but only a few of them (usually 2-3 end-members). This sparsity of the linear combination of end-members leads to the low-rankness of the decomposed coefficient matrix. Therefore, through imposing low-rank regularization to the subspace of HSI instead of the original space, better restoration results would be achieved. Reference [49] fused spectral subspace low-rank learning and superpixel segment, its restoration results proved that subspace does carry stronger low-rank property. By taking advantage of the subspace decomposition, a fast restoration method through using denoiser on subspace is proposed in [50] , but the fixed subspace basis they used also prevents the further rising of restoration results. References [48,51] proposed a Gaussian noise-removal method by combining low-rank subspace decomposition and weighted nuclear norm minimum, and achieved outstanding results. However, this advancement comes with the reduction of usability, as complex parameter settings greatly suppress the real application. Reference [52] modeled mixed denoising problem with a robust l 1 norm, and achieved the unification of complex statistical noise distribution. The number of parameters that need to be manually tuned is reduced, thus its usability and efficiency in real application is further enhanced. However, these studies also show that along with low-rank-based restoration methods, SLR-based methods are also powerless when encountering structured sparse noise and strong Gaussian noise. It is essential to further exploit the high correlation in the spatial dimension of HSI by exploring the NSS prior based on SLR. Meanwhile, the complexity of the parameters should be as simple as possible to make the proposed method easier to be implement. Through the discreet design of a low-dimensional orthogonal dictionary and prior regularization, the joint method of SLR and NSS not only accomplishes precise restoration of HSI, but also reduces the cost of computation. In this way, the NSS property is utilized in subspace, the dilemma of heavy computation cost of NSS-based methods is distinctly compensated. In this paper, a joint method of subspace low-rank learning and non-local 4-D transform filtering, named SLRL4D, is presented for HSI restoration. By assuming that the latent clean HSI exists in the low-dimensional subspace, the NSS prior of the HSI is embedded into the SLR learning architecture. Thus, the high correlation of HSI in both spatial and spectral dimensions is thoroughly explored, and the successful removal of the mixed noises of HSI is realized. The main contributions of this paper are summarized as: 1. The subspace low-rank learning method is designed to restore the clean HSI from the contaminated HSI which is degraded by mixed noises. An orthogonal dictionary with far lower dimensions is learned during the process of alternating updates, thus leading to precise restoration of the clean principal signal of HSI. 2. Based on the full exploration of low-rank property in subspace, BM4D filtering is employed to further explore the NSS property within the subspace of HSI rather than the original space. By preserving the details of the spatial domain more precisely, the complete parameter-free BM4D filtering leads to easier application in real world. 3. Each term in the proposed restoration model is convex after careful design, thus, the convex optimization algorithm based on the alternating direction method of multipliers (ADMM) could be derived to solve and optimize this model. Meanwhile, due to the HSI being decomposed into two sub-matrices of lower dimensions, computational consumption of the proposed SLRL4D algorithm is much lower than other existing restoration methods. 4. Extensive experiments of quantitative analysis and visual evaluation are conducted on several simulated and real datasets, and it is demonstrated that our proposed method not only achieves better restoration results, but also improves reliability and efficiency compared with the latest HSI restoration methods.
The remainder of this paper is organized as follows. Section 2 formulates the preliminary problem of HSI restoration. Then, the combined restoration approach of subspace low-rank learning and nonlocal 4-d transform filtering is elaborated specifically in Section 3. In Sections 4-6, substantial experimental evaluations are conducted to demonstrate that our proposed method is both effective and efficient. Eventually, the conclusion and future work are drawn in Section 7.

Degradation Model
During the process of acquisition and transmission, a succession of noises (e.g., Gaussian, deadline, impulse and stripe noises) frequently degrade the imaging quality of HSI. Denote the observed HSI cube as a 3D tensor Y ∈ R m×n×l with the size of spatial dimension m × n and the number of spectral bands l.
Considering a scenario of mixed noises contamination, the degradation model of HSI restoration is formulated as: where Y ∈ R l×d is the 2D HSI matrix reconstructed from the 3D tensor Y with d pixels and l spectral bands(d = m × n). X ∈ R l×d denotes the latent clean HSI which demands to be recovered, N ∈ R l×d and S ∈ R l×d denote the Gaussian additive noise and sparse noise, respectively. The purpose of HSI denoising is to seek the restoration of clean image X from the observed image Y by isolating S and N. However, it is clear that Equation (1) is an unsolvable ill-conditioned problem. Employing prior knowledge to introduce regularization constraints, the compression of the uncertain solution domain shall be realized: where R 1 (X) and R 2 (S) are the regularization terms that constrain the pure image and sparse noise, respectively. The regularization parameter λ is accountable for the balance of the two regularization terms. Y − X − S 2 F ≤ ε is the fidelity term, which constrains the Gaussian additive noise, and ε is the standard deviation of Gaussian noise.

Low-Rank Regularization
One of the most efficient ways for R 1 (X) is low-rank regularization. This is mainly caused by the high linear correlation of the spectral domain between diverse areas of HSI. The low-rank property of HSI can be effectively explored by enforcing nuclear norm to the original image X. As for the R 2 (S), l 0 -norm is an option that meets the constraint of sparsity, thereby the restoration problem of HSI is reformulated as: In the process of solving, l 1 -norm is employed as the convex relaxation of l 0 -norm, thus, Equation (3) is transformed into a solvable problem.

Non-Local Regularization
Furthermore, the NSS is also a highly effective regularization constraint while exploring the low-rank property of images: Non-local regularization X NL realizes the full advantage of redundant information in the latent clean HSI and achieves the maintenance of image details in the course of restoration. By setting a neighborhood window and a search window, the estimate weight value of the target area is calculated from the area with a similar neighborhood structure.
However, even working favorable for HSI restoration, there still exists the following disadvantages: • For the low-rank regularization, when stripes or deadlines emerge at the same location in HSI, they also appear to have structured low-rank property, hence it is difficult to isolate sparse noise from low-rank signals. Meanwhile, its denoising ability is highly degraded when encountering heavy Gaussian noise.

•
For the non-local regularization, the performance of restoration relies on the selection of search window and the neighborhood window, better performance is often achieved with much higher computational complexity. In practical applications, the excessively high time consumption of processing is fatal.

Subspace Low-Rank Representation
According to the basic principle of manifold learning, high-dimensional data lives in a low-dimensional manifold. Meanwhile, both spatial and spectral domains of HSI are correlated greatly. Therefore, under the linear mixture model, the HSI matrix X can be decomposed into the product of an endmember matrix and an abundance matrix. In specific HSI scenes, each pixel is a mixture of a few pure endmembers, and the number of endmembers is far less than the dimension of original HSI space. Projecting the original image signal into a low-rank subspace, the degraded HSI is revealed as where D ∈ R l×p is the dictionary matrix that spans each subspace of original space X, and A ∈ R p×d is the representation coefficient matrix of X in respect of D. In this equation, p is the dimension of the subspace. Consequently, the subspace low-rank representation is denoted as arg min In this enhanced low-rank restoration model, the nuclear norm is employed to constraint the coefficient matrix A. It has been proven that under the certain dictionary D, integrated HSI could be recovered accurately by the SLR representation model [49]. It inspired that the employment of latent clean HSI as a dictionary might be a good choice. Therefore, the SLR representation model is converted to arg min An additional fidelity item was implemented into Equation (6) by supposing that subspace D equal to clean HSI X. Although certain effects are achieved, the following disadvantage of this strategy still exists.

•
The use of clean HSI X ∈ R l×d as the dictionary D ∈ R l×p will inevitably cause the situation l = p.
Hence the dimension of the coefficients matrix A will increase to l × l immediately. The velocity of calculation and the rate of convergence will be decreased remarkably.

•
In the face of heavy Gaussian noise, this model is not able to perform satisfied restoration results. The main reason is the plentiful sharp feature of the spatial domain in the original HSI are directly ignored.

Proposed SLRL4D Model
To address these challenges, we propose the SLRL4D method as follows. The flowchart of the proposed SLRL4D method is depicted in Figure 2. In our model, an orthogonal dictionary, D ∈ R l×p with p l is learned during the process of HSI restoration. In this manner, the computational cost is much decreased due to the reduction of subspace dimensions. Meanwhile, aiming to achieve the finest preservation of spatial details, the well-known non-local 4-D transform filtering method BM4D is invited into the subspace low-rank learning model. Through superimposing cubes of voxels with similar features in A, the composited hyperrectangle accomplishes successful exploration of the NSS. Under the above regularization, the proposed restoration method SLRL4D is designed as arg min where φ NL (A) is the non-local regularization term which is imposed on the representation coefficient matrix A. The corresponding convex optimization subproblem could be efficiently solved by the BM4D filtering without any parameter. In our design, the following advantages emerge: • The low-rank property of the spectral dimension lives in multiple subspaces and NSS of the spatial dimension is synchronously exploited in the process of restoration. Better restoration results would be achieved consequently.

•
The dimension of subspace p is far lesser than the dimension of original space l which brings far lesser resource consumption and computational time in the practical experimentation. Higher restoration performance would be achieved consequently.

•
Each term of this model is convex, thereby it is solvable with the employment of distributed ADMM algorithm. Easier restoration solution and optimization would be achieved consequently.

Model Optimization
Obviously, attempting to solving Equation (8) directly is extremely difficult. Through introducing two variables L 1 = L 2 = A, the proposed SLRL4D model is translated as the following equivalent form where Γ 1 , Γ 2 and Γ 3 are three augmented Lagrange multipliers, µ is the penalty parameter, and ζ {I} D T D is the indicator function which is transformed from the orthogonal constraint of dictionary D. In this way, ADMM is prone to solve our proposed model by minimizing the above Lagrangian function [54]. The original difficult problem could be solved through alternately optimizing the divided five easily solved subproblems. The main steps of solving Equation (9) are summarized in Algorithm 1. Mixed noises and latent clean HSI are gradually isolated during the process of alternating updates the following five variables: • Update S (Line 4): The subproblem of updating S is given by: The soft-thresholding algorithm is employed to efficiently solve this subproblem: where softTH τ (M) = sign (M) max (0, |M| − τ) is the soft-thresholding operator with threshold value τ. |M| is the absolute value of each element in matrix M.

•
Update D (Line 6): The subproblem of updating D is given by: According to the reduced-rank Procrustes rotation algorithm [49], the solution of this subproblem is formulated as T , Σ and Σ are the left singular vector and right singular vector ofΣ. In this way, the solution to D is given by the D = UV T form.

13:
Update lagrange multipliers Γ i (k+1) , i = 1, 2, 3 14: The subproblem of updating L 1 is given by: Under our design, BM4D filtering could efficiently solve this non-local regularized subproblem at lightning speed: • Update L 2 (Line 10): The subproblem of updating L 2 is given by: Inviting the well-known singular value shrinkage method [56] to solve this nuclear norm subproblem, thus the solution could be formulated as: where D τ (·) is the singular value thresholding operator, which is defined as: The subproblem of updating A is given by: Through the calculation of partial derivative, the solution of this convex subproblem is equivalent to the beneath linear equation: where 3 , and I is the identity matrix, D (k+1)T is the transpose of the matrix D (k+1) .

Experimental Configurations
Aiming to verify the effectiveness and efficiency of our proposed SLRL4D restoration method, extensive experiments were conducted in multiple simulated and real HSI datasets. Main configurations of the experiment environment are Intel Core I7-9700K CPU @4.90 GHz, RTX2060Super GPU with 8GB graphics memory, and 16GB Dual channel RAM. Before the experiment, the gray values of each band were first normalized into the interval [0, 1]. The source code of all the comparison methods is downloaded from the author's homepage or obtained directly from the author. All our experiments were run on MATLAB 2018b and Windows 10. Source code of the proposed SLRL4D, corresponding datasets, and parameter settings are available in https://github.com/cxunhey/SLRL4D.

Datasets
Four different datasets were employed in our experiments, including two simulated datasets and two real datasets. The details of these datasets are listed as follows.

Comparison Methods
Nine different methods were employed in our experiment, and all parameters are manually optimized according to the instructions in the corresponding references. The details of these methods are listed as follows.
• NLR-CPTD [41]: One of the latest HSI restoration methods which combine the low-rank CP tensor decomposition and non-local patches segmentation. It achieves state-of-the-art performance in removing Gaussian-stripes mixed noises removal. • LRTDGS [57]: One of the latest HSI restoration methods that combine the low-rank Tucker tensor decomposition and group sparsity regularization. It achieves state-of-the-art performance in isolating structured sparse noise. • LRTDTV [34]: A combination method of the tensor tucker decomposition and anisotropic TV, the piecewise smooth property of spatial-spectral domains and the global correlation of all bands are successfully explored. • LLRSSTV [30]: Local overlapping patches of HSI are first cropped, then SSTV is embedded in the LRMR framework, therefore, the isolation of clean HSI patch and mixed noises is achieved. • FSSLRL [49]: Superpixel segmentation is creatively introduced into the subspace low-rank learning architecture, and both spatial and spectral correlations are exploited simultaneously. • TDL [58]: A dictionary learning denoising method based on the well-known Tucker tensor decomposition technique, outstanding Gaussian noise removal is achieved with the splitting rate of convergence. • PCABM4D [47]: An enhanced principal component analysis method based on BM4D, BM4D filtering is employed to eliminate the remaining low-energy noisy component after PCA is executed in HSI. • LRMR [21]: A patch-based restoration method under the framework of low-rank matrix recovery, the low-rank property of HSI is explored by dividing the original HSI into several local patches. • BM4D [46]: An improved version of the well-known BM3D denoising algorithm, NSS of HSI is efficiently utilized in voxels instead of pixels.

Simulation Configurations
For simulating the contamination scenario of extreme noise to validate the restoration performance of our proposed method, six different instances of noises were added to the WDC and PU datasets. • Case 4: Deadlines were added to 20 continuous bands, and their widths and number were randomly selected, ranging from [1,3] and [3,10], respectively. The variances and added bands of Gaussian noise were selected in the same manner as in Case 2 still. • Case 5: The stripes were added to 20 continuous bands, whose numbers were also randomly selected, ranging from [3,10]. And the added manner of Gaussian noise were also selected same as in Case 2. • Case 6: Mixed noises composed of multiple diverse types were added to different bands of datasets, that is, Gaussian noise, impulse noise, deadlines and stripes were added in the same manner as Case 2-5 simultaneously.

Visual Evaluation Results
For the visual performance of different methods, representative Case 1 and Case 6 are presented for evaluation. Figures 3 and 4 illustrate the restoration results of WDC dataset. To further evaluate the visual results, zoom-in figures of specific area are enlarged in the bottom right corner.  As shown in Figure 3, the original band 1 of WDC is highly degraded by heavy Gaussian noise, ground information is completely unobservable. BM4D fails to remove noise and distorts the spatial domain of the image. LRMR, PCABM4D and TDL partially remove noise but are still unacceptable. FSSLRL, LLRSSTV, LRTDTV and LRTDGS are powerless in the front of heavy Gaussian noise. The non-local regularized method NLR-CPTD achieved fine removal of Gaussian noise, but there are still some gaps compared to SLRL4D. As the enlarged box shows, SLRL4D performs best in preserving details while the accurate restoration is realized successfully.
Same results are shown in Figure 4. Figure 4b indicates that band 110 of WDC in Case 6 is severely contaminated from mixed noises. Both BM4D, LRMR and PCABM4D failed their mission. TDL failed to completely remove the stripes. Four low-rank based methods, FSSLRL, LLRSSTV, LRTDTV and LRTDGS achieved fine removal of mixed noises, however, our proposed SLRL4D not only achieved the same outstanding mixed noises removal but also achieved the finest detail preservation according to the enlarged box. In addition, Figure 4k shows that the ability of NLR-CPTD in removing stripes partly decreases while heavy Gaussian noise exists.
Restoration results of all methods on PU dataset are shown in Figures 5 and 6. Band 66 of the PU dataset is also degraded by heavy Gaussian noise. Among all methods, our SLRL4D method undoubtedly accomplished the best spatial information restoration as shown in the enlarged boxes of Figure 5. Severe mixed noises corrupted the band 82 of PU in Case 6, as shown in Figure 6b, and no progress was made by BM4D, LRMR, PCABM4D and TDL. Only FSSLRL, LLRSSTV and LRTDGS achieved the removal of mixed noises, but none of them ever transcended SLRL4D in spatial details preservation. In summary, either with Gaussian noise or mixed noises, our SLRL4D method achieved the best results in visual performance for HSI restoration compared with all nine competitors. The superiority of low-rank learning and non-local filtering in the subspace of HSI is revealed.

Quantitative Evaluation Results
Five evaluation metrics with three positive metrics that is, the peak signal-to-noise ratio (PSNR), structure similarity (SSIM), feature similarity (FSIM) and two negative metrics, namely, erreur relative globale adimensionnellede synthese (ERGAS), and mean spectral angle (MSA) are invited to conduct a quantitative evaluation of our proposed method. The higher values of positive metrics indicate better restoration performance, as well as the lower values of negative metrics. Figure 7a-l shows the curves of the above three positive metrics comparison of all restoration methods on simulated WDC and PU datasets in Case 1 and Case 6. For Gaussian noise degraded Case 1 on WDC, SLRL4D achieved the highest PSNR, SSIM and FSIM in most of bands, and overwhelming dominance is achieved in more than half of all bands. Similar results occur in Case 6 on WDC which is heavily contaminated by mixed noises. SLRL4D still achieved optimal results in most of bands. Although the results in few bands are not optimal yet, the average results of all bands still exceed all competing methods. In Case 1 of PU dataset, our SLRL4D method purchased optimal PSNR and SSIM value effortlessly. Meanwhile, the same first-rank results happened in SSIM values of Case 6. As for rest metrics in the PU datasets, gratifying top-ranking results in all bands are also achieved. Meanwhile, the average optimum of all bands is also easily derived by SLRL4D.   Table 1 demonstrates all the five quantitative evaluation metrics of seven competitive methods in WDC and PU datasets under all simulated cases. Due to the poor results of BM4D and LRMR, we have omitted them in this table. The optimal value of metrics among all competitive methods is highlighted in bold. As shown in Case 1 and Case 2, it is obvious that NLR-CPTD achieves second optimal PSNR, SSIM and FSIM values encountering heavy Gaussian noise.
The reason is that self-similarity of HSI is sufficiently utilized by the non-local regularization. However, our SLRL4D still performs better even NLR-CPTD reached competitive results. The superiority of non-local BM4D filtering is revealed compared with patch matching of NLR-CPTD. Meanwhile, as a BM4D filtering based method, PCABM4D achieved only average results. Both above phenomenas indicate that the combination of subspace low-rank learning and non-local BM4D filtering is sincerely powerful for getting rid of heavy Gaussian noise. Case 3-6 show that our proposed method almost achieved global optimal quantitative results in all cases in contrast with every competitor. Thus, the ability to remove mixed noises of SLRL4D is verified convincingly. NLR-CPTD performs delightful in heavy Gaussian noise, but its restoration results on mixed noises case is relativity inferior. LRTDTV and LRTDGS achieved pleasant consequences in mixed noises case, but compared with subspace-based SLRL4D, a conspicuous divergence still exists.

Qualitative Evaluation Results
For further illustrating the ability to preserve spectral signatures while restoring clean HSI of the SLRL4D method, the spectrum difference curves between original and restored HSI are presented in Figures 8-11. As shown in Figure 8b, extensiveness fluctuations occurs at the spectrum of pixel (92,218) in WDC under the contamination of heavy Gaussian noise. It is clear that our SLRL4D achieved most parallel signature of spectral after restoration. Figure 9 shows that in front of intricate mixed noises, BM4D, PCABM4D, TDL and NLR-CPTD disappointed the retrieval of spectral signatures.   TDL [58] FSSLRL [49] LLRSSTV [30] LRTDTV [34] LRTDGS [57] NLR-CPTD [41]  This is unsurprising, particularly given their poor performance on mixed denoising. Combination methods of low-rank and other regularization terms that is, FSSLRL, LLRSSTV and LRTDTV achieved suboptimal spectral signatures preservation, but some middle-scale fluctuations still exist. Both LRTDGS and our SLRL4D makes splendid preservation of spectral signatures, but whether in visual comparison or quantitative analysis, far better results are always accomplished by SLRL4D.
Regarding to the PU dataset, approximate phenomena occurred as illustrated in Figures 10 and 11. In either Case 1 or Case 6, the proposed method performs far better stability on spectrum difference images than any competitive method without even a hitch. Figure 10k-l indicates that while NLR-CPTD shows suboptimal ability on the removal of heavy Gaussian noise with the help of non-local regularization, its ability to maintain spectral correlation is far inferior to SLRL4D. As for mixed noises degraded Case 6, Figure 11 and Table 1 collectively indicate that our SLRL4D not only achieved better quantitative analyzed results than both LRTDTV and LRTDGS, the ability to restore spectral signatures of SLRL4D still exceeded them in the PU dataset. In conclusion, the proposed SLRL4D achieved not only optimal visual and quantitative performance but also transcendental preservation of spectral signatures.

Classification Evaluation Results
In this evaluation, we conducted the classification experiments on the restored Pavia University dataset to further verify the superiority of the proposed method. Representative case 1 and case 6 are employed. Table 2 shows the classification results of restored HSI by every competitor and the proposed SLRL4D, respectively. The well-known support vector machine (SVM) [59][60][61] classier with RBF kernel is used for classification. Three indexes, that is, overall accuracy (OA), average accuracy (AA), and kappa statistic (Kappa), are employed for quantitative evaluation. Before classification, about 10% of samples are randomly selected from the labeled samples for training and the rest samples are for testing. Ten times of Monte Carlo runs are conducted and the averaged results are listed in Table 2. It leads us to conclude that restoration does improve the accuracy of classification, and better restoration delivers higher accuracy. The classification results once again verified that the proposed SLRL4D can accurately and effectively restore the HSI from mixed noise. Table 2. Classification accuracy for Case 1 and Case 6 of Pavia University dataset.

Case
Index Noisy PCABM4D [47] TDL [58] FSSLRL [49] LLRSSTV [30] LRTDTV [34] LRTDGS [57] NLR-CPTD [41] SLRL4D In recent years, deep learning methods have earned a prestigious position that cannot be ignored in the field of computer vision [62][63][64][65][66]. Aiming to further verify the capacity of the proposed SLRL4D to remove heavy Gaussian noise, two state-of-the-art deep learning HSI denoising methods HSID-CNN [15] and HSI-SDeCNN [16] are employed for comparison. Similar to most deep learning methods, HSID-CNN and HSI-SDeCNN are mainly designed to remove Gaussian noise. Therefore, for a fair comparison, in the following experiments, we used the exact same WDC dataset and simulation manner of heavy Gaussian noise as in the article corresponding to HSI-SDeCNN [16], that is, the noise level σ n was set to [50,75,100]. The detailed quantitative evaluation results are shown in Table 3. Again, the optimal results are shown in bold font. It can be seen that the SLRL4D realizes better restoration than both latest deep learning methods from quantitative assessment.

Visual Evaluation Results
The aforementioned discussion demonstrated the outstanding ability in mixed denoising on the simulated datasets. However, noise degradation is far more intricate than simulation in real scenes. Figure 12a shows the representative bands 103, 105, 144, 148 and 210 of the well-known real Urban dataset. In these five bands, the intensity of mixed noises is strengthened gradually. Except for band 103, which is slightly contaminated, the other four bands have been completely contaminated by a variety of heavy mixed noises, attempt to acquire any valuable information is absolutely impossible. For slightly contaminated band 103, PCABM4D, TDL, FSSLRL and LLRSSTV failed to remove the cross stripes, rest three comparison methods achieved the removal of stripes but spatial details contrast of some area are slightly enhanced as observed in the enlarged boxes. Both denoising and details preservation were achieved by our SLRL4D method. For band 105, as the noise density of mixed noises increases, PCABM4D, TDL and NLR-CPTD completely lose their power. Low-rank based methods FSSLRL and LRTDTV cannot deal with heavy stripes. LLRSSTV and LRTDGS achieved fine restoration results, but as shown in the enlarged boxes, distortion happened in the results, and, some areas became too bright in contrast with SLRL4D. For band 144, even LLRSSTV failed to remove the cross stripes. Similar to band 105, LRTDGS caused anamorphoses after restoration as shown in the enlarged box. The proposed SLRL4D successfully achieved delightful preservation of image details. Similar restoration results re-emerged on band 148. The difficulty of restoration was further increased, but SLRL4D still successfully accomplishes its mission without causing any abnormal brightness. As for band 210, SLRL4D also achieved optimal restoration results, the combination of subspace low-rank learning and subspace non-local filtering proved itself again.

Quantitative Evaluation Results
The horizontal mean profiles of bands 103, 105, 144, 148 and 210 on the HYDICE Urban dataset are shown in Figure 13. In each subfig, the horizontal axis represents the bands of urban dataset, and the vertical axis represents the digital number of horizontal mean profiles. As shown in Figure 13a, rapid fluctuations exist in the curve of the original bands due to the contamination of heavy mixed noises in all five bands. Particularly in the last three bands, severe mixed noises caused an abnormal jitter. As shown in Figure 13b-e and h, PCABM4D, TDL, FSSLRL, and LLRSSTV achieved inferior restoration performance. For non-local method NLR-CPTD, fine results are achieved only when the noises level are relativity mild. Figure 13f,g indicates that both LRTDTV and LRTDGS achieved relatively flat horizontal mean profiles curves. However, the blue box shows that compared with these two methods, our proposed method SLRL4D achieved flatter curves in all five bands, which indicates that the cross stripes and other diverse noises of Urban dataset have been eliminated more efficiently. In sum, competitive restoration results of our proposed SLRL4D were verified. To further illustrate the efficiency of our proposed SLRL4D method in real scenarios, a blind non-reference assessment index named Q-metric was employed in our experiment to evaluate the restoration of spatial information. Table 4 lists the average Q values of all bands in both HYDICE Urban and EO-1 Hyperion datasets. A higher value of Q-metric yields a better restoration result.
The optimal values are shown in bold. It can be seen that our proposed method SLRL4D acquires the best restoration indexes compared with all seven competitive methods.

Dataset
Noisy PCABM4D [47] TDL [58] FSSLRL [49] LLRSSTV [30] LRTDTV [34] LRTDGS [57] NLR-CPTD [ The running time of an algorithm can effectively measure its efficiency. Table 5 lists the running time comparison of all methods on two real datasets. Each method is run 10 times and averaged. In order to ensure the fairness of comparison, all methods are conducted in the same experimental environment, and the iteration times of those methods that require iterative solution is set to 50. As shown by the bold values, TDL achieved the fastest running speed of all methods, but its restoration results are also worse than most methods. The proposed SLRL4D achieved a competitive running time while ensuring the denoising performance.

Discussion
In our SLRL4D model, multiple parameters need to be discussed. With the help of BM4D filtering, no parameter of the non-local term requires to be tunning arduously, which makes our algorithm far more efficient and stable than exist non-local methods. Only subspace dimension p, penalty parameter µ, regularization parameters γ and λ, which impact the nuclear norm and sparse norm respectively are left to be discussed. Obviously, less dependence on the parameters is achieved by our model in contrast with the congeneric methods.
6.1. The Impact of Parameters µ, γ and λ Remarkable impacts were caused by the penalty parameter µ, as with the soft-thresholding parameter λ and singular value thresholding parameters γ. More parameters bring more complex tunning, in order to make the method easier to use, we empirically set µ= 10/ max (d, l) √ d + √ l at first [49]. For regularization parameters γ and λ, extensive experiments were conducted in the interval [0. 1,2] and [0.01, 0.1] to evaluate their sensitivity. As shown in Figure 14, mean PSNR, mean SSIM and ERGAS are employed to demonstrate the restoration performance of different parameter combinations on the mixed noises degraded Case 6 of the WDC dataset. It is easy to find that the proposed SLRL4D achieved delightful restoration results on a large scale of parameters, which indicates that our method does not require extremely specific parameter settings. Therefore, we recommend choosing γ and λ within the range of [0.1, 1] and [0.01, 0.05] respectively.

The Dimension of Subspace p
Theoretically, the number of endmembers in an HSI is the perfect candidate for the dimension of subspace p. The HySime signal identification algorithm is a highly effective tool to estimate p. For the simulated Case 6 under WDC dataset, HySime gives p = 3 as the estimated dimension. However, as illustrated in Figure 15, both mean PSNR and MSA achieved the best results at p = 5 instead of p = 3. The primary cause of this phenomenon is that HySime only delivers precise estimation for Gaussian degraded signals. Therefore, for the case of mixed noises, we recommend choosing subspace dimension p within the range of [3,8] according to the specific degradation status itself.
All our experiments select parameters in the above ranges, for instance, the parameters used on simulated case 1-6 of the WDC dataset are listed in Table 6.   Figure 16 shows the rising tendency of the PSNR and SSIM under simulated Case 6 of WDC dataset along with the increase of iterations, both the PSNR and SSIM achieved convergence after rapid growth in the initial few iterations, hence the rapid convergence of the proposed SLRL4D is verified.

Conclusions and Future Work
Focusing on the mixed noises degeneration problem of HSI, in this paper, a joint method of subspace low-rank learning and nonlocal 4-d transform filtering is presented for HSI restoration. In summary, both spatial-spectral correlations of HSI are explored comprehensively through joint low-rank regularization and non-local 4-d filtering in the subspace decomposition framework. In this manner, the low-rank property in the spectral domain and the non-local self-similarity of the spatial domain are exploited simultaneously. In contrast with utilizing prior in original space, our proposed subspace-based method not only obtained better restoration results but also cost lower computational resource. The carefully designed restoration model is efficiently solved by the convex optimization algorithm. Extensive experimental evaluations on multiple simulated and real datasets validate that our proposed HSI restoration method is both effective and efficient. For future work, the tensor representation based subspace low-rank learning methods will be considered to further improve the performance, and more precise non-local self-similarity regularization of spatial domain of HSI will be explored. Meanwhile, using low-rankness for image compression is also worth further study.

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

Abbreviations
The following abbreviations are used in this manuscript:

HSI
Hyperspectral image LRMR Low-rank matrix recovery NSS Non-local self-similarity CP CANDECOMP/PARAFAC BM4D Block-matching and 4D filtering SLR Subspace low-rank ADMM Alternating direction method of multipliers PSNR Peak signal-to-noise ratio SSIM Structure similarity FSIM Feature similarity ERGAS Erreur relative globale adimensionnellede synthese MSA Mean spectral angle