Hyperspectral Unmixing in Presence of Endmember Variability, Nonlinearity or Mismodelling Effects

1 Hyperspectral Unmixing in Presence of Endmember V ariability , Nonlinearity or Mismodelling Ef fects Abderrahim Halimi 1 , Member , IEEE, Paul Honeine 2 , Member , IEEE , Jose Bioucas-Dias 3 , Senior Member , IEEE , Abstract This paper presents three hyperspectral mixture models jointly with Bayesian algorithms for supervised hyperspectral unmixing. Based on the residual component analysis model, the proposed general formulation assumes the linear model to be corrupted by an additi ve term whose e xpression can be adapted to account for nonlinearities (NL), endmember variability (EV), or mismodelling effects (ME). The NL effect is introduced by considering a polynomial expression that is related to bilinear models. The proposed new formulation of EV accounts for shape and scale endmember changes while enforcing a smooth spectral/spatial v ariation. The ME formulation takes into account the effect of outliers and copes with some types of EV and NL. The known constraints on the parameter of each observ ation model are modeled via suitable priors. The posterior distrib ution associated with each Bayesian model is optimized using a coordinate descent algorithm which allows the computation of the maximum a posteriori estimator of the unknown model parameters. The proposed mixture and Bayesian models and their estimation algorithms are validated on both synthetic and real images showing competitiv e results regarding the quality of the inferences and the computational complexity when compared to the state-of-the-art algorithms. Index T erms (1) A. Halimi is with the School of Engineering and Ph ysical Sciences, Heriot-W att Uni versity , Edinbur gh U.K. (e-mail: a.halimi@hw .ac.uk). (2) P . Honeine is with the Normandie Univ , UNIROUEN, UNIHA VRE, INSA Rouen, LITIS, Rouen, France (e-mail: paul.honeine@univ-rouen.fr). (3) J. M. Bioucas-Dias is with the Instituto de T elecomunicações and Instituto Superior Técnico, Uni versidade de Lisboa, Portugal (bioucas@lx.it.pt). This work was supported in part by the HYP ANEMA ANR Project under Grant ANR-12-BS03-003, and in part by the Portuguese Fundacão para a Ciência e T ecnologia (FCT), grant UID/EEA/5008/2013 July 20, 2016 DRAFT 2 Hyperspectral imagery , endmember v ariability , nonlinear spectral unmixing, robust unmixing, mismodelling effect, Bayesian estimation, coordinate descent algorithm, Gaussian process, Gamma Markov random field I . I N T RO D U C T I O N Spectral unmixing (SU) of hyperspectral images has been the subject of intensiv e interest ov er the last two decades. This is a source separation problem consisting of recov ering the spectral signatures (endmembers) of the materials present in the scene, and quantifying their proportions within each hyperspectral image pixel. The linear mixture model (LMM) is the widely used model for SU mainly because of its simplicity . Ho we ver , this model can be inappropriate for some hyperspectral scenarios, namely if there are v olumetric scattering, or terrain relief, or intimate mixtures of materials [1]. Nonlinear mixture models (NLMM) appear then as an alternati v e to better account for those ef fects [2], [3]. There e xists two main families for NLMMs: the first family is signal processing based and seeks to construct flexible models that can represent a wide range of nonlinearities. The second family is the physical based and includes the intimate mixture models [1] and those accounting for multiple scattering such as polynomial [4] or bilinear models [5]–[9]. This paper considers physical based nonlinearity while focusing on polynomial/bilinear models. In both LMM and NLMM, the endmembers are generally assumed fix ed in the whole image. This appears as a clear simplification since in many cases, the endmember spectra vary along the image causing what is kno wn as spectral variability or endmember v ariability (EV) [10], [11]. Spectral variability has been identified as a relev ant source of error in abundance estimation and is attracting growing interest in the hyperspectral community [10]–[12]. Many algorithms have been proposed in the literature to describe this variability and they can be gathered into two main classes. The first class considers each physical material as a set or bundle of spectra that are assumed known [13], [14] or estimated from the data [15], [16]. In this class, we find parametric representation of the endmembers such as the multiplication of each endmember by a pixel dependent constant to account for a change of illumination in the observed scene [17]–[19]. The second class of methods relies on a statistical representation of the endmembers that are assumed to be random vectors with giv en probability distributions. T wo main statistical models of the endmembers hav e been considered in the literature. The normal compositional model (NCM) assuming Gaussian distrib utions for the endmembers [20]–[22] and the Beta compositional model [23] that exploits the physically realistic range July 20, 2016 DRAFT 3 of the endmember reflectances by assigning them a Beta distribution. This paper introduces a spectral smooth EV that can be expressed by a Gaussian distributed endmembers resulting in a model that is closely related to NCM. Mismodelling ef fects (ME) are also a concern when processing hyperspectral images. These effects can be due to the presence of some physical phenomena such as nonlinearity or endmember v ariability . Ho we ver , they can also be due to the propagated errors in the signal processing chain. Indeed, SU is generally performed using three steps: (i) estimating the number of endmembers, (ii) identifying the endmembers using an endmember extraction algorithm (EEA) such as verte x component analysis (VCA) [24], and N-FINDR [25] and (iii) estimating the ab undances under physical non-negati vity and sum-to-one constraints using algorithms such as the fully constraine d least squares [26]. Man y studies consider a supervised unmixing scenario which aims at estimating the abundances while assuming that the two first unmixing steps were successfully implemented [26]–[29]. Howe ver , an error on the estimated number of endmembers or in their spectra may lead to bad abundance estimates. This problem is considered by some ne w rob ust unmixing algorithms that aim at reducing the effect of outliers or mismodelling effects [30], [31]. In this paper, the ME is considered by reducing the effect of smooth spatial/spectral outliers that can be due to the presence of NL, EV or to signal processing chain errors. The first contrib ution of this paper is the introduction of a general formulation for the mixture model to account for three phenomena: (i) nonlinearity , (ii) endmember variability or (iii) mismodelling ef fects. Based on the residual component analysis model [32], the proposed general formulation assumes the linear model to be corrupted by an additi ve term whose expression can be adapted to account for the studied phenomenon. This residual term is expressed as a combination of the abundances or the endmembers depending on the studied EV or NL effects. Indeed, the first model studies NL ef fect which can be modeled by considering a modification to the polynomial term proposed in [28], that depends only on the endmember spectra. The second model considers EV effect by introducing a smooth additiv e de viation of each endmember from kno wn spectra. Thus, contrary to the NCM described in [12], [20], [21], the proposed model takes into account the smooth spectral and spatial v ariation of the endmembers in presence of EV . Moreov er , and to account for ME, a third model is introduced for the residual term while accounting for the spatial and spectral smooth properties of the corrupting term. The proposed formulation is therefore general, cov ers many physical phenomena and can be related to many NL [4], [5], [7]–[9], [33] and EV models July 20, 2016 DRAFT 4 [12], [18], [20], [21]. The second contrib ution of this paper is the introduction of three hierarchical Bayesian models associated with each observ ation model. These hierarchical models introduce prior distributions that enforce kno wn physical constraints on the estimated parameters such as the sum-to-one and positivity of the abundances, positi vity of the nonlinear coef ficients and the smooth spectral behavior of EV and ME. Moreo ver , the spatial correlation of the residual term has been introduced by considering Mark ov random fields [34], [35]. Using the likelihood and the considered prior distrib utions, the joint posterior distrib ution of the unkno wn parameter vector is then deriv ed for each model. The minimum mean square error (MMSE) and maximum a posteriori (MAP) estimators of these parameters cannot be easily computed from the obtained joint posteriors. In this paper , the MAP estimator is e v aluated by considering a coordinate descent algorithm (CD A) [36]–[38] that sequentially updates the abundances, the noise variances and the residual terms. The proposed Bayesian models and estimation algorithms are validated using synthetic and real hyperspectral images. The obtained results are very promising and sho w the potential of the proposed mixture and Bayesian models and their associated inference algorithms. The remainder of this paper in organized as follo ws. Section II introduces the proposed gen- eral formulation for the mixture model and its variants to deal with NL, EV and ME ef fects. The proposed hierarchical Bayesian models and their estimation algorithms are introduced in Sections III and IV. Section V is dev oted to testing and v alidating the proposed techniques using synthetic images with known ground truth. Section VI shows results obtained using real hyperspectral images. Conclusions and future work are finally 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: July 20, 2016 DRAFT 5 N number of pixels R number of endmembers L number of spectral bands D number of the nonlinear coef ficients y i,j ∈ R L × 1 pixel located at the i th row and j th column a i,j ∈ R R × 1 abundance vector of the pixel ( i, j ) e i,j ∈ R L × 1 noise vector of the pixel ( i, j ) γ i,j ∈ R D × 1 nonlinear coefficients of the pixel ( i, j ) c i,j ∈ R illumination coefficient of the pixel ( i, j ) m r ∈ R L × 1 spectrum of the r th fixed endmember c i,j m r ∈ R L × 1 pixel dependent spectrum of the r th endmember k r,i,j ∈ R L × 1 smooth vector for endmember r d i,j ∈ R L × 1 smooth outlier vector φ i,j ∈ R L × 1 residual term of the pixel ( i, j ) Σ ∈ R L × L noise cov ariance matrix σ 2 ∈ R L × 1 diagonal of the noise cov ariance matrix M ∈ R L × R endmember matrix S i,j ∈ R L × R endmember matrix of the pixel ( i, j ) Y ∈ R L × N spectra of the pixels A ∈ R R × N abundance matrix  ∈ R 1 × N v ariance of the residual terms (or their energy) B. Mixing model: nonlinearity , endmember variability and mismodelling effects The proposed formulation is based on a residual component analysis model [32] that is expressed as the sum of a linear model and a residual term. The general observation model for the ( L × 1) pixel spectrum y i,j is giv en by y i,j = R X r =1 a r,i,j s r,i,j + φ i,j ( S i,j , a i,j ) + e i,j = S i,j a i,j + φ i,j ( S i,j , a i,j ) + e i,j , (1) where a i,j is an ( R × 1) vector of abundances associated with the pix el ( i, j ) , R is the number of endmembers, e i,j ∼ N ( 0 , Σ ) is an additi ve centered Gaussian noise with a diagonal cov ariance matrix Σ = diag ( σ 2 ) , σ 2 = ( σ 2 1 , · · · , σ 2 L ) T is an ( L × 1) vector containing the noise variances, L is the number of spectral bands, S i,j ( M ) = S i,j is the endmember matrix July 20, 2016 DRAFT 6 that depends on each pixel to introduce EV effect, M is a fixed endmember matrix that is assumed known (extracted using an EEA) and φ i,j ( S i,j , a i,j ) is a residual term that might depend on the endmembers or the abundances to model NL or EV effects. Due to physical constraints, the abundance vector a i,j = ( a r,i,j , · · · , a R,i,j ) T satisfies the follo wing positivity and sum-to-one (PSTO) constraints a r,i,j ≥ 0 , ∀ r ∈ { 1 , . . . , R } and R X r =1 a r,i,j = 1 . (2) Eq. (1) shows a general model that can be adapted to account for different physical phenom- ena. In this paper , we consider three variants dealing with NL, EV or ME. The NL model is designed to deal with the multiple scattering effect that appears in presence of terrain relief, and/or trees. The EV model accounts for the deviation of the endmembers that is commonly observed in presence of v egetation (such as trees or grass), and shadow . It is common to observe the NL and the EV effects simultaneously when analyzing a scene. Therefore, a ME model has been proposed to account for both effects. The next sections provide details regarding each of these models. 1) Nonlinearity effect: Nonlinear mixing models provide a useful alternati ve for overcom- ing the inherent limitations of the LMM. Indeed, LMM can be inappropriate for some hyper- spectral images, such as those containing trees, vegetation or urban areas. Bilinear/polynomial models hav e sho wn useful results for these scenes by addressing the problem of double scat- tering effects. In addition to the LMM terms, these models consider second order interactions between endmembers and neglects the ef fect of the higher order terms [5], [7], [8]. This paper considers the following polynomial/bilinear nonlinear model y i,j = c i,j M a i,j + φ N L i,j ( M ) + e i,j (3) where the residual component are similar to [28] as follows φ N L i,j ( M ) = c 2 i,j R − 1 X k =1 R X k 0 = k +1 γ ( k,k 0 ) i,j √ 2 m k  m k 0 + R X k =1 γ ( k ) i,j m k  m k ! , (4) with γ i,j =  γ (1) i,j , · · · , γ ( R ) i,j , γ (1 , 2) i,j , · · · , γ ( R − 1 ,R ) i,j  T , ∀ i, j is the ( D × 1) vector of positi ve nonlinear coefficients 1 , D = R ( R +1) 2 ,  denotes the Hadamard (termwise) product, s r,i,j = 1 The nonlinear coefficients represent the additional amount of bilinear interactions between the endmembers, thus, they should be positiv e [2], [3]. July 20, 2016 DRAFT 7 c i,j m r , ∀ i, j , with c i,j a pix el dependent illumination coef ficient. The model (3) generalizes the model [28] by including an EV illumination parameter c i,j that accounts for the main spectral v ariation of endmembers as shown in [18], [19]. Contrary to the RCA model [33], model (3) considers physical positiv e γ i,j (the RCA model [33] can be obtained by marginalizing unconstrained γ i,j as shown in [28]). Note also that (3) generalizes the LMM (obtained when γ i,j = 0 , and c i,j = 1 ∀ , i, j ) and has a polynomial-like form as for the bilinear models (GBM [5], PPNMM [4], Nascimento [7], Fan [8] and Meganem [9] models). Note finally that model (3) (with no illumination v ariation) has been studied in [28] when considering a Marko v chain Monte-Carlo (MCMC) approach and have shown good performance for processing hyper- spectral images. Ho we ver , the MCMC estimation algorithm was computationally expensiv e, and we consider in this paper a faster algorithm based on a coordinate descent algorithm. 2) Endmember variability ef fect: Due to the lo w spatial resolution of h yperspectral images, the image might represent very large scenes. Therefore, it makes sense to assume that the same material (such as ve getation) might differ with respect to (w .r .t.) the image regions resulting in what it is known as EV . This variability introduces a modification in the shape and the scale of the endmembers spectrum in each pixel, i.e., s i,j depends on the pixel location. As an e xample, Fig. 1 (top) presents spectra associated with the T opaz ph ysical element and extracted from the USGS library . Even though these spectra are associated with the same material, they show some differences which is known as EV ef fect. T o highlight this effect, we compute the average spectrum in each spectral band, and assume that EV is obtained by computing the difference between the spectra of Fig. 1 (top) and the av erage spectrum. An example of the obtained differences is sho wn in Fig. 1 (bottom). This figure clearly sho ws that the difference between the spectra (which represents the EV effect) can be approximated by smooth functions (see the dashed lines). Therefore, to account for the shape variability of each endmember , each endmember can be approximated by the sum of a fixed spectrum and a smooth spectral function representing EV . This smooth function can be modeled by a parametric approach such as spline, or a statistical approach as Gaussian process (which will be studied in section III-B3). In this paper , we consider the follo wing model for endmember variability s r,i,j = m r + k r,i,j , (5) where k r,i,j is a smooth spectral function. This model leads to y i,j = M a i,j + φ E V i,j ( a i,j ) + e i,j , (6) July 20, 2016 DRAFT 8 where φ E V i,j ( a i,j ) = R X r =1 a r,i,j k r,i,j . (7) Note that (7) does not consider an illumination parameter c i,j since its effect can be included in the smooth function k r,i,j . Model (7) relates to state-of-the-art models as follows: (i) it generalizes the LMM that can be retrie ved for k r,i,j = 0 L , ∀ i, j , (ii) it generalizes the model [18] by including the effect of shape variability , and (iii) it has a similar formulation as in [39] while accounting for the spectral smoothness of the residuals. Note finally that in the special case where k r,i,j is Gaussian distributed (Gaussian process), the model (6) will improv e the GNCM introduced in [12] by including the smooth behavior of EV (model (6) is closely related to the GNCM that can be obtained by marginalizing s i,j when considering a Bayesian approach). 3) Mismodelling effects (ME) or outliers: The linear model is widely used because of its simplicity . As previously shown, there is a lot of situations where the linear model is not v alid because of the presence of variability , nonlinearity or mismodelling ef fects due to the signal processing chain errors. This section accounts for mismodelling effects or the presence of outliers by considering a residual term that sho ws spatial and spectral correlations. The observ ation model is giv en by y i,j = c i,j M a i,j + φ M E i,j + e i,j , with φ M E i,j = d i,j , (8) where d i,j is a smooth spectral function. Similarly to the previous models, (8) reduces to the LMM when d i,j = 0 L , and c i,j = 1 , ∀ i, j . Note that other models have been introduced in the literature to account for the effect of outliers such as [30] that proposed spatial/spectral correlated outliers by considering discrete MRF and [31] which proposed a positiv e sparse outliers that has no spatial or spectral correlation. Note that (8) can be seen as a special case of the EV model (7) when k r,i,j = k r 0 ,i,j , ∀ r , r 0 and c i,j = 1 , i.e., the same variability is af fecting the dif ferent endmembers. Note also that the NL model (4) reduces to (8) if γ ( k,k 0 ) = γ , ∀ k , k 0 since the spectra m k  m k 0 , ∀ , k , k 0 are generally smooth. Ho we ver , (8) is more flexible since it does not consider the positi vity constraint (mainly to account for EV). The next section introduces the proposed hierarchical Bayesian models associated with each observ ation model. July 20, 2016 DRAFT 9 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 the unkno wn parameters of models (3), (6) and (8) which consider dif ferent residual terms. The common unknown parameters of these models include the ( L × 1 ) diagonal cov ariance matrix of the noise denoted by Σ = diag ( σ 2 ) , the ( R × N ) abundance matrix A , the ( N × 1) vector of illumination v ariability c for the NL (3) and ME (8) models, and finally the (1 × N ) v ariance of the residual terms  (or their energy). The remaining parameters depend on the considered residual term as follows: • the ( D × N ) matrix of nonlinear coefficients Γ for (3) • the ( R × L × N ) matrix of endmember variability coef ficients K for (6) • the ( L × N ) matrix of mismodelling coefficients D for (8). T able I presents the parameters (and their dimensions) related to each algorithm. It also shows when a spatial correlation (sc) or a spectral smoothing (ss) is introduced. A. Likelihood Using the observation model (1), the Gaussian properties of the noise sequence e n , and exploiting independence between the observ ations in different spectral bands, yield the fol- lo wing Gaussian distribution for the likelihood f ( y i,j | a i,j , Σ , S i,j , φ i,j ) ∝ 1 Q L l =1 σ 2 l ! 1 2 × exp ( −  y i,j − S i,j a i,j − φ i,j  T Σ − 1  y i,j − S i,j a i,j − φ i,j  2 ) . (9) The joint likelihood of the observation matrix Y is obtained by considering the independence between the observed pixels as follo ws f ( Y | A , Σ , S , Φ ) ∝ Y i Y j f ( y i,j | a i,j , Σ , S i,j , φ i,j ) . (10) B. P arameter priors This section introduces the prior distributions that we have chosen for the parameters of interest A , c , Γ , K , D , Σ and  . July 20, 2016 DRAFT 10 1) Abundance matrix A : In order to satisfy the constraints (2), the abundance vector should liv e in the following simplex S = ( a i,j   a r,i,j ≥ 0 , ∀ r and R X r =1 a r,i,j = 1 ) . (11) Since there is no additional information about this parameter vector , we propose to assign a uniform prior in the simplex S to the abundances [5], [40]. 2) Prior for c : The variation in illumination is introduced by the parameter c i,j that is pixel-dependent. In man y works, this parameter has been fix ed to the v alue 1 , which represents the absence of illumination variability [5], [26], [27]. In this paper , we allow this parameter to fluctuate around the v alue 1 to include the effect of illumination. Therefore, we assign the follo wing conjugate Gaussian distribution as a prior for c i,j c i,j ∼ N  1 , η 2  , (12) where ∼ means “is distributed according to”, η 2 is a fixed variance that ensures the value of c i,j to be located around 1 , thus, av oiding negati v e values for this parameter in practice ( η 2 = 0 . 01 in the rest of the paper). For simplicity , we denote “ x | θ ∼ ... ”, by “ x ∼ ... ” when the parameter θ is a user fixed parameter . Note finally that the joint prior of c is obtained by assuming a priori independence between the coefficients c i,j , as follows f ( c ) = Q i,j f ( c i,j ) . 3) Prior for the residual terms: a) Nonlinear coef ficients γ i,j : Due to physical constraints, the nonlinear coefficients should satisfy the positi vity constraint. Similarly to [28], γ i,j are assigned the follo wing truncated Gaussian prior γ i,j |  2 i,j ∼ N ( R +) D  0 D ,  2 i,j I D  , (13) where I D denotes the D × D identity matrix and  2 i,j is a variance parameter that is pixel dependent. From (13), it is clear that this variance is related to the strength of the nonlinearities at the pix el ( i, j ) (via the norm || γ i,j || 2 ). Moreov er , as in [28], we assume the nonlinear energies to v ary smoothly from one pixel to another which will be introduced by considering a specific prior for  2 i,j , as explained in Section III-B5. Note finally that the joint prior of Γ is obtained by assuming a priori independence between the nonlinear coefficients, as follows f ( Γ |  ) = Q i,j f  γ i,j |  2 i,j  . July 20, 2016 DRAFT 11 b) EV coefficients k r,i,j : As pre viously explained, k r,i,j is a smooth spectral vector . This property can be satisfied by considering a parametric expression for k r,i,j (such as spline), or a statistical approach such as Gaussian process. In this paper , we assign a Gaussian Markov random field prior [34] for K r = { k r,i,j , ∀ i, j } that ensures both spectral smoothness and spatial correlation between the residuals. The prior is gi ven by K r ∝ exp " − X i,j P ( i 0 ,j 0 ) ∈ ν ( i,j ) || k r,i,j − k r,i 0 ,j 0 || 2 16 β 2 r + k T r,i,j H − 1 k r,i,j 2 α 2 r !# , (14) where || · || denotes the standard l 2 norm, ν ( i, j ) is the neighborhood of the pixel ( i, j ) (an eight neighborhood structure is considered), and H is an ( L × L ) matrix 2 representing the squared-exponential cov ariance function gi ven by H ( `, ` 0 ) = exp h − ( ` − ` 0 ) 2 ( L/ 2) 2 i , which introduces the spectral smoothness on k r,i,j . The conditional distribution of k r,i,j gi ven its neighbors is a conjugate Gaussian distribution that can be expressed as follows k r,i,j | k r,i 0 ,j 0 ∼  N  0 L , α 2 r H  × N  µ r,i,j , β 2 r I L  , (15) for ( i 0 , j 0 ) ∈ ν ( i, j ) and µ r,i,j = 1 / 8 P ( i 0 ,j 0 ) ∈ ν ( i,j ) k r,i 0 ,j 0 is obtained by the av erage of the neighborhood residuals of the pixel ( i, j ) . It is clear from (15) that the first normal (left) introduces the spectral smoothness and the second one (right) introduces the spatial correlations between the residuals. The hyperparameters are fixed by cross-validation, where α 2 r controls the amplitude of the spectral smooth residuals and β 2 r the degree of spatial correlation between residuals. Finally the joint prior of K is obtained by assuming a priori independence between K r , that is f ( K ) = Q R r =1 f ( K r ) . c) Mismodelling coefficients d i,j : The spectral vector d i,j should also be smooth. This property is satisfied by considering a Gaussian prior for d i,j ensuring smoothness as follows d i,j |  2 i,j ∼ N  0 L ,  2 i,j H  , (16) where H is the squared-exponential cov ariance function described in the previous section. As for the nonlinear coef ficients, spatial correlation is introduced by enforcing a smooth variation of the residual energies  d T i,j H − 1 d i,j  . This is achiev ed by considering a specific prior for  2 i,j , as explained in Section III-B5. Finally the joint prior of D is obtained by assuming a priori independence between the mismodelling coefficients f ( D |  ) = Q i,j f  d i,j |  2 i,j  . 2 When processing real images, some bands are removed because of water absorptions and other atmospheric perturbations. These bands should also be removed from the covariance matrix by removing the corresponding columns and rows. July 20, 2016 DRAFT 12 4) Noise variances: The noise variances are assigned a conjugate in verse gamma distri- bution giv en by f  σ 2  = L Y ` =1 f  σ 2 `  , with σ 2 ` ∼ I G ( ϕ ` , ψ ` ) , (17) where σ 2 ` are assumed a priori independent. The hyperparameters ϕ ` and ψ ` are fixed to approximate the Hysime estimated v ariances [41]. Indeed, Hysime estimates the noise by subtracting the spectral correlated (smooth) signal from the observed data. This is in agree- ment with the considered model (1) that represents the observ ations as a sum of a smooth part (linear mixture and residuals) and an independent Gaussian noise. Note finally that the parameters can also be set to ϕ ` = ψ ` = 0 in absence of prior knowledge about σ 2 ` , leading to a noninformativ e Jef freys’ prior . 5) Ener gy of the r esidual term: Due to the spatial organization of hyperspectral images, we expect the energies of the nonlinear coefficients γ i,j and the mismodelling coef ficients d i,j to v ary smoothly from one pixel to another . This behavior is obtained by introducing an auxiliary variable w (of size N row × N col ) and assigning a gamma Marko v random field prior (GMRF) for the couple (  , w ) giv en by (see [28], [35] for more details reg arding this prior) f (  , w | ζ ) = 1 Z ( ζ ) Q ( i,j ) ∈ ν   − 2(4 ζ +1) i,j × Q ( i 0 ,j 0 ) ∈ ν w w 2(4 ζ − 1) i 0 ,j 0 × Q (( i,j ) , ( i 0 ,j 0 )) ∈E exp  − ζ w 2 i 0 ,j 0  2 i,j  , (18) where Z ( ζ ) is a normalizing constant, the partition ν  (resp. ν w ) denotes the collection of v ariables  (resp. w ), the edge set E consists of pairs ( i, j ) representing the connection between the variables and ζ is a fixed coupling parameter that controls the amount of spatial smoothness enforced by the GMRF . This prior ensures that each  2 i,j is connected to four neighbor elements of w and vice-versa. Note that the energies  2 i,j are conditionally independent and that a 2 nd order neighbors (i.e., the spatial correlation) is introduced via the auxiliary variables w . An interesting property of this joint prior is that the conditional prior distrib utions of  and w reduce to conjugate in verse Gamma ( I G ) and gamma ( G ) distributions as follows  2 i,j | w , ζ ∼ I G (4 ζ , 4 ζ ρ 1 ,i,j ( w )) w 2 i,j |  , ζ ∼ G (4 ζ , 1 / (4 ζ ρ 2 ,i,j (  ))) , (19) where July 20, 2016 DRAFT 13 ρ 1 ,i,j ( w ) = ( w 2 i,j + w 2 i +1 ,j + w 2 i,j +1 + w 2 i +1 ,j +1 ) / 4 ρ 2 ,i,j (  ) = (  − 2 i,j +  − 2 i − 1 ,j +  − 2 i,j − 1 +  − 2 i − 1 ,j − 1 ) / 4 . (20) C. P osterior distributions The proposed Bayesian models associated with each observation model are summarized in the directed ac yclic graphs (D A Gs) displayed in Fig. 2. The parameters of interest are Θ NL = ( A , Γ , σ , c ,  , w ) , Θ EV = ( A , K , σ ) , and Θ ME = ( A , D , σ , c ,  , w ) . The joint posterior distribution of these Bayesian models can be computed from the follo wing Bayes’ rule f ( Θ | Y ) ∝ f ( Y | Θ ) f ( Θ ) , (21) with f ( Θ NL ) = f ( A ) f ( σ ) f ( c ) f ( Γ |  ) f (  , w ) f ( Θ EV ) = f ( A ) f ( σ ) f ( K ) f ( Θ ME ) = f ( A ) f ( σ ) f ( c ) f ( D |  ) f (  , w ) , (22) where ∝ means “proportional to” and we hav e assumed a priori independence between the parameters of each model. The MMSE and MAP estimators associated with the posterior (21) are not easy to determine. In this paper , and akin to [38], we propose to e v aluate the MAP estimator by using an optimization technique maximizing the posterior (21) w .r .t. the parameters of interest (or equi valently , minimizing the ne gati ve log-posterior F ( Θ ) = − log [ f ( Θ | Y )] denoted as “cost function” in the rest of this paper). I V . C O O R D I N A T E D E S C E N T A L G O R I T H M S This section describes the optimization algorithms maximizing the posteriors (21) asso- ciated with the NL, EV and ME models w .r .t. the parameters of interest. This provides the MAP estimator of Θ NL , Θ EV , and Θ ME . Because of the large number of parameters to estimate for each model, we propose three coordinate descent algorithms [36]–[38] that sequentially update the different parameters associated with each model, denoted as, CD A- NL, CDA-EV and CD A-ME. For each algorithm and in each step, the posterior distribution is maximized w .r .t. one parameter , the other being fixed. Thus, the algorithm iterativ ely updates each parameter by maximizing its conditional distribution as follows (see also the Appendix): July 20, 2016 DRAFT 14 • Conditional of A : truncated Gaussian distribution (whose one maximum is obtained with SUNSAL-FCLS 3 [27]) • Parameters of the residual term – Conditional of Γ : positiv e truncated Gaussian distrib ution (whose one maximum is obtained with SUNSAL-CLS [27]) – Conditional of K : Gaussian distribution (analytical expression of the mean) – Conditional of D : Gaussian distribution (analytical expression of the mean) • Conditional of  : in verse gamma distribution (analytical expression of the maximum) • Conditional of w : gamma distribution (analytical expression of the maximum) • Conditional of σ 2 : in verse gamma distribution (analytical expression of the maximum) • Conditional of c ME : Gaussian distribution (analytical expression of the mean) • Conditional of c NL : uncommon distribution (see the Appendix). Regarding the sequence generated by the coordinate descent algorithm, the proposition 2.7.1 in [36] asserts that its limit points are stationary points of (21) provided that the maximum of that function w .r .t. Θ along each coordinate is unique. This is easily checked for all the parameters e xcept for c NL . Indeed, the cost function writes as a 4 -order polynomial w .r .t. c NL (leading to 3 possible maxima) and we ha ve chosen the one that maximizes it in the interv al [0 . 2 , 3] . Algo. 1 gathers the main steps of the proposed three algorithms. If a line is in volv ed in a specific algorithm, it will begin by its name. For example, line 11 is only ex ecuted when considering the algorithms CD A-NL and CD A-ME. The cost function is not con vex, thus, the solution obtained might depend on the initial values that should be chosen carefully . Therefore, the ab undances A are initialized with SUNSAL-FCLS [27], the residual terms are initialized by 0 , the noise variance is initialized by HYSIME [41], the illumination coef ficient c is initialized by considering the sum of the abundances that are estimated using only the positi vity constraint with SUNSAL-CLS [27]. W ith these initializations, the proposed algorithm reached minima of “good quality” in the considered simulations (see Sections V and VI). Therefore, the CD A algorithms constitute a good balance between computational ef ficiency (obtained by solving the simple problems associated with each descent step) and the quality of the solution (experimentally observed when considering a good initialization). Finally , we have empirically observed that the algorithm is less sensiti ve to the initialization when beginning the decent with an abundance update (as in Algo. 1). 3 SUNSAL-FCLS satisfies the PSTO constraints while SUNSAL-CLS only ensure the positivity constraint. July 20, 2016 DRAFT 15 Algorithm 1 Coordinate descent algorithm 1: Estimate M using an EEA (VCA [24], N-FINDR [25], ...) 2: Initialization 3: Initialize parameters A (0) (SUNSAL-FCLS), Γ (0) = 0 D,N , K (0) = 0 L,R,N , D (0) = 0 L,N , c (0) (SUNSAL-CLS), Σ (0) (HYSIME),  (0) and t = 1 4: con v = 0 , 5: Parameter update 6: while con v = 0 do 7: Update A ( t ) with SUNSAL-FCLS 8: CD A-NL: Update Γ ( t ) with SUNSAL-CLS 9: CD A-EV : Update K ( t ) by a standard least squares (analytical expression) 10: CD A-ME: Update D ( t ) by a standard least squares (analytical expression) 11: CD A-NL and CD A-ME: Update (  , w ) ( t ) (analytical expression) 12: Update Σ ( t ) (analytical expression) 13: CD A-ME: Update c ( t ) by a standard least squares (analytical expression) 14: CD A-NL: Update c ( t ) by the resolution of a 3 rd order polynomial 15: Set con v = 1 if the con ver gence criteria are satisfied 16: t = t + 1 17: end while 1) Stopping criteria: Algo. 1 is an iterativ e algorithm that requires the definition of some stopping criteria. In this paper , we hav e considered four criteria and the algorithm is stopped if one of them is satisfied. The first criterion compares the new value of the cost function to the previous one and stops the algorithm if the relativ e error between these two v alues is smaller than a giv en threshold, i.e., |F  Θ t +1  − F  Θ t  | ≤ ξ 1 F  Θ t  , (23) where | . | denotes the absolute v alue and the cost function F ( Θ ) = − log [ f ( Θ | Y )] is the negati ve log-posterior . The second criterion ev aluates the new abundance values and stops the algorithm if the follo wing condition is satisfied [42]    A ( t +1) − A ( t )    F ≤ ξ 2     A ( t )    F + ξ 2  . (24) where k . k F denotes the Frobenius norm. The third criterion apply the same condition to the residual term ( Γ or K or D ) while considering a different threshold ξ 3 . The last criterion July 20, 2016 DRAFT 16 is based on a maximum number of iterations T max . These thresholds ha ve been fixed as follo ws ( ξ 1 , ξ 2 , ξ 3 , T max ) = (10 − 5 , 10 − 6 , 10 − 6 , 500) for CD A-NL, ( ξ 1 , ξ 2 , ξ 3 , T max ) = (5 × 10 − 6 , 10 − 4 , 10 − 6 , 500) for CD A-EV and ( ξ 1 , ξ 2 , ξ 3 , T max ) = (10 − 5 , 10 − 6 , 10 − 11 , 500) for CD A- ME. The next sections study the behavior of the proposed algorithms when considering synthetic and real images. 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 D A TA This section ev aluates the performance of the proposed algorithms with synthetic data. The first part introduces the criteria used for the ev aluation of the unmixing quality . In the second part, we ev aluate and compare the performance of the proposed algorithms with the state- of-the-art algorithms when considering different unmixing scenarios. In the third part, we analyze the behavior of the proposed algorithms when v arying the number of endmembers, spectral bands, pixels and the signal-to-noise ratio (SNR). A. Evaluation criteria For synthetic images, the abundances are known and the unmixing quality can be ev al- uated by using the root mean square error: RMSE ( A ) = q 1 N R P N n =1 k a n − ˆ a n k 2 . The unmixing performance can also be ev aluated by considering the reconstruction error: RE = q 1 N L P N n =1 k ˆ y n − y n k 2 and the spectral angle mapper SAM = 1 N P N n =1 arccos  ˆ y T n y n k y n k k ˆ y n k  criteria, where arccos( · ) is the in verse cosine operator and y n , ˆ y n denotes the # n th measured and estimated pixel spectra. B. P erformance of the pr oposed algorithms This section ev aluates the performance of the proposed unmixing algorithms when consid- ering different mixture models. Four synthetic images of size 100 × 100 pixels and L = 207 spectral bands hav e been generated using R = 3 endmembers corresponding to spectral signatures av ailable in the ENVI software library [43]. All images hav e been corrupted by i.i.d. Gaussian noise (with SNR= 25 dB) for a fair comparison with SU algorithms using this assumption. The images have been generated using different mixture models as follows • Linear model: the image I 1 has been generated according to the LMM model with abundances uniformly distrib uted in the simplex S . The illumination coef ficient increases linearly from the left of the image to its right in the interval [0 . 9 , 1 . 15] . July 20, 2016 DRAFT 17 • Nonlinearity: the image I 2 has been generated according to 4 linear/nonlinear models. For this, an image partition into 4 classes has been generated by considering a Potts- Marko v random field (with granularity parameter β = 0 . 8 ). The four spatial classes are associated with the LMM, NL model (3) (with  2 = 0 . 1 ), GBM (with random nonlinear coef ficients in [0 . 8 , 1] ) and PPNMM (with b = 0 . 5 ), respectiv ely . The illumination coef- ficient increases linearly from the left of the image to its right in the interval [0 . 9 , 1 . 15] . Finally , to mak e the SU ev en more challenging, we have considered a highly mix ed scenario by generating the abundance using a Dirichlet distribution whose parameters are selected randomly in the interval [1 , 20] . • Endmember variability: the image I 3 has been partitioned into 4 classes (the same classification map as for I 2 ). In each class, we hav e generated a set of endmembers by considering the model (5) and k r,i,j ∼ N ( 0 L , α 2 H ) , with α 2 = 0 . 005 , ∀ r , i, j . Therefore, the image will be composed by four sets of endmembers, each one associated with a spatial region. T o make the SU ev en more challenging, we ha ve considered a highly mixed scenario by generating the abundance using a Dirichlet distribution whose parameters are selected randomly in the interval [1 , 20] . • Mismodelling effects: the image I 4 has been partitioned into 2 classes (by merging class 1 with 2 , and 3 with 4 of the classification map of I 2 ). Pixels of the first class have been generated according to the ME model (8) with  2 = 0 . 002 . The pixels of class 2 hav e been generated with 4 endmembers to simulate the effect of a bad estimation of the number of endmembers that will be fixed to R = 3 for all unmixing algorithms. The illumination coefficient increases linearly from the left of the image to its right in the interv al [0 . 9 , 1 . 15] . Finally , to make the SU ev en more challenging, we have considered a highly mixed scenario by generating the abundance using a Dirichlet distribution whose parameters are selected randomly in the interval [1 , 20] . These images are processed using different unmixing strategies that are compared to the proposed algorithms. For all algorithms, we ha ve assumed the endmembers to be known and we have considered the ENVI spectra used to create the images. The studied unmixing algorithms are • Linear unmixing: the abundances are estimated using the FCLS algorithm [26] and the July 20, 2016 DRAFT 18 SUNSAL-CLS 4 algorithm [27]. • Nonlinear unmixing: the ab undances are estimated using the MCMC-RCA algorithm [33] and the SKHYPE algorithm [29] • Endmember v ariability: we hav e considered the automated endmember bundles (AEB) algorithm proposed in [16]. In addition to the ENVI spectra, the endmember library includes endmembers extracted using VCA in 10% (resp. 20% ) image subset when processing synthetic (resp. real) images. For each pixel, the R endmembers that provide the smallest RE are selected. The next sections analyse the obtained results w .r .t. each synthetic image. 1) Linear model: T able II sho ws the obtained results when considering the LMM based image I 1 . Considering the abundance RMSE, and as expected, SUNSAL-CLS provides better results than FCLS because of the variation of the illumination parameter c for this image. Indeed, varying c can be seen as a relaxation of the sum-to-one constraint. This variability also af fects RCA-MCMC that presents reduced performance when compared to the other nonlinear based algorithms (CD A-NL and SKhype). The proposed algorithms show v ery good results especially CD A-NL and CD A-ME which validate their use for linear mixture model. Except FCLS, all algorithms provide good reconstruction spectra which is highlighted by the RE and SAM criteria. This table highlights that the LMM based algorithms are the fastest ones (this result is valid for all images) followed by the CD A-NL algorithm which is ≈ 100 times faster than the MCMC based algorithm for this image. Overall, CD A algorithms provide a reduced computational times which sho ws their efficienc y for linear unmixing. Note that CDA-EV requires more computational times mainly because it estimates large matrices K . Finally , the obtained results confirm the good beha vior of the proposed CD A algorithms when processing LMM based images. 2) Nonlinearity: T able III shows the obtained results when considering image I 2 . This table also provides the obtained RMSE associated with each spatial class, i.e., each nonlinear mixture model (the results of class 1 are similar to those obtained for I 1 ). As expected, the linear based algorithms provide poor results which highlights the need for more elaborated algorithms to deal with the NL effect present in the image. Ho wev er, it can be seen that removing the sum-to-one constraint in SUNSAL-CLS improves considerably the performance 4 The SUNSAL algorithm is run with the parameters suggested in [27], i.e., λ = 0 , maximum of iterations = 200 and tolerance = 10 − 4 . July 20, 2016 DRAFT 19 w .r .t. FCLS especially for GBM and PPNMM. This suggest that part of the nonlinearity can be interpreted as a variation of illumination, which is in agreement with [9], [33] that suggest removing this constraint when considering NL models. The best performance are obtained by considering nonlinear-based algorithms. More precisely , CD A-NL is robust to the dif ferent nonlinearities af fecting the data and provides the best RMSEs for I 2 . Note also that CD A-ME sho ws an intermediate performance between the NL models and the linear/EV based models. Regarding the computational times, CD A-NL requires less times than SKhype to achie ve the unmixing, and again outperforms RCA-MCMC (by a factor of 20 ) that shows expensi ve computational cost. These results confirm the good behavior of CD A-NL when considering NL images. They also show that CD A-ME is an intermediate solution between NL and LMM based algorithms since it requires less computational times than NL algorithms while it shows better results than LMM algorithms. 3) Endmember variability: The results associated with I 3 are summarized in T able II. Again, LMM based algorithms suffer from the presence of EV effect. As expected, the best performance in terms of ab undance RMSE and RE are obtained with the CD A-EV algorithm that is especially designed to deal with the EV effect. Similar good performance are obtained by CD A-ME algorithm which confirms its robustness with respect to the observed mixture model (LMM, NL or EV). Note that AEB provides bad RMSE results even though it shows a reduced RE. This result shows that a reduced RE does not ensure a good estimation of the abundances. Note finally that apart from the LMM based algorithms, the CD A algorithms sho w competitiv e computational times when compared to the remaining algorithms. These results confirm the good performance of the CD A-EV and the robustness of CD A-ME when processing images with EV . 4) Mismodelling effects: This section analyses the obtained results when processing image I 4 that is corrupted by mismodelling ef fects (presence of outliers and a fourth endmember). For this scenario, CD A-ME provides the best results in term of abundance RMSEs, RE and SAM as sho wn in T able II. Note that CD A-EV also sho ws interesting results with a reduced computational cost. Nonlinear algorithms suffer from this corruption and show high abundance RMSEs. These results confirm the good performance of CD A-ME in terms of unmixing quality and computational times. T o summarize, the proposed CD A algorithms show reduced computational cost for all images. The best results associated with an image are obtained when considering the corre- sponding CD A algorithm (for example, CD A-EV is best for EV -based images). CD A-ME is July 20, 2016 DRAFT 20 robust to the different mixture models that can af fect hyperspectral images. C. Robustness of the pr oposed algorithms This section ev aluates the robustness of the proposed CD A algorithms when varying the parameters ( R, L, N , SNR ) . Synthetic images hav e been generated according to the proposed RCA models (NL, EV and ME) with the follo wing configurations: • The illumination coefficient increases linearly from the left of the image to its right in the interval [0.9; 1.15]. • The abundances are uniformly distributed in the simplex. • The NL coefficient is fixed to  2 = 0 . 1 in (13), the EV coef ficient is fixed to α 2 = 0 . 005 , ∀ r , i, j in (15), and the ME coefficient is fixed to  2 = 0 . 002 in (16). Each CD A algorithm has been run on its corresponding image, i.e., CD A-NL for the NL image, CD A-EV for the EV image and CD A-ME for the ME image. T able IV reports the obtained RMSE results when v arying ( R , L, N , SNR ) . When not varying, the parameters are fixed to L = 207 bands, N = 10 4 pixels, R = 3 endmembers and SNR = 25 dB. T able IV sho ws that the performance decreases when increasing the number of endmembers R . This result is similar to [12], [44], since increasing R leads to an increase in the number of the estimated parameters which reduces the RMSE performance. In addition, T able IV highlights the decrease in performance when reducing the observations either by reducing the number of bands # L or the number of pixels # N . Note, ho we ver , that the results are still acceptable for L ≥ 52 bands and N ≥ 625 pixels. As expected, the obtained results show that the higher the SNR, the better the RMSE performance is. These results are in good agreement with the literature and confirm the good behavior of the CD A algorithms for dif ferent parameter configurations. Note finally that more results are provided in the technical report [45] and are not shown here for brevity . V I . S I M U L AT I O N R E S U L T S O N R E A L D A TA This section illustrates the performance of the proposed algorithms when applied to tw o real hyperspectral images. The first hyperspectral image has receiv ed much attention in the remote sensing community [5], [40]. This image was acquired ov er Mof fett Field, CA, in 1997 by the A VIRIS. The considered dataset contains 100 × 100 pixels, L = 152 spectral bands (after removing w ater absorption bands) acquired in the interv al 0 . 4 − 2 . 5 µ m, has a spatial resolution of 100 m and is mainly composed of three components: water , soil, and vegetation (see Fig. July 20, 2016 DRAFT 21 3 (left)). For this image, three endmembers were extracted using VCA [24]. The second image, denoted as Madonna, was acquired in 2010 by the Hyspex hyperspectral scanner ov er V illelongue, France (00 03’W and 4257’N). The dataset contains L = 160 spectral bands recorded from the visible to near infrared ( 400 − 1000 nm) with a spatial resolution of 0 . 5 m [46]. It has already been studied in [12], [44] and is mainly composed of forested areas (see Fig. 3 (right)). This image contains 100 × 100 pixels and is composed of R = 4 components: tree, grass, soil and shado w (see Fig. 3 (right)). The Bayesian UsLMM algorithm [47] was used to estimate R = 4 endmembers. A. Unmixing performance The ab undances of each image ha ve been estimated by the considered unmixing algorithms. Fig. 4 sho ws the obtained results for the Moffett image. Note that a white (black) pixel indicates a lar ge (small) proportion of the corresponding materials. Except SUNSAL-CLS, the considered algorithms show similar abundance maps. The behavior of SUNSAL-CLS is due to the presence of a high illumination variation for this image. In addition, some nonlinearity can be interpreted as an illumination variability (as already shown when analysing synthetic data) leading to high values for the parameter c . This beha vior will be further analysed in the next sections. Note that the algorithms estimated similar abundance maps for the Madonna image and for conciseness, the abundances maps are not displayed. The unmixing performance can also be compared by considering the reconstruction errors as shown in T able V. Considering the Mof fett image, the best RE and SAM results are obtained by CD A-ME follo wed by CD A-EV and CD A-NL. These results suggest the presence of both EV and NL in the Mof fett image which requires the use of a robust unmixing algorithm such as CD A-ME (a detailed study of NL and EV will be provided in the next sections). The best performance on the Madonna image are obtained by the CD A-EV algorithm, while we obtain good results when considering the algorithms including EV such as AEB, CD A-ME and SUNSAL-CLS. These results suggest the presence of a higher variability ef fect than nonlinearity in the Madonna image. Note finally that the proposed CD A algorithms sho w a reduced computational cost while pro viding man y details reg arding the physical ef fects corrupting the images as sho wn in the next sections. July 20, 2016 DRAFT 22 B. Residual components This section analyses the obtained residual energies for the considered real images. Fig. 5 sho ws the energies of the difference between the reconstructed signal and the linear model (i.e., || ˆ y i,j − M ˆ a i,j || ) for the Moffett image. For RCA-MCMC, the reconstructed spectrum is not a vailable and we only sho w the estimated nonlinearity le vels. Fig. 5 shows a good agreement between the results of the proposed CD A algorithms where mismodelling effects are mainly located in the coastal zone, and in the right area (composed of ≈ 90% soil and ≈ 10% vegetation). Similarly to RCA-MCMC, the three CD A algorithms show some mismodelling effects in the left zone of the water area. Ho we ver , the energy location shows some dif ferences between RCA-MCMC and CD A-algorithms mainly because of the presence of both NL and EV in the Moffett scene as shown in the previous section (see T able V). Indeed, all CDA algorithms account for EV illumination effect while RCA-MCMC only account for NL effect. The same energies were displayed when considering the Madonna image in Fig. 6. This figure shows a good agreement between the results of RCA-MCMC and CD A algorithms where the high deviation from the linear model is mainly located in the tree and shado w areas (as in [44]). This makes sense since multiscattering effects (i.e., nonlinear interactions) are generally located in these areas. Note finally that for both CD A-NL and CD A-ME, the displayed maps in Figs. 5 and 6 include the effect of the illumination variation (introduced by c i,j ) and that of the residual term φ . These two terms are studied separately in the next section. C. Nonlinearity and illumination As previously noticed, the Moffett image includes both NL and EV ef fects. Fig. 7 separates the residual energies into two components: the first one is related to illumination parameter c and the second one to the term φ . Considering CDA-NL, it can be seen that the nonlinearity captured by φ N L is mainly located in the coastal region and areas comprising multiple physical components (as in [5]). This algorithm also captures a high variation in illumination in the left of the w ater area which has also been captured by RCA-MCMC (see Fig. 5), ho we ver , it has been interpreted as a nonlinear effect. These results sho w that there is a close relation between the effects of illumination v ariation and nonlinearity . This remark can be further justified by considering the results of CD A-ME which captures both EV and NL effects. Fig. 7 (left-bottom) shows that most of the NL effect can be interpreted as July 20, 2016 DRAFT 23 a variation of illumination while the remaining effects can be approximated by a residual smooth term φ M E . Therefore, in presence of nonlinearity , a great unmixing improvement can be obtained by allo wing an illumination v ariation, which also justifies the great improvement of SUNSAL-CLS with respect to FCLS for the Mof fett image (see T able V). Considering Madonna image, Fig. 8 (top) sho ws a reduced NL ef fect and most of the residuals can be interpreted as a variation in illumination. Similarly to CD A-NL, CD A-ME captures an illumination v ariation that is mainly located in the region of trees. In addition, Fig. 8 (right-bottom) sho ws that CD A-ME also captures the shape EV as for CD A-EV . Indeed, CD A-ME is an intermediate algorithm that can deal with illumination variation, NL and EV ef fects. T o conclude, this section shows that most of the NL effect can be captured by varying the illumination coefficient. Moreov er , CD A algorithms capture similar spatial residual effects but giv e a dif ferent interpretation depending on the considered mixture model. Most of these residuals appear in region of intersection between the physical elements and in presence of ve getation such as trees. Finally , CD A-ME is very flexible and can capture different kind of mismodelling effects. D. Endmember variability CD A-EV estimates a set of endmember spectra for each pixel to account for EV . Figs. 9 and 11 compare the estimated spectra with those extracted with AEB, VCA or UsLMM. These figures also show the interval of endmembers obtained with CD A-ME when interpreting the mismodelling term d as to be due to EV . Overall, the estimated spectra are in good agreement with the state-of-the-art algorithms especially for the Moffett image. These figures show that CD A-ME captures effects that are not only due to EV since neg ati ve spectra may appear . Figs. 10 and 12 show the spatial repartition of the captured EV by CD A-EV . For both images, the highest EV is obtained for ve getation (tree and grass) and for regions with multiple components. For the Moffett image, and as an example, the water variability is located in the left area of the water zone while the soil variability is concentrated in the coastal area. For the Madonna image, the tree variability appears in the forest region while the shadow v ariability is located in the shadowed regions. These results sho w the ability of CD A-EV to provide both the spectral behavior and the spatial location of EV . W e conclude this section by providing some comments regarding the proposed algorithms and their use in practice. It has been shown in this section that CD A-ME can capture both NL and EV in presence of vegetation, terrain relief or multiple physical components. In July 20, 2016 DRAFT 24 addition to this flexibility , CD A-ME is generally fast, which encourages its use as a first step to analyze a scene. The obtained results can then be analyzed to highlight the presence of EV or NL. As a second step, CD A-NL or CD A-EV can be used to focus on the phenomenon of interest. CD A-EV is designed to highlight the sensiti ve spectral bands as well as the shape of the spectral variation and the re gions of spatial variation of each materiel. CD A-NL is designed to highlight the region of interaction of the physical components as well as the acti ve bilinear components. Finally , we provide T able VI that summarizes some characteristics and dif ferences between the proposed models. V I I . C O N C L U S I O N S This paper introduced three hyperspectral mixture models and their associated Bayesian algorithms for supervised hyperspectral unmixing. The three models were introduced under a general formulation that can be adapted to account for nonlinearity , endmember variability or mismodelling effects. A hierarchical Bayesian model was proposed for each model to introduce the known constraints on their parameters. Those parameters were estimated using a coordinate descent algorithm that showed a reduced computational cost when compared to state-of-the-art algorithms. The proposed algorithms sho wed good performance when processing synthetic data generated with the linear model or other more sophisticated models. Results on real data confirmed the good performance of the proposed algorithms and sho wed their ability to e xtract different features in the observed scenes. These results showed that most of the nonlinearity can be interpreted by a variation in illumination, that endmember v ariability is mainly located in ve getation areas, and that residual effects (such as nonlinearity or endmember variability) appears in presence of multiple physical components. Future work includes the estimation of the hyperparameters associated with the proposed algorithms, since they control the degree of smoothness of the obtained results (see for instance Figs. 10 and 12). It also includes the estimation of the endmember number using a model selection criteria such as the Akaike information criterion (AIC) [48], the minimum description length (MDL) [49], or other methods such as [41], [50]. Adapting the proposed algorithms to include a dimension reduction step is an interesting issue which should further reduce their computational cost. Finally , detecting the presence of endmember variability and nonlinearity using hypothesis tests is also an interesting issue which would deserve to be in vestigated. July 20, 2016 DRAFT 25 A P P E N D I X C O N D I T I O N A L D I S T R I B U T I O N S This appendix provides the conditional distributions of the estimated parameters where S i,j and φ i,j depend on the considered observation model (NL, EV or ME). The conditional distributions are giv en by                                                        a i,j | Θ \ a i,j ∼ N S  µ a i,j , Λ a i,j  γ i,j | Θ \ γ i,j ∼ N ( R +) D  µ γ i,j , Λ γ i,j  k r,i,j | Θ \ k r,i,j ∼ N  µ k r,i,j , Λ k r,i,j  d i,j | Θ \ γ i,j ∼ N  µ d i,j , Λ d i,j  c M E i,j | Θ \ c i,j ∼ N  µ c i,j , Λ c i,j  f  c N L i,j | Θ \ c i,j  ∝ exp  P 4 m =0 x m c m i,j  σ 2 ` | Θ \ σ 2 ` ∼ I G  ϕ ` + N 2 , ψ ` + P i,j [ y i,j ( ` ) − S i,j ( ` ) a i,j − φ i,j ( ` ) ] 2 2  NL:  2 i,j | Θ \  2 i,j ∼ I G  4 ζ + D 2 , γ > i,j γ i,j 2 + 4 ζ ρ 1 ,i,j ( w )  ME:  2 i,j | Θ \  2 i,j ∼ I G  4 ζ + L 2 , d > i,j H − 1 d i,j 2 + 4 ζ ρ 1 ,i,j ( w )  (25) where                                                              Λ a i,j =  S > i,j Σ − 1 S i,j  − 1 µ a i,j = Λ a i,j S > i,j Σ − 1  y i,j − φ i,j  Λ γ i,j =  c 4 i,j Q > Σ − 1 Q + 1  2 i,j I D  − 1 µ γ i,j = c 2 i,j Λ γ i,j Q > Σ − 1  y i,j − c i,j p i,j  Λ k r,i,j =  1 β 2 r I L + 1 α 2 r H − 1 + a 2 r,i,j Σ − 1  − 1 µ k r,i,j = Λ k r,i,j  a 2 r,i,j Σ − 1 ω r,i,j + 1 β 2 r µ r,i,j  Λ d i,j =  Σ − 1 + 1  2 i,j H − 1  − 1 µ d i,j = Λ d i,j Σ − 1  y i,j − c i,j p i,j  Λ c i,j =  1 η 2 + p > i,j p i,j  − 1 µ c i,j = Λ c i,j h p > i,j Σ − 1  y i,j − φ i,j  + 1 η 2 i (26) with Q =  m 1  m 1 , · · · , m R  m R , √ 2 m 1  m 2 , · · · , √ 2 m R − 1  m R  , p i,j = M a i,j , ω r,i,j = 1 a r,i,j  y i,j − p i,j − P r 0 6 = r a r 0 ,i,j k r 0 ,i,j  , z i,j = Q i,j γ i,j , x 0 = y > i,j Σ − 1 y i,j , x 1 = y > i,j Σ − 1 p i,j , x 2 = 0 . 5  p > i,j Σ − 1 p i,j − 2 z > i,j Σ − 1 y i,j  , x 3 = z > i,j Σ − 1 p i,j and x 4 = 0 . 5 z > i,j Σ − 1 z i,j . July 20, 2016 DRAFT 26 Note that f  w 2 i,j | Θ \ w 2 i,j  is giv en by (19). The independence between the parameters leads to f  A | Θ \ A  = Q i,j f  a i,j | Θ \ a i,j  , f  Γ | Θ \ Γ  = Q i,j f  γ i,j | Θ \ γ i,j  , f  D | Θ \ D  = Q i,j f  d i,j | Θ \ d i,j  , f  c | Θ \ c  = Q i,j f  c i,j | Θ \ c i,j  , f  Σ | Θ \ Σ  = Q ` f  σ 2 ` | Θ \ σ 2 `  , f   | Θ \   = Q i,j f   2 i,j | Θ \  2 i,j  , and f  w | Θ \ w  = Q i,j f  w 2 i,j | Θ \ w 2 i,j  . Note finally that the matrix K is updated iterativ ely using a checkerboard scheme. R E F E R E N C E S [1] B. W . Hapke, “Bidirectional reflectance spectroscopy . I. Theory , ” J. Geophys. Res. , vol. 86, pp. 3039– ˝ U3054, 1981. [2] R. Heylen, M. Parente, and P . Gader , “ A revie w of nonlinear hyperspectral unmixing methods, ” IEEE J. Sel. T opics Appl. Earth Observat. Remote Sens. , vol. 7, no. 6, pp. 1844–1868, June 2014. [3] 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. [4] Y . Altmann, A. Halimi, N. Dobigeon, and J.-Y . T ourneret, “Supervised nonlinear spectral unmixing using a postnonlinear mixing model for hyperspectral imagery , ” IEEE T rans. Image Process. , vol. 21, no. 6, pp. 3017–3025, June 2012. [5] A. Halimi, Y . Altmann, N. Dobigeon, and J.-Y . T ourneret, “Nonlinear unmixing of hyperspectral images using a generalized bilinear model, ” IEEE T rans. Geosci. Remote Sens. , vol. 49, no. 11, pp. 4153–4162, 2011. [6] ——, “Unmixing hyperspectral images using the generalized bilinear model, ” in Proc. IEEE Int. Conf. Geosci. Remote Sens. (IGARSS) , July 2011, pp. 1886–1889. [7] J. M. Bioucas-Dias and J. M. P . Nascimento, “Nonlinear mixture model for hyperspectral unmixing, ” in Pr oc. SPIE Image and Signal Processing for Remote Sensing XV , L. Bruzzone, C. Notarnicola, and F . Posa, Eds., vol. 7477, no. 1. SPIE, 2009, p. 74770I. [8] W . Fan, B. Hu, J. Miller, and M. Li, “Comparativ e study between a new nonlinear model and common linear model for analysing laboratory simulated-forest hyperspectral data, ” International Journal of Remote Sensing , vol. 30, no. 11, pp. 2951–2962, June 2009. [9] I. Meganem, P . Deliot, X. Briottet, Y . Deville, and S. Hosseini, “Linear-quadratic mixing model for reflectances in urban environments, ” IEEE T rans. Geosci. Remote Sens. , vol. 52, no. 1, pp. 544–558, Jan 2014. [10] B. Somers, G. P . Asner , L. Tits, and P . Coppin, “Endmember variability in spectral mixture analysis: A revie w , ” Remote Sensing of Envir onment , vol. 115, no. 7, pp. 1603 – 1616, 2011. [11] A. Zare and K. Ho, “Endmember variability in hyperspectral analysis: Addressing spectral variability during spectral unmixing, ” IEEE Signal Pr ocess. Mag. , vol. 31, no. 1, pp. 95–104, Jan 2014. [12] A. Halimi, N. Dobigeon, and J.-Y . T ourneret, “Unsupervised unmixing of hyperspectral images accounting for endmember variability , ” IEEE T rans. Image Pr ocess. , vol. 24, no. 12, pp. 4904–4917, 2015. [13] 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 Sensing of En vir onment , vol. 65, no. 3, pp. 267–279, Sept. 1998. [14] C. Bateson, G. Asner, and C. W essman, “Endmember bundles: a new approach to incorporating endmember variability into spectral mixture analysis, ” IEEE T rans. Geosci. Remote Sens. , vol. 38, no. 2, pp. 1083–1094, Mar 2000. [15] 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. July 20, 2016 DRAFT 27 [16] 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 Observat. Remote Sens. , vol. 5, no. 2, pp. 396–408, April 2012. [17] J. M. P . Nascimento and J. M. Bioucas-Dias, “Does independent component analysis play a role in unmixing hyperspectral data?” IEEE T rans. Geosci. Remote Sens. , vol. 43, no. 1, pp. 175–187, Jan. 2005. [18] M. A. V eganzones, L. Drumetz, G. T ochon, M. Dalla Mura, A. Plaza, J. M. Bioucas-Dias, and J. C., “A New Extended Linear Mixing Model to Address Spectral V ariability, ” in IEEE W orkshop on Hyperspectral Image and Signal Processing: Evolution in Remote Sensing (WHISPERS 2014) , Lausanne, Switzerland, June 2014. [19] G. Shaw and H. Burke, “Spectral imaging for remote sensing, ” Lincoln Lab . J. , vol. 14, no. 1, pp. 3–28, 2003. [20] 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, June 2010. [21] A. Zare, P . Gader , and G. Casella, “Sampling piece wise conv ex unmixing and endmember extraction, ” IEEE T rans. Geosci. Remote Sens. , vol. 51, no. 3, pp. 1655–1665, March 2013. [22] D. Stein, “ Application of the normal compositional model to the analysis of hyperspectral imagery , ” in Advances in T echniques for Analysis of Remotely Sensed Data, 2003 IEEE W orkshop on , Oct 2003, pp. 44–51. [23] X. Du, A. Zare, P . Gader , and D. Dranishnikov , “Spatial and spectral unmixing using the beta compositional model, ” IEEE J. Sel. T opics Appl. Earth Observat. Remote Sens. , vol. 7, no. 6, pp. 1994–2003, June 2014. [24] J. M. P . 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, April 2005. [25] M. W inter , “Fast autonomous spectral end-member determination in hyperspectral data, ” in Pr oc. 13th Int. Conf. on Applied Geologic Remote Sensing , vol. 2, V ancouver , April 1999, pp. 337–344. [26] 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. [27] J. Bioucas-Dias and M. Figueiredo, “ Alternating direction algorithms for constrained sparse regression: Application to hyperspectral unmixing, ” in Pr oc. IEEE GRSS W orkshop on Hyperspectral Image and SIgnal Pr ocessing: Evolution in Remote Sensing (WHISPERS) , June 2010, pp. 1–4. [28] Y . Altmann, M. Pereyra, and S. McLaughlin, “Bayesian nonlinear hyperspectral unmixing with spatial residual component analysis, ” IEEE T rans. Image Pr ocess. , vol. 1, no. 3, pp. 174–185, Sept. 2015. [29] J. Chen, C. Richard, and P . Honeine, “Nonlinear unmixing of hyperspectral data based on a linear-mixture/nonlinear- fluctuation model, ” IEEE T rans. Signal Pr ocess. , vol. 61, no. 2, pp. 480–492, Jan 2013. [30] Y . Altmann, S. McLaughlin, and A. Hero, “Robust linear spectral unmixing using anomaly detection, ” vol. 1, no. 2, pp. 74–85, June 2015. [31] C. Févotte and N. Dobigeon, “Nonlinear hyperspectral unmixing with robust nonnegativ e matrix factorization, ” IEEE T rans. Image Pr ocess. , vol. 24, no. 12, pp. 4810–4819, June 2015. [32] A. A. Kalaitzis and N. D. Lawrence, “Residual components analysis, ” in Proc. ICML , 2012, pp. 1–3. [33] 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 Pr ocess. , vol. 23, no. 5, pp. 2148–2158, May 2014. [34] H. Rue and L. Held, Gaussian Markov random fields: theory and applications . CRC Press, 2005. [35] O. Dikmen and A. Cemgil, “Gamma markov random fields for audio source modeling, ” IEEE T rans. Audio, Speech, Language Pr ocess. , vol. 18, no. 3, pp. 589–601, March 2010. [36] D. P . Bertsekas, Nonlinear progr amming . Belmont, Massachusetts: Athena Scientific, 1995. July 20, 2016 DRAFT 28 [37] J. Sigurdsson, M. Ulfarsson, and J. Sveinsson, “Hyperspectral unmixing with l q regularization, ” IEEE T rans. Geosci. Remote Sens. , vol. 52, no. 11, pp. 6793–6806, Nov . 2014. [38] 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. , vol. 54, no. 4, pp. 2207–2219, Mar . 2016. [39] 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 Pr ocess. , vol. 64, no. 2, pp. 525–538, Jan 2016. [40] 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. [41] J. M. Bioucas-Dias and J. M. P . Nascimento, “Hyperspectral subspace identification, ” IEEE T rans. Geosci. Remote Sens. , vol. 46, no. 8, pp. 2435–2445, Aug. 2008. [42] K. Madsen, H. B. Nielsen, and O. T ingleff, “Methods for non-linear least squares problems (2nd ed.), ” Richard Petersens Plads, Building 321, DK-2800 Kgs. L yngby , p. 60, 2004. [43] RSI (Research Systems Inc.), ENVI User’ s guide V ersion 4.0 , Boulder , CO 80301 USA, Sept. 2003. [44] Y . Altmann, N. Dobigeon, S. McLaughlin, 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, June 2014. [45] A. Halimi, P . Honeine, and J.Bioucas-Dias, “T echnical report associated with the paper “Hyperspectral unmixing in presence of endmember variability , nonlinearity or mismodelling effects”, ” Université de technologie de T royes, France, T ech. Rep., April 2016. [Online]. A vailable: https://sites.google.com/site/abderrahimhalimi/publications [46] D. Sheeren, M. Fauvel, S. Ladet, A. Jacquin, G. Bertoni, and A. Gibon, “Mapping ash tree colonization in an agricultural mountain landscape: Inv estigating the potential of hyperspectral imagery , ” in Pr oc. IEEE Int. Conf. Geosci. Remote Sens. (IGARSS) , July 2011, pp. 3672–3675. [47] 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 Process. , vol. 57, no. 11, pp. 4355–4368, Nov . 2009. [48] H. Akaike, “ A new look at the statistical model identification, ” IEEE T rans. Autom. Contr . , vol. 19, pp. 716–723, 1974. [49] J. Rissanen, “Modeling by shortest data description, ” Automatica , vol. 14, pp. 465–471, 1978. [50] A. Halimi, P . Honeine, M. Kharouf, C. Richard, and J. Y . T ourneret, “Estimating the intrinsic dimension of hyperspectral images using a noise-whitened eigengap approach, ” IEEE T rans. Geosci. Remote Sens. , vol. 54, no. 7, pp. 3811–3821, 2016. July 20, 2016 DRAFT FIGURES 29 Fig. 1. (T op) USGS spectra showing endmember variability effect. (Bottom) example of dif ference between USGS spectra (continuous lines) and their smooth approximation with Gaussian process (dashed lines). July 20, 2016 DRAFT FIGURES 30 Y c σ Γ A η (  , w ) ϕ ψ ζ (a) NL model Y σ K A α r β r ϕ ψ (b) EV model Y c σ D A η (  , w ) ϕ ψ ζ (c) ME model Fig. 2. DA Gs for the parameter and hyperparameter priors (the user fixed parameters appear in boxes). (a) the NL model, (b) the EV model and (c) the ME model. (a) Moffett image. (b) Madonna image. Fig. 3. Real hyperspectral images. (Left) Moffett image, (right) Madonna image. July 20, 2016 DRAFT FIGURES 31 Fig. 4. Estimated abundance maps with different algorithms for the Moffett image. (Left) vegetation, (middle) water , (right) soil. Fig. 5. Square root of the energies of the difference between the reconstructed signal and the linear model obtained with || ˆ y i,j − M ˆ a i,j || (for RCA-MCMC we show the estimated nonlinearity coefficients) for the Moffett image. July 20, 2016 DRAFT FIGURES 32 Fig. 6. Square root of the energies of the difference between the reconstructed signal and the linear model obtained with || ˆ y i,j − M ˆ a i,j || (for RCA-MCMC we show the estimated nonlinearity coefficients) for the Madonna image. Fig. 7. Residual maps obtained with CDA-NL and CD A-ME algorithms for the Moffett image. (Left) estimated illumination variation | 1 − c i,j | , (right) square root of the energies of the residual terms k φ i,j k . Fig. 8. Residual maps obtained with CDA-NL and CD A-ME algorithms for the Madonna image. (Left) estimated illumination variation | 1 − c i,j | , (right) square root of the energies of the residual terms k φ i,j k . July 20, 2016 DRAFT FIGURES 33 Fig. 9. The estimated R = 3 endmembers of the Moffett image with VCA (continuous red lines), AEB (continuous green lines), CDA-EV (continuous blue lines) and the interval of spectra with CD A-ME (dashed black lines). Fig. 10. Estimated spatial variability maps with CD A-EV for the Moffett image. The # r th map is obtained by computing k k i,j,r k for each pixel. Fig. 11. The estimated R = 4 endmembers of the Madonna image with UsLMM (continuous red lines), AEB (continuous green lines), CD A-EV (continuous blue lines) and the interval of spectra with CDA-ME (dashed black lines). Fig. 12. Estimated spatial variability maps with CDA-EV for the Madonna image. The # r th map is obtained by computing k k i,j,r k for each pixel. July 20, 2016 DRAFT 34 T ABLE I E S T I M A T E D PA R A M E T E R S O F T H E P R O P O S E D A L G O R I T H M S . S PA T I A L C O R R E L ATI O N I S D E N OT E D B Y ( S C ) A N D S P E C T R A L S M O OT H I N G I S D E N OT E D B Y ( S S ) Parameters CD A-NL CD A-EV CD A-ME A R × N R × N R × N Γ D × N - - K - R × L × N (ss, sc) - D - - L × N (ss)  1 × N (sc) - 1 × N (sc) σ 2 L × 1 L × 1 L × 1 c 1 × N - 1 × N T ABLE II R E S U LT S O N S Y N T H E T I C DAT A ( I M AG E I 1 , I 3 A N D I 4 ) . ( × 10 − 2 ) T ime RMSE RE SAM (s) I 1 FCLS 7 . 78 3 . 58 6 . 24 1 . 2 SUNSAL-CLS 3 . 42 2 . 27 5 . 62 0 . 1 (LMM, SKhype 1 . 41 − − 541 K = 1 ) RCA-MCMC 4 . 31 − − 6737 AEB 33 . 8 2 . 31 5 . 62 1507 CD A-NL 1 . 34 2 . 27 5 . 62 45 CD A-EV 3 . 27 2 . 24 5 . 55 142 CD A-ME 1 . 35 2 . 27 5 . 62 88 I 3 FCLS 10 . 22 2 . 91 6 . 0 1 . 3 SUNSAL-CLS 10 . 16 2 . 53 5 . 93 0 . 14 (RCA-EV , SKhype 7 . 69 − − 623 K = 4 ) RCA-MCMC 13 . 33 − − 11025 AEB 19 . 05 2 . 45 5 . 74 1465 CD A-NL 10 . 02 2 . 53 5 . 94 520 CD A-EV 3 . 60 2 . 35 5 . 53 230 CD A-ME 4 . 29 2 . 35 5 . 51 316 I 4 FCLS 12 . 59 5 . 01 7 . 19 1 . 21 SUNSAL-CLS 12 . 65 2 . 91 6 . 94 0 . 12 (RCA-ME, SKhype 11 . 68 − − 374 K = 4 ) RCA-MCMC 18 . 20 − − 6029 AEB 32 . 17 2 . 88 6 . 05 1240 CD A-NL 11 . 85 2 . 88 6 . 88 394 CD A-EV 9 . 74 2 . 36 5 . 71 279 CD A-ME 6 . 07 2 . 29 5 . 56 315 July 20, 2016 DRAFT 35 T ABLE III R E S U LT S O N S Y N T H E T I C DAT A ( I M AG E I 2 ) . RMSE (classes) ( × 10 − 2 ) RMSE RE SAM Time C 1 C 2 C 3 C 4 ( × 10 − 2 ) ( × 10 − 2 ) ( × 10 − 2 ) (s) LMM RCA-NL GBM PPNMM FCLS 10 . 24 44 . 72 15 . 48 23 . 98 24 . 76 15 . 74 10 . 64 1 . 7 SUNSAL-CLS 3 . 84 33 . 81 5 . 68 8 . 45 16 . 55 4 . 17 7 . 57 0 . 07 SKhype 1 . 67 11 . 92 2 . 21 2 . 81 5 . 87 − − 547 RCA-MCMC 5 . 87 6 . 29 5 . 44 3 . 93 5 . 66 − − 9009 AEB 53 . 33 22 . 5 49 . 2 41 . 8 45 . 72 3 . 05 6 . 46 1732 CD A-NL 1 . 62 7 . 27 2 . 16 2 . 89 3 . 86 2 . 86 6 . 16 430 CD A-EV 4 . 84 33 . 27 7 . 12 9 . 01 16 . 59 3 . 34 6 . 64 160 CD A-ME 1 . 63 13 . 55 2 . 27 2 . 76 6 . 61 2 . 89 6 . 17 37 T ABLE IV R M S E O F T H E P RO P O S E D C DA A L G O R I T H M S W H E N V A RY I N G L , N , R A N D S N R . CD A-NL CD A-EV CDA-ME V arying L L = 52 11 . 08 6 . 34 8 . 22 L = 104 9 . 04 4 . 61 5 . 89 L = 207 7 . 58 3 . 67 4 . 26 V arying N N = 169 21 . 00 9 . 93 9 . 52 N = 625 7 . 71 4 . 26 4 . 32 N = 2500 7 . 71 3 . 74 4 . 29 N = 10000 7 . 58 3 . 67 4 . 26 SNR= 15 R = 3 14 . 52 8 . 09 9 . 05 R = 6 19 . 21 9 . 67 11 . 11 SNR= 25 R = 3 7 . 57 3 . 62 4 . 26 R = 6 13 . 94 5 . 95 7 . 52 T ABLE V R E S U LT S O N R E A L I M A G E S . Moffett image Madonna image RE SAM T ime RE SAM T ime ( × 10 − 3 ) ( × 10 − 2 ) (s) ( × 10 − 3 ) ( × 10 − 2 ) (s) FCLS 13 . 6 12 . 7 2 6 . 3 2 . 7 1 SUNSAL-CLS 6 . 4 9 . 1 0 . 1 6 . 2 2 . 6 0 . 1 SKhype − − 241 − − 131 RCA-MCMC − − 11160 − − 13206 AEB 9 . 7 9 . 7 324 6 . 2 2 . 6 1841 CD A-NL 5 . 7 9 . 0 364 6 . 4 2 . 7 1588 CD A-EV 3 . 8 5 . 57 519 6 . 0 2 . 5 1277 CD A-ME 2 . 9 3 . 25 259 6 . 2 2 . 6 91 July 20, 2016 DRAFT 36 T ABLE VI C H A R A C T E R I S T I C S O F T H E P R O P O S E D M O D E L S . “ P O S .” S TAN D S F O R P O S I T I V I T Y , “ I L L U M I N .” F O R I L L U M I N A T I O N , “ C O R R E L .” F O R C O R R E L ATI O N , ( + + ) B E S T R E S U LT S , A N D ( + ) G O O D R E S U LT S . Residuals Illumin. T ime Pos. Spatial Spectral coeff. c NL X Correl. X X + energies EV X Correl. Correl. X + values values ME X Correl. Correl. X ++ energies values July 20, 2016 DRAFT