IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. ?, NO. ?, 1 Implementation strategies for hyperspectral unmixing using Bayesian source separation Fr´ed´eric Schmidt, Albrecht Schmidt, Erwan Tr´eguier, Ma¨el Guiheneuf, Sa¨ıd Moussaoui and Nicolas Dobigeon, Member, IEEE, I. INTRODUCTION In visible and near infrared hyperspectral imaging, each image recorded by the sensor is the solar light reflected and diffused back from the observed planet surface and atmo- sphere at a particular spectral band. Under some assumptions related to surface and atmosphere properties – e.g., Lam- bertian surface, no intimate mixture, no diffusion terms in the atmosphere, homogeneous geometry in the scene – each measured spectrum – i.e., each pixel of the observed image for several spectral bands – is modeled as a linear mixture of the scene component spectra ( endmembers ) [1]–[3]. In this model, the weight of each component spectrum is linked to its abundance in the surface area which corresponds to the underlying pixel. The main goal of hyperspectral unmixing is to identify the components of the imaged surface and to estimate their respective abundances [4], [5]. By considering P pixels of an hyperspectral image acquired in L frequency bands, the observed spectra are gathered Fr´ed´eric Schmidt was with ESA, ESAC, Villanueva de la Canada, Madrid, Spain. He is now with University Paris-Sud, Laboratoire IDES, UMR8148, Orsay, F-91405, CNRS, Orsay, F-91405 (e-mail: frederic.schmidt@u-psud.fr). Albrecht Schmidt, Erwan Tr´eguier and Ma¨el Guiheneuf are with ESA, ESAC, Villanueva de la Canada, Madrid, Spain (e-mail: Albrecht.Schmidt@esa.int, mael.guiheneuf@gmail.com, erwan.treguier@sciops.esa.int). Sa¨ıd Moussaoui is with IRCCyN, UMR CNRS 6597, ´ Ecole Centrale Nantes, France (e-mail: said.moussaoui@irccyn.ec-nantes.fr). Nicolas Dobigeon is with the University of Toulouse, IRIT/INP- ENSEEIHT/T´eSA, 2 rue Charles Camichel, BP 7122, 31071 Toulouse cedex 7, France (e-mail: Nicolas.Dobigeon@enseeiht.fr). in a P × L data matrix X , potentially ignoring spatiality. Each row of this matrix contains a measured spectrum at a pixel with spatial index p = 1 , . . . , P . According to the linear mixing model, the p th spectrum, 1 ≤ p ≤ P , can be expressed as a linear combination of R pure spectra of the surface components. Using matrix notations, this linear spectral mixing model can be written as X ≈ AS (1) where non-negative matrices A ∈ R P × R + and S ∈ R R × L + approximate X ∈ R P × L + in the sense that 1 $_{2}$| | AS − X || 2 is minimized ($_{R}$ ·×· + denotes the space of matrices with only non-negative entries of respective dimensions). The rows of matrix S now contain the pure surface spectra of the R components, and each element a$_{pr}$ of matrix A corresponds to the abundance of the r th component in pixel with spatial index p . For a qualitative and quantitative description of the observed scene composition, the estimation problem consists of finding matrices S and A that allow one to explain the data matrix X and have a coherent physical interpretation. This approach casts the hyperspectral unmixing as a source separation prob- lem under a linear instantaneous mixing model [6]. Source separation is a statistical multivariate data processing problem whose aim is to recover unknown signals (called sources ) from noisy and mixed observations of these sources [7], [8]. This problem has been studied in-depth in recent years, starting with pioneer work more than 15 years ago [9], [10]. From a statistical point of view, the problem is also related to Principal Component Analysis (PCA) and k-means clustering (see [11] for an overview). Also note that the factorization AS is not uniquely defined. For instance, for any matrices Z ∈ R R × . + such that ZZ − 1 = I , then AZZ − $^{1}$S = ( AZ )( Z − $^{1}$S ) = A $^{′}$S′ is a solution as well; this holds even if the minimization is able to find a global minimum. However, when solving this separation problem with hyperspectral data, several constraints can be considered to reduce the set of admissible solutions. A first hard constraint is the non-negativity of the elements of both matrices S and A since they correspond to pure spectra and abundances of the surface components, respectively. A second constraint that may be imposed is the sum-to-one (additivity) constraint of the abundances. Indeed the abundance weights correspond to proportions and therefore should sum to unity. Several algorithms have been proposed in the literature to solve fully constrained unmixing problems, i.e., handling both of the constraints imposed on the spectra and abundances. Specifically, an iterative algorithm called ICE (Iterated Con- Abstract —Bayesian Positive Source Separation (BPSS) is a useful unsupervised approach for hyperspectral data unmixing, where numerical non-negativity of spectra and abundances has to be ensured, such in remote sensing. Moreover, it is sensible to impose a sum-to-one (full additivity) constraint to the estimated source abundances in each pixel. Even though non-negativity and full additivity are two necessary properties to get physically interpretable results, the use of BPSS algorithms has been so far limited by high computation time and large memory requirements due to the Markov chain Monte Carlo calculations. An implementation strategy which allows one to apply these algorithms on a full hyperspectral image, as typical in Earth and Planetary Science, is introduced. Effects of pixel selection, the impact of such sampling on the relevance of the estimated component spectra and abundance maps, as well as on the computation times, are discussed. For that purpose, two different dataset have been used: a synthetic one and a real hyperspectral image from Mars. Index Terms —Hyperspectral imaging, source separation, Bayesian estimation, implementation strategy, computation time. arXiv:1001.0499v2 [astro-ph.EP] 22 Jul 2010