Table Of ContentJOURNALOFLATEXCLASSFILES,VOL.14,NO.8,AUGUST2015 1
Denoising Hyperspectral Image with Non-i.i.d.
Noise Structure
Yang Chen, Xiangyong Cao, Qian Zhao, Deyu Meng, Member, IEEE, Zongben Xu
Abstract—Hyperspectral image (HSI) denoising has been at- The simplest denoising way is to utilize the traditional 2-D
tracting much research attention in remote sensing area due or 1-D denoising methods to reduce noise in the HSI pixel
to its importance in improving the HSI qualities. The existing
by pixel [17] or band by band [12], [14], [30], [44], [50].
HSI denoising methods mainly focus on specific spectral and
However, this kind of processing method ignores the correla-
7 spatialpriorknowledgeinHSIs,andshareacommonunderlying
1 assumption that the embedded noise in HSI is independent and tions between different spectral bands or adjacent pixels, and
0 identicallydistributed(i.i.d.).Inrealscenarios,however,thenoise often results in a relatively low-quality result. Recently, more
2 existed in a natural HSI is always with much more compli- efficient methods have been proposed to handle this issue.
cated non-i.i.d. statistical structures and the under-estimation
b the main idea is to elaborately encode the prior knowledge
to this noise complexity often tends to evidently degenerate the
e on the structure underlying a natural HSI, especially the
robustness of current methods. To alleviate this issue, this paper
F
attempts the first effort to model the HSI noise using a non- characteristic across the spatial and spectral dimensionality.
1 i.i.d. mixture of Gaussians (NMoG) noise assumption, which is E.g., Othman and Qian [38] made an initial attempt to this
finelyinaccordancewiththenoisecharacteristicspossessedbya issue by designing a hybrid spatial-spectral derivative-domain
] naturalHSIandthusiscapableofadaptingvariousnoiseshapes
wavelet shrinkage model based on the dissimilarity of the
V
encountered in real applications. Then we integrate such noise
signal regularity existing along the space and spectrum of
C modelingstrategyintothelow-rankmatrixfactorization(LRMF)
model and propose a NMoG-LRMF model in the Bayesian a natural HSI. And then Chen et al. [5] proposed another
.
s framework. A variational Bayes algorithm is designed to infer approach to encoding both spatial and spectral knowledge by
c the posterior of the proposed model. All involved parameters combining the bivariate wavelet thresholding with principal
[
can be recursively updated in closed-form. Compared with the component analysis. To further enhance the denoising capa-
1 current techniques, the proposed method performs more robust bility, Yuan et al. [51] employed a spectral-spatial adaptive
v beyondthestate-of-the-arts,assubstantiatedbyourexperiments
8 implemented on synthetic and real noisy HSIs. total variation. Later, Chen et al. [9] proposed a spatially
adaptive weighted prior by combining the smoothness and
9 Index Terms—Hyperspectral image denoising, Non-i.i.d. Noise
0 modeling, Low-rank matrix factorization. discontinuity preserving properties along the spectral domain.
0 By further considering the spatial and spectral dependencies,
0 Zhong and Wang [58] proposed a multiple-spectral-band CRF
. I. INTRODUCTION
2 (MSB-CRF) model in a unified probabilistic framework.
HYPERSPECTRALimage(HSI)iscapturedfromsensors
0 Besides, by explicitly treating HSI data as a tensor, a
7 over various bands and contains abundant spatial and
series of methods that expanding wavelet-based method from
1 spectral knowledge across all these bands beyond the tradi-
2-D to 3-D has been proposed. E.g., Dabov [11] designed
: tional gray-scale or RGB images. Due to its preservation of
v VBM3D method by applying the concepts of grouping and
i full-bands information under a real scene, it has been widely
X collaborativefilteringtovideodenoising.Then,Letexieretal.
used in military and civilian aspects such as terrain detection,
[28] proposed a generalized multi-dimensional Wiener filter
r mineral exploration, pharmaceutical counterfeiting, vegetation
a for denoising hyperspectral image. Similarly, Chen et al. [6]
and environmental monitoring [10], [16], [40], [42], [52],
extended Sendur and Selesnick’s bivariate wavelet threshold-
[56]. In real cases, however, a HSI is always corrupted by
ing from 2-D image denoising to 3-D data cube denoising.
noisesduetoequipmentlimitations,suchassensorsensitivity,
For better denoising results, Chen et al. [7] proposed a new
photon effects and calibration error [45]. Besides, due to the
signal denoising algorithm by using neighbouring wavelet
limited radiance energy and generally narrow band width, the
coefficients, which considered both translation-invariant (TI)
energy captured by each sensor might be low and thus the
andnon-TIversions.Later,asanextensionofBM3Dmethod,
shot noise and thermal noise then always happen inevitably.
Maggioni et al. [32] presented BM4D model. Meanwhile,
These noise severely degrades the quality of the imagery
another type of method that based on tensor decomposition
and limits the performance of the subsequent processing, i.e.,
has appeared. Karami et al. [26] developed a Genetic Kernel
classification[43],unmixing[2],[41]andtargetdetection[46],
Tucker Decomposition (GKTD) algorithm to exploit both
on data. Therefore, it is a critical preprocessing step to reduce
spectral and the spatial information in HSIs. To address the
the HSI noise [20], [39] to a general HSI processing task.
uniquenessofmultipleranksofTuckerDecomposition,Liuet
Yang Chen, Xiangyong Cao, Qian Zhao, Deyu Meng (corresponding al. [31] proposed PARAFAC method that make the number
author), Zongben Xu are with the School of Mathematics and Statistics of estimated rank reduced to one. Later, Peng et al. [39]
and Ministry of Education Key Lab of Intelligent Networks and Net-
proposed a decomposable nonlocal tensor dictionary learning
work Security, Xi’an Jiaotong University, Xi’an 710049, China (e-mail:
[email protected]). (TDL) model, which fully considered the non-local similarity
JOURNALOFLATEXCLASSFILES,VOL.14,NO.8,AUGUST2015 2
overspaceandtheglobalcorrelationacrossspectrum.Among
these methods, BM4D and TDL achieved the state-of-the-art
performance in more general noisy MSI cases. （a）
Most of current HSI denoising methods have mainly con-
sidered the HSI prior spectral and spatial knowledge into
= = = = =
their methods, while only use L loss term to rectify the
2
deviation between the reconstruction and the original HSI. In
（b）
the viewpoint of statistical theory, the utilization of such loss
termimpliesthattheHSIdatanoisefollowsani.i.d.Gaussian
distribution. However, in real scenarios, HSI noises generally + + + + +
have more complicated statistical structures. This means that
such easy loss term is too subjective to reflect the real cases （c）
and inclines to degenerate the performance of the methods in
more complex while more realistic non-i.i.d. noise case .
The idea of considering more complex HSI noise beyond 6000 10000 3000 3000 4000
only Gaussian in HSI denoising has been attracting attention （d）4000 5000 2000 2000 23000000
very recently in the framework of low-rank matrix analysis. 2000 1000 1000 1000
Since adjacent HSI bands usually exhibit strong correlations, -00.1 0 0.1 -00.1 0 0.1 0 -0.2 0 0.2 0 -0.2 0 0.2 0 -0.2 0 0.2
band8 band118 band189 band190 band191
byreshapingeachbandasavectorandstackingallthevectors
Fig.1. (a)DifferentbandsofimagesintheoriginalHSI.(b)RestoredHSIs
into a matrix, the reshaped HSI matrix can be assumed to
obtainedbyourproposedmethod.(c)Extractednoisebytheproposedmethod.
be with low rank. Various low-rank matrix models have been (d)Histogramsofthenoiseinallbands.
presented in different noise scenes in recent decades. Along
thisline,theclassicalLow-rankMatrixFactorization(LRMF)
model is presented by K.Mitra et al. [37] and T.Okatani et to further improve the capability of current HSI denoising
al. [36] for Gaussian noise, and its global optimal solution techniques especially under real complicated scenarios.
can be directly obtained by using Singular Value Decom- Tomakethispointclearer,let’strytoanalyzesomeevident
position (SVD) [18]. To add more robustness, the L -norm noise characteristics possessed by a HSI collected from real
1
LRMF [13], [15], [25], [27] is generally used and many cases. Fig. 1 presents a real HSI case for auxiliary under-
algorithms have been designed to solve this model, such as standing. Firstly, an image in a band is generally collected
L Wiberg [15], CWM [35] and RegL1ALM [57]. This L - under the similar sensor setting, and thus the i.i.d. distribution
1 1
normLRMFmodelactuallyassumesani.i.d.Laplaciannoises is a rational assumption for the noise over the elements in
embeddedindata.Tohandlemorecomplexnoisecases,Meng the band. This can be evidently observed from the band-
andDelaTorre[34]modeledthenoiseasmoreadaptablei.i.d. noise of the HSI image shown in Fig. 1. Secondly, due to
mixture of Gaussians (MoG) distributions to represent noise. the difference of the sensor sensitivity for collecting images
SuchnoisemodelingideawasfutherextendedtotheBayesian located in various HSI bands, the noise of different-band-
frameworkbyChenetal.[8],toRPCAbyZhaoetal.[55]and images always depicts evident distinctions [19]. From Fig.
to the tensor factorization by Chen et al. [49] . Very recently, 1, it is easy to see that images located in some bands are
Cao et al. [?], [4] modeled the noise as a more general i.i.d. very clean, while are extremely noisy in some others [16],
mixture of Exponential Power (MoEP) distribution, which [53]. Thirdly, albeit different, the noise distributions along
achievescompetitiveperformanceinHSIdenoisingunderreal different bands have certain correlation. E.g., along adjacent
noise scenarios. Besides, some other attempts have also been bands, the noise tends to be more or less similar since images
proposed to model the noise as a combination of sparse and on neighboring spectrums are generally collected under small
Gaussiannoise[1],andZhangetal.[54]alsoutilizedthisidea deviation of sensor settings and wavelength [58]. From Fig.
for HSI denoising task. Later, He et al. [21] further enhance 1, it is easy to see the noise similarity for images located in
the capability of the method by adding a TV regularization to 189-191bands.Accordingly,therealHSInoisedistributionis
low-rank reconstruction. Besides, based on the mixture noise generallynon-i.i.d.andhasamorecomplicatedconfigurations
assumption,Heet.al [19]proposedanoise-adjustedlow-rank than the i.i.d. noise assumption of the current HSI denoising
methods and Wang et al. [47] proposed a GLRR denoising techniques.Suchdeviationinclinestomaketheirperformance
methodforthereconstructionofthecorruptedHSIs.Allthese degenerate under more practical cases, which will be clearly
approaches also achieve good performance on HSI denoising observed in our experiments.
in mixed noise cases. To address this issue, this paper proposes a new noise
From Gaussian noise assumption to MoEP noise assump- modelingframework,bycarefullydesigningnoisedistribution
tion, such advancements make the model more adaptive to structure to make it possibly faithfully deliver the real HSI
variousHSInoisesencounteredinpractice.However,allofthe noise configurations. Specifically, we model the noise of
aforementionedLRMFmodelsjustsimplyassumethattheHSI each HSI band with different Mixture of Gaussian (MoG)
noiseisi.i.d.,whichismoreorlessdeviatedfromthepractical distributions (i.e., parameters of MoG are different). Besides,
scenarios, where the noises in a HSI is generally with non- MoG parameters across different bands are encoded under a
i.i.d. configurations. In this sense, there is still a large room similartop-levelpriordistribution,representingthecorrelation
JOURNALOFLATEXCLASSFILES,VOL.14,NO.8,AUGUST2015 3
betweennoisedistributionsacrossdifferentbands.Inthisway,
=1,…,
the non-i.i.d. noise structure under a practical HSI image can
be more properly encoded and an appropriate HSI denoising
effect is thus to be expected.
In this study, we embed such noise modeling framework
into the low-rank matrix factorization (LRMF) realization,
which easily assumes a low-rank structure (both spatial and j=1,…,
spectral) of the to-be-recovered HSI. Our experimental results
show that, even under this simple setting of HSI structure,
such noise modeling idea can greatly help enhance the HSI Fig.2. GraphicalmodelofNMoG-LRMF.Yij denotestheith HSIelement
denoising ability beyond previous techniques and obtain the in its jth band. ui and vj are columns of low-rank matrix U and V,
respectively, generated from a Gaussian distribution with precise γ, with
state-of-the-art performance in real HSI denoising tasks. Gammapriordistributionwithhyper-parametersξ0andδ0.Theportioninthe
The rest of the paper is organized as follows: The proposed redboxcorrespondstothenoiseencodingpart.µj,τj arethek-dimentional
vectorsrepresentingthemeanandvarianceofallMoGcomponentsinthejth
model and the corresponding variational inference algorithm
band,withhyper-parametersβ0,m0,d,c0,wherediswithhyper-parameters
are presented in Section 2. Experimental results are shown η0,λ0. zij is the hidden variable generated from Multinomial distribution
in Section 3. Finally, conclusions are drawn in Section 4. withparameterπj,withhyper-parametersα0.
Throughout the paper, we denote scalars, vectors, matrices,
tensors as non-bold letters, bold lower case letters, bold upper
τ are mean and precision of the kth Gaussian component in
case letters and decorated letters, respectively. jk
thejth band,respectively.NotethatMoGparametersπ ,µ
jk jk
and τ are different, implying different MoG distributions
jk
II. NON-I.I.D.MOGMETHODFORHSIDENOISING
across different bands.
In this section, we first introduce our non-i.i.d. noise en- Considering the noise property (iii), we further provide the
coding strategy and then propose a non-i.i.d. MoG LRMF hypothesis that MoG parameters µ s and τ s of all bands are
j j
(NMoG-LRMF) model by using this strategy in the LRMF generated from a two-level prior distribution:
model for HSI denoising. Finally, the corresponding varia-
µ ,τ ∼N(µ |m ,(β τ )−1)Gam(τ |c ,d),
tionalinferencealgorithmfortheproposedmodelisdesigned. jk jk jk 0 0 jk jk 0 (3)
d∼Gam(d|η ,λ ),
0 0
A. Non-i.i.d. noise encoding whereGam(·)representstheGammadistribution.Inthisway,
Let’s represent a given HSI as a matrix Y ∈RN×B, where the correlation between noise distributions among different
bandsisthenrationallyencoded.Throughintroducingalatent
N andB meanthespatialandspectraldimensionsoftheHSI,
variable z , we can equivalently rewrite Eq. (2) as the
respectively. By assuming that noise is addictive, we have the ijk
following form:
following expression:
K
Y =L+E, (1) e ∼ (cid:89)N(e |µ ,τ−1)zijk,
ij ij jk jk
where L ∈ RN×B denotes the clean HSI and E ∈ RN×B k=1 (4)
denotes the embedded noise. zij ∼Multinomial(zij|πj),
In the previous section, we have summarized three intrinsic π ∼Dir(π |α ),
j j 0
properties possessed by HSI noise: (i) Noise of an image
where Multinomial(.) and Dir(.) represent the Multinomial
locatedinthesamebandarei.i.d.;(ii)Noiseofimageslocated
and Dirichlet distributions, respectively. Then, Eq. (3) and (4)
indifferentbandsarenon-identical;(iii)Noisedistributionsin
together encode the noise structure embedded in a HSI. Fig.
different bands are correlated.
2 shows the graphical model for noise encoding within the
Based on the above properties (i) and (ii) of HSI noise
red box. All involved parameters can be inferred from data,
structure, we model noise located in each band as an inde-
as introduced in the following Section II.C.
pendent MoG distribution while assume that the parameters
for MoG distributions in different bands are different. The
MoG is utilized due to its universal approximation capability B. NMoG-LRMF model
for any continuous densities [33], which has been extensively
For the prior structure of the underlying clean HSI matrix
used and studied in [34], [55].
L, we readily employ the low-rank structure to encode it.
Denotee astheelementlocatedinith rowandjth column
ij Specifically, let’s consider the following LRMF model:
of the noise matrix E, and as aforementioned, the noise
distribution located in the jth band can then be modeled as: Y =UVT +E, (5)
K where L = UVT implies the low-rank structure underlying
(cid:88)
p(eij)= πjkN(eij|µjk,τj−k1), (2) the clean HSI.
k=1 For most deterministic LRMF model, the rank r of matrix
where π is the mixing proportion with π ≥ 0 and L is fixed. By modeling the problem into certain generative
jk jk
(cid:80)K π =1,K istheGaussiancomponentnumber,µ and model [1], [55], the rank can also be adaptively learned from
k=1 jk jk
JOURNALOFLATEXCLASSFILES,VOL.14,NO.8,AUGUST2015 4
data. Specifically, suppose the columns of U and V are form inference equation for each component of (11)1.
generated from Gaussian distribution. For l = 1,2,··· ,R, Estimation of noise component: We first list the updating
where R is a preset larger value beyond the true rank r, then equation for the noise components involved in Eq. (11). The
posterior distribution of mean and precision for noise in each
u ∼N(u |0,γ−1I ), v ∼N(v |0,γ−1I ), (6)
·l ·l l N ·l ·l l B band is updated by the following equation:
where u·l and v·l are the lth columns of U and V, respec- (cid:89) 1
tively. IN(IB) denotes the N×N(B×B) identity matrix. γl q∗(µj,τj)= N(µjk|mjk,β τ )Gam(τjk|cjk,djk), (12)
jk jk
is the precision of u and v with prior as follows: k
·l ·l
where the parameters in the above equation can be calculated
γ ∼Gam(γ |ξ ,δ ). (7)
l l 0 0 from data in the following way:
Note that each column pair u·l and v·l of U, V has the 1 (cid:88)
m = {β m (cid:104)z (cid:105)(Y −(cid:104)µ (cid:105))}, (13)
same sparsity profile characterized by the common precision jk βjk 0 0 i ijk ij jk
variable γ . It has been validated that such a modeling could
l
olefadauttoomlaragtiecaplrleycicsoionnduvcatliunegsloofws-ormanekγelss,tiamnadteheonfceLis[1c]a.pable βjk =β0+(cid:88)i(cid:104)zijk(cid:105), cjk =c0+ 12(cid:88)i(cid:104)zijk(cid:105), (14)
Combining Eqs. (2)-(7) together, we can construct the full d =(cid:104)d(cid:105)+ 1{(cid:88) (cid:104)z (cid:105)(cid:104)(Y −u vT)2(cid:105)+β m2
NMoG-LRMF Bayesian model. And the graphical model jk 2 i ijk ij i· j· 0 0
(15)
representationofthismodelisshowninFig.2.Thegoalturns − 1 [(cid:88) (cid:104)z (cid:105)(Y −u vT)+β2m2]2}.
to infer the posterior of all involved variables: βjk i ijk ij i· j· 0 0
The update equation of the latent variable Z is:
p(U,V,Z,µ,τ,π,γ,d|Y), (8)
where Z = {zij}N×B, µ = {µjk}B×K, τ = {τjk}B×K, q∗(zij)=(cid:89)k(cid:37)zijijkk, (16)
π =(π ,··· ,π ), γ =(γ ,··· ,γ ).
1 B 1 R where the involved parameters are calculated by
(cid:88)
(cid:37) =ρ / ρ ,
ijk ijk ijk
k √
C. Variational Inference lnρ =(cid:104)lnπ (cid:105)−ln 2π+(cid:104)lnτ (cid:105)/2 (17)
ijk jk jk
We use variational Bayes (VB) method [3] to infer the −(cid:104)τjk(Yij −µjk−ui·vjT·)2(cid:105)/2.
posterior.Specifically,VBaimstouseavariationaldistribution
Similarly, theupdate equationfor themixing proportion π
j
q(θ) to approximate the true posterior p(θ|D), where θ over the jth band can be written as:
denotes the set of parameters and D denotes the observed
data. To achieve this goal, we need to solve the following q∗(πj)=(cid:89)kπjαkjk−1, (18)
optimization problem:
(cid:80)
where α =α + (cid:104)z (cid:105).
(cid:90) p(θ|D) jk 0k i ijk
minKL(q(θ)||p(θ|D))=− q(θ)ln{ }dθ, (9) The update equation on the hyper-parameter d is:
q∈C q(θ)
q∗(d)=Gam(d|η,λ), (19)
where KL(q||p) represents the KL divergence between two
distribution q and p, and C denotes the constraint set of prob- where η =η0+c0KB and λ=λ0+(cid:80)j,k(cid:104)τjk(cid:105).
abilitydensitiestomaketheminimizationtractable.Assuming Estimation of low-rank component.Completingthecom-
C is the distribution family which can be factorized with ponentupdatefornoise,wethenestimatetheposterioroflow-
(cid:81)
respecttosomedisjointgroups:q(θ)= q (θ ),theclosed- rankcomponentu (i=1,··· ,N)andv (j =1,··· ,B)as:
i i i i· j·
form optimal solution q∗(θ ) can be obtained by
i i q∗(u )=N(u |µ ,Σ ), (20)
i· i· ui· ui·
exp{(cid:104)lnp(θ,D)(cid:105) }
qi∗(θi)= (cid:82) exp{(cid:104)lnp(θ,D)θ(cid:105)\θi }, (10) where
θ\θi (cid:88)
µ ={ (cid:104)z (cid:105)(cid:104)τ (cid:105)(Y −(cid:104)µ (cid:105))(cid:104)v (cid:105)}Σ ,
where (cid:104)·(cid:105) denotes the expectation and θ\θi denotes the set ui· j,k ijk jk ij jk j· ui·
ofTθhewnitwheθciarnemanoavlyetdi.callyinferallthefactorizeddistributions Σui· ={(cid:88)j,k(cid:104)zijk(cid:105)(cid:104)τjk(cid:105)(cid:104)vjT·vj·(cid:105)+(cid:104)Γ(cid:105)}−1.
involved in Eq. (8). Suppose that the approximation of poste- q∗(v )=N(v |µ ,Σ ), (21)
rior distribution (8) possesses a factorized form as follows: j· j· vj· vj·
(cid:89) (cid:89) where
q(U,V,Z,µ,τ,π,γ,d)= q(u ) q(v )
i· j· (cid:88)
µ ={ (cid:104)z (cid:105)(cid:104)τ (cid:105)(Y −(cid:104)µ (cid:105))(cid:104)u (cid:105)}Σ ,
(cid:89) (cid:89) i (cid:89) j (11) vj· (cid:88)i,k ijk jk ij jk i· vj· (22)
q(zij) q(µj,τj)q(πj) q(γl)q(d), Σvj· ={ i,k(cid:104)zijk(cid:105)(cid:104)τjk(cid:105)(cid:104)uTi·ui·(cid:105)+(cid:104)Γ(cid:105)}−1,
ij j l
where u (v ) are the ith (jth) row of matrix U (V),
i· j· 1Inference details to obtain these equations can be referred to in
respectively. According to Eq. (10), we can get the closed- http://dymeng.gr.xjtu.edu.cn.
JOURNALOFLATEXCLASSFILES,VOL.14,NO.8,AUGUST2015 5
where Γ = diag((cid:104)γ(cid:105)). For γ which controls the rank of U Besides, the performance of TDL [39] and BM4D [32] are
l
and V, we have : also compared, and both methods represent the state-of-the-
art methods for HSI denoising by considering HSI priors. All
q∗(γ )=Gam(γ |ξ ,δ ), (23) experiments were implemented in Matlab R2014b on a PC
l r l l
with 4.0GHz CPU and 31.4GB RAM.
where
ξ =ξ +(m+n)/2,
l 0 A. Simulated HSI denoising experiments
(cid:88) (cid:88)
δ =δ + (cid:104)u2(cid:105)/2+ (cid:104)v2(cid:105)/2.
l 0 il jl Inthisexperiment,wefocusontheperformanceofNMoG-
i j
LRMR in HSI denoising with syntheic noise. Two HSIs were
mandnrepresenttheimagesizeamongthespatialdimension.
employed:WashingtonDCMall2withsizeof1208×307×191
As discussed by Babacan et al. [1], some γ s tend to be
l and RemoteImage3 provided by Liu et. al. [29] with size of
very large during the inference process and the corresponding
205 × 246 × 96. After cropping the main part of HSI and
rows will be removed from U and V. The low-rank purpose
deleting some evident visual contaminative spectral channels,
can thus be rationally conducted. In all our experiments, we
Washington DC Mall and RemoteImage are resized to 200×
just automatically infer the rank of the reconstructed matrix
200×160 and 200×200×89, respectively. The gray value
through throwing away those comparatively very large γ as
l of each band are normalized into [0,1].
previous literatures did [55].
Real-world HSIs are usually contaminated by several dif-
The proposed variational inference method for the NMoG-
ferent types of noise, including the most common Gaussian
LRMF model can then be summarized in Algorithm 1.
noise, impulse noise, dead pixels or lines, and stripes [54]. In
ordertosimulatetheserealHSInoisescenarios,weaddedsix
Algorithm 1 NMoG-LRMF algorithm
kinds of noises to the original HSI data.
Input: the original HSI matrix Y ∈ RN×B, Gaussian component
(1)I.i.d.Gaussiannoise:Entriesinallbandswerecorrupted
number K, and maximum iteration number.
by zero-mean i.i.d. Gaussian noise N(0,σ2) with σ =0.05.
Output: U =U , V =V .
opt t opt t
Initialization: Parameters (m ,β ,c ,d ,η ,λ ) in noise prior. (2)Non-i.i.d.Gaussiannoise:Entriesinallbandswerecor-
0 0 0 0 0 0
Low-rank components U ,V and parameters in model prior rupted by zero-mean Gaussian noise with different intensity.
0 0
(ξ0,δ0); t=1. The signal noise ratio (SNR) value of each band is generated
1: while not coverged do fromuniformdistributionwithvalueintherangeof[5,10]dB.
2: UpdateapproximateposteriorofnoisecomponentZt,πt by
(3) Gaussian + Stripe noise: All bands were corrupted by
Eq. (16)-(18).
3: Updateapproximateposteriorofnoisecomponentµt,τt by GaussiannoiseasCase(2).Besides,40bandsinDCmalldata
Eq. (12)-(15). (30bandsinRemoteImagedata)wererandomlychosentoadd
4: Update approximate posterior of noise component dt by Eq. stripe noise. The number of stripes in each band is from 20
(19).
to 40.
5: Update approximate posterior of low-rank component
Ut,Vt by Eq. (20)-(22). (4) Gaussian + Deadline noise: Each band was contam-
6: Updateapproximateposteriorofparametersinnoisecompo- inated by Gaussian noise as Case (2). 40 bands in DCmall
nent γt by Eq. (23). data (30 bands in RemoteImage data) were chosen randomly
7: t←t+1. to add deadline noise. The number of deadline is from 5 to
8: end while
15.
(5)Gaussian+Implusenoise:Allbandswerecorruptedby
Setting of hyper-parameters: We set all the hyper-
Gaussian noise as Case (2). 40 bands in DCmall (30 bands in
parametersinvolvedinourmodelinanon-informativemanner
RemoteImage data) were randomly chosen to added impluse
to make them possibly less affect the inference of posterior
noise with different intensity, and the percentage of impluse
distributions[3].Throughoutourexperiments,wesetm as0,
0 is from 50% to 70%.
andβ ,c ,d ,η ,λ ,ξ ,δ asasmallvalue10−3.Ourmethod
0 0 0 0 0 0 0 (6)Mixturenoise:Eachbandwasrandomlycorruptedbyat
performs stably well under such easy settings.
least one kind of noise mentioned in case (2)-(5).
Three criteria were utilized to measure performance: (1)
III. EXPERIMENTALRESULTS MPSNR[24]:Meanofpeaksignal-to-noiseratio(PSNR)over
all bands between clean HSI and recovered HSI. (2) MSSIM
In this section, to evaluate the performance of the proposed
[48]: Mean of structural similarity (SSIM) between clean HSI
NMoG-LRMF method, we conducted a series of experiments
andrecoveredHSIoverallbands.(3)Time:Timecostofeach
on both synthetic and real HSI data. Compared methods
method used to complete the denoising process.
includeLRMR[54]andLRTV[21],consideringdeterministic
The parameters of competing methods are set as follows:
Gaussian noise and sparse noise. Meanwhile, five representa-
the block size in LRMR is 20×20×B and the step size
tive low-rank matrix analysis methods considering different √
is 10. For LRTV, (cid:15) = (cid:15) = 10−8 and λ = 1/ NB. For
kinds of i.i.d. noise distributions were also considered for 1 2
TDLd, we set the block size as 6×6×B and the step size
comparison, including PMoEP [4] (assuming i.i.d. mixture of
is 2. For BM4D, we set the block size as 8 × 8 × 8 and
Exponential Power noise), MoG-RPCA [55] (assuming i.i.d.
MoG noise), RegL1ALM [57], CWM [35] (assuming i.i.d. 2http://engineering.purdue.edu/∼biehl/MultiSpec/hyperspectral.html
Laplacenoise)andSVD[18](assumingi.i.d.Gaussiannoise). 3http://peterwonka.net/Publications/code/LRTC Package Ji.zip
JOURNALOFLATEXCLASSFILES,VOL.14,NO.8,AUGUST2015 6
TABLEI
PERFORMANCECOMPARISONOFALLCOMPETINGMETHODSONDCMALLDATAWITHSYNTHETICNOISE.
NoisyHSI SVD RegL1ALM CWM MoG-RPCA PMoEP LRMR LRTV TDL BM4D NMoG
i.i.d.GaussianNoise
MPSNR 26.02 44.76 41.19 42.89 45.26 44.26 42.48 43.83 46.08 39.33 45.27
MSSIM 0.513 0.983 0.965 0.973 0.985 0.980 0.972 0.977 0.989 0.948 0.985
time - 0.290 221.2 851.1 59.1 201.7 5367.0 512.3 44.9 965.2 77.8
Noni.i.d.GaussianNoise
MPSNR 40.96 43.34 47.15 45.00 46.15 44.09 48.26 47.19 41.80 45.50 51.11
MSSIM 0.890 0.973 0.985 0.980 0.981 0.973 0.987 0.986 0.903 0.983 0.990
time - 0.292 212.0 800.2 121.5 1672.3 778.6 478.9 821.2 940.8 106.7
Gaussian+StripeNoise
MPSNR 44.72 46.14 48.89 48.74 48.54 48.08 51.11 49.06 44.94 46.66 54.16
MSSIM 0.941 0.989 0.992 0.992 0.992 0.990 0.994 0.993 0.942 0.985 0.997
time - 0.515 238.5 782.6 307.6 1505.5 414.4 683.2 843.1 369.5 334.7
Gaussian+Deadlinenoise
MPSNR 44.36 44.08 48.95 47.81 48.65 46.95 50.34 48.95 44.84 46.02 54.04
MSSIM 0.938 0.980 0.993 0.990 0.991 0.986 0.991 0.993 0.943 0.979 0.997
time - 0.291 177.6 690.4 136.2 2004.7 670.8 321.2 772.3 307.5 143.6
Gaussian+ImpluseNoise
MPSNR 43.05 44.05 49.16 47.51 49.89 46.79 51.05 49.11 43.65 45.77 52.96
MSSIM 0.925 0.981 0.993 0.991 0.993 0.986 0.994 0.993 0.932 0.977 0.997
time - 0.293 181.3 567.4 146.5 1640.9 743.6 501.3 695.3 1028.6 150.2
MixtureNoise
MPSNR 42.93 43.32 48.75 44.82 47.94 44.26 48.77 48.54 43.22 43.97 50.38
MSSIM 0.913 0.977 0.991 0.982 0.989 0.978 0.986 0.991 0.917 0.968 0.996
time - 0.289 176.5 688.7 136.3 940.3 708.2 358.9 653.8 325.4 142.0
TABLEII
PERFORMANCECOMPARISONOFALLCOMPETINGMETHODSONREMOTEIMAGEDATAWITHSYNTHETICNOISE.
NoisyHSI SVD RegL1ALM CWM MoG-RPCA PMoEP LRMR LRTV TDL BM4D NMoG
i.i.d.GaussianNoise
MPSNR 26.02 38.19 36.63 36.78 38.66 36.72 37.46 37.89 39.81 38.98 38.65
MSSIM 0.591 0.956 0.937 0.939 0.964 0.956 0.946 0.949 0.971 0.971 0.964
time - 0.247 105.9 331.6 65.7 38.1 409.1 211.9 47.7 357.7 69.6
Noni.i.d.GaussianNoise
MPSNR 25.17 34.36 35.51 35.37 36.41 35.67 35.95 35.73 29.74 36.82 38.26
MSSIM 0.539 0.902 0.921 0.922 0.942 0.928 0.926 0.927 0.729 0.952 0.962
time - 0.253 105.4 329.3 100.2 301.9 398.7 220.2 151.0 327.8 101.1
Gaussian+StripeNoise
MPSNR 29.08 37.48 39.31 38.64 39.70 38.62 40.13 39.36 29.76 37.03 40.68
MSSIM 0.710 0.955 0.968 0.964 0.974 0.965 0.971 0.968 0.731 0.938 0.983
time - 0.279 108.9 375.8 88.6 299.8 606.4 238.3 1969.9 160.4 112.3
Gaussian+Deadlinenoise
MPSNR 27.97 33.37 39.08 37.94 39.58 38.20 38.74 38.72 29.98 34.02 40.43
MSSIM 0.693 0.907 0.967 0.961 0.972 0.963 0.958 0.965 0.772 0.888 0.983
time - 0.299 108.9 361.7 119.6 1205.9 725.1 378.9 273.5 683.6 127.5
Gaussian+ImpluseNoise
MPSNR 28.20 35.00 41.89 39.48 41.15 39.15 41.32 40.63 30.63 36.97 42.90
MSSIM 0.660 0.917 0.983 0.971 0.981 0.969 0.977 0.977 0.745 0.951 0.990
time - 0.236 100.0 317.2 110.1 202.8 473.1 251.6 190.7 445.3 118.1
MixtureNoise
MPSNR 25.84 31.55 37.93 35.72 38.68 36.09 37.76 38.31 26.93 30.29 39.32
MSSIM 0.590 0.868 0.958 0.936 0.966 0.940 0.944 0.962 0.633 0.782 0.979
time - 0.283 122.0 416.4 69.5 650.7 696.4 260.8 312.2 86.7 105.0
the step size is 4. For NMoG-LRMF method, The component methods.Ini.i.d.Gaussiancase,insteadofonlyusingthesim-
number K in each band was fixed as 1 in Case (1)-(2), and 3 plelow-rankpriorinourmethod,multiplecompetingmethods,
in Case (3)-(6). The rank of all low-rank based methods is set like TDL and BM4D, utilize more useful HSI priors in their
as 5 in DCmall data experiment and 4 in RemoteImage data model, and thus tend to have relatively better performance.
experiment.Allparametersinvolvedinthecompetingmethods While on more complex but more practical complicated non-
wereempiricallytunedorspecifiedassuggestedbytherelated i.i.d. noise cases, the advantage of the proposed method is
literatures to guarantee their possibly good performance. All evident.Thiscanbeeasilyexplainedbythebetternoisefitting
competingmethodsrunwith20randominitializationsineach capability of the proposed method, i.e., it can more properly
noise case, and the average result is reported. extractnoisesembeddedinHSI,whichthennaturallyleadsto
The results of all competing methods in DCmall and its better HSI recovery performance.
RemoteImage HSI data are shown in Table I and Table II, Furthermore, it also can be seen that the computational
respectively. The superiority of the proposed method can be cost of the proposed method is with almost similar order of
easilyobserved,exceptinthei.i.d.Gaussiannoisecase,which magnitude with other competing methods, except the known
complies with the basic noise assumption of conventional SVD, which we use the mature toolkit in Matlab and can be
JOURNALOFLATEXCLASSFILES,VOL.14,NO.8,AUGUST2015 7
SVD RegL1ALM CWM MoG RPCA PMoEP LRMR LRTV TDL BM4D NMoG LRMF
45 45 45 45 50 45
45 40 40 40 45 40
40 35 40 35
PSNR35 40 PSNR3305 PSNR3305 PSNR2350 PSNR233505 PSNR2350
35 25 25 20 20 20
30 20 40 60 80 100 120 140 160 20 20 40 60 80 100 120 140 160 20 20 40 60 80 100 120 140 160 15 20 40 60 80 100 120 140 160 20 40 60 80 100 120 140 160 15 20 40 60 80 100 120 140 160
Band number R Band number Band number Band number Band number Band number
N30
S
00000.....099999.888999468249 20 40 60 80 100 12P0 14205160 000...00999..7899955589 20 40 60 80 100 120 140 160 000...00999..7899955589 20 40 60 80 100 120 140 160 00..9968 20 40 60 80 100 120 140 160 0.981 20 40 60 80 100 120 140 160 00..9968 20 40 60 80 100 120 140 160
1 20 1 1 1 1 1
SSIM0000....99992468 15 SSIM00..089.559 20 40SSIM000...000789...555789 60 Band 8nSSIM0u000...00789..55589mber 100 1SSIM0002...00789..555890 140 SSIM00001....678960
0.9 20 40 60Band n80umber100 120 140 160 0.8 20 40 60Band n80umber100 120 140 160 0.65 20 40 60Band n80umber100 120 140 160 0.7 20 40 60Band n80umber100 120 140 160 0.7 20 40 60Band n80umber100 120 140 160 20 40 60Band n80umber100 120 140 160
(a) I.i.d. Gaussian noise (b) Non i.i.d. Gaussian noise (c) Gaussian + Stripe noise (d) Gaussian+ Deadline noise (e) Gaussian + Impluse noise (f) Mixture noise
Fig.3. EachcolumnshowstheaveragePSNRandSSIMmeasurementsamong20initializationsofallmethodsundercertaintypeofnoiseinDCmalldata:
(a)Gaussiannoise.(b)Gaussian+Stripenoise.(c)Gaussian+Deadlinenoise.(d)Gaussian+Implusenoise.(e)Mixturenoise.Thedemarcatedareaofthe
subfigureindicatesthecurvelocalityonalargerscale.
SVD RegL1ALM CWM MoG RPCA PMoEP LRMR LRTV TDL BM4D NMoG LRMF
42 42 40 40 50 40
40 40
38 45 38 35 45 35
PSNR333246 40 PSNR333246 PSNR3305 PSNR2350 PSNR334050 PSNR2350
2380 20 40 60 3805 2380 20 40 60 80 25 20 40 60 80 20 20 40 60 80 2205 20 40 60 80 20 20 40 60 80
Band number R Band number Band number Band number Band number Band number
N30
0000...00999...9678995555781 20 40 60 PS228005 00..000099....789999556878 20 40 60 80 000..0099...97899555781 20 40 60 80 0000...00999.0..967899.55559781 10 20 30 40 50 60 70 80 0.09.5911 20 40 60 80 0000...999.96781 10 20 30 40 50 60 70 80
SSIM0.9 15 SSIM000...0899.8249 20 40SSIM0.008..589 60 8SSIM000..0078..5578 100 1SSIM200..780 140 SSIM0001...67860
0.85 0.86 0.75 Band nu0.65mber 0.6 0.5
0.8 20 Ban4d0 number60 80 0.84 20 Ban4d0 number60 80 0.7 20 Ban4d0 number60 80 0.6 10 20 30Ba4n0d num5b0er 60 70 80 0.5 20 Ban4d0 number60 80 0.4 10 20 30Ba4n0d num5b0er 60 70 80
(a) I.i.d. Gaussian noise (b) Non i.i.d. Gaussian noise (c) Gaussian + Stripe noise (d) Gaussian+ Deadline noise (e) Gaussian + Impluse noise (f) Mixture noise
Fig.4. EachcolumnshowstheaveragePSNRandSSIMmeasurementsamong20initializationsofallmethodsundercertaintypeofnoiseinRemoteImage
data:(a)Gaussiannoise.(b)Gaussian+Stripenoise.(c)Gaussian+Deadlinenoise.(d)Gaussian+Implusenoise.(e)Mixturenoise.
implemented very efficiently. Considering its better capability better reconstruction effect both visually and quantitatively,
in fitting much wider range of noises than current methods, it and both finely removes the unexpected noises, and better
shouldberationaltosaythattheproposedmethodisefficient. recovers HSI details like edges and textures.
WefurthershowthePSNRandSSIMmeasurementsacross
all bands of the HSI under six types of noise settings in two B. Real HSI denoising experiments
experimentsinFig.3andFig.4,respectively.Fromthefigures,
Inthissection,weevaluatetheperformanceoftheproposed
it is easy to see that TDL obtains the best PSNR and SSIM
method on two real HSI datasets, Urban dataset4 and EO-1
values across all bands in the i.i.d. Gaussian noise. In other
Hyperion dataset5.
more complex noise cases, NMoG-LRMF achieves the best
The similar competing methods as the last section were
PSNR and SSIM values across almost all bands. This verifies
compared in the real HSI experiments. Like the synthetic
the robustness of the proposed method over entire HSI bands.
experiments, all involved parameters have been finely tuned
Figs.5and6givetherestorationresultsoftwotypicalbands
or directly suggested by the original references to promise
in DCmall and RemoteImage HSIs, respectively. The original
a possibly well performance of all competing methods. The
HSIs are corrupted by mixture of Gaussian, stripes, deadline
component number of each band in the proposed method is
and impluse noise. It can be observed that the competing
easily set as 3 throughout all experiments. The gray value of
methodsSVD,CWM,PMoEP,LRMR,andBM4Dcanhardly
HSI were normalized into [0,1].
remove this complex mixture noise from the images. We can
The first Urban dataset is of the size 307 × 307 × 210,
also see that although RegL1ALM, MoG-RPCA and LRTV
and some bands are seriously polluted by atmosphere and
have a better denoising performance, the restored HSI still
water.Weuseallofdatawithoutremovinganybandstomore
relativelyblurmanydetailsascomparedwiththegroundtruth,
which are also reflected from their PSNR and SSIM values in 4http://www.tec.army.mil/hypercube
the band. Comparatively, our NMoG-LRMF method achieves 5http://datamirror.csdb.cn/admin/dataEO1Main.jsp
JOURNALOFLATEXCLASSFILES,VOL.14,NO.8,AUGUST2015 8
(a) Original HSI (b) Noisy HSI (c) SVD (d) RegL1ALM (e) CWM (f) MoG-RPCA
PSNR=15.19, SSIM=0.414 PSNR=31.29, SSIM=0.944 PSNR=20.36, SSIM=0.680 PSNR=31.19, SSIM=0.955 PSNR=30.71, SSIM=0.958
(g)PMoEP (h) LRMR (i) LRTV (j) TDL (k) BM4D (l) NMoG-LRMF
PSNR=30.31, SSIM=0.951 PSNR=24.67, SSIM=0.845 PSNR=25.80, SSIM=0.847 PSNR=16.45, SSIM=0.483 PSNR=23.07, SSIM=0.780 PSNR=32.86 SSIM=0.972
Fig. 5. Restoration results of band 75 under mixture noise in DCmall data: (a) Original HSI, (b) Noisy HSI, (c) SVD, (d) RegL1ALM, (e) CWM, (f)
MoG-RPCA,(g)PMoEP,(h)LRMR,(i)LRTV,(j)TDL,(k)BM4D,(l)NMoG.Thisfigureshouldbebetterseenbyzoomingonacomputerscreen.
(a) Original HSI (b) Noisy HSI (c) SVD (d) RegL1ALM (e) CWM (f) MoG-RPCA
PSNR=12.62, SSIM=0.112 PSNR=18.41, SSIM=0.326 PSNR=16.74, SSIM=0.241 PSNR=25.13, SSIM=0.712 PSNR=32.67, SSIM=0.902
(g)PMoEP (h) LRMR (i) LRTV (j) TDL (k) BM4D (l) NMoG-LRMF
PSNR=32.63, SSIM=0.907 PSNR=22.06, SSIM=0.576 PSNR=16.88, SSIM=0.244 PSNR=13.16, SSIM=0.126 PSNR=15.77, SSIM=0.226 PSNR=34.06, SSIM=0..944
Fig.6. Restorationresultsofband86undermixturenoiseinRemoteImagedata:(a)OriginalHSI,(b)NoisyHSI,(c)SVD,(d)RegL1ALM,(e)CWM,(f)
MoG-RPCA,(g)PMoEP,(h)LRMR,(i)LRTV,(j)TDL,(k)BM4D,(l)NMoG.
（a）Original HSI （b）SVD （c）RegL1ALM （d）CWM （e）MoG-RPCA （f）PMoEP
（g）LRMR （h）LRTV （i）TDL （j）BM4D （k）NMoG
Fig.7. Restorationresultsofallmethodsonband103inUrbanHSIdata.(a)OriginalHSI,(b)SVD,(c)RegL1ALM,(d)CWM,(e)MoG-RPCA,(f)PMoEP,
(g)LRMR,(h)LRTV,(i)TDL,(j)BM4D,(k)NMoG.
JOURNALOFLATEXCLASSFILES,VOL.14,NO.8,AUGUST2015 9
（a）Original HSI （b）SVD （c）RegL1ALM （d）CWM （e）MoG-RPCA （f）PMoEP
（g） （h） （i） （j）
（g）LRMR （h）LRTV （i）TDL （j）BM4D （k）NMoG
Fig.8. Restorationresultsofallmethodsonband139inUrbanHSIdata.(a)OriginalHSI,(b)SVD,(c)RegL1ALM,(d)CWM,(e)MoG-RPCA,(f)PMoEP,
(g)LRMR,(h)LRTV,(i)TDL,(j)BM4D,(k)NMoG.
（a）Original HSI （b）SVD （c）RegL1ALM （d）CWM （e）MoG-RPCA （f）PMoEP
（g）LRMR （h）LRTV （i）TDL （j）BM4D （k）NMoG
Fig.9. Restorationresultsofallmethodsonband207inUrbanHSIdata.(a)OriginalHSI,(b)SVD,(c)RegL1ALM,(d)CWM,(e)MoG-RPCA,(f)PMoEP,
(g)LRMR,(h)LRTV,(i)TDL,(j)BM4D,(k)NMoG.
rationally verify the robustness of the proposed method in the represents the mean digital number value of each row. As
presence of such heavy noises. shown in Fig. 10(a), due to the existence of mixed noise,
there are rapid fluctuations in the curve. After the restoration
Figs. 7, 8 and 9 present bands 103, 139 and 207 of the re-
processing, the fluctuations are more or less suppressed. It is
storedimagesobtainedbyallcompetingmethods,respectively.
easytoobservethattherestorationcurveofourNMoG-LRMF
From Fig. 7-9 (a), we can see that the original HSI bands are
method provides evidently smoother curves.
contaminatedbycomplexstructuralnoise,likethestripenoise.
It can be easily observed from Fig. 7-9 (b)-(f) that the LRMF ThesecondHyperiondatasetisofthesize3128×256×242
methods with different i.i.d. noise distribution assumptions andsomebandsaresoseriouslypollutedthatallsignaturesof
cannotfinelyrestoreacleanHSI.Thiscanbeeasilyexplained pixel is equal 0. Thus we only use a subset of its 198 bands
by the fact that such structural noises is obvious non-i.i.d. in this experiment after removing those zero signatures bands
andthedeviationbetweentherealnoiseconfigurationandthe 1-7, 58-76, and 224-242, which are evidently outliers with-
encoded knowledge in the model then naturally degenerate out any intrinsically useful information. We further spatially
theperformanceofthecorrespondingmethods.Comparatively, croppedasquareareafromtheHSIwithsize256×256×198
albeitconsideringlesspriorknowledgeontheto-be-recovered withrelativelyevidentnoisestospecificallytestifythedenois-
HSI, our method can still achieve a better recovery effect ing capability of a utlized method.
in visualization in its better restoration of texture and edge
Figs. 11, 12 and 13 show bands 100, 144 and 197 of
details and less preservation of structural noises due to its
the restored HSI images by all competing methods. From
powerful noise modeling ability. This further substantiates the
the figures, it can be easily observed that SVD and TDL
robustness of the proposed method in practical scenarios.
methodshavelittleeffectsinremovingthenoise;RegL1ALM
Then we give some quantitative comparison by showing and CWM methods can only partially eliminate the noise and
the horizontal mean profiles of band 207 in Urban dataset MoG-RPCA and PMoEP methods have a better performance
before and after restoration in Fig. 10. The horizontal axis in than SVD, RegL1ALM and CWM due to their better noise
the figure represents the row number, and the vertical axis fitting ability brought by the MoG and MoEP models. They,
JOURNALOFLATEXCLASSFILES,VOL.14,NO.8,AUGUST2015 10
0.4 0.4 0.4 0.4 0.4 0.4
0.35 0.35 0.35 0.35 0.35 0.35
0.3 0.3 0.3 0.3 0.3 0.3
0.25 0.25 0.25 0.25 0.25 0.25
0.2 0.2 0.2 0.2 0.2 0.2
0.15 0.15 0.15 0.15 0.15 0.15
0.1 0.1 0.1 0.1 0.1 0.1
0.05 0.05 0.05 0.05 0.05 0.05
0 50 10B0and 15n0umbe2r00 250 300 0 50 10B0and 15n0umbe2r00 250 300 0 50 10B0and 15n0umbe2r00 250 300 0 50 10B0and 15n0umbe2r00 250 300 0 50 10B0and 15n0umbe2r00 250 300 0 50 10B0and 15n0umbe2r00 250 300
（a）Original HSI （b）SVD （c）RegL1ALM （d）CWM （e）MoG-RPCA （f）PMoEP
0.4 0.4 0.4 0.4 0.4
0.35 0.35 0.35 0.35 0.35
0.3 0.3 0.3 0.3 0.3
0.25 0.25 0.25 0.25 0.25
0.2 0.2 0.2 0.2 0.2
0.15 0.15 0.15 0.15 0.15
0.1 0.1 0.1 0.1 0.1
0.05 0.05 0.05 0.05 0.05
0 50 10B0and 15n0umbe2r00 250 300 0 50 10B0and 15n0umbe2r00 250 300 0 50 10B0and 15n0umbe2r00 250 300 0 50 10B0and 15n0umbe2r00 250 300 0 50 10B0and 15n0umbe2r00 250 300
（g）LRMR （h）LRTV （i）TDL （j）BM4D （k）NMoG
Fig.10. Horizontalmeanprofilesofband207inUrbanHSIdata.(a)OriginalHSI,(b)SVD,(c)RegL1ALM,(d)CWM,(e)MoG-RPCA,(f)PMoEP,(g)
LRMR,(h)LRTV,(i)TDL,(j)BM4D,(k)NMoG.
（a）Original HSI （b）SVD （c）RegL1ALM （d）CWM （e）MoG-RPCA （f）PMoEP
（g）LRMR （h）LRTV （i）TDL （j）BM4D （k）NMoG
Fig. 11. Restoration results of all methods on band 100 in EO-1 Hyperion data. (a) Original HSI, (b) SVD, (c) RegL1ALM, (d) CWM, (e) MoG-RPCA,
(f)PMoEP,(g)LRMR,(h)LRTV,(i)TDL,(j)BM4D,(k)NMoG.
however, still miss lots of details in their HSI recovery due 43 0.995 0.98
42 37.5 0.979
tbouttihoenisruimndpelricHitSiIm. pArlolpoetrhei.ri.cdo.maspseutminpgtimonesthoodnsnaolissoecdainsntroi-t MPSNR344901 MSSIM0.09.8959 MPSNR363.57 MSSIM000...999777678
ahcahvieevceonasisdaetirsefdacmtooryreHHSSIIdpenriooirsiknngorwesleudltgseevinenthtehiorumghodtheelsy. 33780 2 Com4pone(nat 6)N umber8 10 0.980 2Com4ponen(tb 6N)u mber8 10 353.560 2 Com4ponen(ta 6N)u mber8 10 00..9977450 2Com4ponen(tb 6N)u mber8 10
Comparatively, the superiority of the proposed method can be Fig.15. StabilitytestoftheproposedmethodtoGaussiancomponentnumber
easilyobservedinbothdetailpreservationandnoiseremoving. K on DC Mall data in Gaussian + impluse noise and mixture noise cases.
(a) MPSNR tendency curve with respect to K. (b) MSSIM tendency curve
Fig.14showstheverticalmeanprofilesofband100before
withrespecttoK.
andafterrestoration.AsshowninFig.14(a),duetotheexis-
tence of mixed noise especially the stripes disorderly located
intheimage,thereareevidentfluctuationsovervariousplaces method tends to be stable and not very sensitive to the choise
of the curve. After the restoration processing, all competing of K value. Actually, in all our real experiments, we just
methodsshowevidentnon-smoothnessoverthecorresponding simply set K as 3, and our method can consistently perform
curves, except that obtained from the proposed one. This well throughout all our experiments.
implies that the HSI obtained by proposed NMoG-LRMF
method can better preserve such prior knowledge possessed IV. CONCLUSION
by real HSIs, as also substantiated by Fig. 10. In this paper, we initially propose a strategy to model the
HSInoiseusinganon-i.i.d.noiseassumption.Thenweembed
such noise modeling strategy into the low-rank matrix factor-
C. Effect of component number on denoising performance
ization (LRMF) model and propose a non-i.i.d LRMF model
Inthissection,weexaminethesensitivityofNMoG-LRMF under the Bayesian framework. A variational Bayes algorithm
method to the setting of the Gaussian component number K. is presented to infer the posterior of the proposed model.
We run NMoG-LRMF with 20 initializations on the DC Mall Compared with the current state-of-the-arts techniques, the
data in Gaussian + impluse noise and mixture noise cases, proposedmethodperformsmorerobustduetoitscapabilityof
respectively, with K varying from 1 to 10. The results are adapting to various noise shapes encountered in applications,
shown in Fig. 15. It can be easily observed that after K which is substantiated by our experiments implemented on
is larger than 2, the denoising performance of the proposed synthetic and real noisy HSIs.
