Nonlinear unmixing of hyperspectral images: models and algorithms

1 Nonlinear Unmixing of Hyperspectral Images: Models and Algorithms Nicolas Dobigeon (1) , Jean-Yves T ourneret (1) , C ´ edric Richard (2) , Jos ´ e C. M. Bermudez (3) , Stephen McLaughlin (4) , and Alfred O. Hero (5) (1) Uni versity of T oulouse, IRIT/INP-ENSEEIHT/T ´ eSA, T oulouse, France { Nicolas.Dobigeon,Jean-Yves.Tourneret } @enseeiht.fr (2) Uni versit ´ e de Nice Sophia-Antipolis, OCA, Laboratoire Lagrange, Nice Cede x 2, France cedric.richard@unice.fr (3) Uni versity of Santa Catarina, Florian ´ opolis, SC, Brazil j.bermudez@ieee.org (4) School of Engineering and Physical Sciences, Heriot-W att Univ ersity , Edinbur gh, U.K. s.mclaughlin@hw.ac.uk (5) Uni versity of Michigan, Department of EECS, Ann Arbor , USA. hero@umich.edu Abstract When considering the problem of unmixing hyperspectral images, most of the literature in the geoscience and image processing areas relies on the widely used linear mixing model (LMM). Howe ver , the LMM may be not valid and other nonlinear models need to be considered, for instance, when there are multi-scattering effects or intimate interactions. Consequently , over the last few years, sev eral significant contributions hav e been proposed to ov ercome the limitations inherent in the LMM. In this paper , we present an overvie w of recent advances in nonlinear unmixing modeling. 2 I . M OT I V A T I O N F O R N O N L I N E A R M O D E L S Spectral unmixing (SU) is widely used for analyzing hyperspectral data arising in areas such as: remote sensing, planetary science chemometrics, materials science and other areas of micro- spectroscopy . SU provides a comprehensi ve and quantitative mapping of the elementary materials that are present in the acquired data. More precisely , SU can identify the spectral signatures of these materials (usually called endmembers ) and can estimate their relativ e contributions (or abundances ) to the measured spectra. Similar to other blind source separation tasks, the SU problem is naturally ill-posed and admits a wide range of admissible solutions. As a consequence, SU is a challenging problem that has received considerable attention in the remote sensing, signal and image processing communities [1]. Hyperspectral data analysis can be supervised , when the endmembers are known, or unsupervised , when they are unkno wn. Irrespectiv e of the case, most SU approaches require the definition of the mixing model underlying the observations. A mixing model describes in an analytical fashion how the endmembers combine to form the mixed spectrum measured by the sensor . The abundances parametrize the model. Gi ven the mixing model, SU boils down to estimating the in v erse of this formation process to infer the quantities of interest, namely the endmembers and/or the abundances, from the collected spectra. Unfortunately , defining the direct observation model that links these meaningful quantities to the measured data is a non-trivial issue, and requires a thorough understanding of complex physical phenomena. A model based on radiativ e transfer (R T) could accurately describe the light scattering by the materials in the observed scene [2], but would lead to very complex unmixing problems. Fortunately , in voking simplifying assumptions can lead to exploitable mixing models. When the mixing scale is macroscopic and each photon reaching the sensor has interacted with just one material, the measured spectrum y p ∈ R L in the p th pixel can be accurately described by the linear mixing model y p = R X r =1 a r,p m r + n p (1) where L is the number of spectral bands, R is the number of endmembers present in the image, m r is the spectral signatures of the r th endmember , a r,p is the abundance of the r th material in the p th pixel and n p is an additiv e term associated with the measurement noise and the modeling error . The abundances can be interpreted as the relativ e areas occupied by the materials in a giv en image pixel [3]. Thus it is natural to consider additional constraints regarding the abundance coef ficients a r,p    a r,p ≥ 0 , ∀ p, ∀ r P R r =1 a r,p = 1 , ∀ p (2) 3 In that case, SU can be formulated as a constrained blind source separation problem, or constrained linear regression, depending on the prior knowledge a v ailable regarding the endmember spectra. Due to the relati ve simplicity of the model and the straightforward interpretation of the analysis results, LMM-based unmixing strategies predominate in the literature. All of these techniques have been shown to be very useful whene ver the LMM represents a good approximation to the actual mixing. There are, howe ver , practical situations in which the LMM is not a suitable approximation [1]. As an illustrative example, consider a real hyperspectral image, composed of L = 160 spectral bands from the visible to near infrared, acquired in 2010 by the airborne Hyspex hyperspectral sensor over V illelongue, France. This image, with a spatial resolution of 0 . 5 m, is represented in Fig. 1 (top). From primary inspection and prior knowledge coming from av ailable ground truth, the 50 × 50 pixel region of interest depicted in Fig. 1 (bottom, right) is kno wn to be composed of mainly R = 3 macroscopic components (oak tree, chestnut tree and an additional non-planted-tree component). When considering the LMM to model the interactions between these R = 3 components, all the observed pixels should lie in a 2 -dimensional linear subspace, that can be easily identified by a standard principal component analysis (PCA). Conv ersely , if nonlinear effects are present in the considered scene, the observed data may belong to a 2 -dimensional nonlinear manifold. In that case, more complex nonlinear dimension reduction procedures need to be considered. The accuracy of these dimension reduction procedures in representing the dataset into a 2 -dimensional linear or nonlinear subspace can be ev aluated thanks to the av erage reconstruction error (ARE), defined as ARE = v u u t 1 LP N X n =1 k y n − ˆ y n k 2 (3) where y n are the observ ed pixels and ˆ y n the corresponding estimates. Here we contrast two ap- proaches, a locally linear Gaussian process latent variable model (LL-GPL VM) introduced in [4] and PCA. When using PCA to represent the data, the obtained ARE is 8 . 4 × 10 − 3 while using the LL- GPL VM, the ARE is reduced to 7 . 9 × 10 − 3 . This demonstrates that the inv estigated dataset should be preferably represented in a nonlinear subspace, as clearly demonstrated in Fig. 1 (bottom, left), where the nonlinear simplex identified by the fully constrained LL-GPL VM has been represented as blue lines. Consequently , more comple x mixing models need to be considered to cope with nonlinear interac- tions. These models are e xpected to capture important nonlinear ef fects that are inherent characteristics of hyperspectral images in se veral applications. They hav e proven essential to un veil meaningful information for the geoscience community [5]–[10]. Sev eral approximations to the R T theory hav e been proposed, such as Hapke’ s bidirectional model [3]. Unfortunately , these models require highly non-linear and inte gral formulations that hinder practical implementations of unmixing techniques. T o ov ercome these difficulties, se veral physics-based approximations of Hapke’ s model hav e been 4 Fig. 1. T op: real hyperspectral Madonna data acquired by the Hyspex hyperspectral scanner ov er V illelongue, France. Bottom, right: the region of interest shown in true colors. Bottom, left: Representation of the N = 2500 pixels (black dots) of the data and boundaries of the estimated nonlinear simplex (blue lines). proposed, mainly in the spectroscopy literature (e.g., see [3]). Howe v er , despite their wide interest, these approximations still remain difficult to apply for automated hyperspectral imaging. In particular, for such models, there is no unsupervised nonlinear unmixing algorithm able to jointly extract the endmembers from the data and estimate their relativ e proportions in the pixels. Meanwhile, sev eral approximate but exploitable nonlinear mixing models have been recently proposed in the remote sensing and image processing literatures. Some of them are similarly moti v ated by ph ysical ar guments, such as the class of bilinear models introduced later . Others exploit a more flexible non-linear math- ematical model to improve unmixing performance. Developing effecti ve unmixing algorithms based on nonlinear mixing models represents a challenge for the signal and image processing community . Supervised and unsupervised algorithms need to be designed to cope with nonlinear transformations that can be partially or totally unknown. Solving the nonlinear unmixing problem requires innov ative approaches to existing signal processing techniques. 5 (a) (b) (c) Fig. 2. (a) Linear mixing model: the imaged pixel is composed of two materials. (b) Intimate mixture: the imaged pixel is composed of a microscopic mixture of sev eral constituents. (c) Bilinear model: the imaged pixel is composed of two endmembers, namely tree and soil. In addition to the individual contribution of each material, bilinear interactions between the tree and the soil reach the sensor . More than 10 years after Kesha v a and Mustard’ s comprehensiv e revie w article on spectral unmixing [11], this article provides an updated revie w focusing on non-linear unmixing techniques introduced in the past decade. In [11], the problem on nonlinear mixtures was thoroughly addressed but, at that time, very few algorithmic solutions were av ailable. Capitalizing on almost one decade of advances in solving the linear unmixing problem, scientists from the signal and image processing communities hav e developed, and continue to do so, automated tools to extract endmembers from nonlinear mixtures, and to quantify their proportions in nonlinearly mixed pixels. The paper is organized as follo ws. The principal nonlinear mixing models are presented in the next section. Then the most popular nonlinear unmixing algorithms are revie wed. Model-based and model-free algorithms are considered. Existing solutions for supervised and unsupervised unmixing are also discussed. At the end of this paper , we present some recent strategies for detection of nonlinear mixtures in hyperspectral data. Finally , challenges and future directions for hyperspectral unmixing are reported in the conclusions. I I . N O N L I N E A R M O D E L S In [1], it is explained that linear mixtures are reasonable when two assumptions are wholly fulfilled. First the mixing process must occur at a macroscopic scale [12]. Secondly , the photons that reach the sensor must interact with only one material, as is the case in checkerboard type scenes [13]. An illustration of this model is depicted in Fig. 2(a) for a scene composed of two materials. When one of these two assumptions does not hold, dif ferent nonlinear effects may occur . T wo families of nonlinear models are described in what follows. A. Intimate mixtures The first assumption for linear mixtures is a macroscopic mixing scale. Howe ver , there are common situations when interactions occur at a microscopic lev el. The spatial scales in volv ed are typically 6 smaller than the path length follo wed by the photons. The materials are said to be intimately mixed [3]. Such mixtures ha ve been observ ed and studied for some time, e.g., for imaged scenes composed of sand or mineral mixtures [14]. The y hav e been advocated for analyzing mixtures observed in laboratory [15]. Based on R T theory , sev eral theoretical frame works ha ve been deriv ed to accurately describe the interactions suffered by the light when encountering a surface composed of particles. An illustration of these interactions is represented in Fig. 2(b). Probably the most popular ap- proaches dealing with intimate mixtures are those of Hapke in [3] since the y in volv e meaningful and interpretable quantities that hav e physical significance. Based on these concepts, se veral simplified nonlinear mixing models have been proposed to relate the measurements to some physical character- istics of the endmembers and to their corresponding abundances (that are associated with the relativ e mass fractions for intimate mixtures). In [2], the author deri ves an analytical model to express the measured reflectances as a function of parameters intrinsic to the mixtures, e.g., the mass fractions, the characteristics of the indi vidual particles (density , size) and the single-scattering albedo. Other popular approximating models include the discrete-dipole approximation [16] and the Shkuratov’ s model [17] (interested readers are in vited to consult [3] or the more signal processing-oriented papers [18], [19]). Howe v er these models also strongly depend on parameters inherent to the experiment since it requires the perfect knowledge of the geometric positioning of the sensor with respect to the observed sample. This dependency upon external parameters makes the in version (i.e., the estimation of the mass fractions from the collected spectra) very dif ficult to implement and, obvious ly , even more challenging in a unsupervised scenario, i.e., when the spectral signatures of the materials are unkno wn and need to be also recovered. More generally , it is worth noting that the first requirement of having a macroscopic mixing scale is intrinsically related to the definition of the endmembers. Indeed, defining a pure material requires specification of the spatial or spectral resolution, which is application dependent. Consider , for example, a simple scene composed of 3 materials A , B and C . It is natural to e xpect retriev al of these components individually when analyzing the scene. Howe ver , in other circumstances, one may be interested in the material components themselves, for instance, A 1 , A 2 , B 1 , B 2 , C 1 and C 2 if we assume that each material is composed of 2 constituents. In that case, pairs of subcomponents combine and, by performing unmixing, one might also be interested in reco v ering each of these 6 components. Con v ersely , it may be well known that the material A can never be present in the observed scene without the material B . In such case, unmixing would consist of identifying the couple A + B and C , without distinguishing the subcomponent A from the subcomponent B . This issue is frequently encountered in automated spectral unmixing. In each scenario, it is clear that when more details are desired, the mixtures should not occur at a macroscopic scale. T o circumvent this dif ficulty in defining the mixture scale, it makes sense to associate pure components with indi vidual instances 7 whose resolutions hav e the same order of magnitude than the sensor resolution. F or example, a patch of sand of spatially homogeneous composition can be considered as a unique pure component. In that case, most of the interactions occurring in most of the scenes of interest can be reasonably assumed to occur at a macroscopic level, at least when analyzing airborne and spaceborne remotely sensed images. B. Bilinear models Another type of nonlinear interaction occurs at a macroscopic scale, in particular in so-called multilayer ed configurations. One may encounter this nonlinear model when the light scattered by a giv en material reflects off other materials before reaching the sensor . This is often the case for scenes acquired ov er forested areas, where there may be many interactions between the ground and the canopy . An archetypal example of this kind of scene is sho wn in Fig. 2(c). Se veral models hav e been proposed to analytically describe these interactions. They consist of including powers of products of reflectance. Howe v er they are usually employed such that interactions of orders greater than two are neglected. The resulting models are known as the family of the bilinear mixing models. Mathematically , for most of these bilinear models, the observed spectrum y p ∈ R L in L spectral bands for the i th pixel is approximated by the following expansion y p = R X r =1 a r,p m r + R − 1 X i =1 R X j = i +1 β i,j,p m i  m j + n p . (4) where  stands for the termwise (Hadamard) product m i  m j =      m 1 ,i . . . m L,i            m 1 ,j . . . m L,j      =      m 1 ,i m 1 ,j . . . m L,i m L,j      In the right-hand side of (4), the first term, also found in (1), summarizes the linear contribution in the mixture while the second term models nonlinear interactions between the materials. The coefficient β i,j,p adjusts the amount of nonlinearities between the components m i and m j in the p th pixel. Se veral alternati ves for imposing constraints on these nonlinear coef ficients hav e been suggested. Similarly to [10], Nascimento and Dias assume in [20] that the (linear) ab undance and nonlinearity coefficients obey          a r,p ≥ 0 , ∀ p, ∀ r β i,j,p ≥ 0 , ∀ p, ∀ i 6 = j P R r =1 a r,p + P R − 1 i =1 P R j = i +1 β i,j,p = 1 . (5) It is worth noting that, from (5), this Nascimento model (NM), also used in [21], can be interpreted as a linear mixing model with additional virtual endmembers. Indeed, considering m i  m j as a pure 8 component spectral signature with corresponding abundance β i,j,p , the model in (5) can be rewritten y p = ˜ R X s =1 ˜ a s,p ˜ m s + n p with the positivity and additivity constraints in (2) where    ˜ a s,p , a r,p , ˜ m s , m r s = 1 , . . . , R ˜ a s,p , β i,j,p , ˜ m s , m i  m j s = R + 1 , . . . , ˜ R and ˜ R = 1 2 R ( R + 1) . This NM reduces to the LMM when ˜ a s,p = 0 for s = R + 1 , . . . , ˜ R . Con v ersely , in [9], Fan and his co-authors hav e fixed the nonlinearity coefficients as functions of the (linear) abundance coefficients themselves: β i,j,p = a i,p a j,p ( i 6 = j ). The resulting model, called the Fan Model (FM) in what follows, is thus fully described by the mixing equation y p = R X r =1 a r,p m r + R − 1 X i =1 R X j = i +1 a i,p a j,p m i  m j + n p (6) subject to the constraints in (2). One ar gument to explain the direct relation between the abundances and the nonlinearity coefficients is the following: if the i th endmember is absent in the p th pixel, then a i,p = 0 and there are no interactions between m i and the other materials m j ( j 6 = i ). More generally , it is quite natural to assume that the quantity of nonlinear interactions in a given pixel between two materials is directly related to the quantity of each material present in that pixel. Howe ver , it is clear that this model does not generalize the LMM, which can be a restrictive property . More recently , to alleviate this issue, the generalized bilinear model (GBM) has been proposed in [22] by setting β i,j,p = γ i,j,p a i,p a j,p y p = R X r =1 a r,p m r + R − 1 X i =1 R X j = i +1 γ i,j,p a i,p a j,p m i  m j + n p . (7) where the interaction coefficient γ i,j,p ∈ (0 , 1) quantifies the nonlinear interaction between the spectral components m i and m j . This model has the same interesting characteristic as the FM: the amount of nonlinear interactions is governed by the presence of the endmembers that linearly interact. In particular , if an endmember is absent in a pixel, there is no nonlinear interaction supporting this endmember . Ho we ver , it also has the significant adv antage of generalizing both the LMM when γ i,j,p = 0 and the FM when γ i,j,p = 1 . Having γ i,j,p > 0 indicate that only constructiv e interactions are considered. For illustration, synthetic mixtures of R = 3 spectral components have been randomly generated according to the LMM, NM, FM and GBM. The resulting data set are represented in the space spanned by the three principal eigen vectors (associated with the three lar gest eigen values of the sample cov ariance matrix of the data) identified by a principal component analysis in Fig. 3. These plots illustrate an interesting property for the considered dataset: the spectral signatures of the pure 9 components are still extremal points, i.e., vertices of the clusters, in the cases of FM and GBM mixtures contrary to the NM. In other words, geometrical endmember extraction algorithms (EEAs) and, in particular , those that are looking for the simplex of lar gest volume (see [23] for details), may still be valid for the FM and the GBM under the assumption of weak nonlinear interactions. All these bilinear models only include between-component interactions m i  m j with i 6 = j but no within-component interactions m i  m i . Finally , in [24], the authors deriv ed a nonlinear mixing model using a R T model applied to a simple canyon-lik e urban scene. Successive approximations and simplifying assumptions lead to the following linear -quadratic mixing model (LQM) y p = R X r =1 a r,p m r + R X i =1 R X j = i β i,j,p m i  m j + n p (8) with the positivity and additivity constraints in (2) and β i,j,p ∈ (0 , 1) . This model is similar to the general formulation of the bilinear models in (4), with the noticeable difference that the nonlinear contribution includes quadratic terms m i  m i . This contribution also shows that it is quite legitimate to include the termwise products m i  m j as additional components of the standard linear contrib ution, which is the core of the bilinear models described in this section. C. Other appr oximating physics-based models T o describe both macroscopic and microscopic mixtures, [25] introduces a dual model composed of two terms y p = R X r =1 a r,p m r + a R +1 ,p R R X r =1 f r,p w r ! + n p . The first term is similar to the one encountered in LMM and comes from the macroscopic mixing process. The second one, considered as an additional endmember with abundance a R +1 ,p , describes the intimate mixture by the a verage single-scattering albedo [2] e xpressed in the reflective domain by the mapping R ( · ) . On the use of geometrical LMM-based EEAs to identify nonlinearly mixed endmembers The first automated spectral unmixing algorithms, proposed in the 1990’ s, were based on geometrical concepts and were designed to identify endmembers as pure pixels (see [1] and [23] for comprehensiv e revie ws of geometrical linear unmixing methods). It is worth noting that this class of algorithms does not e xplicitly rely on the assumption of pixels coming from linear mixtures. They only search for endmembers as extremal points in the hyperspectral dataset. Provided there are pure pixels in the analyzed image, this might indicate that some of these geometrical approaches can be still valid for nonlinear mixtures that preserve this property , such as the GBM and the FM as illustrated in Fig. 3. 10 Fig. 3. Clusters of observ ations generated according to the LMM, the NM, the FM and the GBM (blue) and the corresponding endmembers (red crosses). Altmann et al. hav e proposed in [26] an approximating model able to describe a wide class of nonlinearities. This model is obtained by performing a second-order e xpansion of the nonlinearity defining the mixture. More precisely , the p th observed pixel spectrum is defined as a nonlinear transformation g p ( · ) of a linear mixture of the endmember spectra y p = g p R X r =1 a r,p m r ! + n p (9) where the nonlinear function g p is defined as a second order polynomial nonlinearity parameterized by the unique nonlinearity parameter b p g p : (0 , 1) L → R L x 7→  x 1 + b p x 2 1 , . . . , x L + b p x 2 L  T (10) This model can be rewritten y p = Ma p + b p ( Ma p )  ( Ma p ) + n p where M = [ m 1 , . . . , m R ] and a p = [ a 1 ,p , . . . , a R,p ] T . The parameter b p tunes the amount of nonlinearity present in the p th pixel of the image and this model reduces to the standard LMM when b p = 0 . It can be easily sho wn that this polynomial post-nonlinear model (PPNM) includes 11 bilinear terms m i  m j ( i 6 = j ) similar to those defining the FM, NM and GBM, as well as quadratic terms m i  m i similar to the LQM in (8). This PPNM has been shown to be sufficiently fle xible to describe most of the bilinear models introduced in this section [26]. D. Limitations a pixel-wise nonlinear SU Having revie wed the above physics-based models, an important remark must be made. It is impor- tant to note that these models do not take into account spatial interactions from materials present in the neighborhood of the targeted pixel. It means that these bilinear models only consider scattering ef fects in a giv en pixel induced by components that are present in this specific pixel. This is a strong simplifying assumption that allo ws the model parameters (abundance and nonlinear coef ficients) to be estimated pix el-by-pixel in the in version step. Note howe ver that the problem of taking adja- cency effects into account, i.e., nonlinear interactions coming from spectral interference caused by atmospheric scattering, has been addressed in an unmixing context in [27]. I I I . N O N L I N E A R U N M I X I N G A L G O R I T H M S Significant promising approaches hav e been proposed to nonlinearly unmix hyperspectral data. A wide class of nonlinear unmixing algorithms rely explicitly on a nonlinear physics-based parametric model, as detailed earlier . Others do not require definition of the mixing model and rely on very mild assumptions regarding the nonlinearities. For these two classes of approaches, unmixing algorithms hav e been considered under two different scenarios, namely supervised or unsupervised, depending on the a v ailable prior kno wledge on the endmembers. When the endmembers are known, supervised algorithms reduce to estimating the abundance coefficients in a single supervised in version step. In this case, the pure spectral signatures present in the scene must have been previously identified. For instance, they use prior information or suboptimal linear EEA. Indeed, as pre viously noted, when considering weakly nonlinearly mixed data, the LMM-based EEA may produce good endmember estimates when there are pure pixels in the dataset. In contrast, an unsupervised unmixing algorithm jointly estimates the endmembers and the abundances. Thus the unmixing problem becomes even more challenging, since a blind source separation problem must be solved. A. Model-based parametric nonlinear unmixing algorithms Gi ven a nonlinear parametric model, SU can be formulated as a constrained nonlinear regression or a nonlinear source separation problem, depending on whether the endmember spectral signatures are known or not. When dealing with intimate mixtures, some authors hav e proposed con v erting the measured reflectance into a single scattering albedo average; since this obeys a linear mixture, the mass fractions associated with each endmember can be estimated using a standard linear unmixing 12 algorithm. This is the approach adopted in [15] and [18] for known and unknown endmembers, re- specti vely . T o av oid the functional in version of the reflectance measurements into the single scattering albedo, a common approach is to use neural-networks (NN) to learn this nonlinear function. This is the strategy followed by Guilfoyle et al. in [28], for which sev eral improv ements have been proposed in [29] to reduce the computationally intensi v e learning step. In these NN-based approaches, the endmembers are assumed to be known a-priori, and are required to train the NN. Other NN-based algorithms have been studied in [30]–[33]. For the bilinear models introduced pre viously , supervised nonlinear optimization methods have been de veloped based on the assumption that the endmember matrix M is known. When the observed pixel spectrum y p is related to the parameters of interest θ p (a vector containing the abundance coefficients as well as any other nonlinearity parameters) through the function ϕ ( M , · ) , unmixing the pixel y p consists of solving the following minimization problem ˆ θ p = argmin θ k y p − ϕ ( M , θ ) k 2 2 . (11) This problem raises two major issues: i) the nonlinearity of the criterion resulting from the underlying nonlinear model ϕ ( · ) and ii) the constraints that have to be satisfied by the parameter vector θ . Since the NM can be interpreted as a linear mixing model with additional virtual endmembers, estimation of the parameters can be conducted with a linear optimization method as in [20]. In [9], [34] dedicated to FM and GBM, the authors propose to linearize the objectiv e criterion via a first-order T aylor series expansion of ϕ ( · ) . Then, the fully constrained least square (FCLS) algorithm of [35] can be used to estimate parameter vector θ . An alternativ e algorithmic scheme proposed in [34] consists of resorting to a gradient descent method where the step-size parameter is adjusted by a constrained line search procedure enforcing the constraints inherent to the mixing model. Another strategy initially introduced in [22] for the GBM is based on Monte Carlo approximations, dev eloped in a fully Bayesian statistical frame work. The Bayesian setting has the great advantage of providing a conv enient way to include the parameter constraints within the estimation problem, by defining appropriate priors for the parameters. This strategy has been also considered to unmix the PPNM [26]. When the spectral signatures M inv olv ed in these bilinear models need also to be identified in addition to the abundances and nonlinearity parameters, more ambitious unmixing algorithms need to be designed. In [36], the authors differentiate the NM to implement updating rules that generalize the SPICE algorithm introduced in [37] for the linear mixing model. Con versely , NMF-based iterati ve algorithms hav e been adv ocated in [38] for the GBM defined in (7), and in [24] for the LQM described in (8). More recently , an unsupervised version of the Bayesian PPNM-based unmixing algorithm initially introduced in [26] has been in vestigated in [39]. Adopting a geometrical point-of-view , Heylen and Scheunders propose in [40] an integral formula- 13 tion to compute geodesic distances on the nonlinear manifold induced by the GBM. The underlying idea is to deriv e an EEA that identifies the simplex of maximum volume contained in the manifold defined by the GBM-mixed pixels. B. Model-fr ee nonlinear unmixing algorithms When the nonlinearity defining the mixing is unkno wn, the SU problem becomes ev en more challenging. In such cases, when the endmember matrix M is fixed a priori, a classification approach can be adopted to estimate the ab undance coefficients which can be solved using support vector machines [41], [42]. Con versely , when the endmember signatures are not known, a geometrical- based unmixing technique can be used, based on graph-based approximate geodesic distances [43], or manifold learning techniques [44], [45]. Another promising approach is to use nonparametric methods based on reproducing kernels [46]–[51] or on Gaussian processes [4] to approximate the unkno wn nonlinearity . These two later techniques are described below . Nonlinear algorithms operating in reproducing kernel Hilbert spaces (RKHS) ha ve receiv ed con- siderable interest in the machine learning community , and have prov ed their efficienc y in solving nonlinear problems. K ernel-based methods ha ve been widely considered for detection and classifica- tion in hyperspectral images. Surprisingly , nonlinear unmixing approaches operating in RKHS hav e been inv estigated in a less in-depth way . The algorithms deriv ed in [46], [47] were mainly obtained by replacing each inner product between endmember spectra in the cost functions to be optimized by a kernel function. This can be vie wed as a nonlinear distortion map applied to the spectral signature of each material, independently of their interactions. This principle can be extremely efficient in solving detection and classification problems as a proper distortion can increase the detectability or separability of some patterns. It is howe ver of little physical interest in solving the unmixing problem because the nonlinear nature of the mixtures is not only governed by individual spectral distortions, but also by nonlinear interactions between the materials. In [48], a new kernel-based paradigm was proposed to take the nonlinear interactions of the endmembers into account, when these endmembers are assumed to be a priori known. It solves the optimization problem min ψ θ ∈H L X ` =1 [ y `,p − ψ θ ( m λ ` )] 2 + µ k ψ θ k 2 H (12) where m λ ` is the vector of the endmember signatures at the ` -th frequency band, namely , m λ ` = [ m `, 1 , . . . , m `,R ] T , with H a gi ven functional space, and µ a positive parameter that controls the trade-of f between regularity of the function ψ θ ( · ) and fitting. Again, θ is a vector containing the abundance coef ficients as well as any other nonlinearity parameters. It is interesting to note that (12) is the functional counterpart to (11), where ψ θ ( · ) defines the nonlinear interactions between the 14 endmembers assumed to be kno wn in [48]. Clearly , this strategy may fail if the functional space H is not chosen appropriately . A successful strategy is to define H as an RKHS in order to exploit the so- called kernel trick. Let κ ( · , · ) be the reproducing kernel of H . The RKHS H must be carefully selected via its kernel in order to make it flexible enough to capture wide classes of nonlinear relationships, and to reliably interpret a v ariety of experimental measurements. In order to extract the mixing ratios of the endmembers, the authors in [48] focus their attention on partially linear models, resulting in the so-called K-HYPE SU algorithm. More precisely , the function ψ θ ( · ) in problem (12) is defined by an LMM parameterized by the abundance vector a , combined with a nonparametric term, ψ θ ( m λ ` ) = a > m λ ` + ψ nlin ( m λ ` ) (13) possibly subject to the constraints in (2), where ψ nlin can be any real-valued function of an RKHS denoted by H nlin . This model generalizes the standard LMM, and mimics the PPNM when H nlin is defined to be the space of polynomial functions of de gree 2 . Remember that the latter is induced by the polynomial kernel κ ( m λ ` , m λ ` 0 ) = ( m > λ ` m λ ` 0 ) q of degree q = 2 . More complex interaction mechanisms can be considered by simply changing κ ( m λ ` , m λ ` 0 ) . By virtue of the reproducing kernel machinery , the problem can still be solved in the framew ork of (12). Another strategy introduced in [4] considers a kernel-based approach for unsupervised nonlinear SU based on a nonlinear dimensionality reduction using a Gaussian process latent variable model (GPL VM). In this work, the authors have used a particular form of kernel which extends the gen- eralized bilinear model in (7). The algorithm proposed in [4] is unsupervised in the sense that the endmembers contained in the image and the mixing model are not known. Only the number of endmembers is assumed to be known. As a consequence, the parameters to be estimated are the kernel parameters, the endmember spectra and the abundances for all image pixels. The main advantage of GPL VMs is their capacity to accurately model many different nonlinearities. GPL VMs construct a smooth mapping from the space of fractional abundances to the space of observed mixed pixels that preserves dissimilarities. This strategy has been also considered in [51] by Nguyen et al. , who solve the so-called pre-image problem [52] studied in the machine learning community . In the SU context, it means that pixels that are spectrally dif ferent hav e different latent v ariables and thus dif ferent abundance vectors. Howe ver , preserving local distances is also interesting: spectrally close pix els are expected to hav e similar abundance vectors and thus similar latent variables. Sev eral approaches hav e been proposed to preserve similarities, including back-constraints and locally linear embedding. For illustration, a small set of experiments has been conducted to ev aluate some of the model-based and model-free algorithms introduced above. First, 4 synthetic images of size 50 × 50 hav e been generated by mixing R = 3 endmember spectra (i.e., green grass, oliv e green paint and galvanized steel metal) extracted from the spectral libraries provided with the ENVI softw are [53]. These 4 15 Mixing models – with pure pixels Mixing models – w/o pure pixels LMM PPNM GBM FM LMM PPNM GBM FM Model-based algo. LMM N-FINDR + FCLS 1 . 42 14 . 1 7 . 71 13 . 4 3 . 78 13 . 2 6 . 83 9 . 53 unsupervised MCMC 0 . 64 12 . 4 5 . 71 8 . 14 0 . 66 10 . 9 4 . 21 3 . 92 PPNM Geodesic + GBA 1 . 52 10 . 3 6 . 04 12 . 1 4 . 18 6 . 04 4 . 13 3 . 74 unsupervised MCMC 0 . 39 0 . 73 1 . 32 2 . 14 0 . 37 0 . 81 1 . 38 2 . 25 GBM Geodesic + GBA 2 . 78 14 . 3 6 . 01 13 . 0 4 . 18 11 . 1 5 . 02 1 . 45 FM Geodesic + GBA 13 . 4 21 . 8 9 . 90 3 . 40 12 . 2 18 . 1 7 . 17 4 . 97 Geodesic + K-HYPE 2 . 43 9 . 71 5 . 23 11 . 3 2 . 44 5 . 92 3 . 18 2 . 58 T ABLE I A B U N D A N C E R N M S E S ( × 10 − 2 ) F O R V A R I O U S L I N E A R / N O N L I N E A R U N M I X I N G S C E N A R I O S . images hav e been generated according to the standard LMM (1), the GBM (7), the FM (6) and the PPNM (9), respectively . For each image, the abundance coef ficient vectors a p , [ a 1 ,p , . . . , a 3 ,p ] ( p = 1 , . . . , 2500 ) hav e been randomly and uniformly generated in the admissible set defined by the constraints (2). W e hav e also considered the more challenging scenario defined by the assumption that there is no pure pixel (by imposing a r,p < 0 . 9 , ∀ r , ∀ p ). The nonlinearity coefficients are uniformly drawn in the set [0 , 1] for the GBM. The PPNM-parameters b p , p = 1 . . . , P hav e been generated uniformly in the set [ − 0 . 3 , 0 . 3] . F or both scenario (i.e., with of without pure pix els), all images hav e been corrupted by an additiv e independent and identically distributed (i.i.d) Gaussian noise of v ariance σ 2 = 10 − 4 , which corresponds to an a verage signal-to-noise ratio of 20 dB (note that the usual SNR for most of the spectro-imagers are not belo w 30 dB). V arious unmixing strategies hav e been implemented to recov er the endmember signatures and then estimate the abundance coefficients. For supervised unmixing, the N-FINDR algorithm [54] and its nonlinear geodesic-based counterpart [43] hav e been used to extract the endmembers from linear and nonlinear mixtures, respectiv ely . Then, dedicated model-based strategies were used to recover the abundance fractions. The fully constrained least square (FCLS) algorithm [35] was used for linear mixtures. Gradient-based algorithms (GBA) were used for nonlinear mixtures. The GBA are detailed in [55], [34] and [9] for the PPNM, GBM and FM, respectiv ely . For comparison with supervised unmixing, and to ev aluate the impact of having no pure pixels in these images, joint estimations of endmembers and abundances was implemented using the Markov chain Monte Carlo techniques detailed in [56] and [39] for the LMM and PPNM images, respectively . Finally , the model-free supervised K-HYPE algorithm detailed in [48] was also coupled with the nonlinear EEA in [43]. The performance of these unmixing strategies has been e v aluated in term of abundance estimation error measured by RNMSE = v u u t 1 RP N X n =1 k a p − ˆ a p k 2 16 where a p is the n th actual abundance v ector and ˆ a p its corresponding estimate. The results are reported in T able I. These results clearly show that the prior knowledge of the actual mixing model underlying the observ ations is a clear adv antage for abundance estimation. Ho we ver , in the absence of such kno wledge, using an inappropriate model-based algorithm may lead to poor unmixing results. In such cases, as advocated before, PPNM seems to be sufficiently flexible to provide reasonable estimates, whatev er the mixing model may be. Otherwise, one may prefer to resort to model-free based strategy such as K-HYPE. I V . D E T E C T I N G N O N L I N E A R M I X T U R E S Consideration of nonlinear ef fects in hyperspectral images can provide more accurate results in terms of endmember and abundance identification. Howe ver , working with nonlinear models generally requires a higher computational complexity than approaches based on the LMM. Thus, unmixing linearly mixed pixels using nonlinear models should be av oided. Consequently , it is of interest to de vise techniques to detect nonlinearities in the mixing process before applying any unmixing method. Linearly mixed pixels can then be unmixed using linear unmixing techniques, leaving the application of more in v olved nonlinear unmixing methods to situations where the y are really necessary . This section describes approaches that ha ve been recently proposed to detect nonlinear mixing in hyperspectral images. A. Detection using a polynomial post-nonlinear model (PPNM) One interesting approach for nonlinearity detection is to assume a parametric nonlinear mixing model that can model dif ferent nonlinearities between the endmembers and the observations. A model that has been successfully applied to this end is the PPNM (9) studied in [26], [55]. PPNM assumes the post-nonlinear mixing described in (9) with the polynomial nonlinearity g p defined in (10). Hence, the nonlinearity is characterized by the parameter b p for each pixel in the scene. This parameter can be estimated in conjunction with the ab undance vector a p and the noise variance σ 2 . Denote as s 2 ( a p , b p , σ 2 ) the variance of the maximum likelihood estimator ˆ b p of b . Using the properties of the maximum likelihood estimator , it makes sense to approximate the distribution of ˆ b p by the following Gaussian distribution ˆ b p ∼ N  b p , s 2 ( a p , b p , σ 2 )  . The nonlinearity detection problem can be formulated as the following binary hypothesis testing problem    H 0 : y p is distributed according to the LMM (1) H 1 : y p is distributed according to the PPNM (9) (14) 17 Hypothesis H 0 is characterized by b p = 0 whereas nonlinear models ( H 1 ) correspond to b p 6 = 0 . Then, (14) can be rewritten as    H 0 : ˆ b p ∼ N (0 , s 2 0 ) H 1 : ˆ b p ∼ N ( b p , s 2 1 ) (15) where s 2 0 = s 2 ( a p , 0 , σ 2 ) and s 2 1 = s 2 ( a p , b p , σ 2 ) with b p 6 = 0 . Detection can be performed using the generalized likelihood ratio test. This test accepts H 1 (resp. H 0 ) if the ratio T , ˆ b 2 p /s 2 0 is greater (resp. lo wer) than a threshold η . As sho wn in [55], the statistic T is approximately normally distributed under the two hypotheses. Consequently , the threshold η can be explicitly related to the probability of false alarm (PF A) and the probability of detection (PD), i.e., the power of the test. Howe v er , this detection strategy assumes the prior knowledge of the v ariances s 2 0 and s 2 1 . In practical applications, Altmann et al. show that hav e proposed to modify the previous test strategy as follows [55] ˆ T = ˆ b 2 p ˆ s 2 0 H 1 ≷ H 0 η ∗ (16) where ˆ s 2 0 can be calculated as ˆ s 2 0 = CCRLB(0; ˆ a p , ˆ σ 2 ) (17) In (17), CCRLB is the constrained Cram ´ er-Rao lo wer-bound [57] on estimates of the parameter vector θ = [ a T p , b p , σ 2 ] T under H 0 , and  ˆ a p , ˆ σ 2  is the MLE of  a p , σ 2  . The performance of the resulting test is illustrated in Fig. 4 which shows the pixels detected as linear (red crosses) and nonlinear (blue dots) when generated according to various mixing models (LMM, FM, GBM and PPNM). B. Robust model-fr ee detection The detector discussed in the previous section assumes a specific nonlinear mixing model under the alternativ e hypothesis. Ho we v er , there are situations where the actual mixing does not obey any av ailable model. It is also possible that there is insufficient information to opt for any existing nonlinearity model. In these cases, it is interesting to address the problem of determining whether an observed pixel is a linear function of endmembers or results from a generic nonlinear mixing. One may consider the LMM (1) and the hyperplane P defined by P : ( z p     z p = Ma p , R X r =1 a r,p = 1 ) . (18) In the noise-free case, the hyperplane P lies in an ( R − 1) -dimensional subspace embedding all observ ations distrib uted according to the LMM. On the other hand, consider the general nonlinear mixing model y p = Ma p + µ p + n p (19) 18 Fig. 4. Pixels detected as linear (red crosses) and nonlinear (blue dotted) for the four subimages generated according the LMM, FM, GBM, and PPNM. Black lines depict the simplex corresponding to the noise-free case LMM. where µ p is an L × 1 deterministic vector that does not belong to P , i.e., µ p / ∈ P and a p satisfies the constraints (2). Note that a similar nonlinear mixing model coupled with a group-sparse constraint on µ p has been explicitly adopted in [58], [59] to make more robust the unmixing of hyperspectral pixels. In (19), µ p can be a nonlinear function of the endmember matrix M and/or the abundance vector a p and should be denoted as µ p ( M , a p ) [60]. Howe ver , the arguments M and a p are omitted here for bre vity . Giv en an observation vector y p , the detection of nonlinear mixtures can be formulated as the follo wing binary hypothesis testing problem    H 0 : y p is distributed according to the LMM (1) H 1 : y p is distributed according to the model (19) . Using the statistical properties of the noise n p , we obtain E[ y p |H 0 ] = Ma p ∈ P whereas E[ y p |H 1 ] = Ma p + µ p / ∈ P . As a consequence, it makes sense to consider the squared Euclidean distance δ 2 ( y p ) = min z p ∈P k y p − z p k 2 (20) between the observed pixel y p and the hyperplane P to decide which hypothesis ( H 0 or H 1 ) is true. As shown in [60], the test statistic δ 2 ( y p ) is distributed according to χ 2 distribution under the two hypotheses H 0 and H 1 . The parameters of this distribution depend on the known matrix M , the noise v ariance σ 2 and the nonlinearity vector µ p . If σ 2 is kno wn, the distribution of δ 2 ( y p ) is perfectly 19 kno wn under H 0 and partially known under H 1 . In this case, one may employ a statistical test that does not depend on µ p . This test accepts H 1 (resp. H 0 ) if the ratio T , δ 2 ( y p ) /σ 2 is greater (resp. lo wer) than a threshold η . As in the PPNM-based detection procedure, the threshold η can be related to the PF A and PD through closed-form expressions. In particular , it is interesting to note that the PD is intrinsically related to a non-Euclidean norm of the residual component µ p (see [60, Eq. (11)]), which is unfortunately unknown in most practical applications. If the noise variance σ 2 is unkno wn, which is the case in most practical applications, one can replace σ 2 with an estimate ˆ σ 2 , leading to T ∗ , δ 2 ( y p ) ˆ σ 2 H 1 ≷ H 0 η (21) where η is the threshold computed as pre viously indicated. The PF A and PD of the test (21) are then explicitly obtained using cumulativ e distribution functions of the χ 2 distribution. It was shown in [60] that the better the estimation of σ 2 , the closer the distributions of T and T ∗ and thus the closer the performances of the two corresponding tests. Se veral techniques can be used to estimate σ 2 . For instance, ˆ σ 2 has been estimated in [60] through an eigen-analysis of the sample cov ariance matrix of a set of pixels assumed to share the same v ariance. The v alue of ˆ σ 2 was determined as the av erage of the smallest eigenv alues of the sample cov ariance matrix. The accuracy of the estimator is a function of the number of eigen v alues considered. It was shown in [60] that a PF A smaller (resp. larger) than P ∗ F A is obtained if ˆ σ 2 > σ 2 (resp. ˆ σ 2 < σ 2 ). V . C O N C L U S I O N S A N D O P E N C H A L L E N G E S T o ov ercome the intrinsic limitations of the linear mixing model, se v eral recent contributions hav e been made for modeling of the ph ysical processes that underly h yperspectral observations. Some models attempt to account for between-material interactions affecting photons before they reach the spectro-imager . Based on these models, se veral parametric algorithms hav e been proposed to solve the resulting nonlinear unmixing problem. Another class of unmixing algorithms attempts to avoid the use of any rigid nonlinear model by using nonparametric machine learning-inspired techniques. The price to pay for handling nonlinear interactions induced by multiple scattering effects or intimate mixtures is the computational complexity and a possible de gradation of unmixing performance when processing large hyperspectral images. T o overcome these difficulties, one possible strategy consists of detecting pixels subjected to nonlinear mixtures in a pre-processing step. The pix els detected as linearly mixed can then benefit from the huge and reliable literature dedicated to the linear unmixing problem. The remaining pixels (detected as nonlinear) can then be the subject of particular attention. This paper has described dev elopments methods in nonlinear mixing for hyperspectral imaging. Se veral important se veral interesting challenges remain. First of all, better integration of algorithmic approaches and physical models ha ve the potential to greatly improv e non-linear unmixing perfor - 20 mance. By fully accounting for complex R T effects, such as scattering, dispersion, and beam interac- tion depth, a physical model can guide the choice of simplified mathematical and statistical models. Preliminary results hav e been recently communicated in [61], based on in situ measurements coupled with simulation tools (e.g., ray-tracing techniques). A second challenge is to develop unmixing models that take heterogeneity of the medium into account. Heterogeneous regions consist of combinations of linear , weakly non-linear and strongly non-linear pixels. The detection strategies detailed abov e might be one solution to tackle this problem since they are able to locate the areas where a non-linear model may outperform a linear model and vice-v ersa. Another approach adopted in [58], [59], which works well when there are only a few non-linear subregions, consists of using a statistical outlier approach to identify the non-linear pixels. Moreov er , as any nonlinear blind source separation problem, deriving flexible unsupervised unmixing algorithms is still a major challenge, especially if one wants to go one step further than a crude pixel-by-pixel analysis by exploiting spatial information inherent to these images. Finally , we observe that the presence of nonlinearity in the observed spectra is closely related to the number R of endmembers which is usually unkno wn. For e xample, in analogy to k ernelization in machine learning, after non-linear transformation, a nonlinear mixture of R components can often be represented as a linear mixture of ˜ R endmembers, with ˜ R > R . Recent advances in manifold learning and dimensionality estimation are promising approaches to the non-linear unmixing problem. A C K N O W L E D G M E N T S Part of this work has been funded by the Hypanema ANR Project n ◦ ANR-12-BS03-003 and the MADONN A project supported by INP T oulouse, France. Some results obtained in this paper result from fruitful discussions during the “CIMI W orkshop on Optimization and Statistics in Image Processing” (T oulouse, June 24-28, 2013). The authors are grateful to Y oann Altmann (Univ ersity of T oulouse) for sharing his numerous MA TLAB codes and for stimulating discussions. 21 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 ov erview: Geometrical, statistical, and sparse regression-based approaches, ” IEEE J. Sel. T opics Appl. Earth Observations and Remote Sens. , vol. 5, no. 2, pp. 354–379, April 2012. [2] B. W . Hapke, “Bidirectional reflectance spectroscopy . I. Theory , ” J. Geophys. Res. , vol. 86, pp. 3039–3054, 1981. [3] B. Hapke, Theory of Reflectance and Emittance Spectroscopy . Cambridge, UK: Cambridge Univ . Press, 1993. [4] Y . Altmann, N. Dobigeon, S. McLaughlin, and J.-Y . T ourneret, “Nonlinear spectral unmixing of hyperspectral images using Gaussian processes, ” IEEE T r ans. Signal Process. , vol. 61, no. 10, pp. 2442–2453, May 2013. [5] T . Ray and B. Murray , “Nonlinear spectral mixing in desert vegetation, ” Remote Sens. Envir onment , vol. 55, pp. 59–64, 1996. [6] J.-P . Combe, P . Launeau, V . Carr ` ere, D. Despan, V . M ´ l ` eder , L. Barill ´ e, and C. Sotin, “Mapping microphytobenthos biomass by non-linear inv ersion of visible-infrared hyperspectral images, ” Remote Sens. Envir onment , vol. 98, no. 4, pp. 371–387, 2005. [7] W . Liu and E. Y . Wu, “Comparison of non-linear mixture models: sub-pixel classification, ” Remote Sens. En vir onment , vol. 18, no. 2, pp. 145–154, Jan. 2005. [8] K. Arai, “Nonlinear mixture model of mixed pixels in remote sensing satellite images based on Monte Carlo simulation, ” Advances in Space Research , vol. 41, no. 11, pp. 1715–1723, 2008. [9] 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, ” Int. J. Remote Sens. , vol. 30, no. 11, pp. 2951–2962, June 2009. [10] B. Somers, K. Cools, S. Delalieux, J. Stuckens, D. V . der Zande, W . W . V erstraeten, and P . Coppin, “Nonlinear hyperspectral mixture analysis for tree cover estimates in orchards, ” Remote Sens. Env ir onment , vol. 113, pp. 1183– 1193, Feb. 2009. [11] N. Keshav a and J. F . Mustard, “Spectral unmixing, ” IEEE Signal Pr ocess. Mag. , vol. 19, no. 1, pp. 44–57, Jan. 2002. [12] R. B. Singer and T . B. McCord, “Mars: Large scale mixing of bright and dark surface materials and implications for analysis of spectral reflectance, ” in Pr oc. 10th Lunar and Planetary Sci. Conf. , March 1979, pp. 1835–1848. [13] R. N. Clark and T . L. Roush, “Reflectance spectroscop y: Quantitati v e analysis techniques for remote sensing applications, ” J. Geophys. Res. , vol. 89, no. 7, pp. 6329–6340, 1984. [14] D. B. Nash and J. E. Conel, “Spectral reflectance systematics for mixtures of powdered hypersthene, labradorite, and ilmenite, ” J. Geophys. Res. , vol. 79, pp. 1615–1621, 1974. [15] J. F . Mustard and C. M. Pieters, “Photometric phase functions of common geologic minerals and applications to quantitativ e analysis of mineral mixture reflectance spectra, ” J. Geophys. Res. , vol. 94, no. B10, pp. 13,619–13,634, Oct. 1989. [16] B. T . Draine, “The discrete-dipole approximation and its application to interstellar graphite grains, ” Astr ophysical Journal , vol. 333, pp. 848–872, Oct. 1988. [17] Y . Shkuratov , L. Starukhina, H. Hoffmann, and G. Arnold, “ A model of spectral albedo of particulate surfaces: implication to optical properties of the Moon, ” Icarus , vol. 137, no. 2, pp. 235–246, 1999. [18] J. M. P . Nascimento and J. M. Bioucas-Dias, “Unmixing hyperspectral intimate mixtures, ” in Proc. SPIE Image and Signal Pr ocessing for Remote Sensing XVI , L. Bruzzone, Ed., vol. 74830. SPIE, Oct. 2010, p. 78300C. [19] R. Close, P . Gader , A. Zare, J. W ilson, and D. Dranishnikov , “Endmember extraction using the physics-based multi- mixture pixel model, ” in Pr oc. SPIE Imaging Spectrometry XVII , S. S. Shen and P . E. Lewis, Eds., vol. 8515. San Diego, California, USA: SPIE, Aug. 2012, pp. 85 150L–14. 22 [20] J. M. P . Nascimento and J. M. Bioucas-Dias, “Nonlinear mixture model for hyperspectral unmixing, ” in Proc. SPIE Image and Signal Pr ocessing for Remote Sensing XV , L. Bruzzone, C. Notarnicola, and F . Posa, Eds., vol. 7477, no. 1. SPIE, 2009, p. 74770I. [21] N. Raksuntorn and Q. Du, “Nonlinear spectral mixture analysis for hyperspectral imagery in an unknown en vironment, ” IEEE Geosci. and Remote Sensing Lett. , vol. 7, no. 4, pp. 836–840, 2010. [22] A. Halimi, Y . Altmann, N. Dobigeon, and J.-Y . T ourneret, “Nonlinear unmixing of hyperspectral images using a generalized bilinear model, ” IEEE T r ans. Geosci. and Remote Sensing , vol. 49, no. 11, pp. 4153–4162, Nov . 2011. [23] W .-K. Ma, J. M. Bioucas-Dias, P . Gader , T .-H. Chan, N. Gillis, A. Plaza, A. Ambikapathi, and C.-Y . Chi, “Signal processing perspective on hyperspectral unmixing, ” IEEE Signal Process. Mag. , 2013, this issue. [24] I. Meganem, P . D ´ eliot, X. Briottet, Y . Deville, and S. Hosseini, “Linear -quadratic mixing model for reflectances in urban environments, ” IEEE T r ans. Geosci. and Remote Sensing , 2013, to appear . [25] R. Close, P . Gader, J. W ilson, and A. Zare, “Using physics-based macroscopic and microscopic mixture models for hyperspectral pixel unmixing, ” in Pr oc. SPIE Algorithms and T echnologies for Multispectral, Hyperspectral, and Ultraspectr al Imagery XVIII , S. S. Shen and P . E. Lewis, Eds., vol. 8390. Baltimore, Maryland, USA: SPIE, April 2012, pp. 83 901L–83 901L–13. [26] Y . Altmann, A. Halimi, N. Dobigeon, and J.-Y . T ourneret, “Supervised nonlinear spectral unmixing using a post- nonlinear mixing model for hyperspectral imagery , ” IEEE T rans. Image Process. , vol. 21, no. 6, pp. 3017–3025, June 2012. [27] D. Burazerovic, R. Heylen, B. Geens, S. Sterckx, and P . Scheunders, “Detecting the adjacency ef fect in hyperspectral imagery with spectral unmixing techniques, ” IEEE J. Sel. T opics Appl. Earth Observations and Remote Sens. , vol. 6, no. 3, pp. 1070–1078, June 2013. [28] K. J. Guilfoyle, M. L. Althouse, and C.-I. Chang, “ A quantitativ e and comparativ e analysis of linear and nonlinear spectral mixture models using radial basis function neural networks, ” IEEE T r ans. Geosci. and Remote Sensing , vol. 39, no. 8, pp. 2314–2318, Aug. 2001. [29] Y . Altmann, N. Dobigeon, S. McLaughlin, and J.-Y . T ourneret, “Nonlinear unmixing of hyperspectral images using radial basis functions and orthogonal least squares, ” in Pr oc. IEEE Int. Conf. Geosci. Remote Sens. (IGARSS) , V ancouver , Canada, July 2011, pp. 1151–1154. [30] J. Plaza, P . Martinez, R. Perez, and A. Plaza, “Nonlinear neural network mixture models for fractional abundance estimation in A VIRIS hyperspectral images, ” in Pr oc. XIII N ASA/Jet Pr opulsion Laboratory Airborne Earth Science W orkshop , Pasadena, CA, USA, 2004. [31] J. Plaza, A. Plaza, R. Perez, and P . Martinez, “ Automated generation of semi-labeled training samples for nonlinear neural network-based abundance estimation in hyperspectral data, ” in Pr oc. IEEE Int. Conf. Geosci. Remote Sens. (IGARSS) , 2005, pp. 345–350. [32] J. Plaza, A. Plaza, R. P ´ erez, and P . Martinez, “Joint linear/nonlinear spectral unmixing of hyperspectral image data, ” in Pr oc. IEEE Int. Conf. Geosci. Remote Sens. (IGARSS) , 2007, pp. 4037–4040. [33] J. Plaza, A. Plaza, Rosa, Perez, and P . Martinez, “On the use of small training sets for neural network-based characterization of mix ed pixels in remotely sensed hyperspectral images, ” P attern Recognition , v ol. 42, pp. 3032–3045, 2009. [34] A. Halimi, Y . Altmann, N. Dobigeon, and J.-Y . T ourneret, “Unmixing hyperspectral images using the generalized bilinear model, ” in Pr oc. IEEE Int. Conf. Geosci. Remote Sens. (IGARSS) , V ancouver , Canada, July 2011, pp. 1886– 1889. [35] 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. and Remote Sensing , vol. 29, no. 3, pp. 529–545, March 2001. 23 [36] P . Gader, D. Dranishnikov , A. Zare, and J. Chanussot, “ A sparsity promoting bilinear unmixing model, ” in Proc. IEEE GRSS W orkshop Hyperspectral Image SIgnal Pr ocess.: Evolution in Remote Sens. (WHISPERS) , Shanghai, China, June 2012. [37] A. Zare and P . Gader, “Sparsity promoting iterated constrained endmember detection in hyperspectral imagery , ” IEEE Geosci. and Remote Sensing Lett. , vol. 4, no. 3, pp. 446–450, July 2007. [38] N. Y oko ya, J. Chanussot, and A. Iwasaki, “Generalized bilinear model based nonlinear unmixing using semi- nonnegati v e matrix factorization, ” in Pr oc. IEEE Int. Conf. Geosci. Remote Sens. (IGARSS) , Munich, Germany , 2012, pp. 1365–1368. [39] 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 Process. , 2013, submitted. [40] R. Heylen and P . Scheunders, “Calculation of geodesic distances in nonlinear mixing models: Application to the generalized bilinear model, ” IEEE Geosci. and Remote Sensing Lett. , vol. 9, no. 4, pp. 644–648, July 2012. [41] J. Plaza, A. Plaza, P . Martinez, and R. Perez, “Nonlinear mixture models for analyzing laboratory simulated-forest hyperspectral data, ” Pr oc. SPIE Image and Signal Processing for Remote Sensing IX , vol. 5238, pp. 480–487, 2004. [42] P .-X. Li, B. W u, and L. Zhang, “ Abundance estimation from hyperspectral image based on probabilistic outputs of multi-class support vector machines, ” in Proc. IEEE Int. Conf. Geosci. Remote Sens. (IGARSS) , 2005, pp. 4315–4318. [43] R. Heylen, D. Burazerovic, and P . Scheunders, “Non-linear spectral unmixing by geodesic simplex volume maximiza- tion, ” IEEE J. Sel. T opics Signal Pr ocess. , vol. 5, no. 3, pp. 534–542, June 2011. [44] H. Nguyen, C. Richard, P . Honeine, and C. Theys, “Hyperspectral image unmixing using manifold learning methods. deriv ations and comparativ e tests, ” in Pr oc. IEEE Int. Conf. Geosci. Remote Sens. (IGARSS) , Munich, Germany , July 2012. [45] G. Licciardi, X. Ceamanos, S. Doute, and J. Chanussot, “Unsupervised nonlinear spectral unmixing by means of NLPCA applied to hyperspectral imagery , ” in Pr oc. IEEE Int. Conf. Geosci. Remote Sens. (IGARSS) , Munich, Germany , July 2012, pp. 1369–1372. [46] J. Broadwater , R. Chellappa, A. Banerjee, and P . Burlina, “Kernel fully constrained least squares abundance estimates, ” in Pr oc. IEEE Int. Conf. Geosci. Remote Sens. (IGARSS) , Barcelona, Spain, July 2007, pp. 4041–4044. [47] J. Broadwater and A. Banerjee, “ A comparison of kernel functions for intimate mixture models, ” in Pr oc. IEEE GRSS W orkshop Hyperspectral Image SIgnal Process.: Evolution in Remote Sens. (WHISPERS) , Aug. 2009, pp. 1–4. [48] 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. 60, no. 2, pp. 480–492, Jan. 2013. [49] ——, “Nonlinear abundance estimation of hyperspectral images with L1-norm spatial regularization, ” IEEE T rans. Geosci. and Remote Sensing , 2013, submitted. [50] X. Li, J. Cui, and L. Zhao, “Blind nonlinear hyperspectral unmixing based on constrained kernel nonnegati ve matrix factorization, ” Signal, Image and V ideo Pr ocessing , pp. 1–13, 2012. [51] N. H. Nguyen, J. Chen, C. Richard, P . Honeine, and C. Theys, “Supervised nonlinear unmixing of hyperspectral images using a pre-image methods, ” EAS Publications Series , vol. 59, pp. 417–437, 2013. [52] P . Honeine and C. Richard, “Preimage problem in kernel-based machine learning, ” IEEE Signal Processing Magazine , vol. 28, no. 2, pp. 77–88, March 2011. [53] RSI (Research Systems Inc.), ENVI User’ s guide V ersion 4.0 , Boulder , CO 80301 USA, Sept. 2003. [54] M. W inter , “Fast autonomous spectral end-member determination in hyperspectral data, ” in Proc. 13th Int. Conf. on Applied Geologic Remote Sensing , vol. 2, V ancouv er , April 1999, pp. 337–344. [55] Y . Altmann, N. Dobigeon, and J.-Y . T ourneret, “Nonlinearity detection in hyperspectral images using a polynomial post-nonlinear mixing model, ” IEEE T r ans. Image Pr ocess. , vol. 22, no. 4, pp. 1267–1276, April 2013. 24 [56] 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 r ans. Signal Process. , vol. 57, no. 11, pp. 4355–4368, Nov . 2009. [57] J. D. Gorman and A. O. Hero, “Lower bounds for parametric estimation with constraints, ” IEEE T r ans. Inf. Theory , vol. 36, no. 6, pp. 1285 –1301, Nov . 1990. [58] 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 r ans. Image Pr ocess. , 2013, submitted. [59] N. Dobigeon and C. F ´ evotte, “Robust nonne gati ve matrix factorization for nonlinear unmixing of h yperspectral images, ” in Pr oc. IEEE GRSS W orkshop Hyperspectral Imag e SIgnal Pr ocess.: Evolution in Remote Sens. (WHISPERS) , Gainesville, FL, June 2013. [60] Y . Altmann, N. Dobigeon, J.-Y . T ourneret, and J. C. M. Bermudez, “ A robust test for nonlinear mixture detection in hyperspectral images, ” in Proc. IEEE Int. Conf. Acoust., Speech, and Signal Pr ocessing (ICASSP) , V ancouver , Canada, June 2013, pp. 2149–2153. [61] L. Tits, W . Delabastita, B. Somers, J. Farifteh, and P . Coppin, “First results of quantifying nonlinear mixing effects in heterogeneous forests: A modeling approach, ” in Pr oc. IEEE Int. Conf. Geosci. Remote Sens. (IGARSS) , 2012, pp. 7185–7188.