I Introduction
Hyperspectral imaging employs an imaging spectrometer to collect hundreds of spectral bands ranging from ultraviolet to infrared wavelengths for the same area on the surface of the Earth. It has a wide range of applications including environmental monitoring, military surveillance, mineral exploration, among numerous others [1], [2]. Due to various factors, e.g., thermal electronics, dark current, and stochastic error of photocounting in imaging process, hyperspectral images (HSIs) are inevitably corrupted by severe noise during the acquisition process. This greatly degrades the visual quality of the HSIs and further affects the precision of previous listed applications. Hence, the task of removing the noise in hyperspectral imagery is a valuable research topic and has received much attention in the past decades.
A natural way for restorazing HSIs is to regard each band as a graylevel image and then apply traditional 2D or 1D denoising methods to remove noise bandbyband. See, e.g., [3, 4, 5, 6]. However, this kind of method ignores the correlations among all the spectral bands or spatial pixels, and thus usually could not provide satisfactory results. Aiming at taking account of such correlations, a variety of studies have been conducted in the literature. For example, Othman and Qian [7] proposed a hybrid spatialsepectral wavelet shrinkage method to take advantage of dissimilarity of the signal regularity in both the spatial and spectral domains of the HSIs. Zhong and Wang [8] designed a multiple spectralband conditional random field method, which can simultaneously model and utilize the spatial and spectral dependences in a unified probabilistic framework. In [9], an efficient HSI denoising procedure is proposed by designing a spectralspatial adaptive total variation model, where the spectral noise differences and spatial information differences were both considered. Additionally, several advanced techniques in traditional image processing, including the nonlocal similarity [10], wavelet shrinkage [11], anisotropic diffusion [12], have been adopted for HSI restoration in the recent years.
Actually, all the aforementioned methods are designed to remove only one or two types of noise, i.e., Gaussian noise, impulse noise, or hybrid Gaussianimpulse noise. In realworld scenarios, however, HSIs are usually corrupted by a combination of several different types of noise during the acquisition process, e.g., Gaussian noise, impulse noise, dead lines, stripes, and many others. Although several methods based on lowrank matrix modeling have been proposed for removing mixed noise in HSIs, e.g., [13, 14, 15, 16, 17], modeling a HSI cube as a matrix is not able to utilize the finer spatialandspectral information, leading to suboptimal restoration results in some heavy noisy situations. See Fig. 1 for an instance. It can be easily to observe that, for the HSI band polluted by severe noise (e.g., band 108), the stateoftheart lowrank matrix method, e.g., LRTV proposed by [17], could not provide satisfactory restoration as for moderately noisy HSI band (e.g., band 103).
Recently, a number of studies, both practically and theoretically, have demonstrated the advantages of direct tensor modeling techniques over matricization techniques in dealing with highorder tensor data. See, e.g., [18, 19, 20, 21, 22] among others. Motivated by such studies, we propose in this paper a novel mixed noise removal approach using direct tensor modeling techniques to fully exploit the spatialandspectral priors underlying the clean HSI part and characterize the intrinsic structures of the heavy noise part. To highlight our contributions, we shall go over related work on HSI restoration using different lowrank modeling techniques and contrast our innovations with the existing literature.
Ia Related Work
In the past few years, various approaches based on lowrank matrix approximation (LRMA) have been proposed for HSI restoration, and can be represented as stateoftheart techniques. Motivated by the idea of robust principle component analysis (RPCA) [23], Zhang et al. [13] explored the lowrank property by lexicographically ordering a patch of the HSI into a 2D matrix and modeled the nonGaussian noise (including impulse noise, dead lines, and stripes) as a sparse part. The socalled “Go Decomposition” (GoDec) [24] algorithm was used to enforce rank and cardinality
constraints for the lowrank and sparse parts, respectively. Considering that the noise intensity in different bands is different, a noiseadjusted iterative LRMA method on the basis of patchwise randomized singular value decomposition was proposed in
[15]. Encouraged by the powerfulness of total variation (TV) regularization in various image restoration tasks, He et al. [17] integrated the bandbyband TV regularization into a rankconstrained RPCA model to explore the spatial piecewise smoothness of the HSI, and as a result enhanced the capability of the LRMA technique for HSI restoration. Similarly, Wu et al. [25] proposed a HSI mixed denoising model by combining the bandbyband TV regularization with the widely used weighted nuclear norm minimization (WNMM) [26, 27]. Aiming at giving better approximation to the lowrank assumption of HSI data, a weighted Schatten norm regularization was introduced by Xie et al. [28] into the LRMA framework. To exploit both the local similarity within a HSI patch and the nonlocal similarity across patches in a group simultaneously, a novel group lowrank representation model was consider in [16]. As stated before, though this kind of LRMA methodology has been an increasingly useful technique in HSI restoration, it fails to fully exploit the prior knowledge on the intrinsic structures of HSI cube after vectorizing the HSI bands.
Despite the efficiency of tensor methods on removing Gaussian noise or hybrid Gaussianimpulse noise in HSIs, e.g., [29, 30, 31, 32, 33, 34], to the best of our knowledge, only two studies [35, 36] based on tensor techniques have been conducted to remove mixed noise in HSIs. Specifically, in [35]
, the tensor nuclear norm introduced by Liu et al.
[18] was used to acquire lowrank property of HSI cube, and the mixednorm was used to impose sparsity property of the outliers and nonGaussian noise; the work
[36] integrated the structure tensor TV [37] into the WNNM model, and demonstrated outperformance over bandbyband TVregularized WNMM method [25], because of the utilization of finer spatial structure information.Basically, our work is related to the aforementioned works, however, there exits significant differences between our work and other ones. Fig. 2 illustrates the frameworks of the popular matrix modeling idea especially LRMA and our direct tensor modeling idea in dealing with noisy HSI cubes. It can be seen that, the matricization techniques should preliminarily vectorize all HSI bands at the cost of losing spatial structures of the HSI cube, whereas direct tensor modeling techniques could deliver more faithful underlying information of the HSI cube without destroying the spatial structures. Despite their connection and similarity, compared with the work [35], we adopt an anisotropic spatialspectral TV to further exploit the local piecewise smoothness in both spatial and spectral domains; compared with the work [36], lowrank Tucker decomposition is used to further characterize the spatial correlation in each HSI band, and the proposed spatialspectral TV is more simple than the structure tensor TV. Additionally, the use of Frobenius norm term for dealing with very large Gaussian noise was not consider in [35, 36]. A detailed motivation of our work can be found in Section III.
IB Our contributions
In this paper, we mainly focus on HSI mixed noise removal. The main contributions of this paper are summarized as follows.

The lowrank Tucker decomposition is applied to separate the clean HSI from the raw observation corrupted by complex nose, which can be fully capable of exploiting the global spatialandspectral correlation among the HSI bands. The classical HOOI algorithm is used to achieve such Tucker decomposition efficiently, without bringing extra computational burden.

A designed spatialspectral total variation (SSTV) regularizer is incorporated into the lowrank Tucker decomposition framework. It is shown that the SSTV regularization has the ability of enhancing the spatial information and thus preserves the spectral signatures of HSIs very well.

The wellknown augmented Lagrange multiplier (ALM) method is used for solving the TV regularized low rank tensor decomposition model. And extensive studies on simulated and real data experiments are carried out to illustrate that the proposed method clearly improves the HSI restoration results over some other popular techniques, in terms of both the quantitative evaluation and the visual comparison.
Notations  Explanations 
, , ,  Tensor, matrix, vector, scalar. 
or  Mode matricization of tensor , obtained by arranging the mode fibers as the columns of the resulting matrix of size . 
Multilinear rank, where ,  
Inner product of tensor and .  
Frobenius norm of tensor .  
Mode multiplication of and U with the matrix representation . 
IC Organization of The Paper
This paper is organized as follows. Section II introduces some notations and preliminaries of tensors, which will be used for presenting our procedure. In Section III, the anisotropic SSTV regularized lowrank tensor decomposition model and its motivations are introduced. We then develop an efficient ALM algorithm for solving the proposed model. In Section IV, extensive experiments on both simulated and real datasets are carried out to illustrate the merits of our model. Finally, we conclude this paper with some discussions on future research in Section V.
Ii Notation and Preliminaries
It is known that a tensor can be seen as a multiindex numerical array, and its order is defined as the number of its modes or dimensions. A realvalued tensor of order is denoted by and its entries by . We then can consider an vector as a tensor of order one, and an matrix as a tensor of order two. Following [38], we shall provide a brief introduction on tensor algebra.
The inner product of two samesized tensors is defined as . Then the corresponding Frobenius norm is defined as . The socalled mode matricization of a tensor is denoted as , where the tensor element maps to the matrix element satisfying with . The mode multiplication of a tensor with a matrix is denoted by and elementwise, we have The multilinear rank is defined as an array where , .
The tensor notations used in this work are summarized in Table I. Interested readers are referred to [38] for a more detailed introduction.
Iii HSI Restoration Via TV Regularized Lowrank Tensor Decomposition
Iiia Motivation
In many real situations, the observed HSI data are contaminated by a mixture of several different kinds of noise [13, 15, 17]. As a result, a noise HSI cube denoted by a three order tensor ,, where each matrix represents th band, with the height and width , and denotes the number of bands, can be described as
(1) 
where and are with the same size of
, which represent the clean HSI cube and the mixed noise term, respectively. Now the objective of HSI restoration is to estimate
from the observed by exploiting the structures of the clean HSI and the noise terms .Following the studies of [17, 28], we divide the noise term into two subterms as is shown in Fig.2, i.e., the Gaussian noise term and the sparse noise term including stripes, impulse noise, and dead pixels, leading to the following degradation model:
(2) 
As such, the Frobenius norm and the norm can be naturally used to model such two noise terms and respectively.
It is well known that each spectral signature can be represented by a linear combination of a small number of pure spectral endmembers, which means that the mode3 matricization can be factorized as , where is the socalled endmember matrix, and is regarded as the abundance matrix. As stated in [13], the number of endmembers is relatively small, that is, or , which means that only a small fraction of singular values are greater than zero as shown in the singular value curve of in Fig. 3. Besides, like other real images, there have certain correlations in both two spatial modes as also shown in Fig. 3, where the singular value curves of and have obvious decaying trends. Therefore, we can utilize lowrank Tucker decomposition to characterize the aforementioned spatialandspectral correlation.
TV regularization has been widely used to explore the spatial piecewise smooth structure for tackling HSI restoration task. See, e.g., [17], [25]. One can find in Fig. 4 an illustrative example of edge detection for demonstrating the spatial smoothness of one HSI. Actually, as presented in the third row of Fig. 4, there also exits local smooth structure of a HSI along with its spectral mode. That is, most of the difference value between adjacency bands in spectral domain nearly equal to 0. It is obvious that the commonly used bandbyband TV regularizer neglects such spectral smoothness, which motivates us to design a new SSTV regularizer to fully explore the piecewise smooth structure in both spatial and spectral domains. It should be also pointed that the works [39, 40] are very different from our work, although they both utilized SSTV regularizers for destriping and denoising tasks, respectively.
IiiB Lowrank Tensor Decomposition with Anisotropic SSTV
As previously stated, exploiting the prior knowledge is a key consideration for HSI mixednoise removal. Based on the discussion of the above Section II.A, by combining the lowrank and TV properties in both spatial and spectral models, we introduce a TV regularized lowrank tensor decomposition (LRTDTV for short) model, that is,
(3) 
where is the Tucker decomposition with core tensor and factor matrices of rank , and the SSTV term is defined as
where is the th entry of , and is the weight along the th mode of that controls its regularization strength. It is not hard to see that the above problem (3) is a nonconvex optimization problem due to the nonconvexity of Tucker decomposition. Thus, it would be possible to find a good local solution by using the popular augmented Lagrange multiplier (ALM) method [41] that we shall show in the next subsection.
It is worthy noting that the proposed model can fully capture the spatial and spectral information of the HSI, and thus is expected to have a strong ability of mixed noise removal. Specifically, the Tucker decomposition could make use of the spectral similarity of all the pixels and the certain correlations in both two spatial modes. Once the sparse noise term, including impulse noise, dead lines and stripes, is detected by the norm term, the designed SSTV norm term is used to characterize the piecewise smooth structure in both spatial and spectral domains and, as a result, help to remove the Gaussian noise. The Frobenius norm term is further used to model Gaussian noise and can enhance the performance of the proposed model in some heavy Gaussian noise situations.
IiiC Optimization Procedure
By introducing some auxiliary variables , we rewrite (3) as the following equivalent minimization problem:
(4) 
where is the socalled weighted threedimensional difference operator and , , are the firstorder difference operators respect to three different directions of HSI cube. Based on the ALM methodology, the above problem (4) can be transformed into minimizing the following augmented Lagrangian function:
(5) 
under the constraints , and , where is the penalty parameter, and are the Lagrange multipliers. Therefore, we can alternative optimize the augmented Lagrangian function (5) over one variable, while fixing the others. Specifically, in the th iteration, variables involved in the model (3) can be updated as follows:
1) Update . Extracting all terms containing from the augmented Lagrangian function (5), we need to solve:
which can be easily transformed into the following equivalent problem:
By using the classic HOOI algorithm [38], we can easily get the , and . Then can be updated as follows:
(6) 
2) Update . Extracting all terms containing from the augmented Lagrangian function (5), we can deduced:
Thus, optimizing this problem can be treated as solving the following linear system:
where indicates the adjoint operator of . Thanks to the blockcirculant structure of the matrix corresponding to the operator , it can be diagonalized by the 3D FFT matrix. Therefore, we have
(7) 
where fftn and ifftn indicate fast Fourier transform and its inverse transform respectively, is the elementswise square, and the division is also performed elementwisely.
3) Update . Extracting all terms containing from (5), we can get
By introducing the socalled softthresholding operator:
where and , then we can update as
(8) 
4) Update . Similarly, we should consider
By using the softthresholding operator introduced before, then the solution of above problem can be expressed as
(9) 
5) Update . Extracting all items related to from (5), it can be easily deduced:
Thus through a simple calculation, one can get its solution as follows:
(10) 
6) Updating Multipliers. Based on the ALM algorithm, the multipliers are updated by the following equations:
(11) 
Summarizing the aforementioned discussion, we now arrive at obtaining an ALM method to solve the proposed LRTDTV model (3), as presented in Algorithm 1.
In the LRTDTV solver, the inputs are the noisy image , desired rank for Tucker decomposition, the stopping criteria , and the regularized parameters , , and . Considering the fact that these three parameters have certain proportional relationship, we simply set and then tune and . More details and discussions would be presented in Section IV.C. For another important parameter , we first initialize it as and then update it as in each iteration. This strategy of determining the variable has been widely used in the ALMbased methods [41, 42], which can facilitate the convergence of the algorithm.
IiiD The LRTDTV Approximate Model
In the LRTDTV model (3
), we separate the mixture noise into two parts namely the sparse noise and the Gaussian noise. Owing to the TV regularization has certain effect to remove the Gaussian noise, when the the variance of Gaussian noise is not heavy, we can simplify (
3) as follows:(12) 
Similarly, we can derive an ALMbased optimization procedure. For simplicity, we directly give the closedform solution of each subproblem in the following.
Solving subproblem:
We need to consider the following optimization problem:
By using HOOI algorithm, we can get the estimation of the core tensor , and . Then can be updated as follows:
(13) 
Solving subproblem:
(14) 
Solving subproblem:
(15) 
Solving subproblem:
(16) 
The multipliers are updated by the following equations:
(17) 
Similarly, we could get an efficient ALM method to solve the LRTDTV approximate model (12), that is, sequentially update (1316) until certain convergence condition is satisfied. Note that all the parameters used for updating (1316) can be set as the same values that used in Algorithm 1.
Iv Experimental Results and Discussion
Both simulated and real image data experiments were carried out to demonstrate the effectiveness of the proposed LRTDTV method for HSI restoration. To thoroughly evaluate the performance of LRTDTV, we implemented seven different popular HSI restoration methods for comparison, i.e., the nuclear norm minimazition (NNM) based LRMA [23], the weighted nuclear norm minimazition (WNNM) based LRMA [26], the GoDec based lowrank matrix recovery (LRMR) [13], the weighted schatten norm minimazition (WSNM) based LRMA[28], the TVregularized lowrank matrix factorization (LRTV)[17], the block matching 4D filtering (BM4D) [43], and the decomposable nonlocal tensor dictionary learning (TDL) [33]. These methods represent the stateoftheart HSI restoration methods, especially WSNM and LRTV, and their implementation codes can be directly obtained from the authors’ websites.
In all the following experiments, the parameters in those compared methods were manually adjusted according to their default strategies. While for our LRTDTV solver, we would like to present a detailed discussion of parameter selection in Section IV.C. In addition, to facilitate the numerical calculation and visualization, all the bands of the HSI are normalized into and then will be stretched to original level after restoration.
Iva Simulated Data Experiments
The synthetic data were firstly considered by [13], which were generated using the ground truth of the Indian Pines dataset [44]. The size of the synthetic HSI was , and the reflectance values of all the voxels in the HSI were linearly mapped to .
To simulate noisy HSI data, we added several types of noise to the original synthetic HSI data under six different cases to test the performance of all the compared restoration methods, both in visual quality and quantitative perspective. We list the details in the following:
Case 1): For different bands, the noise intensity was equal. And the same distribution of zeromean Gaussian noise was added to each band. The variances of the Gaussian noise was 0.1.
Case 2): In this case, the Gaussian noise was added to each band just as in Case 1). In addition, some deadlines were added from band 91 to band 130, with the number of stripes was randomly selected from 3 to 10, and the width of the stripes was randomly generated from 1 to 3.
Case 3): In this case, the noise intensity was equal for different bands. And the same distribution of zeromean Gaussian noise and the same percentage of impulse noise were added to each band. The variance of the Gaussian noise was 0.075, and the percentage of the impulse noise was 0.15.
Case 4): In this case, the Gaussian and impulse noise were added just like in Case 3). Besides, the deadlines were added from band 91 to band 130, with the number of stripes was randomly selected from 3 to 10, and the width of the deadlines was randomly generated from 1 to 3.
Case 5): In this case, the noise intensity was different for different bands. That is, different variance zeromean Gaussian noise was added to each band, with the variance value being randomly selected from 0 to 0.2, and different percentages of impulse noise were added, which were randomly selected from 0 to 0.2. In addition, the deadlines noise were added to some band just as in Case 4).
Case 6): In this case, the Gaussian noise, impulse noise, and deadlines were added just as in Case 5). In addition, some stripes were added from band 161 to band 190, with the number of stripes being randomly selected from 20 to 40.
IvA1 Visual quality comparison
In terms of visual quality, two representative bands of restored HSIs in two typical cases obtained by different methods are shown in Fig. 5 and Fig. 6, respectively. More precisely, Fig. 5 shows the restoration results of band 36, which is corrupted by both the Gaussian noise and impulse noise in Case 3). And Fig. 6 shows the restoration results of band 116, which is severely corrupted by a mixture of Gaussian noise, impulse noise and deadline in Case 5).
It is evident that, for both cases, the original clean HSIs have undergone great changes as shown in Fig. 5(b) and Fig. 6(b) respectively. Aiming at better visual comparison, we marked a same subregion of each subfigure in Figs. 5 and 6 by a green box and then enlarged it in a red box. Several observations can be easily made from both Figs. 5 and 6. Firstly, all the compared methods, to some extent, could remove such mixednoise. Secondly, the proposed LRTDTV method performs best among all the compared methods, effectively removing the noise while preserving the essential structures of original HSIs. For instance, compared to red dashed boxes in Figs. 6(ci), the red dashed box in Fig. 6(j) produced by our LRTDTV method preserves the edge clearer and sharper. Thirdly, the two TVregularized methods, i.e., LRTV and LRTDTV, perform much better than other competing methods, and our proposed LRTDTV outperforms LRTV to a certain extent. This demonstrates the power of both TV regularization and directly low rank tensor modeling techniques for HSI restoration.
IvA2 Quantitative comparison
To better understand the overall performance of different methods, the mean peak signaltonoise ratio (MPSNR) and mean structural similarity index (MSSIM) and Erreur Relative Globale Adimensionnelle de Synthese (ERGAS) are adopted to objectively evaluate the restoration quality:
where and denote the PSNR and SSIM values for the th band, and denote the th band of the reference image the restoration image, respectively.
Table II documents the restoration results by all the compared methods in terms of aforementioned three indices. Because the Gaussian noise is not the dominant noise in Case 2)6), so we adopt the approximation model (15) in such cases except Case 1) where the Gaussian noise dominates other noises. It is clear that the proposed LRTDTV enjoys a superior performance over the other popular approaches. We also show in Fig 7. that the PSNR and SSIM values of each band in all the simulated data cases. It is easy to see that, the TVregularized methods including LRTV and the proposed LRTDTV can get much higher SSIM and PSNR values than other ones for almost every band. In addition, except that the SSIM values of LRTV are slightly higher than those of LRTDTV in Case 1), LRTDTV achieves the best performance among all the methods in terms of both SSIM and PSNR .


Noise case  Evaluation index  Noise  NNM  WNNM  LRMR  BM4D  TDL  WSNM  LRTV  LRTDTV 


Case 1)  MPSNR(dB)  19.99  30.97  32.58  36.20  38.44  38.05  37.32  38.68  40.76 
MSSIM  0.3672  0.8727  0.8420  0.9311  0.9763  0.9674  0.9453  0.9853  0.9804  
ERGAS  233.99  69.02  57.65  36.85  29.04  30.72  32.37  28.47  23.02  


Case 2)  MPSNR(dB)  19.34  30.81  32.46  35.67  35.10  33.22  36.12  38.04  40.54 
MSSIM  0.3592  0.8715  0.8415  0.9291  0.9391  0.8743  0.9402  0.9818  0.9895  
ERGAS  257.88  70.06  58.39  39.84  110.40  112.35  44.89  49.18  23.44  


Case 3)  MPSNR(dB)  13.07  31.61  32.36  36.40  28.59  27.50  38.14  39.54  41.08 
MSSIM  0.1778  0.8889  0.8786  0.9345  0.8401  0.9076  0.9551  0.9866  0.9910  
ERGAS  520.53  64.36  59.71  36.04  90.29  102.08  29.52  35.22  21.98  


Case 4)  MPSNR(dB)  12.92  31.45  32.29  35.76  27.02  26.54  36.63  38.75  40.72 
MSSIM  0.1748  0.8878  0.8777  0.9316  0.8073  0.8365  0.9485  0.9826  0.9906  
ERGAS  529.82  65.45  60.10  39.99  128.11  120.54  46.32  55.44  22.90  


Case 5)  MPSNR(dB)  13.80  29.62  31.33  33.72  27.95  24.34  35.02  36.54  38.83 
MSSIM  0.2038  0.8633  0.8445  0.8951  0.8201  0.5931  0.9285  0.9742  0.9859  
ERGAS  500.68  81.16  66.39  50.16  121.95  158.75  51.33  72.52  28.66  


Case 6)  MPSNR(dB)  13.73  29.54  29.97  33.42  27.53  23.34  33.88  36.35  38.63 
MSSIM  0.2022  0.8612  0.8431  0.8918  0.8060  0.5583  0.9261  0.9736  0.9852  
ERGAS  504.37  82.18  82.48  52.62  127.28  174.28  53.29  72.05  29.82  

To further compare the performances of all the restoration methods, it would be necessary to show the spectral signatures before and after restoration. As such, Fig. 8 shows one simple example namely the spectral signatures of pixel (20, 20) in case 5). Combining with the ERGAS values in Table II, it can be clearly seen that the proposed LRTDTV method produces the best spectral signature among all the compared methods.
IvB Real Data Experiments
Two realworld HSI data sets were used in our experiments, i.e., the Hyperspectral Digital Imagery Collection Experiment (HYDICE) urban data set [45] and the Airborne Visible/Infrared Imaging Spectrometer (AVIRIS) Indian Pines data set [44]. Before conducting the restoration process, the gray values of each HSI band were normalized to [0, 1]. We have implemented the proposed two models (3) and (12), and indeed found that they performed similar in these real data experiments. Therefore, we only present the restoration results of the approximation model (12) for visual comparison.
1) HYDICE Urban Data Set: The size of original image is . As shown in Fig. 9, the full urban image is polluted by stripes, deadlines, the atmosphere, water absorption, and other unknown noise. In this experiment, NNM, WNNM, BM4D, and TDL were implemented using the default values of parameters. For other competing methods, we manually adjusted their parameters to give best visual results. For instance, the two parameters of LRMR, i.e., and , were set as 3 and 9500, respectively.
In this set of experiments, we observed that NNM and WNNM have similar performance, while WNNM performed slightly better than NNM. Therefore, for simplicity, we don’t show the visual restoration results of NNM in the following.
Figs. 10 and 11 show the restorations of band 109 and band 207, respectively. It is not hard to see that several lowrank matrix methods including NNM, WNNM, LRMR and WSNM cannot effectively remove the stripes. In addition, some structure of objects are not clearly differentiated. This is mainly because of the fact that stripes and deadlines exist in the same place from band 104 to 110, and exist in the same place from band 199 to 210. That is, in the low rank and sparse decomposition, the stripes are more likely regarded as the lowrank part, which is assumed to be the clean image. As for BM4D, though it is a tensor filtering based method, it mainly assumes that the noise is subject to Gaussian noise, so once the HSIs suffer from heavy sparse noise, or some structure noise including stripes, this method cannot perform well. As for TDL, it loses its effectiveness when HSIs suffer from heavy mixture corrupted noise since it was also designed for removing Gaussian noise. Similar to before, since LRTV uses the spatial TV regularization to exploit the spatial information, it can remove this kind of stripes and preserve the edge information to a certain extent. By combining the spectral TV and spatial TV into a unified SSTV regularization, and utilizing the lowrank Tucker decomposition to encode the global spectralandspatial correlation, our proposed LRTDTV method can better remove the complex mixednoise and preserve the spatial texture information compared with other competing methods.
We then show in Fig. 12 that the horizontal mean profiles of band 109 before and after restoration. As shown in Fig 12(a), due to the existence of mixednoise, especially stripes and deadlines, there are rapid fluctuations in the curve before processing the restoration. After restoration, the fluctuations are more or less suppressed by all the methods. Here, one can see that the curve of proposed LRTDTV method is most stable, which is in accordance with the visual results presented in Fig. 10.
2) AVIRIS Indian Pines Data Set: This data set was acquired by the NASA AVIRIS instrument over the Indian Pines test site in Northwestern Indiana in 1992, and it has pixels and bands. The full urban image showed in Fig. 13 is mainly corrupted by deadlines, atmosphere, water absorption, and others. Similar to HYDICE data experiment, we implemented WNNM, BM4D and TDL using the default parameters, and the rest using the finetuned parameters.
We selected two typical bands, namely band 108 and band 220, to present the performance of all the compared methods in Figs. 14 and 15. Since the TDL method completely failed to restore these two band, we excluded it in this visual comparison. Here, it is easy to observe that lowrank matrix methods including NNM, WNNM, LRMR and WSNM can more or less remove some noise, but when the noise is heavy, such methods lose their utility and even occur some degradation of the gray value. See Fig. 14(f) for an example of WSNM. As for BM4D, it would oversmooth the image and distort the structure as shown in Fig. 14(e), and thus fails to give satisfactory results. Compared to the aforementioned nonTVregularized methods, LRTV and the proposed LRTDTV can remove lots of noise, but LRTV cannot preserve edges and local detail information as well as LRTDTV. More detailed visual comparison results can be seen in such red boxes. In a word, for this data set, the proposed LRTDTV can still get best performance for removing such heavy mixednoise.
We also present the vertical mean profiles of band 108 before and after restoration in Fig. 16. Similar to Fig. 12, it can again be clearly observed that LRTDTV gives the best curve among all the restored vertical mean profiles.
To sum up, extensive experiments on both simulated and real data demonstrate the clear superiority of LRTDTV over other popular methods, because more useful structures are exploited.
IvC Discussion
Basically, we introduced two LRTDTV models, one is the general model (3), the other is the approximate model (12). When the noise HSI was corrupted by very heavy Gaussian noise, the general model (3) would perform better than the approximate model (12) because of the use of the additional Frobenius norm term for modeling the Gaussian noise. In fact, the noise of real HSI is usually a mixture of several types of noise, and compared to other noise terms, the extent of Gaussian noise was not severely. Thus, considering the fact that the TV regularization has the ability of removing Gaussian noise, we would like to use the approximate model (12) for simplicity.
In the LRTDTV model, there exist several parameters need to be carefully identified. Specifically, for the Frobenius term regularization parameter , we set it as the reciprocal to the variance of Gaussian noise in general model, just as WNNM did; For the TV regularization parameter in two models, similar to many other works, it can be fixed as constant 1. Then, we need to consider how to select the term regularization parameter , the rank s among three dimensions of HSI cube, and the weights of SSTV. In the following, we shall provide the reasons for choosing such parameters in our experiments. Additionally, we present an empirical analysis of convergence of our proposed LRTDTV solver. All the results were based on the simulated data experiment in Case 5).
1) Sensitivity Analysis of Parameter : It is easy to see that is the parameter used to restrict the sparsity of the sparse noise. As stated in the RPCA model[23], the sparsity regularization parameter was set to , which was good enough to guarantee the existence of an optimal solution. Owing to bring the new TV penalty, our model is different from RPCA, we thus set as where is a tuning parameter. Fig. 17 shows the restoration results as varied in the set . It can be clearly observed from this figure that the results of the LRTV solver are relatively stable in terms of both MPSNR and MSSIM values, with the value of was changed from 10 to 25. Therefore, we suggest the use of in all the simulated data experiments.
2) Effectiveness of the Rank Constraint: In the LRTDTV optimization, we adopted the Tucker decomposition to encode the lowrank prior. So we should give the estimated ranks along the three modes before running the algorihm. For two spatial modes, we empirical adopted the eighty percent of size to ultilize the lowrank property adequately. Fig. 18 presents the MPSNR and MSSIM values of the LRTDTV solver with different rankconstrained values of spectral mode. It can be easily observed that the MPSNR and MSSIM values first increase and then decrease with the growth of the estimated rank value. This inspires us to use HSI subspace estimation method (e.g., HySime [46]) to estimated the range of values for the rank, and then choose the best desired rank value among candidates for spectral mode. As such, in all our simulation experiments, the value of spectral rank is set to 10.
3) Effectiveness of the weight of SSTV: For hyperspectral image, the spectral characteristics is a very important feature, which can be encoded by spectral TV prior. This prior can help us to remove the noise, especially for some structure noise whose distribution was different on adjacent bands as observed. Then, based on the statistical analysis of the TV value along three modes shown in Fig. 4, the weight of spatial TV is default as 1, and we only change the weight of SSTV along the spectral mode from 0 to 1. The influence of the weight of spectral TV is presented in Fig. 19, which clearly shows the benefit of tuning the weight of SSTV.
4) Convergence of the LRTDTV Solver: Fig. 20 presents the MPSNR and MSSIM gains versus the iteration number of the LRTDTV solver. Here, we can observe that, as the number of iteration increasing to a relatively large value, the relative changes of MPSNR and MSSIM converge to zero. This clearly illustrates the convergent behavior of the proposed method, which further encourages us to use it for more practical situations.
V Conclusion
In this paper, we have proposed a novel tensorbased approach for removing mixed noise in HSIs. Specifically, the low rank tensor Tucker decomposition is utilized to describe the global spatialandspectral correlation among all HSI bands, and a spatialspectral total variation (SSTV) regularization is applied to characterize the piecewise smooth structure in both spatial and spectral domain of HSIs. Besides the commonlyused norm to detect the sparse noise, we have also considered an additional Frobenius norm term for modeling the heavy Gaussian noise that may exist in some practical situations. Though the resulting nonconvex tensor optimization problem seems to be difficult to solve, we have designed an efficient algorithm based on the augmented Lagrange multiplier method. A series of simulated and real data experiments have been conducted to demonstrate the superior performance of the proposed method over some popular methods in terms of both the quantitative evaluation and the visual comparison.
In the future, we are interested in conducting the following studies. Firstly, incorporate the noise modeling idea of [47]
into our SSTV regularized low rank tensor decomposition framework to further enhance its capability for removing more complex noise in some realword scenarios. Secondly, extend the deep learning ideas of
[48, 49] to design a deep tensor architecture to learn the multilinear structure of clean HSI and identify the noise structures of observered HSI. Finally, design a new SSTV regularizer to exploit the gradient similarity among all HSI bands to further characterize the smooth structure, and as a result, improve the noise removal capability of the proposed procedure.Acknowledgment
This work was supported by the National Natural Science Foundation of China (Grant Nos. 11501440, 61603292, 61673015, 61373114).
References
 [1] A. F.H. Goetz, “Three decades of hyperspectral remote sensing of the earth: a personal view,” Remote Sensing of Environment, vol. 113, pp. S5–S16, 2009.
 [2] R. Willett, M. Duarte, M. Davenport, and et al, “Sparsity and structure in hyperspectral imaging: Sensing, reconstruction, and target detection,” IEEE Signal Process. Mag., vol. 31, no. 1, pp. 116–26, 2014.
 [3] A. A. Green, M. Berman, P. Switzer, and et al, “A transformation for ordering multispectral data in terms of image quality with implications for noise removal,” IEEE Transactions on Geoscience and Remote Sensing, vol. 26, no. 1, pp. 65–74, 1988.
 [4] M. Elad and M. Aharon, “Image denoising via sparse and redundant representations over learned dictionaries,” IEEE Transactions on Image Processing, vol. 15, no. 12, pp. 3736–3745, 2006.
 [5] K. Dabov, A. Foi, V. Katkovnik, and et al, “Image denoising by sparse 3d transformdomain collaborative filtering,” IEEE Transactions on Image Processing, vol. 16, no. 8, pp. 2080–2095, 2007.

[6]
J. Mairal, F. Bach, J. Ponce, and et al, “Nonlocal sparse models for image
restoration,” in
Proceedings of IEEE Conference on Computer Vision and Pattern Recognition
, 2009, pp. 2272–2279.  [7] H. Othman and S.E. Qian, “Noise reduction of hyperspectral imagery using hybrid spatialspectral derivativedomain wavelet shrinkage,” IEEE Transactions on Geoscience and Remote Sensing, vol. 44, no. 2, pp. 397–408, 2006.
 [8] P. Zhong and R. Wang, “Multiplespectralband crfs for denoising junk bands of hyperspectral imagery,” IEEE Transactions on Geoscience and Remote Sensing, vol. 51, no. 4, pp. 2260–2275, 2013.
 [9] Q. Yuan, L. Zhang, and H. Shen, “Hyperspectral image denoising employing a spectral–spatial adaptive total variation model,” IEEE Transactions on Geoscience and Remote Sensing, vol. 50, no. 10, pp. 3660–3677, 2012.
 [10] Y. Qian and M. Ye, “Hyperspectral imagery restoration using nonlocal spectralspatial structured sparse representation with noise estimation,” IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, vol. 6, no. 2, pp. 499–515, 2013.

[11]
G. Chen and S.E. Qian, “Denoising of hyperspectral imagery using principal component analysis and wavelet shrinkage,”
IEEE Transactions on Geoscience and Remote Sensing, vol. 49, no. 3, pp. 973–980, 2011.  [12] Y. Wang, R. Niu, and X. Yu, “Anisotropic diffusion for hyperspectral imagery enhancement,” IEEE Sensors Journal, vol. 10, no. 3, pp. 469–477, 2010.
 [13] H. Zhang, W. He, L. Zhang, and et al, “Hyperspectral image restoration using lowrank matrix recovery,” IEEE Transactions on Geoscience and Remote Sensing, vol. 52, no. 8, pp. 4729–4743, 2014.
 [14] H. Song, G. Wang, and K. Zhang, “Hyperspectral image denoising via lowrank matrix recovery,” Remote sensing letters, vol. 5, no. 10, pp. 872–881, 2014.
 [15] W. He, H. Zhang, L. Zhang, and et al, “Hyperspectral image denoising via noiseadjusted iterative lowrank matrix approximation,” IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, vol. 8, no. 6, pp. 3050–3061, 2015.
 [16] M. Wang, J. Yu, J.H. Xue, and et al, “Denoising of hyperspectral images using group lowrank representation,” IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, vol. 9, no. 9, pp. 4420–4427, 2016.
 [17] W. He, H. Zhang, L. Zhang, and et al, “Totalvariationregularized lowrank matrix factorization for hyperspectral image restoration,” IEEE Transactions on Geoscience and Remote Sensing, vol. 54, no. 1, pp. 178–188, 2016.
 [18] J. Liu, P. Musialski, and P. Wonka, “Tensor completion for estimating missing values in visual data,” IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 34, no. 1, pp. 208–20, 2013.
 [19] M. Yuan and C.H. Zhang, “On tensor completion via nuclear norm minimization,” Foundations of Computational Mathematics, vol. 16, no. 4, pp. 1031–68, 2016.
 [20] W. Cao, Y. Wang, J. Sun, and et al., “Total variation regularized tensor rpca for background subtraction from compressive measurements,” IEEE Transactions on Image Processing, vol. 25, no. 9, pp. 4075–90, 2016.

[21]
A. Anandkumar, P. Jain, Y. Shi, and et al, “Tensor vs. matrix methods: robust
tensor decomposition under block sparse perturbations,” in
Proceedings of the 19th International Conference on Artificial Intelligence and Statistics
, 2016.  [22] C. Lu, J. Feng, Y. Chen, and et al., “Tensor robust principal component analysis: exact recovery of corrupted lowrank tensors via convex optimization,” in Proceedings of IEEE Conference on Computer Vision and Pattern Recognition, 2016.
 [23] E. Candès, X. Li, Y. Ma, and et al, “Robust principal component analysis?” Journal of the ACM, vol. 58, no. 3, pp. 1–39, 2011.

[24]
T. Zhou and D. Tao, “Godec: Randomized lowrank & sparse matrix decomposition
in noisy case,” in
Proceedings of the 28th International Conference on Machine Learning (ICML11)
, 2011, pp. 33–40.  [25] Z. Wu, Q. Wang, Z. Wu, and et al, “Total variationregularized weighted nuclear norm minimization for hyperspectral image mixed denoising,” Journal of Electronic Imaging, vol. 25, no. 1, pp. 013 037–013 037, 2016.
 [26] S. Gu, L. Zhang, W. Zuo, and et al, “Weighted nuclear norm minimization with application to image denoising,” in Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, 2014, pp. 2862–2869.
 [27] Y. Peng, J. Suo, Q. Dai, and et al, “Reweighted lowrank matrix recovery and its application in image restoration,” IEEE Transactions on Cybernetics, vol. 44, no. 12, pp. 2418–2430, 2014.
 [28] Y. Xie, Y. Qu, D. Tao, and et al, “Hyperspectral image restoration via iteratively regularized weighted schatten norm minimization,” IEEE Transactions on Geoscience and Remote Sensing, vol. 54, no. 8, pp. 4642–4659, 2016.
 [29] N. Renard, S. Bourennane, and J. BlancTalon, “Denoising and dimensionality reduction using multilinear tools for hyperspectral images,” IEEE Geoscience and Remote Sensing Letters, vol. 5, no. 2, pp. 138–142, 2008.
 [30] A. Karami, M. Yazdi, and A. Z. Asli, “Noise reduction of hyperspectral images using kernel nonnegative tucker decomposition,” IEEE Journal of Selected Topics in Signal Processing, vol. 5, no. 3, pp. 487–493, 2011.
 [31] X. Liu, S. Bourennane, and C. Fossati, “Denoising of hyperspectral images using the parafac model and statistical performance analysis,” IEEE Transactions on Geoscience and Remote Sensing, vol. 50, no. 10, pp. 3717–3724, 2012.
 [32] X. Guo, X. Huang, L. Zhang, and et al, “Hyperspectral image noise reduction based on rank1 tensor decomposition,” ISPRS journal of photogrammetry and remote sensing, vol. 83, pp. 50–63, 2013.
 [33] Y. Peng, D. Meng, Z. Xu, and et al, “Decomposable nonlocal tensor dictionary learning for multispectral image denoising,” in Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, 2014, pp. 2949–2956.
 [34] Q. Xie, Q. Zhao, D. Meng, and et al, “Multispectral images denoising by intrinsic tensor sparsity regularization,” in Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, 2016, pp. 1692–1700.
 [35] C. Li, Y. Ma, J. Huang, and et al, “Hyperspectral image denoising using the robust lowrank tensor recovery.” Journal of the Optical Society of America. A, Optics, image science, and vision, vol. 32, no. 9, pp. 1604–1612, 2015.
 [36] Z. Wu, Q. Wang, J. Jin, and et al, “Structure tensor total variationregularized weighted nuclear norm minimization for hyperspectral image mixed denoising,” Signal Processing, vol. 131, pp. 202–219, 2017.
 [37] S. Lefkimmiatis, A. Roussos, P. Maragos, and M. Unser, “Structure tensor total variation,” SIAM Journal on Imaging Sciences, vol. 8, no. 2, pp. 1090–1122, 2015.
 [38] T. G. Kolda and B. W. Bader, “Tensor decompositions and applications,” Siam Review, vol. 51, no. 3, pp. 455–500, 2009.
 [39] Y. Chang, L. Yan, H. Fang, and et al, “Anisotropic spectralspatial total variation model for multispectral remote sensing image destriping,” IEEE Transactions on Image Processing, vol. 24, no. 6, pp. 1852–1866, 2015.
 [40] C. Jiang, H. Zhang, L. Zhang, and et al, “Hyperspectral image denoising with a combined spatial and spectral weighted hyperspectral total variation model,” Canadian Journal of Remote Sensing, vol. 42, no. 1, pp. 53–72, 2016.
 [41] Z. Lin, M. Chen, and Y. Ma, “The augmented lagrange multiplier method for exact recovery of corrupted lowrank matrices,” arXiv preprint arXiv:1009.5055, 2010.
 [42] D. P. Bertsekas, Constrained optimization and Lagrange multiplier methods. Academic press, 2014.
 [43] M. Maggioni, V. Katkovnik, K. Egiazarian, and A. Foi, “Nonlocal transformdomain filter for volumetric data denoising and reconstruction.” IEEE Transactions on Image Processing, vol. 22, no. 1, pp. 119–33, 2012.
 [44] https://engineering.purdue.edu/~biehl/MultiSpec/hyperspectral.html, [Online; accessed 19April2017].
 [45] http://www.tec.army.mil/hypercube, [Online; accessed 19April2017].
 [46] J. M. BioucasDias and J. M. Nascimento, “Hyperspectral subspace identification,” IEEE Transactions on Geoscience and Remote Sensing, vol. 46, no. 8, pp. 2435–2445, 2008.
 [47] X. Chen, Z. Han, Y. Wang, and et al, “A general model for robust tensor factorization with unknown noise,” arXiv preprint arXiv:1705.06755, 2017.
 [48] K. Gregor and Y. LeCun, “Learning fast approximations of sparse coding,” in Proceedings of the 27th International Conference on Machine Learning (ICML10), 2010, pp. 399–406.
 [49] P. Sprechmann, A. M. Bronstein, and G. Sapiro, “Learning efficient sparse and low rank models,” IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 37, no. 9, pp. 1821–1833, 2015.
Comments
There are no comments yet.