Unsupervised Unmixing of Hyperspectral Images Accounting for Endmember Variability

1 Unsupervised Unmixing of Hyperspectral Images Accounting for Endmember V ariability Abderrahim Halimi, Nicolas Dobigeon and Jean-Yves T ourneret Uni versity of T oulouse, IRIT/INP-ENSEEIHT/TéSA, T oulouse, France {Abderrahim.Halimi,Nicolas.Dobigeon,Jean-Yves.Tourneret}@enseeiht.fr Abstract This paper presents an unsupervised Bayesian algorithm for hyperspectral image unmixing accounting for endmember variability . The pixels are modeled by a linear combination of endmembers weighted by their corresponding abundances. Howe ver , the endmembers are assumed random to take into account their variability in the image. An additiv e noise is also considered in the proposed model generalizing the normal compositional model. The proposed algorithm exploits the whole image to provide spectral and spatial information. It estimates both the mean and the cov ariance matrix of each endmember in the image. This allows the beha vior of each material to be analyzed and its variability to be quantified in the scene. A spatial segmentation is also obtained based on the estimated abundances. In order to estimate the parameters associated with the proposed Bayesian model, we propose to use a Hamiltonian Monte Carlo algorithm. The performance of the resulting unmixing strategy is ev aluated via simulations conducted on both synthetic and real data. I . I N T RO D U C T I O N Hyperspectral imaging is a remote sensing technology that collects 3 dimensional data cubes composed of 2 D spatial images acquired in numerous contiguous spectra bands. Due to the limited spatial resolution of the observed image, each pix el generally consists of sev eral physical elements that are linearly [1], [2] or nonlinearly [3]–[5] mixed. Spectral unmixing (SU) consists of decomposing the pixel spectrum to recov er these materials, kno wn as endmembers , and estimating the corresponding proportions or abundances [6]. The linear mixture model (LMM) has receiv ed great interest in Part of this work has been funded by the Hypanema ANR Project n ◦ ANR-12-BS03-003. 2 the literature and has been used intensi v ely for SU. The unmixing is generally performed using two distinct steps: (i) identifying the endmembers using an endmember extraction algorithm (EEA) such as verte x component analysis (VCA) [7], pixel purity index (PPI) [8] and N-FINDR [9], (ii) estimating the abundances under ph ysical non-ne gati vity and sum-to-one constraints using algorithms such as the fully constrained least squares [2]. Some algorithms also tackle the SU problem in an unsupervised manner , i.e., by jointly estimating the endmembers and the abundances. This is generally achie ved under a statistical frame work using optimization techniques [10] or Markov chain Monte Carlo (MCMC) simulation methods [6], [11]. The unsupervised algorithms generally provide more sophisticated results and appear to be less sensiti v e to the absence of pure pixels [3]. The previous described algorithms provide one endmember spectrum for each physical component present in the image (see Fig. 1(a)). This appears as a clear simplification since in many cases, the endmember spectra v ary along the image causing what is known as spectral variability . Spectral v ariability has been identified as one of the most profound sources of error in abundance estimation and is kno wing growing interest in the hyperspectral community [12], [13]. Many algorithms have been proposed in the literature to describe this variability by considering each endmember as a finite set or as a statistical distribution. Some deterministic approaches represent each physical material as a set or bundle of spectra (see Fig. 1(b)). One can distinguish between algorithms assuming a known spectral library [14], [15] and those estimating it from the data [16], [17]. SU resulting from these approaches is generally sensitiv e to the quality of the a v ailable or e xtracted endmember libraries. There are also statistical approaches assuming that each endmember is a random vector with a giv en distribution (see Fig. 1(c)). Statistical approaches provide a parametric representation of the endmembers and thus can estimate endmembers that are not present in the observed data. This property makes these algorithms more robust in absence of pure pixels [18]–[20]. A more detailed discussion about these algorithms, their adv antages and challenges is a v ailable in [12], [13]. The main contrib ution of this paper is the consideration of endmember variability under a Bayesian frame work. Any endmember is considered via a probability distribution to model its variability . T wo main approaches have been considered in the literature assuming the endmembers are random vectors: the Beta compositional model [19] and the normal compositional model (NCM) [18], [20], [21]. This paper considers a generalization of the NCM model characterized by Gaussian variability for the endmembers (as for the NCM) and an additiv e Gaussian noise modeling fitting errors (which was not present in the NCM). Moreover , the proposed model considers a dif ferent mean and covariance matrix for each endmember to analyze each component separately . These parameters are both estimated to generalize the works of [18] and [20] that estimated the endmembers means and covariances, respecti vely . Moreov er , the endmember fluctuation with respect to the spectral bands is quantified by considering non identically distributed endmember v ariances. 3 (a) (b) (c) Fig. 1: Simplex representation for (a) endmembers without v ariability , (b) endmembers as a finite set (or bundle) and (c) endmembers as a distribution. Another important point concerning hyperspectral unmixing is the spatial correlation between pixels. Indeed, ev en if many algorithms consider a pixel-by-pixel context, recent studies have shown the interest of considering spatial information to improve the unmixing quality [22]–[24]. W ithin a Bayesian framew ork, this spatial correlation can be introduced using Markov random fields (MRFs) as already shown in [22], [23], [25]. In this work, a Potts model is considered since it has already shown good performance when processing hyperspectral images [22], [23]. The image is then segmented into regions sharing similar abundance characteristics. Note that this segmentation was also achieved in [23] and [10] by considering Gaussian and Dirichlet distributions for the abundances. This paper proposes an unsupervised Bayesian algorithm to estimate the parameters associated with endmembers and abundances. In addition to the abundance Dirichlet priors, it assumes appropriate prior for the remaining parameters/hyperparameters to satisfy the kno wn physical constraints. The joint posterior distribution of the proposed Bayesian model is then deriv ed. Howe ver , the classical minimum mean square error (MMSE) and maximum a posteriori (MAP) estimators cannot be easily computed from this joint posterior . A classical way of alle viating this problem is to generate samples distributed according to the posterior using MCMC methods. This goal is achiev ed in this paper using a Gibbs sampler coupled with a Hamiltonian Monte Carlo (HMC) method. HMC is well adapted for problems with a large number of parameters to be estimated [26]. Moreover , this method presents good mixing properties when compared to the classical Metropolis-Hasting algorithm. This paper considers a constrained-HMC (CHMC) that has been introduced in [26, Chap. 5] and successfully used for hyperspectral SU in [11]. This CHMC accounts for inequality constraints which is required to satisfy the physical constraints related to the proposed SU problem. 4 The paper is structured as follows. The unmixing problem considered in this study is formulated in Section II. The different components of the proposed Bayesian model are studied in Section III. Section IV introduces the Gibbs sampler and the CHMC method which will be used to generate samples asymptotically distributed according to the joint posterior of the unknown parameters and hyperparameters. Section V analyzes the performance of the proposed algorithm when applied to syn- thetic images. Results on real hyperspectral images are presented in Section VI whereas conclusions and future works are reported in Section VII. I I . P RO B L E M F O R M U L A T I O N A. Notations The variables used in this paper are described as follows N number of pixels R number of endmembers L number of spectral bands K number of spatial classes Y ∈ R L × N spectra of the pixels A ∈ R R × N abundance matrix T ∈ R R − 1 × N reparameterized abundance matrix M ∈ R L × R endmember means Σ ∈ R R × L matrix containing the diagonal of endmember covariances Ψ ∈ R 1 × N noise variances K ∈ R L × N matrix whose rows equals Ψ z ∈ R 1 × N labels C ∈ R R × K Dirichlet parameters B. Mixing model and endmember variability This section introduces the proposed mixture model. The classical LMM assumes the pixel spectrum y n , n ∈ { 1 , · · · , N } , where N is the number of pixels in the image, is a linear combination of R deterministic endmembers s r , r ∈ { 1 , · · · , R } , corrupted by an additiv e noise as follows y n = R X r =1 a rn s r + e n = S a n + e n (1) with e n ∼ N  0 L , ψ 2 n I L  (2) where R is the number of endmembers, y n is an ( L × 1) vector representing the n th observed pixel, L is the number of spectral bands, 0 L is an ( L × 1 ) vector of 0 , I L is the ( L × L ) identity matrix, 5 a n = [ a 1 n , · · · , a Rn ] T is the ( R × 1) abundance vector of the n th pixel, S = [ s 1 , · · · , s R ] is an ( L × R ) matrix of endmembers and e n is a centered additiv e, independent and identically distributed Gaussian noise. The endmembers are generally variable in the observed image due to en vironmental conditions or inherent variability [13]. In this paper , we introduce a model taking into account this variability . The proposed model can be seen as a generalization of the NCM model (GNCM) since it introduces an additional residual Gaussian noise e as follows y n = R X r =1 a rn s rn + e n = S n a n + e n (3) with s rn ∼ N  m r , diag  σ 2 r  (4) where S n = [ s 1 n , · · · , s Rn ] , σ 2 r =  σ 2 r 1 , · · · , σ 2 rL  is the variance vector of the r th endmember and M = [ m 1 , · · · , m R ] is the ( L × R ) matrix containing the endmember means of the image. The main difference between model (3) and the LMM used in [6] is that the endmember matrix S n depends on each observed pixel in order to introduce the spectral variability . Each physical element is then represented in a giv en pixel by an endmember s rn that has its own Gaussian distribution whose v ariances σ 2 r change from one band to another . This allows the GNCM to capture the spectral v ariations of each physical element with respect to each spectral band. The GNCM also includes an additional Gaussian noise e n ∼ N  0 L , ψ 2 n I L  (that is independent from the v ariables s 1 n , · · · , s Rn ) whose goal is to make the proposed model more robust with respect to mismodeling. Moreover , we consider that the endmember v ariability is the main source of randomness in the observed pixel, which is ensured by assigning a very sparse prior to the noise variance (see Eq. (19)). Note finally that the proposed model reduces to the NCM for ψ 2 n = 0 , ∀ n . Thus, it generalizes the model of [18], [20] by considering a non-isotropic covariance matrix for each endmember . C. Abundance r eparametrization Since the abundance vector a n usually represents spatial coverage of the material in a giv en pixel, it should satisfy the physical positi vity and sum-to-one constraints a rn ≥ 0 , ∀ r ∈ { 1 , . . . , R } and R X r =1 a rn = 1 . (5) Ho we ver , in order to transform the sum to one constraint into an inequality constraint (which will be handled more easily in the algorithm), we propose the following reparametrization a rn = r − 1 Y k =1 t kn ! ×    1 − t rn , if r < R 1 , if r = R . (6) 6 The transformation (6) has been introduced in [27] and has sho wn interesting properties for hyperspec- tral unmixing in [11]. Its main advantage is to express the positivity and the sum to one constraints for the abundances as follows 0 < t rn < 1 , ∀ r ∈ 1 , · · · , R − 1 (7) which will be easily handled in the sampling procedure developed in this paper (see Sections III and IV). I I I . H I E R A R C H I C A L B A Y E S I A N M O D E L This section introduces a hierarchical Bayesian model for unsupervised hyperspectral SU accounting for spectral v ariability . The unkno wn parameters of this model include the ( L × R ) endmember mean matrix M , the ( R × L ) matrix Σ gathering the endmember v ariances (with Σ r,l = σ 2 rl ), the ( R − 1) × N reparameterized abundance matrix T (whose n th column is T : n = t n ), the ( 1 × N ) label vector z and the ( 1 × N ) vector Ψ containing the noise v ariances (with Ψ n = ψ 2 n ). A. Likelihood Using the observation model (3), the Gaussian properties of both the noise sequence e n and the endmembers, and exploiting independence between the observations in dif ferent spectral bands, yield f ( y n | T , M , Σ , z , Ψ ) ∝ 1 Q L l =1 Ω ln ! 1 2 exp  − 1 2 Λ T : n [( y n − M a n )  ( y n − M a n )]  (8) where Ω = Σ T ( A  A ) + K is an ( L × N ) matrix, A = [ a 1 , · · · , a N ] is an ( R × N ) abundance matrix, K = 1 L ⊗ Ψ is an ( L × N ) matrix whose ro ws are equal to Ψ , 1 L is an ( L × 1 ) vector of 1 , Λ is an ( L × N ) matrix with Λ ln = 1 Ω ln ,  denotes the Hadamard (termwise) product and ⊗ denotes the Kroneck er product. Note that the abundance v ector a n ( t n ) has been denoted as a n in (8) for bre vity . Moreo v er , contrary to the LMM, Eq. (8) shows that the elements 1 of Ω depend jointly on the pixel ab undances and on the pixel index # n . This property was also satisfied by the NCM model as previously shown in [18], [20]. Note finally that the joint likelihood of the observation matrix Y can be obtained by exploiting independence between the observed pixels f ( Y | T , M , Σ , z , Ψ ) ∝ N Y n =1 f ( y n | T , M , Σ , z , Ψ ) . (9) B. P arameter priors This section introduces the prior distributions that we have chosen for the parameters of interest z , T (or A ), M , and Σ . 1 The matrix Ω gathers the noise and endmember v ariances. 7 1) Classification prior modeling: Many recent works related to hyperspectral imaging have been considering spatial correlation between the image pixels to segment the image into homogeneous regions with similar abundances [10], [23]. In this paper, we propose to exploit this correlation by di viding the observed image into K classes sharing the same abundance properties [23]. Each pixel is assigned to a specific class by using a latent label variable z n that takes its value into a finite set { 1 , · · · , K } . The whole set of random variables { z n } n =1 , ··· ,N forms a random field. The correlation between neighboring pixels is then introduced by considering a Markov random field prior for z n as follo ws f  z n | z \ n  = f  z n | z ν ( n )  (10) where ν ( n ) denotes the pixel neighborhood as in [23] (a four neighborhood structure will be consid- ered in the rest of the paper), z ν ( n ) = { z i , i ∈ ν ( n ) } and z \ n = { z i , i 6 = n } . As in [22], [23], [28], this paper considers a Potts-Marko v model which is appropriate for hyperspectral image segmentation. The prior of z is then obtained using the Hammersley-Clif ford theorem f ( z ) = 1 G ( β ) exp   β N X n =1 X n 0 ∈ ν ( n ) δ ( z n − z n 0 )   (11) where β > 0 is the granularity coefficient, G ( β ) is a normalizing (or partition) constant and δ ( . ) is the Dirac delta function. The parameter β controls the degree of homogeneity of each region in the image. It is assumed known a priori in this paper . Howe ver , it could be also included within the Bayesian model and estimated using the strategy described in [29]. 2) Abundance matrix T : In order to satisfy the constraints (5), the abundance vector should li ve in the following simplex S S = ( a n   a rn ≥ 0 , ∀ r and R X r =1 a rn = 1 ) . (12) Thus, a natural choice for the prior of a n is a uniform distribution on S [5], [30]. Howe ver , we want to define a prior enforcing strong correlations for close pixels. Therefore, we propose to assign a Dirichlet prior to the abundances of the k th class of the image with Dirichlet parameters c k = ( c 1 k , · · · , c Rk ) T as follows a n | z n = k , c k ∼ Dir ( c k ) , for n ∈ I k (13) where Dir ( . ) denotes the Dirichlet distribution, and n ∈ I k means that y n belongs to the k th class (which is also equiv alent to z n = k ). This prior allows the data to be located in se veral dif ferent clusters inside the simplex [10]. Note that assigning a Dirichlet prior for a n corresponds to a beta distribution prior for the coefficient t rn as shown in [11], [27] t rn | z n = k , C r : R ,k ∼ B e R X i = r +1 c ik , c rk ! , for n ∈ I k (14) 8 where C = [ c 1 , · · · , c K ] is an R × K matrix containing the Dirichlet parameters. The prior associated with the vector t n is finally obtained by assuming prior independence between its elements leading to f ( t n | z n = k , c k ) = Γ  P R i =1 c ik  Q R i =1 Γ ( c ik ) 1 [0 , 1] R − 1 ( t n ) R − 1 Y r =1 t P R i = r +1 c ik − 1 rn (1 − t rn ) c rk − 1 (15) for n ∈ I k , where 1 [0 , 1] R − 1 ( . ) is the indicator of the set [0 , 1] R − 1 . 3) Endmember means: The endmember mean matrix M contains reflectances that should satisfy the following constraints [11] 0 < m rl < 1 , ∀ r ∈ { 1 , · · · , R } , ∀ l ∈ { 1 , · · · , L } . (16) Moreov er , it makes sense to assume that the reflectances are close to estimates computed using an EEA. Therefore, we choose a truncated Gaussian prior for each endmember as follows [11], [20] m r ∼ N [0 , 1] L  f m r ,  2 I l  (17) where f m r denotes an estimated endmember (resulting from an EEA such as VCA 2 ) and  2 is a v ariance term defining the confidence that we have on this estimated endmember f m r . 4) Endmember variances: The absence of kno wledge about the endmember variances can be considered by choosing a Jeffre ys distribution for the parameters σ 2 rl , i.e., f ( Σ : l ) ∝ R Y r =1 1 σ 2 rl 1 R +  σ 2 rl  (18) where we hav e assumed prior independence between the endmember variances. 5) Noise variance prior: T o av oid identifiability problems, the noise effect should be smaller than the effect of endmember variability . This can be achie ved by choosing an exponential prior f  ψ 2 n | λ  = λ exp  − λψ 2 n  1 R +  ψ 2 n  (19) where λ is a large coef ficient imposing sparsity for ψ n ( λ = 10 7 in our simulations). W e furthermore assume prior independence between the random v ariables ψ 2 n , ∀ n ∈ { 1 , · · · , N } . Note that the estimation of ψ 2 n can be remo ved from the proposed Bayesian algorithm without changing significantly the estimation performance (see Section V -D). This paper presents a general formulation allowing the noise effect to be remov ed by setting to zero the noise variance. 2 W e consider in this paper the VCA algorithm e ven if other algorithms such as N-FINDR [9] and pixel purity index (PPI) [8] could also be in vestig ated. 9 C. Hyperparameter priors 1) Dirichlet parameters: The Dirichlet parameters c k are assigned the following conjugate prior [31] f ( c k | z n = k ) =   Γ  P R r =1 c rk  Q R r =1 Γ ( c rk )   γ exp − α R X r =1 c rk + Rα ! R Y r =1 1 R + ( c rk ) (20) where α and γ are fixed constants that have been chosen to ensure a non-informativ e prior (flat distribution). D. P osterior distribution The parameters of the proposed Bayesian model are included in the vector θ = { θ p , θ h } where θ p = { T , M , Σ , z , Ψ } (parameters) and θ h = { C } (hyperparameters). This Bayesian model is summarized in the directed acyclic graph (D A G) displayed in Fig. 2. The joint posterior distribution of the unknown parameter/h yperparameter vector θ can be computed from the following hierarchical structure f ( θ p , θ h | Y ) ∝ f ( Y | θ p , θ h ) f ( θ p , θ h ) (21) where f ( Y | θ p , θ h ) = f ( Y | θ p ) has been defined in (9) and f ( θ p , θ h ) is the joint prior of the unkno wn parameters. Assuming prior independence between the parameters yields f ( θ p , θ h ) = f ( θ p | θ h ) f ( θ h ) = f ( T | C ) f ( M ) f ( Σ ) f ( z ) f ( Ψ ) f ( C ) . (22) The joint posterior distribution f ( θ p , θ h | Y ) can be computed up to a multiplicativ e constant after replacing (9) and (22) in (21). Unfortunately , it is difficult to obtain closed form expressions for the standard Bayesian estimators associated with (21). In this paper , we propose to use MCMC methods to generate samples asymptotically distributed according to (21) and to build estimators of θ from these generated samples. Due to the large number of parameters to be sampled, we use an HMC algorithm which improves the mixing properties of the sampler and reduces the required number of iterations to approximate the target distribution [26]. The parameters are finally estimated using the minimum mean square error (MMSE) estimator for { T , M , Σ , Ψ , C } and the maximum a posteriori (MAP) estimator for the labels z . The next section defines the proposed sampling procedure based on a hybrid Gibbs sampler including a CHMC method. I V . H Y B R I D G I B B S A L G O R I T H M The principle of the Gibbs sampler is to generate samples according to the conditional distrib utions of the target distribution (here the posterior (21)) [32]. When a conditional distribution cannot be 10 Y A Σ z Ψ C β λ M α γ  2 f M Fig. 2: D A G for the parameter and hyperparameter priors (the fixed parameters appear in boxes). Note that the dashed box defines the statistical distribution of the endmember matrix S . sampled directly , sampling techniques such as the Metropolis-Hasting (MH) algorithm can be applied. In this paper , we consider HMC as the proposal strategy since it provides better mixing property than independent or random walk MH moves especially for high-dimensional problems. The next section describes the CHMC algorithm followed by the description of the sampling procedure of the conditional distributions. A. Constrained Hamiltonian Monte Carlo method HMC is used to sample the high dimensional parameter vector of the proposed Bayesian model. It exploits the gradient of the target distribution to improve the quality of the generated samples. Denoting as f ( q ) (resp. q ) the distribution (resp. d-dimensional variable) to be sampled from, HMC defines the Hamiltonian function after introducing a Gaussian momentum variable p as follows H ( p , q ) = U ( q ) + K ( p ) (23) where U ( q ) = − log [ f ( q )] is the potential ener gy related to the tar get distrib ution f ( q ) and K ( p ) = 1 2 p T p is the momentum energy which results from an independent centered Gaussian distribution for p [11]. The ev olution of the ( q , p ) samples is determined using the partial deriv ati ves of the Hamiltonian referred to as Hamiltonian equations [26]. For computer implementations, these equations should be discretized which can be done using the leapfrog method that ensures volume preserv ation and rev ersibility of the chains. This leapfrog discretization scheme moves the samples by an  stepsize, i.e., from the n th state ( q n , p n ) to the ( n + 1) th state  q ( n +1) , p ( n +1)  using N L iteration steps defined 11 by p ( i,n +1 / 2) = p ( i,n ) −  2 ∂ U ∂ q T h q ( i,n ) i (24) q ( i,n +1) = q ( i,n ) +  p ( i,n +1 / 2) (25) p ( i,n +1) = p ( i,n +1 / 2) −  2 ∂ U ∂ q T h q ( i,n +1) i . (26) The resulting samples are accepted with probability ρ giv en by ρ = min n 1 , exp h H ( q n , p n ) − H  q ( n +1) , p ( n +1) io . (27) This procedure ensures the resulting samples to be asymptotically distributed according to the target distribution. In the presence of inequality constraints ( q ( i,n ) ∈ [ q l , q u ] ), we adopt the procedure presented in [11] and [26, Chap. 5]. This procedure replaces a sample that violates the constraints at each leapfrog iteration by its symmetric to the bound (see [11] for more details). For example, the candidate q ( i,n ) = q u + h with 0 < h < ( q u − q l ) will be replaced by q ( i,n ) = q u − h (and similarly q ( i,n ) = q l − h will be replaced by q ( i,n ) = q l + h ) when a constraint is not satisfied. B. Sampling the abundance matrix T It can be shown that the N vectors t n , n ∈ { 1 , · · · , N } are a posteriori independent leading to f ( T | Y , M , Σ , C ) = K Y k =1 Y n ∈I k f ( t n | z n = k , y n , M , Σ , c k ) . (28) Moreov er , using the likelihood (8) and the prior (15) leads to the following conditional distribution f ( t n | z n = k , y n , M , Σ , c k ) ∝ 1 Q L l =1 Ω ln ! 1 2 exp  − 1 2 Λ T : n [( y n − M a n )  ( y n − M a n )]  × 1 [0 , 1] R − 1 ( t n ) R − 1 Y r =1 t P R i = r +1 c ik − 1 rn (1 − t rn ) c rk − 1 (29) for n ∈ I k . The conditional distribution (29) is not easy to sample. Ho we ver , the CHMC frame work is well suited for sampling the independent vectors t n , n ∈ { 1 , · · · , N } in an effecti v e parallel procedure that reduces the computational cost. Moreov er , the small size of these vectors (of size ( R − 1) × 1 ) improv es the con ver gence of the sampler . Note that the CHMC requires the definition of the potential energy U ( t n ) = − log [ f ( t n | z n = k , y n , M , Σ , c k )] giv en by U ( t n ) = U 1 + U 2 + U 3 (30) 12 with U 1 = 1 2 Λ T : n [( y n − M a n )  ( y n − M a n )] U 2 = − R X r =1 ( R X i = r +1 c ik − 1 ! log ( t rn ) + ( c rk − 1) log (1 − t rn ) ) U 3 = 1 2 L X l =1 log ( Ω ln ) . (31) Note finally that the deriv ati v es of U with respect to the variable of interest t n (that are required for the CHMC steps) are provided in the appendix. C. Sampling the mean endmember matrix M Straightforward computations using the posterior distribution (21) yield f ( M | Y , T , Σ ) = L Y l =1 f ( M l : | Y l : , T , Σ : l ) (32) where f ( M l : | Y l : , T , Σ : l ) ∝ exp  − 1 2 [( Y l : − M l : A )  ( Y l : − M l : A )] Λ T l :  × exp − || M l : − f M l : || 2 2  2 ! 1 [0 , 1] R ( M l : ) . (33) Equation (32) results from the independence between the columns of the matrix M (vectors of small size R × 1 ). This interesting property promotes the use of a parallel CHMC algorithm for sampling T . The potential energy V associated with the conditional distribution of M l : is giv en by V ( M l : ) = 1 2 [( Y l : − M l : A )  ( Y l : − M l : A )] Λ T l : + || M l : − f M l : || 2 2  2 . (34) The deriv ati v es of V with respect to M l : are provided in the appendix. D. Sampling the variance of the endmember matrix Considering (21) yields the following conditional distribution for matrix Σ containing the end- member variances f ( Σ | Y , T , M ) = L Y l =1 f ( Σ : l | Y l : , T , M l : ) (35) with f ( Σ : l | Y l : , T , M l : ) ∝ 1 Q N n =1 Ω ln ! 1 2 exp  − 1 2 [( Y l : − M l : A )  ( Y l : − M l : A )] Λ T l :  × R Y r =1 1 σ 2 rl 1 R +  σ 2 rl  . (36) 13 Sampling from (36) can again be performed using a CHMC algorithm (as in Sections IV -B and IV -C). The potential energy associated with the vector Σ : l is W ( Σ : l ) = W 1 + W 2 + W 3 (37) with W 1 = 1 2 [( Y l : − M l : A )  ( Y l : − M l : A )] Λ T l : W 2 = R X r =1 log  σ 2 rl  W 3 = 1 2 N X n =1 log ( Ω ln ) . (38) The deriv ati v es of W with respect to Σ : l are provided in the appendix. E. Sampling the labels The conditional distribution associated with the discrete random variable z n is giv en by f ( z n = k | t n , c k ) ∝ f ( t n | z n = k , c k ) exp   β X n 0 ∈ ν ( n ) δ ( k − z n 0 )   (39) where f ( t n | z n = k , c k ) has been defined in (15). Sampling from this conditional distribution is classically performed by drawing a discrete value in the finite set { 1 , · · · , K } with the probabilities (39). F . Sampling the noise variance Ψ Considering (21) yields the following conditional distribution for the noise variance matrix Ψ f ( Ψ | z , T , Y , M , Σ , c ) = N Y n =1 f  ψ 2 n | z n = k , t n , y n , M , Σ , c k  (40) with f  ψ 2 n | z n = k , t n , y n , M , Σ : l , c k  ∝ 1 Q L l =1 Ω ln ! 1 2 exp  − 1 2 Λ T : n [( y n − M a n )  ( y n − M a n )]  × exp  − λψ 2 n  1 R +  ψ 2 n  (41) for n ∈ I k . This distribution is sampled using a parallel CHMC procedure with the following potential energy H  ψ 2 n  = U 1 + U 3 + λψ 2 n . (42) 14 G. Sampling the Dirichlet coefficients Using (21) and (22), it can be easily sho wn that the conditional distrib ution of c k | T , z n ∈I k is given by f ( c k | T , z n ∈I k ) ∝ Y n ∈I k        Γ  P R r =1 c rk  Q R r =1 Γ ( c rk )   γ +1 exp − α R X r =1 c rk + Rα ! R Y r =1 a c rk − 1 rn      (43) for k ∈ { 1 , · · · , K } . This distribution is also sampled using a CHMC procedure. The corresponding potential energy is giv en by P ( c k ) = P 1 + P 2 (44) with P 1 = ( γ + 1) X n ∈I k " − log Γ R X r =1 c rk ! + R X r =1 log Γ ( c rk ) # P 2 = X n ∈I k " α R X r =1 c rk − Rα − R X r =1 log  a c rk − 1 rn  # . (45) V . S I M U L A T I O N R E S U LT S O N S Y N T H E T I C DA T A This section ev aluates the performance of the proposed algorithm with synthetic data. It is di vided into four parts whose objecti ves are: 1) introduce the criteria used for the ev aluation of the unmixing quality , 2) present the different parameters that are estimated in the proposed unmixing approach, 3) analyze the behavior of the proposed algorithm as a function of the number of endmembers and the size of the image, 4) compare the proposed strategy with other state-of-the-art algorithms from the literature. A. Evaluation criteria Abundances and endmembers are known for synthetic images. In this case, the quality of the unmixing strategy can be measured by comparing the estimated and actual abundances by using the av erage root mean square error (aRMSE) defined by aRMSE ( A ) = v u u t 1 N R N X n =1 k a n − ˆ a n k 2 (46) where || · || denotes the standard l 2 norm such that || x || 2 = x T x . The mean of the r th estimated endmember can be compared with the actual one by using RMSE ( m r ) or the spectral angle mapper SAM ( m r ) defined as follows RMSE ( m r ) = 1 √ L k ˆ m r − m r k , SAM ( m r ) = arccos ˆ m T r m r k m r k k ˆ m r k ! (47) 15 where arccos( · ) is the in verse cosine operator . Moreover , the global endmember error is ev aluated by the av eraged RMSE (aRMSE) and av eraged SAM (aSAM) given by aRMSE ( M ) = v u u t 1 R R X r =1 [ RMSE ( m r )] 2 , aSAM ( M ) = 1 R R X r =1 SAM ( m r ) . (48) Note finally that the RE and SAM criteria can also be e v aluated for the # p th measured and estimated pixel spectra y n , ˆ y n as follows RE = v u u t 1 N L N X n =1 k ˆ y n − y n k 2 , SAM = 1 N N X n =1 arccos ˆ y T n y n k y n k k ˆ y n k ! . (49) B. P erformance of the pr oposed algorithm This section considers a 50 × 50 synthetic image generated according to (3) with R = 3 endmembers (construction concrete, green grass and micaceous loam) that hav e been extracted from the ENVI software library [33]. The considered endmember v ariances depend on the spectral bands as sho wn in Fig. 3 (dashed lines). This image contains K = 3 classes whose label maps hav e been generated Fig. 3: Actual endmember v ariances (dassed line) and estimated v ariances by the proposed UsGNCM (continuous line) for the considered R = 3 endmembers. using (11) with β = 1 . 5 (see Fig. 4). The ab undances of each class share the same Dirichlet parameters (that are reported in T able I) leading to the observed pixels displayed in Fig. 5. Note that the generated abundances hav e been truncated ( a r < 0 . 9 , ∀ r ) to av oid the presence of pure pixels in the image. Finally , we have considered a noise v ariance equal to 10 − 7 (note that the noise v ariance has to be 16 smaller than the endmember v ariances). The proposed unsupervised GNCM-based algorithm, denoted by UsGNCM, has been run using N bi = 11000 burn-in iterations and N MC = 12000 iterations. Fig. T ABLE I: Actual and estimated Dirichlet parameters in each spatial class. Dirichlet parameters c 1 k c 2 k c 3 k ˆ c 1 k ˆ c 2 k ˆ c 3 k k = 1 15 15 1 14.97 14.85 1.00 k = 2 1 8 8 1.05 8.24 8.19 k = 3 3 1 3 3.12 1.02 3.03 4 (right) displays the estimated classification map obtained with the proposed algorithm. This map is in a very good agreement with the ground truth shown in Fig. 4 (left). Note that the Dirichlet parameters used in this simulation correspond to three distinguishable classes that are well separated using the proposed algorithm. The obtained classification results can also be observed with the data projected in the plane associated with the two most discriminant principle components as shown in Fig. 5. The proposed algorithm also allows the Dirichlet parameters to be estimated accurately as sho wn in T able I. A significant advantage of the proposed algorithm is its ability to estimate the endmember means and variances. Fig. 5 shows the estimated endmembers obtained using the VCA algorithm (diamonds) [7], the UsLMM algorithm (circles) [6] and the proposed UsGNCM approach (triangles). Contrary to the VCA algorithm that provides bad endmember estimates because of the absence of pure pixels in the image, both UsLMM and UsGNCM strategies yield good endmember estimations. As explained before, the good performance of the UsGNCM algorithm can be explained by the fact that it is able mitigate the endmember v ariability . Fig. 6 displays the endmember means (continuous lines), the endmember distributions (colored areas in Figs. 6(a), (b) and (c)) and the associated v ariability interv als defined by mean ± 3 σ (Fig. 6 (d)). Fig. 3 displays the actual and estimated endmember v ariances for the three endmembers that are clearly in good agreement. These results show the good performance of the proposed approach that fully exploits the spatial (segmentation map, abundances and noise v ariances) and spectral (endmember means and variances) correlations. The ne xt section studies the robustness of the proposed approach with respect to the number of endmembers and the image size (number of pixels). C. P erformance as a function of the number of endmembers and the image size The UsGNCM algorithm estimates many parameters which might require a lot of observations in order to obtain acceptable performance. The first part of this section deals with this problem by 17 Fig. 4: Actual (left) and estimated (right) classification maps of a synthetic image. Fig. 5: Classified projected pixels (colored crosses), actual endmembers (red stars), endmembers esti- mated by VCA (black diamonds), endmembers estimated by UsLMM (cyan circle) and endmembers estimated by UsGNCM (blue triangles). analyzing the proposed algorithm when varying the number of observed pixels. The considered image has been generated using the three endmembers considered in Section V -B with abundances uniformly distributed in the simplex S defined by the positivity and sum-to-one constraints (the corresponding Dirichlet parameters are c rk = 1 , ∀ r , ∀ k ). Fig. 7 shows the obtained aRMSE ( A ) , RE and SAM when varying the size of the observed image. As expected, the unmixing performance improves by 18 (a) (b) (c) (d) Fig. 6: Actual endmembers (crosses) and endmember means estimated by UsGNCM (continuous lines). The estimated endmember distributions are represented in (a), (b), (c) by colored areas. The bottom-right figure (d) shows the endmembers estimated by UsGNCM ± 3 σ (dashed lines). increasing the number of observations. This figure also shows that the aRMSE ( A ) con ver ges to a constant value for √ N > 50 while RE and SAM continue to improv e when increasing N . Note, ho we ver , that the obtained results are quite good for N ≥ 100 . The second part of this section analyzes the behavior of UsGNCM with respect to the number of endmembers. T able II shows the obtained aRMSE ( A ) , aRMSE ( M ) and aSAM ( M ) criteria for R = { 3 , 4 , 5 , 6 } . The considered endmembers are construction concrete, green grass, micaceous loam, oli v e green paint, bare red brick, and galv anized steel metal. These spectra have been extracted from the spectral libraries provided with the ENVI software [33]. As previously , the images associated with R = { 3 , 4 , 5 , 6 } hav e been generated with abundances uniformly distrib uted in the simple x S . As expected, increasing the number of parameters (i.e., increasing R ) reduces the estimation performance. Howe v er , the obtained results 19 are still acceptable confirming the rob ustness of UsGNCM with respect to the number of endmembers R . Fig. 7: UsGNCM performance for different numbers of pixels. T ABLE II: UsGNCM performance for different number of endmembers. aRMSE ( A ) aRMSE ( M ) aSAM ( M ) ( × 10 − 2 ) ( × 10 − 2 ) ( × 10 − 2 ) R = 3 0.5734 0.1826 0.4306 R = 4 0.542 0.2115 0.5073 R = 5 0.8053 0.2790 0.6837 R = 6 1.4049 0.8404 1.6550 D. Comparison with state-of-the-art algorithms This section ev aluates the performance of the proposed UsGNCM algorithm for different images. All images ha v e been constructed using R = 3 endmembers with truncated ab undances (with a i < 0 . 9 , ∀ i ∈ 1 , · · · , R ) to av oid the presence of pure pixels. The remaining parameters hav e been defined as follo ws 20 • the image I 1 has been generated according to the GNCM model with K = 1 class and ab undances uniformly distributed in the simplex S . The endmember variances ha ve been adjusted as in Fig. 3. The noise variance is ψ 2 n = 10 − 7 . • the image I 2 is the GNCM image used in Section V -B. • the image I 3 has been generated according to the LMM model with K = 3 classes and the Dirichlet parameters of T able I. The noise v ariances v ary linearly with respect to the spectral bands with ψ 2 l = 10 − 4  4 L − 1 l + L + 3 L − 1  , for l ∈ [1 , · · · , L ] . These images are processed using dif ferent unmixing strategies that are compared to the proposed UsGNCM algorithm. More precisely , we hav e considered the following unmixing algorithms • VCA+FCLS: the endmembers are extracted from the whole image using VCA and the ab undances are estimated using the FCLS algorithm [2]. • UsLMM: the unsupervised Bayesian algorithm of [6] is used to jointly estimate the endmembers and abundances. • AEB: this is the automated endmember bundles algorithm proposed in [17]. W e consider a 10% image subset and the VCA algorithm to extract the endmembers. For each pixel, the 3 endmembers that provide the smallest RE are selected. • UsNCM: the proposed unmixing strategy with ψ n = 0 (i.e., the additive noise e n of (3) is remov ed). Note that the resulting algorithm reduces to the NCM model. The first two algorithms provide one estimate for each endmember while the other algorithms estimate endmember variability . Note that the UsNCM is introduced to study the effect of the additive noise. T able III reports the quality of the estimated abundances and endmembers by unmixing the three images with the different algorithms. This table sho ws bad performance for VCA+FCLS and AEB algorithms which is mainly due to the absence of pure pixels in the considered images. The UsLMM provides good results for the three images. Howe ver , it appears to be sensitiv e to the variation of endmember/noise variances with respect to the spectral band and to the spatial correlations between adjacent pixels. Indeed, the UsLMM did not consider spatial correlation which leads to a performance reduction when processing the images I 2 and I 3 . Note also that the UsLMM algorithm provides one estimate for each endmember and does not take into account the spatial variability of endmembers in the processed images. The best performance is generally obtained by the proposed UsNCM and UsGNCM strategies that provide almost similar results. Howe v er , the UsGNCM algorithm is more robust than UsNCM when processing the LMM image I 3 . Moreover , the UsGNCM provides the best endmember estimates as highlighted by the criteria ARE and ASAM. These results confirm the superiority of the proposed approach in presence of endmember variability , spatial correlation 21 between pixels and in absence of pure pixels in the observed scene. T ABLE III: Results on synthetic data. Criteria ( × 10 − 2 ) aRMSE RMSE RMSE RMSE SAM SAM SAM aRMSE aSAM ( A ) ( m 1 ) ( m 2 ) ( m 3 ) ( m 1 ) ( m 2 ) ( m 3 ) ( M ) ( M ) image I 1 VCA+FCLS 4.78 2.29 1.97 2.31 6.14 1.71 5.74 2.20 4.53 UsLMM 0.52 0.25 0.10 0.15 0.77 0.20 0.33 0.18 0.43 (GNCM, AEB 3.73 2.33 1.80 2.55 4.92 2.04 7.75 2.25 4.90 K = 1 ) UsNCM 0.48 0.09 0.10 0.21 0.30 0.20 0.43 0.14 0.31 UsGNCM 0.48 0.07 0.09 0.18 0.25 0.19 0.40 0.12 0.28 image I 2 VCA+FCLS 3.71 2.89 2.98 2.08 9.54 5.60 5.09 2.68 6.74 UsLMM 0.76 0.21 0.40 0.73 0.67 0.86 1.28 0.49 0.94 (GNCM, AEB 9.46 3.48 4.67 4.37 7.96 13.26 4.94 4.20 8.72 K = 3 ) UsNCM 0.56 0.14 0.11 0.28 0.39 0.21 0.69 0.19 0.43 UsGNCM 0.48 0.14 0.11 0.21 0.46 0.26 0.51 0.16 0.41 image I 3 VCA+FCLS 9.51 3.63 6.21 2.61 12.20 7.73 5.60 4.42 8.51 UsLMM 1.01 0.75 0.18 0.34 2.62 0.30 0.73 0.49 1.22 (LMM, AEB 9.30 5.98 4.67 4.61 16.05 5.86 10.84 5.13 10.92 K = 3 ) UsNCM 0.86 0.45 0.44 0.55 1.55 0.91 0.99 0.48 1.15 UsGNCM 0.74 0.20 0.29 0.47 0.70 0.56 0.98 0.34 0.74 V I . S I M U L AT I O N R E S U LT S O N R E A L D A T A A. Description of the Hyperspectral Data This section illustrates the performance of the proposed UsGNCM algorithm when applied to a real hyperspectral data set. The real image used in this section was acquired in 2010 by the Hyspex hyperspectral scanner over V illelongue, France (00 03’W and 4257’N). The dataset contains L = 160 spectral bands recorded from the visible to near infrared with a spatial resolution of 0 . 5 m [34]. The proposed unmixing algorithm has been applied to two subimages: scene #1 of size 50 × 50 which is composed of R = 4 components: tree, grass, soil and shadow (see Fig. 8 (right)), and scene #2 of size 31 × 31 which is composed of R = 3 components: grass, road and ditch (see Fig. 8 (left)). B. Endmember V ariability The proposed UsGNCM algorithm can estimate both the endmember means and variances. Fig. 9 compares the endmember estimates of this algorithm with those obtained with VCA and UsLMM 22 Fig. 8: Real Madonna image and the considered subimages shown in true colors. (Right) scene 1, (left) scene 2 when considering scene #1 . The estimated endmembers are globally in good agreement. Note that VCA pro vides a different shado w endmember because it estimates the endmember as the purest pix el in the image while UsLMM and UsGNCM estimate both the abundances and endmembers resulting in a better shadow estimate (lower amplitude). Moreover , the proposed algorithm provides endmember distributions (blue lev el areas in Fig. 9) which measure the endmember variability in the considered image. It can be seen that the higher relati ve variation is obtained for the shadow spectrum because of its lo w amplitude. Moreover , the v ariation is more pronounced for high spectral bands ( l > 80 ) which is in agreement with the results presented in [11]. Fig. 10 shows the obtained endmembers when considering scene #2 . This figure presents similar results between UsGNCM and UsLMM, especially for capturing spectral components having lo w amplitudes as for ditch. 23 Fig. 9: The R = 4 endmembers estimated by VCA (continuous red lines), UsLMM (continuous black lines), UsGNCM (continuous blue lines) and the estimated endmember distribution (blue lev el areas) for the Madonna image. C. Abundance Estimation and Image Classification The fraction maps of scene #1 estimated by the proposed method are shown in Fig. 11 (bottom). Note that a white (black) pixel indicates a large (small) proportion of the corresponding materials. These pictures are in good agreement with the FCLS and UsLMM results shown in Fig. 11 (top) and (middle), respectiv ely . Note that the compared algorithms also provide similar abundance maps when considering scene #2 . Howe ver , these results are not presented here for brevity . In addition to unmixing, UsGNCM also provides a spatial segmentation of the considered scenes as shown in Fig. 12(a) for scene #1 and Fig. 13(a) for scene #2 . These classifications clearly highlight the area of each physical element in the scene. Indeed, for scene #1 we have 4 classes that represent tree, shado w , soil and grass zones while for scene #2 we hav e 3 classes representing road, ditch and grass areas. T able IV finally reports the estimated Dirichlet parameters and the number of pixels for each spatial class when considering scene #1 . These parameters suggest a highly non uniform distrib ution 24 Fig. 10: The R = 3 endmembers estimated by VCA (continuous red lines), UsLMM (continuous black lines), UsGNCM (continuous blue lines) and the estimated endmember distribution (blue lev el areas) for the Madonna image. ov er the simplex which promote the use of the proposed approach. Fig. 11: Abundance maps estimated by FCLS (top), UsLMM (middle) and the proposed UsGNCM (bottom) for the Madonna image. 25 (a) Classification map. (b) Noise variances. Fig. 12: Estimated maps with the UsGNCM algorithm for the scene #1 of Madonna image. (a) Classification map and (b) noise variances. (a) Classification map. (b) Noise variances. Fig. 13: Estimated maps with the UsGNCM algorithm for the scene #2 of Madonna image. (a) Classification map and (b) noise variances. D. Residual Components The proposed algorithm also provides a measure of the noise variance for each observed pixel. This parameter brings an information about pixels that are inaccurately described by a linear formulation, i.e., allows modeling errors to be quantified. Fig. 12(b) shows the obtained noise variances for the scene #1 . This figure shows a higher error in the shadow area and around trees, i.e., for regions where possible interactions between physical components might occur (e.g., tree/soil) resulting in a more complex model than the proposed linear one. The noise variances associated with the scene #2 are sho wn in Fig. 13(b). This figure shows a higher error near the ditch area which might be 26 T ABLE IV: Estimated Dirichlet parameters for the Madonna image. Dirichlet parameters number of ˆ c 1 k ˆ c 2 k ˆ c 3 k ˆ c 4 k pixels k = 1 7.8767 2.8933 1.0139 5.1277 613 k = 2 2.8914 7.6524 1.3115 1.7289 318 k = 3 12.3176 16.1875 21.2009 21.1454 445 k = 4 25.7654 26.4822 17.0927 49.9141 1124 due to the presence of nonlinearities as explained in [11]. Note finally that both Fig. 12(b) and Fig. 13(b) highlight the presence of regular vertical patterns that hav e also been observed in [35] and were associated with a sensor defect or other miscalibration problems. V I I . C O N C L U S I O N S This paper introduced a Bayesian model for unsupervised unmixing of hyperspectral images ac- counting for spectral variability . The proposed algorithm was based on a generalization of the normal compositional model and includes an additiv e Gaussian noise for modeling errors. This algorithm esti- mated the endmembers of the scene, their variabilities provided by their variances and the correspond- ing ab undances. The observ ed image w as also spatially segmented into re gions sharing homogeneous abundance characteristics. The physical constraints of the abundances were ensured by choosing a Dirichlet distribution for each spatial class of the image. Due to the complexity of the resulting joint posterior distribution, a Mark ov chain Monte Carlo procedure based on a Gibbs algorithm was used to sample the posterior of interest and to approximate the Bayesian estimators of the unkno wn parameters using the generated samples. The sampling was achiev ed using an Hamiltonian Monte Carlo method which is well suited for problems with a lar ge number of parameters. The proposed algorithm showed good performance when processing data presenting endmember variability , spatial correlation between pixels and in absence of pure pixels in the observed scene. UsGNCM fully exploits both the spatial dimension (segmentation, abundance and noise estimation) and the spectral dimension (estimation of endmember means and variances). Future work includes the study of endmember variability with nonlinear mixing models. This point is an interesting issue that is currently under in vestig ation. 27 A P P E N D I X D E R I V A T I V E S O F T H E P OT E N T I A L F U N C T I O N S The deriv ati v e of U with respect to t n is giv en by ∂ U ∂ t n = ∂ U 1 ∂ a n ∂ a n ∂ t n + ∂ U 2 ∂ t n + ∂ U 3 ∂ a n ∂ a n ∂ t n (50) with ∂ U 1 ∂ a n = − [ Λ : n  ( y n − M a n )] T M + 1 2 [( y n − M a n )  ( y n − M a n )] T  ∂ Λ : n ∂ a n  T  ∂ Λ ln ∂ a n  T = − 2 diag ( a n ) Σ : ,l Ω 2 ln ∂ U 3 ∂ a n = a T n  [ ΣΛ : n ] T ∂ U 2 ∂ t rn = − P R i = r +1 c ik − 1 t rn + c rk − 1 1 − t rn , ∀ r ∈ { 1 , · · · , R − 1 } (51) and ∂ a rn ∂ t in =            0 if i > r a rn t in − 1 if i = r a rn t in if i < r . (52) The deriv ati v e of V with respect to M l : is giv en by ∂ V ∂ M l : = − [ Λ l :  ( Y l : − M l : A )] A T + 1  2  M l : − f M l :  . (53) The deriv ati v es of W with respect to Σ : l are giv en by ∂ W 1 ∂ Σ : l = − 1 2  ( Y l : − M l : A )  ( Y l : − M l : A ) Ω l :  Ω l :  ( A  A ) T ∂ W 2 ∂ Σ 2 rl = ∂ W 2 ∂ σ 2 rl = 1 σ 2 rl , ∀ r ∈ { 1 , · · · , R } ∂ W 3 ∂ Σ : l = 1 2 h Λ l : ( A  A ) T i (54) The deriv ati v es of H with respect to ψ 2 n is giv en by ∂ T ∂ ψ 2 n = ∂ U 1 ∂ ψ 2 n + ∂ U 3 ∂ ψ 2 n + λ (55) with ∂ U 1 ∂ ψ 2 n = − 1 2 L X l =1 ( y n − M a n )  ( y n − M a n ) Ω : n  Ω : n ∂ U 3 ∂ ψ 2 n = − 1 2 L X l =1 Λ : n (56) 28 The deriv ati v e of P with respect to c rk is giv en by ∂ P ∂ c rk = ∂ P 1 ∂ c rk + ∂ P 2 ∂ c rk (57) with ∂ P 1 ∂ c rk = ( γ + 1) X n ∈I k " − Υ R X r 0 =1 c r 0 k ! + Υ ( c rk ) # ∂ P 2 ∂ c rk = X n ∈I k [ α − log ( a rn )] (58) where Υ denotes the polygamma function, i.e., the deriv ati ve of the log-gamma function. 29 R E F E R E N C E S [1] N. K esha va and J. F . Mustard, “Spectral unmixing, ” IEEE Signal Pr ocess. Mag. , pp. 44–57, Jan. 2002. [2] D. C. Heinz and C. -I Chang, “Fully constrained least-squares linear spectral mixture analysis method for material quantification in hyperspectral imagery , ” IEEE T rans. Geosci. Remote Sens. , vol. 29, no. 3, pp. 529–545, March 2001. [3] J. Bioucas-Dias, A. Plaza, N. Dobigeon, M. Parente, Q. Du, P . Gader, and J. Chanussot, “Hyperspectral unmixing ov erview: Geometrical, statistical, and sparse regression-based approaches, ” IEEE J. Sel. T opics Appl. Earth Observat. Remote Sens. , vol. 5, no. 2, pp. 354–379, April 2012. [4] N. Dobigeon, J.-Y . T ourneret, C. Richard, J. Bermudez, S. McLaughlin, and A. Hero, “Nonlinear unmixing of hyperspectral images: Models and algorithms, ” IEEE Signal Pr ocess. Mag . , vol. 31, no. 1, pp. 82–94, Jan 2014. [5] A. Halimi, Y . Altmann, N. Dobigeon, and J.-Y . T ourneret, “Nonlinear unmixing of hyperspectral images using a generalized bilinear model, ” IEEE T r ans. Geosci. Remote Sens. , vol. 49, no. 11, pp. 4153–4162, 2011. [6] N. Dobigeon, S. Moussaoui, M. Coulon, J.-Y . T ourneret, and A. O. Hero, “Joint Bayesian endmember extraction and linear unmixing for hyperspectral imagery , ” IEEE T rans. Signal Pr ocess. , v ol. 57, no. 11, pp. 4355–4368, Nov . 2009. [7] J. M. Nascimento and J. M. Bioucas-Dias, “V ertex component analysis: A f ast algorithm to unmix hyperspectral data, ” IEEE Tr ans. Geosci. Remote Sens. , v ol. 43, no. 4, pp. 898–910, April 2005. [8] J. Boardman, “ Automating spectral unmixing of A VIRIS data using con vex geometry concepts, ” in Summaries 4th Annu. JPL Airborne Geoscience W orkshop , vol. 1. W ashington, D.C.: JPL Pub., 1993, pp. 11–14. [9] M. W inter , “F ast autonomous spectral end-member determination in hyperspectral data, ” in Pr oc. 13th Int. Conf. on Applied Geologic Remote Sens. , v ol. 2, V ancouver , April 1999, pp. 337–344. [10] J. Nascimento and J. Bioucas-Dias, “Hyperspectral unmixing based on mixtures of Dirichlet components, ” IEEE T rans. Geosci. Remote Sens. , vol. 50, no. 3, pp. 863–878, March 2012. [11] Y . Altmann, N. Dobigeon, S. McLaughlin, and J.-Y . T ourneret, “Unsupervised post-nonlinear unmixing of h yperspectral images using a Hamiltonian Monte Carlo algorithm, ” IEEE T r ans. Image Pr ocess. , vol. 23, no. 6, pp. 2663–2675, June 2014. [12] B. Somers, G. P . Asner , L. Tits, and P . Coppin, “Endmember variability in spectral mixture analysis: A revie w , ” Remote Sens. Envir on. , vol. 115, no. 7, pp. 1603 – 1616, 2011. [13] A. Zare and K. Ho, “Endmember v ariability in hyperspectral analysis: Addressing spectral v ariability during spectral unmixing, ” IEEE Signal Pr ocess. Mag . , vol. 31, no. 1, pp. 95–104, Jan 2014. [14] D. Roberts, M. Gardner, R. Church, S. Ustin, G. Scheer , and R. Green, “Mapping Chaparral in the Santa Monica Mountains Using Multiple Endmember Spectral Mixture Models, ” Remote Sens. Envir on. , v ol. 65, no. 3, pp. 267–279, Sept. 1998. [15] C. Bateson, G. Asner, and C. W essman, “Endmember bundles: a ne w approach to incorporating endmember variability into spectral mixture analysis, ” IEEE T rans. Geosci. Remote Sens. , vol. 38, no. 2, pp. 1083–1094, Mar 2000. [16] M. Goenaga, M. T orres-Madronero, M. V elez-Reyes, S. V an Bloem, and J. Chinea, “Unmixing analysis of a time series of hyperion images over the Gunica dry forest in Puerto Rico, ” IEEE J. Sel. T opics Appl. Earth Observat. Remote Sens. , vol. 6, no. 2, pp. 329–338, April 2013. [17] B. Somers, M. Zortea, A. Plaza, and G. Asner, “ Automated extraction of image-based endmember b undles for improv ed spectral unmixing, ” IEEE J . Sel. T opics Appl. Earth Observat. Remote Sens. , v ol. 5, no. 2, pp. 396–408, April 2012. [18] O. Eches, N. Dobigeon, C. Mailhes, and J.-Y . T ourneret, “Bayesian estimation of linear mixtures using the normal compositional model. Application to hyperspectral imagery , ” IEEE T rans. Image Pr ocess. , vol. 19, no. 6, pp. 1403– 1413, June 2010. [19] A. Zare, P . Gader, D. Drashnikov , and T . Glenn, “Beta compositional model for hyperspectral unmixing, ” in Proc. IEEE GRSS WHISPERS , Gainesville, USA, Jun. 2013. 30 [20] A. Zare, P . Gader , and G. Casella, “Sampling piecewise con v ex unmixing and endmember extraction, ” IEEE T r ans. Geosci. Remote Sens. , vol. 51, no. 3, pp. 1655–1665, March 2013. [21] D. Stein, “ Application of the normal compositional model to the analysis of hyperspectral imagery , ” in Proc. IEEE W orkshop Advances in T echniques for Analysis of Remotely Sensed Data , Oct 2003, pp. 44–51. [22] N. Bali and A. Mohammad-Djafari, “Bayesian approach with hidden Markov modeling and mean field approximation for hyperspectral data analysis, ” IEEE T rans. Image Pr ocess. , v ol. 17, no. 2, pp. 217–225, Feb . 2008. [23] O. Eches, N. Dobigeon, and J.-Y . T ourneret, “Enhancing hyperspectral image unmixing with spatial correlations, ” IEEE Tr ans. Geosci. Remote Sens. , v ol. 49, no. 11, Nov 2011. [24] J. Chen, C. Richard, and P . Honeine, “Nonlinear estimation of material abundances in hyperspectral images with ` 1 -norm spatial regularization, ” IEEE T rans. Geosci. Remote Sens. , vol. 52, no. 5, pp. 2654–2665, May 2014. [25] R. Rand and D. Keenan, “Spatially smooth partitioning of hyperspectral imagery using spectral/spatial measures of disparity , ” IEEE T rans. Geosci. Remote Sens. , v ol. 41, no. 6, pp. 1479–1490, June 2003. [26] S. Brooks, A. Gelman, G. L. . Jones, and X.-L. Meng, Handbook of Markov Chain Monte Carlo . ser . Chapman & Hall/CRC Handbooks of Modern Statistical Methods. T aylor & Francis, 2011. [27] M. J. Betancourt, “Cruising the simplex: Hamiltonian Monte Carlo and the Dirichlet distribution, ” in ArXiv e-prints , Oct. 2013. [28] Y . Altmann, N. Dobigeon, S. McLaughlin, and J.-Y . T ourneret, “Residual component analysis of hyperspectral images: Application to joint nonlinear unmixing and nonlinearity detection, ” IEEE T rans. Image Process. , vol. 23, no. 5, pp. 2148–2158, May 2014. [29] M. Pereyra, N. Dobigeon, H. Batatia, and J.-Y . T ourneret, “Estimating the granularity coefficient of a Potts-Markov random field within a Markov Chain Monte Carlo algorithm, ” IEEE T rans. Imag e Pr ocess. , vol. 22, no. 6, pp. 2385– 2397, June 2013. [30] N. Dobigeon, J.-Y . T ourneret, and C.-I Chang, “Semi-supervised linear spectral unmixing using a hierarchical Bayesian model for hyperspectral imagery , ” IEEE T rans. Signal Pr ocess. , vol. 56, no. 7, pp. 2684–2695, July 2008. [31] Z. Ma, “Bayesian estimation of the Dirichlet distribution with expectation propagation, ” in Proc. EUSIPCO , Bucharest, Romania, Aug. 2012. [32] C. P . Robert, The Bayesian Choice: fr om Decision-Theor etic Motivations to Computational Implementation , 2nd ed., ser . Springer T e xts in Statistics. New Y ork: Springer-V erlag, 2007. [33] RSI (Research Systems Inc.), ENVI User’ s guide V ersion 4.0 , Boulder , CO 80301 USA, Sept. 2003. [34] D. Sheeren, M. Fauv el, S. Ladet, A. Jacquin, G. Bertoni, and A. Gibon, “Mapping ash tree colonization in an agricultural mountain landscape: In vestigating the potential of hyperspectral imagery , ” in Proc. IEEE Int. Conf. Geosci. Remote Sens. (IGARSS) , July 2011, pp. 3672–3675. [35] C. Févotte and N. Dobigeon, “Nonlinear hyperspectral unmixing with robust nonnegati ve matrix factorization, ” in ArXiv e-prints , March 2014.