Hyperspectral Image Denoising with Partially Orthogonal Matrix Vector Tensor Factorization

In this paper, scalar, vector, matrix and tensor are denoted by lowercase letters $x$, boldface lowercase letters $\mathbf{x}$ and boldface capital letters $\mathbf{X}$, and calligraphic letters $\mathcal{X}$, respectively.

The $\ell_{1}$ norm and the Frobenius norm of a matrix $\mathbf{X}$ are defined as:

and

respectively, where $x_{i_{1},i_{2}}$ is the $(i_{1},i_{2})$-th entry of $\mathbf{X}$.

The inner product of two tensors $\mathcal{X},\mathcal{Y}\in\mathbb{R}^{I_{1}\times I_{2}\times\cdots\times I_{N}}$ is defined as:

The cross product of two tensors $\mathcal{X}\in\mathbb{R}^{I_{1}\times I_{2}\times\cdots\times I_{N}}$ and $\mathcal{Y}\in\mathbb{R}^{J_{1}\times J_{2}\times\cdots\times J_{M}}$ is defined as:

The mode-$n$ product of an $N$-th order tensor $\mathcal{X}\in\mathbb{R}^{I_{1}\times I_{2}\times\cdots\times I_{N}}$ with a matrix $\mathbf{U}\in\mathbb{R}^{J_{n}\times I_{n}}$ is a new tensor of size $I_{1}\times\cdots\times I_{n-1}\times J_{n}\times I_{n+1}\times\cdots\times I_{N}$, which can be denoted as follows:

The mode-$n$ unfolding of an $N$-th order tensor $\mathcal{X}\in\mathbb{R}^{I_{1}\times I_{2}\times\cdots\times I_{N}}$ is expressed as $\mathbf{A}_{(n)}$. The mode-$n$ unfolding operator arranges the $n$-th mode of $\mathcal{A}$ as the row while the rest modes as the column of the mode-$n$ unfolding matrix. Mathematically, the elements of $\mathbf{A}_{(n)}$ satisfy

where $j=\overline{i_{1},\cdots,i_{n-1},i_{n+1},\cdots,i_{N}}$.

For a third-order hyperspectral data $\mathcal{X}\in\mathbb{R}^{I\times J\times K}$, the total variation of HSI is denoted by [45, 46].:

where $\mathfrak{D}(\cdot)=\left[\mathfrak{D}_{h}(\cdot);\mathfrak{D}_{v}(\cdot);\mathfrak{D}_{z}(\cdot)\right]$ is a three-dimensional difference operator, $\mathfrak{D}_{h}$, $\mathfrak{D}_{v}$, $\mathfrak{D}_{z}$ are first-order finite difference operators along three different ways of hyperspectral data, which can be defined as:

The MVTF decomposes a tensor $\mathcal{X}\in\mathbb{R}^{I\times J\times K}$ into the sum of several component tensors, and each component tensor can be written in the cross product form of a matrix $\mathbf{G}_{r}\in\mathbb{R}^{I\times J}$ and a vector $\mathbf{c}_{r}\in\mathbb{R}^{K}$, which can be denoted as [39]:

where $\mathcal{G}\in\mathbb{R}^{I\times J\times R}$, $\mathbf{G}_{r}$ is the $r$-th slice of $\mathcal{G}$ and $\mathbf{C}\in\mathbb{R}^{K\times R}$.

In linear spectral mixture model of hyperspectral image, $\mathbf{c}_{r}$ can be modeled as the $r$-th endmember spectrum, and abundance matrix $\mathbf{G}_{r}$ represents the spatial information of $r$-th endmember. Due to the strong spatial correlation of HSI, each abundance matrix is always low rank.

Fixing other variables to update $\mathcal{G}$, the optimization model can be rewritten as:

where $\mathcal{G}(:,:,r)=\mathbf{G}_{r},r=1,\cdots,R$. The optimization problem (16) can be equivalent to

where $\mathcal{M}=(\mathcal{X}+\frac{\Lambda_{4}}{\beta_{4}})\times_{3}\mathbf{C}^{\text{T}}$. The detailed solutions of $\mathcal{G}$ are concluded in Algorithm 1, where SVD means singular value decomposition and the soft thresholding operator of $x$ is defined as:

Fixing $\mathbf{C}$ in equation (III), we can obtain the following equation:

Solving problem (III-B) is equal to :

Letting $\mathbf{M}=\mathbf{G}_{(3)}(\Lambda_{4(3)}^{\text{T}}+\beta_{4}\mathbf{X}_{(3)}^{\text{T}})$ and performing SVD on $\mathbf{M}$ as $\mathbf{M}=\mathbf{U}\bm{\Sigma}\mathbf{V}^{\text{T}}$, the solution of $\mathbf{C}$ can be obtained by

Fixing other variables to update $\mathcal{X}$, the optimization model (III) can be rewritten as:

Using the objective function in (III-C) to derive the $\mathcal{X}$, we can obtain the solution of $\mathcal{X}$ as follows:

Similarly, the optimization model (III) with respect to $\mathcal{Z}$ can be rewritten as:

The solution of this optimization problem can be transformed into the solution of the following linear system:

Considering the block circulant structure of the operator $\mathfrak{D}^{*}\mathfrak{D}$, it can be transformed into the Fourier domain and fast calculated. The fast computation of $\mathcal{Z}$ can be written as:

where $\mathcal{M}=\beta_{2}\mathcal{X}-\Lambda_{2}+\mathfrak{D}^{*}(\beta_{3}\mathbf{L}+\Lambda_{3})$; $\operatorname{fftn}$ and $\operatorname{ifftn}$ are 3D fast Fourier transform and 3D fast inverse Fourier transform, respectively.

To update $\mathbf{L}$, we fix all the rest variables in equation (III), and the sub-problem can be formulated as follows:

It can be transformed into:

The solution of $\mathbf{L}$ can be obtained by equation (18) as follows:

By fixing all the rest variables in equation (III), we can update $\mathcal{S}$ by the sub-problem:

Similarly to solving $\mathbf{L}$ above, this optimization model can be solved by:

Fixing the other variables in equation (III), the optimization model to update $\mathcal{N}$ can be formulated as follows:

It can be easily solved by:

The Lagrangian multiplier updating scheme can be as follows:

We summarize the algorithm for the SRLRTR with partially orthogonal MVTF for HSI denoising in Algorithm 2.

The convergence condition is reached when the relative error between two successive tensors $\mathcal{X}$ is smaller than a threshold, i.e. ${\left\|\mathcal{X}^{(k)}-\mathcal{X}^{(k+1)}\right\|_{\text{F}}^{2}}/{\left\|\mathcal{X}^{(k+1)}\right\|_{\text{F}}^{2}}\leq\varepsilon$ where $\mathcal{X}^{(k)}$ is the recovered image in $k$-th iteration and $\varepsilon$ is the threshold.

The main complexity of the proposed algorithm comes from the updating of $\mathcal{G}$. For a hyperspectral image $\mathcal{X}\in\mathbb{R}^{I\times J\times K}$, the most time-consuming part is from calculating $\mathcal{G}$, which needs to compute the SVD of $\mathcal{M}_{(:,:,r)}\in\mathbb{R}^{I\times J},r=1,\cdots,R$. Assuming $I=J$, the computational complexity is $O(I^{3})$. Therefore, the total complexity is $O(PRI^{3})$, where $P$ is the number of iterations.

1. Washington DC Mall (WDC): The Washington DC Mall dataset¹¹1https://engineering.purdue.edu/ biehl/MultiSpec/hyperspectral.html is provided with the permission of Spectral Information Technology Application Center of Virginia. The size of each scene is $1208\times 307$. The dataset includes 210 bands in the spectral range from 400 nm to 2400 nm. We use the 191 bands by removing bands where the atmosphere is opaque. And $256\times 256$ pixels are selected from each image, forming the HSI with the size of $256\times 256\times 191$. The false color image (R: 29, G: 19, B: 9) of WDC is shown in Fig. 3(a) .

2. Pavia University (PaviaU): This scene²²2https://rslab.ut.ac.ir/data is acquired by ROSIS sensor that covers the Pavia, northern Italy. The spatial size of the PaviaU is $610\times 340$ and there are 103 channels, ranging from 430 nm to 860 nm. And $300\times 300$ pixels are considered, forming the HSI with the size of $300\times 300\times 103$. The false color image (R: 29, G: 19, B: 9) of PaviaU is shown in Fig. 3(b).

3. EO-1 Hyperion datasets (EO-1): The EO-1 Hyperion hyperspectral dataset³³3https://www.usgs.gov/centers/eros/science/ covers an agricultural area of the state of Indiana, USA. It contains $1000\times 1000$ pixels and 242 bands spanning 350nm-2600nm. After removing water absorption bands, 166 bands are retained, and $200\times 200$ pixels are selected from original scene, forming the HSI with size of $200\times 200\times 166$. It is mainly contaminated by Gaussian noise, stripes and dead lines.The false color image (R: 2, G: 132, B: 136) of EO-1 is shown in Fig. 3(c).

4. HYDICE Urban Dataset (Urban): This scene is acquired from the real Hyperspectral Digital Imagery Collection Experiment (HYDICE) [47]. It contains $307\times 307$ pixels and there are 210 wavelengths ranging from 400 nm to 2500 nm. It is mainly contaminated by stripes and dead lines. The false color image (R: 2, G: 109, B: 207) of Urban is shown in Fig. 3(d).

To compare with the proposed method, we select 9 state-of-the-art algorithms and one representative wavelet denoising algorithm. The detailed descriptions are as follows:

• •

  SNLRSF [48]: a novel subspace method based on nonlocal low-rank and sparse constraints for mixed noise reduction of HSI.

• •

  NLRCPTD  [34]: a method based on nonlocal low-rank regularized CP decomposition with rank determination for HSI denoising.

• •

  3DATVLR  [12]: a method based on 3D anisotropic total variation and low rank constraints for mixed-noise removal of HSI.

• •

  LRTDTV  [31]: a method based on total variation and low rank Tucker decomposition for HSI mixed noise removal.

• •

  RSLRNTF  [49]: a novel method based on matrix-vector tensor factorization with low rank and nonnegative constrains for HSI denoising.

• •

  SSTVLRTF [50]: a method based on spatial-spectral total variation and low rank t-SVD factorization to remove mixed noise in HSI.

• •

  TVLASDA [51]: a method based on spectral difference-induced total variation and low rank approximation for HSI denoising.

• •

  HyRes [52]: a method based on sparse low-rank model, which automatically tunes parameter.

• •

  BM4D [53]: a wavelet method based on nonlocal self-similarities, which achieves great success in nature image denoising.

The parameters of the comparison algorithms are set according to the values in its corresponding literatures.

In addition, to evaluate the performance quantitatively, four metrics are used to evaluate the image denoising quality, including MPSNR, MSSI, ERGAS and CPU time. The details description are as follows:

• •

  $\text{MPSNR}=\frac{1}{K}\sum_{k=1}^{K}\text{PSNR}_{k}$

• •

  $\text{MSSIM}=\frac{1}{K}\sum_{k=1}^{K}\text{SSIM}_{k}$

• •

  $\text{ERGAS}=\sqrt{\frac{1}{K}\sum_{k=1}^{K}{\left\|\hat{\mathcal{X}}(:,:,k)-\mathcal{X}(:,:,k)\right\|_{\text{F}}^{2}}/{\mu_{k}^{2}}}$

where $K$ is number of spectral bands. $\mu_{k}$ is the average value of $\mathcal{X}(:,:,k)$ and $\hat{\mathcal{X}}(:,:,k)$ is the recovered data. All the experiments are conducted using MatLab 2019b on a computer with 2.4 GHz Quad-Core Intel Core i5 processor and 16 GB RAM.

WDC and PaviaU hyperspectral images are used as simulated datasets in this experiment. To testify the proposed method, we add Gaussian noise, salt-and-pepper noise, dead line and strip noise to the simulated ground-truth datasets. Four groups of experiments are performed, and the experimental settings are as follows:

Case 1: Gaussian noise + impulse noise

In this group of experiments, the noise intensities of different bands are the same. In each frequency band, the same zero mean Gaussian noise and the same percentage of impulse noise are added. The variance of Gaussian noise is 0.2, and the percentage of impulse noise is 0.2.

Case 2: Gaussian noise +dead lines

In this group of experiments, Gaussian noise is added to each band in the same way where the noise level is equal to 0.15. In addition, a dead line is added to the spectral segment from band 41 to band 100 with the number of stripes randomly chosen from 3 to 10. The width of the stripes were randomly generated from one line to three lines.

Case 3: Gaussian noise + impulse noise+dead lines

In this group, the noise intensities of different bands are not the same. Gaussian noise with zero mean and 0.05 variance, is added to each band. The percentage of impulse noise is 0.1. The dead lines are added randomly between 41 and 100 segments, the number of stripes is selected randomly from 3 to 10, and the width ranges randomly from one line to three lines.

Case 4: Gaussian noise + impulse noise+dead lines+stripe

Based on the 3rd group of experiments, this group randomly adds stripe noise between band 101 and band 190 and the number of random stripe ranges from 20 to 40 for WDC data. For PaviaU, the dead lines are added randomly between 41 and 60 segments, the number of stripes is selected randomly from 3 to 10, and the width ranges randomly from one line to three lines. In addition, stripe noise is added between band 61 and band 100 when the number of random stripe ranges from 20 to 40.

In the proposed method, we set the parameters $\lambda_{\text{TV}}=0.0002,\lambda_{\mathcal{S}}=0.02,\lambda_{\mathcal{N}}=0.1,\lambda_{\mathcal{G}}=0.1$ and $R=5$ for WDC and PaviaU according to the highest value of MPSNR. The detailed information is put on the discussion part.

Table I shows the results of experiments based on simulated data in terms of MPSNR, MSSIM, ERGAS and CPU time. It can be observed, in most cases, the proposed algorithm is superior to all other methods in terms of the MPSNR, MSSIM and ERGA. Specifically, our algorithm shows advantage on removing mixed Gaussian noise and impulse noise in case 1. It means that our algorithm can deal with Gaussian noise and impulse noise well. Similarly, in case 3 and case 4, when the dead line and stripe noise are added, our proposed method is also superior to existing state-of-art methods in terms of recovery quality. This may imply the superiority of the proposed method in processing mixed noise in different conditions. However, in case 2 where it contains Gaussian noise and dead lines, SNLRSF, NLRCPTD and HyRes performs better.

To explain the recovery results more vividly, Fig. 4 and Fig. 5 present comparative results for the 121st band of WDC HSI in case 1 and the 85th band of PaviaU HSI in case 3, respectively. Fig. 4 shows the recovered images by different methods on mixed Gaussian noise and impulse noise condition. We could observe that the resolution of recovered image in our proposed one is clearer than that of others. In addition, Fig. 5 shows the recovered images on mixed Gaussian noise, impulse noise and dead lines condition. We could see that recovered images of BM4D and RSLRTF still suffer from the dead lines noise.

To give the details of results on each spectral segment, Fig. 6 and Fig. 7 provide the PSNR of each spectral segment in case 4 for WDC and PaviaU data, respectively. In case 4, we add different kinds of noises in different spectral segments. From Fig. 6, we can see the PSNR values change smooth in the 1-41 spectral segments where Gaussian noise and impulse noise are added while the PSNR values fluctuate in 50-190 spectral segments where dead line and stripe noise are added. In addition, for each scene, the recovery performance of the proposed one is better than that of others with respect to PSNR. It indicates that the proposed method is superior in removing mixed noises. In Fig. 7, the PSNR values of the proposed method is better than the others in almost every spectral segment. It means that it can remove noise and retain the structural characteristics of HSI better.

In conclusion, the proposed method performs better in the denoising of Gaussian noise and impulse noise. In the processing of stripe noise, it can also have a leading position among the compared methods.

We choose EO-01 and Urban datasets in this experiment. Both scenes are seriously contaminated by the Gaussian noise, stripes and dead lines. In this case, we set the parameters $\lambda_{\text{TV}}=0.00001,\lambda_{\mathcal{S}}=0.013,\lambda_{\mathcal{N}}=0.1,\lambda_{\mathcal{G}}=0.1,R=2$ for Urban and EO-01data.

Fig. 8, Fig. 9 and Fig 10 show the denoising resulting of EO-1 hyperspectral image in bands 68, 96 and 132. It can be observed in Fig. 8, the original data are heavily contaminated by stripes noise. Our proposed algorithm can recover the HSI while local details and structural information of the HSI are preserved. However, HyRes and BM4D fail to recover the HSI in stripes noise. This is mainly because they cannot model the sparse noise. Similarly, HyRes and BM4D performs badly in Fig. 9 and Fig. 10. In addition, it can be observed that the proposed method obtains better performance in terms of image resolution.

Fig. 11, Fig. 12 and Fig. 13 show the 109th, 151st and 207th spectrum segments of Urban hyperspectral image after denoising. Fig. 11 shows that the TV based methods including LRTDTV, 3DATVLR, SSTVLRTF, TVLASDS and our proposed one have a good recovery performance on stripe noise condition. Meanwhile, the nonlocal based methods such as SNLRSF and NLRCPTD can recover the image. However, with the noise increase as shown in Fig. 12, some TV based methods including LRTDTV, 3DATVLR and SSTVLRTF are over-smoothness to recover the HSI. Instead, SNLRSF, TVLASDS and our proposed can successfully recover the clean image from the severely noise condition. Specially, as shown in Fig. 13, the scene is severely polluted by serious Gaussian noise, stripe noise and dead line, only SNLRSF, TVLASDS and our proposed method have a good recovery performance.

In addition, table II shows the CPU time of our algorithm and state-of-art algorithms on real data. It is noticed that due to the poor performance of the BM4D and HyRes methods in the above experiments, we do not consider the CPU time of them for comparison. As seen in table II, 3DATVLR is the fastest one compared with others, followed by our proposed one. But the recovery performance of 3DATVLR is worse than that of SRLRTR. Meanwhile, compared with TVLASDS and SNLRSF which have a good performance on severely noisy condition, the proposed one is faster than them. It implies that the partially orthogonal MVTF method can be more efficient and suitable to explore HSI information.

We introduce a new smooth and robust low rank partially orthogonal MVTF method for denoising of HSI. The partially orthogonal MVTF is crucial to model the clean HSI for the proposed algorithm. Compared with traditional MVTF including orthogonal MVTF [54] and nonnegative MVTF [39], the proposed one is different in two aspects.

1. The traditional orthogonal MVTF in [54] is defined as:

where $\mathbf{G}_{r}=\mathbf{A}_{r}\mathbf{D}_{r}\mathbf{B}_{r}^{\text{T}}$ is the SVD of $\mathbf{G}_{r}$ and $\mathbf{c}_{r}$ is under the unit-norm constraint. In this case, the matrices $\mathbf{A}_{r}$ and $\mathbf{B}_{r}$ are orthogonal. Compared with traditional orthogonal MVTF, the proposed one constrains $\mathbf{C}^{\text{T}}\mathbf{C}=\mathbf{I}_{R}$, $\mathbf{C}(:,r)=\mathbf{c}_{r},r=1,\cdots,R$, resulting in partially orthogonal MVTF. The partially orthogonal MVTF has a good physical interpretation under the assumption of linear spectral mixture model for HSI, when $\mathbf{c}_{r}$ is considered as the $r$-th endmember spectrum, and endmember spectrums are uncorrelated.

2. The nonnegative MVTF in [39] is denoted by:

where $\mathbf{G}_{r}$ is of rank $L_{r},r=1,\cdots,R$, and can be factorized as $\mathbf{A}_{r}\mathbf{B}_{r}^{\text{T}}$ with predefined rank $L_{r}$. In HSI, the abundance matrix is low rank.

This low rank model has shown its superiority on HSI unmixing [39], where $\mathbf{c}_{r}$ is considered as an endmember, and matrix $\mathbf{G}_{r}$ is the corresponding abundance map. In [49], RSLRNTF uses this model for HSI denoising. In these cases, the ranks $L_{r}$, $r=1,\cdots,R$ are required to be defined in advance. In practice, accurate determination of the rank of abundance matrix $\mathbf{G}_{r}$ of the $r$-th endmember is impossible. Therefore, we directly minimize the rank of each abundance matrix under the assumption that endmember spectrums are irrelevant, resulting in low rank partially-orthogonal MVTF. In addition, compared with RSLRNTF [49], our proposed method shows a better recovery performance in experiments.

There are four parameters in equation (III), but we only need to tune three parameters because the parameters represent the relative weights of different terms in objective function. In this case, we fix $\lambda_{\mathcal{G}}$ and analyze the impacts of parameters $\lambda_{\mathcal{S}},\lambda_{\mathcal{N}},\lambda_{\text{TV}}$. In addition, the parameter $R$ in our algorithm has its physical meaning, which represents the number of endmembers (also called pure materials). In the following, we will take the experiments in case 4 on WDC datasets as an example and address how to choose these parameters.

1. The impact of parameters $\lambda_{\text{TV}}$, $\lambda_{\mathcal{S}}$ and $\lambda_{\mathcal{N}}$: For tuning the parameters $\lambda_{\text{TV}}$, $\lambda_{\mathcal{S}}$ and $\lambda_{\mathcal{N}}$, we set $\lambda_{\mathcal{G}}=0.1$. The other three parameters are related to the weights of TV regularizer, sparse noise (i.e., impulse noise, stripes, and dead lines), and Gaussian noise, respectively. Fig. 15 shows the MPSNR values as a function with respect to $\lambda_{\text{TV}}$ and $\lambda_{\mathcal{S}}$ for the proposed algorithm with $\lambda_{\mathcal{N}}$ chosen from the set {0.01, 0.018, 0.032, 0.056, 0.1}. $\lambda_{\text{TV}}$ is selected from the set {0.0001, 0.00017, 0.00046, 0.0013, 0.0036, 0.01} and $\lambda_{\mathcal{S}}$ is selected from the set {0.001, 0.0028, 0.0077, 0.02, 0.05, 0.1}. It is obvious that the smaller $\lambda_{\text{TV}}$ and the larger $\lambda_{\mathcal{N}}$, the value of MPSNR becomes larger. When $\lambda_{\mathcal{S}}=0.02$, the proposed method can reach the peak of MPSNR.

2. The impact of parameters $R$:

The parameter $R$ in our algorithm represents the number of endmembers. Fig. 16 shows the MPSNR value as a function of $R$, where $R$ changes from 2 to 15 with step size 1. From it, we can see when $R=5$, the value of MPSNR arrives at its peak.