A Hierarchical Bayesian Model Accounting for Endmember Variability and Abrupt Spectral Changes to Unmix Multitemporal Hyperspectral Images

1 A Hierarchical Bayesian Model Accounting for Endmember V ariability and Abrupt Spectral Changes to Unmix Multitemporal Hyperspectral Images Pierre-Antoine Thouvenin, Member , IEEE , Nicolas Dobigeon, Senior Member , IEEE and Jean-Yves T ourneret, Senior Member , IEEE Abstract —Hyperspectral unmixing is a blind source separation problem which consists in estimating the reference spectral signatures contained in a hyperspectral image, as well as their relati ve contribution to each pixel according to a given mixture model. In practice, the process is further complexified by the inherent spectral variability of the observed scene and the possible presence of outliers. More specifically , multi-temporal hyperspectral images, i.e., s equences of h yperspectral images acquired over the same area at different time instants, are likely to simultaneously exhibit moderate endmember variability and abrupt spectral changes either due to outliers or to significant time intervals between consecutive acquisitions. Unless properly accounted f or , these two perturbations can significantly affect the unmixing process. In this context, we propose a new unmixing model for multitemporal h yperspectral images accounting for smooth temporal variations, construed as spectral variability , and abrupt spectral changes interpr eted as outliers. The proposed hierarchical Bayesian model is inferr ed using a Markov chain Monte-Carlo (MCMC) method allo wing the posterior of interest to be sampled and Bayesian estimators to be approximated. A comparison with unmixing techniques from the literatur e on synthetic and real data allows the interest of the proposed approach to be appreciated. Index T erms —Hyperspectral imagery, multitemporal images, endmember variability , Markov chain Monte-Carlo (MCMC) methods. I . I N T R O D U C T I O N A CQUIRED in hundreds of contiguous spectral bands, hyperspectral (HS) images ha ve received an increasing interest due to the significant spectral information they conv ey , which is somewhat mitigated by their lower spatial resolution in remote sensing applications. This limitation, combined with possibly complex interactions between the incident light and the observed materials, implies that the observed spectra are This work was supported in part by the Hypanema ANR Project no. ANR-12-BS03-003, by the MapIn vPlnt ERA-NET MED Project no. ANR-15- NMED-0002-02, by the Thematic T rimester on Image Processing of the CIMI Labex under Grant ANR-11-LABX-0040-CIMI within the Program ANR-11- IDEX-0002-02 and by the Direction G ´ en ´ erale de l’Armement, French Ministry of Defence. The authors are with the Univ ersity of T oulouse, IRIT/INP-ENSEEIHT , 31071 T oulouse, France. (e-mail: { pierreantoine.thouv enin, Nicolas.Dobigeon, Jean-Yves.T ourneret } @enseeiht.fr This paper has supplementary downloadable material a vailable at http://ieeexplore.ieee.or g., provided by the authors. The material includes complementary illustrations to the main paper , as well as further details on the values chosen for the hyperparameters inv olved in the experiments. Contact pierreantoine.thouvenin@enseeiht.fr for further questions about this work. mixtures of several signatures corresponding to distinct mate- rials. Spectral unmixing then consists in identifying a limited number of reference spectral signatures composing the data – referred to as endmembers – and their abundance fractions in each pixel according to a predefined mixture model. The choice of a specific model generally reflects the practitioners’ prior kno wledge on the en vironmental factors possibly affect- ing the acquisitions, such as declivity or multiple reflections. T raditionally , a linear mixing model (LMM) is adopted since it is appropriate to describe hyperspectral data when the declivity of the scene and microscopic interactions between the observed materials are negligible [ 1 ]. Depending on the applications, various models ha ve also been in vestig ated to capture higher order interactions (i.e., nonlinearities) between the incident light and the observ ed materials (see [ 2 ], [ 3 ] for recent reviews on this topic). Howev er, varying acquisition conditions, such as local illumination variations or the natural ev olution of the scene, may significantly alter the shape and the amplitude of the acquired signatures [ 4 ], [ 5 ], thus affecting the extracted endmembers. Endmember variability has hitherto been extensiv ely considered within a single HS image, either in a deterministic [ 6 ]–[ 9 ] or a statistical setting [ 10 ]–[ 12 ]. Recent works also considered temporal variability by ex- ploiting the possibilities offered by multitemporal HS (MTHS) images [ 13 ], [ 14 ]. From a hyperspectral unmixing perspective, MTHS images, i.e., sequences of HS images acquired ov er the same area at different time instants, can be of interest to exploit information redundancy between consecuti ve images (e.g., through features exhibiting moderate or smooth temporal variations as in [ 15 ], [ 16 ]) while allowing the endmember temporal e v olution to be characterized. For instance, MTHS hav e been recently exploited to improve endmember unmixing results [ 13 ], [ 14 ], [ 17 ] and used in a change detection problem in volving two HS images [ 18 ], [ 19 ]. Even though the approaches proposed in [ 13 ], [ 14 ], [ 17 ] specifically allo w smooth temporal variations of some of the mixture parameters to be considered, they do not account for abrupt spectral changes either due to outliers or to possibly significant time intervals between two consecuti ve images. In practice, such situations can be reasonably expected, depend- ing on the acquisition dates and possible climatic hazards, e.g., when vegetation or water is present in the observ ed scene. Unless specifically accounted for , this situation fre- 2 quently observ ed in real datasets has a significant impact on the recovered endmembers, which motiv ates the present work. Inspired by [ 13 ], [ 20 ], [ 21 ], and based on an original interpretation of the unmixing problem under study , our contri- bution consists in jointly accounting for smooth endmember variations – construed as temporal endmember v ariability – and abrupt changes interpreted as outliers (e.g., significant variability within a single image or presence of non-linearities) using a carefully designed hierarchical Bayesian model. More precisely , we focus our analysis on scenes in which mostly the same materials are expected to be observed from an image to another . In this context, using the endmembers extracted from the reference scene as a starting point to unmix the whole MTHS image constitutes a reasonable attempt to generalize the analyses previously conducted for a single image. On the one hand, the endmembers identified in each single image can in fine be considered as time-varying instances of reference signatures shared by the different images, thus justifying the use of a modified version of the perturbed linear mixing model (PLMM) proposed in [ 14 ]. This formulation will notably allo w smooth spectral variations occurring over time to be captured, leading to competitive results when compared to methods analyzing the images individually . On the other hand, the signatures corresponding to materials appearing in only a few images, which induce abrupt spectral changes, can be regarded as outliers with respect to the commonly shared endmembers. This paper studies a new Bayesian model allowing both spectral variability and presence of outliers to be considered in the unmixing of MTHS images. The resulting unmixing task is solved using a Markov chain Monte-Carlo (MCMC) allowing the posterior of interest to be sampled and Bayesian estimators to be approximated. The paper is organized as follows. The mixing model considered in this paper is introduced in Section II , and the associated hierarchical Bayesian model is de veloped in Section III . Section IV in vestig ates a Gibbs sampler to solve the resulting mixed integer non-linear problem. The performance of the proposed approach on synthetic and real data is studied in Sections V and VI . In particular , the results obtained with the proposed method are compared to those of the VCA/FCLS algorithm [ 22 ], [ 23 ], the SISAL/FCLS algorithm [ 24 ], the algorithm associated with the robust LMM (RLMM) proposed in [ 25 ] and the MTHS optimization method [ 14 ]. Finally , Section VII concludes this work and outlines further research perspectiv es. I I . P RO B L E M S T AT E M E N T W e consider a sequence of HS images acquired at T different time instants o ver the same area, where mostly the same materials are expected to be observed ov er time. In the following, at most R endmembers are assumed to be shared between the T images composing the sequence, where R is a priori known. Since the observed instances of a giv en endmember can be reasonably expected to vary from an image to another , we propose to account for smooth endmember spectral variations via a modified version of the perturbed linear mixing model (PLMM) proposed in [ 9 ], [ 14 ]. Inspired by the total least squares problem [ 26 ], the PLMM consists in representing each pixel y nt by a linear combination of the R endmembers – denoted by m r – af fected by an additi ve error term dm r,t accounting for temporal endmember v ariability . Howe ver , this model sho ws notable limitations when the v ector y n,t is affected by abrupt changes. Consequently , this paper in vestigates a new unmixing model jointly accounting for endmember variability and abrupt changes possibly affecting MTHS images. T o this end, the proposed model is a general- ized PLMM, which includes an additional term x n,t to capture significant de viations from the LMM, i.e., significant spatial variability or non-linearities within each image [ 20 ], [ 25 ]. The resulting observation model can thus be written y n,t = R X r =1 a r,n,t  m r + dm r,t  + x n,t + b n,t (1) for n = 1 , . . . , N and t = 1 , . . . , T , where y n,t denotes the n th image pix el at time t , m r is the r th endmember , a r,n,t is the proportion of the r th endmember in the n th pixel at time t , dm r,t denotes the perturbation of the r th endmember at time t , and x n,t denotes the contribution of outliers in the n th pix el at time t . Finally , b n,t represents an additiv e noise resulting from the data acquisition and the modeling errors. The so-called robust PLMM can be written Y t = ( M + dM t ) A t + X t + B t (2) where Y t = [ y 1 ,t , . . . , y N ,t ] is an L × N matrix containing the pixels of the t th image, M denotes an L × R matrix containing the endmembers that are common to all the images of the sequence, A t is an R × N matrix composed of the abundance vectors a n,t , dM t is an L × R matrix whose columns contain the variability inherent to the t th image, X t is an L × N matrix whose columns are the outliers present in the image t , and B t is an L × N matrix accounting for the noise at time t . The constraints considered to reflect physical considerations are A t  0 R,N , A T t 1 R = 1 N , ∀ t ∈ { 1 , . . . , T } M  0 L,R , M + dM t  0 L,R , ∀ t ∈ { 1 , . . . , T } X t  0 L,N , ∀ t ∈ { 1 , . . . , T } (3) where  denotes a term-wise inequality . Note that the outlier term X t is intended to describe abrupt changes due for instance to the appearance of one or se veral new endmembers that were not present in the reference image. This justifies the corresponding non-negati vity constraint, similar to the one imposed on the other endmembers. Note ho we ver that dif ferent phenomena not considered in this work, possibly represented by the terms X t , can induce a decrease in the total reflectance, e.g., shado wing ef fects or some nonlinearities as detailed in [ 27 ]. T o address this case, the non-negati vity constraint on the outlier terms X t should be removed. Giv en the mixture model ( 2 ), the unmixing problem con- sidered in this work consists in inferring the ab undances A t , the endmembers M , the v ariability dM t and the outliers X t from the observations Y t , t = 1 , . . . , T . In the next section, this problem is tackled in a Bayesian framew ork to easily incorporate all the prior kno wledge av ailable on the mixture parameters. 3 I I I . B AY E S I A N M O D E L This section details the specific structure imposed on the parameters to be inferred via appropriate prior distrib utions. Note that dependencies with respect to constant parameters are omitted in the following paragraphs to simplify the notations. A. Likelihood Assuming the additiv e noise b n,t is distributed according to a Gaussian distribution b n,t ∼ N ( 0 L , σ 2 t I L ) , the observation model ( 2 ) leads to y n,t | M , dM t , A t , X t , σ 2 t ∼ N  ( M + dM t ) a n,t + x n,t , σ 2 t I L  . In addition, assuming prior independence between the pix els within each image and between the images Y t themselves, the likelihood function of all images Y ¯ = [ Y 1 , . . . , Y T ] is p ( Y ¯ | Θ) ∝ T Y t =1 ( σ 2 t ) − N L/ 2 × exp  − 1 2 σ 2 t k Y t − ( M + dM t ) A t − X t k 2 F  (4) where the underline notation stands for the overall set of the corresponding parameters, k·k F is the Frobenius norm, Θ = { Θ p , Θ h } and Θ p = { M , dM ¯ , A ¯ , X ¯ , σ 2 , Z } , Θ h = { Ψ 2 , s 2 , β } (5) denote the parameters and hyperparameters whose priors are defined in the follo wing paragraphs. Note that the indepen- dence assumption between the observed images conditionally on the unknown parameters is justified by the fact that the sequence of images has been acquired by possibly dif ferent sensors at different time instants. Remark. The proposed method can easily accommodate dif- ferent structures for the noise covariance matrix (e.g., diagonal or full covariance matrix) in case the correlation between the spectral bands is significant (see, e.g., [ 28 ]). Howe ver , this modification would increase the computational and memory cost of the estimation algorithm introduced in Section IV . B. P ar ameter priors 1) Abundances: W e propose to promote smooth temporal variations of the abundances between successi ve time instants for pixels that are not classified as outliers. T o this end, we first introduce the binary latent v ariables z t ∈ { 0 , 1 } N to describe the support of the outliers (i.e., z n,t = 0 in the absence of outliers in the pixel ( n, t ) , 1 otherwise). W ith this notation, we introduce a new ab undance prior defined for n = 1 , . . . N as a n, 1 | z n, 1 = 0 ∼ U S R (6) a n,t | z n,t = 1 ∼ U f S R , for t = 1 , . . . , T (7) p  a n,t | z n,t = 0 , A ¯ \{ a n,t }  ∝ exp  − 1 2 ε 2 n ×  [ T 1 n,t 6 = ∅ ] k a n,t − a n,τ 1 n,t k 2 2   1 S R ( a n,t ) , for t ≥ 2 (8) where U S R denote the uniform distribution on the set S R , 1 S R is the indicator function of the set S R , [ P ] denotes the Iverson bracket applied to the logical proposition P , i.e., [ P ] =  1 , if P is true; 0 , otherwise and S R = { x ∈ R R | ∀ i, x i ≥ 0 and x T 1 R = 1 } (9) f S R = { x ∈ R R | ∀ i, x i ≥ 0 and x T 1 R ≤ 1 } (10) T 1 n,t = { τ < t | z n,τ = 0 } , τ 1 n,t = max τ ∈ T 1 n,t τ . (11) By con vention, we set T 1 n,t = ∅ when t = 1 . T o be more explicit, consider an image at time t and a pixel n within this image which is not corrupted by outliers (i.e., z n,t = 0 ). For t = 1 , a uniform distribution defined in the unit simplex is selected to reflect the absence of specific prior knowledge while accounting for the related constraints in ( 3 ). For t > 1 , smooth v ariations of a n,t are promoted via a one- dimensional Gaussian Markov field [ 12 ], [ 29 ] penalizing the Euclidean distance between a n,t and the abundance of the last corresponding outlier-free pixel in the preceding images of the sequence, i.e., at time instant τ 1 n,t . On the contrary , when outliers are present in the pixel ( n, t ) ( x n,t = 1 ), the usual abundance sum-to-one constraint is relaxed ( a T n,t 1 R ≤ 1 ) so that the prior allows cases in which the linear model does not exhaustiv ely describe the data to be addressed. Note that the a priori independence assumptions between the abundance vectors a n,t (conditionally to the labels z n,t ) is reasonable from a physical point of vie w , since the y can ev olve independently from a pix el to another . In the following, the joint abundance prior is denoted by p ( A ¯ | Z ) = N Y n =1  J n Y j =1 t j : z n,t j =1 p ( a n,t j | z n,t j = 1)  ×  I n Y i =1 t i : z n,t i =0 p ( a n,t i | a n,t i − 1 , z n,t i = 0)  (12) with I n = ] { t : z n,t = 0 } , J n = T − I n and ] denotes the cardinal operator . Note that the e vents [ z n,t = 0] and [ x n,t = 0 L ] (respectiv ely [ z n,t = 1] and [ x n,t 6 = 0 L ] ) are equiv alent, which allows p ( A ¯ | X ¯ ) to be defined. In the following paragraph, the latent variables z n,t are assigned a specific prior to reflect the fact that outliers are a priori assumed to represent a limited number of pixels within the sequence of image. 2) Outliers X ¯ and label maps Z : Similarly to [ 20 ], out- liers are a priori assumed to be spatially sparse. Different approaches hav e been proposed in the literature to include this prior knowledge, either relying on the ` 1 penalty (such as the LASSO [ 30 ]) or on mixtures of probability distrib utions in volving a Dirac mass at zero and a continuous probabil- ity distribution [ 31 ] (such as the Bernoulli-Laplace [ 32 ] or Bernoulli-Gaussian distributions [ 33 ], [ 34 ], e xtensively used in 4 the literature [ 35 ]–[ 37 ]). In this work, we propose to assign the following prior to the outliers x n,t to promote spatial sparsity p ( x n,t | z n,t , s 2 t ) = (1 − z n,t ) δ ( x n,t ) + z n,t N R L + ( 0 L , s 2 t ) (13) where N R L + denotes a Gaussian distribution truncated to the set R L + . Note that z n,t = 1 if an outlier is present in the corresponding pixel, and 0 otherwise. The proposed prior notably allows outliers to be a priori described by a truncated Gaussian distrib ution when z n,t = 1 , since the outliers x n,t are mainly due to the appearance of new endmembers (i.e., that were not present in the reference image). W ith this context in mind, we further propose to promote spatial correlations between the outliers’ support, since ne w materials are likely to appear in multiple contiguous pixels. The binary label maps z t ∈ R N ( t = 1 , . . . , T ) are consequently modeled as Ising-Markov random fields [ 12 ], [ 38 ], [ 39 ], for which the Hammersley-Clif ford theorem yields p ( z t | β t ) = 1 C ( β t ) exp  β t N X n =1 X k ∈V ( n ) δ ( z n,t − z k,t )  (14) where V ( n ) denotes the 4-neighbourhood of the pixel n , and C ( β t ) is the partition function [ 40 ]. In practice, the outlier terms x n,t can be assumed to be a priori independent conditionally on z n,t (since the values of outliers are not a priori correlated over space and time from a physical point of vie w). A similar assumption can be made on the labels z t , which leads to p ( X ¯ | Z , s 2 ) = Y n,t p ( x n,t | z n,t , s 2 t ) (15) p ( Z | β ) = Y t p ( z t | β t ) (16) with Z ∈ R N × T and β ∈ R T . Note that the prior ( 13 ) leads to the following result, which will be useful to sample the label maps in Section IV -E p ( x n,t | z \ n,t , s 2 t , β t ) = (1 − ω n,t ) δ ( x n,t ) + ω n,t N R L + ( 0 L , s 2 t I L ) where z \ n,t denotes the label map z t whose n th entry has been removed, and ω n,t = 1 C exp  β t X k ∈V ( n ) δ (1 − z k,t )  . (17) with C = P 1 i =0 exp  β t P k ∈V ( n ) δ ( i − z k,t )  . 3) Endmembers: A non-informativ e prior is adopted for the endmember matrix M to reflect the absence of specific prior knowledge about the spectral signatures contained in the image. More precisely , as in pre vious studies related to hyperspectral unmixing [ 12 ], [ 20 ], we consider the follo wing truncated multiv ariate Gaussian distribution m r ∼ N R L + ( 0 L , ξ I L ) , for r = 1 , . . . , R (18) where ξ is set to a sufficiently large value to ensure an uninformativ e prior (e.g., ξ = 1 ). Assuming the endmembers m r are independent (which is physically reasonable since the endmembers characterize different materials), the joint prior for the endmembers can be written as p ( M ) = R Y r =1 p ( m r ) . (19) In addition, the endmembers can be a priori assumed to li ve in a subspace of dimension K = R − 1 [ 41 ] whose practical determination can be performed by a principal component analysis (PCA) or a robust PCA (rPCA) [ 42 ]. This dimen- sionality reduction step is essentially aimed at reducing the computational complexity of the proposed approach. More explicitly , the PCA applied to the original data Y ¯ leads to a decomposition which can be expressed as [ 41 ] m r = Ue r + ˇ y , ˇ y = ( I L − UU T ) ¯ y , U T U = I K (20) where U denotes a basis of the subspace of dimension K and ¯ y denotes the average spectral signature obtained from Y ¯ . Note that using an rPCA would result in similar expres- sions (modulo a simple change of notations). The projected endmembers e r are then assigned the following truncated Gaussian prior , which ensures the non-ne gati vity of the end- members e r ∼ N E r ( 0 K , ξ I K ) , for r = 1 , . . . , R (21) with E r = [ e − 1 ,r , e + 1 ,r ] × . . . × [ e − K,r , e + K,r ] (22) e − k,r = max ` ∈U + k  − ˇ y ` + P j 6 = k u `,j e j,r u `,k  (23) e + k,r = min ` ∈U − k  − ˇ y ` + P j 6 = k u `,j e j,r u `,k  (24) U − k = { r : u k,r < 0 } , U + k = { r : u k,r > 0 } . (25) 4) Endmember variability: W e consider a prior for the vectors dm r,t (associated with the endmember variability) promoting smooth temporal variations while accounting for the term-wise non-ne gati vity of the observ ed endmembers (i.e., m r + dm r,t  0 L,R ), expressed as dm r, 1 | m r ∼ N I r ( 0 L , ν I L ) (26) dm `,r,t | m `,r , dm `,r, ( t − 1) , ψ 2 `,r ∼ N I `,r  dm `,r, ( t − 1) , ψ 2 `,r  (27) for ` = 1 , . . . , L , r = 1 , . . . , R , t = 1 , . . . , T , where I r = I 1 ,r × . . . × I L,r and I `,r = [ − m `,r , + ∞ ) . Assuming a priori independence between the dif ferent endmember variabilities (since the variabilit y can be independent from a material to another), the joint v ariability prior can finally be expressed as p ( dM ¯ | M , Ψ 2 ) = R Y r =1  p ( dm r, 1 | m r ) × T Y t =2 p ( dm r,t | m r , dm r, ( t − 1) , ψ 2 r )  . (28) 5) Noise variance: A non-informati ve inv erse-gamma con- jugate prior is selected for the noise variance σ 2 t ∼ I G ( a σ , b σ ) (29) for t = 1 , . . . , T , with a σ = b σ = 10 − 3 in order to ensure a weakly informative prior . The noise variances σ 2 t can be 5 Y ¯ M ξ dM ¯ ν Ψ 2 a Ψ b Ψ σ 2 a σ b σ A ¯ ε 2 X ¯ Z s 2 β a s b s Fig. 1. Directed acyclic graph associated with the proposed Bayesian model (fixed parameters appear in boxes). assumed to be a priori independent (giv en the absence of a priori correlation between the noise in dif ferent images), thus leading to p ( σ 2 ) = Y t p ( σ 2 t ) . (30) C. Hyperparameters In order to complete the description of the proposed hier- archical Bayesian model, we consider the following generic priors for the dif ferent hyperparameters. Note that the a priori independence assumptions made in this section are more of a computational nature, i.e., aimed at simplifying the estimation procedure detailed in the next section. (i) Non-informativ e conjugate in verse-gamma priors for the variability variances Ψ 2 and the outlier variances s 2 , i.e., for ` = 1 , . . . , L , r = 1 , . . . , R and t = 1 , . . . , T ψ 2 `,r ∼ I G ( a Ψ , b Ψ ) , s 2 t ∼ I G ( a s , b s ) (31) where I G ( a Ψ , b Ψ ) denotes the in verse gamma distribu- tion and a Ψ = b Ψ = a s = b s = 10 − 3 . Classical indepen- dence assumptions for the different hyperparameters lead to p ( Ψ 2 ) = Y `,r p ( ψ 2 `,r ) , p ( s 2 ) = Y t p ( s 2 t ) . (32) (ii) A uniform prior for the granularity parameter of a Potts- Markov random field ( a fortiori of an Ising-Marko v random field). Pre vious studies ha ve shown that it is rea- sonable to constrain the granularity parameter to belong to the interval [0 , 2] [ 43 ], leading to β t ∼ U [0 , 2] , for t = 1 , . . . , T . (33) Assuming the granularity parameters are a priori inde- pendent for different time instants finally yields p ( β ) = Y t p ( β t ) . (34) D. Joint posterior distribution Applying Bayes’ theorem, the joint posterior distribution of the parameters of interest is giv en by p (Θ | Y ¯ ) ∝ p ( Y ¯ | Θ) p ( A ¯ | X ¯ ) p ( X ¯ | Z , s 2 ) p ( s 2 ) × p ( Z | β ) p ( β ) p ( dM ¯ | M , Ψ 2 ) p ( M ) p ( Ψ 2 ) p ( σ 2 ) . (35) The complexity of the proposed Bayesian model summa- rized in the directed ac yclic graph of Fig. 1 and its resulting posterior ( 35 ) prevent a simple computation of the maximum a posteriori (MAP) or minimum mean square (MMSE) es- timators. For instance, the optimization problem associated with the determination of the MAP estimator of Θ is clearly complex, since the ne gati ve log-posterior is non-conv ex and parameterized by mixed continuous and discrete variables. In this context, classical matrix factorization techniques such as [ 25 ] cannot be used efficiently . An MCMC method is consequently adopted to sample the posterior ( 35 ) and to build estimators of the parameters in volv ed in the proposed Bayesian model using the generated samples. I V . G I B B S S A M P L E R This section studies a Gibbs sampler, which is guaranteed to produce samples asymptotically distributed according to the target distribution ( 35 ). This sampler described in Algo. 1 consists in generating samples distributed according to the conditional distribution of each parameter of interest. Sec- tion IV -A introduces the proposed sampling method, and the conditional distributions of all the parameters of interest (see Fig. 1 ) are detailed in the following paragraphs. A. Bayesian inference and parameter estimation The main steps of the proposed Gibbs sampler are summa- rized in Algo. 1 . Similarly to [ 20 ], the sequence { Θ ( q ) } N MC q = N bi +1 generated by the proposed sampler (i.e., after N bi burn-in iterations) is used to approximate the MMSE estimators of the different unkno wn parameters M , A t , dM t and X t by replacing the expectations by empirical a verages. c M MMSE ' 1 N MC − N bi N MC X q = N bi +1 M ( q ) (36) b A MMSE t ' 1 N MC − N bi N MC X q = N bi +1 A ( q ) t (37) d dM MMSE t ' 1 N MC − N bi N MC X q = N bi +1 dM ( q ) t (38) b X MMSE t ' 1 N MC − N bi N MC X q = N bi +1 X ( q ) t . (39) This choice is justified by the fact that MAP estimators computed using MCMC algorithms are often less accurate when the number of unknown parameters is relativ ely large [ 44 ]. Finally , the following marginal maximum a posterior (mMAP) estimator is considered for the label maps b z mMAP n,t = arg max z n,t ∈{ 0 , 1 } p  z n,t | y n,t , Θ \{ z n,t }  . (40) It is approximated by b z mMAP n,t '  0 , if ] { q > N bi : z ( q ) n,t = 0 } ≤ N MC − N bi 2 1 , otherwise . (41) 6 Algorithm 1: Proposed hybrid Gibbs sampler . Data: N bi , N MC , M (0) , A ¯ (0) , dM ¯ (0) , X ¯ (0) , σ 2(0) , Z (0) , β (0) , s 2(0) , Ψ 2(0) . begin for q = 1 to N MC do Sample the endmembers M ( q ) , cf. § IV -C ; Sample the variability terms dM ¯ ( q ) , cf. § IV -D ; Sample the abundances A ¯ ( q ) , cf. § IV -B ; Sample the labels and the outliers ( Z ( q ) , X ¯ ( q ) ) , cf. § IV -E ; Sample the outlier variances s 2( q ) , cf. § IV -F ; Sample the noise variances σ 2( q ) , cf. § IV -G ; Sample the variability variances Ψ 2( q ) , cf. § IV -H ; Sample the granularity parameters β ( q ) , cf. § IV -I ; Result: n M ( q ) , dM ¯ ( q ) , A ¯ ( q ) , Z ( q ) , X ¯ ( q ) , σ 2( q ) , Z ( q ) , β ( q ) , s 2( q ) , Ψ 2( q ) o N MC q =1 . B. Sampling the ab undances A ¯ The likelihood function ( III-A ) combined with the prior giv en in Section III-B1 leads to the following conditional distribution for the abundances a n,t | y n,t , Θ \{ a n,t } ∼ N S R ( µ ( A ) n,t , Λ n,t ) (42) Λ − 1 n,t = 1 σ 2 t M T t M t + 1 ε 2 n  [ T 1 n,t 6 = ∅ ] + [ T 2 n,t 6 = ∅ ]  I R (43) M t , M + dM t (44) µ ( A ) n,t = Λ n,t  1 σ 2 t M T t ( y n,t − x n,t ) + 1 ε 2 n  [ T 1 n,t 6 = ∅ ] a n,τ 1 n,t + [ T 2 n,t 6 = ∅ ] a n,τ 2 n,t   (45) where N S R ( µ , Λ ) denotes a Gaussian distribution truncated to the set S R and T 2 n,t = { τ > t | z n,τ = 0 } , τ 2 n,t = min τ ∈ T 2 n,t τ (46) with the conv ention T 2 n,t = ∅ if t = T . Samples distrib uted according to the above truncated mul- tiv ariate Gaussian distributions can be generated by a Gibbs sampler described in [ 45 , Section IV .B.] [ 46 ], by an Hamil- tonian Monte-Carlo procedure [ 47 ], [ 48 ] or by the general method recently proposed in [ 49 ]. In this work, the Gibbs sampler [ 45 , Section IV .B.] has been adopted to sample the parameters of interest. Note that the ab undance v ectors a n,t can be sampled in parallel to accelerate the algorithm. C. Sampling the endmember s M Combining ( III-A ) and the endmember prior giv en in Sec- tion III-B3 leads to m `,r | Y ¯ , Θ \{ m `,r } ∼ N [ b `,r , + ∞ ) ( µ ( M ) `,r , κ 2 `,r ) (47) b `,r = max n 0 , max t ( − dm `,r,t ) o (48) µ ( M ) `,r = κ 2 `,r X t 1 σ 2 t  e y `,t − e x `,t − e m l, \ r A \ r,t − g dm `,t A t  e a T r,t (49) κ 2 `,r = " X n,t a 2 r,n,t σ 2 t + 1 ξ # − 1 (50) where g dm `,t is the ` th row of dM t , e m `, \ r is the ` th row of M whose r th entry has been remo ved, A \ r,t denotes the matrix A t without its r th ro w and e a r,t is the r th ro w of A t . Samples distributed according to the above truncated Gaussian distributions can be efficiently generated using the algorithm described in [ 50 ]. When using a PCA as a preprocessing step ( 20 ), the projected endmembers e r , for r = 1 , . . . , R have a truncated multiv ariate Gaussian distribution [ 41 ] e r | Y ¯ , Θ \{ e r } ∼ N E r ( µ ( E ) r , Λ r ) (51) where E r , µ ( E ) r and Λ r hav e been reported in Appendix A . Note that the ro ws of M (resp. of the projected endmember matrix E ) can be sampled in parallel to decrease the compu- tational time required by the algorithm. D. Sampling the variability terms dM ¯ Similarly , the likelihood function ( III-A ) and the prior gi ven in Section III-B4 lead to dm `,r,t ∼ N [ − m `,r , + ∞ ) ( µ ( dM ) `,r,t , η 2 `,r,t ) (52) with 1 η 2 `,r,t = 1 σ 2 t X n a 2 r,n,t + 1 ν  t = 1  + 1 ψ 2 `,r  1 +  1 < t < T   (53) µ ( dM ) `,r,t =  1 σ 2 t  e y `,t − g dm `, \ r,t A \ r,t − e m ` a n,t − x `,n,t  ˜ a T r,t + 1 ψ 2 `,r  [ t < T ] dm `,r, ( t +1) + [ t > 1] dm `,r, ( t − 1)  # η 2 `,r,t (54) where g dm `, \ r,t denotes the ` th row of dM t whose r th element has been remov ed, e m ` is the ` th ro w of M and A \ r,t is the matrix A t without its r th row . The rows of each variability matrix dM t can be sampled in parallel to reduce the computational time of the sampler . E. Sampling the label maps Z and the outliers X ¯ According to ( III-A ) and Section III-B2 , the outliers admit the follo wing group-sparsity promoting conditional distribu- tions p ( x n,t | y n,t , Θ \{ z n,t , x n,t } ) = (1 − w n,t ) δ ( x n,t ) + w n,t N R L + ( µ ( X ) n,t , ϑ 2 t I L ) (55) which are mixtures of a Dirac mass at 0 and of truncated multiv ariate Gaussian distrib utions, where w n,t = ˜ w n,t ˜ w n,t + (1 − ω n,t ) , ϑ 2 t = σ 2 t s 2 t σ 2 t + s 2 t (56) ˜ w n,t = ω n,t ( s 2 t ) L/ 2 ( ϑ 2 t ) L/ 2 exp  1 2 ϑ 2 t k µ ( X ) n,t k 2 2  (57) µ ( X ) n,t = s 2 t σ 2 t + s 2 t  y n,t − ( M + dM t ) a n,t  . (58) 7 In practice, the labels z n,t are first sampled according to a Bernoulli distribution to select one of the two models for x n,t , with probability P [ z n,t = 1 | y n,t , Θ \{ z n,t , x n,t } ] = w n,t . Note that the labels z n,t can be sampled in parallel by using a checkerboard scheme. In addition, the outliers x n,nt can be sampled in parallel to decrease the computational time. F . Sampling the outlier variances s 2 According to Sections III-B2 and III-C , we can easily identify the conditional law of s 2 t for t = 1 , . . . , T as the following in verse gamma distribution s 2 t | Θ \{ s 2 t } ∼ I G  a s + ] { n : z n,t = 1 } L 2 , b s + 1 2 k X t k 2 F  . (59) G. Sampling the noise variances σ 2 Using Sections III-B5 and III-C , we obtain for t = 1 , . . . , T σ 2 t | Y t , Θ \{ σ 2 t } ∼I G  a σ + LN 2 , b σ + 1 2 k Y t − ( M + dM t ) A t − X t k 2 F  . (60) H. Sampling the variability variances Ψ 2 Similarly , Section III-B4 and III-C lead to ψ 2 `,r | Θ \{ ψ 2 `,r } ∼I G  a Ψ + T − 1 2 , b Ψ + 1 2 T X t =2 ( dm `,r,t − dm `,r,t − 1 ) 2  . (61) I. Sampling the gr anularity parameters β t Provided square images are considered, the partition func- tions C ( β t ) have the closed-form expressions [ 40 ], [ 51 ] e C ( β t ) = 1 2 log(2 sinh β t ) + 1 2 N N X n =1 acosh ∆ n ( β t ) + β t ∆ n ( β t ) = v ( β t ) − C n , v ( β t ) = cosh 2 β t sinh β t e C ( β t ) = 1 N log C ( β t ) , C n = cos  2 n − 1 2 N π  . The exact partition function can then be used to sample the parameters β t using Metropolis-Hastings steps. In this work, new v alues of the granularity parameters have been proposed by the following Gaussian random w alk β ∗ t = β ( q ) t + ε t , ε t ∼ N  0 , σ 2 β ( t )  (62) where the parameters σ 2 β ( t ) are adjusted during the b urn-in iterations to yield acceptance rates in the interv al [0 . 4 , 0 . 6] . J. Computational complexity Assuming elementary arithmetic operations and scalar pseudo-random number generations are O (1) operations, the ov erall computational complexity is dominated by matrix products needed to compute the parameters related to the (a) Endmember 1 (b) Endmember 2 (c) Endmember 3 Fig. 2. Endmembers ( m r , red lines) and their variants affected by variability ( m r + dm r,t , blue dotted lines) used to generate the synthetic mixtures with R = 3 . Signatures corresponding to different time instants are represented in a single figure to better appreciate the variability introduced in the data. conditional distribution of the v ariability vectors. Since R  L  N and T  L , the per-iteration computational cost of the proposed algorithm is O ( LR 2 N T ) per iteration. As detailed in the preceding paragraphs, many parameters can be sampled in parallel to reduce the computational time of the proposed algorithm. In comparison, the computational complexity of VCA is O ( R 2 N ) [ 22 ] per image, and the per iteration complexity of the others algorithms for a single image are respectively: O ( N 2 ) for FCLS [ 23 ], O ( RN ) for SISAL [ 24 ], O ( LRN ) for rLMM [ 25 ] and O ( R 2 ( L + N )) for OU [ 14 ]. V . E X P E R I M E N T S W I T H S Y N T H E T I C D A TA The proposed method has been applied to an MTHS image composed of 10 acquisitions of size 50 × 50 with L = 413 bands. The first scenario deals with the appearance of a ne w material in specific re gions of a fe w images. T o this end, 4 images out of the 10 hav e been corrupted by spatially sparse outliers, corresponding to a ne w endmember extracted from a spectral library . Each image of the sequence corresponds to a linear mixture of 3 endmembers affected by smooth time-varying variability , and the synthetic ab undances v ary smoothly from one image to another . First, so-called reference ab undance maps corresponding to the first time instant have been generated (e.g., for Fig. 4 , we hav e taken the abundance maps obtained by VCA/FCLS on the widely studied Moffett dataset [ 41 ]). Then, the abundance maps corresponding to the remaining time instants have been generated by multiplying the reference maps with trigonomet- ric functions to ensure a sufficiently smooth temporal ev olu- tion. For the first dataset composed of R = 3 endmembers, the reference maps associated with the first two endmembers hav e been respectively multiplied by cos  π 100 + t 48 π 100  and sin  π 100 + t 48 π 100  , with t ∈ { 1 , . . . , T } . The temporal e volution of the last abundance map has finally been obtained by lev er- aging the sum-to-one condition. W ith these abundance maps, the contribution of a gi ven endmember (assumed to punctually disappear) has been replaced at specific time instants by a new endmember signature in pixels originally corresponding to its highest abundance coefficients (e.g., abov e 0.8). The mixtures have finally been corrupted by an additive white Gaussian noise to ensure a resulting signal-to-noise ratio (SNR) between 25 and 30 dB . Similarly , two complementary scenarii in volving 5 HS images, of size 100 × 100 , composed of 6 and 9 endmembers, hav e been considered to analyze the performance of the method in the presence of a larger 8 T ABLE I F I XE D PAR A M E TE R S , A N D I N I T I A L V A L U E S A S S O C I A T ED I N T H E E X PE R I M EN T S T O P A R AM E T E RS L ATE R I N F E R R E D F R O M T H E M OD E L . Synthetic data Real data Fixed parameters ε 2 n 10 − 3 10 − 2 ξ 1 1 ν 10 − 3 10 − 5 a s , a Ψ , a σ 10 − 3 10 − 3 b s , b Ψ , b σ 10 − 3 10 − 3 N bi 350 450 N MC 400 500 Initial values σ 2 t 10 − 4 10 − 4 s 2 t 5 × 10 − 3 5 × 10 − 3 ψ 2 `,r 10 − 3 10 − 2 β t 1.7 1.7 number of endmembers. Note that a larger image size has been considered for these two datasets to reflect the fact that a larger number of endmembers is expected to be observed in larger scenes. In addition, the images of this experiment do not satisfy the pure pixel assumption to assess the proposed method in challenging situations. Controlled spectral v ariability has been introduced by using the product of reference endmembers with randomly gen- erated piece wise-af fine functions as in [ 9 ], where different affine functions ha ve been generated for each endmember at each time instant. T ypical instances of the signatures used in this experiment are depicted in Fig. 2 . The robustness of the proposed method to moderate spatial v ariability , i.e., endmember variability occurring within single images, has also been ev aluated. The corresponding results can be found in the technical report [ 52 , Appendix E] due to space constraints. A. Compar ed methods The results of the proposed algorithm have been compared to those of several unmixing methods from the literature, some of which are specifically designed to unmix a single HS image. In the follo wing lines, the most relev ant implementation details specific to each method are briefly recalled. 1) VCA/FCLS (no v ariability , single image): the endmem- bers are first extracted on each image using the verte x component analysis (VCA) [ 22 ], which requires pure pixels to be present. The abundances are then estimated for each pixel by solving a fully constrained least squares problem (FCLS) using the alternating direction method of multipliers (ADMM) [ 23 ]. Note that the estimates provided by the VCA algorithm vary from one run to another , gi ven its stochastic nature; 2) SISAL/FCLS (no v ariability , single image): the end- members are extracted on each image by the simplex identification via split augmented Lagrangian (SISAL) [ 24 ], and the abundances are estimated for each pixel by FCLS. The tolerance for the stopping rule has been set to 10 − 3 ; 3) RLMM (no variability , single image): the unmixing method associated with the robust linear mixing model (RLMM) proposed in [ 25 ] has been applied to each image of the series independently . The algorithm has been initialized with SISAL/FCLS, and the regularization parameter specific to this method is set as in [ 25 ]; 4) OU: the endmembers are estimated using the online un- mixing (OU) algorithm introduced in [ 14 ] with endmem- bers initialized by the output of VCA applied to the first image of the sequence. The abundances are initialized by FCLS, and the variability matrices are initialized with all their entries equal to 0 . The other parameters are set to the same values as those gi ven in [ 14 , T able I]; 5) Proposed approach: the endmembers are initialized with VCA applied to the first image of the sequence, within which the observed materials are well represented (i.e., with sufficiently high abundance coefficients for each material). In this context the VCA algorithm, which requires pure pixels to be present in the data, has been observed to yield relev ant results for the initialization. Howe ver , other endmember e xtraction techniques might be used to initialize the proposed algorithm if needed. The abundances are initialized by FCLS, and the variability matrices and label maps are initialized with all their entries equal to 0 (i.e., the images are a priori assumed to contain no outlier). The v alues chosen for the other parameters are summarized in T able I . Further details on these values can be found in the supplementary material provided by the authors. Performance assessment has been conducted in terms of (i) endmember estimation through the average spectral angle mapper (aSAM) aSAM( M ) = 1 R R X r =1 arccos  m T r b m r k m r k 2 k b m r k 2  ; (63) (ii) abundance and variability estimation through the global mean square errors (GMSEs) GMSE( A ) = 1 T R N T X t =1 k A t − b A t k 2 F (64) GMSE( dM ) = 1 T LR T X t =1 k dM t − d dM t k 2 F ; (65) (iii) quadratic reconstruction error (RE) RE = 1 T LN T X t =1 k Y t − b Y t k 2 F (66) where b Y t is the matrix composed of the pix els recon- structed with the estimated parameters. B. Results The endmembers estimated by the proposed algorithm are compared to those of VCA/FCLS, SISAL/FCLS, RLMM and OU in Fig. 3 , whereas the corresponding abundance maps are displayed in Fig. 4 . Note that the abundance maps and the endmembers obtained for the mixtures of 6 and 9 endmembers are included in a separate technical report [ 52 , Appendix D] due to space constraints (see Figs. 17–24 for R = 6 , and Figs. 26–37 for R = 9 ). The unmixing performance of each method, reported in T able II , leads to the following conclusions. 9 (a) VCA (b) SISAL (c) RLMM (d) OU (e) Proposed Fig. 3. Second endmember ( m 2 , red lines) and its variants affected by variability ( m 2 + dm 2 ,t , blue dotted lines) recovered by the different methods from the synthetic mixtures with R = 3 . Due to space restrictions, the signatures extracted for the other two endmembers have been included in the associated supplementary material. Signatures corresponding to different time instants are represented on a single figure to better appreciate the v ariability recovered from the data. The spectra represented in black correspond to signatures corrupted by outliers. • Endmember estimation: the proposed method shows an interesting robustness with respect to spatially sparse outliers in the sense that the estimated signatures (Figs. 3e , see the supplementary material for the two other endmembers) are v ery close to the corresponding ground truth (Fig. 2 ). In comparison, the shape of the end- members recovered by VCA, SISAL and RLMM and the variability extracted by OU are significantly affected by outliers, as ex emplified in Figs. 3a , 3b , 3c and 3d respectiv ely . These qualitativ e results are confirmed by the quantitativ e performance measures of each method provided in T able II . Note that the endmembers recovered by the SISAL and RLMM methods are very sensitive to the VCA initialization, as illustrated by the similarity between the signatures estimated by these methods (Figs. 3a to 3c ). • Abundance estimation: the abundance maps estimated by FCLS, RLMM and SISAL reflect the high sensiti vity of VCA (used to initialize SISAL and RLMM) to the presence of outliers (see the figures delineated in red in Fig. 4 ). On the contrary , the abundances recovered by OU and the proposed approach are much closer to the ground truth. These observations are confirmed by the abundance estimation performance reported in T able II . The proposed abundance smoothness prior appears to mit- igate the errors induced by the presence of outliers as can be seen in Fig. 6 . More precisely , for images corrupted by outliers, the abundance coefficients estimated by the proposed approach are closer to the ground truth than the results of the other methods. A more detailed v ersion of T able II , along with a complementary figure illustrating the interest of the proposed abundance prior can be found in the supplementary material. • Overall performance: the performance measures re- ported in T able II are globally fav orable to the proposed approach. It is important to mention that the price to pay with the good performance of the proposed method t = 1 TRUE t = 2 ∗ t = 3 t = 4 t = 5 ∗ 0 0.5 1 VCA 0 0.5 1 SISAL 0 0.5 1 RLMM 0 0.5 1 OU 0 0.5 1 MCMC 0 0.5 1 Fig. 4. Abundance map of the first endmember reco vered by the different methods (in each row) at the first fiv e time instants (giv en in column) for the experiment with R = 3 [the different rows correspond to the true abundances, VCA/FCLS, SISAL/FCLS, RLMM, OU and the proposed method]. The images delineated in red show that sev eral methods are highly sensitive to the presence of outliers, and the time instants represented with ∗ denote images containing outliers. Due to space restrictions, the ab undance maps obtained at each time instant for each endmember have been included in the supplementary material. is its computational complexity , which is common with MCMC methods. As a complementary output, the proposed algorithm is able to recover the location of the outliers within each image, as illustrated in Fig. 5 . Up to a few false detections, the estimated labels are very close to the ground truth. The label errors observed for t = 7 , 8 and 9 partly result from the different abundance constraints considered when an outlier is detected or not. V I . E X P E R I M E N T W I T H R E A L D A T A A. Description of the dataset W e consider a real sequence of A VIRIS HS images acquired ov er the Lake T ahoe region (California, United States of America) between 2014 and 2015 1 . The scene of interest ( 100 × 100 ), composed of a lake and a nearby field, has been unmixed with R = 3 endmembers based on the results of the noise-whitened eigengap algorithm (NWEGA) [ 53 ] applied to each image of the series (see T able III ). This choice is further supported by results obtained from a pre vious analysis conducted on the same dataset [ 54 , Appendix E]. For R = 4 and 5 , the signatures of water , soil and ve getation were split 1 The images from which the scene under study is extracted are freely av ailable from the online A VIRIS flight locator tool at http://a viris.jpl.nasa. gov/alt locator/ . 10 Estimated labels t = 1 Ground truth t = 2 t = 3 t = 4 t = 5 t = 6 t = 7 t = 8 t = 9 t = 10 Fig. 5. Outlier labels z t estimated for each image of the synthetic dataset with 3 endmembers (the different rows correspond to the true labels, and the estimated labels) [0 in black, 1 in white]. T ABLE II S I MU L A T I ON R E S U LTS O N S Y N T HE T I C DAT A ( A S A M ( M ) I N (° ) , G M SE ( A ) × 10 − 2 , G MS E ( dM ) × 10 − 4 , R E × 10 − 4 , T I M E I N ( S ) ) . aSAM( M ) GMSE( A ) GMSE( dM ) RE time R = 3 VCA/FCLS 6.07 2.32 / 3.91 1 SISAL/FCLS 5.07 1.71 / 2.28 2 RLMM 5.13 2.04 / 0.31 463 OU 1.90 0.42 3.22 2.61 98 Proposed 2.03 0.15 1.85 2.00 2530 R = 6 VCA/FCLS 3.81 1.57 / 3.09 2 SISAL/FCLS 5.76 0.91 / 4.49 3 RLMM 2.73 1.26 / 0.29 1453 OU 2.74 0.38 3.70 1.13 420 Proposed 1.48 0.16 2.84 0.51 8691 R = 9 VCA/FCLS 3.74 0.65 / 6.83 4 SISAL/FCLS 5.91 0.36 / 5.56 5 RLMM 2.48 0.54 / 0.31 1447 OU 6.08 0.47 2.19 0.89 1024 Proposed 2.23 0.15 8.38 0.82 17151 t 1 2 3 4 5 6 7 8 9 10 a 1 , 939 ,t 0.1 0.2 0.3 0.4 0.5 0.6 TRUE VCA SISAL RLMM OU MCMC Fig. 6. Evolution over time of the abundance associated with the first endmember in a giv en pixel. The similarity between the recovered result and the ground truth illustrates the relev ance of the proposed abundance prior to mitigate the errors induced by the presence of outliers in the image (time instants 2, 5, 6 and 10). into two or more components by the different algorithms, suggesting R = 3 is more appropriate for this study . Note that prior studies led in [ 14 ] rev ealed that this dataset contains outliers (area delineated in red in Fig. 7e ). After remo ving the seemingly corrupted bands and the water absorption bands, 173 out of the 224 spectral bands were finally exploited. The initial parameters used for the proposed algorithm are giv en in T able I . The other methods have been run with the same parameters as in Section V . Note that the VCA results reported in this section are representative of those obtained o ver multiple runs (no significant differences hav e been observed from one run to another). B. Results In the absence of any ground truth, the performance of the unmixing methods is assessed in terms of RE (T able IV ) while taking into account the consistency of the estimated abundance maps reported in Figs. 9 , 10 and 11 . More precisely , the abundances associated with the ve getation area are expected to be very high for t = 1 , 3 , 5 (corresponding to Figs. 7a , 7c and 7e ) where the ve getation visually appears to be sufficiently irrigated (hence well represented). On the contrary , the abundance coefficients are supposed to be much lower for t = 2 , 4 , 6 (corresponding to Figs. 7b , 7d and 7f ), where the ve getation is visually drier or almost absent. Concerning the presence of water in the bottom left-hand corner of the images, the latent variables introduced in Section III-B2 are e xpected to reflect the abrupt v ariations in the presence of water observed at t = 3 , 4 and 5 . These observations, combined with the extracted signatures (Fig. 8 ) and the estimated abundances (Figs. 9 to 11 ) lead to the follo wing comments. • Endmember estimation: the signature recovered for the soil by VCA, SISAL and RLMM at time t = 5 shows an amplitude which is significantly greater than the ampli- tude of the signatures extracted at the other time instants, and a shape incompatible with what can be expected based on physical considerations (see the black signatures in Figs. 8a , 8d and 8g ). This is a clear indication that outliers are present in the corresponding image. A similar observation can be made for the vegetation signature obtained by VCA, SISAL and RLMM at time t = 5 . On the contrary , the endmembers recovered by OU and the proposed approach are much more consistent from this point of vie w . • Abundance estimation: the estimated ab undances glob- ally reflect the previous comments made on the extracted endmembers. Notably , the abundance coef ficients esti- mated at t = 5 by VCA, SISAL and RLMM (delineated in red in Figs. 9 to 11 ) are visually inconsistent with the temporal e volution of the materials observed in the true color composition gi ven in Fig. 7 . More explicitly , the soil is not supposed to be concentrated on a fe w pixels as suggested by the corresponding ab undance maps in Fig. 9 . Similarly , the water is not supposed to be present in high proportions in all the pixels of the image as indicated 11 (a) 04/10/2014 (b) 06/02/2014 (c) 09/19/2014 (d) 11/17/2014 (e) 04/29/2015 (f) 10/13/2015 Fig. 7. Scenes used in the experiment, given with their respectiv e acquisition date. The area delineated in red in Fig. 7e highlights a region known to contain outliers (this observation results from a previous analysis led on this dataset in [ 14 ]). T ABLE III E N DM E M B ER N U MB E R R E ST I M AT E D B Y N WE G A [ 53 ] O N E AC H I M AG E O F T H E R EA L DAT A S E T . 04/10/2014 06/02/2014 09/19/201 11/17/2014 04/29/2015 10/13/2015 NWEGA 3 3 3 4 3 4 T ABLE IV S I MU L A T I ON R E S U LTS O N R E A L DA TA ( R E × 10 − 4 ) . RE time (s) R = 3 VCA/FCLS 45.05 1 SISAL/FCLS 1.65 2 rLMM 2.51 390 OU 2.50 508 Proposed 0.34 23608 in Fig. 10 . These results, in contradiction with Fig. 7 , suggest that outliers are present at t = 5 . In addition, the abundance maps estimated at t = 4 and 6 by FCLS for the water and the ve getation (delineated in green in Figs. 9 and 10 ) suggest that the water contribution has been split into two spectra. The corresponding signatures are represented in green in Figs. 8a and 8c . On the contrary , the results reported for OU and the proposed method are consistent with the e xpected ev olution of water and ve getation over time (abundance values close to 1 at time t = 1 , 3 , 5 , lower values at time t = 2 , 4 , 6 ). Finally , the ve getation abundance maps estimated by the proposed method globally presents a better contrast than those obtained with OU (Fig. 11 ). The pre vious comments, along with the lower reconstruction error reported in T able IV , suggest that the proposed approach is robust to spatially sparse outliers while allowing smooth temporal variations to be exploited. Indeed, the pixels cor- responding to abrupt v ariations of the water signature have been properly detected. Furthermore, the outliers pre viously detected in this dataset [ 14 ] for t = 5 (highlighted in red in Fig. 7e ) are well captured by the latent v ariables Z (see Fig. 12 ). In addition, the spatial distribution of the estimated outlier labels (Fig. 12 ) is in agreement with the results of the RLMM (in terms of the spatial distribution of the outlier energy) and with the non-linearity detector [ 55 ] applied to each image of the sequence with the SISAL-estimated endmembers (see Fig. 13 ). Concentrated on re gions where non-linear ef fects can be reasonably expected, the activ e latent variables Z tend to capture the spatial distribution of the non-linearities possibly occurring in the observ ed scene. V I I . C O N C L U S I O N A N D F U T U R E W O R K This paper introduced a Bayesian model accounting for both smooth and abrupt variations possibly occurring in multitemporal hyperspectral images. The adopted model was specifically designed to handle datasets in which mostly the same materials were expected to be observed at dif ferent time instants, thus allo wing information redundancy to be exploited. An MCMC algorithm was derived to solve the resulting unmixing problem in order to precisely assess the performance of the proposed approach on multitemporal HS images of moderate size (i.e., moderate spatial and temporal dimensions). This algorithm was used to sample the posterior of the proposed hierarchical Bayesian model and to use the generated samples to build estimators of the unkno wn model parameters. Giv en its computational cost, the proposed approach is not intended to be applied to lar ge datasets, for which different unmixing methods can provide a rougher analysis at a smaller computational cost. The proposed approach is rather meant to be used as a complementary tool to carry out an in-depth anal- ysis of scenes of moderate size. Future research perspectiv es include the use of relaxation methods to the Ising field to tackle similar problems with online optimization techniques, and the dev elopment of distrib uted unmixing procedures to ef ficiently unmix lar ger datasets. Designing unmixing methods scaling with the problem dimension while simultaneously accounting for temporal and spatial endmember variability is another interesting prospect. A P P E N D I X A S A M P L I N G T H E P RO J E C T E D E N D M E M B E R S E When using a PCA as a preprocessing step, the projected endmembers e r , for r = 1 , . . . , R , are distributed according to the following truncated Gaussian distrib utions e r | Y ¯ , Θ \{ e r } ∼ N E r ( µ ( E ) r , Λ r ) (67) with E r = [ c 1 ,r , d 1 ,r ] × . . . × [ c K,r , d K,r ] , and for k = 1 , . . . , K c k,r = max ` ∈U + k  − ˇ y ` + P j 6 = k u `,j e j,r + b `,r u `,k  d k,r = min ` ∈U − k  − ˇ y ` + P j 6 = k u `,j e j,r + b `,r u `,k  b `,r = min n 0 , min t ( dm `,r,t ) o 12 (a) Soil (VCA) (b) W ater (VCA) (c) V eg. (VCA) (d) Soil (SISAL) (e) W ater (SISAL) (f) V eg. (SISAL) (g) Soil (RLMM) (h) W ater (RLMM) (i) V e g. (RLMM) (j) Soil (OU) (k) W ater (OU) (l) V eg. (OU) (m) Soil (Prop.) (n) W ater (Prop.) (o) V e g. (Prop.) Fig. 8. Endmembers ( m r , red lines) and their variants affected by variability ( m r + dm r,t , blue dotted lines) recovered by the different methods from the real dataset depicted in Fig. 7 . The spectral gaps in the recov ered signatures correspond to the lo w SNR bands which have been remo ved prior to the unmixing procedure. Signatures corresponding to different time instants are represented in a single figure to better appreciate the variability recovered from the data. The spectra represented in black correspond to signatures corrupted by outliers, while those giv en in green represent endmembers which have been split into several components by the associated estimation procedure. Λ − 1 r = " 1 ξ + X n,t a 2 r,n,t σ 2 t i I R − 1 (68) µ ( E ) r = Λ r U T  X t,n 1 σ 2 t  y n,t − x n,t − dM t a n,t − ˇ y a r,n,t − X j 6 = r a j,n,t m j  a r,n,t  . R E F E R E N C E S [1] J. M. Bioucas-Dias, A. Plaza, N. Dobigeon, M. Parente, Q. Du, P . Gader , and J. Chanussot, “Hyperspectral unmixing overview: Geometrical, statistical, and sparse regression-based approaches, ” IEEE J . Sel. T opics Appl. Earth Observ . in Remote Sens. , v ol. 5, no. 2, pp. 354–379, Apr. 2012. [2] N. Dobigeon, J.-Y . T ourneret, C. Richard, J. C. M. Bermudez, S. McLaughlin, and A. O. Hero, “Nonlinear unmixing of hyperspectral images: Models and algorithms, ” IEEE Signal Process. Mag. , vol. 31, no. 1, pp. 89–94, Jan. 2014. t = 1 VCA t = 2 t = 3 t = 4 t = 5 t = 6 0 0.5 1 SISAL 0 0.5 1 RLMM 0 0.5 1 OU 0 0.5 1 MCMC 0 0.5 1 Fig. 9. Soil abundance map recovered by the different methods (in each row) at each time instant (given in column) for the experiment on the real dataset [the dif ferent rows correspond to VCA/FCLS, SISAL/FCLS, RLMM, OU, and the proposed method]. The images delineated in red suggest that some of the methods are particularly sensitive to the presence of outliers. t = 1 VCA t = 2 t = 3 t = 4 t = 5 t = 6 0 0.5 1 SISAL 0 0.5 1 RLMM 0 0.5 1 OU 0 0.5 1 MCMC 0 0.5 1 Fig. 10. W ater abundance map recov ered by the different methods (in each row) at each time instant (gi ven in column) for the experiment on the real dataset [the different ro ws correspond to VCA/FCLS, SISAL/FCLS, RLMM, OU, and the proposed method]. On the one hand, the images delineated in red suggest that some of the methods are particularly sensitive to the presence of outliers. On the other hand, the images delineated in green represent the abundance maps associated with signatures which ha ve been split into two components by the corresponding unmixing procedures. [3] R. Heylen, M. Parente, and P . Gader , “ A review of nonlinear hyper- spectral unmixing methods, ” IEEE J. Sel. T opics Appl. Earth Observ . in Remote Sens. , vol. 7, no. 6, pp. 1844–1868, Jun. 2014. [4] B. Somers, G. P . Asner, L. Tits, and P . Coppin, “Endmember variability in spectral mixture analysis: A revie w , ” Remote Sens. En vir onment , vol. 115, no. 7, pp. 1603–1616, Jul. 2011. [5] A. Zare and K. C. Ho, “Endmember variability in hyperspectral im- agery , ” IEEE Signal Pr ocess. Mag. , vol. 31, no. 1, pp. 95–104, Jan. 2014. [6] B. Somers, M. Zortea, A. Plaza, and G. Asner, “ Automated extraction of image-based endmember bundles for impro ved spectral unmixing, ” IEEE J. Sel. T opics Appl. Earth Observ . in Remote Sens. , vol. 5, no. 2, pp. 396–408, Apr . 2012. [7] M. A. V eganzones, L. Drumetz, G. T ochon, M. D. Mura, A. Plaza, J. M. Bioucas-Dias, and J. Chanussot, “ A new extended linear mixing model 13 t = 1 VCA t = 2 t = 3 t = 4 t = 5 t = 6 0 0.5 1 SISAL 0 0.5 1 RLMM 0 0.5 1 OU 0 0.5 1 MCMC 0 0.5 1 Fig. 11. V egetation abundance map reco vered by the different methods (in each row) at each time instant (given in column) for the experiment on the real dataset [the different ro ws correspond to VCA/FCLS, SISAL/FCLS, RLMM, OU, and the proposed method]. The images delineated in red suggest that some of the methods are particularly sensitiv e to the presence of outliers. RLMM 0 0.5 1 MCMC 0 0.5 1 t = 1 Label maps t = 2 t = 3 t = 4 t = 5 t = 6 Fig. 12. mMAP estimates of the label maps recov ered by the proposed approach, displayed at each time instant (the different rows correspond to: the estimated label map (pixels detected as outliers appear in white), the outlier energy map re-scaled in the interv al [0 , 1] obtained by the proposed method, and by RLMM). t = 1 t = 2 t = 3 t = 4 t = 5 t = 6 Fig. 13. Non-linearity maps estimated by the detector [ 55 ] applied to each image with the SISAL-extracted endmembers, with a probability of false alarm of 10 − 3 (pixels detected as non-linearities appear in white). to address spectral variability , ” in Proc. IEEE GRSS W orkshop Hyper- spectral Ima ge Signal Pr ocess.: Evolution in Remote Sens. (WHISPERS) , Lausanne, Switzerland, Jun. 2014. [8] T . Uezato, R. J. Murphy , A. Melkumyan, and A. Chlingaryan, “ A novel spectral unmixing method incorporating spectral variability within endmember classes, ” IEEE T rans. Geosci. Remote Sens. , vol. 54, no. 5, pp. 2812–2831, May 2016. [9] P .-A. Thouvenin, N. Dobigeon, and J.-Y . T ourneret, “Hyperspectral unmixing with spectral variability using a perturbed linear mixing model, ” IEEE T rans. Signal Process. , vol. 64, no. 2, pp. 525–538, Jan. 2016. [10] 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 Process. , vol. 19, no. 6, pp. 1403–1413, Jun. 2010. [11] X. Du, A. Zare, P . Gader , and D. Dranishniko v , “Spatial and spectral unmixing using the beta compositional model, ” IEEE J. Sel. T opics Appl. Earth Observ . in Remote Sens. , vol. 7, no. 6, pp. 1994–2003, Jun. 2014. [12] A. Halimi, N. Dobigeon, and J.-Y . T ourneret, “Unsupervised unmixing of hyperspectral images accounting for endmember variability , ” IEEE T rans. Image Process. , vol. 24, no. 12, pp. 4904–4917, Dec. 2015. [13] A. Halimi, N. Dobigeon, J.-Y . T ourneret, S. McLaughlin, and P . Honeine, “Unmixing hyperspectral images accounting for temporal and spatial endmember variability , ” in Pr oc. Eur opean Signal Process. Conf. (EU- SIPCO) , Nice, France, Sep. 2015, pp. 1686–1690. [14] P .-A. Thouvenin, N. Dobigeon, and J.-Y . T ourneret, “Online unmixing of multitemporal hyperspectral images accounting for spectral variability , ” IEEE T rans. Image Pr ocess. , vol. 25, no. 9, pp. 3979–3990, Sep. 2016. [15] A. Halimi, C. Mailhes, and J.-Y . T ourneret, “Nonlinear regression using smooth Bayesian estimation, ” in Pr oc. IEEE Int. Conf. Acoust., Speech, and Signal Pr ocessing (ICASSP) , Brisbane, Australia, Apr . 2015, pp. 2634–2638. [16] A. Halimi, C. Mailhes, J.-Y . T ourneret, and H. Snoussi, “Bayesian estimation of smooth altimetric parameters: Application to conv entional and delay/Doppler altimetry , ” IEEE T rans. Geosci. Remote Sens. , v ol. 54, no. 4, pp. 2207–2219, Apr . 2016. [17] S. Henrot, J. Chanussot, and C. Jutten, “Dynamical spectral unmixing of multitemporal hyperspectral images, ” IEEE Tr ans. Image Pr ocess. , vol. 25, no. 7, pp. 3219–3232, Jul. 2016. [18] A. Ert ¨ urk and A. Plaza, “Informative Change Detection by Unmixing for Hyperspectral Images, ” IEEE Geosci. Remote Sens. Lett. , vol. 12, no. 6, pp. 1252–1256, Jun. 2015. [19] S. Liu, L. Bruzzone, F . Bov olo, and P . Du, “Unsupervised multitemporal spectral unmixing for detecting multiple changes in hyperspectral im- ages, ” IEEE T rans. Geosci. Remote Sens. , vol. 54, no. 5, pp. 2733–2748, May 2016. [20] Y . Altmann, S. McLaughlin, and A. O. Hero, “Robust linear spectral unmixing using anomaly detection, ” IEEE T rans. Comput. Imag. , vol. 1, no. 2, pp. 74–85, Jun. 2015. [21] C. Chenot, J. Bobin, and J. Rapin, “Robust sparse blind source separa- tion, ” IEEE Signal Process. Lett. , vol. 22, no. 11, pp. 2172–2176, Nov . 2015. [22] J. M. Nascimento and J. M. Bioucas-Dias, “V ertex component analysis: a fast algorithm to unmix hyperspectral data, ” IEEE T rans. Geosci. Remote Sens. , vol. 43, no. 4, pp. 898–910, Apr . 2005. [23] J. M. Bioucas-Dias and M. A. T . Figueiredo, “ Alternating direction al- gorithms for constrained sparse regression: Application to hyperspectral unmixing, ” in Pr oc. IEEE GRSS W orkshop Hyperspectral Image Signal Pr ocess.: Evolution in Remote Sens. (WHISPERS) , Reykjavik, Iceland, Jun. 2010. [24] J. M. Bioucas-Dias, “ A v ariable splitting augmented Lagrangian ap- proach to linear spectral unmixing, ” in Pr oc. IEEE GRSS W orkshop Hyperspectral Ima ge Signal Pr ocess.: Evolution in Remote Sens. (WHIS- PERS) , Grenoble, France, Aug. 2009. [25] C. F ´ evotte and N. Dobigeon, “Nonlinear hyperspectral unmixing with robust nonnegativ e matrix factorization, ” IEEE T rans. Image Pr ocess. , vol. 24, no. 12, pp. 4904–4917, Dec. 2015. [26] G. H. Golub and C. F . V . Loan, “ An analysis of the total least squares problem, ” SIAM J . Numer . Anal. , vol. 17, no. 6, pp. 883 – 893, Dec. 1980. [27] R. Heylen and P . Scheunders, “ A multilinear mixing model for nonlinear spectral unmixing, ” IEEE T rans. Geosci. Remote Sens. , vol. 54, no. 1, pp. 240–251, 2016. [28] N. Dobigeon, J.-Y . T ourneret, and A. O. Hero III, “Bayesian linear unmixing of hyperspectral images corrupted by colored Gaussian noise with unknown co variance matrix, ” in Pr oc. IEEE Int. Conf . Acoust., Speech, and Signal Pr ocessing (ICASSP) , Las V egas, USA, Mar . 2008, pp. 3433–3436. [29] V . Mazet, S. Faisan, S. A w ali, M.-A. Gav eau, and L. Poisson, “Unsu- pervised joint decomposition of a spectroscopic signal sequence, ” Signal Pr ocess. , vol. 109, pp. 193–205, Apr. 2015. [30] R. Tibshirani, “Regression shrinkage and selection via the LASSO, ” J. Roy . Stat. Soc. Ser . B , vol. 58, no. 1, pp. 267–288, 1996. [31] J. P . V ila and P . Schniter, “Expectation-Maximization Gaussian-mixture approximate message passing, ” IEEE T rans. Signal Pr ocess. , v ol. 61, no. 19, pp. 4658–4672, Oct. 2013. [32] N. Dobigeon, A. O. Hero, and J.-Y . T ourneret, “Hierarchical Bayesian sparse image reconstruction with application to MRFM, ” IEEE T rans. Image Pr ocess. , vol. 18, no. 9, pp. 2059–2070, Sep. 2009. [33] J. J. Kormylo and J. M. Mendel, “Maximum likelihood detection and estimation of Bernoulli-Gaussian processes, ” IEEE T rans. Inf . Theory , vol. 28, no. 3, pp. 482–488, May 1982. 14 [34] M. Lavielle, “Bayesian deconvolution of Bernoulli-Gaussian processes, ” Signal Process. , vol. 33, no. 1, pp. 67–79, Jul. 1993. [35] S. Bourguignon and H. Carfantan, “Bernoulli-Gaussian spectral analysis of une venly spaced astrophysical data, ” in Pr oc. IEEE-SP W orkshop Stat. and Signal Pr ocessing (SSP) , Bordeaux, France, Jul. 2005, pp. 811–816. [36] C. Bazot, N. Dobigeon, and J.-Y . T ourneret, “Bernoulli-Gaussian model for gene expression analysis, ” in Proc. IEEE Int. Conf. Acoust., Speech, and Signal Pr ocessing (ICASSP) , Prague, Czech Republic, May 2011, pp. 5996–5999. [37] L. Chaari, J.-Y . T ourneret, and C. Chaux, “Sparse signal recovery using a Bernoulli generalized Gaussian prior, ” in Proc. Eur opean Signal Process. Conf. (EUSIPCO) , Nice, France, Sep. 2015, pp. 1711–1715. [38] O. Eches, N. Dobigeon, and J. Y . T ourneret, “Enhancing hyperspectral image unmixing with spatial correlations, ” IEEE T rans. Geosci. Remote Sens. , vol. 49, no. 11, pp. 4239–4247, Nov . 2011. [39] 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 Tr ans. Image Pr ocess. , vol. 23, no. 5, pp. 2148–2158, Jun. 2014. [40] J.-F . Giovannelli, “Estimation of the Ising field parameter thanks to the exact partition function, ” in Proc. IEEE Int. Conf. Ima ge Pr ocessing (ICIP) , Hong-Kong, China, Sep. 2010, pp. 1441–1444. [41] 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 Tr ans. Signal Pr ocess. , vol. 57, no. 11, pp. 4355–4368, Nov . 2009. [42] E. J. Cand ` es, X. Li, Y . Ma, and J. Wright, “Rob ust principal component analysis?” journal of ACM , vol. 58, no. 1, pp. 1–37, 2009. [43] J.-F . Giovannelli, “Ising field parameter estimation from incomplete and noisy data, ” in Pr oc. IEEE Int. Conf. Image Pr ocessing (ICIP) , Brussels, Belgium, Sep. 2011, pp. 1853–1856. [44] A. Doucet, S. J. Godsill, and C. P . Robert, “Marginal maximum a posteriori estimation using Markov chain Monte Carlo, ” Stat. Comput. , vol. 12, no. 1, pp. 77–84, Jan. 2002. [45] N. Dobigeon and J.-Y . T ourneret, “Efficient sampling according to a multiv ariate Gaussian distribution truncated on a simplex, ” IRIT/ENSEEIHT/T ´ eSA, France, T ech. Rep., Mar . 2007. [Online]. A vailable: http://www .enseeiht.fr/ ∼ dobigeon/papers/ Dobigeon T echReport 2007b .pdf [46] Y . Altmann, S. McLaughlin, and N. Dobigeon, “Sampling from a multiv ariate Gaussian distrib ution truncated on a simplex: a revie w , ” in Pr oc. IEEE-SP W orkshop Stat. and Signal Pr ocessing (SSP) , Gold Coast, Australia, Jul. 2014, pp. 113–116. [47] Y . Altmann, N. Dobigeon, and J.-Y . T ourneret, “Unsupervised post- nonlinear unmixing of hyperspectral images using a Hamiltonian Monte Carlo algorithm, ” IEEE T rans. Image Pr ocess. , vol. 23, no. 6, pp. 2663– 2675, Jun. 2014. [48] A. Pakman and L. Paninski, “Exact Hamiltonian Monte Carlo for trun- cated multiv ariate Gaussians, ” journal of Computational and Graphical Statistics , vol. 23, no. 2, pp. 518–542, 2014. [49] Z. I. Bote v , “The normal law under linear restrictions: Simulation and estimation via minimax tilting, ” J. Roy . Stat. Soc. Ser . B , vol. 78, no. 3, 2016, to appear . [50] N. Chopin, “Fast simulation of truncated Gaussian distributions, ” Statis- tics and Computing , vol. 21, no. 2, pp. 275–288, 2011. [51] L. Onsager, “ A two-dimensional model with an order-disorder transi- tion, ” Phys. Rev . , vol. 65, no. 3 & 4, pp. 117–149, Feb. 1944. [52] P .-A. Thouvenin, N. Dobigeon, and J.-Y . T ourneret, “ A hierarchical Bayesian model accounting for endmember variability and abrupt spectral changes to unmix multitemporal hyperspectral images – complementary results and supporting material, ” Univ ersity of T oulouse, IRIT/INP-ENSEEIHT , T ech. Rep., May 2017. [Online]. A vailable: http://thouvenin.perso.enseeiht.fr/papers/Thouvenin2017TR.pdf [53] A. Halimi, P . Honeine, M. Kharouf, C. Richard, and J.-Y . T ourneret, “Estimating the intrinsic dimension of hyperspectral images using a noise-whitened eigen-gap approach, ” IEEE T rans. Geosci. Remote Sens. , vol. 54, no. 7, pp. 3811–3821, Jul. 2016. [54] P .-A. Thouvenin, N. Dobigeon, and J.-Y . T ourneret, “Online unmixing of multitemporal hyperspectral images accounting for spectral variability - complementary results and supporting material, ” Univ ersity of T oulouse, IRIT/INP-ENSEEIHT , T ech. Rep., Oct. 2015. [Online]. A vailable: http://thouvenin.perso.enseeiht.fr/papers/Thouvenin TR 2015.pdf [55] Y . Altmann, N. Dobigeon, J.-Y . T ourneret, and J. C. M. Bermudez, “ A robust test for nonlinear mixture detection in hyperspectral images, ” in Pr oc. IEEE Int. Conf. Acoust., Speech, and Signal Pr ocessing (ICASSP) , V ancouver , Canada, Jun. 2013, pp. 2149–2153. Pierre-Antoine Thouvenin (S’15–M’17) received the state engineering de- gree in electrical engineering from ENSEEIHT , T oulouse, France, and the M.Sc. degree in signal processing from the National Polytechnic Institute of T oulouse (INP T oulouse), both in 2014. At the time of writing this paper, he was working towards the Ph.D. degree within the Signal and Communications Group of the IRIT Laboratory , T oulouse, France. His research interests include statistical modeling, optimization techniques and hyperspectral unmixing. Nicolas Dobigeon (S’05–M’08–SM’13) received the state engineering degree in electrical engineering from ENSEEIHT , T oulouse, France, and the M.Sc. degree in signal processing from the National Polytechnic Institute of T oulouse (INP T oulouse), both in June 2004, as well as the Ph.D. degree and Habilita- tion ` a Diriger des Recherches in Signal Processing from the INP T oulouse in 2007 and 2012, respectively . He was a Post-Doctoral Research Associate with the Department of Electrical Engineering and Computer Science, Uni versity of Michigan, Ann Arbor, MI, USA, from 2007 to 2008. Since 2008, he has been with the National Polytechnic Institute of T oulouse (INP-ENSEEIHT , University of T oulouse) where he is currently a Professor. He conducts his research within the Signal and Communications Group of the IRIT Laboratory and he is also an affiliated faculty member of the T elecom- munications for Space and Aeronautics (T ´ eSA) cooperative laboratory . His current research interests include statistical signal and image processing, with a particular interest in Bayesian in verse problems with applications to remote sensing, biomedical imaging and genomics. Jean-Yves T ourneret (SM’08) recei ved the ing ´ enieur degree in electrical engineering from the Ecole Nationale Sup ´ erieure d’Electronique, d’Electrotechnique, d’Informatique, d’Hydraulique et des T ´ el ´ ecommunications (ENSEEIHT) de T oulouse in 1989 and the Ph.D. degree from the National Polytechnic Institute from T oulouse in 1992. He is currently a professor in the univ ersity of T oulouse (ENSEEIHT) and a member of the IRIT laboratory (UMR 5505 of the CNRS). His research activities are centered around statistical signal and image processing with a particular interest to Bayesian and Marko v chain Monte-Carlo (MCMC) methods. He has been in volv ed in the organization of several conferences including the European conference on signal processing EUSIPCO’02 (program chair), the international conference ICASSP’06 (plenaries), the statistical signal processing workshop SSP’12 (international liaisons), the International W orkshop on Computational Advances in Multi-Sensor Adaptiv e Processing CAMSAP 2013 (local arrangements), the statistical signal processing workshop SSP’2014 (special sessions), the w orkshop on machine learning for signal processing MLSP’2014 (special sessions). He has been the general chair of the CIMI workshop on optimization and statistics in image processing hold in T oulouse in 2013 (with F . Malgouyres and D. Kouam ´ e) and of the International W orkshop on Computational Advances in Multi-Sensor Adaptive Processing CAMSAP 2015 (with P . Djuric). He has been a member of different technical committees including the Signal Processing Theory and Methods (SPTM) committee of the IEEE Signal Processing Society (2001-2007, 2010-present). He has been serving as an associate editor for the IEEE Transactions on Signal Processing (2008- 2011, 2015-present) and for the EURASIP journal on Signal Processing (2013-present).