Hyperspectral Image Mixed Noise Removal Using Subspace Representation and Deep CNN Image Prior

Abstract

The ever-increasing spectral resolution of hyperspectral images (HSIs) is often obtained at the cost of a decrease in the signal-to-noise ratio (SNR) of the measurements. The decreased SNR reduces the reliability of measured features or information extracted from HSIs, thus calling for effective denoising techniques. This work aims to estimate clean HSIs from observations corrupted by mixed noise (containing Gaussian noise, impulse noise, and dead-lines/stripes) by exploiting two main characteristics of hyperspectral data, namely low-rankness in the spectral domain and high correlation in the spatial domain. We take advantage of the spectral low-rankness of HSIs by representing spectral vectors in an orthogonal subspace, which is learned from observed images by a new method. Subspace representation coefficients of HSIs are learned by solving an optimization problem plugged with an image prior extracted from a neural denoising network. The proposed method is evaluated on simulated and real HSIs. An exhaustive array of experiments and comparisons with state-of-the-art denoisers were carried out.

1. Introduction

Hyperspectral cameras measure the radiation arriving at a sensor with high spectral resolution over a sufficiently broad spectral band such that the acquired spectrum can be used to uniquely characterize and identify any given material [1]. Hyperspectral imaging plays an important role in remote sensing and has been used in a wide array of applications, such as the identification of various minerals in mining and oil industries [2], monitoring the development and health of crops in agriculture [3], detecting the development of cracks in pavements in civil engineering [4], and so on.

A hyperspectral image (HSI) is a three-dimensional data cube, where the first two dimensions represent the spatial domain and the third dimension represents the spectral domain. In contrast to multispectral imaging, hyperspectral cameras capture electromagnetic information in hundreds of narrow spectral bands, instead of multiple wide spectral bands. However, due to the decrease of the width of spectral bands, hyperspectral cameras receive fewer photons per band and tend to acquire images with a lower signal-to-noise ratio (SNR). The decreased SNR reduces the reliability of measured features or information extracted from HSIs [1]. Therefore, hyperspectral image denoising is a fundamental inverse problem before further applications.

The degradations linked with various mechanisms also result in different types of noise, such as Gaussian noise, Poissonian noise, impulse noise, dead-lines/stripes, and cross-track illumination variation. In this paper, we focus on the discussion of additive and signal-independent noise (namely, Gaussian noise, impulse noise, and dead-lines/stripes) and attack hyperspectral mixed noise composed of these additive noises.

The hyperspectral image denoising problem is usually studied by exploiting the characteristics of HSIs and noise. Some models for clean HSIs (i.e., priors, regularizers, constraints) have been demonstrated to work well in hyperspectral image denoising problems [5,6]. For example, the spectral–spatial adaptive hyperspectral total variation (SSAHTV) denoising method [7] uses a vector total variation (TV) regularization, which can promote localized step gradients within the image bands and align the “discontinuities” across the bands. Because of the very high spectral–spatial correlation, an HSI admits a sparse representation on a given basis or frame [8,9]; e.g., the Fourier basis [10], wavelet basis [11], Discrete Cosine basis, and data-adaptive basis [8]. For example, a 3D wavelet-based hyperspectral denoising method was introduced in [11], which took advantage of a fundamental property of HSIs; i.e., 3D wavelet coefficients of HSIs are sparse or compressible. The sparsity of coefficients means that the majority of coefficients are zero or close to zero. The compressibility of coefficients means that its elements have fast-decaying tails. In the 3D wavelet-based denoising method, the ${\ell }_1$ norm, jointly with a quadratic data fidelity term, is used to promote sparsity on the wavelet coefficients.

The low rankness of HSIs in the spectral domain is also a widely used image prior for solving the denoising problem [12,13]. By minimizing the rank of the estimated image [14], one can remove a bulk of the Gaussian noise in observations due to the fact that Gaussian noise is full-rank. The non-convex rank constraint is usually relaxed by minimizing the nuclear norm of HSIs, as done in the low-rank matrix recovery (LRMR) method [15] and in the weighted sum of weighted tensor nuclear norm minimization (WSWTNNM) method [16]. The low-rank structure of HSIs is also exploited by representing the spectral vectors of the clean image in an orthogonal subspace in [8,17,18,19,20] for Gaussian noise removal and in [21,22] for mixed noise removal. Subspace representation is an explicit low-rank representation in the sense that the rank is constrained by the dimension of subspace.

The self-similarity of HSIs also has been exploited for HSI denoising problems mainly in two ways: (a) the non-local (generalized) mean: for each patch, seek for its similar patches in the image and produce a patch estimate based on the found patches [17]; (b) dictionary learning: express each patch as a sparse representation in a given dictionary, which may be learned from the data. Patch-based learning is arguably the state-of-the-art in HSI denoising [8,18].

Recently, deep learning techniques also have been developed for image denoising. Aiming at single-band images or RGB images, representative networks are feed-forward denoising convolutional neural networks (DnCNNs) [23], a fast and flexible denoising convolutional neural network (FFDNet) [24], and a convolutional blind denoising network (CBDNet) [25]. Some neural networks have been conceived for HSI; for example, a spatial–spectral gradient network (SSGN) [26], a CNN-based HSI denoising method HSI-DeNet [27], and a novel deep spatio-spectral Bayesian posterior (DSSBP) framework [28]. As the performance of deep learning based denoisers highly depends on the quality and quantity of training data, a challenge of deep-denoisers is the lack of real HSIs that can be used as training data or the simulation of pairs of clean–noisy images close to real HSIs. However, once a network is trained properly, deep-denoisers are much faster than the traditional machine learning (ML)-based denoisers. One potential solution is to incorporate deep-denoisers that have been well-trained using vast amounts of RGB images into the HSI denoising framework. This paper studies this possibility in the sense that we incorporate the well-known single-band deep denoiser, FFDNet, into a mixed noise removal framework derived using a traditional ML technique. These fall into the area of research called the plug-and-play technique [29,30] or regularization by denoising (RED) framework [31].

1.1. Related Work

The subspace representation of spectral vectors in HSIs has been successfully used to remove Gaussian noise by regularizing the representation coefficients of HSIs. We refer to representative work on global and nonlocal low-rank factorizations (GLF) [17], fast hyperspectral denoising (FastHyDe) [8], and non-local meets global (NGmeet) [18]. The idea of regularizing the subspace representation coefficients of HSIs underpins state-of-the-art Gaussian denoisers and has been extended to address mixed noise. The challenge of this extension lies in the estimation of the subspace.

A spectral subspace can be estimated from an observed HSI by simply carrying out the singular value decomposition (SVD) of observed image matrix when noise is i.i.d. Gaussian. When noise is a mixture of Gaussian noise, stripes, and impulse noise, the spectral subspace is usually estimated jointly with subspace coefficients of the HSI; for example, in double-factor-regularized low-rank tensor factorization (LRTF-DFR) [21] and subspace-based nonlocal low-rank and sparse factorization (SNLRSF) [22]. The joint estimation of the subspace and the corresponding coefficients of the HSI usually produce poor estimates of the subspace when HSI is affected by severe mixed noise. To sidestep the joint estimation, one strategy is to estimate the subspace and the corresponding coefficients of the HSI separately, which motivated us to develop the L1HyMixDe [32] method.

The proposed method in this paper is also based on two approaches: (a) learning subspace and coefficients of the HSI separately and (b) regularizing representation coefficients. However, the new work is distinct from L1HyMixDe in terms of the method of learning the subspace and the regularization imposed on the subspace representation coefficients. The main difference in subspace learning is that L1HyMixDe performs median filtering band by band, exploiting spatial correlation, whereas the new work estimates a coarse image using Hampel filtering pixel by pixel, exploiting spectral correlation. The Hampel filtering is more effective in outlier removal procedures. Furthermore, as anomalous target detection is an important task in hyperspectral imaging, it occurs that both anomalous targets and sparse noise are sparsely distributed in the spatial domain. If only the spatial information is considered—e.g., L1HyMixDe performs median filtering band by band, exploiting spatial correlation—then anomalous targets may be wrongly detected as sparse noise and would not be represented in the estimated subspace. This motivated us to develop a new method exploiting spectral information; i.e., spectral signatures of anomalous targets are smooth while spectral signatures of materials corrupted by sparse noise have unusual jumps in the value of the pixels. Furthermore, L1HyMixDe regularizes subspace representation coefficients with a non-local image prior, while this work adopts a more powerful deep CNN image prior. Compared with using non-local patch-based image priors, the new method using a deep-learning-based image prior is much faster as long as the deep network has been well trained. A computational efficient HSI mixed noise removal method is of importance in practice.

1.2. Contributions

The work aims to estimate a clean HSI from observations corrupted by mixed noise (containing Gaussian noise, impulse noise, and dead-lines/stripes) by exploiting two main characteristics of hyperspectral data, namely low-rankness in the spectral domain and high correlation in the spatial domain. Contributions of this work are summarized as follows:

• Instead of estimating the subspace basis and the corresponding coefficients of HSIs jointly and iteratively, we decouple the estimation of the subspace basis and the corresponding coefficients. A new subspace learning method, which works independently from coefficient estimation and is robust to mixed noise, is proposed.
• An image prior extracted from a state-of-the-art neural denoising network, FFDNet, is seamlessly embedded within our HSI mixed noise removal framework, which is a successful combination of the traditional machine learning technique and deep learning technique.

This paper is organized as follows. Section 2 formally introduces a mixed noise removal method—HySuDeep. Section 3 extends the proposed method to address mixed noise containing non-i.i.d. Gaussian noise. Section 4 and Section 5 show and analyze the experimental results of the proposed method and the comparison methods. Finally, we present a conclusion of this paper in Section 6.

2. Problem Formulation

Some notations and tensor operations used in this paper and their definitions are provided in

Table 1. Notations And Definitions.

Notation	Definition
$\mathcal{X}\in {\mathbb{R}}^{I_1\times I_2\times I_3}$	Three-dimensional tensor (calligraphic letter)
$\mathbf{X}$	Matrix (boldface capital letter)
$\mathbf{x}$	Vector (boldface lowercase letter)
x	Scalar (italic lowercase letter)
Order N	Number of dimensions in a tensor
Mode n	nth dimension of a tensor
Mode-3 vectors of $\mathcal{X}$	$I_3$-dimensional vectors obtained from $\mathcal{X}$ by varying the 3rd index while keeping 1st and 2nd indices fixed.
Mode-3 slices of $\mathcal{X}$	Matrices obtained from $\mathcal{X}$ by varing every index but the 3rd index.
$\mathcal{X}(:,:,i)$	ith mode-3 slice of $\mathcal{X}$, a matrix obtained by fixing the mode-3 index of $\mathcal{X}$ to be i.
${\mathbf{X}}_{\left(3\right)}\in {\mathbb{R}}^{I_3\times (I_1*I_2)}$	Mode-3 unfolding of $\mathcal{X}$. A tensor can be unfolded into a matrix by rearranging its mode-3 vectors, which are the column
	vectors of ${\mathbf{X}}_{\left(3\right)}$.
$\mathcal{X}{\times }_3\mathbf{E}$	Tensor matrix multiplication. The mode-3 product of a tensor $\mathcal{X}\in {\mathbb{R}}^{I_1\times I_2\times I_3}$ by a matrix $\mathbf{E}\in {\mathbb{R}}^{J_n\times I_3}$ is a tensor
	$\mathcal{Y}\in {\mathbb{R}}^{I_1\times I_2\times J_n}$, denoted as $\mathcal{Y}=\mathcal{X}{\times }_3\mathbf{E}$, which is corresponding to a matrix multiplication, ${\mathbf{Y}}_{\left(3\right)}=\mathbf{E}{\mathbf{X}}_{\left(3\right)}$.
${\parallel \mathcal{X}\parallel }_F$	The definition of Frobenius norm of a matrix is extended to a tensor as follows: ${\parallel \mathcal{X}\parallel }_F=\sqrt{{\sum }_{i_1,\cdots ,i_N}{\left|x_{i_1,\cdots ,i_N}\right|}^2}$

.

2.1. Observation Model

Let $\mathcal{X}\in {\mathbb{R}}^{r\times c\times n_b}$ denote an underlying clean HSI with $r\times c$ pixels and $n_b$ bands. Assuming that noise is additive, we can write an observation model as

$$\mathcal{Y}=\mathcal{X}+\mathcal{S}+\mathcal{G},$$

(1)

where $\mathcal{Y},\mathcal{S},\mathcal{G}\in {\mathbb{R}}^{r\times c\times n_b}$ are observed HSI data, sparse noise, and Gaussian noise, respectively. Elements in $\mathcal{G}$ are assumed to be i.i.d. Gaussian with zero-mean and variance ${\sigma }^2$. Non-i.i.d. noise is discussed in Section 3.

2.2. Subspace Representation of HSIs

The high correlation in the spectral domain of HSIs implies that spectral vectors (i.e., mode-3 vectors) of an HSI $\mathcal{X}$ can be well represented in a low-dimensional subspace [1,8,19]. We take a cropped Washington DC Mall image (shown in Figure 1a) as an example. To analyze the correlation of pixel values in the spectral domain, we obtained the autocorrelation matrix of its spectral vectors by computing $\mathbf{E}\left[{\mathbf{X}}_{\left(3\right)}{\mathbf{X}}_{\left(3\right)}^T\right]$. From Figure 1b, we can see that the first 100 bands are highly correlated. Due to the high spectral correlation, the autocorrelation matrix can be approximated by a low-rank matrix. The low-rankness may also be inferred in Figure 1c, where the first two singular values (the largest) dominate all other singular values, which decrease to near-zero rapidly. From Figure 1d,e, we can see that two spectral vectors of size 191 can be well approximated by vectors in a subspace with a dimension of 10.

Based on the subspace representation of spectral vectors, we write

$$\mathcal{X}=\mathcal{Z}{\times }_3\mathbf{E}+\Xi ,$$

(2)

with $\mathbf{E}\in {\mathbb{R}}^{n_b\times p}$ ($p\ll n_b$), $\mathcal{Z}\in {\mathbb{R}}^{r\times c\times p}$, and $\Xi \in {\mathbb{R}}^{r\times c\times n_b}$. Matrix $\mathbf{E}$ holds an orthogonal basis for the signal subspace, and the entries of $\mathcal{Z}$ are representation coefficients of $\mathcal{X}$ with respect to $\mathbf{E}$. Hereafter, mode-3 slices of $\mathcal{Z}$ are termed the eigenimages. Tensor $\Xi$ represents residuals of a low-dimensional representation of $\mathcal{X}$. When SVD is adopted, most of the image variation can be compactly represented in a low-dimensional subspace, and the residual tensor contains a small amount of image energy. Thus, we omit $\Xi$ when reconstructing the image; that is, $\mathcal{X}=\mathcal{Z}{\times }_3\mathbf{E}$. Figure 2 illustrates the subspace representation of a hyperspectral image. The observation model (1) can be written as

$$\mathcal{Y}=\mathcal{Z}{\times }_3\mathbf{E}+\mathcal{S}+\mathcal{G}.$$

(3)

The subspace representation of spectral vectors has been demonstrated to be a powerful tool for solving some hyperspectral inverse problems, such as hyperspectral image super-resolution [33], hyperspectral image segmentation [34], hyperspectral image classification [35], and hyperspectral unmixing [36]. In this paper, we propose a hyperspectral image mixed noise removal method using subspace representation and deep CNN Image prior (HySuDeep), the flowchart of which is illustrated in Figure 3.

However, for the hyperspectral mixed noise removal problem, the challenge of exploiting subspace representation lies in the estimation of the subspace from observations corrupted not only by Gaussian noise but also by sparse noise.

2.3. Subspace Learning against Mixed Noise

When noise in an observed HSI is approximately Gaussian distributed, the subspace of the underlying clean image can be estimated easily from the noisy image by performing SVD on the noisy image matrix; that is, $\{\mathbf{U},$$\Sigma ,\mathbf{V}\}=\mathrm{SVD}\left({\mathbf{Y}}_{\left(\mathbf{3}\right)}\right)$. The estimated subspace is spanned by k eigenvectors corresponding to the k largest eigenvalues; that is, $\widehat{\mathbf{E}}={\left[\mathbf{U}\right]}_{:,1:k}$, where the subscript $\left(:,1:k\right)$ means extracting the first k columns from a matrix.

However, when noise is a mixture of Gaussian noise and stripe noise, performing SVD on a noisy image matrix is not an optimal way to learn the spectral subspace. SVD aims to find the most important dimensions in the data by finding the directions of maximum variance, and the addition of i.i.d. Gaussian noise increases the variance of each dimension uniformly (thus, it does not change the order of singular values and consequently does not change the estimation of subspace), whereas the addition of stripes increases the variance of each dimension non-uniformly (thus, it does change the order of singular values). Experimental analysis is given in Section 4.3.

One of our contributions is to estimate the spectral subspace from observations corrupted by mixed noise and thus propose an efficient denoising method, which does not need to estimate spectral subspace iteratively.

2.3.1. Outlier Removal Using Hampel Filtering

Since the existence of sparse noise significantly disturbs the identification of signal subspace, we may estimate a coarse image by removing the sparse noise from observations. An outlier removal step applying Hampel filtering to spectral vectors of the observed image is introduced below.

Given a spectral vector of observed HSI, denoted as $\mathbf{y}={[y_1,y_2,\cdots ,y_{n_b}]}^T\in {\mathbb{R}}^{n_b}$, and a sliding window of length q, we define the point-to-point median and median absolute deviation (MAD) estimates for the samples in the ith band using

$$\begin{array}{cc}\hfill m_i=& \mathrm{median}(y_{i-q},y_{i-q+1},y_{i-q+2},\cdots ,y_i,\cdots ,\hfill \\ \hfill & y_{i+q-2},y_{i+q-1},y_{i+q}),\hfill \end{array}$$

(4)

and

$$\mathrm{MAD}=\mathrm{median}\left(\right|y_{i-q}-m_i|,\cdots ,|y_{i+q}-m_i\left|\right),$$

(5)

respectively. The standard deviation (denoted as ${\sigma }_{i,MAD}$) is estimated by scaling the local median absolute deviation (MAD) by a factor of $\kappa =1/\left(\sqrt{2}{\mathrm{erf}}^{-1}(1/2)\right)\approx 1.4826$:

$${\sigma }_{i,MAD}=\kappa \mathrm{MAD}.$$

(6)

The Hampel filtering is a variation of the three-sigma rule of statistics, which is robust against outliers. If a sample $y_i$ is such that

$$|y_i-m_i|>n_{\sigma }{\sigma }_{i,MAD}$$

(7)

for a given threshold $n_{\sigma }$, then the Hampel filtering declares $y_i$ an outlier and replaces it with $m_i$; that is,

$${\tilde{y}}_i=\left\{\begin{array}{cc}{m}_i\hfill & \mathrm{if}|y_i-m_i|>n_{\sigma }{\sigma }_{i,MAD}.\hfill \\ y_i\hfill & \mathrm{otherwise}.\hfill \end{array}\right.$$

(8)

Here, we simply set $n_{\sigma }$ to 3.

The coarse image $\tilde{\mathcal{Y}}$ is obtained by applying the Hampel filtering (8) to each spectral vector of the observed image. We remark that although the Hampel filtering can detect the positions of impulse noise and stripes well, it may not produce exact spectral estimates for the polluted pixels. Thus, the goal of performing the Hampel filtering on the observed image is to obtain a coarse image for subspace learning.

2.3.2. Subspace Identification

Given the coarse image $\tilde{\mathbf{Y}}$ (mainly containing Gaussian noise, not sparse noise), the subspace basis can be directly obtained by performing SVD on $\tilde{\mathcal{Y}}$, as $\tilde{\mathcal{Y}}$ mainly contains Gaussian noise which is independent and identically distributed (i.i.d.). That is, the signal subspace is approximately spanned by the first p left-singular vectors of $\tilde{\mathbf{Y}}$:

$$\widehat{\mathbf{E}}=\tilde{\mathbf{U}}(:,1:p),$$

(9)

where $\tilde{\mathbf{U}}\in {\mathbb{R}}^{n_b\times n_b}$ is an orthogonal matrix and $\{\tilde{\mathbf{U}},\tilde{\Sigma },\tilde{\mathbf{V}}\}=\mathrm{SVD}\left({\tilde{\mathbf{Y}}}_{\left(3\right)}\right)$ with singular values in $\tilde{\Sigma }$ ordered by non-increasing magnitude. The dimension of subspace, p, can be automatically estimated by the HySime method [37].

2.4. Fast Eigenimage Learning

2.4.1. Objective Function

Given the estimate of spectral subspace and the noise assumption (i.e., a mixture of i.i.d. Gaussian noise and sparse noise), we estimate the unknown variables, $\mathcal{Z}$ and $\mathcal{S}$, in (3) by solving

$$\{\widehat{\mathcal{Z}},\widehat{\mathcal{S}}\}\in arg\underset{\mathcal{Z},\mathcal{S}}{min}\frac{1}{2}\left|\right|\mathcal{Y}-\mathcal{Z}{\times }_3\widehat{\mathbf{E}}-{\mathcal{S}\left|\right|}_F^2+{\lambda }_1\varphi \left(\mathcal{Z}\right)+{\lambda }_2{\left|\right|\mathcal{S}\left|\right|}_1.$$

(10)

The first term on the above optimization problem is derived from the observation model (3). The third term is a regularizer with respect to matrix $\mathcal{S}$. The ${\ell }_1$ norm of $\mathcal{S}$, given by ${\parallel \mathcal{S}\parallel }_1={\sum }_{i=1}^r{\sum }_{j=1}^c{\sum }_{k=1}^{n_b}\left|\mathcal{S}(i,j,k)\right|$, promotes the sparsity of $\mathcal{S}$. Finally, ${\lambda }_1,{\lambda }_2\ge 0$ are regularization parameters that trade off the importance of the respective regularizers. If the $\varphi$ is a convex function, then the optimization (10) is a convex problem.

The second term in the optimization problem (10) is a regularization expressing prior information tailored to spatially correlated eigenimages [8,17]. We use the plug-and-play technique [29,30] to assign a prior for eigenimages; rather than investing efforts in designing new regularization for eigenimages, we can plug an image prior extracted from an off-the-shelf denoiser conceived for natural images. A deep image prior extracted from a CNN-based denoising network, FFDnet, is employed, and thus an explicit definition of the function $\varphi (·)$ is not given in optimization (10). More detailed analysis with regard to the deep prior $\varphi (·)$ can be found in Section 2.4.3.

Let $\mathcal{A}=$ concat $\left(\right[\mathcal{Z},\mathcal{S}],3)$ be a $r\times c\times (p+n_b)$ tensor that concatenates the eigenimages $\mathcal{Z}$ and the sparse noise $\mathcal{S}$ along the third dimension. The optimization problem (10) can be rewritten as

$$\begin{array}{cc}\widehat{\mathcal{A}}\hfill & \in arg\underset{\mathcal{A}}{min}\frac{1}{2}{\parallel \mathcal{Y}-\mathcal{A}{\times }_3[\widehat{\mathbf{E}},{\mathbf{I}}_{n_b}]\parallel }_F^2+\hfill \\ \hfill & {\lambda }_1\varphi \left(\mathcal{A}{\times }_3[{\mathbf{I}}_p,{\mathbf{0}}_{(p\times n_b)}]\right)+{\lambda }_2{\parallel \mathcal{A}{\times }_3[{\mathbf{0}}_{(n_b\times p)},{\mathbf{I}}_{n_b}]\parallel }_1,\hfill \end{array}$$

(11)

where ${\mathbf{I}}_a$ denotes an identity matrix of size a and ${\mathbf{0}}_{(a\times b)}$ is a zero matrix of size $a\times b$.

2.4.2. Solver

The optimization problem (11) can be solved by the constrained-split augmented Lagrangian shrinkage algorithm (C-SALSA) [38,39]. C-SALSA is an instance of an alternating direction method of multipliers (ADMM) [40] that was developed to solve the constrained optimization formulation of regularized image restoration. By using variable splitting, we can convert the original optimization into a constrained one:

$$\begin{array}{cc}\hfill \underset{\mathcal{A},{\mathcal{V}}_1,{\mathcal{V}}_2,{\mathcal{V}}_3}{min}& \frac{1}{2}\parallel \mathcal{Y}-{\mathcal{V}}_1{\parallel }_F^2+{\lambda }_1\varphi \left({\mathcal{V}}_2\right)+{\lambda }_2{\parallel {\mathcal{V}}_3\parallel }_1\hfill \\ \hfill \mathrm{s}.\mathrm{t}.& {\mathcal{V}}_1=\mathcal{A}{\times }_3[\hat{\mathbf{E}},{\mathbf{I}}_{n_b}]\hfill \\ \hfill & {\mathcal{V}}_2=\mathcal{A}{\times }_3[{\mathbf{I}}_p,{\mathbf{0}}_{(p\times n_b)}]\hfill \\ \hfill & {\mathcal{V}}_3=\mathcal{A}{\times }_3[{\mathbf{0}}_{(n_b\times p)},{\mathbf{I}}_{n_b}].\hfill \end{array}$$

The augmented Lagrangian function of the above optimization is

$$\begin{array}{cc}L(\mathcal{A},\hfill & {\mathcal{V}}_1,{\mathcal{V}}_2,{\mathcal{V}}_3,{\mathcal{D}}_1,{\mathcal{D}}_2,{\mathcal{D}}_3)=\hfill \\ \hfill & \frac{1}{2}\parallel \mathcal{Y}-{\mathcal{V}}_1{\parallel }_F^2+{\lambda }_1\varphi \left({\mathcal{V}}_2\right)+{\lambda }_2{\parallel {\mathcal{V}}_3\parallel }_1\hfill \\ \hfill & +\frac{{\mu }_1}{2}{\parallel {\mathcal{V}}_1-\mathcal{A}{\times }_3[\widehat{\mathbf{E}},{\mathbf{I}}_{n_b}]-{\mathcal{D}}_1\parallel }_F^2\hfill \\ \hfill & +\frac{{\mu }_2}{2}{\parallel {\mathcal{V}}_2-\mathcal{A}{\times }_3[{\mathbf{I}}_p,{\mathbf{0}}_{(p\times n_b)}]-{\mathcal{D}}_2\parallel }_F^2\hfill \\ \hfill & +\frac{{\mu }_3}{2}{\parallel {\mathcal{V}}_3-\mathcal{A}{\times }_3[{\mathbf{0}}_{(n_b\times p)},{\mathbf{I}}_{n_b}]-{\mathcal{D}}_3\parallel }_F^2,\hfill \\ \hfill & +\mathrm{constant}\hfill \end{array}$$

(12)

where ${\mu }_1,{\mu }_2,{\mu }_3>0$ are C-SALSA penalty parameters.

The application of C-SALSA to (13) leads to Algorithm 1.  A MATLAB demo of the proposed HySuDeep method will available at https://github.com/LinaZhuang (accessed on 8 October 2021).

Algorithm 1HySuDeep for mixed noise containing i.i.d. Gaussian noise

1: Set $k=0$, choose ${\mu }_1,{\mu }_2,{\mu }_3>0$, ${\mathcal{V}}_{1,0}$, ${\mathcal{V}}_{2,0}$, ${\mathcal{D}}_{1,0}$, ${\mathcal{D}}_{2,0}$. 2: repeat 3: $$\begin{gathered} \hspace{0.33em}{\mathcal{A}}_{k+1}=arg{min}_{\mathcal{A}}\frac{{\mu }_1}{2}\parallel {\mathcal{V}}_{1,k}-\mathcal{A}{\times }_3[\widehat{\mathbf{E}},{\mathbf{I}}_{n_b}]-{\mathcal{D}}_{1,k}{\parallel }_F^2+\frac{{\mu }_2}{2}\parallel {\mathcal{V}}_{2,k}-\mathcal{A}{\times }_3[{\mathbf{I}}_p,{\mathbf{0}}_{(p\times n_b)}]- \\ {\mathcal{D}}_{2,k}{\parallel }_F^2+\frac{{\mu }_3}{2}{\parallel {\mathcal{V}}_{3,k}-\mathcal{A}{\times }_3[{\mathbf{0}}_{(n_b\times p)},{\mathbf{I}}_{n_b}]-{\mathcal{D}}_{3,k}\parallel }_F^2 \end{gathered}$$   4: $\hspace{0.33em}{\mathcal{V}}_{1,k+1}=arg{min}_{{\mathcal{V}}_1}\frac{1}{2}\parallel \mathcal{Y}-{\mathcal{V}}_1{\parallel }_F^2+\frac{{\mu }_1}{2}{\parallel {\mathcal{V}}_1-{\mathcal{A}}_{k+1}{\times }_3[\widehat{\mathbf{E}},{\mathbf{I}}_{n_b}]-{\mathcal{D}}_{1,k}\parallel }_F^2$   5: $\hspace{0.33em}{\mathcal{V}}_{2,k+1}=arg{min}_{{\mathcal{V}}_2}{\lambda }_1\varphi \left({\mathcal{V}}_2\right)+\frac{{\mu }_2}{2}{\parallel {\mathcal{V}}_2-{\mathcal{A}}_{k+1}{\times }_3[{\mathbf{I}}_p,{\mathbf{0}}_{(p\times n_b)}]-{\mathcal{D}}_{2,k}\parallel }_F^2$   6: $\hspace{0.33em}{\mathcal{V}}_{3,k+1}=arg{min}_{{\mathcal{V}}_3}{\lambda }_2\parallel {\mathcal{V}}_3^T{\parallel }_1+\frac{{\mu }_3}{2}{\parallel {\mathcal{V}}_3-{\mathcal{A}}_{k+1}{\times }_3[{\mathbf{0}}_{(n_b\times p)},{\mathbf{I}}_{n_b}]-{\mathcal{D}}_{3,k}\parallel }_F^2$   7: $\hspace{0.33em}{\mathcal{D}}_{1,k+1}={\mathcal{D}}_{1,k}-({\mathcal{V}}_{1,k+1}-{\mathcal{A}}_{k+1}{\times }_3[\widehat{\mathbf{E}},{\mathbf{I}}_{n_b}])$   8: $\hspace{0.33em}{\mathcal{D}}_{2,k+1}={\mathcal{D}}_{2,k}-({\mathcal{V}}_{2,k+1}-{\mathcal{A}}_{k+1}{\times }_3[{\mathbf{I}}_p,{\mathbf{0}}_{(p\times n_b)}])$   9: $\hspace{0.33em}{\mathcal{D}}_{3,k+1}={\mathcal{D}}_{3,k}-({\mathcal{V}}_{3,k+1}-{\mathcal{A}}_{k+1}{\times }_3[{\mathbf{0}}_{(n_b\times p)},{\mathbf{I}}_{n_b}])$   10: $\hspace{0.33em}k\to k+1$ 11: until stopping criterion is satisfied.

The optimizations with regard to $\mathcal{A}$ and ${\mathcal{V}}_1$ (on lines 3 and 4 of Algorithm 1) are quadratic problems, whose solutions are given by

$$\begin{array}{cccc}{\mathcal{A}}_{k+1}& =& & (\mu 1({\mathcal{V}}_1-{\mathcal{D}}_1){\times }_3{[\widehat{\mathbf{E}},{\mathbf{I}}_{n_b}]}^T+\mu 2({\mathcal{V}}_2-{\mathcal{D}}_2){\times }_3\hfill \\ & & & {[{\mathbf{I}}_p,{\mathbf{0}}_{(p\times n_b)}]}^T+\mu 3({\mathcal{V}}_3-{\mathcal{D}}_3){\times }_3{[{\mathbf{0}}_{(n_b\times p)},{\mathbf{I}}_{n_b}]}^T){\times }_3\hfill \\ & & & (\mu 1{[\widehat{\mathbf{E}},{\mathbf{I}}_{n_b}]}^T[\widehat{\mathbf{E}},{\mathbf{I}}_{n_b}]+\mu 2{[{\mathbf{I}}_p,{\mathbf{0}}_{(p\times n_b)}]}^T\hfill \\ & & & {[{\mathbf{I}}_p,{\mathbf{0}}_{(p\times n_b)}]+\mu 3{[{\mathbf{0}}_{(n_b\times p)},{\mathbf{I}}_{n_b}]}^T[{\mathbf{0}}_{(n_b\times p)},{\mathbf{I}}_{n_b}])}^{-1}\hfill \end{array}$$

(13)

(for line 3) and

$${\mathcal{V}}_{1,k+1}={(1+{\mu }_1)}^{-1}\left[\mathcal{Y}+{\mu }_1({\mathcal{A}}_{k+1}{\times }_3[\widehat{\mathbf{E}},{\mathbf{I}}_{n_b}]+{\mathcal{D}}_{1,k})\right].$$

(14)

(for line 4).

Line 6 of Algorithm 1 is a proximity operator of the ${\ell }_1$ norm applied to ${\mathcal{V}}_{3,k}^{\prime }={\mathcal{A}}_{k+1}{\times }_3[{\mathbf{0}}_{(n_b\times p)},{\mathbf{I}}_{n_b}]+{\mathcal{D}}_{3,k}$. It can be solved by an element-wise soft-threshold function; that is,

$$\begin{array}{cc}\hfill {\mathcal{V}}_{3,k+1}=& arg\underset{{\mathcal{V}}_3}{min}{\lambda }_2\parallel {\mathcal{V}}_3{\parallel }_1+\frac{{\mu }_3}{2}{\parallel {\mathcal{V}}_3-{\mathcal{V}}_{3,k}^{\prime }\parallel }_F^2\hfill \\ \hfill =& [\mathrm{thresh}({\mathcal{V}}_{3,k}^{\prime }(i,j,k),\lambda /{\mu }_3),i=1,\cdots ,r,\hfill \\ \hfill & j=1,\cdots ,c,k=1,\cdots ,n_b],\hfill \end{array}$$

(15)

and threshold $(x,\tau )$ returns the soft thresholding of x (where $\tau$ is the threshold value) [41]:

$$x\mapsto \frac{max\left(\right|x|-\tau ,0)}{max\left(\right|x|-\tau ,0)+\tau }x.$$

Line 5 of Algorithm 1 is a proximity operator of $\varphi$ applied to ${\mathcal{V}}_{2,k}^{\prime }={\mathcal{A}}_{k+1}{\times }_3[{\mathbf{I}}_p,{\mathbf{0}}_{(p\times n_b)}]+{\mathcal{D}}_{2,k}$. As the scaled Lagrange multiplier ${\mathcal{D}}_{2,k}$ tends to 0 as the C-SALSA iteration goes on, ${\mathcal{V}}_{2,k}^{\prime }$ tends to be close to ${\mathcal{Z}}_{k+1}$. Considering that orthogonal projection is a decorrelation transformation and mode-3 slices of $\mathcal{Z}$ tend to be decorrelated, it is reasonable to decouple $\varphi (·)$ with respect to the mode-3 slices; that is,

$$\varphi \left({\mathcal{V}}_{2,k}^{\prime }\right)=\sum _{i=1}^p{\varphi }_i\left({\mathcal{V}}_{2,k}^{\prime }(:,:,i)\right),$$

(16)

where ${\mathcal{V}}_{2,k}^{\prime }(:,:,i)$ denotes the ith mode-3 slice. The solution of the subproblem on line 5 can be decoupled with regard to mode-3 slices of ${\mathcal{V}}_{2,k}^{\prime }$ and can be written as

$${\mathcal{V}}_{2,k+1}={\mathsf{\Psi }}_{{\lambda }_1\varphi /{\mu }_2}\left({\mathcal{V}}_{2,k}^{\prime }\right)=\left[\begin{array}{c}{\mathsf{\Psi }}_{{\lambda }_1{\varphi }_1/{\mu }_2}\left({\mathcal{V}}_{2,k}^{\prime }(:,:,1)\right)\\ ⋮\\ {\mathsf{\Psi }}_{{\lambda }_1{\varphi }_p/{\mu }_2}\left({\mathcal{V}}_{2,k}^{\prime }(:,:,p)\right)\end{array}\right]$$

(17)

where

$${\mathsf{\Psi }}_{\lambda {\varphi }_i}\left(\mathbf{U}\right)=arg\underset{\mathbf{X}}{min}\frac{1}{2}{\parallel \mathbf{X}-\mathbf{U}\parallel }_F^2+\lambda {\varphi }_i\left(\mathbf{X}\right)$$

(18)

is the so-called Moreau proximity operator (MPO) [42] of $\varphi$ or the denoising operator.

2.4.3. Plug-and-Play Prior, $\varphi (·)$

Instead of investing great effort in conceiving a new eigenimage regularization ${\varphi }_i\left(\mathbf{X}\right)$, we make use of the plug-and-play (PnP) prior [8,29,30]. The PnP tool is a powerful framework that solves an inverse image problem by directly using an existing regularizer from a state-of-the-art image denoiser. The PnP framework has been applied to hyperspectral image fusion, spectral unmixing, hyperspectral image inpainting, and anomaly detection, and produces surprisingly good image recovery results. Inspired by these successful applications of PnP, we plug the prior of a deep denoising network, FFDnet [24], into (18), leading to

$${\widehat{\mathcal{V}}}_{2,k}(:,:,i)\leftarrow \mathrm{FFDNet}\left({\mathcal{V}}_{2,k}^{\prime }(:,:,i)\right),$$

(19)

where FFDNet $(·)$ represents the network and outputs a denoised image. It is worth noting that other state-of-the-art image denoisers can also be plugged in (18). We choose FFDNet in this work based on the following considerations: (a) deep-learning-based FFDNet is much faster than other machine learning-based denoisers as long as it has been well trained and (b) FFDNet is flexible and able to address images with various noise levels without retraining the network.

Since the denoiser FFDNet, plugged into the subproblem with regard to ${\mathcal{V}}_2$, is not a proximity operator, we do not have a theoretical convergence guaranteed for the implemented variant of C-SALSA. However, impressive performances using plugged non-convex regularizations have been observed in a number of inverse problems [19,43,44]. The numerical convergence of the proposed HySuDeep is analyzed in Section 4.5.

2.5. HSI Recovery

After estimating $\widehat{\mathcal{Z}}$, the denoised image is obtained as follows:

$$\widehat{\mathcal{X}}=\widehat{\mathcal{Z}}{\times }_3\widehat{\mathbf{E}}.$$

(20)

3. Model Extension to Non-i.i.d. Gaussian Noise

Gaussian noise in the observation model (1) is assumed to be i.i.d., meaning that the model is simplified and we can focus on the discussion of the image structure. Let ${\mathbf{C}}_{\lambda }=\mathbb{E}\left[{\mathbf{n}}_i{\mathbf{n}}_i^T\right]$ (where ${\mathbf{n}}_i$ is a generic column of unfolded mode-3 matrix ${\mathbf{N}}_{\left(3\right)}$) define the covariance matrix of the spectral noise. We have ${\mathbf{C}}_{\lambda }={\sigma }^2{\mathbf{I}}_{n_b}$ when the noise is i.i.d.

However, we remark that Gaussian noise in real HSIs tends to be non-i.i.d.; that is, it is pixel-wise independent but band-wise dependent. Before implementing denoising according to model (1), we need to convert the non-i.i.d. scenario into an i.i.d. scenario by implementing data whitening:

$$\overline{\mathcal{Y}}=\mathcal{Y}{\times }_3\sqrt{{\mathbf{C}}_{\lambda }^{-1}},$$

(21)

where $\sqrt{{\mathbf{C}}_{\lambda }}$ and ${\mathbf{C}}_{\lambda }^{-1}$ denote the square root of ${\mathbf{C}}_{\lambda }$ and the inverse of ${\mathbf{C}}_{\lambda }$, respectively.

The estimate of the covariance matrix, ${\mathbf{C}}_{\lambda }$, is challenging when noise in observations is a mixture of Gaussian and sparse noise. If the HSI is only corrupted by Gaussian noise, then the noise covariance matrix can be estimated by exploiting the spectral correlation of HSI. The spectrally uncorrelated components extracted from the HSI are considered as Gaussian noise, whose covariance matrix can be computed easily. This problem has been studied deeply, and we refer readers to representative works [37,45]. However, when the noise in the image is a mixture of Gaussian noise and sparse noise, the spectrally uncorrelated components extracted from the HSI are mixed noise, instead of single Gaussian noise. To solve this problem, we propose a method to split two kinds of noise and estimate ${\mathbf{C}}_{\lambda }$ in two steps, as listed below:

• Estimate a coarse image, $\tilde{\mathcal{Y}}$, by removing the sparse noise from observation. A outlier removal step applying Hampel filtering to spectral vectors of the observed image is given in $\left(8\right)$.
• We apply linear regression to each band of the image $\tilde{\mathcal{Y}}$; i.e., each band is represented as a linear combination of the remaining bands [37]. That is,

  $${\left[{\tilde{\mathbf{Y}}}_{\left(3\right)}^T\right]}_{:,i}={\left[{\tilde{\mathbf{Y}}}_{\left(3\right)}^T\right]}_{{\alpha }_i}{\mathbf{\beta }}_i+{\mathbf{\xi }}_i,$$

  (22)

  where ${\tilde{\mathcal{Y}}}_{\left(3\right)}^T$ denotes the transpose of the mode-3 unfolding matrix ${\tilde{\mathbf{Y}}}_{\left(3\right)}$, the subscript ${[·]}_{:,i}$ means extracting the ith column from a matrix, a matrix with the subscript ${[·]}_{{\alpha }_i}$ means the matrix including all columns except ith column, ${\mathbf{\beta }}_i\in {\mathbb{R}}^{n_b-1}$ denotes regression coefficients, and ${\widehat{\mathbf{\xi }}}_i\in {\mathbb{R}}^n$ denotes regression error.
  The regression coefficients ${\mathbf{\beta }}_i$ can be estimated by the least squares method; i.e.,

  $$\widehat{\mathbf{\beta }}=arg\underset{\mathbf{\beta }}{min}{\parallel {\left[{\tilde{\mathbf{Y}}}_{\left(3\right)}^T\right]}_{:,i}-{\left[{\tilde{\mathbf{Y}}}_{\left(3\right)}^T\right]}_{{\alpha }_i}{\mathbf{\beta }}_i\parallel }_F^2.$$

  (23)
  Given $\widehat{\mathbf{\beta }}$, the regression error, ${\widehat{\mathbf{\xi }}}_i$, is computed by

  $${\widehat{\mathbf{\xi }}}_i={\left[{\tilde{\mathbf{Y}}}_{\left(3\right)}^T\right]}_{:,i}-{\left[{\tilde{\mathbf{Y}}}_{\left(3\right)}^T\right]}_{{\alpha }_i}{\widehat{\mathbf{\beta }}}_i.$$

  (24)
  The regression errors are taken as a coarse estimate of the Gaussian noise; thus, its covariance matrix, ${\mathbf{C}}_{\lambda }=\mathrm{diag}\left([{\sigma }_1^2,{\sigma }_2^2,\cdots ,{\sigma }_{n_b}^2]\right)$, can be estimated by

  $${\sigma }_i^2=\mathrm{var}\left({\widehat{\mathbf{\xi }}}_i\right),$$

  (25)

  where $\mathrm{var}(·)$ is a function computing the variance of the elements of the input vector.

Given the estimate of ${\mathbf{C}}_{\lambda }$, we can convert band-dependent Gaussian noise to i.i.d. Gaussian noise via (21). Here, we analyze the impact of image conversion on Gaussian noise. Given the image conversion,

$$\overline{\mathcal{Y}}=\mathcal{X}{\times }_3\sqrt{{\mathbf{C}}_{\lambda }^{-1}}+\mathcal{S}{\times }_3\sqrt{{\mathbf{C}}_{\lambda }^{-1}}+\mathcal{G}{\times }_3\sqrt{{\mathbf{C}}_{\lambda }^{-1}},$$

(26)

we can compute the spectral covariance of the Gaussian noise in the converted image as follows:

$$\begin{array}{cc}\hfill {\overline{\mathbf{C}}}_{\lambda }& =\mathbb{E}\left[\sqrt{{\mathbf{C}}_{\lambda }^{-1}}{\mathbf{n}}_i{\left(\sqrt{{\mathbf{C}}_{\lambda }^{-1}}{\mathbf{n}}_i\right)}^T\right]=\sqrt{{\mathbf{C}}_{\lambda }^{-1}}\mathbb{E}\left[{\mathbf{n}}_i{\mathbf{n}}_i^T\right]{\sqrt{{\mathbf{C}}_{\lambda }^{-1}}}^T\hfill \\ \hfill & =\sqrt{{\mathbf{C}}_{\lambda }^{-1}}{\mathbf{C}}_{\lambda }{\sqrt{{\mathbf{C}}_{\lambda }^{-1}}}^T=\mathbf{I}.\hfill \end{array}$$

(27)

From (27), we can see that after image conversion, the Gaussian noise is i.i.d. and standard distributed, which meets the noise assumption in model (1); thus, the converted image, $\overline{\mathcal{Y}}$, can be denoised as discussed in Section 2. Finally, we reconstruct the clean image as follows:

$$\widehat{\mathcal{X}}=\widehat{\overline{\mathcal{X}}}{\times }_3\sqrt{{\mathbf{C}}_{\lambda }},$$

(28)

where $\widehat{\overline{\mathcal{X}}}$ is the estimated clean version of image $\overline{\mathcal{Y}}$.

The pseudocode in Algorithm 2 shows how HySuDeep is implemented to reduce mixed noise for an HSI. Given an HIS of size r (rows) $\times c$ (columns) $\times n_b$ (bands) with subspace dimension p ($p\ll n_b$), the computational complexity of obtaining a coarse image $\tilde{\mathcal{Y}}$ in line 2 and ${\mathbf{C}}_{\lambda }$ in line 3 is $\mathcal{O}(r*c*n_b*q*logq)$ and $\mathcal{O}(r*c*n_b^3)$, respectively. The Gaussian noise-whitening in line 4, its inverse transformation in line 8, and the image reconstruction step in line 7 have the same computational complexity; that is, $\mathcal{O}(r\times c\times n_b^2)$. The estimation of the spectral subspace in line 5 and eigenimages in line 6, respectively, cost $\mathcal{O}(r^2*c^2*n_b)$ and $\mathcal{O}(r*c*{(p+n_b)}^2+{(p+n_b)}^3+p*d)$. Here, d represents the computational complexity of denoising an eigenimage using FFDNet. Consequently, the overall computational complexity of HySuDeep is $\mathcal{O}(r*c*n_b^3+r^2*c^2*n_b+p*d)$.

Algorithm 2HySuDeep for mixed noise containing non-i.i.d. Gaussian noise

1: Input a noisy HSI, $\mathcal{Y}$, and parameter ${\lambda }_2$. 2: Estimate a coarse image $\tilde{\mathcal{Y}}$ via (8). 3: Estimate ${\mathbf{C}}_{\lambda }$ from $\tilde{\mathcal{Y}}$ via (25). 4: Obtain an whitened image $\overline{\mathcal{Y}}$ using (21). 5: Denoise the whitened image $\overline{\mathcal{Y}}$ using Algorithm 1, and obtain $\widehat{\overline{\mathcal{X}}}$. 6: Compute $\widehat{\mathcal{X}}=\widehat{\overline{\mathcal{X}}}{\times }_3\sqrt{{\mathbf{C}}_{\lambda }}$, an estimate of the clean HSI.

4. Experiments with Simulated Images

Experiments were carried out on three simulated datasets to evaluate the performance of the proposed HySuDeep method compared with six state-of-the-art denoising methods for HSI mixed noise removal.

4.1. Simulation of Noisy Images and Comparisons

4.1.1. Simulation of Noisy Images

Three public hyperspectral datasets (shown in Figure 4a–c) were employed to simulate the noisy images. Following the same procedure in [32], we generated noisy HSIs of Pavia University data (of size 310 (rows) × 250 (columns) × 87 (bands)), a subregion of Washington DC Mall data (of size 150 (rows) × 200 (columns) × 191 (bands)), and Terrain data (of size 500 (rows) × 307 (columns) × 166 (bands)).

For the first two datasets, we removed bands that were severely corrupted by water vapor in the atmosphere. To obtain relatively clean images, we projected spectral vectors of each image onto a subspace spanned by eight principal eigenvectors of each image. The projection of each image was considered to be a clean image. To simulate noisy HSIs, we added four kinds of additive noise into images as follows:

Case 1 (Gaussian non-i.i.d. noise): Synthetic data with Gaussian noise. The noise in ith pixel vector ${\mathbf{g}}_i\sim \mathcal{N}(\mathbf{0},{\mathbf{D}}^2)$, where $\mathbf{D}$ is a diagonal matrix with diagonal elements sampled from a uniform distribution $U(0,0.01)$.

Case 2 (Gaussian noise + stripes): Synthetic data with Gaussian noise (described in case 1) and oblique stripe noise randomly affecting 30% of the bands and, for each band, about 10% of the pixels, at random.

Case 3 (Gaussian noise + “Salt and Pepper” noise): Synthetic data with Gaussian noise (described in case 1) and “Salt and Pepper” noise with a noise density of 0.5%, meaning approximately 0.5% of elements in $\mathcal{X}$ are affected.

Case 4 (Gaussian noise + stripes + “Salt and Pepper” noise): Synthetic data with Gaussian noise (described in case 1), random oblique stripes (described in case 2) and “Salt and Pepper” noise (described in case 3).

To evaluate the impact of denoising methods on the hyperspectral unmixing task, we simulated a clean semi-real HSI based on the publically available Terrain image (Figure 4c) following generation steps in [36]. A MATLAB demo used to generate a semi-real HSI is available at https://github.com/LinaZhuang/NMF-QMV_demo, accessed on 8 October 2021. The original Terrain image has a size of 500 (rows) × 307 (columns) × 166 (bands) and is mainly composed of soil, tree, grass, and shadow. The number of endmembers is empirically set to 5, as performed in [36,46,47]. Briefly, a clean Terrain image was synthesized based on a linear mixing model [1]; i.e., $\mathcal{X}=\mathcal{A}{\times }_3\mathbf{M}$, where $\mathbf{M}$ and $\mathcal{A}$ are endmembers and abundances, respectively, estimated from the original Terrain image. Next, we generated a noisy Terrain HSI by adding the Gaussian noise and oblique stripe noise (as described in case 2)—i.e., $\mathcal{Y}=\mathcal{A}{\times }_3\mathbf{M}+\mathcal{S}+\mathcal{G}$—yielding an MPSNR of 28.90 dB.

4.1.2. Comparisons

To thoroughly evaluate the performance of the proposed method, six state-of-the-art HSI denoising methods were selected for comparison: the fast hyperspectral denoising method (FastHyDe) [8], the noise-adjusted iterative low-rank matrix approximation method (NAIRLMA) [48], the spatio-spectral total variation based method (SSTV) [49], the low-rank matrix recovery method (LRMR) [15], double-factor-regularized low-rank tensor factorization (LRTF-DFR) [21], and the ${\ell }_1$-norm based hyperspectral mixed noise denoising method (L1HyMixDe) [32]. These compared methods were carefully selected. The FastHyDe served as a benchmark to see whether mixed noise could be simply addressed by Gaussian denoisers. The NAIRLMA, SSTV, and LRMR are from highly cited papers working on hyperspectral mixed noise based on low-rank and sparse representations. The LRTF-DFR and L1HyMixDe methods are subspace representation methods and remove noise by filtering the subspace coefficients of HSIs. All experiments were implemented in MATLAB (R2016a) on Windows 10 with an Intel Core i7-7700HQ 2.8-GHz processor and 24 GB RAM. A MATLAB demo of this work will be available at https://github.com/LinaZhuang, (accessed on 8 October 2021) for the sake of reproducibility.

The size of the sliding window in the Hampel filter, q, was fixed to 7. The dimension of the subspace, p, for FastHyDe, LRTF-DFR, L1HyMixDe, and HySuDeep methods was estimated by HySime [37] automatically. The other parameters of FastHyDe, NAILRMA, LRTF-DFR, and L1HyMixDe were set as suggested in their original references. We fine-tuned the parameters $\lambda$ and $\mathsf{\mu }$ in SSTV and parameters r and s in LRMR to obtain the optimal results.

For quantitative assessment, the peak signal-to-noise ratio (PSNR) index, the structural similarity (SSIM) index, and the feature similarity (FSIM) index of each band were calculated. The mean PSNR (MPSNR), mean SSIM (MSSIM), and mean FSIM (MFSIM) of the proposed and compared methods on Pavia University Data and Washington DC Mall Data are presented in

Table 2. Performance of the proposed and comparison methods on Pavia University Data and Washington DC Mall Data.

Cases	Indexes	Noisy	FastHyDe	NAILRMA	SSTV	LRMR	LRTF-DFR	L1HyMixDe	HySuDeep
			[8]	[48]	[49]	[15]	[21]	[32]
Pavia University data
Case 1	MPSNR	34.29	50.50	47.25	45.73	43.72	43.77	47.62	49.08
	MSSIM	0.8221	0.9981	0.9949	0.9929	0.9867	0.9906	0.9971	0.9978
	MFSIM	0.9272	0.9991	0.9977	0.9966	0.9945	0.9949	0.9986	0.999
	Time (s)	-	10	185	282	107	141	114	31
Case 2	MPSNR	27.77	40.94	32.33	35.06	33.58	43.43	47.26	48.55
	MSSIM	0.7451	0.9087	0.9236	0.9584	0.9176	0.9902	0.9969	0.9975
	MFSIM	0.8683	0.9569	0.9519	0.9702	0.9461	0.9948	0.9985	0.9989
	Time (s)	-	9	182	268	67	141	156	33
Case 3	MPSNR	22.56	39.66	37.03	43.85	43.37	43.72	47.90	48.93
	MSSIM	0.6892	0.9781	0.9556	0.9876	0.9858	0.9905	0.9971	0.9977
	MFSIM	0.8780	0.9870	0.9799	0.9948	0.9942	0.9948	0.9986	0.999
	Time (s)	-	9	123	250	87	141	112	32
Case 4	MPSNR	19.29	31.83	28.20	33.59	33.64	43.32	47.43	48.31
	MSSIM	0.6258	0.8249	0.8531	0.9463	0.9153	0.9901	0.9968	0.9972
	MFSIM	0.8267	0.9103	0.9169	0.9656	0.9455	0.9947	0.9985	0.9987
	Time (s)	-	10	193	274	89	140	114	38
Washington DC Mall data
Case 1	MPSNR	19.13	42.04	37.22	26.78	27.15	36.96	39.52	41.17
	MSSIM	0.7781	0.9985	0.9962	0.9674	0.9710	0.9961	0.9981	0.9983
	MFSIM	0.8634	0.9980	0.9931	0.9630	0.9706	0.9959	0.9978	0.999
	Time (s)	-	6	218	228	55	145	40	29
Case 2	MPSNR	12.52	30.58	20.56	24.08	25.39	34.12	36.65	40.05
	MSSIM	0.6851	0.8608	0.8440	0.9512	0.9636	0.9909	0.9946	0.9978
	MFSIM	0.7973	0.9368	0.8940	0.9547	0.9663	0.9948	0.9956	0.9987
	Time (s)	-	5	196	234	37	145	64	29
Case 3	MPSNR	7.29	29.78	24.94	26.70	27.29	36.62	39.90	40.75
	MSSIM	0.6421	0.9818	0.9642	0.9667	0.9707	0.9959	0.9982	0.9981
	MFSIM	0.8052	0.9849	0.9713	0.9625	0.9706	0.9958	0.9978	0.9988
	Time (s)	-	5	137	229	40	144	42	30
Case 4	MPSNR	3.84	20.34	15.82	24.03	24.94	34.03	36.81	39.07
	MSSIM	0.5676	0.8029	0.7958	0.9506	0.9617	0.9910	0.9947	0.997
	MFSIM	0.7508	0.8886	0.8663	0.9544	0.9649	0.9947	0.9956	0.998
	Time (s)	-	5	178	207	43	144	52	30

, where we have highlighted the best results in bold.

4.2. Mixed Noise Removal

We compare the proposed method HySuDeep with FastHyDe [8], NAILRMA [48], SSTV [49], LRMR [15], LRTF-DFR [21], and L1HyMixDe [32] methods in terms of MPSNR, MSSIM, and MFSIM (provided in Table 2 and

Table 3. Denoising performance of the proposed and comparison methods on the Terrain image and unmixing performance (in terms of NMSEA${}_{\mathbf{M}}$ and NMSES${}_{\mathcal{A}}$) of “VCA + FCLS” applied to the denoised images.

	Noisy	FastHyDe	NAILRMA	SSTV	LRMR	LRTF-DFR	L1HyMixDe	HySuDeep
MPSNR	28.90	48.94	40.72	35.58	39.84	45.84	50.84	51.91
MSSIM	0.7296	0.9933	0.9843	0.9168	0.9831	0.9932	0.9975	0.9987
MFSIM	0.8719	0.9961	0.9882	0.9548	0.9830	0.9953	0.9984	0.9989
Time (s)	-	20	587	951	129	572	562	137
NMSEA${}_{\mathbf{M}}$	0.10	0.19	0.17	0.20	0.05	0.06	0.03	0.04
NMSES${}_{\mathcal{A}}$	0.52	0.48	0.47	0.48	0.46	0.51	0.15	0.09

Note: NMSEA${}_{\mathbf{M}}$ and NMSES${}_{\mathcal{A}}$ denote root mean square error of endmembers and abundances, respectively.

). The band-wise PSNR is depicted in Figure 5 for quantitative assessment.

Among the compared methods, only FastHyDe is shown to address pure Gaussian noise, not mixed noise. We include FastHyDe to evaluate whether mixed noise can be reduced by a simple Gaussian denoiser. As shown in Table 2, FastHyDe outperformed other methods in case 1 (where HSIs are corrupted by only Gaussian noise). However, when the images were contaminated by mixed noise in Cases 2–4, the results of FastHyDe show that heavy mixed noise cannot be reduced well by a Gaussian denoiser. The existence of mixed noise in HSIs calls for effective denoising techniques.

To address mixed noise, an additive term is usually introduced to model sparse noise so that the entries in the data fidelity term, $(\mathcal{Y}-\mathcal{X}-\mathcal{S})$, follow Gaussian distributions. The proposed HySuDeep, NAIRLMA, SSTV, LRMR, and LRTF-DFR methods fall in this line of research. The critical differences among these methods are the regularizations imposed on the HSI, $\mathcal{X}$. SSTV minimizes the total variation of the HSI in the spectral–spatial domain. NAILRMR and LRMR impose spectral low-rankness in spatial patches of the image. LRTF-DFR and the proposed HySuDeep both enforce the low-rankness of spectral vectors by subspace representation. To exploit the spatial correlation of eigenimages, LRTF-DFR minimizes the spatial difference (i.e., total variation) of eigenimages, while HySuDeep uses a CNN-regularization (i.e., a deep image prior extracted from a CNN network).

The proposed method uniformly yields the best results in Cases 2–4, as shown in Table 2 and Table 3. The main reason is that, compared with SSTV, NAILRMR, LRMR, LRTF-DFR, and L1HyMixDe methods, our method uses a more efficient spatial regularization term. Although the CNN regularization works like a black box, we can see its superiority over TV regularization (used in SSTV and LRTF-DFR) and patch-based regularization (used in NAILRMA, LRMR, and L1HyMixDe) from Table 2 and Table 3. Furthermore, HySuDeep is similar to the L1HyMixDe method in terms of using subspace representations of spectral vectors and imposing a regularization on representation coefficients. Table 2 and Table 3 show that HySuDeep achieves better results than L1HyMixDe. The difference is caused by the accuracy of spectral subspace learning. Compared with L1HyMixDe, HySuDeep elaborates an outlier removal operator, which improves the estimates of subspace $\widehat{\mathbf{E}}$ (see discussion in Section 4.3).

For visual comparison, we display the 37th band of Pavia University data and the 126th band of Washington DC Mall data in Figure 6 and Figure 7, respectively. For case 1 (Gaussian noise), all methods can reduce noise significantly. As shown in Figure 6 and Figure 7, SSTV and LRMR are able to remove light stripes, but still leave some wide stripes. Heavy noise still remains in the results of NAILRMA and FastHyDe. LRTF-DFR, L1HyMixDe, and HySuDeep methods visually yield comparable results in Figure 6 and Figure 7. To show their differences, we generated and presented the residual images in Figure 8a,b, where we can see the results of proposed HySuDeep contain less residual compared with LRTF-DFR and L1HyMixDe.

4.3. Subspace Learning against Mixed Noise

One of the contributions of this paper is the introduction of an orthogonal subspace identification method that is robust to mixed noise. We conducted experiments using simulated images to compare the proposed subspace learning method with representative subspace learning methods: namely, SVD, RPCA, L1HyMixDe [32], and HySime [37]. Results are reported in

Table 4. Quantitative comparison of representation ability of subspaces learned by different methods.

		SVD	RPCA	L1HyMixDe	HySime	HySuDeep
		[51]	[52]	[32]	[37]
Pavia University data
Case 1	${\gamma }_1$	0.9994	0.9993	0.9995	1.0000	0.9996
	${\gamma }_2$	0.1789	0.1816	0.1753	0.0954	0.1500
Case 2	${\gamma }_1$	0.9854	0.9868	0.9994	0.9881	0.9996
	${\gamma }_2$	0.2552	0.2526	0.1599	0.2270	0.1314
Case 3	${\gamma }_1$	0.9993	0.9991	0.9995	0.9997	0.9996
	${\gamma }_2$	0.1095	0.1102	0.1066	0.0909	0.1019
Case 4	${\gamma }_1$	0.9855	0.9867	0.9993	0.9869	0.9997
	${\gamma }_2$	0.2355	0.2338	0.1526	0.2109	0.1271
Washington DC Mall data
Case 1	${\gamma }_1$	0.9998	0.9998	0.9998	1.0000	0.9998
	${\gamma }_2$	0.0772	0.0775	0.0751	0.0466	0.0630
Case 2	${\gamma }_1$	0.9826	0.9838	0.9997	0.9878	0.9997
	${\gamma }_2$	0.2738	0.2700	0.0978	0.2426	0.0585
Case 3	${\gamma }_1$	0.9997	0.9997	0.9998	0.9998	0.9998
	${\gamma }_2$	0.0483	0.0480	0.0456	0.0406	0.0441
Case 4	${\gamma }_1$	0.9825	0.9834	0.9997	0.9868	0.9997
	${\gamma }_2$	0.2342	0.2321	0.0883	0.2110	0.0600

Note: ${\gamma }_1=\frac{\parallel \mathcal{X}{\times }_3{\widehat{\mathbf{E}}}^T{\times }_3\widehat{\mathbf{E}}{\parallel }_F^2}{{\parallel \mathcal{X}\parallel }_F^2}$ and ${\gamma }_2=\frac{\parallel \mathcal{N}{\times }_3{\widehat{\mathbf{E}}}^T{\times }_3\widehat{\mathbf{E}}{\parallel }_F^2}{{\parallel \mathcal{N}\parallel }_F^2}$.

, where we show two metrics, namely ${\gamma }_1=\parallel \mathcal{X}{\times }_3{\widehat{\mathbf{E}}}^T{\times }_3\widehat{\mathbf{E}}{\parallel }_F^2/{\parallel \mathcal{X}\parallel }_F^2$ and ${\gamma }_2=\parallel \mathcal{N}{\times }_3{\widehat{\mathbf{E}}}^T{\times }_3\widehat{\mathbf{E}}{\parallel }_F^2/{\parallel \mathcal{N}\parallel }_F^2$, measuring the relative power of the clean spectral vectors and noise lying in the estimated subspace, respectively.

In Case 1, two images were corrupted only by Gaussian noise, and values of ${\gamma }_1$ for all methods are higher than 0.9993, implying that all the methods can estimate a subspace representing image signal very well. Case 3 shows results similar to Case 1; that is, all the methods achieve relatively high ${\gamma }_1$. We conclude that when noise is Gaussian distributed or a mixture of Gaussian and “Salt and Pepper” noise, signal subspace learning is not a challenging task and we can simply use SVD. The reason is that Gaussian noise and “Salt and Pepper” noise are randomly and uniformly distributed in each channel; thus, the noise increases singular values in the direction of each eigenvector almost uniformly and does not change the order of the singular values of clean images. The signal subspace is approximately spanned by p singular vectors of the noisy image corresponding to the largest p singular values.

However, subspace learning from noisy images in Case 2 and Case 4 is challenging and we cannot simply use SVD because the stripe noise usually exist in specific channels, instead of being uniformly distributed in all channels. The stripe noise will significantly increase the variances of those channels affected by stripe noise. Due to the fact that SVD tends to learn a subspace representing information in channels with high variances, SVD is not suitable for Case 2 and Case 4. This inference is consistent with the results in Table 4, where SVD obtains the worst results in terms of ${\gamma }_1$ and ${\gamma }_2$ in case 2 and case 4. If we only focus on the table rows highlighted with gray color, the proposed learning method performs better than others. Among the compared methods, RPCA learns a subspace representing low-rank signal and excluding sparse noise. However, stripe noise is also low-rank [50], and RPCA cannot separate an image signal and stripes using only a low-rank regularization. HySime is conceived based on an assumption of Gaussian noise and is clearly not suitable for mixed noise. L1HyMixDe performs median filtering band by band, exploiting spatial correlation, whereas HySuDeep estimates a coarse image using Hampel filtering pixel by pixel, exploiting spectral correlation. As shown in gray rows of Table 4, HySuDeep uniformly provides the best performances.

4.4. Analysis of Regularization Parameters

There are two regularization parameters—namely, ${\lambda }_1$ and ${\lambda }_2$—in the objective function of the proposed HySuDeep method. Controlling the trade-off between Gaussian noise reduction and detail preservation, the parameter ${\lambda }_1$ is related to the standard deviation of Gaussian noise, which can be estimated via (25).

The second parameter ${\lambda }_2$ controls the sparsity of estimation of sparse noise $\mathcal{S}$. The setting of ${\lambda }_2$ should depend on the intensity of sparse noise. Figure 9 depicts the denoising performance of HySuDeep in the Pavia University image and in the Washington DC Mall image as a function of the value of parameter ${\lambda }_2$. When the value of ${\lambda }_2$ is set to {1, 3, 5, 7, 9}, denoising results are acceptable in both images. Therefore, we simply set ${\lambda }_2$ to 3 in all experiments in this paper.

4.5. Numerical Convergence of the HySuDeep

Since CNN-regularization is incorporated into (10), the proposed HySuDeep cannot be considered as a convex optimization problem; thus, its theoretical convergence is not guaranteed. However, its numerical convergence is systematically observed when the augmented Lagrangian parameters are set to ${\mathsf{\mu }}_i=1(i=1,2,3)$. As shown in Figure 10, the relative change converges to near zero after 15 iterations, implying the convergence of the proposed method can be numerically guaranteed.

4.6. Application in Hyperspectral Unmixing

Hyperspectral denoising is usually implemented as a preprocessing step for subsequent applications. This subsection takes hyperspectral unmixing as an example to evaluate whether the image denoising step has a positive impact on the performance of subsequent applications. Spectral unmixing mainly involves two stages: (i) identifying materials present in the scene (termed endmembers) and (ii) estimating the fraction (or abundance) of each material in each pixel. The hyperspectral unmixing task was chosen considering that the endmember extraction step can be very sensitive to sparse noise.

We first performed denoising on the Terrain image using FastHyDe, NAILRMA, SSTV, LRMR, LRTF-DFR, L1HyMixDe, and HySuDeep. The mean PSNR (MPSNR), mean SSIM (MSSIM), and mean FSIM (MFSIM) of the proposed and compared methods on Terrain image are presented in Table 3, where we have highlighted the best results in bold. Then, we unmixed the spectra of the denoised images by using vertex component analysis (VCA) [53] to estimate the endmembers and used fully constrained least squares (FCLS) [54] to estimate the abundances. Two metrics were computed: the normalized mean square error (NMSE) of endmembers $\mathbf{A}$ and abundances $\mathbf{S}$, denoted as ${\mathrm{NMSE}}_{\mathrm{A}}$ and ${\mathrm{NMSE}}_{\mathrm{S}}$, respectively. As reported in Table 3, LRMR, LRTF-DFR, L1HyMixDe, and HySuDeep obtain a lower NMSE of endmembers than the counterpart image without denoising processing. For images denoised by FastHyDe, NAILRMA, and SSTV, the errors in the estimates of endmembers directly result in high errors in the estimates of abundances. Among LRMR, LRTF-DFR, L1HyMixDe, and HySuDeep, the proposed HySuDeep leads to the lowest NMSE of abundances.

5. Experimental Results for Real Images

We evaluate the performance of the proposed method on two real HSI datasets, as shown in Figure 4d,e.

5.1. Hyperion Cuprite Dataset

The Cuprite HSI was captured at Cuprite, NV, USA, by the Hyperion sensor, which divides the spectrum from 355 nm to 2577 nm into 242 channels with a spectral resolution of 10 nm. The spatial resolution of the image is 30 m. A subregion of size $240\times 178$ pixels with 177 spectral channels (after removing very low SNR channels) was cropped for the test. Four bands are shown in Figure 11, where bands 93, 133, and 134 were corrupted by severe stripes, and band 156 was affected by a dead-line. The full-resolution images of Figure 11 will be available at https://github.com/LinaZhuang, (accessed on 8 October 2021). All methods (except NAILRMA) are able to remove the dead-line in band 156. For severe stripes in the other three bands, we can see that FastHyDe, LRMR, L1HyMixDe, and HySuDeep achieved good restoration results while obvious stripes remained within the results of other methods.

5.2. Tiangong-1 Dataset

The Tiangong-1 dataset was acquired over an area of Qinghai Province, China in May 2013, by a sensor placed in the Tiangong-1 imager, which has a 75-band push broom scanner with nominal bandwidths of 23 nm short wave infrared (SWIR), covering from 800 nm to 2500 nm. A subregion image of size $351\times 253$ pixels was tested. Three bands displaying strong noise are shown in Figure 12, where band 1 has severe cross-track illumination error, and bands 36 and 65 contain obvious stripes. The full-resolution images of Figure 12 will be available at https://github.com/LinaZhuang (accessed on 8 October 2021). Comparing the images before denoising and after denoising, we can see that FastHyDe, L1HyMixDe, and HySuDeep can correct the illumination error in band 1. In bands 36 and 65, FastHyDe, LRMR, L1HyMixDe, and HySuDeep can alleviate stripe noise. If we focus on band 65, our HySuDeep method obtains the visually best result considering the restoration of the areas in the red circles.

6. Conclusions

This paper introduces a hyperspectral mixed noise removal method, HySuDeep, by exploiting two important characteristics of HSIs. HySuDeep takes advantage of the spectral low-rankness of HSIs by representing clean spectral vectors in a low-dimensional subspace, which significantly improves the conditioning of the denoising problem. The spatial correlation of HSIs is exploited by adding CNN regularization for eigenimages. Although the network was trained using grayscale images acquired from commercial cameras, without using remote sensing images, the network still achieves impressive performance for HSI denoising. The reason is that both kinds of images are natural images, sharing the same properties, such as local and non-local similarity, and piece-wise smoothness. Therefore, the image prior learned by the network from grayscale images is also applicable to HSIs. Experimental results on both simulated and real HSIs show the superiority of the proposed method for mixed noise in HSIs.