Mathematical Problems in Engineering

Mathematical Problems in Engineering / 2021 / Article

Research Article | Open Access

Volume 2021 |Article ID 5535169 | https://doi.org/10.1155/2021/5535169

Shuo Wang, Zhibin Zhu, Ruwen Zhao, Benxin Zhang, "Hyperspectral Image Denoising via Nonconvex Logarithmic Penalty", Mathematical Problems in Engineering, vol. 2021, Article ID 5535169, 18 pages, 2021. https://doi.org/10.1155/2021/5535169

Hyperspectral Image Denoising via Nonconvex Logarithmic Penalty

Academic Editor: Ana C. Teodoro
Received05 Jan 2021
Revised22 Mar 2021
Accepted08 Jun 2021
Published24 Jun 2021

Abstract

Hyperspectral images (HSIs) can help deliver more reliable representations of real scenes than traditional images and enhance the performance of many computer vision tasks. However, in real cases, an HSI is often degraded by a mixture of various types of noise, including Gaussian noise and impulse noise. In this paper, we propose a logarithmic nonconvex regularization model for HSI mixed noise removal. The logarithmic penalty function can approximate the tensor fibered rank more accurately and treats singular values differently. An alternating direction method of multipliers (ADMM) is also presented to solve the optimization problem, and each subproblem within ADMM is proven to have a closed-form solution. The experimental results demonstrate the effectiveness of the proposed method.

1. Introduction

A hyperspectral image (HSI) consists of multiple discrete bands at specific frequencies. It can deliver additional information that the human eye fails to capture for real scenes and has been attracting much interest from researchers in a wide range of application fields, such as land use analysis, environmental monitoring, and field surveillance [13]. However, HSIs always suffer from various degradations, such as Gaussian noise, impulse noise, and random noise, which can affect the subsequent image processing, such as unmixing, classification, and target detection [4, 5]. Improving the HSI quality merely through a hardware scheme is unsustainable and impractical. Therefore, it is natural to introduce image processing-based approaches to obtain a high-quality HSI before subsequent applications.

The numerical optimization algorithm plays an important role in HSI denoising, such as Liu et al. [6] proposed a two-step wavelet-domain estimation method to extract the noise map, and Lu et al. [7] presented some representative high-order variational models in the context of image denoising. From the perspective of the prior data format, we classify the existing HSI restoration methods into three categories: (1) 1D vector-based sparse representation methods [813]; (2) 2D matrix-based low-rank matrix recovery (LRMR) methods [1421]; and, (3) 3D tensor-based approximation methods [2233]. Although the existing works have made significant progress in HSI restoration, there are still drawbacks that need to be improved, such as when the multidimensional HSI data is transformed into a vector or matrix, it usually breaks the spectral-spatial structural correlation. The tensor low-rankness characterization of an HSI is expected to explore the global correlation and preserve the intrinsic structural information. The tensors’ recovery under a limited number of measurements is an important problem that arises in a variety of applications, such as computer vision [3436]. Based on low tubal-rank tensor recovery, Zhang et al. [37] proposed an HSI mixed noise removal model. However, the framework of the tensor singular value decomposition (t-SVD) lacks flexibility for handling different correlations along the different modes of HSIs, leading to suboptimal denoising performance. Then, Zheng et al. [38] proposed an HSI mixed noise removal model with tensor fibered rank, which is based on the mode-k t-SVD. Moreover, Zheng et al. [38] also proposed a three-directional tensor nuclear norm (3DTNN) as its convex relaxation to provide an efficient numerical solution and a three-directional log-based tensor nuclear norm (3DLogTNN) as its nonconvex relaxation to promote the low rankness of the solution. Compared to 3DTNN, 3DLogTNN has two advantages. First, it is a closer approximation to the fibered rank than 3DTNN. Second, by using the sum of the log function of singular values, 3DLogTNN can better approximate to the fibered rank than 3DTNN.

It is well known that suitable nonconvex penalty functions induce sparsity among the singular values more effectively. However, the use of nonconvex penalty functions will lead to nonconvex optimization problems. Then, it suffers from numerous issues such as spurious local minima in the subproblem, for example, 3DLogTNN in [38]. To avoid the intrinsic difficulties, we introduce a new nonconvex logarithmic regularization model, which allows the use of nonconvex penalty function while maintaining convexity of the subproblem within ADMM. Also, the new model can provide a good approximation for the fibered rank and preserve the major information well.

The rest of the paper is structured as follows. Section 2 presents some notations, explains t-SVD and defines mode-k t-SVD. Section 3 introduces the proposed ADMM based on the parametric penalty function for HSI denoising. The experimental results and analysis are reported in Section 4. Finally, the conclusion is given in Section 5.

2. Brief Overview of Tensor Singular Value Decomposition

In this section, we first describe the notations used throughout the paper and then introduce the tensor decomposition proposed in [3941] and mode-k t-SVD proposed in [38].

2.1. Notation and Indexing

For a third-order tensor , denote as the horizontal, lateral, and frontal slices, respectively. denote the mode-1, mode-2, and mode-3 fibers, respectively; Figures 1 and 2 show the horizontal, lateral, and frontal slides, denoted by , respectively, of a third-order tensor . Suppose ; we adopt the definition of the Frobenius norm of a tensor and the definition of the norm of a tensor, .

Let denote , that is, the frontal slice of . Then,

Similarly, , and

It can be seen in Figure 3, for , let ; then, are the frontal slices of tensor , where is computed by applying the Fast Fourier Transform (FFT) along each tube of , i.e., and .

2.2. t-SVD

In this section, we exploit the proposed t-SVD. A t-SVD interprets third-order tensors as linear operators on the space of oriented matrices [39]. Based on a t-SVD, O. Semerci exploited the decomposition, completion, and recovery of multilinear data [42], and Zhang et al. [34] proposed novel methods for completion and denoising of multilinear data and, as an application, considered 3D and 4D (color) video data completion and denoising.

Definition 1. (t-product). For and , the t-product is a tensor of size . For and ,The t-product is analogous to matrix multiplication except that circular convolution replaces the multiplication operations between the elements, which are represented by tubes.

Theorem 1. (t-SVD). For , the t-SVD of is given bywhere and are orthogonal tensors, is a rectangular diagonal tensor, and denotes the t-product.

Figure 4 illustrates the decomposition for the 3D case. Additionally, one can obtain this decomposition by computing matrix SVDs in the Fourier domain, see Algorithm 1.

However, when handling different correlations along different modes of an HSI, the t-SVD and the induced tubal rank lack flexibility. This inflexible HSI characterization usually does not have good denoising effects. Zheng et al. [38] proposed a novel tensor decomposition by generalizing the t-SVD to the mode-k t-SVD, which achieves a more flexible and accurate HSI low-rankness characterization.

Input: , .
for to do
end for
for to do
;
end for
for to do
;
end for
2.3. Mode-k t-SVD

In this section, we introduce the mode-k t-SVD and the related definitions.

For a third-order tensor , the mode-k block circulation operation is denoted aswhere is the mode-k slice of .

The mode-k block diagonalization operation and its inverse operation are defined as

The mode-k block vectorization operation and its inverse operation are defined as

Definition 2. (Mode-k t-product). For and , the mode-1 t-product is a tensor of size :For and , the mode-2 t-product is a tensor of size :For and , the mode-3 t-product is a tensor of size :

Definition 3. (Mode-k identity tensor). is the mode-k identity tensor, whose first mode-k slice is an identity matrix and other mode-k slices are all zeros.

Definition 4. (Mode-k conjugate transpose). For is the mode-k conjugate transpose of , which is obtained by transposing each of the mode-k slices and then reversing the order of transposed mode-k slices 2 through .

Definition 5. (Mode-k diagonal tensor). For , each of its mode-k slices is a diagonal matrix, and then, is a mode-k diagonal.

Definition 6. (Mode-k orthogonal tensor). Ifwhere the tensor is mode-k orthogonal.

Definition 7. (Tensor mode-k permutation). For , is the mode-k permutation of , whose mode-3 slice is the mode-k slice of , and its inverse operation is .

Theorem 2. (Mode-k t-SVD). The factorization of iswhere and are the mode-k orthogonal tensors and is the mode-k diagonal tensor.

Theorem 2 is proven in [38] (Th. 2).

Definition 8. (Tensor fibered rank). For , is the fibered rank of , whose element is the mode-k tensor fibered rank . , where the element of is the rank of the mode-k slice of .
The mode-k t-SVD can be efficiently obtained by computing a series of matrix SVDs in the Fourier domain and and achieves more flexible and accurate HSI low-rankness characterization (see Algorithm 2).

Input: , k
for to do
end for
output:

3. HSI Denoising Model and Its ADMM

In this section, we show our new denoising model and ADMM for solving the proposed model is also derived.

3.1. The Logarithmic Penalty Function

This article mainly proposes logarithmic penalty function [43] that serves as the model for the penalty function developed in the HSI denoising model below and is designed to have less bias. The logarithmic penalty is given bywhere controls the degree of nonconvexity of the penalty function. This function satisfies the following conditions:(A1)(A2)(A3)(A4)

The proximity operator associated with the nonconvex function is

In [4446], the authors prove that, for , the function in equation (14) is convex. Therefore, the proximity operator finds an optimal solution for convex minimization problem (14). The proximity operator associated with logarithmic penalty equation (13) is a continuous nonlinear threshold function with as the threshold value [44], namely,and is given by

3.2. HSI Denoising Model

Zhang et al. proposed an HSI mixed noise removal model based on low tubal-rank tensor recovery (LRTR) [37]. It can address the mixed noise in HSIs and decompose a noisy HSI into three parts, i.e., a low-tubal-rank part (the clean HSI), a Gaussian noise part, and a sparse noise part. Zheng [38] proposed an HSI denoising model based on a low-fibered-rank prior, and it is formulated aswhere is the observed noisy HSI, is Gaussian noise, is sparse noise, is the clean hyperspectral signal, which has the low-fibered-rank tensor property, and and are regularization parameters.

The Gaussian noise model (large degree of freedom) corresponds to a dense noise type [47, 48]. should not be low in rank. Otherwise, cannot be recovered from random noise. Due to the stochastic nature of Gaussian noise, it is assumed that there is no correlation (or a weak correlation) among noise components. Thus, the rank of is normally full and much larger than the rank of . Therefore, the low-fibered-rank component can be recovered from equation (17) by properly choosing the tuning parameters and .

Directly minimizing the tensor fibered rank is NP-hard. Based on logarithmic penalty equation (13), we propose that the following 3DNLogATNN HSI denoising model is commonly formulated:where and is the weight of the fibered rank. are set as in the HIS denoising model:where is the number of singular values of , is the singular values of , and is the mode-k slice of .

3.3. ADMM for Solving Model (18)

We use the ADMM to solve equation (18). By casting the three auxiliary variables , equation (18) can be rewritten as follows:

The augmented Lagrangian function of equation (18) iswhere and are the Lagrange multipliers and and are positive scalars. Now we can solve the problem within the ADMM framework.

With the other parameters fixed, can be updated by solvingwhich is equivalent to the following subproblem for :

To solve equation (23), we can rewrite it aswhere , and . From equation (16),where . Therefore,

Though model (20) is based on a nonconvex penalty function, the parameters of which are constrained to ensure the convexity of the subproblem .

With other parameters fixed, can be updated by solvingwhich is equivalent to the following subproblem:

It has the following closed-form solution:

With the other parameters fixed, can be updated by solvingwhich is equivalent to the following subproblem:and it has the following closed-form solution:

With the other parameters fixed, can be updated by solvingwhich is equivalent to the following subproblem:

Using the tensor soft thresholding operator, the following solution can be obtained [38]:

The Lagrange multipliers and can be updated by solving

Hence, the proposed algorithm for HSI denoising is summarized in Algorithm 3.

Input: The noisy HSI , parameters .
, and .
while not converged do
Update with equation (26), .
Update with equation (29).
Update with equation (32).
Update with equation (35).
Update with equation (36), .
Update with equation (37).
Let ; and .
Check the convergence condition:
.
end while
Output: Denoised HIS .

4. Experiment Results

To validate the effectiveness of the proposed method for HSI denoising, we perform experiments on simulated data and compare the experimental results both quantitatively and visually. The Washington DC Mall data, Pavia City Center data, and the Indian Pines data are used. In our experiments, the Washington DC Mall data uses only a subimage (191 bands and size of each band is ). The Pavia City Center data uses only a subimage (80 bands and size of each band is ). And, the synthetic data by Zhang et al. [37] was generated using the ground truth of the Indian Pines dataset; the size of the synthetic HSI was . In addition, to facilitate the numerical calculation and visualization, all the bands of the HSI are normalized into [0, 1], and they will be stretched to the original level after restoration.

The three evaluation measures are the mean peak signal-to-noise ratio (MPSNR), mean structure similarity (MSSIM), and spectral angle mapping (SAM). The three metrics are defined as follows to measure the quality of the denoised result:where represents the image size in the space, represents the original image, represents the distortion image, is a pixel that can achieve the maximum value, and is the number of PSNR:where is a constant for is the same as and represent the and standard deviations, respectively, represent the original image and the distorted image, respectively, and represent the jth local window contents, and is the number of local windows:where is the number of wavelengths, and represent the spectrum vectors, and

The PSNR and structural similarity index measure (SSIM) are two conventional perceptual quality indexes (PQIs) in image processing and computer vision. They evaluate the similarities between the target image and the reference image based on the mean square error (MSE) and structural consistency. The larger these two measures are, the closer the target HSI is to the reference HSI. The SAM is a physically based spectral classification method that uses an N-dimensional angle to match pixels to reference spectra. Different from the former two measures, the smaller the SAM is, the more similar the target HSI is to the reference HSI.

Real HSIs usually include several different types of noise. To simulate real-noise scenarios, we consider the variance of the Gaussian noise and the variance of the impulse noise . We use statistical structures to simulate different types of noise, including independent and identically distributed (i.i.d.) and non-i.i.d. noise, which are elaborated as follows:(1)Case 1 (non-i.i.d. Gaussian noise): all the bands of the test dataset are contaminated by zero-mean Gaussian noise with different intensities. The variance in the Gaussian noise is randomly changed from U (0.1, 0.2) and U (0.55, 0.65).(2)Case 2 (non-i.i.d. impulse noise): in this case, all bands are contaminated by impulse noise with different ratios. The ratios of impulse noise are randomly changed from U (0.35, 0.45), U (0.45, 0.55), and U (0.55, 0.65).(3)Case 3 (i.i.d. Gaussian + i.i.d. impulse noise): in this case, all bands are corrupted by zero-mean Gaussian noise and impulse noise. The variance in the Gaussian noise is 0.3, and the ratio of the impulse noise is 0.1; the variance in the Gaussian noise is 0.1, and the ratio of the impulse noise is 0.4; the variance in the Gaussian noise is 0.3, and the ratio of the impulse noise is 0.5.(4)Case 4 (non-i.i.d. Gaussian + i.i.d. impulse noise): in this case, all bands are corrupted by zero-mean Gaussian noise and impulse noise with different intensities. The variance in the Gaussian noise is randomly changed from U (0.3, 0.4), and the ratio of the impulse noise is 0.2; the variance in the Gaussian noise is randomly changed from U (0.2, 0.3), and the ratio of impulse noise is 0.3; the variance in the Gaussian noise is randomly changed from U (0.4, 0.5), and the ratio of impulse noise is 0.1.(5)Case 5 (i.i.d. Gaussian + non-i.i.d. impulse noise): in this case, all bands are corrupted by zero-mean Gaussian noise and impulse noise with different intensities. The variance in the Gaussian noise is 0.1, and the ratios of impulse noise are randomly varied from U (0.5, 0.6); the variance in the Gaussian noise is 0.3, and the ratios of impulse noise are randomly varied from U (0.3, 0.4); the variance in the Gaussian noise is 0.2, and the ratios of impulse noise are randomly varied from U (0.4, 0.5).(6)Case 6 (non-i.i.d. Gaussian + non-i.i.d. impulse noise): in this case, all bands are corrupted by zero-mean Gaussian noise and impulse noise with different intensities. The variance in the Gaussian noise is randomly changed from U (0.2, 0.3), U (0.1, 0.2), and U (0.4, 0.5), and the ratios of impulse noise are randomly varied from U (0.2, 0.3), U (0.3, 0.4), U (0.1, 0.2), and U (0.4, 0.5).

Parameter setting: we analyze the parameters involved in the proposed method on HSIs’ Washington DC Mall, Pavia City Center, and Indian Pines, i.e., the weight , the regularization parameters and , the threshold parameter , the penalty parameter , and a constant in 3DLogATNN. In all the following experiments, the parameters in these compared methods were manually adjusted according to their default strategies.

The regularization parameter for 3DLogATNN: it is easy to see that is the parameter used to restrict the sparsity of the Gaussian noise. We set , where is the standard deviation of Gaussian noise and is a tuning parameter. The results were based on the simulated data experiment in case 1-1. Figure 5 shows the restoration results as varied in the set {0.0001, 0.0005, 0.0009, 0.001, 0.002, 0.003, 0.004, 0.005, 0.007, 0.009, 0.11, 0.5, 0.15}. It can be clearly observed from this figure that the results of the 3DLogATNN solver are relatively stable in terms of both MPSNR and MSSIM values, with the value of changing from 0.002 to 0.003. Therefore, we suggest the use of in all the simulated data experiments.

The regularization parameter for 3DLogATNN: it is easy to see that is the parameter used to restrict the sparsity of the impulse noise. We set , whereand is a tuning parameter. The results were based on the simulated data experiment in case 2-1. Figure 6 shows the restoration results as varied in the set {0.05, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9}. It can be clearly observed from this figure that the results of the 3DLogATNN solver are relatively stable in terms of both MPSNR and MSSIM values, with the value of changing from 0.4 to 0.5. Therefore, we suggest the use of in all the simulated data experiments.

The constant for 3DLogATNN: the parameter controls the degree of nonconvexity of the penalty function. The results were based on the simulated data experiment in case 2-1. Figure 7 shows the restoration results as varied in the set {0.01, 0.04, 0.05, 0.06, 0.1, 0.5, 1, 5, 10}. It can be clearly observed from this figure that the results of the 3DLogATNN solver are relatively stable in terms of both MPSNR and MSSIM values, with the value of changing from 0.05 to 0.06. Therefore, we suggest the use of in all the simulated data experiments.

We adjust the parameters to achieve the best visual result, and the parameter setting is presented in Table 1.


MethodData

3DTNNWashington DC Mall
Pavia City Center
Indian Pines

3DLogTNNWashington DC Mall50
Pavia City Center80
Indian Pines50

3DLogATNNWashington DC Mall0.05
Pavia City Center
Indian Pines

Compared with the state-of-the-art methods, including TRPCA + BM4D [36, 49], LRMR [37], LRTR [37], LRTDTV [50], and NMoG [51], on low-rank matrix/tensor approximation and noise modeling, the extensive experimental results demonstrate that the 3DTNN and 3DLogTNN [38] methods are better at removing the mixed noise. Therefore, the denoising results of the proposed method are quantitatively and visually compared with two state-of-the-art HSI denoising methods, i.e., 3DTNN and 3DLogTNN. The denoising results of all the methods in six cases are shown in Table 2. Three typical bands of the denoised HSIs in the mixture noise case obtained with different methods are shown in Figures 810. Figure 8 shows the denoising results at band 71 of the Washington DC Mall HSI, Figure 9 shows the denoising results at band 52 of the Pavia City Center HSI, and Figure 10 shows the denoising results at band 28 of the Indian Pines HSI. It can be seen that the proposed 3DLogATNN can effectively remove the mixed noise and preserve the detailed information of the original image. The proposed method obtains the best visual quality by removing all the mixture noise and preserving the details well. Table 2 shows that the 3DLogATNN method converges faster than the 3DLogTNN-based method on all the Washington DC Mall, Pavia City Center data, and Indian Pines data, and our method outperforms the compared ones for the Pavia City Center data. Besides, 3DLogATNN is more stable than the other two algorithms, because it can get the best results in most cases.


CaseDataGaussian noiseImpulse noiseIndicatorsNoise3DTNN3DLogTNN3DLogATNN

Case 1-1Washington DC MallMPSNR16.547230.775232.1420
MSSIM0.36330.93870.9449
SAM36.39415.89696.7688
Time (s)310.708207.411

Case 1-2Washington DC MallMPSNR4.382022.990518.1073
MSSIM0.03990.71280.3943
SAM67.767931.642811.0287
Time (s)241.556371.838

Cases 1–3Pavia City CenterMPSNR16.709530.338031.5492
MSSIM0.27420.90370.8933
SAM34.63293.92145.7779
Time (s)117.293883.271

Cases 1–4Pavia City CenterMPSNR4.446920.819110.4980
MSSIM0.02190.36130.0785
SAM68.250219.934752.9761
Time (s)112.223120.206

Cases 1–5Indian PinesMPSNR16.671929.290330.9744
MSSIM0.27050.90460.8793
SAM16.93013.03292.5550
Time (s)226.048163.207

Cases 1–6Indian PinesMPSNR4.435819.73427.7503
MSSIM0.04270.07340.6672
SAM49.19096.653538.7423
Time (s)232.497137.487

Case 2-1Washington DC MallMPSNR8.474733.518342.1350
MSSIM0.08220.94740.9965
SAM48.82567.42831.3193
Time (s)235.137332.130

Case 2-2Washington DC MallMPSNR7.530519.530033.5233
MSSIM0.05760.46010.8926
SAM50.045225.859314.2450
Time (s)376.93292.934

Case 2-3Washington DC MallMPSNR7.458318.892724.8326
MSSIM0.05650.42860.6793
SAM50.136927.041022.7390