# Hyperspectral Image Denoising via Nonconvex Logarithmic Penalty

Authors: Shuo Wang; Zhibin Zhu; Ruwen Zhao; Benxin Zhang
Journal: Mathematical Problems in Engineering
 (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.


---

## Body

## 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 [1–3]. 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 [8–13]; (2) 2D matrix-based low-rank matrix recovery (LRMR) methods [14–21]; and, (3) 3D tensor-based approximation methods [22–33]. 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 [34–36]. 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. Section2 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 [39–41] and mode-k t-SVD proposed in [38].


### 2.1. Notation and Indexing

For a third-order tensorA∈ℝn1×n2×n3, denote Ak,:,:,A:,k,:,andA:,:,k as the kth horizontal, lateral, and frontal slices, respectively. A:,i,j,Ai,:,j,andAi,j,: denote the i,jth mode-1, mode-2, and mode-3 fibers, respectively; Figures 1 and 2 show the horizontal, lateral, and frontal slides, denoted by Ai,:,:,A:,j,:,andA:,:,k, respectively, of a third-order tensor A. Suppose A=aijk∈ℝn1×n2×n3; we adopt the definition of the Frobenius norm of a tensor AF=∑i=1n1∑j=1n2∑k=1n3aijk2 and the definition of the l1 norm of a tensor, A1=∑i=1n1∑j=1n2∑k=1n3aijk.Figure 1
Fibers of a 3rd-order tensor. (a) Mode-1 (column) fiber:x:jk. (b) Mode-2 (row) fiber: xi:k. (c) Mode-3 (tube) fiber:xij:.
(a)(b)(c)Figure 2
Slices of a 3rd-order tensor. (a) Horizontal slices:xi::. (b) Lateral slices: x:j:. (c) Frontal slices: x::k.
(a)(b)(c)LetAi denote A:,:,i, that is, the ith frontal slice of A. Then,(1)unfoldA=A1A2⋮An3∈ℝn1n3×n2.Similarly,Ai∈ℝn1n3×n2, and(2)foldA1A2⋮An3=A∈ℝn1×n2×n3.It can be seen in Figure3, for A∈ℝn×n×3, let Ai=A:,:,i; then, Di are the frontal slices of tensor D, where D is computed by applying the Fast Fourier Transform (FFT) along each tube of A, i.e., Di=fftA,,i and A=ifftDi,,i.Figure 3
Examples of the unfold and fold operators.


### 2.2.t-SVD

In this section, we exploit the proposedt-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 G∈ℝn1×n2×n3 and K∈ℝn2×n4×n3, the t-product D=G∗K is a tensor of size n1×n4×n3. For i=1,2,…,n1 and j=1,2,…,n4,(3)Di,j,:=∑k=1n2Gi,k,:∗Kk,j,:.
Thet-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 D∈ℝn1×n2×n3, the t-SVD of D is given by(4)D=U∗S∗VT,where U∈ℝn1×n1×n3 and V∈ℝn2×n2×n3 are orthogonal tensors, S∈ℝn1×n2×n3 is a rectangular diagonal tensor, and ∗ denotes the t-product.Figure4 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, thet-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.Figure 4
Thet-SVD of a tensor.Algorithm 1:t-SVD.
Input: D∈ℝn1×n2×⋯×nN, w¯=n3n4,…,nN.fori=3 to NdofftD,,i⟶D¯end forfori=1 to w¯doU,S,V=SVDD¯:,:,iU¯:,:,i=U;S¯:,:,i=S;V¯:,:,i=V ;end forfori=3 to NdoifftU¯,,i⟶U;ifftS¯,,i⟶S;ifftV¯,,i⟶V;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 tensorA∈ℝn1×n2×n3, the mode-k block circulation operation is denoted as(5)bcircA,k≔Ak1Aknk⋯Ak2Ak2Ak1⋯Ak3⋮⋮⋱⋮AknkAk^nk−1⋯Ak1,where Aki is the ith mode-k slice of A.The mode-k block diagonalization operation and its inverse operation are defined as(6)bdiagA,k≔Ak1Ak2⋱Aknk,bdfoldbdiagA,k,k≔A.The mode-k block vectorization operation and its inverse operation are defined as(7)bvecA,k≔Ak1Ak2⋮Aknk,bvfoldbvecA,k,k≔A.Definition 2.
(Mode-k t-product). For A∈ℝn1×n2×n3 and ℬ∈ℝn1×n3×n4, the mode-1 t-product is a tensor of size n1×n2×n4:(8)C=A∗1ℬ=bvfoldbcircA,1bvecℬ,1,1⇔C:,j,t=∑s=1n3A:,j,s∗ℬ:,s,t.
ForA∈ℝn1×n2×n3 and ℬ∈ℝn4×n2×n1, the mode-2 t-product is a tensor of size n4×n2×n3:(9)C=A∗2ℬ=bvfoldbcircA,2bvecℬ,2,2⟺Ci,:,t=∑s=1n1As,:,t∗ℬi,;,s.
ForA∈ℝn1×n2×n3 and ℬ∈ℝn2×n4×n3, the mode-3 t-product is a tensor of size n1×n4×n3:(10)C=A∗3ℬ=bvfoldbcircA,3bvecℬ,3,3⟺Ci,j,:=∑s=1n1Ai,s,:∗ℬs,j,:.Definition 3.
(Mode-k identity tensor). ℐk∈ℝn1×n2×n3 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 A∈ℝn1×n2×n3,ATk∈ℝn2×n1×n3 is the mode-k conjugate transpose of A, which is obtained by transposing each of the mode-k slices and then reversing the order of transposed mode-k slices 2 through nk.Definition 5.
(Mode-k diagonal tensor). For A∈ℝn1×n2×n3, each of its mode-k slices is a diagonal matrix, and then, A is a mode-k diagonal.Definition 6.
(Mode-k orthogonal tensor). If(11)ATk∗kA=A∗kATk=ℐk,where the tensor A∈ℝn1×n2×n3 is mode-k orthogonal.Definition 7.
(Tensor mode-k permutation). For A∈ℝn1×n2×n3, A¯k≔permuteA,k is the mode-k permutation of A, whose ith mode-3 slice is the ith mode-k slice of A, and its inverse operation is A≔ipermuteA¯k,k.Theorem 2.
(Mode-k t-SVD). The factorization of A∈ℝn1×n2×n3 is(12)A=Uk∗kSk∗kVkTk,k=1,2,and3,where Uk and Vk are the mode-k orthogonal tensors and Sk∈ℝn1×n2×n3 is the mode-k diagonal tensor.Theorem2 is proven in [38] (Th. 2).Definition 8.
(Tensor fibered rank). ForA∈ℝn1×n2×n3, rankfA is the fibered rank of A, whose kth element is the mode-k tensor fibered rank rankfkA. rankfkA=maxrankmkA, where the ith element of rankmkA is the rank of the ith mode-k slice of A¯k=fftA,,k.
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).Algorithm 2: Mode-k t-SVD for three-way tensors.
Input: A∈ℝn1×n2×n3, kfftA,,k⟶A¯kfori=1 to nkdoU,S,V=SVDA¯kki;U⟶U¯kki;S⟶S¯kki;V⟶V¯kkiend forifftU¯k,,k⟶Uk;ifftS¯k,,k⟶Sk;ifftV¯k,,k⟶Vkoutput: Uk,Sk,Vk


## 2.1. Notation and Indexing

For a third-order tensorA∈ℝn1×n2×n3, denote Ak,:,:,A:,k,:,andA:,:,k as the kth horizontal, lateral, and frontal slices, respectively. A:,i,j,Ai,:,j,andAi,j,: denote the i,jth mode-1, mode-2, and mode-3 fibers, respectively; Figures 1 and 2 show the horizontal, lateral, and frontal slides, denoted by Ai,:,:,A:,j,:,andA:,:,k, respectively, of a third-order tensor A. Suppose A=aijk∈ℝn1×n2×n3; we adopt the definition of the Frobenius norm of a tensor AF=∑i=1n1∑j=1n2∑k=1n3aijk2 and the definition of the l1 norm of a tensor, A1=∑i=1n1∑j=1n2∑k=1n3aijk.Figure 1
Fibers of a 3rd-order tensor. (a) Mode-1 (column) fiber:x:jk. (b) Mode-2 (row) fiber: xi:k. (c) Mode-3 (tube) fiber:xij:.
(a)(b)(c)Figure 2
Slices of a 3rd-order tensor. (a) Horizontal slices:xi::. (b) Lateral slices: x:j:. (c) Frontal slices: x::k.
(a)(b)(c)LetAi denote A:,:,i, that is, the ith frontal slice of A. Then,(1)unfoldA=A1A2⋮An3∈ℝn1n3×n2.Similarly,Ai∈ℝn1n3×n2, and(2)foldA1A2⋮An3=A∈ℝn1×n2×n3.It can be seen in Figure3, for A∈ℝn×n×3, let Ai=A:,:,i; then, Di are the frontal slices of tensor D, where D is computed by applying the Fast Fourier Transform (FFT) along each tube of A, i.e., Di=fftA,,i and A=ifftDi,,i.Figure 3
Examples of the unfold and fold operators.


## 2.2.t-SVD

In this section, we exploit the proposedt-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 G∈ℝn1×n2×n3 and K∈ℝn2×n4×n3, the t-product D=G∗K is a tensor of size n1×n4×n3. For i=1,2,…,n1 and j=1,2,…,n4,(3)Di,j,:=∑k=1n2Gi,k,:∗Kk,j,:.
Thet-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 D∈ℝn1×n2×n3, the t-SVD of D is given by(4)D=U∗S∗VT,where U∈ℝn1×n1×n3 and V∈ℝn2×n2×n3 are orthogonal tensors, S∈ℝn1×n2×n3 is a rectangular diagonal tensor, and ∗ denotes the t-product.Figure4 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, thet-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.Figure 4
Thet-SVD of a tensor.Algorithm 1:t-SVD.
Input: D∈ℝn1×n2×⋯×nN, w¯=n3n4,…,nN.fori=3 to NdofftD,,i⟶D¯end forfori=1 to w¯doU,S,V=SVDD¯:,:,iU¯:,:,i=U;S¯:,:,i=S;V¯:,:,i=V ;end forfori=3 to NdoifftU¯,,i⟶U;ifftS¯,,i⟶S;ifftV¯,,i⟶V;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 tensorA∈ℝn1×n2×n3, the mode-k block circulation operation is denoted as(5)bcircA,k≔Ak1Aknk⋯Ak2Ak2Ak1⋯Ak3⋮⋮⋱⋮AknkAk^nk−1⋯Ak1,where Aki is the ith mode-k slice of A.The mode-k block diagonalization operation and its inverse operation are defined as(6)bdiagA,k≔Ak1Ak2⋱Aknk,bdfoldbdiagA,k,k≔A.The mode-k block vectorization operation and its inverse operation are defined as(7)bvecA,k≔Ak1Ak2⋮Aknk,bvfoldbvecA,k,k≔A.Definition 2.
(Mode-k t-product). For A∈ℝn1×n2×n3 and ℬ∈ℝn1×n3×n4, the mode-1 t-product is a tensor of size n1×n2×n4:(8)C=A∗1ℬ=bvfoldbcircA,1bvecℬ,1,1⇔C:,j,t=∑s=1n3A:,j,s∗ℬ:,s,t.
ForA∈ℝn1×n2×n3 and ℬ∈ℝn4×n2×n1, the mode-2 t-product is a tensor of size n4×n2×n3:(9)C=A∗2ℬ=bvfoldbcircA,2bvecℬ,2,2⟺Ci,:,t=∑s=1n1As,:,t∗ℬi,;,s.
ForA∈ℝn1×n2×n3 and ℬ∈ℝn2×n4×n3, the mode-3 t-product is a tensor of size n1×n4×n3:(10)C=A∗3ℬ=bvfoldbcircA,3bvecℬ,3,3⟺Ci,j,:=∑s=1n1Ai,s,:∗ℬs,j,:.Definition 3.
(Mode-k identity tensor). ℐk∈ℝn1×n2×n3 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 A∈ℝn1×n2×n3,ATk∈ℝn2×n1×n3 is the mode-k conjugate transpose of A, which is obtained by transposing each of the mode-k slices and then reversing the order of transposed mode-k slices 2 through nk.Definition 5.
(Mode-k diagonal tensor). For A∈ℝn1×n2×n3, each of its mode-k slices is a diagonal matrix, and then, A is a mode-k diagonal.Definition 6.
(Mode-k orthogonal tensor). If(11)ATk∗kA=A∗kATk=ℐk,where the tensor A∈ℝn1×n2×n3 is mode-k orthogonal.Definition 7.
(Tensor mode-k permutation). For A∈ℝn1×n2×n3, A¯k≔permuteA,k is the mode-k permutation of A, whose ith mode-3 slice is the ith mode-k slice of A, and its inverse operation is A≔ipermuteA¯k,k.Theorem 2.
(Mode-k t-SVD). The factorization of A∈ℝn1×n2×n3 is(12)A=Uk∗kSk∗kVkTk,k=1,2,and3,where Uk and Vk are the mode-k orthogonal tensors and Sk∈ℝn1×n2×n3 is the mode-k diagonal tensor.Theorem2 is proven in [38] (Th. 2).Definition 8.
(Tensor fibered rank). ForA∈ℝn1×n2×n3, rankfA is the fibered rank of A, whose kth element is the mode-k tensor fibered rank rankfkA. rankfkA=maxrankmkA, where the ith element of rankmkA is the rank of the ith mode-k slice of A¯k=fftA,,k.
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).Algorithm 2: Mode-k t-SVD for three-way tensors.
Input: A∈ℝn1×n2×n3, kfftA,,k⟶A¯kfori=1 to nkdoU,S,V=SVDA¯kki;U⟶U¯kki;S⟶S¯kki;V⟶V¯kkiend forifftU¯k,,k⟶Uk;ifftS¯k,,k⟶Sk;ifftV¯k,,k⟶Vkoutput: Uk,Sk,Vk


## 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 by(13)ϕx;a=1alog1+ax,where a>0 controls the degree of nonconvexity of the penalty function. This function satisfies the following conditions:(A1)
ϕx;a∈ℂ2ℝ+∗,ϕx;a∈C0ℝ+,ℝ+≔t∈ℝ:t≥0,andℝ+∗≔t∈ℝ:t>0(A2)
ϕ′x;a>0and∀t∈ℝ+∗(A3)
ϕ″x;a≤0and∀t∈ℝ+∗(A4)
supt∈R+∗ϕ′x;a=1andinft∈R+∗ϕ″x;a=ϕ″0+;a=−aThe proximity operatorΘ:ℝ⟶ℝ associated with the nonconvex function ϕx;a is(14)Θy;λ¯,a=proxϕy;λ¯,a=argminx∈ℝfx=12y−x2+λ¯ϕx;a.In [44–46], the authors prove that, for 0<a≤1/λ¯, the function f 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,(15)proxλ¯ϕy;λ¯,a=0,∀y<λ¯,and is given by(16)proxλ¯ϕy;λ¯,a=y2−12a+y2+12a2−λ¯a,y>λ¯,0,y≤λ¯.


### 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 as(17)minA.N.SrankfA+λ1NF2+λ2S1s.t.Y=A+N+S,where Y is the observed noisy HSI, N is Gaussian noise, S is sparse noise, A∈ℝn1×n2×n3 is the clean hyperspectral signal, which has the low-fibered-rank tensor property, and λ1 and λ2 are regularization parameters.The Gaussian noise model (large degree of freedom) corresponds to a dense noise type [47, 48]. N should not be low in rank. Otherwise, A 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 N is normally full and much larger than the rank of A. Therefore, the low-fibered-rank component A can be recovered from equation (17) by properly choosing the tuning parameters λ1 and λ2.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:(18)minA.N.S∑k=13τkℱkA+λ1NF2+λ2S1s.t.Y=A+N+S,where τk≥0k=1,2,and3 and ∑k=13τk=1 is the weight of the fibered rank. ℱkA are set as LogATNNkA in the HIS denoising model:(19)LogATNNkA=∑i=1nk∑j=1m1alog1+aσjA¯kki,nk=n1,n2,andn3,where m is the number of singular values of A¯, σjA¯kki is the jth singular values of A¯kki, and A¯kki is the ith mode-k slice of A¯=fftA,,k.


### 3.3. ADMM for Solving Model (18)

We use the ADMM to solve equation (18). By casting the three auxiliary variables Xkk=1,2,and3, equation (18) can be rewritten as follows:(20)minA.N⋅S∑k=13τkℱkXk+λ1NF2+λ2S1s.t.Y−A+N+S=0A−Xk=0k=1,2,and3.The augmented Lagrangian function of equation (18) is(21)Lζk,ρXk,A,N,S,W,Q=∑k=13τkℱkXk+A−Xk,Wk+ζk2A−XkF2+λ1NF2+λ2S1+Y−A+N+S,Q+ρ2Y−A+N+SF2,where Wk and Q are the Lagrange multipliers and ζkk=1,2,3 and ρ are positive scalars. Now we can solve the problem within the ADMM framework.With the other parameters fixed,X can be updated by solving(22)minXLζk,ρXk,A,N,S,W,Q=minX∑k=13τkℱkXk+A−Xk,Wk+ζk2A−XkF2,which is equivalent to the following subproblem for Xk∈ℝn1×n2×n3k=1,2,and3:(23)Xkl+1=argminXkτkℱkXk+ζk2Al−Xk+WklζkF2.To solve equation (23), we can rewrite it as(24)argminXλ¯logATNNkX+12X−ℳF2,where ℳ=A+Wk/ξk=U∗kS∗kVTk,S¯k=fftS,,k, and λ¯=τk/ζk. From equation (16),(25)S^k=θS¯ki,j,r=S¯ki,j,r2−12a+S¯ki,j,r2+12a2−λ¯asignS¯ki,j,r,S¯ki,j,r≥λ¯,0,S¯ki,j,r<λ¯,where S˜k=ifftS^k,,kandX=Dℳ,k=U∗kS˜∗kV⊤k,λ¯>0and0<a≤1/λ¯. Therefore,(26)Xkl+1=DAl+Wklζk,k.Though model (20) is based on a nonconvex penalty function, the parameters of which are constrained to ensure the convexity of the subproblem Xk.With other parameters fixed,A can be updated by solving(27)argminyLζk,ρXk,A,N,S,W,Q=∑i=13A−Xk,Wk+ζk2A−XkF2+Y−A+N+S,Q+ρ2Y−A+N+SF2,which is equivalent to the following subproblem:(28)Al+1=argminy∑k=13ζk2A−Xkl+1+WklζkF2+ρ2Y−A+N+S+QlρF2.It has the following closed-form solution:(29)Al+1=∑k=13ζkXkl+1−Wkl/ζk+ρY−Nl−Sl+Ql/ρ∑k=13ζk+ρ.With the other parameters fixed,N can be updated by solving(30)argminNLζk,ρXk,A,N,S,W,Q=λ1NF2+Y−A+N+S,Q+ρ2Y−A+N+SF2,which is equivalent to the following subproblem:(31)Nl+1=argminNλ1NF2+ρ2Y−Al+1+N+Sl+QlρF2,and it has the following closed-form solution:(32)Nl+1=ρY−Al+1−Sl+Ql/ρ2λ1+ρ.With the other parameters fixed,S can be updated by solving(33)argminsLζk,ρXk,A,N,S,ℳ,Q=λ2S1+Y−A+N+S,Q+ρ2Y−A+N+SF2,which is equivalent to the following subproblem:(34)Sl+1=argminSλ2S1+ρ2Y−Al+1+Nl+1+S+QlρF2.Using the tensor soft thresholding operator, the following solution can be obtained [38]:(35)Sl+1=signY−Al+1−Nl+1+QlρmaxY−Al+1−Nl+1+Qlρ−λ2ρ,0.The Lagrange multipliersWk and Q can be updated by solving(36)Wkl+1=Wkl+ζkAl+1−Xkl+1,(37)Ql+1=Ql+ρY−Al+1−Nl+1−Sl+1.Hence, the proposed algorithm for HSI denoising is summarized in Algorithm3.Algorithm 3: Algorithm for HSI denoising with the ADMM framework.
Input: The noisy HSI Y, parameters τk,ζkk=1,2,and3,λ1,λ2,ρ,t=1.2,l=0.X0=0,N0=0,S0=0,Xk0=0,Wk0=0, and Q0=0.while not converged doUpdateXkl+1 with equation (26), k=1,2,and3.UpdateAl+1 with equation (29).UpdateNl+1 with equation (32).UpdateSl+1 with equation (35).UpdateWkl+1 with equation (36), k=1,2,and3.UpdateQl+1 with equation (37).Letζk=tζk,k=1,2,and3;ρ=tρ; and l=l+1.Check the convergence condition:Al+1−AlF/AlF<10−4.end whileOutput: Denoised HIS X.


## 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 by(13)ϕx;a=1alog1+ax,where a>0 controls the degree of nonconvexity of the penalty function. This function satisfies the following conditions:(A1)
ϕx;a∈ℂ2ℝ+∗,ϕx;a∈C0ℝ+,ℝ+≔t∈ℝ:t≥0,andℝ+∗≔t∈ℝ:t>0(A2)
ϕ′x;a>0and∀t∈ℝ+∗(A3)
ϕ″x;a≤0and∀t∈ℝ+∗(A4)
supt∈R+∗ϕ′x;a=1andinft∈R+∗ϕ″x;a=ϕ″0+;a=−aThe proximity operatorΘ:ℝ⟶ℝ associated with the nonconvex function ϕx;a is(14)Θy;λ¯,a=proxϕy;λ¯,a=argminx∈ℝfx=12y−x2+λ¯ϕx;a.In [44–46], the authors prove that, for 0<a≤1/λ¯, the function f 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,(15)proxλ¯ϕy;λ¯,a=0,∀y<λ¯,and is given by(16)proxλ¯ϕy;λ¯,a=y2−12a+y2+12a2−λ¯a,y>λ¯,0,y≤λ¯.


## 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 as(17)minA.N.SrankfA+λ1NF2+λ2S1s.t.Y=A+N+S,where Y is the observed noisy HSI, N is Gaussian noise, S is sparse noise, A∈ℝn1×n2×n3 is the clean hyperspectral signal, which has the low-fibered-rank tensor property, and λ1 and λ2 are regularization parameters.The Gaussian noise model (large degree of freedom) corresponds to a dense noise type [47, 48]. N should not be low in rank. Otherwise, A 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 N is normally full and much larger than the rank of A. Therefore, the low-fibered-rank component A can be recovered from equation (17) by properly choosing the tuning parameters λ1 and λ2.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:(18)minA.N.S∑k=13τkℱkA+λ1NF2+λ2S1s.t.Y=A+N+S,where τk≥0k=1,2,and3 and ∑k=13τk=1 is the weight of the fibered rank. ℱkA are set as LogATNNkA in the HIS denoising model:(19)LogATNNkA=∑i=1nk∑j=1m1alog1+aσjA¯kki,nk=n1,n2,andn3,where m is the number of singular values of A¯, σjA¯kki is the jth singular values of A¯kki, and A¯kki is the ith mode-k slice of A¯=fftA,,k.


## 3.3. ADMM for Solving Model (18)

We use the ADMM to solve equation (18). By casting the three auxiliary variables Xkk=1,2,and3, equation (18) can be rewritten as follows:(20)minA.N⋅S∑k=13τkℱkXk+λ1NF2+λ2S1s.t.Y−A+N+S=0A−Xk=0k=1,2,and3.The augmented Lagrangian function of equation (18) is(21)Lζk,ρXk,A,N,S,W,Q=∑k=13τkℱkXk+A−Xk,Wk+ζk2A−XkF2+λ1NF2+λ2S1+Y−A+N+S,Q+ρ2Y−A+N+SF2,where Wk and Q are the Lagrange multipliers and ζkk=1,2,3 and ρ are positive scalars. Now we can solve the problem within the ADMM framework.With the other parameters fixed,X can be updated by solving(22)minXLζk,ρXk,A,N,S,W,Q=minX∑k=13τkℱkXk+A−Xk,Wk+ζk2A−XkF2,which is equivalent to the following subproblem for Xk∈ℝn1×n2×n3k=1,2,and3:(23)Xkl+1=argminXkτkℱkXk+ζk2Al−Xk+WklζkF2.To solve equation (23), we can rewrite it as(24)argminXλ¯logATNNkX+12X−ℳF2,where ℳ=A+Wk/ξk=U∗kS∗kVTk,S¯k=fftS,,k, and λ¯=τk/ζk. From equation (16),(25)S^k=θS¯ki,j,r=S¯ki,j,r2−12a+S¯ki,j,r2+12a2−λ¯asignS¯ki,j,r,S¯ki,j,r≥λ¯,0,S¯ki,j,r<λ¯,where S˜k=ifftS^k,,kandX=Dℳ,k=U∗kS˜∗kV⊤k,λ¯>0and0<a≤1/λ¯. Therefore,(26)Xkl+1=DAl+Wklζk,k.Though model (20) is based on a nonconvex penalty function, the parameters of which are constrained to ensure the convexity of the subproblem Xk.With other parameters fixed,A can be updated by solving(27)argminyLζk,ρXk,A,N,S,W,Q=∑i=13A−Xk,Wk+ζk2A−XkF2+Y−A+N+S,Q+ρ2Y−A+N+SF2,which is equivalent to the following subproblem:(28)Al+1=argminy∑k=13ζk2A−Xkl+1+WklζkF2+ρ2Y−A+N+S+QlρF2.It has the following closed-form solution:(29)Al+1=∑k=13ζkXkl+1−Wkl/ζk+ρY−Nl−Sl+Ql/ρ∑k=13ζk+ρ.With the other parameters fixed,N can be updated by solving(30)argminNLζk,ρXk,A,N,S,W,Q=λ1NF2+Y−A+N+S,Q+ρ2Y−A+N+SF2,which is equivalent to the following subproblem:(31)Nl+1=argminNλ1NF2+ρ2Y−Al+1+N+Sl+QlρF2,and it has the following closed-form solution:(32)Nl+1=ρY−Al+1−Sl+Ql/ρ2λ1+ρ.With the other parameters fixed,S can be updated by solving(33)argminsLζk,ρXk,A,N,S,ℳ,Q=λ2S1+Y−A+N+S,Q+ρ2Y−A+N+SF2,which is equivalent to the following subproblem:(34)Sl+1=argminSλ2S1+ρ2Y−Al+1+Nl+1+S+QlρF2.Using the tensor soft thresholding operator, the following solution can be obtained [38]:(35)Sl+1=signY−Al+1−Nl+1+QlρmaxY−Al+1−Nl+1+Qlρ−λ2ρ,0.The Lagrange multipliersWk and Q can be updated by solving(36)Wkl+1=Wkl+ζkAl+1−Xkl+1,(37)Ql+1=Ql+ρY−Al+1−Nl+1−Sl+1.Hence, the proposed algorithm for HSI denoising is summarized in Algorithm3.Algorithm 3: Algorithm for HSI denoising with the ADMM framework.
Input: The noisy HSI Y, parameters τk,ζkk=1,2,and3,λ1,λ2,ρ,t=1.2,l=0.X0=0,N0=0,S0=0,Xk0=0,Wk0=0, and Q0=0.while not converged doUpdateXkl+1 with equation (26), k=1,2,and3.UpdateAl+1 with equation (29).UpdateNl+1 with equation (32).UpdateSl+1 with equation (35).UpdateWkl+1 with equation (36), k=1,2,and3.UpdateQl+1 with equation (37).Letζk=tζk,k=1,2,and3;ρ=tρ; and l=l+1.Check the convergence condition:Al+1−AlF/AlF<10−4.end whileOutput: Denoised HIS X.


## 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 is256×256). The Pavia City Center data uses only a subimage (80 bands and size of each band is 200×200). 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 145×145×224. 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:(38)PSNR=10∗log10L2MN∑x=1M∑y=1NI^x,y−Ix,y2,MPSNR=1B∑x=1BPSNRi,where M×N represents the image size in the space, I^x,y represents the original image, Ix,y represents the distortion image, L is a pixel that can achieve the maximum value, and B is the number of PSNR:(39)SSIM=2μxμy+C12σxy+C2μx2+μy2+C1σx2+σy2+C2,MSSIMX,Y=1Q∑j=1QSSIMxj,yj,where C1 is a constant for μx2+μy2,C2 is the same as C1,σx and σy represent the x and y standard deviations, respectively, XandY represent the original image and the distorted image, respectively, xj and yj represent the jth local window contents, and Q is the number of local windows:(40)SAMv,v′=cos−1∑i=1Zvivi′∑i=1Zvi2∑i=1Zvi′2,where Z is the number of wavelengths, v and v′ represent the spectrum vectors, and(41)v=v1,v2,…,vn,v′=v1′,v2′,…,vn′.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 anN-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τk, the regularization parameters λ1 and λ2, the threshold parameter ςk=τk/ξ, the penalty parameter ρ=1/meanξ, and a constant a in 3DLogATNN. In all the following experiments, the parameters in these compared methods were manually adjusted according to their default strategies.The regularization parameterλ1 for 3DLogATNN: it is easy to see that λ1 is the parameter used to restrict the sparsity of the Gaussian noise. We set λ1=C/β, where β is the standard deviation of Gaussian noise and C is a tuning parameter. The results were based on the simulated data experiment in case 1-1. Figure 5 shows the restoration results as C 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 C changing from 0.002 to 0.003. Therefore, we suggest the use of C=0.002 in all the simulated data experiments.Figure 5
Sensitivity analysis of parameterλ1 for 3DLogATNN. (a) Change in the MPSNR value. (b) Change in the MSSIM value.
(a)(b)The regularization parameterλ2 for 3DLogATNN: it is easy to see that λ2 is the parameter used to restrict the sparsity of the impulse noise. We set λ2=Bλ, where(42)λ=τ1maxn2,n3n1+τ2maxn3,n1n2+τ3maxn1,n2n3,and B is a tuning parameter. The results were based on the simulated data experiment in case 2-1. Figure 6 shows the restoration results as B 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 B changing from 0.4 to 0.5. Therefore, we suggest the use of B=0.5 in all the simulated data experiments.Figure 6
Sensitivity analysis of parameterλ2 for 3DLogATNN. (a) Change in the MPSNR value. (b) Change in the MSSIM value.
(a)(b)The constanta for 3DLogATNN: the parameter a>0 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 a 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 a changing from 0.05 to 0.06. Therefore, we suggest the use of a=0.05 in all the simulated data experiments.Figure 7
Sensitivity analysis of parametera for 3DLogATNN. (a) Change in the MPSNR value. (b) Change in the MSSIM value.
(a)(b)We adjust the parameters to achieve the best visual result, and the parameter setting is presented in Table1.Table 1
Parameter settings in the HSI denoising methods for simulated data.
MethodDataτkξλ1λ2εa3DTNNWashington DC Mall1,1,0.001/2.001102,102,1020.002/βλ——Pavia City CenterIndian Pines3DLogTNNWashington DC Mall1,1,0.001/2.001105,105,1050.005/β0.011λ50—Pavia City Center80Indian Pines503DLogATNNWashington DC Mall1,1,0.001/2.001103,103,1030.002/β0.5λ—0.05Pavia City CenterIndian PinesCompared 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 8–10. 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.Table 2
Performance comparison of competing methods. The best result in each case is highlighted in bold.
CaseDataGaussian noiseImpulse noiseIndicatorsNoise3DTNN3DLogTNN3DLogATNNCase 1-1Washington DC Mallβ∈U0.1,0.2—MPSNR16.547230.775232.142032.5883MSSIM0.36330.93870.94490.9508SAM36.39415.53805.89696.7688Time (s)119.688310.708207.411Case 1-2Washington DC Mallβ∈U0.55,0.65—MPSNR4.382022.990518.107324.9584MSSIM0.03990.71280.39430.7436SAM67.76799.891231.642811.0287Time (s)241.556371.838202.24Cases 1–3Pavia City Centerβ∈U0.1,0.2—MPSNR16.709530.338031.549232.9023MSSIM0.27420.90370.89330.9288SAM34.63293.92145.77793.6840Time (s)41.121117.293883.271Cases 1–4Pavia City Centerβ∈U0.55,0.65—MPSNR4.446920.819110.498024.9550MSSIM0.02190.36130.07850.6598SAM68.250219.934752.97616.1381Time (s)112.223120.20671.473Cases 1–5Indian Pinesβ∈U0.1,0.2—MPSNR16.671929.290330.974431.7479MSSIM0.27050.90460.87930.8438SAM16.93013.03292.55502.4113Time (s)90.68226.048163.207Cases 1–6Indian Pinesβ∈U0.55,0.65—MPSNR4.435819.73427.750323.2606MSSIM0.04270.73540.07340.6672SAM49.19096.653538.74235.9968Time (s)120.883232.497137.487Case 2-1Washington DC Mall—δ∈U0.35,0.45MPSNR8.474733.518342.135049.3435MSSIM0.08220.94740.99650.9992SAM48.82567.42831.31930.7387Time (s)235.137332.130227.661Case 2-2Washington DC Mall—δ∈U0.45,0.55MPSNR7.530519.530033.523345.1362MSSIM0.05760.46010.89260.9981SAM50.045225.859314.24501.0228Time (s)263.654376.93292.934Case 2-3Washington DC Mall—δ∈U0.55,0.65MPSNR7.458318.892724.832642.2501MSSIM0.05650.42860.67930.9964SAM50.136927.041022.73901.3989Time (s)239.041337.120224.460Cases 2–4Pavia City Center—δ∈U0.35,0.45MPSNR8.854241.260045.921353.6297MSSIM0.04980.99560.99800.9995SAM45.63000.36930.32040.1519Time (s)91.438140.218112.862Cases 2–5Pavia City Center—δ∈U0.55,0.65MPSNR7.064918.43637.292139.0625MSSIM0.02320.28960.02510.9888SAM46.839726.478546.65341.1001Time (s)111.60676.873112.652Cases 2–6Indian Pines—δ∈U0.35,0.45MPSNR9.277839.063737.347338.3288MSSIM0.07740.97300.97630.9784SAM34.41061.03030.84210.7242Time (s)149.143213.263150.303Cases 2–7Indian Pines—δ∈U0.55,0.65MPSNR7.488114.118913.506336.9272MSSIM0.04020.16360.15420.9704SAM40.390221.608622.92390.8861Time (s)150.468229.334165.124Cases 3–1Washington DC Mallβ=0.3δ=0.1MPSNR9.000025.592627.843527.3793MSSIM0.10870.81290.86240.8470SAM53.89058.58457.84968.5841Time (s)137.651317.362198.887Case 3-2Washington DC Mallβ=0.1δ=0.4MPSNR8.159124.527330.400030.5960MSSIM0.07510.72260.92130.9226SAM49.859312.86866.38626.4664Time (s)202.923332.368162.639Case 3-3Washington DC Mallβ=0.3δ=0.5MPSNR5.700916.337014.905520.9115MSSIM0.03590.29830.24060.6499SAM55.956721.925528.763513.2614Time (s)220.647367.041200.073Case 3-4Pavia City Centerβ=0.3δ=0.1MPSNR9.097425.799713.131926.8088MSSIM0.06420.74480.13270.8029SAM53.28585.400942.71535.6316Time (s)51.95593.24965.679Cases 3–5Pavia City Centerβ=0.1δ=0.4MPSNR8.471128.496330.214631.0972MSSIM0.04430.86150.86130.9098SAM46.80834.26066.46143.9908Time (s)63.319116.12763.407Cases 3–6Indian Pinesβ=0.1δ=0.4MPSNR8.908325.914830.128430.7026MSSIM0.07190.68170.86440.8164SAM35.51255.20592.81852.7058Time (s)134.676234.822165.194Cases 3–7Indian Pinesβ=0.3δ=0.1MPSNR9.217623.622326.862627.1730MSSIM0.09710.81050.78070.7415SAM34.23245.55984.07143.9805Time (s)102.124224.456150.277Cases 4–1Washington DC Mallβ∈U0.3,0.4δ=0.2MPSNR7.153322.789324.898425.2959MSSIM0.07020.71240.80040.8008SAM57.150611.23789.47889.3893Time (s)220.148370.379226.554Cases 4–2Washington DC Mallβ∈U0.2,0.3δ=0.3MPSNR7.682922.130125.147725.7826MSSIM0.07370.68540.81460.8178SAM53.481912.56299.58059.3802Time (s)200.652338.836165.677Case 4-3Washington DC Mallβ∈U0.4,0.5δ=0.1MPSNR6.237723.72769.608025.6329MSSIM0.06020.72830.11280.8029SAM61.581410.016353.28279.1467Time (s)190.534323.047172.353Case 4-4Pavia City Centerβ∈U0.2,0.3δ=0.3MPSNR7.860223.587324.627325.6363MSSIM0.04210.59110.69430.7801SAM51.64279.91426.06695.9457Time (s)113.882150.16091.212Case 4-5Pavia City Centerβ∈U0.3,0.4δ=0.2MPSNR7.253523.519324.336825.1019MSSIM0.03950.57040.66720.7372SAM56.230110.13076.01745.9217Time (s)92.428139.00687.948Cases 4–6Indian Pinesβ∈U0.2,0.3δ=0.3MPSNR8.197022.191526.003626.5240MSSIM0.06980.78020.77710.7409SAM37.69257.05874.58084.3126Time (s)117.283230.989152.424Cases 4–7Indian Pinesβ∈U0.4,0.5δ=0.1MPSNR6.357220.838015.048623.8612MSSIM0.05780.75110.20080.6976SAM43.24857.033319.41835.5472Time (s)106.580226.007133.670Cases 5–1Washington DC Mallβ=0.1δ∈U0.5,0.6MPSNR6.876616.589310.670528.4585MSSIM0.04470.28550.11300.8819SAM51.257927.348843.52167.8981Time (s)213.342317.644172.836Cases 5–2Washington DC Mallβ=0.3δ∈U0.3,0.4MPSNR6.698719.960512.141323.0269MSSIM0.05570.55400.15810.7539SAM55.593815.500241.954911.3455Time (s)200.605323.034166.227Cases 5–3Pavia City Centerβ=0.1δ∈U0.5,0.6MPSNR7.214320.501226.637429.0107MSSIM0.02620.36580.79000.8651SAM47.569119.06505.42844.5068Time (s)105.725137.456103.408Case 5-4Pavia City Centerβ=0.3δ∈U0.3,0.4MPSNR6.902621.095422.976723.7700MSSIM0.03170.42850.61310.6875SAM53.563612.63936.36476.2442Time (s)88.2136.35586.321Case 5-5Indian Pinesβ=0.1δ∈U0.5,0.6MPSNR7.621417.678727.658528.9263MSSIM0.04540.25270.74420.7765SAM39.873614.34204.16423.2634Time (s)143.077243.403159.432Case 5-6Indian Pinesβ=0.2δ∈U0.4,0.5MPSNR7.617720.638526.286326.3867MSSIM0.05280.51390.77330.7446SAM39.74189.11724.49094.3622Time (s)106.580226.007133.670Cases 6–1Washington DC Mallβ∈U0.2,0.3δ∈U0.2,0.3MPSNR8.112723.302226.335126.2174MSSIM0.08450.74640.83780.8400SAM53.204011.01238.81038.8638Time (s)196.620340.366167.766Cases 6–2Washington DC Mallβ∈U0.1,0.2δ∈U0.3,0.4MPSNR8.241823.923127.940228.4975MSSIM0.08060.72560.88020.8831SAM50.583112.21937.82577.9288Time (s)199.073341.498167.085Cases 6–3Washington DC Mallβ∈U0.4,0.5δ∈U0.1,0.2MPSNR5.957622.78229.129424.5838MSSIM0.05440.68960.09880.7766SAM61.417011.302553.70849.8922Time (s)212.886326.244173.614Cases 6–4Pavia City Centerβ∈U0.1,0.2δ∈U0.4,0.5MPSNR7.684923.140425.833427.7378MSSIM0.03390.53150.75190.8332SAM58.469512.97265.99275.2728Time (s)90.948136.67789.459Case 6-5Pavia City Centerβ∈U0.4,0.5δ∈U0.1,0.2MPSNR5.949822.301324.013924.3767MSSIM0.03010.45130.62790.6621SAM61.324713.43686.25486.2623Time (s)94.386139.48285.760Case 6-6Indian Pinesβ∈U0.1,0.2δ∈U0.4,0.5MPSNR8.059121.639327.984628.3006MSSIM0.05730.48070.81670.7662SAM38.32178.56983.65723.5710Time (s)139.625232.800155.648Case 6-7Indian Pinesβ∈U0.4,0.5δ∈U0.1,0.2MPSNR6.043620.293213.170123.4520MSSIM0.05260.73670.15510.6912SAM44.35147.740123.63095.9105Time (s)113.345236.529137.036Figure 8
Denoising results of the Washington DC Mall dataset with mixture noise.
(a)(b)Figure 9
Denoising results of the Pavia City Center dataset with mixture noise.
(a)(b)Figure 10
Denoising results of the Indian Pines dataset with mixture noise.
(a)(b)


## 5. Conclusion

In this paper, we present a new 3DLogATNN method for HSI denoising by mode-k t-SVD. The logarithmic penalty function is introduced in 3DLogATNN, which enables it to extract low-rank and sparse components more accurately from a degraded HSI contaminated by several types of noise. In addition, the ADMM-based algorithm is applied to effectively solve the proposed HSI denoising model. The experiments have substantiated the superiority of the proposed method over state-of-the-art methods.



---

*Source: 5535169-2021-06-24.xml*