Abstract
The acquired hyperspectral images (HSIs) are affected by a mixture of several types of noise, which often suffer from information missing. Corrupted HSIs limit the precision of the subsequent processing. In this paper, the weighted Total Variation-regularized Hybrid Model of CP and Tucker (WTV-HMCT) is proposed to accurately identify the intrinsic structures of the clean HSIs. By jointly minimizing CP rank and Tucker rank in the low-rank tensor approximation, WTV-HMCT fully exploits the high-dimensional structure correlations of HSI. To ensure the piecewise smoothness of the recovered image, the hybrid low-rank tensor decomposition approach integrates the weighted spatial spectral total variation regularization for the separation of the noise-free HSI and mixed noise. By the Alternating Direction Method of Multipliers (ADMM), the optimization model is transformed into two subproblems. Finally, an efficient proximal alternating minimization algorithm is developed to optimize the proposed hybrid low-rank tensor decomposition efficiently. The experimental results show that the proposed model effectively handles Gauss noise, striping noise, and mixed noise and that it outperforms the most advanced methods in terms of evaluation metrics and visual evaluation.
1. Introduction
In the process of acquisition and transmission of HSIs, HSIs are usually contaminated by a variety of noises, such as Gaussian noise, stripes, deadlines, impulse noise, and their hybrids [1], which make subsequent analysis and application of HSIs difficult. Therefore, as a preprocessing step, HSI restoration [2–6] is the main challenge faced by many researchers.
The high-dimensional HSI is composed of hundreds of separate images banded together. The classical HSI restoration algorithm is based on low-rank matrix approximation [6], image denoising via spatiospectral total variation [7], and the restoration method with nuclear norm minimization [8]. However, conventional HSI restoration denoising approaches can only investigate the structural characteristics of each pixel or band individually, ignoring the significant correlations between all spectral bands. As a result, the quality of their restored images is substantially low. To solve this issue, many academics have proposed HSI restoration algorithms based on low-rank tensor decomposition [9–13]. Furthermore, low-rank tensor factorization was applied to spectral correlation representations for HSI noise reduction to eliminate the calculation of singular value decomposition and enhance the efficiency [11–13]. These approaches efficiently explore the spatial structure throughout the denoising process, but they fail to capture the high-order low-rank structure of HSI.
1.1. Related Work
Tensors have a different form of decomposition from matrix singular value decomposition, which has applications in many scenarios, such as image dynamic enhancement [14], medical image encryption [15], rib segmentation [16], atrial fibrillation detection [17], and video segmentation [18]. At the moment, tensor decomposition primarily includes Tucker model decomposition [13] and CANDECOMP/PARAFAC (CP) model decomposition [19–21]. Tensor’s singular value decomposition (T-SVD) [9], tensor train (TT) [22], and tensor ring (TR) are the other tensor decompositions [23]. T-SVD is an extension of SVD applied to matrix filling data, in which a third-order tensor is decomposed into a singular tensor with normal property, and the number of channels with data in the diagonal tensor is defined as the tensor rank [9], and the tensor is then recovered by minimizing the tensor’s Tubal rank. The TT decomposition of a tensor sequence is similar to SVD in that it decomposes a higher-order tensor into a set of third-order singular tensors. Tensor sequence-based approaches include rank minimization and variational energy minimization based on ALS to recover tensors. A tensor tree is a hierarchical Tucker model that decomposes a tensor to produce a subset of modules [24]. The TR decomposition corresponds to the TT decomposition, which decomposes a higher-order tensor into a set of cyclic contraction third-order tensors.
As a result, the Tucker model and the CP model are two important tensor decomposition approaches, each of which focuses on revealing the tensor structure in its way. To enhance the overall performance of HSI image reconstruction, two decomposition algorithms are employed to constrain low-rank tensors and more accurately characterize low-rank features.
1.2. Our Contributions
Contributions to this article are as follows: (1)To represent the low-rank structure and sparsity of HSI, a Hybrid Low-Rank Tensor Decomposition (WTV-HMCT) model is presented. The low-rank hybrid tensor approximate decomposition model is used to break down the CP model and Tucker Model to fully use the high-dimensional structure information of HSI(2)By incorporating the weighted spatial total variation into the hybrid tensor low-rank decomposition model, the restored image is guaranteed to have a good local smoothing structure as well as additional spectral complementary information to reduce spectral distortion(3)Using the alternating direction method of multipliers (ADMM), the optimization model is transformed into two subproblems. The rank tensor updating approach is utilized for the two main subproblems, respectively. Experiments demonstrate that this method can effectively deal with Gauss noise, strip noise, and mixed noise and that it outperforms the most advanced methods in terms of evaluation index and visual assessment
1.3. Organization of the Paper
This paper is organized as follows. To facilitate our presentation, we first introduce some notations and preliminaries of tensors in Section 2. In Section 3, the TV-regularized low-rank hybrid tensor decomposition model is introduced. We then develop an efficient ADMM algorithm for solving the proposed model. In Section 4, extensive experiments on both simulated and real datasets are carried out to illustrate the merits of our model. We then conclude this paper with some discussions on future research in Section 5.
2. Low-Rank Tensor Decomposition
2.1. Symbolic and Data Representation
This subsection introduces various tensor-related notes and preliminaries. Tensors are multidimensional data arrays expressed by upper case calligraphic letters, such as the elements of tensor represented by or . Matrix data is represented in capital bold letters, and vector data is represented in lowercase bold letters. In addition, scalars are represented in lowercase letters, such as . An -dimensional tensor of a real number is expressed as , and a three-dimensional tensor contains rows, columns, and tube fibers, defined as , , and , respectively. The two-dimensional portion of a three-dimensional tensor is referred to as a slice and is represented by all subscripts except two subscripts. The lateral, horizontal, and front slips of the three-dimensional tensor are represented by , , and ; is the expansion of the tensor’s mode , which is formed by arranging the mode fibers into columns.
2.2. Tucker Decomposition for a Third-Order Tensor
For a third-order tensor, the Tucker decomposition yields three-factor matrices , , , and a kernel tensor , each factor matrix in mode is called the tensor’s fundamental matrix or principal component in each mode, so Tucker decomposition is also called higher-order PCA, higher-order SVD, etc. Therefore, the Tucker decomposition can be used in data dimension reduction, feature extraction, and tensor subspace learning. For example, the low-rank Tucker decomposition of an HSI image can remove the hyperspectral image noise, we can also use the tensor subspace to select the features of the HSI image and use the Tucker decomposition to compress the data.
The Tucker decomposition represents a tensor as a core tensor multiplied by a matrix along with each mode. For the third-order tensor , the Tucker decomposition is defined by
Factor matrices are usually orthogonal and can be regarded as principal components along with the corresponding mode. The Tucker decomposition is expressed in a matrix form. For example, the expansion of the third-order Tucker decomposition is the matrix form of the Tucker Decomposition:
The rank uniqueness of the tensor cannot be guaranteed by the Tucker decomposition, so additional constraints need to be added. The core tensor is usually required to be as simple as possible so that the principal components of each mode are as orthogonal or sparse as possible.
2.3. CP Decomposition for Third-Order Tensors
CP decomposition is a different form of tensor representation from Tucker decomposition, where a tensor is expressed as the sum of a finite number of rank tensors. A third-order tensor can be decomposed into , where , , and are the matrices of the corresponding vectors in the three rank-one tensor, also known as a factor matrix. According to the factor matrix, the CP decomposition of the third-order tensor can be written as an expansion:
Decomposed by a slice, CP decomposition of a third-order tensor is sometimes written as
Rank decomposition is a special CP decomposition, , similar to the SVD of the matrix. There is no direct method to solve any given rank of a tensor, which is proved to be an NP-hard problem. Therefore, we generally use the tensor’s low-rank approximation to restore the tensor, but we need to satisfy the following conditions: , denote the k-rank of the matrix : Any column is linearly independent of the largest . Compared with the SVD of the matrix, the uniqueness of the rank decomposition of higher-order tensors does not need the guarantee of orthogonality.
For the calculation of CP decomposition, there is no uniform perfect solution to the problem of how many rank-one tensors (components) can be decomposed into the best one. From the effect, the general approach is to start from rank-1, using ALS (alternate Least Square) algorithm, and gradually increase the rank until the error range. If you have a strong application background and prior information, you can specify the value of the rank in advance. Because the ALS algorithm cannot guarantee convergence to a minimum point, or even to a stable point, it can only find a point where the objective function is no longer descending. The initialization of the algorithm can be random or the factor matrix can be initialized to the corresponding outspread singular vector, such as which is initialized to the first left singular vectors . From the two different methods of CP decomposition and Tucker decomposition, it is also easy to see that CP decomposition is a special form of Tucker decomposition: If the core tensor is diagonal, then Tucker decomposition will degenerate into CP decomposition.
3. Proposed WTV-HLRTD Model for HSI Reconstruction
3.1. Problem Formulation
The hyperspectral image with spatial size (spatial length and spatial width) and spectral band number is denoted as , and the HSI image to be recovered is usually degraded by additive mixed noise, mainly including Gauss noise, impulse noise, and band noise. Therefore, the HSI mixed noise degradation model is defined by where represents the noisy HSI image, represents the noise-free HSI image, represents sparse noise, which mainly consists of impulse noise, fringe noise, and truncation noise, and is gauss noise; the four components are third-order tensors of the same size is the space size of each band and is the number of spectral bands, the main task of HSI reconstruction is to restore the latent image from the observed noise image.
3.2. Total Variational Regularization and Low-Rank Tensor Recovery
Recovering a clean image from a noised image is a difficult and ill-posed problem. To solve this problem, unknown variables need to be regularized by prior information to improve the recovery effect. Therefore, the optimal model for removing the mixed noise problem of hyperspectral images is formulated as:
are the prior knowledge regularization constraints describing the denoised hyperspectral image and the sparse noise , respectively. is the regularization parameter, and is the variance of gauss noise density. Typically, is used to constrain the sparse regularization of the sparse noise , such as . According to the linear mixture model, hyperspectral images have a strong correlation in spectral dimension, which shows that has a low rank prior, Therefore, different types of low-rank matrices are used to explore the prior information . All spectral bands are vectorized using the low-rank matrix recovery method. All spectral bands are vectorized using the low-rank matrix recovery method. However, the correlation of spectral-spatial structure is not considered. Since the hyperspectral image is a third-order tensor, it is more reasonable to model it with tensor tools to preserve the details and spatial structure. To simultaneously aggregate the spatial-spectral correlation of hyperspectral images in 3D, a valid Tucker decomposition can be used to constrain the clean image . By choosing suitable regularization parameters, the regularization model (6) is equivalent to the following optimization problem under the constraint of low-order Tucker decomposition, which can be transformed by where and the factorization matrices , , and are orthogonal in space and spectral dimensions.
According to [20], a single band of the hyperspectral image can be considered a grayscale image, so it has local smoothness of the image structure in a spatial dimension. Therefore, TV regularization [20] is introduced into the low-rank matrix tensor decomposition framework using both spatial and spectral prior information. The associated TV normalized low-rank tensor model is mathematically formulated as where is the different low-rank matrix tensor hyperspectral image approximation, is the TV regularization of the hyperspectral image, and and are the regularization parameters, respectively. Due to the effectiveness of TV regularization constraints, many improved versions of TV regularization have been proposed, such as BANDTV [25], space-spectral TV [26], and 3D spectral space Crossover TV [27]. With these extended full-variation regularization constraints, anisotropic space-spectral TV regularization based on the low-rank Tucker decomposition method [1] has achieved excellent results in hyperspectral image reconstruction.
By combining the hybrid low-rank tensor decomposition with SSTV regularization and considering the global spatial-spectral correlation and spatial-spectral smoothness of the hyperspectral image, satisfactory hyperspectral image restoration results can be obtained. However, this still poses some problems. The frequency band TV and the TV regularization model based on the spatial spectrum make full use of the sparse prior of the spatial spectral difference image and usually use the convex norm to describe the sparse prior. In general, sparsity is an effective constraint to promote the smooth structure of each band segment. However, the sparse prior only describes the number of nonzero elements, ignoring the local structure of nonzero elements. Therefore, to make up for this deficiency, space-spectral TV regularization is implemented in the restoration tensor to better describe the prior and improve the restoration effect of the hyperspectral image.
3.3. Description of the WTV-HMCT Model
By integrating CP decomposition and the Tucker model of the tensor, this form is complementary to the characteristic of low rank and sparsity, which can be expressed as where is the rank of tensor based on CP decomposition, is the rank of core tensor derived from tensor based on Tucker decomposition, and is the regularization parameter. This comprehensive consideration of the new measure makes the tensor have two internal sparse configurational nuclear tensors and low rank along with each tensor module in the tensor subspace; thus, the limitations of the Tucker decomposition on sparse elements and the CP decomposition on low-rank features are solved.
In HSI reconstruction, TV regularization based on the spatial domain is widely used to ensure the segmentation and smooth structure of the image. There is also a strong local smooth structure in the spectral domain. In [14, 25–27], because the spectrum is smooth, a 3D TV constrained optimization model is proposed to restore the image. To describe the edge and structure more accurately, the weighted space-spectral domain TV (WTV) is utilized to reserve a better smooth structure.
Therefore, the WTV-HMCT model is described as follows: where is the space and the spectral mode TV, is the weight parameter, and is the two differential operators in the two-dimensional space and and is defined as where represent the position coordinates in 3D, respectively. It is worth noting that the proposed WTV-HMCT model can capture the HSI of spatial and spectral information well and has a strong mixed noise removal ability. In particular, the spectral similarity and structural correlation of all pixels in the two spatial modes can be fully utilized by combining CP model decomposition with the rank of the Tucker Model, the rank estimation of CP model decomposition has good recognition ability for sparse noise items including impulse noise, dead-line noise, and fringe noise, and the WTV norm term is designed to characterize the piecewise smooth structure, which can remove gauss noise in both spatial and spectral domains.
3.4. Decomposition and Optimization of the WTV-HMCT Model
In the WTV-HMCT optimization model (11), because and are nonconvex, it is difficult to solve the rank of the tensor directly, which makes the objective function difficult to be solved. Therefore, the CP kernel norm and the Tucker kernel norm are used as convex proxies for functions and . The approximate convex replacement functions RC and RT are solved, and the model is transformed into the form of where and are the CP and Tucker nuclear norms, respectively. The WTV-HMCT model is a convex optimization problem, and the ADMM algorithm is used to solve the corresponding model iteratively. To solve the convex optimization problem (11), two tensor variables and are introduced and reformulated as equivalent global consistency problems, as shown in where . Then, the ADMM algorithm is used to solve the optimization model. Based on the ADMM algorithm, the constraints of model (13) are transformed into the optimization Lagrangian shown in
According to the optimization and Lagrangian, the model can be decomposed into four subproblems. (1)Solving subproblem, by fixing other variables, the subproblem is shown in
Tucker factorization factors can be well solved by the high-order orthogonal iterative HOOI Algorithm [28]. When you get the factorization factor , then update by . (2)To solve the subproblem, by fixing other variables, the subproblem is shown in
The problem in (16) is solved by using the tensor kernel norm minimization approach to low-rank CP.
The model described in (16) is reformulated as the equivalent form of
The iterative solutions are obtained by using ADMM-BCD (alternate direction method of multipliers and block coordinate descent (ADMM-BCD)).
When the th iteration is carried out, can be directly solved as
Similarly, and can be obtained by the same method.
can be calculated according to the regularized soft threshold algorithm and can be rewritten as
Update where is the contraction operator, and can be obtained by . After getting and , is obtained by (3)Solve the subproblem by fixing other variables, the subproblem of solving can be formulated as
The subproblem (22) can be solved precisely by using the soft threshold contraction operator and can be directly solved as where . (4)Solve subproblem , by fixing other variables; the subproblem can be rewritten as
Similarly, by applying the soft threshold contraction operator, the subproblem in (24), as shown in (25), can be solved precisely: (5)Update multiplier . According to the ADMM algorithm, update Lagrange multiplier be obtained by
In addition, the weight in (12) is updated by
4. Experiment
To verify the performance of the proposed method, this section will conduct experiments based on simulated data to evaluate the performance of the proposed WTV-HMCT model for hyperspectral image denoising. For illustrating the superiority of the combination of the low-rank hybrid tensor decomposition model and SSTV regularization, our methods can contrast the performance of their corresponding five most advanced hyperspectral image restoration methods. These methods include low rank matrix factorization and TV regularization methods (LRTV) [25], low-rank tensor decomposition and TV regularization of anisotropic spatial spectrum (LRTDTV) [25], tensor sparse recovery based on Kronecker expression (KBR) [27], BM4D [28], and WGLRTD [29]. The code for the comparison method is available from the author’s home page. In addition, the choice of parameters for these models in all experiments is carefully determined by the author’s code or the suggestions in the author’s paper to obtain the best performance.
Each band’s pixel values were scaled in the [0, 1] range before the experiment began. All experiments were run in Matlab R2018B on a laptop with an Intel Core i7-8750 h CPU running at 2.20 GHz and 32GB of RAM.
4.1. Simulation Experiment Environment
(1)Experimental environment: To verify the robustness of the proposed WTV-HMCT method, we select two clear HSI data sets for simulation. The first data set is the Washington DC Mall (WDC) at http://lesson.weebly.com/hyperspectral-data-set.html, which selects a subimage dataset of size . The second dataset is a hyperspectral Pavia City data center set (PAC) with a size of pixels.(2)Noise type description: To simulate the complex noise in the real scene, three different types of noise are added to the two original and clean hyperspectral data sets. A detailed description of these noise types is given below.(i)Case 1: Add zero mean Gauss noise to all bands, with a noise variance of 0.15 for each band.(ii)Case 2: Gauss mixed impulse noise is added to each band and Gauss noise is added in the same way. In addition, different percentages are added to each band and the percentages are randomly selected within the range [0, 0.2].(iii)Case 3: A mixture of Gauss noise, impulse noise, deadline noise, and striping noise is added to the hyperspectral image data. Gauss noise and impulse noise are added in the same way as (2). In addition, 20% of the bands were randomly selected from all bands, and the number of bands ranged from 3 to 10.
4.2. HSI Experiments and Analysis
The reconstruction performance of all comparison methods is analyzed and evaluated from three aspects: visual effect, quantity, and quality.
4.2.1. Visual Comparison
To achieve visual comparison, two representative noise cases 1 and 3 are chosen to compare the performance of the different methods. Figure 1 shows the result of band 36 data recovery from the WDC dataset with noise case 1. To make a better visual contrast, a random region in the image was enlarged deliberately, and the enlarged block was placed in the lower right corner of the image to observe the restoration effect of each algorithm model.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
In Figure 1, the recovery results of all comparison methods are compared and analyzed. The results show that LRTDTV, KBR, and BM4D cannot eliminate the noise in the enlarged images, while LRTV can remove the noise, it appears serious oversmoothing and loses a lot of detailed information in the images. On the contrary, WGLRTD and WTV-HMCT have better restoration effects and obvious noise elimination effects. Compared with other methods, the WTV-HMCT model represented in Figure 1(c) can eliminate all mixed noise and preserve the edges and details of the image effectively, and it shows the performance of the proposed WTV-HMCT method in recovering the WDC Dataset is improved clearly.
Figure 2 shows the results of six algorithmic models recovering 30 bands of PAC 6 data contaminated by noise case 3. After visual comparison, it is not difficult to find that the results of LRTV, LRTDTV, and BM4D still contain part of the noise. The WGLRTD removes the noise but loses the detail of the image, resulting in a pattern and verification. WTV-HMCT designed in this section can get rid of all noise and preserve the original details of the image and get a good visual effect. In short, from the two data sets of Figures 1 and 2, WTV-HMCT for image restoration achieved obvious visual effect and better restoration effect.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
4.2.2. Quantitative Comparison
The comparison of the restored visual effect between Figures 1 and 2 shows the effectiveness of the WTV-HMCT designed in this paper. In this subsection, four objective quantitative evaluation indicators prove the good performance of the WTV-HMCT model in all simulation experiments. These indicators are the average peak signal-to-noise ratio (MPSNR) of all bands, average structural similarity (MSSIM) of all bands, average characteristic structural similarity (MFSIM) of all bands, and Spectral Angle Mapper (SAM).
In this experiment, MPSNR, MSSIM, and MFSIM were evaluated by averaging all bands of the two datasets. The results of the four indicators are shown in Tables 1 and 2. The higher values represented the better MPSNR, MSSIM, and MFSIM of the restoration results. Concerning the SAM index, the lower values indicate a better recovery.
The experimental results show that the performance of WTV-HMCT is better than other methods based on tensor decomposition, thanks to the combination of hybrid low-rank tensor decomposition and WTV regularization, as well as the joint minimization of CP rank and Tucker rank. Compared with other algorithms, our WTV-HMCT obtains the best evaluation index in all cases with four indicators.
The sparse structure and low-rank features of the tensor cannot be fully explored due to the traditional decomposition method, which is compensated by the WTV-HMCT model. In particular, all indicators of the WDC data set have been improved, and this method can achieve better results than all other comparison methods because the model makes full use of the same local smoothing property in all spatial and spectral regions, the norm can be used to capture the difference image effectively. To sum up, the robustness and effectiveness of the WTV-HMCT method are verified by the quantitative comparison of three different noise conditions in two noise-free data sets.
4.2.3. Qualitative Comparison
In the discussion above, we quantitatively analyzed the effectiveness of the algorithm by averaging the evaluation metrics. Figure 3 shows the PSNR and SSIM evaluation index curves for case 2 at each band over both datasets.
(a)
(b)
(c)
(d)
Furthermore, Figures 3(a) and 3(c) show the PSNR value of each band for the PAC and WDC data sets. It is evident that the results obtained by the proposed WTV-HMCT method achieve much higher PSNR for almost every band and indicate the robustness of the proposed method. Figures 3(b) and 3(d) show the SSIM value of each band for the PAC and WDC data sets; the proposed WTV-HMCT method can better retain most of the details of each band. Therefore, the WTV-HMCT achieves the highest PSNR and SSIM values in all bands of WDC and PAC data sets.
In summary, extensive simulation experiments demonstrated that the proposed WTV-HMCT method outperforms the existing denoising approaches.
5. Conclusion
Previous strategies of hyperspectral image reconstruction based on low-rank tensor decomposition could not effectively capture the high-dimensional structure of the tensor, making adaption to different mixed noise hyperspectral image denoising challenging. Furthermore, the lack of constraints on the local smoothness of the spectrum results in spectral curve distortion and distortion. In this paper, a hybrid tensor decomposition (WTV-HMCT) is proposed to describe the low-rank structure and sparsity of hyperspectral images. The low-rank hybrid tensor decomposition model, on the one hand, is proposed by integrating the Tucker model rank into the CP model rank. The weighted spatial and spectral total variation regularization is introduced into the low-rank hybrid tensor decomposition model, which constrains the spatial structure information and local smoothness of the spectral information so that the mixed noise is eliminated. The model is solved by using the joint alternating direction multiplier and block coordinate descent optimization method, which ensures the restored image has a good local smooth structure. From several data experiments, it is shown that the proposed WTV-HMCT method is more effective than the state-of-the-art denoising methods.
In the future, we hope to combine the proposed WTV-HMCT with the deep learning model to learn a more appropriate regularization and improve the ability to remove mixed noise.
Data Availability
The HSI data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This work was supported in part by the Major Research Platform Construction Project Sichuan Province of China under Grant 2019JDPT0014, the Scientific Research Project of Ministry of Education under Grant 2021KSA01008, the Scientific Research Project of Sichuan Higher Education Association under Grant GJXHXXH21-YB-27, and the Science and Technology Support Project of Panzhihua City of Sichuan Province of China under Grant 2021CY-S-6.