Noise Reduction in Hyperspectral Imagery : Overview and Application

Hyperspectral remote sensing is based on measuring the scattered and reflected electromagnetic signals from the Earth’s surface emitted by the Sun. The received radiance at the sensor is usually degraded by atmospheric effects and instrumental (sensor) noises which include thermal (Johnson) noise, quantization noise, and shot (photon) noise. Noise reduction is often considered as a preprocessing step for hyperspectral imagery. In the past decade, hyperspectral noise reduction techniques have evolved substantially from two dimensional bandwise techniques to three dimensional ones, and varieties of low-rank methods have been forwarded to improve the signal to noise ratio of the observed data. Despite all the developments and advances, there is a lack of a comprehensive overview of these techniques and their impact on hyperspectral imagery applications. In this paper, we address the following two main issues; (1) Providing an overview of the techniques developed in the past decade for hyperspectral image noise reduction; (2) Discussing the performance of these techniques by applying them as a preprocessing step to improve a hyperspectral image analysis task, i.e., classification. Additionally, this paper discusses about the hyperspectral image modeling and denoising challenges. Furthermore, different noise types that exist in hyperspectral images have been described. The denoising experiments have confirmed the advantages of the use of low-rank denoising techniques compared to the other denoising techniques in terms of signal to noise ratio and spectral angle distance. In the classification experiments, classification accuracies have improved when denoising techniques have been applied as a preprocessing step.


Introduction
Remote sensing has been substantially influenced by hyperspectral imaging in the past decades [1].Hyperspectral cameras provide contiguous electromagnetic spectra ranging from visible over near-infrared to shortwave infrared spectral bands (from 0.3 µm to 2.5 µm).The spectral signature is the consequence of molecular absorption and particle scattering, allowing to distinguish between materials with different characteristics.Hyperspectral remote sensing applications include agriculture, Remote Sens. 2018, 3, 482 2 of 28 environmental monitoring, weather prediction, military [2], food industry [3], biomedical [4], and forensic research [5].
A hyperspectral image (HSI) is a three dimensional (3D) datacube in which the first two dimensions represent spatial information and the third dimension represents the spectral information of a scene.Figure 1 shows an illustration of a hyperspectral datacube.Hyperspectral spaceborne sensors capture data in several narrow spectral bands, instead of a single wide spectral band.In this way, hyperspectral sensors can provide detailed spectral information from the scene.However, since the width of spectral bands significantly decreases, the received signal by the sensor also decreases.This leads to a trade-off between spatial resolution and spectral resolution.Therefore, to improve the spatial resolution of hyperspectral images, airborne imagery has been widely used.Further information about the different types of hyperspectral sensors is given in [2].In this review, we focus on the hyperspectral cameras which provide the reflectance from a scanned scene.In real word HSI applications, the observed HSI is degraded by different sources, related to the applied imaging technology, system, environment, etc., and therefore, the noise free HSI needs to be estimated.When the observed signal is degraded by noise sources, the estimation task is referred to as "denoising".
The received radiance at the remote sensing hyperspectral camera is degraded by atmospheric effects and instrumental noises.The atmospheric effects should be compensated to provide the reflectance.Instrumental (sensor) noise includes thermal (Johnson) noise, quantization noise and shot (photon) noise which cause corruption in the spectral bands by varying degrees.These corrupted bands degrade the efficiency of the HSI analysis techniques and therefore they are often removed from the data before any further processing.Alternatively, HSI denoising can be considered as a preprocessing step in HSI analysis to improve the signal to noise ratio (SNR) of HSI.
Figure 2 illustrates the dynamics of the important subject of hyperspectral image denoising in the hyperspectral community.The reported numbers include both scientific journal and conference papers on this particular subject using "hyperspectral" and "(denoising, restoration, or noise reduction)" as the main keywords used in the abstracts.In order to highlight the increase in this number, the time period has been split into a number of equal time slots (i.e., 1998-2001, 2002-2005, 2006-2009, 2010-2013, 2014-2017 (October 1st)).The exponential growth in the number of papers reveals the popularity of this subject.In this review, we have two main objectives: (1) giving an overview of the denoising techniques which have been developed for HSI and compare their performance in terms of improving the SNR; (2) demonstrating the effect of these techniques when applied as preprocessing for classification applications.

Notation
In this paper, the number of bands and pixels in each band are denoted by p and n = (n 1 × n 2 ), respectively.Matrices are denoted by bold and capital letters, column vectors by bold letters, the element placed in the ith row and jth column of matrix X by x ij , the jth row by x T j , and the ith column by x (i) .The identity matrix of size p × p is denoted by I p .X stands for the estimate of the variable X.The Frobenius norm and the Kronecker product are denoted by .F and ⊗, respectively.Operator vec vectorizes a matrix and vec −1 in the corresponding inverse operator.tr(X) denotes the trace of matrix X.In Table 1, the different symbols and their definition are given.
Table 1.The different symbols and their definition.

Symbols Definition
x i the ith entry of the vector x x ij the (i, j)th entry of the matrix X x (i) the ith column of the matrix X x T j the jth row of the matrix X x 0 l 0 -norm of the vector x, i.e., the number of nonzero entries.x 1 l 1 -norm of the vector x, obtained by ∑ i |x i |.
x 2 l 2 -norm of the vector x, obtained by ∑ i x 2 i .X 0 l 0 -norm of the matrix X, i.e., the number of nonzero entries.
X 1 l 1 -norm of the matrix X, obtained by ∑ i,j x ij .

X F
Frobenius-norm of the matrix X, obtained by ∑ i,j x 2 ij .X * Nuclear-norm of the matrix X, obtained by ∑ i σ i (X), i.e., the sum of the singular values.X the estimate of the variable X. tr(X) the trace of the matrix X.

Dataset Description
The datasets that are used in this paper are described below.

Houston
This hyperspectral dataset was captured by the Compact Airborne Spectrographic Imager (CASI) over the University of Houston campus and the neighboring urban area in June 2012.This dataset is composed of 144 bands ranging from 0.38 µm to 1.05 µm and the spatial resolution is 2.5 m.The image contains 349 × 1905 pixels.The available groundtruth covers 15 classes of interest.Table 2 gives information on all 15 classes including the number of training and test samples.Figure 3 shows a three-band false color composite image and its corresponding and already-separated training and test samples.

III. DIM
The increasing sp benefits precision pa the memory capacity conventional signal p spatial dimension of per-pixel, the data vo bands.The data volu series hyperspectral mental changes.The the data will easily e personal computers.M ratio between the spe samples is high, hig from the well-known Dimensionality reduc eliminating statistical keeping as much spe used in hyperspectral can represent most of

Trento Data
The second dataset was captured over a rural area in the south of the city of Trento, Italy.This dataset is of size 600 by 166 pixels with a spatial resolution of 1 m with 63 spectral bands (ranging from 402.89 to 989.09 nm) captured by the AISA Eagle sensor.The available groundtruth covers six classes of interest including Building, Woods, Apple trees, Roads, Vineyard, and Ground.Figure 4 illustrates a false color composite representation of the hyperspectral data and the corresponding training and test samples.Table 3 provides information about the corresponding number of training and test samples.The third dataset was taken by the Airborne Visible/Infrared Imaging Spectrometer (AVIRIS) sensor over an agricultural area located at northwestern Indiana.This dataset is composed of 220 bands ranging from 400 nm to 2500 nm.The image contains 145 × 145 pixels with a spatial resolution of 20 m per pixel.The available groundtruth covers 16 classes of interest.Table 4 provides detailed information about all 16 classes and the corresponding number of training and test samples.Figure 5 presents a three-band false color composite image and its corresponding training and test samples.The last dataset is the Washington DC mall which is an airborne dataset captured over the Washington DC Mall in August 1995 using the HYDICE sensor (Hyperspectral Digital Imagery Collection Experiment).The sensor provides 210 bands in the 0.4-2.4µm spectral region where each band contains 1280 lines with 307 pixels.After removing noisy bands, the available dataset contains 191 bands.The reference ground truth contains 7 classes of interest, given in Table 5. Figure 6 shows the test and training samples, and a false color image of the dataset using bands 60, 27 and 17.

Hyperspectral Modeling
A hyperspectral image can be represented by 1 dimensional (1D), 2 dimensional (2D), or 3 dimensional (3D) models.In 1D, 2D, and 3D models the HSI is treated as a combination of spectral pixel-vectors, spectral bands and as a whole cube, respectively.In other words, in 1D and 2D models, spatial and spectral correlations are ignored, respectively.However, in 3D models, both spatial and spectral correlations are taken into account.Using the matrix form we can represent the observed degraded HSI as a combination of a true unknown signal and additive noise: where H is an n × p matrix containing the vectorized observed image at band i in its ith column, X is the true unknown signal which needs to be estimated, and N is an n × p matrix representing the noise.
Other models may also be considered for HSI [6], however, model ( 1) is widely used in the literature.Model (1) can be generalized by ( [7]): where A is an n × n and M is a p × r matrix (1 ≤ r ≤ p).These are two dimensional and one dimensional projection matrices, respectively.W (n × r) is the unknown projected HSI.A and M are often selected to decorrelate the signal and noise in the HSI spatially and spectrally, respectively, and they can be known or unknown (in HSI denoising they are usually assumed to be known).Model ( 2) is a 3D model (see Appendix A for more details), however 2D and 1D models can be obtained as special cases of model (2).If M = I, then model ( 2) becomes a 2D model.If A = I, then model ( 2) becomes a 1D model, and if A = M = I, then model ( 2) is equivalent to model (1) which is also a 1D model.Assuming model (2) and A and M known, the HSI denoising task is to estimate W, and the HSI is restored by X = A ŴM T .
As an example, consider a 3D model obtained by using a 3D wavelet basis: where where the signal is only projected spatially (i.e., 2D wavelet applied on each band separately), and if D 2 = I, a 1D wavelet model is obtained: where the signal is only projected spectrally (i.e., 1D wavelet applied on each spectral pixel separately).
Another example is the 1D model that is widely used for spectral unmixing: which is a special case of model ( 2) with A = I and M = E.In (6) the HSI is projected spectrally by E, a matrix of endmembers, and W contains the abundance maps.

Hyperspectral Denoising
Assuming model ( 1), the denoising task is to estimate the original (unknown) signal X.This can be done by penalized least squares optimization: where the first term is the fidelity term, φ(X) is the penalty term, and λ determines the tradeoff between both terms.Equivalently, model (1) can be solved by solving the constrained minimization problem: where is a small number to relax the exact solution to the problem: Remote Sens. 2018, 3, 482 9 of 28 Note that the penalty method turns the constrained minimization problem (8) into the unconstrained minimization problem (7).Another constrained formulation of the denoising problem is to minimize the fidelity term subject to some constraint on the penalty term: where K is the upper bound of the constraint.The penalty term is usually selected based on the chosen model, the prior knowledge, and characteristics of the data.For instance, if we use sparsifying bases in model (2) for A and M, such as wavelet bases (model ( 3)) then it is better to use penalties which promote sparsity such as 1 (or 0 ) in (7): Note that function φ can also be a combination of multiple penalty terms.

HSI Denoising Challenges
HSI denoising is a delicate task and needs specific attention compared to denoising of other images due to the importance of preserving spectral information.The high spectral correlation provides a great advantage for the denoising task, however, oversmoothing causes loss of valuable spectral information.In the next sections, we point out the main challenges related to the development of HSI denoising algorithms.

Hyperspectral Model and Parameter Selection
The selection of the HSI model requires attention.For instance, in model ( 2), A and M need to be carefully selected for an improved performance of HSI denoising.Additionally, the selection of the model parameters is not a trivial task.For instance, over-or under-estimating λ in (7) yields either information loss or poor denoising performance.Providing ad hoc or experimental strategies for these estimations may make denoising algorithms unreliable and problematic to be used as a preprocessing step in HSI analysis.It is worth mentioning that selecting more complicated models or penalties makes the parameter selection task much harder.
In [7], a model and parameter selection criterion is given for a general model of the form (2) where A and M are orthogonal matrices, N is Gaussian noise, and the hyperspectral signal is given by X = A ŴM T , where Ŵ is given by and B = [b ij ] = A T HM. λ is the tuning parameter and B has rank r.Since the estimated signal X = A ŴM T , the performance of denoising techniques is highly dependent on the selection of those parameters.In denoising applications, it is often of interest to select the models and the corresponding parameters based on the minimization of the mean squared error (MSE), Unfortunately, in real world applications the true signal X is unknown and thus it is impossible to compute the MSE.In [8], an unbiased estimator of the MSE, called Stein's unbiased risk estimator (SURE), was derived for deterministic signals with Gaussian noise.The general form of SURE is given by where . ., σ 2 p , and σ p is the noise standard deviation in band p.In [7], a model and parameter selection criterion was proposed, based on the estimate of MSE by using SURE for X = A ŴM T where Ŵ is given by ( 12): where I is the indicator function.The main advantage of ( 15) is that model parameters can be selected based on the MSE estimator.As can be seen, ( 15) is only dependent on the observed signal (H).Therefore, (15) lets us select the model (in the form of model ( 2)) and the models parameters (r and λ) w.r.t. the minimum of the estimation of the MSE.Equation ( 15) is called hyperspectral SURE (HySURE), and is proposed in [9] in the context of spectral unmixing to determine the subspace dimension (or the number of endmembers), r, in the absence of the noise free signal X (The Matlab code online available in [10]).Additionally, a model selection criterion that is not dependent on the unknown signal (such as ( 15)) provides an instrument to compare the denoising techniques without using simulated (noisy) HSI and by only using the observed HSI itself.

Spectral Distortion and Band-Wise Normalization Issues
Spectral information in HSIs is of great importance in HSI analysis.Therefore, it is essential that HSI denoising techniques preserve the spectral information.Both signal and noise variances are varying throughout the hyperspectral bands, which makes the noise estimation and denoising task very challenging.To cope with this issue, some denoising techniques use band-wise normalization to obtain spectral bands of similar scale.Band-wise normalization causes spectral distortion and it is not recommended for HSI.One way to deal with varying signal and noise variances is to define the model parameters to be variable w.r.t. the spectral bands.For instance, in model ( 3), ( 11) can be rewritten as where the tuning parameter λ is defined to vary w.r.t. the spectral bands (the columns of W) to cope with noise power that varies between spectral bands.

Noise Variance Estimation
Denoising techniques, and particularly the model parameter selection criteria, are often highly reliable on the estimation of the noise variance.One of the most common techniques used for HSI noise parameter estimation is multiple linear regression (MLR) [11,12].MLR was proposed in [13] as a noise estimation technique which assumes that each band is a linear combination of the other bands and therefore can be estimated by using least squares estimation.The main reason of the success of MLR is the high spectral correlation of the pixels.This technique does not take into account spatial information.On the other hand, conventional noise variance estimation techniques, such as the median estimator applied on the wavelet coefficients [14] only take into account the spatial correlations.Therefore, it is of interest to investigate variance noise estimation techniques which exploit both spectral and spatial correlations.

Dominant Noise Type Investigation
Assuming mixed noise scenarios, it is of interest to investigate the dominant noise type within an HSI.Additionally, noise estimation in such scenarios is not a trivial task and therefore it is an open question if it is more efficient to estimate the mixed noises simultaneously or hierarchically.

HSI Denoising as a Preprocessing Step
Despite the considerable progress in HSI restoration techniques, they have usually not been used in HSI analysis as a preprocessing step.This might be due to several reasons such as the computational cost, efficiency, reliability, and automation of the algorithms.The main goal of the denoising-based preprocessing stage is to improve the SNR of the observed dataset.It is of interest to investigate the contribution of the various HSI restoration approaches as a preprocessing step for further HSI analysis, such as change detection, resolution enhancement, classification or unmixing.In this paper, we will address this matter for the classification application.
1.5.6.Computational Cost HSI restoration approaches need to be computationally efficient to be useful as a preprocessing step in real-world applications.Fast computing techniques such as parallel computing and GPU programing may be considered for the fast implementation of HSI restoration approaches in the future.Particularly, fast computing techniques can considerably speed up the patch-wise or pixel-wise HSI denoising techniques.

HSI Datasets
Usually in benchmark datasets, corrupted and noisy spectral bands are removed.To evaluate the performance of HSI denoising as a preprocessing technique for further HSI analysis, the access to the complete datasets may be required.

Hyperspectral Noise Assumptions
The presence of different noise sources in a HSI makes its modeling and the denoising task very challenging.Therefore, HSI denoising approaches often consider either of the following noise types or a mixture of them: 1.6.1.Signal Independent Noise Thermal noise and quantization noise in HSI are modeled by signal independent Gaussian additive noise [15,16].Usually, noise is assumed to be uncorrelated spectrally, i.e., having a diagonal noise covariance matrix [16,17].The Gaussian assumption has been broadly used in hyperspectral analysis since it considerably simplifies the analysis and the noise variance estimation.1.6.2.Signal Dependent Noise Shot (photon) noise in HSI is modeled by the Poisson distribution for which the noise variance is signal dependent.The noise variance estimation under this assumption is more challenging than in the signal independent case [18].

Sparse Noise
Impulse noises such as salt and pepper noise, missing pixels, missing lines and other outliers often exist in the acquired HSI, and are usually due to a malfunctioning of the sensor.In this review, we categorize them as sparse noise due to their sparse characteristics.Sparsity techniques or sparse and low-rank decomposition techniques are used to remove sparse noise from the signal.In [19], impulse noise is removed by using an 1 -norm for both penalty and data fidelity terms in the proposed minimization problem.

Pattern Noise
Hyperspectral imaging systems may also induce artifacts in hyperspectral images, usually referred to as pattern noise.For instance, in push-broom imaging systems, the target is scanned line by line and the image lines are acquired in different wavelengths by an area-array detector (usually, a charged coupled device (CCD)).This line by line scanning causes an artifact called striping noise which is often due to calibration errors and sensitivity variations of the detector [20].Striping noise reduction (also referred to as destriping in the literature) for push-broom scanning techniques has been widely studied in the remote sensing literature [21,22] in particular for HSI [20,[23][24][25].

HSI Denoising Overview
During the past couple of years, a considerable amount of research has been devoted to hyperspectral image denoising.Conventional denoising methods based on 2D modeling and convex optimization techniques were not efficient for HSI because these ignore the spectral information.The highly correlated spectral bands in HSI have been found very useful to improve HSI denoising.As a result, HSI denoising techniques have evolved to methods that incorporate spectral information.These HSI denoising approaches can be categorized in four main groups, which will be treated below.

3D Model-Based and 3D filtering Approaches
3D model-based HSI denoising techniques utilize model (2), where both the spatial and spectral projections using matrices A and M, respectively are applied.The projection matrices A and M are usually selected to decorrelate the signal spatially and spectrally, respectively, and for this either dictionaries or bases are used.HSI restoration techniques based on 3D filtering are categorized in this group.Those methods try to decorrelate the noise from the signal in all 3 dimensions (spatial and spectral).In [26], the discrete Fourier transform (DFT) was used to decorrelate the signal in the spectral domain, and the 2D discrete wavelet transform (2D DWT) for denoising the signal in the spatial domain.3D wavelet shrinkage was applied for multispectral image denoising in [27].In [28], 2D bivariate wavelet shrinkage (2D BiShrink) [29] was extended to 3D for the purpose of HSI denoising.3D non-local means filtering (NLM) [30] was exploited for HSI denoising in [31].In [32], 2D DWT and principal component analysis (PCA) were used to decorrelate the noise and the signal spatially and spectrally, receptively.A 3D (blockwise) nonlocal sparse denoising method [33] was presented in [34] where the minimization problem contained a group lasso penalty and a dictionary consisting of the 3D DCT and the 3D DWT.To solve the minimization problem, the accelerated proximal gradient method was used.In [35], the HSI was treated as a 3D datacube and a HSI denoising technique was proposed that uses sparse analysis regularization and the undecimated wavelet transform (UWT), where the function φ in Equation ( 7) is the 3D undecimated wavelet transform: where U 2 and U 1 are 2D and 1D UWTs.
In [17,35], the advantages of (orthogonal and undecimated) 3D wavelets over 2D ones were discussed for HSI denoising.In [36], a new 3D model was proposed for HSI, given by X = D 2 WV T , where V contains the spectral eigenvectors of H and is given by the singular value decomposition (SVD): SVD(H) = Ũ SV T .To estimate the true signal, the 1 penalized least squares optimizer was used: Additionally, SURE was used for the selection of the regularization parameter.In [7], it was shown that for the 1 penalized least squares method of (11), 3D models outperform 2D models.Another important observation that was confirmed in [7], is the advantage of using spectral eigenvectors for the spectral projection.It was demonstrated that a 1D model that projects the data on the spectral eigenvectors (i.e., in model ( 2), A = I and M contains the spectral eigenvectors) outperforms even 3D models that use wavelet bases.

Spectral and Spatial-Spectral Penalty-Based Approaches
In HSI, the spectral bands are highly correlated.Some denoising approaches have proposed penalties which exploit this high spectral correlation.These methods usually use either model (1) or model (2) where M = I and define the function φ in (7) for spectral penalties.It is worth mentioning that when using a spectral projection matrix M, the spectral bands and therefore the spectral penalties are decorrelated.
In [37], it was assumed that X = D 2 W where D 2 contains 2D wavelet bases, and a group of 2 penalties on the 2D wavelet coefficients were proposed: This method was improved in [7] for the purpose of HSI denoising by taking into account the spectral noise variance and solving the obtained minimization problem by using the alternating direction method of multipliers (ADMM): Note that ( 19) is separable and has a closed form solution while (20) is not separable and needs to be solved iteratively by using convex optimization algorithms.
To exploit the redundancy and high correlation in the spectral bands in HSI, a penalized least squares method using a first order spectral roughness penalty (FOSRP) was proposed for HSI denoising in [38]: where R p is a (p − 1) × p difference matrix given by Assuming X = D 2 W, to exploit the multiresolution analysis (MRA) property of wavelets, the following penalty function was applied on the 2D wavelet coefficients: where λ varies with the decomposition level of the wavelet coefficients, l (1 ≤ l ≤ L).It was shown that ( 21) is separable and has a closed form solution.Additionally, SURE was utilized to select the tuning parameters, yielding an automatic and fast HSI denoising technique.The advantage of the FOSRP compared to the group 2 and 1 penalties in the wavelet domain was confirmed in [7,38].
In [39], it was shown that the use of a combination of the FOSRP and the group lasso penalty: outperforms the use of each penalty solely.Total variation (TV) [40] is a widely used and efficient regularization technique for denoising in image processing.In TV denoising, the function φ in (7) represents the total variation of the signal X.In a HSI, the penalty term can account for spatial and/or spectral variations.Cubic total variation (CTV) was proposed in [41]: where D h and D v are the matrix operators for calculating the first order vertical and horizontal differences, respectively, of a vectorized image.For an image of size β determines the weight of the spectral difference w.r.t. the spatial one.CTV exploits the gradient in the spectral direction and as a consequence improves the denoising results compared to band by band TV denoising.In [42], an adaptive version of CTV was applied for preserving texture and edges simultaneously: where ω defines spatial weights on pixels (see [42] for more detail about the selection of ω).In [43], a spatial-spectral HSI denoising approach was developed where a method of spectral derivation was proposed to concentrate the noise in the low frequencies, after which noise is removed by applying the 2D DWT in the spatial domain and the 1D DWT in the spectral domain.A spatial-spectral penalty was proposed in [44] which was based on five derivatives, one along the spectral direction and the rest applied in the spatial domain for the four neighborhood pixels.

Low-Rank Model-Based Approaches
Low-rank (LR) modeling has been widely used in HSI analysis and applications such as dimensionality reduction, feature extraction, unmixing, compression etc. Due to the redundancy in the spectral bands, the LR models often assume a much lower spectral dimension than the one provided by the HSI cameras, i.e., in model (2) r p.In [7], model (11) and HySURE [9] were applied and it was shown that the low-rank model outperforms the full-rank one.
A low-rank representation technique called Tucker3 decomposition [45] was used for hyperspectral image denoising in [46].The HSI was described by a third order tensor and the rank of the decomposition was estimated by minimizing a Frobenius norm.In [47], a similar idea was exploited by applying a higher spectral reduction.In [48], a genetic algorithm (GA) was developed for choosing the rank of the Tucker3 decomposition.This work was followed by [49], in which a kernelized version (using Gaussian radial basis functions) of the Tucker3 decomposition was proposed.A multidimensional Wiener filter on a Tucker3 decomposition of HSI was proposed in [50], where the flattening of the HSI was performed by estimating the main direction corresponding to the smallest rank.
Another low-rank modeling for HSI denoising is Parallel Factor Analysis (PARAFAC) [51].In [36], the low-rank version of model ( 18) was applied, using X = D 2 WV T where V and W are low-rank matrices (i.e., r p): where the penalty is applied on the reduced matrix W of size n × r.Recently, an automatic hyperspectral restoration technique, called HyRes was proposed in [52].HyRes used (27) and HySURE to select the model parameters which led to a parameter free technique.
In [7,53], sparse reduced-rank restoration (SRRR) using both synthesis and analysis undecimated wavelets was proposed.Assuming the low-rank model X = S 2 WV T , where W and V are low rank matrices and S 2 contains 2D synthesis undecimated wavelets, the SRRR optimizer is given by: Assuming the low-rank model X = FV T , where F and V are low-rank matrices, the SRRR is given by: where U 2 contains 2D analysis undecimated wavelets.Assuming the same model, a low-rank TV regularization was proposed in [7,54]: Finally, in [12], a wavelet-based reduced-rank regression was proposed, where in the low-rank model X = D 2 WV T , both W and V are unknown variables.For the simultaneous estimation of the two unknown matrices, an orthogonality constraint was added, which led to a non-convex optimization problem:

Approaches Making the Mixed Noise Assumption
All previous methods inherently assume signal independent additive Gaussian noise as the main source of noise.Other methods take mixed noises into consideration for HSI modeling and denoising, where the HSI X in model ( 1) is assumed to be corrupted by a mixture of the different noise sources as described in Section 1.6.

Mixed Signal Dependent and Signal Independent Noises
Noise models including a mixture of signal dependent noise (N SD ) and signal independent noise (N SI ) (N = N SI + N SD ) were proposed in [11,34,55].In these models, two parameters need to be estimated which are the variances of N SI and N SD , modeled by a Gaussian and a Poisson distribution, respectively.In [34], a 3D (block-wise) non-local sparse denoising method is proposed.The minimization problem uses a group Lasso penalty and a dictionary consisting of a 3D discrete cosine transform (3D-DCT) and a 3D discrete wavelet transform (3D-DWT), and is solved by using the accelerated proximal gradient method.In [55], N SI and N SD are removed sequentially.Maximum likelihood is used to estimate the two parameters of the noise model, and MLR is used for an initial estimation of the noise.

Mixed Signal Independent and Striping Noises
This assumption is common in techniques proposed for striping noise removal.The noise model is given by N = N SI + N Str where N Str is the striping noise.N Str depends on the signal level and the position of detectors of the acquisition array in the cross-track direction (either i or j) [23].In [20] a striping noise removal method was proposed by assuming that the striping noise contains higher spatial frequencies than the surface radiance.Then the striping noise frequencies were detected and removed using low-pass filtering.A low-rank technique was proposed in [24] which uses a regularized cost function to preserve the spatial structures of subimages from spectral bands.In [23], a subspace-based approach is proposed to restore a HSI corrupted by striping noise and signal independent noise.The noise parameters were estimated by the least squares method.A moment matching (MM) technique [56,57] has been adapted for HSI striping noise removal in [25] by exploiting the spectral correlations.

Mixed Signal Independent and Sparse Noises
A widely used mixed noise assumption for HSI denoising is N = N SI + N SP where N SP is sparse noise, described in Section 1.6.As a result, model ( 1) is rewritten as where X is assumed to have a low rank.To estimate the low rank matrix X and the sparse marix N SP simultaneously, a joint minimization problem is obtained in the form: Common choices for φ 1 and φ 2 are the nuclear norm and the 1 sparsity norm [58], leading to: This mixture assumption was used in [59] where HSI was denoised by solving where the upper bound of the rank r of X and the cardinality K of N SP were assumed to be known variables.That method was improved in [60], by taking into account the changes of the noise variance throughout the spectral bands.In [61], a denoising method was presented by adding a TV penalty to the denoising criterion (34): where X HTV denotes the norm of the sum of total variations on the spectral bands: In [62], a weighted Schatten p-norm is defined: and used to induce the low-rank property on X: where the Gaussian noise matrix N SI is also estimated in the minimization problem.In [63], a patchwise approach was proposed that exploits the nonlocal similarity across patches, using (35).Recently, a low-rank and sparse restoration technique was proposed ( [64]), using a low-rank constraint applied on the spectral difference matrix: Some methods cope with the mixed sparse and Gaussian noise without enforcing the low-rank property.The denoising method proposed in [65] utilizes a TV penalty, called spatio-spectral total variation (SSTV): where A TV-based method was proposed in [66], leading to the following minimization problem: where X CrTV was given by and ω defines spatial weights on pixels.For more detail regarding the selection of ω, we refer to [66].

Comparison of HSI Denoising Techniques
In this section, different HSI denoising methods are compared qualitatively and quantitatively on a simulated and a real dataset.
All the results in this section are means over 10 experiments (adding random Gaussian noise) and the error bars show the standard deviations.Wavelab Fast (A fast wavelet toolbox developed for HSI analysis) [67] was used for the implementation of the wavelet transforms.For all the experiments performed in this paper, Daubechies wavelets were used with 2 and 10 coefficients for spectral and spatial bases, respectively.Five decomposition levels were used for the filter banks.The Matlab codes for FORPDN and HyRes are online available in [68,69], respectively.
To evaluate the restoration results for the simulated dataset, SNR and MSAD (Mean Spectral Angle Distance) are used.The output SNR, in decibels is given by: while the noise input level for the whole cube is given by: MSAD, in degrees is given by: Figure 7a shows the comparison of the HSI denoising techniques based on SNR.The figure shows the level of the denoised hyperspectral signal SNR out compared to the level of the noise added (in dB) to the simulated dataset SNR in .The results are shown for varying SNR in from 5 to 45 dB with increments of 5 dB.The blue line shows the original noise levels, the performance of the HSI denoising methods is compared based on the gain obtained w.r.t. the original noise levels.As can be seen, 2D-Wavelet generates the lowest SNR out and for SNR in ≥ 35 dB, the gain is close to negligible.3D-Wavelet outperforms the 2D version.One can further notice that FORPDN ouperforms 3D-Wavelet consistently for all SNR in and outperforms LRMR when SNR in ≤ 15 dB, and SNR in ≥ 35 dB, while when 15 dB ≤ SNR in ≤ 35 dB, LRMR and FORPDN perform similarly.The results also show that HyRes and NAILRMA, both low-rank denoising techniques, considerably outperform all the other methods.They are both designed to cope with varying noise power throughout the spectral bands.HyRes generates the highest SNR out .
Since the spectral information is higly valuable in HSI analysis, we also use MSAD to compare the spectral preservation of the HSI denoising techniques.Figure 7b plots the MSAD of the different HSI restoration techniques w.r.t. the input noise power.The results are shown in logarithmic scale for a better visual representation.It can be seen that 2D-Wavelet produces the highest MSAD.In this experiment, 3D-Wavelet improves over LRMR when SNR in ≤ 15 dB, and SNR in ≥ 35 dB, while LRMR and 3D-Wavelet perform similarly when 15 dB ≤ SNR in ≤ 35 dB.The results also show that FORPDN outperforms 2D-Wavelet, 3D-Wavelet and LRMR consistently for all SNR in .It can also be seen that, also based on MSAD, HyRes and NAILRMA outperform the other techniques.A highly corrupted band from the simulated dataset was selected for a visual comparison of the restoration methods in Figure 8. 2D-Wavelet shows a very poor performance, which is not surprising, since denoising is applied on each spectral band individually and the information is highly corrupted in that specific band.3D-Wavelet considerably improves the visual quality, due to the incorporation of the information from the other bands by 3D modeling and filtering.FORPDN, NAILRMA and HyRes all perform very well.The weak performance of LRMR on this band is due to the fact that it is not designed to cope with the variation of the noise variance throughout the spectral bands.
We also applied the denoising methods on a real dataset.

Classification Application
In this section, we investigate the effect of the HSI denoising techniques as a preprocessing step for HSI classification.In order to evaluate the performance of different denoising approaches, we have applied three well-known classifiers including support vector machines (SVM) [70], random forest (RF) classifiers [71], and extreme learning machines (ELM) [72].
An SVM tries to separate training samples belonging to different classes by locating maximum margin hyperplanes in the multidimensional feature space where the samples are mapped [73].SVMs were originally introduced to solve linear classification problems.However, they can be generalized to nonlinear decision functions by considering the so-called kernel trick.A kernel-based SVM (using the Radial Basis Function (RBF) kernel) projects the pixel vectors into a higher dimensional space where the available samples are linearly separable and estimates maximum margin hyperplanes in this new space in order to improve the linear separability of data.
RF [71] is an ensemble method (a collection of tree-like classifiers) based on decision trees for classification.Ensemble classifiers run several (an ensemble of) classifiers which are individually trained, after which the individual results are combined through a voting process.Ideally, an RF classifier should be an independent and identically distributed randomization of weak learners.The RF classifies an input vector by running down each decision tree (a set of binary decisions) in the forest (the set of all trees).Each tree leads to a unit vote for a particular class and the forest chooses the eventual classification label based on the highest number of votes.
ELMs have been developed to train single layer feedforward neural networks (SLFN).Traditional gradient-based learning algorithms assume that all the parameters of the feedforward network, including weight and bias, need to be tuned.In [74], it was shown that the input weights w i and the hidden layer biases b i of the network can be initialized randomly in the beginning of the learning process and the hidden layer H can remain unchanged during the learning process.Therefore, by fixing the input weights w i and the hidden layer biases b i , one can train the SLFN in a similar manner to find a least-squares solution α of the linear system Hα = Y where α is the weight vector which connects the ith hidden node and the output nodes.In contrast with the traditional gradient-based approach, an ELM obtains not only the smallest training error but also the smallest norm of the output weights.Detailed information about the aforementioned pixel-wise classifiers can be found in [75,76] where the performances of the classifiers have been critically compared.
In order to compare classification accuracies obtained by different classifiers, three metrics have been applied: (1) Overall accuracy (OA), which is the number of correctly classified pixels, divided by the number of test samples; (2) Average accuracy (AA), which is the average value of the classification accuracies of all available classes; Kappa coefficient (K), which is a statistical measurement of agreement between the final classification map and the ground-truth map.

Setup
In the case of the SVM, the RBF kernel (as mentioned above) is employed.The optimal hyperplane parameters C (parameter that controls the amount of penalty during the SVM optimization) and γ (spread of the RBF kernel) were selected in the range of C = 10 −2 , 10 −1 , ..., 10 4 and γ = 10 −3 , 10 −2 , ..., 10 4 using five-fold cross validation.In the case of the RF, the number of trees is set to 300.The value of the prediction variable is set approximately to the square root of the number of input bands.In the case of the ELM, the regularization parameter was selected in the range of C = 1, 10 1 , ..., 10 5 using five-fold cross validation.

Results
In this section, the above-mentioned classifiers (i.e., ELM, SVM, and RF) have been applied to the four datasets (i.e., Houston, Trento, Indian Pines, and Washington DC).With reference to Tables 6-9, the following points can be observed: 1.In general, the prior use of denoising approaches improves the performance of the subsequent classification technique compared to the use of the input data without denoising step.The improvements reported in the cases of Trento (up to 19.4% in OA using ELM) and Indian Pines (up to 26.31% in OA using ELM) confirms the importance of the use of denoising techniques as a preprocessing step for HSI classification.2. The use of denoising approaches is clearly advantageous for raw data (i.e., before performing any preprocessing or low SNR bands removal).For instance, the classification accuracies have been considerably improved for the Indian Pines dataset while the amount of improvement is not that significant for the Houston and Washington DC data whose low SNR bands had been already eliminated.3.In general, ELM and SVM demonstrate a superior performance for the classification of denoised datasets compared to RF. Figure 10 shows several classification maps obtained by applying ELM on the Indian Pines dataset, without denoising and denoised by using 2D-Wavelet, 3D-Wavelet, FORPDN, LRMR, NAILRMA and HyRes.As can be seen, the classification map obtained by ELM on the raw data dramatically suffers from salt and pepper noise.This issue can be partially addressed using all the denoising approaches investigated in this paper.In particular, LRMR and NAILRMA considerably reduce the salt and pepper noise and produce homogeneous regions in the classification maps.

Summary, Conclusion and Future Challenges
In the past decade, HSI denoising has been considerably evolved.Conventional denoising methods based on 2D modeling and convex optimization were not efficient enough for HSI due to the ignorance of spectral information.Therefore, advanced denoising techniques have been developed to take into account the HSI characteristics.The high correlation between spectral bands in HSI has been found very useful for HSI denoising and therefore denoising techniques have been developed to exploit spectral information.3D model-based and 3D filtering approaches [26,35] have been suggested to use spectral information by treating the HSI as a 3D datacube.Spectral and spectral-spatial penalty-based approaches [38,39,42,44] have been developed to incorporate spectral information.Moreover, the advantages of low-rank modeling have been explored.Many techniques have been proposed based on low-rank modeling [7,12,36,46,49,51,53,54].Low-rank based techniques have also been utilized for mixed noise removal in HSIs [59][60][61][62].
In this paper, the state-of-the-art and the recent developments in the area of HSI denoising have been presented.In addition, this paper provided a background for HSI modeling and denoising, in which HSI denoising challenges and the different noise types have been discussed.Experimental results presented in this paper provide a comparative study over different generations of denoising techniques and confirmed the advantage of the low-rank techniques over the other ones.The gain in SNR is up to 20 dB for very low input SNR i.e., SNR in = 5 dB.In more detail, based on the experimental results, one can conclude that conventional band by band denoising techniques have provided the poorest performance compared to 3D filtering and spectral denoising approaches.Additionally, the effects of denoising as a preprocessing step for HSI classification have been investigated on four datasets.The experimental results have shown that the performances of the denoising techniques are consistent for the three classifiers and four datasets used.In general, from the classification experiments one can conclude that exploiting the denoising techniques has improved the classification accuracies.For Trento and Indian Pines datasets the improvements of the OA are very substantial compared to the spectral classification without the prior use of a denoising approach.The improvements in the OA obtained using ELM to classify the Indian Pines and Trento datasets reach to 26.31% and 19.4%, respectively.Despite the developments in HSI denoising, several important challenges remain to be addressed.Future challenges include; further investigation of the contribution of the HSI denoising approaches as a preprocessing step for further HSI analysis such as unmixing, change detection, resolution enhancement etc., exploring a model selection criterion which is not restricted to the Gaussian noise model, noise parameter estimation, investigating of the influences of the different noise types and indication of the dominant noise, and incorporating high performance computing techniques to obtain computationally efficient implementations.Investigating the performances of the HSI denoising approaches for denoising hyperspectral image sequences [77] and thermal hyperspectral images [78] (where emissivity becomes important and therefore other type of noise i.e., dark current noise can be encountered) are also of interest.vec(XR T n 2 ) = (R n 2 ⊗ I n 2 )vec(X) = D h x (i) .
where D h x (i) is the first order horizontal differences of X. Analogously, it can be shown that D v x (i) contains the first order vertical differences of X.

Figure 1 .
Figure 1.Left: Hyperspectral data cube.Right: The reflectance of the material within a pixel.

Figure 2 .
Figure 2. The number of journal and conference papers that appeared in IEEE Xplore on the vibrant topic of hyperspectral image denoising within different time periods.

Fig. 3 .Fig. 4 .Fig. 4 :
Fig. 3. Statistics on papers related to hyperspectral image and signal processing published in IEEE journals during 2009-2012 and 2013-2016.The size of each pie is proportional to the number of papers.The total number of papers for each time period is shown at the center of each pie chart.

Figure 3 .
Figure 3. Houston-from top to bottom-a color composite representation of the hyperspectral data using bands 70, 50, and 20, as R, G, and B, respectively; training samples; test samples; and the corresponding color bar.

Figure 4 .
Figure 4. Trento-from top to bottom-a color composite representation of the hyperspectral data using bands 40, 20, and 10, as R, G, and B, respectively; training samples; test samples; and the corresponding color bar.

Figure 5 .
Figure 5. Indian Pines-(a) a color composite representation of the hyperspectral data; (b) test samples; (c) training samples; and (d) the corresponding color bar.

Figure 6 .
Figure 6.Washington DC Mall-(a) a color composite representation of the hyperspectral data; (b) training samples; (c) test samples; and the corresponding color bar.

Figure 7 .
Figure 7.Comparison of the performances of the studied HSI restoration methods applied on the simulated dataset w.r.t.different levels of input Gaussian noise.(a) SNR (dB); (b) MSAD (logarithmic scale).

Figure 9 Figure 8 .Figure 9 .
Figure 8. Visual comparison of the performances of the studied HSI restoration methods applied on the simulated dataset.

Table 2 .
Houston-Number of training and test samples.

Table 3 .
Trento-Number of training and test Samples.

Table 4 .
Indian Pines-Number of training and test Samples.

Table 5 .
Washingto DC Mall-Number of training and test samples.
2 and D 1 represent 2D and 1D wavelet bases, respectively.Note that D 2 (D 2 = D 1 ⊗ D 1 , see Appendix B) and D 1 project the signal spatially and spectrally, respectively.If D 1 = I, we obtain a 2D wavelet model:

Table 7 .
Trento-Classification accuracies obtained by different denoising approaches before using ELM, SVM, and RF.The metrics AA an OA are reported in percentage, the Kappa coefficient is unitless.

Table 8 .
Indian Pines-Classification accuracies obtained by different denoising approaches before using ELM, SVM, and RF.The metrics AA an OA are reported in percentage, the Kappa coefficient is unitless.

Table 9 .
Washington DC Mall-Classification accuracies obtained by different denoising approaches before using ELM, SVM, and RF.The metrics AA an OA are reported in percentage, the Kappa coefficient is unitless.