Abstract
For the hyperspectral image (HSI) denoising problem, a symmetric proximal alternating direction method multiplier (spADMM) is proposed to solve the sparse optimization problem which cannot be solved accurately by traditional ADMM. The proposed method finds a high-quality recovery method using the traditional low-rank Tucker decomposition method, which can fully take into account the overall spatial and spectral correlation between HSI bands by using the Tucker decomposition. By choosing appropriate steps to update the Lagrange multipliers twice, it makes the selection and use of variables more flexible and better for solving sparsity problems. To maintain stability, we also add appropriate proximity terms to solve the problem during the computation. Experiments have shown that the spADMM has better results than the traditional ADMM. The final experimental results on the dataset demonstrate the effectiveness of the method.
1. Introduction
The hyperspectral image consists of discrete frequencies with distinct bands. The ability of hyperspectral images to provide additional information to real scenes cannot be captured by the human eye. A wide range of applications, including classification [1, 2], anomaly detection [3], and agricultural production [4, 5], have been explored by researchers. Hyperspectral images are always subject to random noise caused by streak noise, photon effects, low spatial resolution and blurring, and so on. There are various degradation phenomena. Improving hyperspectral image quality by improving hardware solutions is not economically viable and impractical. Therefore, before subsequent applications, it is naturally necessary to introduce methods based on image processing to obtain high-quality HSIs. An equation representing HSI recovery can be written as follows mathematically:where is an image with low spatial resolution, is the original, highly spatially resolved image, stands for the sparse noise, stands for the Gaussian noise, and is the spatially linear degradation operator.
First proposed by Rudin et al., the total variational regularization (TV) method can be viewed here [6]. It can effectively preserve the spatial sparsity as well as segmental smoothness of images and is widely used for image processing tasks such as super-resolution, denoising, and magnetic resonance. As a result, many scholars have used the TV method over the past few years. A spectral-spatial adaptive total variational (SSAHTV) denoising method was proposed by Yuan et al. for the denoising of spatial noise and spectral information differences [7]. He et al. combined the advantages of total variational and low rank and proposed a total variational regularized low-rank matrix decomposition [8]. In a subsequent study, Wei et al. proposed the spatial-spectral full-variance regularized local low-rank matrix recovery method (LLRSSTV), in which each directional anisotropy of HSI is exploited as total variational [9]. Kong et al. developed a tensor-weighted kernel norm minimization, which models nonlocal self-similarity as well as global smoothness in spatial-spectral data [10]. Moreover, other methods of total variational denoising for hyperspectral images are proposed [11, 12].
Low-rank matrix factorization is becoming increasingly popular as a tool for image analysis. A low-rank model describes the discovery and exploitation of low-dimensional structural problems in high-dimensional data. A hyperspectral image usually exhibits strong correlation between adjacent spectral bands, as well as between adjacent spatial pixels, both of which indicate that hyperspectral images have a low-rank structure. A variety of hybrid noise removal methods based on low-rank matrix modeling are presented. A robust principal component analysis (RPCA) method for low-rank matrix approximation (LRMA) was proposed by Cands et al. [13]. On the basis of the LRMR method, He et al. developed the iterative low-rank matrix approximation method (NAILRMA) [14]. HSI denoising has widely used this model and its extensions [15, 16]. Reference [8] proposed TV-regularized low-rank matrix factorization (LRTV). However, LRTV does not provide satisfactory recovery results for the moderately noisy HSI band. In addition, Sun et al. proposed the low rank of spectral difference image (LRRSDS) method using low-rank representation on the hyperspectral first-order difference operator [17]. Low-rank representation on spectral difference image (LRRSDS) were proposed using the hyperspectral first-order difference operator as a low-rank representation for spectral difference images [18]. However, some of the above-given algorithms tend to use chunking to process images, which tends to ignore spatial information and cannot maintain both spectral and spatial properties.
Tensor decomposition and tensor-rank for low-rank tensors denoising methods developed in recent years can maintain better spectral-spatial information. Since the Schatten-p parametrization is closer to the low-rank components, Xie et al. used tensor SVD decomposition to obtain matrix eigenvalues and proposed the weighted schatten p-norm low-rank matrix approximation algorithm (WSN-LRMA) [16]. Zheng et al. proposed double factor regularized low-rank tensor factorization (LRTF-DFR) by performing double factor constrained optimization of the decomposed tensor [19]. Xue et al. proposed a nonlocal low-rank regularized CP tensor decomposition (NLR-CPTD) using two a priori information, the global correlation across spectrum (GCS), and nonlocal self-similarity (NSS) [20]. Then, they proposed a new multilayer sparsity-based tensor decomposition (MLSTD) based on CP decomposition [21]. Recently, this team proposed another method to encode the sparsity of the general tensor using the Tucker decomposition based on the Laplace scale mixture (LSM) modeling of the trivial transformation (TLT) of the factor subspace a priori [22]. Zeng et al. combined the spatial-spectral full variational regularization ( SSTV) with local t-SVD tensor decomposition to propose a hyperspectral image recovery model [23]. Chen et al. proposed weighted group sparsity-regularized low-rank tensor decomposition model (LRTDGS) using the group sparse [24]. Recoveries based on low-rank matrices, a hybrid noise removal method for HSIs was developed by Wang et al. [25]. First, the method separated the pure HSI from the noise-contaminated raw observations using low-rank Tucker decomposition, and the correlation between the hyperspectral bands was exploited to its fullest potential. Then, SSTV regularizer was applied to a low-rank Tucker decomposition framework.
According to the studies mentioned above, most of the TV-regularized low-rank tensor decomposition model were solved using the alternating direction method of multipliers (ADMMs). However, the separable structure of the objective function is not fully utilized by ADMM. Consequently, the objective function’s special properties cannot be fully exploited. This paper uses the symmetric multiplicative alternating direction method with the addition of the proximity operator (spADMM) to compute a better image recovery model based on the literature [25]. The main feature of spADMM is that its pairwise variables are symmetrically updated after each update of the separated original variables, which makes the choice of pairwise variables more flexible than the traditional ADMM and can make full use of the sparse decomposition in the model. The proximity point operator ensures the stability of the algorithm. Thus, the method obtains better optimal sparse and low-rank solutions for image processing.
This study contributes the following:(1)We introduce a symmetric ADMM for tensor decomposition models with TV regularization that solves low-rank tensor decomposition models(2)In an iterative process, we add the proximity operator to limit the steps and obtain better convergence by adding the proximity operator(3)Based on experimental findings, the spADMM has better performance than the traditional ADMM, and the new computational method outperforms the comparison algorithm by a significant margin
The other parts of this paper are organized as follows. The low-rank tensor decomposition with SSTV regularized model is introduced in Section 2. In Sections 3-4, the performance of the proposed method is evaluated through simulated and real data. Finally, in Section 5, the conclusions of the paper and future work are presented.
2. Proposed HSI Denoising Model and Algorithm
2.1. Notations and Preliminaries
We first give the symbols to be used throughout the paper [26–28]. A third-order tensor, denoted by ∈ , is considered in this paper. The three indices are called ways or modes, and n-mode is the n-th dimension of this tensor. is a real manifold. The elements of this tensor can be expressed as with . The tensor has a norm equal to the square root of the sum of the squares, which is
A matrix Frobenius norm is denoted by for a matrix. The inner product of two tensors of the same size is given by
The following definitions are . Tensor with a matrix has an n-mode product denoted by . Elementally, this corresponds to
Readers interested in a more detailed introduction can refer to [29].
2.2. SSTV with Anisotropic Low-Rank Tensor Decomposition
TV regularization is one of the effective methods to process images. Based on TV, denoising algorithms have become a popular denoising method due to their effectiveness in maintaining spatial segmentation and edge information smoothing [12, 30]. In [15], the hyperspectral image was recovered using a spatial-spectral TV method. An iterative TV architecture for reconfiguring HSI was proposed by Kuitating et al. [31]. Using the SSAHTV model, Yuan et al. proposed a HSI recovery algorithm that considered both spectral noise differences and differences in spatial information [7]. Recently the combination of low-rank constrained methods as well as TV regularization is often used for hyperspectral image denoising in [8, 32, 33]. This spectral smoothing is ignored by TV regularizers, so a new SSTV regularizer is designed to explore the segmental smoothing structure both spatially and spectrally [25].
Gaussian noise is often modeled using the Frobenius parametrization, which can improve the model’s performance in some situations. The parametric term can be used to filter sparse noise, such as impulse noise. There is a need to estimate the underlying HSI from its observation, which is weakened by a variety of types of noise. The model in [25] has TV-regularized low-rank tensor decomposition combined with low-rank tensor decomposition in spatial and spectral modes:where Tucker decomposition using a core tensor and factor matrix of rank .
The SSTV term iswhere is its regularization strength along the th mode, and is the entry in .
Next, we introduce a new ADMM algorithm with a proximal point optimization procedure to solve (5) in the following.
2.3. Procedure for Optimization
The optimization problem in (5) can be solved using a variety of methods. Here, we propose a new ADMM method for efficient solution. As a first step, we convert (5) in [25] into the equivalent problem:where as for an HSI cube, is the weighted three-dimensional difference operator, and represents the first-order difference operator.
To solve the optimization problem in (7), let and . The augmented Lagrangian model of (7) can be rewritten by introducing some auxiliary variables aswhere is a penalty parameter and , and are the Lagrange multipliers.
First, we show that solving the model (7) with the classical ADMM. Within the framework of ADMM ([25]), the variables involving the model (7) are alternatively updated as
Now, mainly inspired by [34–36], we use a generalized ADMM with proximity point operator, which can guarantee the step size in a wider convergence domain. Following is a brief description of our approach:
Next, we consider how to get the solution of subproblems in (10). The subproblem is
Then, the above problem can be solved by the following formulation:
Let , where , and can be updated by using the classic algorithm for higher-order orthogonal tensors (HOOI) [29]. Then, can be updated with the following equation:
The subproblem is
To solve (14), the subproblem can be solved aswhere the adjoint operator of is indicated by . The three-dimensional (3-D) FFT matrix can be used to diagonalize . As a result, we have the following equation:where the fast 3-D Fourier transform and its inverse transform are denoted as fftn and ifftn, respectively. The division is carried out in an elemental manner, and the square of the element is .
The subproblem is
Soft threshold operators can be used to solve the subproblem (17):where the soft threshold operator is defined asand , .
The subproblem is
The subproblem (20) can be solved using the soft threshold operator
The subproblem is
The following closed form solutions exist:
As mentioned above, a new ADMM is utilized to solve the model as shown in the Algorithm 1.
|
For model (7), we propose a new iterative algorithm different from [25]. The new iterative algorithm is obtained by adding appropriate proximity point terms and pairwise variables to the conventional ADMM. For the purpose of verifying the method’s effectiveness, the following experiments have been conducted.
3. Experiment
On HSIs, we evaluate the performance of the proposed spADMM. It is necessary to normalize the testing HSIs band by band to [0, 1]. The methods of comparison are as follows:(i)Through TV regularization, LRTV [8] eliminates Gaussian noise by combining the low-rank constraint with TV regularization. Efficient removal of sparse noise using low-rank constraints.(ii)LRMR [15] is used to separate the target HSI into a series of full-band cubes. Each full-band cube is then transformed into a matrix. The expansion matrix has a low rank as a result.(iii)Using the Tucker decomposition, the LRTDTV [25] considers global correlations in the HSI and uses the SSTV to describe the segmented smoothing prior as well as the spatial and spectral modes.(iv)TLR- SSTV [23] combines the advantages of global regularization and local LTR models to propose regularized local LTR model for hyperspectral recovery.(v)SLRL4D [37] is a method based on subspace low-rank learning and nonlocal four-dimensional transform filtering. Both spectral and spatial correlations are explored by performing low-rank learning and nonlocal 4-d transform filtering on the subspace.(vi)LRTDGS [24] combined low-rank Tucker tensor decomposition and group sparse regularization for the HSI restoration.(vii)FRCTR-PnP [38] introduced a fibered rank constrained tensor restoration framework with an embedded plug-and-play (PnP)-based regularization.
Each parameter of the methods is set according to the suggestions in their articles or according to the code written by the authors.
In this paper, three evaluation metrics are mean peak signal-to-noise ratio (MPSNR), mean structural similarity (MSSIM), and error relative global error (ERGAS) for quantitative evaluation to measure the quality of denoising results. The following are the definitions of these three metrics. First, MPSNR is given by the following equation:where the size of the image in space is represented by , the original image is represented by , the distorted image is represented by , is the maximum pixel value attainable, and is the PSNR number. MSSIM is defined aswhere is a constant for ; is similar to ; and represent the standard deviation of and ; the original image and the distorted image is represented by and ; the jth local contents is represented by and ; and the number of local windows is .
ERGAS is defined aswhere the reference image and the recovered image of bands are denoted by and , and the number of bands is denoted by .
In image processing and computer vision, the PSNR and SSIM are two conventional perceptual quality indexes (PQIs). The closer the target HSI is to the reference HSI, the greater the difference between the two. The ERGAS is able to measure the global radiometric distortion between two images. ERGAS differs from the previous two measurements in that the smaller its value, HSI with better spectral quality, so the index provides a quick and accurate picture of the overall spectral similarity between the two images.
3.1. Simulated Experiments
In this section, we use three subimages as test data. For Pavia City Center, the dataset size is . For Simulated Indian Pines, the dataset size is . For Cuprite, the dataset size is .
The mixing of various types of noise often degrades HSIs. So to simulate realistic noise scenarios, we consider various combinations of Gaussian noise with varying standard deviations, impulse noise with varying intensities, and stripes noise.(1)Case 1. (i.i.d. Gaussian noise): the test data set is contaminated with independent identically distributed zero-mean Gaussian noise of different intensities, where the variance of the Gaussian noise takes a value of 0.2;(2)Case 2. (i.i.d. Gaussian noise + Stripes with different scales): all bands are affected by identically distributed zero-mean Gaussian noise and deadlines. The variance of the Gaussian noise takes the value of 0.4. And, there are some deadlines added in the 60–80 band, where the number of stripes is generated randomly from 3 to 10, and the width is generated randomly from 1 to 3;(3)Case 3. (i.i.d. Gaussian noise + i.i.d. impulse noise): all frequency bands are corrupted by independent zero-mean Gaussian noise and impulse noise. The variance of the Gaussian noise takes the value of 0.4. And, the impulse noise takes the value of 0.02;(4)Case 4. (Non-i.i.d. Gaussian noise + i.i.d. impulse noise): nonindependent zero-mean Gaussian noise and impulse noise are present in all bands in this case. The variance of the Gaussian noise takes the value U (0.3.0.4). And, the impulse noise takes the value 0.015;(5)Case 5. (Non-i.i.d. Gaussian noise + deadlines): nonindependent zero-mean Gaussian noise and deadlines affect all bands in this case. The variance of the Gaussian noise takes the value of U (0.4.0.5). And, there are some deadlines added in the 60–80 band, where the number of stripes is generated randomly from 3 to 10, and the width is generated randomly from 1 to 3;(6)Case 6. (Non-i.i.d. Gaussian noise + Non-i.i.d. impulse noise): this case involves nonindependent identically distributed zero-mean Gaussian noise and nonindependent identically distributed impulse noise across all bands. The variance of the Gaussian noise takes the value U (0.2.0.3) and the impulse noise takes the value U (0.02.0.03);(7)Case 7. (Non-i.i.d. Gaussian noise + Non-i.i.d. impulse noise + Stripes with different scales): the nonindependent identically distributed zero-mean Gaussian noise, the nonindependent identically distributed impulse noise, and the deadline corrupt all bands. The variance of the Gaussian noise takes the value U (0.4.0.5), the impulse noise takes the value U (0.01.0.02). Stripes are added to bands 61 to 80, with the number of stripes being randomly selected between 20 and 40.
The results for the different HSI data are presented in this section. MPSNR, MSSIM, and ERGAS in Pavia City Center are listed in Table 1. MPSNR, MSSIM, and ERGAS in Simulated Indian Pines are listed in Table 2. MPSNR, MSSIM, and ERGAS in Cuprite are listed in Table 3. The best results are highlighted in bold. The denoising results for band 73 of the Pavia City Center HSI are shown in Figure 1. The denoising results for band 61 of Simulated Indian Pines HSI are shown in Figure 2. The denoising results for band 61 of Cuprite HSI are shown in Figure 3.



For Pavia City Center, the MPSNR, MSSIM and ERGAS values of the proposed method and other comparison methods under different noise environments are given in Table 1. Calculated results show that MPSNR and MSSIM are improved, while ERGAS is reduced with noise suppression, indicating that the method proposed in the paper is effective. FRCTR-PnP and SLRL4D have better denoising effect in Cases 1 and 6. It can be seen that the denoising effect of the proposed spADMM is better in the case of large Gaussian noise. For calculation methods with similar effects, such as FRCTR-PnP and SLRL4D, our proposed method has shorter calculation time.
For different noise scenarios, Table 2 gives MPSNR, MSSIM, and ERGAS values for the proposed method and other compared methods in Simulated Indian Pines. Calculated results show that MPSNR and MSSIM are improved, while ERGAS is reduced with noise suppression, indicating that the method proposed in the paper is effective. LRTDGS performs better in partial noise cases, but our proposed algorithm is only second to its performance. Compared to other methods, the proposed spADMM is more effective.
We give the recovery results for Cuprite in Table 3. From the quantitative evaluation of recovery, FRCTR-PnP, SLRL4D, and our proposed method achieved good experimental results compared to other comparative methods, but our proposed method had the least calculation time.
Figures 1–3 show that the LRMR method is applicable except for the edge regions. As LRMR is a patch-based denoising method, it may lose multidimensional and detailed information. In selected bands, LRTV, and LRTDTV can remove most of the noise. However, LRTV and LRTDTV cannot keep edge information and local detail information. LRTDGS is not able to recover local details well, especially some edges in Pavia City Center. TLR- SSTV is not able to fully remove the noises. FRCTR-PnP and SLRL4D outperform Simulated Indian Pines in image denoising at Pavia City Center. FRCTR-PnP, SLRL4D, and our proposed SPADMM method all perform better in Cuprite data. It is evident from the experimental results that the proposed spADMM can still achieve the better results in most cases.
3.2. Real Data Experiments
In this section, we adopt real Indian Pines data collected by the AVIRIS sensor in north-western Indiana. For Indian Pines, the dataset size is . Comparisons are made in real HSI experiments using the same competitive methods as in simulated experiments. Before the denoising process, the grayscale values of each HSI band are normalized to [0, 1]. The real Indian Pines data mainly contain Gaussian noise and sparse noise.
The Figure 4 shows the grayscale images of the 146th band in the real Indian Pines data and the different visual effects. LRMR, LRTV, LRTDTV, LRTDGS, TLR- SSTV, FRCTR-PnP, SLRL4D, and the proposed spADMM all have good denoising capabilities. However, the proposed spADMM is better in recovering local details especially at the edges.

4. Discussion
4.1. Parameter Analysis
A parameter analysis of the proposed method is presented in this paper, including the regularization parameters and , the weights , and the parameters , , , , , , , and . Simulation data from Case 1 are used for the analysis in Pavia City Center data.
Regularization parameter : In impulse noise, controls its sparsity. We set , where height and width of the HSI band are and , respectively. Figure 5 shows the result of HSI recovery with taking values of trials in the set {1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100}. Figure 5 shows that PSNR and SSIM are highest when when . Therefore, we recommend using in the simulated data experiments.

Regularization parameter : Gaussian noise is limited by , a parameter used to limit its sparsity. Figure 6 shows the recovery results, with taking values in the set {10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 110, 120, 130, 140, 150 } for the test. According to the figure, the method produces the highest PSNR and SSIM values when . As a result, simulations with simulated data should use .

For the TV regularization parameter , the Figure 7 shows the recovery results for selected in the set {0.5, 1, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0, 2.1, 2.2, 2.3, 2.4, 2.5, 2.6, 2.7, 2.8, 2.9, 3}. By observation, the test sets the parameter to 2.3.

The robustness of , , , , , , , and are analyzed in Figure 8. Specifically, the method obtains the highest PSNR at , , , , , , , and values.

Similar to weights in [25], weights , and are chosen for SSTV regularization. The values of and are set to 1, and the values of are adjusted in the interval [0, 1].
4.2. Comparative Analysis of spADMM and ADMM
Our proposed spADMM performs symmetric updates after each update of the separated original variables, which makes the selection of paired variables more flexible than the traditional ADMM and makes full use of the sparse decomposition in the model. In Figure 9, we compare the number of nonzero singular values obtained by the spADMM and the traditional ADMM in Cuprite and Pavia City Center. The spADMM has got fewer nonzero singular values than the traditional ADMM. The spADMM has advantages to guarantee the optimal low-rank solution than the traditional ADMM. In Figure 10, we calculated the value of in Cuprite and Pavia City Center. It can be clearly observed that the values of spADMM are closer to the original image in most bands, which represent the spADMM further improves the sparsity compared to the traditional ADMM.


(a)

(b)
4.3. Complexity Analysis
In this method, we mainly consider the time spent on and since and other parameters are solved using simple algebraic calculations. Let the size of the noise hyperspectral image is . For the update of , several SVDs need to be executed using the HOOI algorithm to obtain an estimate of . The computational complexity is as in HOOI, where is the number of iterations. For a matrix of size , the computational complexity of SVD is . According to the computational complexity of matrix multiplication, we obtain the computational complexity of as . The computational complexity of FFT is , where D is the data size. Thus, the main cost of updating the subproblem of is . The above-given data show that the complexity of the proposed algorithm is .
4.4. Convergence Analysis
We analyze the convergence numerically of the proposed spADMM algorithm in a test at Simulated Indian Pines, Pavia City Center, and Cuprite. A plot of the relative change value R of the hyperspectral reconstructed images versus the number of iterations of the spADMM algorithm is shown in Figure 11. It is evident that R converges to zero as the number of iterations increases. This shows that the proposed method spADMM is numerical convergent.

5. Conclusion
For solving SSTV regularized low-rank tensor decomposition models, we use a new symmetric ADMM incorporating proximity operators. The different types of noise considered in the test images in the numerical experiments: impulse noise, Gaussian noise, and stripes noise of different proportions. The proposed method outperforms some of the comparable methods in terms of removing mixing noise and finely preserving the HSI bands based on extensive experimental results. In future work, low-rank representation methods can make use of the diversity of tensor decomposition, which can be used to improve the performance of HSI spatial domains and explore their nonlocal self-similarity.
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 (Grant nos. 61967004, 11901137, 11961011, 72061007, and 62171147), the Guangxi Key Laboratory of Cryptography and Information Security (Grant nos. GCIS201621 and GCIS201927), the Guangxi Key Laboratory of Automatic Detecting Technology and Instruments (Grant nos. YQ20113 and YQ20114), and the Basic Ability Promotion Project for Young and Middle-Aged Teachers in Universities of Guangxi under Grant no. 2019KY0253.