Contributions to Hyperspectral Unmixing Sina Nakhostin 2017 Contributions to Hyperspectral Unmixing Sina Nakhostin 2017 Joint Anomaly detection and spectral Unmixing

Keywords: Hyperspectral Unmixing Sina Nakhostin 2017 Planetary Hyperspectral Unmixing, Anomaly Detection, Manifold Learning, Nonnegative Matrix Factorization, Overcomplete Dictionary, Kernel based Learning. Contributions to Hyperspectral Unmixing Sina Nakhostin 2017 4 Supervised Planetary Unmixing with Optimal Transport Optimal Transport, Spectral Unmixing, Endmember Variability, Wassertein Distance, Bregman Projection

Hyperspectral images are commonly used in the context of planetary exploration, especially for the analysis of the composition of planets. As several instruments have been sent throughout the Solar System, a huge quantity of data is getting available for the research community. Among classical problems in the analysis of hyperspectral images, a crucial one is unsupervised non-linear spectral unmixing, which aims at estimating the spectral signatures of elementary materials and determining their relative contribution at a sub-pixel level. While the unmixing problem is well studied for earth observation, some of the traditional problems encountered with earth images are somehow magnified in planetary exploration. Among them, large image sizes, strong non-linearities in the mixing (often different from those found in the earth images) and presence of anomalies are usually impairing the unmixing algorithms. This chapter presents a new method that scales favorably with the problem posed by this analysis. It performs an unsupervised unmixing jointly with anomaly detection capacities, and has a global linear complexity. Non linearities are handled by decomposing the hyperspectral data on an overcomplete set of spectra, combined with a specific sparse projection, which guarantees the interpretability of the analysis. A theoretical study is proposed on synthetic datasets, and results are presented over the challenging 4-Vesta asteroid dataset.

2.1

In the ideal case, where there is no effect of noise, also with the assumption of pure endmembers being present in the data, the abundance coefficients as the result of inversion part of linear unmixing will lay within a convex polygon (n-1 simplex) whose vertexes are consist of the endmembers. Illustration courtesy of Dobigeon et al. . .

2.2

From left to right: linear mixing model(LMM), intimate mixing model (IMM) and bilinear mixing model (BMM). Courtesy of [START_REF] Dobigeon | Nonlinear unmixing of hyperspectral images: Models and algorithms[END_REF]. . . . . . .

3.1

Overcomplete endmember dictionary for non-linear unmixing . . . . .

3.2

Illustration of difficulties in presence of outliers. First Row: endmember extraction using SAGA without anomalies: in this situation extracted endmembers are a reliable representation of the manifold ; Second Row: endmember extraction using SAGA in presence of anomalies: in this situation extracted endmembers (in red) do not correspond to real borders of the data manifold ; Third Row: endmember extraction and anomaly detection In presence of anomalies using SAGA+: here extracted endmembers (in red) by are in accordance with borders of the manifold embedding data while the anomalies (in cyan) are properly selected.

Evolution of SOPE (Sum Of Projection Errors

) when all data in green are represented through a single endmember uniformly represented on a meshgrid. One observes a decrease in SOPE when candidates get closer to the nominal distribution. . . . . . . . . . . .

3.12

Results on real 4-Vesta asteroid. Top: Vesta-4 mosaic with the spatial position of the extracted endmembers. The red points are detected anomalies and the 10 green points with labels are selected endmember prototypes. Bottom: The respective estimated endmember signatures. . . . . . . . . . . . . . . . . . . . . . . . . . .

3.13

Top: laboratory spectra of mixtures involving olivine and orthopyroxene (spectra from the RELAB library with numbers referring respectively to the olivine and orthopyroxene contents (% olivine / % orthopyroxene). Bottom: comparison with the three extracted endmembers P0, P5 and P9. The shift in absorption center at 2µm is due to an increase in Ca and/or Fe in the pyroxene, compared to the pure orthopyroxene composition. . . . . . . . . . . . . . . . . . . . .

3.15

Abundance maps and respective spectral signatures extracted by SAGA+ (with anomaly detection procedure). Using an anomaly detection procedure, one observes that signatures whit magnitudes higher than 1 and those with sharp peaks are not selected as endmembers. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

5.

12 Dictionary of the manually extracted endmembers. Based on expert observations, pixel vectors related to highly distinct materials within the scene have been extracted to be used as a dictionary of endmembers. In order to account for material variability, for each material two representatives endmembers were selected in order to build-up a dictionary of 24 signatures. . . . . . . . . . . . . . . . . . . . . . . . .

5.13

Using the exact simplex volume maximization (E-SiVM) procedure a dictionary of 24 atoms with maximum variability among their signature profiles has been generated. One can observe even larger variability among endmember signatures while comparing this dictionary to the one extracted manually. . . Hyperspectral (HS) remote sensing is the scientific field which is emerged as the synergy of airborne/satellite imagery and spectroscopy. Hyperspectral sensors are categorized as passive sensors. Unlike active sensors who perceive and measure the reflection of their own emitted signal, passive sensors perceive the reflected signal of another source. In case of hyperspectral remote sensing, this source is usually the natural light emitted from the sun and reflected by the surface or atmosphere of the astronomical bodies. Figure 1.1 depicts this difference.

Domains of application

The range of applications of these sensors is considerably large and concerns for example astronomy, geoscience, agriculture, archeology, surveillance and environ-ment monitoring where the HS acquisition senses each object not through the three RGB wavelength bands but through hundreds of bands covering large portions of the electromagnetic range, which is depicted in Figure 1.2. Another domain which in recent years has been highly benefited from the advances of hyperspectral images is the medical imaging. In the following a brief look is taken on the diverse domains where this technology can be used. In many applications therefore, hyperspectral data combines more than one captor.

HSI for Environmental monitoring

Desertification, land degradation in arid, semi-arid, and dry sub-humid regions, is emerging as a global environmental dilemma [START_REF] Yang | Desertification assessment in china: An overview[END_REF]. This has called for accurate assessment of the status, change, and trend of desertification for controlling this global threat. Identification, mapping and monitoring of the salt-affected areas have been exploiting remote sensing data during the recent decades [START_REF] Farifteh | Assessing salt-affected soils using remote sensing, solute modelling, and geophysics[END_REF][START_REF] Aghakouchak | Aral sea syndrome desiccates lake urmia: Call for action[END_REF]. In Figure 1.3 one can follow the chronological process of the shrinkage in the size of a salt lake. As other HSI applications in environmental monitoring one can name oil spill detection [START_REF] Brekke | Oil spill detection by satellite remote sensing[END_REF], water pollution monitoring [START_REF] Thiemann | Lake water quality monitoring using hyperspectral airborne data-a semiempirical multisensor and multitemporal approach for the mecklenburg lake district, germany[END_REF] and estimation of land surface temperature (LST) [START_REF] Weng | Estimation of land surface temperature-vegetation abundance relationship for urban heat island studies[END_REF]. In Figure 1.4 a heat map generated according the satellite remote sensing images is depicted. In chapter 5, a number of algorithms developed during this PhD work are going to be used for generation of urban surface abundance maps which by themselves help to estimate LST maps of urban areas.

HSI for medical analysis

During the recent years hyperspectral imagery has increasingly been used for the early detection of different sorts of cancer using noninvasive or minimally invasive techniques [START_REF] Akbari | Cancer detection using infrared hyperspectral imaging[END_REF]. In order to show the riveting role of HSI in shaping the future of medicine one can refer to a recent European project HELICoiD [START_REF] Salvador | Helicoid: Interdisciplinary and collaborative project for real-time brain cancer detection: Invited paper[END_REF]. The aim of HELICoiD was to demonstrate the use of hyperspectral imaging for real-time identification of tumor margins during surgery, helping the surgeon to extract the entire tumor and to spare as much of the healthy tissue as possible. In [START_REF] Kabwama | Intra-operative hyperspectral imaging for brain tumour detection and delineation: Current progress on the helicoid project[END_REF] and [16] for example new manifold embedding techniques are proposed for delineation and segmentation of the tumor effected tissues of brain, the result of which is presented in Figure 1.5. By exploiting the spectral depth in HS images of healthy and cancerous tissue, machine learning algorithms are being trained to detect differences in their spectral signature. The learned model then being applied to identify cancerous and healthy tissue in hyperspectral images acquired during surgery. It can be imagined that in the future, this information could be available for surgeon in real time by overlaying conventional images with a color-map that indicates the likelihood of that particular area being cancerous.

HSI in ART

The need to conserve historical artworks is not new. Many of the methods of art preservation which have developed during 19th century, have been seen to actually contribute to, rather than prevent, the deterioration of paintings. Today it is hoped that by using more advance scientific methods and technology, analysis of paintings, identification of conservation issues and detection of additions and forgeries become faster and more accurate. In general major problems related to art and painting preservation are listed as : Traditional scientific methods include chemical analysis of paint samples to determine the molecular nature of the pigments and varnishes. X-ray technology is one often utilized method for observation of beneath surface layers and discovering the intentions and processes of the artist. Despite their relative success, many of the methods used either give limited information, are invasive to some degree, or both. HSI turns to be a useful method in material identification and mapping on works of art and paintings [START_REF] Kubik | Chapter 5 hyperspectral imaging: A new technique for the non-invasive study of artworks[END_REF][START_REF] Mounier | Hyperspectral imaging, spectrofluorimetry, fors and xrf for the non-invasive study of medieval miniatures materials[END_REF][START_REF]Hyperspectral imaging applied to the analysis of goya paintings in the museum of zaragoza (spain)[END_REF]. By using hyperspectral imaging systems it has become possible to produce false color infrared images that can show restored zones of the painting and therefore help to guide the choice for the pigment analysis of significant zones of the original paintings. Also in a recent work [START_REF] Dooley | Mapping of egg yolk and animal skin glue paint binders in early renaissance paintings using near infrared reflectance imaging spectroscopy[END_REF] it showed possible to identify which pigment binders had been used in specific paintings (See Figure 1.6 for an illustration).

• Pigment Identification

HSI for Archeology

Another domain of application for hyperspectral images is Archeology. In a very similar fashion as to the Underdrawings of the old paintings, historical sites and landscapes might be palimpsests, or layered, while mixing of past and present occupation layers can obscure imagery analysis from space [START_REF] Sarah | Satellite remote sensing for archaeology[END_REF]. Archeological sites, ancient road systems etc. therefore could be well buried under tons of thousands of years of silt deposit or even under the current existing cities and towns. Remote sensing hyperspectral data could be used along other ground-based equipments in order to draw more accurate maps of old archaeological landscapes [START_REF] Rosa | Remote hyperspectral imagery as a support to archaeological prospection[END_REF]. In [START_REF] Ur | Corona satellite photography and ancient road networks: a northern mesopotamian case study[END_REF] an old road system in northern Mesopotamian has been partially detected using the remote sensing imagery.

HSI in Astrophysics

From planetary surface analysis [START_REF] Prakash Chauhan | Hyperspectral remote sensing of planetary surfaces: an insight into composition of inner planets and small bodies in the solar system[END_REF] to atmospheric studies to detection of galaxies, hyperspectral images are being widely used in different astrophysical studies. In [START_REF] Combe | Analysis of omega/mars express data hyperspectral data using a multiple-endmember linear spectral unmixing model (melsum): Methodology and first results[END_REF] and [START_REF] Le Mouelic | An iterative least squares approach to decorrelate minerals and ices contributions in hyperspectral images: Application to cuprite (earth) and mars[END_REF] the hyperspectral data of Observatoire pour la Minéralogie, l'eau, les glaces et l'activité (OMEGA) which is mount on board of MARS EXPRESS and has been monitoring the surface of the planet Mars since 2003 is analyzed. Using a Multiple-Endmember Linear Spectral Unmixing Model (MELSUM), the existence of different surface materials including CO 2 and H 2 O has been studied. For this a reference mineral library composed of 24 different mineral fingerprints (endmember signatures) were used. In Figure 1.7 The global atlas maps produced by this model are illustrated. As for the endmember selection process a systematic exploration of whole set of combinations of four components of the endmember signatures were deployed by manually excluding the maps whit negative values. Although the utilization of this linear method was positively evaluated, it is believed that more subtleties might be captured in case non-linear models of mixing were considered.

In [START_REF] Moussaoui | On the decomposition of mars hyperspectral data by ica and bayesian positive source separation[END_REF] independent component analysis (ICA) and Bayesian Positive Source Separation (BPSS) methods were exploited instead to fuse the spatial and spectral information of the same data and enhance the process of endmember estimation and abundance map estimation. The idea is to do a rough classification of pixels based on ICA which provides a reduced set of signatures, which by itself is being used by BPSS to do endmember estimation and abundance estimation. In another study [START_REF] Wh Bakker | Processing omega/mars express hyperspectral imagery from radiance-at-sensor to surface reflectance[END_REF] OMEGA data is being used to come up with a processing methodology for converting OMEGA data from radiance-at-sensor to surface reflectance, performing an atmospheric correction and filtering the systematic and random noise. Also the presence of iron-rich clays on some areas is reported within this work. As another example of the applications of HS remote sensing for planetary data analysis one can refer to NASA's DAWN satellite whose mission is to orbit around two important protoplanets of the asteroid belt (4-VESTA and CERES), for which it utilizes a sensor composed of 432 bands. Through these bands each material acquires a fingerprint which is called spectral signature of that material. In chapters 3 and 4 applications of the methods developed during this thesis will be presented on a hyperspectral image from the 4-VESTA asteroid.

Ground-based HSI

The Multi Unit Spectroscopic Explorer (MUSE) instrument installed at the Very Large Telescope (VLT) of the European Southern Observatory (ESO) is a ground based integral-field spectrograph which has combined very high spectral and spatial resolution characteristics in an unprecedented way into one instrument. In Figure 1.8 an illustration of the instrument together with an abstract data cube generated by MUSE is visible. Because of this unique property MUSE is mainly meant to be used for extra-galactic studies including globular clusters, planetary nebulae and Lyman-alpha emitters also supermassive black holes. [START_REF] Bourguignon | Sparsity-based spatial-spectral restoration of muse astrophysical hyperspectral data cubes[END_REF] and [START_REF] Shen | Spatial-spectral unmixing of hyperspectral data for detection and analysis of astrophysical sources with the muse instrument[END_REF] are among the first efforts in utilizing state-of-the-art spectral analysis procedures in detection and unmixing of MUSE data. In [START_REF] Diener | The muse-wide survey: detection of a clustering signal from lyman α emitters in the range 3< z< 6[END_REF] a clustering analysis of a sample of 238 Lyman α emitters is presented at redshift 3 ≤ z ≤ 6 from the MUSE-Wide survey. In [START_REF] Meillier | Selfi: an object-based, bayesian method for faint emission line source detection in muse deep field data cubes[END_REF] a new Bayesian method optimized for detection of faint galaxies in MUSE deep fields is presented. The list of publications as the result of research on MUSE data is considerably growing since its down in 2014. One can access an updating list of research related to MUSE in a dedicated web page1 .

Multispectral vs Hyperspectral: literal vs non-literal

Hyperspectral remote sensors can register images in the format of 3D data cubes. These data cubes are being characterized based on their spatial resolution, spectral resolution and their signal-to-noise ratio. The spatial resolution is the measure of the area being represented through one pixel (thus providing so-called literal information) and spectral resolution is the width of the electromagnetic field being represented through one band of the cube (therefore providing so-called non-literal information). The signal-to-noise ratio (SNR) reflects spectrometer's radiance response and relates to its acquired image quality. The large number of spectral bands in hyperspectral images (HSI) differentiates them from multispectral images (MSI) where the spectral representation of pixels never exceeds a few bands. This is illustrated in Figure 1.9. Therefore both thematics and methodologies associated with HS data are specific and different than the ones used with MSI. In the case of MSI, because of the very low spectral resolution, in practice no or very limited information w.r.t the existing subpixel materials is retrievable, therefore, except a few interesting attempts for MSI subpixel classification( [START_REF] Kopriva | Blind multispectral image decomposition by 3d nonnegative tensor factorization[END_REF] and [START_REF] Moussa Sofiane Karoui | Blind spatial unmixing of multispectral images: New methods combining sparse component analysis, clustering and non-negativity constraints[END_REF]), processing of these types of data has been historically focused mainly on the spatial information and spatial correlation among pixels. Spatial pattern classification is therefore one of the most studied problems that could be addressed using MSI and correlation-based methods which have been devised upon them. These are statistical inference methods where the independent and identically distributed (i.i.d) assumption does not hold upon pixels therefore, by somehow defining a measure of spatial distance between the pixels, one can assign n total pixels of the image to different classes. The measure of closeness is usually defined through a weighting matrix which then operates on the pixels in the vicinity of the current pixel under study. Thanks to this characteristic, correlation-based methods [START_REF] Gao | Hyperspectral Image Classification by Using Pixel Spatial Correlation[END_REF] have been extensively used for different spatial classification problems for MSI analysis.

With hyperspectral images, thanks to their rich spectral information, it is possible to assess the composition of the various materials embedded in a single pixel. Problems including unmixing associated with anomaly detection and endmember extraction are within this category.

HSI from the view of Pigeon-Hole Principle

In the pigeon-hole principle we have that for D places and objects, where 0 > D > , there will be at least one place with more than one occupying object. It is possible to see the relationship between different targets/endmembers and spectral bands of HSI through the pigeon-hole principle. Assume different HSI bands as the places with total number of D, and different targets of interest as the objects. In case of D > , one band can be used to detect one target. Since in HSI there are hundreds of contiguous bands available, we can therefore distinguish hundreds of spectrally distinguished targets. For this however, three requirements must meet: First) The number of spectral bands should be larger than the number of endmembers (D > ).

Second) Number of interesting targets/endmembers ( ) should be known and Third) No two distinguished endmember signatures to be accommodated by a single band.

Spectral Unmixing

The current thesis focuses mainly on one specific problem in the field of hyperspectral remote sensing called "Spectral Unmixing". Spectral Unmixing has been an active area of research during the last years and found its application in domains including but not limited to remote sensing, audio signal processing and chemistry. In audio signal processing the term "blind source separation" and in the chemical sciences "chemical agent detection" is used instead to refer to the same problem. Despite their very high spectral resolution, hyperspectral images are indeed known to be of relatively low spatial resolution. This low resolution is a relative notion and is due to technological limitations of the HSI captors, also is negatively proportionates with the distance of the captor from the scene of view. As a consequence the values of pixels within an image are likely to be mixtures of different materials in the scene. In Figure 1.10 one can observe the process of mixing of three materials within a sensor (contributing to one pixel of an image) acquisition. Hyperspectral Unmixing then can be defined as an inverse procedure that consists of identifying in each pixel the amount of pure elements that contributed to the pixel's mixture. The total number of pure elements (also called endmembers) and the number of pure elements included in one pixel are then two informations tricky to retrieve.

The simplest situation is when both the total number and type of existing materials within the scene are known, associated with a linear mixing process assumption. Many different unmixing algorithms rely on this prior assumption (detailed in chapter 2). Though efficient in some situations, this linearity assumption however does not necessarily hold in many real world scenarios. Due to many physical factors issued from the interactions of light and the surface materials before it reaches the captor, or characteristics of the imaging instrument, the registered pixels might not be linear mixtures of the existing components. Therefore a number of approaches has been considered to bring this non-linearity into account, yielding to more complex strategies. A major categorization among the unmixing algorithms can be explained as:

• Methods assuming to have a priori the endmembers. These methods are called Supervised unmixing and consists in estimating the abundances for each (known) endmember contributing to each pixel.

• Methods that do not know the endmembers. These Unsupervised unmixing approaches aim both at estimating endmembers and associated abundances in each pixel.

In addition to these issues, it is important to outline that in practice, we rarely have an exact spectral representation of the pure materials. Each pure material usually has several different representations depending on the conditions of acquisition, its precise composition, etc. Nuances among different spectral signatures of the same endmember can be either issued from the atmospheric condition in which signatures acquired or simply due to different levels of purity of the material.

Endmember Variability

Pure materials are indeed rarely found in the nature and they are usually contaminated with other elements or compounds. In particular case of planetary surfaces, mineral composition, mineral mixing, grain size and moisture content are important factors influencing the spectral signature of the endmembers. This phenomena in the literature is referred to as endmember variability. One way to address the issue of endmember variability is through overcomplete dictionary learning techniques.

These approaches aim at estimating several spectral signatures for only one endmember. Through constructing dictionaries of overcomplete atoms (endmembers inside dictionary), with the number of atoms much larger than the true number of endmembers within the scene, these techniques try to increase the accuracy of the analysis (unmixing, classification etc.) by having a more robust representations of the actual endmembers.

Noise and Anomaly

In addition to the earlier mentioned challenges, many domains of hyperspectral image analysis are concerned by important defects of noises. The nature of such noise is variable and depends on the atmospheric conditions at the time of image acquisition, captors precision and/or physical malfunctioning of the captors at the onset of image registration (spectrally and spatially). In general the nature of noise inside hyperspectral imagery can be divided in two categories:

• Noises which affect specific bands of the image.

• Noise which affect certain pixels of the image, across the bands of the image cube. In this thesis the later type of noise is brought under consideration (henceforth being referred as anomaly). Figure 1.11 illustrates the presence of anomaly within a hyperspectral image from 4-Vesta asteroid. Anomaly detection received not much attention during the early years because of :

1. the very poor spatial and spectral resolution of earlier multi-spectral images which made anomalies either embedded or mixed with other materials within a pixel.

2. the very spatial character of the image processing techniques that designed to take advantage of spatial correlation to perform the analysis, brought considerable limitations. Due to the spatial scarcity of anomalies and even worse their existence in a sub-pixel scale, those spatial domain-based methods mostly failed in capturing anomalies.

By rapid improvements in spectral and spatial resolutions of the hyperspectral images, many subtle materials which were unknown and undetectable by visual assessment, became detectable. In many cases including astronomical studies and planetary surface observations, these subtleties are most interesting targets and contain crucial information. According to [START_REF] Chang | Hyperspectral data processing: algorithm design and analysis[END_REF] anomalies are characterized by three characteristic identifiers :

• Their existence which usually occur with low probability (w.r.t nominal pixels) and without any a priori knowledge.

• Their presence which in many cases can be under/over estimated as they can be presented as subpixel targets.

• Their population which by nature should not be large compared to the total number of pixels.

Unfortunately no sound number (probability, ratio, ...) is available in practice. Indeed the size of population to be considered as anomaly is not well defined in most scenarios. Moreover, the question of defining a suitable distance between anomaly and nominal data, as well as between anomalies of different nature is tricky. Lastly, how an unsupervised global algorithm can reliably account for these issues while performing other image analysis tasks like spectral unmixing? The design of unmixing approaches robust to anomalies is then a tricky task that is addressed in this PhD thesis.

Contributions of this PhD

One of the contributions of this thesis is a new method that scales favorably with the problem posed by planetary exploration. It performs an unsupervised non-linear unmixing jointly with anomaly detection capacities, and has a global linear complexity. Non linearities are handled by analyzing the data in a Reproducing Kernel Hilbert Space (RKHS) while endmember variability is managed by decomposing the hyperspectral data on an overcomplete set of spectra, combined with a specific sparse projection. A theoretical study is proposed on synthetic datasets, and results are presented over the challenging 4-Vesta asteroid dataset. In another attempt, a new supervised unmixing algorithm is presented and focused on the exploitation of Optimal Transport. Optimal Transport consists in estimating a plan that transports a spectrum onto another with minimal cost, enabling to compute an associated distance (Wasserstein distance) that can be used as an alternative metric to compare hyperspectral data. This is exploited for spectral unmixing where abundances in each pixel are estimated on the basis of their projections in a Wasserstein sense onto known endmembers. In this work an over-complete dictionary is also used to deal with internal variability between endmembers, while a regularization term, also based on Wasserstein distance, is used to promote prior proportion knowledge in the endmember groups. The last part of this PhD focuses on the application of the above mentioned algorithms in the framework of the CATUT (French acronym for Urban Temperatures Mapping) project. The aim of this project is to estimate urban land surface temperatures (LST) using mainly hyperspectral data. The long term objective is to specify the characheristics of a future hyperspectral satellite sensor able to estimate LST. The process is first to estimate the nature of the ground and then, to estimate the associated temperature using emissivity models associated to each material. Because of physical properties, it is impossible to provide such a satellite with the very high spatial resolution required to identify all material inside an urban scene. As a consequence, hyperspectral unmixing is of prime importance here. To this end, LST estimation is mentioned as one of the goals of the Dual-use European Security IR Experiment 2008 (DESIREX 2008) campaign in Madrid Spain which carried out by European space agency (ESA) through several missions between June 23rd and 6th of July. The data used in this report, are the result of this campaign. Figure 1.12 is the classification map of a part of the city Madrid based on one of the airborne images of these series. This PhD is structured as follows :

• Chapter two comes with a review of state-of-the-art. First a formal definition of the mixing process is given. Then by introducing the possible existing mixing models (linear and non-linear), different procedures for addressing the spectral unmixing are presented. In particular supervised and unsupervised methods of unmixing are explained. Further the issue of endmember estimation is explained together with the problem of endmember variability are presented.

We also discuss about existing hyperspectral measures and their relationships. This chapter will be finished by giving an overview on the state-of-the-art of hyperspectral anomaly detection procedures.

• Chapter three is devoted to an unsupervised kernel-based method of anomaly detection and endmember extraction for unmixing. This is the first main contribution of this PhD.

• Chapter four details a supervised model of unmixing based on optimal transport theory. This is the second main contribution of this PhD.

• Chapter five is dedicated to a detailed application of the introduced methods within the context of CATUT project. Within this chapter we first give an outline of the goals of the project together with the procedure of data acquisition. Then by referring to the characteristics and limitations of the data we explain the data preprocessing steps necessarily to conduct the unmixing procedures. At the end the results of experiments are presented together with some notes on their analysis. This chapter will complete this thesis.

• Chapter Six comes with the conclusion and also gives indicators to the possible future extensions of the current algorithms. As it said earlier due to a number of reasons the values of sensors at the time of acquisition could be mixtures of two or more pure materials. In such cases the classification of hyperspectral pixels (assigning one label to each pixel) turns less representative. This calls for a way to find and retrieve information which is hidden inside each pixel value regarding the types of materials being mixed inside a pixel and their ratio of dominance. Following the introduction provided in 1.1.3 on the issue of spectral unmixing, within this chapter we provide a formal definition of mixing processes and introduce the latest achievements in the field of spectral unmixing. In particular linear vs nonlinear and supervised vs unsupervised methods of unmixing will be under focus. Then we will give an overview of the existing spectral measures and distances used for spectral comparison. In the last section of this chapter different approaches for hyperspectral anomaly detection are being introduced. The aim of spectral data analysis is to provide comprehensive descriptions of hyperspectral measurements. In the scope of spectral unmixing it consists of detection of endmember signatures and quantification of their respective abundances over the scene of study. The unmixing of a dataset X = [x 1 , . . . , x N ] (where x i is of dimension (1 × D) and represents the spectral signature in a D dimensional space) is an inverse procedure which approximates X with the product of two (smaller in size) matrices X ≈ EG where E is of dimension (D × ) whose columns contain the endmembers and G is of dimension ( × N ) and represents the projections of each pixel onto this endmember basis. Here we call E the endmemers matrix and G the abundance matrix. The concept and definition of endmember is highly dependent on the application. The abundance coefficients usually come with non-negativity (g j ≥ 0, j = 1, . . . , ) and additivity ( j=1 g j = 1) principles, due to the general assumption that they quantify the presence of each endmember within the pixels. These two constraints however, are also relative and not necessarily hold for all different types of mixtures. In some cases for example fractional masses of the materials could be used for interpretation of abundance coefficients. We can write the approximation of the value of one hyperspectral pixel x i as:

x i ≈ ϕ(E, g i ), (2.1)
where the approximation sign stands for the measurement noise or errors inherent to the function ϕ(•, •). This way, the solution of the spectral unmixing (SU) will be the following minimization :

min E,G N i=1 D i [x i |ϕ(E, g i )], s.t g i,j ≥ 0, j = 1, . . . , , (2.2)
where D(•|•) measures the discrepancy between the observed measurement x i and the prediction function ϕ(E, g i ). In many cases the approximation in 2.1 assumes the existence of additive i.i.d noise centered and normally distributed, with covariance matrix Σ. By this assumption in mind we can take an optimization view and solve a weighted least squares problem.

D i (u|v) = ||u -v|| 2 Σ , ∀i (2.3)
where ||u|| 2 Σ = u Σ -1 u is the Mahalanobis norm which is induced by Σ. Regardless of the type of framework (geometrical, statistical etc.) used for the formulation of the SU problem, the most important issue roots in the choice of ϕ(•, •) which relates the observed pixel with the endmember basis matrix E. The choice of this function reveals crucial information regarding the underlying mixing processes. In the ideal case, where there is no effect of noise, also with the assumption of pure endmembers being present in the data, the abundance coefficients as the result of inversion part of linear unmixing will lay within a convex polygon (n-1 simplex) whose vertexes are consist of the endmembers. Illustration courtesy of Dobigeon et al.

(non)Linear mixtures & (Un)supervised methods

A large number of algorithms have been proposed to deal with HS unmixing (see [START_REF] Bioucas-Dias | An overview on hyperspectral unmixing: Geometrical, statistical, and sparse regression based approaches[END_REF][START_REF] Bioucas-Dias | Hyperspectral unmixing overview: Geometrical, statistical, and sparse regression-based approaches[END_REF] for a comprehensive survey). These algorithms, aim either at estimating endmembers and the associated mixing abundances (unsupervised methods) or they focus only on estimating abundance coefficients on the assumption that endmembers are known (supervised methods). As for unsupervised scenarios the estimation of E and G, depending on the underlying approach can be either simultaneous as in [START_REF] Chris Hq Ding | Convex and semi-nonnegative matrix factorizations[END_REF] or sequential as in [START_REF] Courty | Saga: sparse and geometry-aware non-negative matrix factorization through non-linear local embedding[END_REF]. The unmixing methods can also be classified based on their a priori assumption on Linearity or Non-linearity of the mixing procedure.

In linear mixing model (LMM), the spectral signature representing each pixel is assumed to be a linear combination of endmembers. It gives also a good eyesight to look SU as a soft multi-class classification, as each spectral signature is linearly associated to the endmembers. In Figure 2.1 one can see the distribution of the abundance coefficients as the result of linear mixing assumption. when the additivity constrain does not hold, the coefficients instead belong to a convex cone defined by the endmembers. One of the earliest unmixing procedures is independent component analysis (ICA) [START_REF] Hyvärinen | Independent component analysis[END_REF] which comes with linearity assumption. In [START_REF] Lennon | Spectral unmixing of hyperspectral images with the independent component analysis and wavelet packets[END_REF], an enhanced version of ICA based on analysis of the spectra issued from wavelet packets transformation is introduced.

Though unmixing algorithms with linear assumption, show relatively good balance between computational complexity and acceptable accuracy, the underlying linear hypothesis is violated in many practical cases [START_REF] Bioucas-Dias | An overview on hyperspectral unmixing: Geometrical, statistical, and sparse regression based approaches[END_REF]. for example because of multiple scattering phenomenon [START_REF] Jie Chen | Nonlinear unmixing of hyperspectral data based on a linear-mixture/nonlinear-fluctuation model[END_REF] that occurs due to the interaction of light with several materials (in macroscopic scale). It can also be due to intimate mixing of materials which refers to the simultaneous interaction of light with all the endmembers (in microscopic scale) [START_REF] Dobigeon | Nonlinear unmixing of hyperspectral images: Models and algorithms[END_REF]. Non-linear techniques have thus been developed to cope with this difficulty including intimate mixtures, bilinear models and other physicsbased nonlinear mixing models in order to give an approximate formulation of the real mixing process. In Figure 2.2 three different mixing scenarios are illustrated. For intimate mixtures, a number of models based on the modeling and inversion of radiative transfer [START_REF] Hapke | Theory of reflectance and emittance spectroscopy[END_REF] have been derived. The main advantage of such approaches is the interpretability of the mixing process, since they involve real physical quantities such as grain size, roughness, single-scattering albedo, etc. However, such physical quantities, most of the time depend on unknown parameters, related to photometry and optical constants which yield a difficult estimation process. To relax the dependency on physical parameters, analytical formulations based on bilinear models have been successfully applied to represent the multiple scattering phenomenon which is the interactions of light with several materials. According to BMM, the value of pixel x i is being approximated w.r.t the basis matrix E, this time using two parameters α as in LMM for calculation of the linear term and β for calculation of the non-linear term as:

x i ≈ j=1 α (j,i) e j + -1 k=1 h=k+1 β (k,h,i) e k e h , (2.4)
where stands for the Hadamard product and e j is the jth endmember in E. Efficient associated techniques can be found in [START_REF] José | Nonlinear mixture model for hyperspectral unmixing[END_REF], the general bilinear model [START_REF] Halimi | Nonlinear unmixing of hyperspectral images using a generalized bilinear model[END_REF], linear-quadratic mixing [START_REF] Meganem | Linear-quadratic mixing model for reflectances in urban environments[END_REF] where recently in [START_REF] Zohra Benhalouche | Hyperspectral and multispectral data fusion based on linearquadratic nonnegative matrix factorization[END_REF] is used for the fusion of highspectral and high-spatial resolution HSI. Also a comparison of these models in a forestry context can be seen in [START_REF] Fan | Comparative study between a new nonlinear model and common linear model for analysing laboratory simulated-forest hyperspectral data[END_REF]. To deal with nonlinearity, kernel-based algorithms which project the input data in a feature space of (possibly very) high dimension where linearity is improved also offer interesting alternatives [START_REF] Jie Chen | Nonlinear unmixing of hyperspectral data based on a linear-mixture/nonlinear-fluctuation model[END_REF][START_REF] Camps | Kernel-based methods for hyperspectral image classification[END_REF].

It worths mentioning that although many unmixing algorithms assume one type of non-linearity in the mixing process, some of them (see [START_REF] Plaza | Nonlinear mixture models for analyzing laboratory simulated-forest hyperspectral data[END_REF] as an example) do not take any assumption on the type of nonlinearity in this process (as those working on a feature space through the kernel-trick, as detailed in the following paragraph). The unsupervised kernel-based method, introduced in this thesis belongs to this family of methods. In the following we discuss about endmember and abundance estimation.

Estimating endmembers and abundances

To estimate the pure elements, also called endmemders, a number of methods assume that at least one pure instance per endmember is included in the data [START_REF] Chang | A fast iterative algorithm for implementation of pixel purity index[END_REF][START_REF] Michael | N-findr: An algorithm for fast autonomous spectral endmember determination in hyperspectral data[END_REF][START_REF] Nascimento | Vertex component analysis: a fast algorithm to unmix hyperspectral data[END_REF]. Various techniques exist to extract them. They are either based on a pixel purity index [START_REF] Chang | A fast iterative algorithm for implementation of pixel purity index[END_REF] or on the identification of the simplex that encompasses data [START_REF] Craig | Minimum-volume transforms for remotely sensed data[END_REF][START_REF] Hollósi | Convex constraint analysis: a natural deconvolution of circular dichroism curves of proteins[END_REF]. This latter idea exploits the property that a linear combination of some points is included in the simplex, issued from these points; as a consequence extreme points of the simplex correspond to endmembers. N-Finder [START_REF] Michael | N-findr: An algorithm for fast autonomous spectral endmember determination in hyperspectral data[END_REF] is for example, issued from this representation since the maximum volume of a simplex is defined by the purest pixels in a dataset. According to the iterative version of this method, it is assumed that the number of endmembers is known a priori. By applying a dimensionality reduction technique like PCA [START_REF] Wold | Principal component analysis[END_REF] or MNF [START_REF] Andrew A Green | A transformation for ordering multispectral data in terms of image quality with implications for noise removal[END_REF], the dimensionality of the original space of data is being reduced to -1. The iterative part of the algorithm starts by randomly selecting signatures from the data, taking them as endmembers, then within an outer loop of k iteration the algorithm tries to find the volume of the simplex specified by those endmembers. In particular for the -simplex S{e 1 , . . . , e }, the volume of S is being calculated as:

V(e 1 , . . . , e ) = det 1 1 . . . 1 e 1 e 2 . . . e ( -1)! . (2.5)
Within an inner loop all the data samples are being added to S and the volume of the new simplex will be re-calculated in order to reach the maximum volume of the simplex.

Techniques based on non-negative matrix factorization (NMF) [START_REF] Jia | Constrained nonnegative matrix factorization for hyperspectral unmixing[END_REF][START_REF] Miao | Endmember extraction from highly mixed data using minimum volume constrained nonnegative matrix factorization[END_REF][START_REF] Sajda | Recovery of constituent spectra using nonnegative matrix factorization[END_REF][START_REF] Yokoya | Generalized bilinear model based nonlinear unmixing using semi-nonnegative matrix factorization[END_REF] exploit also this property and the simplex is sometimes computed through geodesic distances [START_REF] Heylen | Calculation of geodesic distances in nonlinear mixing models: Application to the generalized bilinear model[END_REF] or in Reproducing Kernel Hilbert Space (RKHS) [START_REF] Courty | Saga: sparse and geometry-aware non-negative matrix factorization through non-linear local embedding[END_REF]. This kernel-trick procedure enables to project data in a feature space (generally of high dimension), where the linearity assumption is more relevant. This efficient strategy has enabled to tackle a number of unsupervised, non-linear unmixing processes (see [START_REF] Dobigeon | Nonlinear unmixing of hyperspectral images: Models and algorithms[END_REF][START_REF] Jie Chen | Nonlinear unmixing of hyperspectral data based on a linear-mixture/nonlinear-fluctuation model[END_REF][START_REF] Broadwater | Kernel methods for unmixing hyperspectral imagery[END_REF][START_REF] Broadwater | A comparison of kernel functions for intimate mixture models[END_REF][START_REF] Broadwater | Kernel fully constrained least squares abundance estimates[END_REF]). Finally let us point out VCA (Vertex Component Analysis [START_REF] Nascimento | Vertex component analysis: a fast algorithm to unmix hyperspectral data[END_REF]) which is based on the idea that affine transformation of a simplex is also a simplex. According to this method, endmembers are being estimated iteratively by projecting data onto the direction which is orthogonal to the subspace spanned by the previously found endmembers, until the preset number of endmembers is reached.

Similar to N-finder, VCA knows the number of endmembers and starts by dimensionality reduction in the preprocessing step but the initial endmembers are being calculated as A = [e 0 , 0, . . . , 0], where e 0 = (0, 0, . . . , 1) . Then in each iteration the volume of a convex hull is being calculated after adding new orthogonal vector to the current basis :

f k = I -A (k-1) A (k-1) # w k / I -A (k-1) A (k-1) # w k . (2.6)
Although efficient, all these techniques assume the presence of pure elements in the image which is not guaranteed in practice. By relaxing the assumption of existing pure pixels, the minimum volume(MV) algorithms try to estimate the closest simplex encompassing data, making the estimation a much more difficult non-convex optimization procedure. MVC-NMF [START_REF] Miao | Endmember extraction from highly mixed data using minimum volume constrained nonnegative matrix factorization[END_REF], SISAL [START_REF] José | A variable splitting augmented lagrangian approach to linear spectral unmixing[END_REF], MVSA [START_REF] Li | Minimum volume simplex analysis: A fast algorithm to unmix hyperspectral data[END_REF] etc. are examples of this trend of SU algorithms. In the next part we consider a new approach for addressing the issue of non-pure pixels.

Endmember variability & Overcomplete dictionaries

The problem of endmember quality has been the focus of several recent studies [START_REF] Zare | Endmember variability in hyperspectral analysis: Addressing spectral variability during spectral unmixing[END_REF][START_REF] Revel | Inertia-constrained pixel-by-pixel nonnegative matrix factorisation: a hyperspectral unmixing method dealing with intra-class variability[END_REF]. As a matter of fact in many cases endmembers are either estimated up to some uncertainties or characterized not only by a single spectrum but rather by a set of different spectra. Using only one spectrum to decompose the information is then hazardous and likely to decrease the efficiency of the unmixing process. In order to cope with these issues, overcomplete dictionaries of hyperspectral signatures have revealed their efficiency [START_REF] Bieniarz | Joint sparsity model for multilook hyperspectral image unmixing[END_REF][START_REF] Bieniarz | On the use of overcomplete dictionaries for spectral unmixing[END_REF]. The general idea is to construct a set/dictionary of spectral signatures larger than the effective number of endmembers to deal with the internal heterogeneity inside a particular material (which should correspond to only one endmember in classical approaches). In some literature, methods where overcomplete dictionaries are used in SU are referred as semi-supervised methods. In order to handle the large number of signatures in dictionaries, the minimization procedure which is used in inversion phase is modified to account for sparsity constraint or penalization while the additivity constraint is being relaxed. Therefore, with a quadratic cost function, the SU consist in solving a sparse LS regression :

min g i ||x i -Mg i || 2 2 s.t ||g i || 0 ≤ λ , g j,i ≥ 0 ∀j, (2.7)
where λ is the sparsity level. The ||g i || 0 therefore, controls the permitted number of endmembers from the dictionary M to be involved in the mixture. Greedy algorithms are needed for solving these NP-hard problems. Well-know techniques are including but not limited to Orthogonal matching pursuit [START_REF] Chandra Pati | Orthogonal matching pursuit: Recursive function approximation with applications to wavelet decomposition[END_REF] and SUNSAL [START_REF] José | Alternating direction algorithms for constrained sparse regression: Application to hyperspectral unmixing[END_REF].

The overcomplete dictionary principle will be explored in the inversion phase of both supervised and unsupervised methods of unmixing in this thesis. Before embarking upon the inversion part we first introduce some distance measures that are commonly used in unsupervised SU procedures.

Distance Measures for HSI

Hyperspectral vector-based measures are mainly used in applications like detection, discrimination, classification and identification. In the case of unsupervised scenarios with vectorial data in particular they are the only means of evaluation and comparison. In this section therefore we introduce the major mathematical measures commonly used in analysis of HSI. In the following we introduce a number of vector-based measures which share some inherent similarities. Although these measures were not used in this thesis, because they show very interesting and promising properties, also the fact that their utilization is less emphasized in the remote sensing literature, we decided to give the definition of the sample spectral correlation based measures and their major subcategories.

Signature vector-based measures

Vector-based spectral measures are mainly used in applications like detection, discrimination and less frequently in classification. In the absence of prior class information, these measures are mainly used to compare two spectral signatures. These measures are effective mostly when the spectral vectors to be compared are believed to be the true representatives of the materials/endmembers which are assigned to. For two spectral signatures x i and x j of size (1 × D) we define:

Euclidean Distance (ED) ED(x i , x j ) = ||x i -x j || = D l=1 (x il -x jl ) 2 , (2.8)
which is the most well known distance measure defined over hermitian spaces.

Spectral Angle Mapper (SAM)

SAM(x i , x j ) = cos -1 x i • x j || x i |||| x j || , (2.9)
where || x i || is the norm of x i . SAM measures the angle between two vector variables. It is also referred as the cosine distance. When the angle between two signatures is small SAM and ED give similar results.

Orthogonal Projection Divergence (OPD)

OPD(x i , x j ) = x i P ⊥ x j x i + x j P ⊥ x i x j 1/2 , (2.10)
where

P ⊥ x k = I -x k (x k x k ) -1
x k and I (D×D) is an identity matrix. This measure works based on the orthogonal subspace projection algorithm [START_REF] Joseph | Hyperspectral image classification and dimensionality reduction: An orthogonal subspace projection approach[END_REF], and finds the residuals of orthogonal projection of x i and x j onto each other. From this perspective OPD and SAM also find similarities, as both of them are projective methods.

Spectral Information Divergence (SID)

SID(x i , x j ) = D(x i ||x j ) + D(x j ||x i ), (2.11)
where D(x i ||x j ) is the relative entropy (Kullback-Leibler information divergence). SID measures discrepancy between two signature vectors x i and x j in terms of their corresponding probability mass functions. while OPD calculates the divergence of one signature projected onto the other, SID can be though of as the stochastic version of OPD where orthogonal projection is replaced by information divergence.

Mixed SID -SAM Also By fusion of SAM and SID in [START_REF] Du | New hyperspectral discrimination measure for spectral characterization[END_REF] it is shown that:

SID -SAM = SID(x i , x j ) × tan(SAM(x i , x j )) (2.12)
and

SID -SAM = SID(x i , x j ) × sin(SAM(x i , x j )) (2.13)
both show better discriminability in some cases [START_REF] Du | New hyperspectral discrimination measure for spectral characterization[END_REF].

Correlation-weighted measures

The idealistic assumption based on which the vector-based measures operate (presented in 2.2.1) is not true in many cases. In particular in mixed signature discrimination and identification where spectral vectors are usually mixed with a number of other signatures. Another situation is the subsample target discrimination and identification where the interesting target/endmember is embedded in a single signature and mixed with other signatures. In such situations, the precision of vector-based measures deteriorates. In order to ameliorate this problem correlation-weighted measures extend the vector-based measures by including the sample spectral correlation to account for spectral variability of the materials. Based on the way these models incorporate the correlation information, the correlation-weighted measures are being categorized into measures weighted by a priori correlation and measures weighted by a posteriori correlation.

HS measures weighted by a priori correlation

By incorporation of a priori correlation, it is meant that there is a known correlation available prior to the comparison of signatures (identification or discrimination). The orthogonal subspace projection (OSP) approach introduced in [START_REF] Joseph | Hyperspectral image classification and dimensionality reduction: An orthogonal subspace projection approach[END_REF] is exploiting this a priori knowledge. It is possible to decorrelate signature vector x i , by orthogonally projecting it onto the basis formed by the columns of E which means the space < E > ⊥ orthogonal to the space linearly spanned by the columns of E. Different identification and discrimination measures are forged based on the OSP principle

P ⊥ E = I -E(E E) -1 E , (2.14)
where E is the set of endmember signatures.

HS measures weighted by a posteriori correlation

As availability of the a priori knowledge on the existing endmembers E is not always guarantied, an approximation of this knowledge is shown [START_REF] Chang | Hyperspectral data processing: algorithm design and analysis[END_REF] to be achievable by the inverse of the sample spectral correlation matrix R -1 . In this respect sample spectral correlation/covariance can be used to come up with a posteriori correlation weighted measures.

Inversion

Once endmembers known, an optimization procedure is often used to estimate the abundances in each pixel. It is called inversion and usually takes the form ĝi = arg min g ||x i -ϕ(E, g)|| 2 2 where E represents the endmembers and ϕ(E, •) is a (non)linear function that links the vector of abundances g to the current pixel x i to unmix. As we mentioned in 2.1.1 various mixing processes ϕ and penalization norms have already been proposed (see for example [START_REF] Halimi | Nonlinear unmixing of hyperspectral images using a generalized bilinear model[END_REF][START_REF] José | Nonlinear mixture model for hyperspectral unmixing[END_REF][START_REF] Jie Chen | Nonlinear unmixing of hyperspectral data based on a linear-mixture/nonlinear-fluctuation model[END_REF]). In cases where the discrepancy measure takes the form of euclidean distance, this optimization procedure boils down to constrained least squares (CLS) regression. The difficulty with the imposed constraints (resulting in impossibility of a closed-form solution) gave raise to other forms of optimization including fully constrained least squares (FCLS) [START_REF] Heinz | Fully constrained least-squares based linear unmixing [hyperspectral image classification[END_REF] and later on their kernelized versions [START_REF] Broadwater | Kernel fully constrained least squares abundance estimates[END_REF]. In the case of FCLS we formulate the minimization problem as :

ĝi = min g (x i -Êg) (x i -Êg) s.t g ≥ 0 and j=1 g j = 1 (2.15)
Because of the existence of the non-negativity inequality constraints there are no closed-form solution for the optimization of FCLS. In [START_REF] Heinz | Fully constrained least-squares based linear unmixing [hyperspectral image classification[END_REF] an iterative solution is given which simultaneously enforces the non-negativity and sum-to-one constraints on the abundance values.

For this, in order to account for the sum-to-one constraint a surrogate of the endmember matrix E and the datum x i is calculated as:

Ê = δE 1 (2.16) and xi = δx i 1 (2.17)
Where 1 = (1, . . . , 1) and δ controls the impact of the sum-to-one constraint.

Since g ≥ 0 is a set of inequalities, it is impossible for Lagrange multiplier to reach optimal solutions, therefore a positive constraint constant vector c is being defined to take care of non-negativity constraint. FCLS can be derived by forming a Lagrangian J as:

J = 1 2 (x i -Êg) (x i -Êg) + λ(g -c), (2.18)
where c = (c 1 , . . . , c ) is a self defined positive constraint constant vector and c i > 0, ∀i. Then iterating trough the two equations:

ĝ = ( Ê Ê) -1 Ê xi -( Ê Ê) -1 λ (2.19) and λ = Ê (x i -Êĝ), (2.20)
which can be used to find the optimal solution ĝ and the Lagrange multiplyer vector λ = (λ 1 , . . . , λ ).

It should be noticed that in many real-world scenarios the inverse problem is ill-posed and this bounds the efficiency of matrix factorization models like FCLS where the strong assumption of linearity is hold for the mixing process. This is why other discrepancy measures, including KL-divergence are being used instead of the ED. Also care must be taken that least squares does not produce sparse coordinates, while in most cases it is highly probable that the pixel is a mixture of few number of endmembers λ where λ . Therefore methods which account for this sparsity have been gained attention during the past years. In [START_REF] Lu | Manifold regularized sparse nmf for hyperspectral unmixing[END_REF], taking the so-called manifold assumption, a graph regularizer is added to sparsity constrained NMF in order to account for the sparsity of solutions as well as the geometry of the manifold embedding data. For this a nearest neighbor graph is created with data points as its vertices together with a weighting matrix based on the heat kernel. The optimization will be an extension of L 1/2 N M F introduced in [START_REF] Qian | Hyperspectral unmixing via l_{1/2} sparsity-constrained nonnegative matrix factorization[END_REF]:

min (E,G) = 1 2 ||X -EG|| 2 F + λ||G|| 1/2 + µ 2 T r(GΛG ), (2.21)
where µ ≥ 0 is a regularization parameter and Λ is a symmetric matrix calculated based on the weighting matrix. Similar to FCLS an iterative procedure with the use of Lagrange multipliers solves the minimization. As for the initialization of E and G, it is proposed to use VCA (for initialization of E) and FCLS (for initialization of G). Similar to FCLS, during the updates, 2.16 and 2.17 will be calculated based on E and x i , in order to guaranty the non-negativity and sum-to-one constraints.

In the next section we give an overview on the current methods of anomaly detection for HSI.

Anomaly detection

As mentioned in chapter 1, anomalies are likely to be present in hyperspectral data. Although noise can more or less be efficiently managed with some dimensionality reduction and denoising techniques, anomalies, because of their coherent structure, still affect their precision. All techniques, relying on endmember estimation directly from the image are sensitive to anomalies. For example, strategies based on extraction of a simplex will embed anomalies that are most of the time outside the simplex [START_REF] Li | Minimum volume simplex analysis: A fast algorithm to unmix hyperspectral data[END_REF][START_REF] Chan | A simplex volume maximization framework for hyperspectral endmember extraction[END_REF][START_REF] Nascimento | Vertex component analysis: a fast algorithm to unmix hyperspectral data[END_REF]. The identification of these existing anomalies is then of prime importance.

Many anomaly detection techniques rely on statistical approaches where the goal is to identify among data, points with strongly different statistics w.r.t the majority of data [START_REF] Matteoli | A tutorial overview of anomaly detection in hyperspectral images[END_REF][START_REF] Stein | Anomaly detection from hyperspectral imagery[END_REF], these latter being viewed in a local context (local neighborhood of the current test pixel) or global one (the whole scene). On this basis, a number of approaches have been proposed as for example the use of a Gaussian Mixture Model (GMM) [START_REF] Veracini | Fully unsupervised learning of gaussian mixtures for anomaly detection in hyperspectral imagery[END_REF] or Gaussian Markov Random Field [START_REF] Schweizer | Hyperspectral imagery: clutter adaptation in anomaly detection[END_REF] to model the nominal distribution of data.

Let us finally outline that other approaches try to represent data in other basis to better identify outliers, as in [START_REF] Mohsen | Improvement of anomoly detection algorithms in hyperspectral images using discrete wavelet transform[END_REF] where a Discrete Wavelet Transform (DWT) is used to represent points whereas a kernel-PCA is exploited in [START_REF] Chapel | Anomaly detection with score functions based on the reconstruction error of the kernel pca[END_REF] and [START_REF] Hoffmann | Kernel-PCA for novelty detection[END_REF].

By this overview on different concepts and tools which we encounter in the realm of remote sensing and hyperspectral data analysis we finish this chapter and in the next one an original algorithm of spectral unmixing and anomaly detection which exploits many of the above mentioned concepts and tools is being presented.

In this chapter we first introduced the concept of spectral mixing and the known mathematical models used for explaining different types of mixtures. In particular we made the difference between linear and non-linear mixing models clear. We also provided an introduction to the existing and frequently used methods of spectral unmixing, their assumptions and different properties that make each suitable for addressing certain types of mixtures. Then we introduced different metrics that are frequently used for measuring difference between spectral signatures. In the next part we explained the inversion procedures for calculating the abundance coefficients corresponding to each endmember. At the end we introduced the common problem of anomaly and the existing algorithms for dealing with anomalies within hyperspectral data.

Part II Kernel methods for Unmixing

Chapter Summary

Endmember extraction and anomaly detection on synthesized 2D data

In this chapter a nonlinear unsupervised kernel-based model of unmixing associated with an anomaly detection procedure is presented. This model is meant to address the ubiquitous problem of noisy datasets where the existence of noise distorts the correct representation of the manifold on which the data exists. The heart of this representation lays on using a kernel to project data from its original space (where data sits on a non-simplical manifold) onto a feature space where the projected data is embedded by a simplical manifold. The unmixing of a dataset X = [x 1 , . . . , x n ] (where x i is of dimension (1 × n) and represents the spectral signature in a n dimensional space) is an inverse procedure which approximates X with the product of two (smaller in size) matrices X ≈ EG where E is of dimension (D × ) whose columns contain the endmembers and G is of dimension ( × N ) and represents the projections of each pixel onto this endmember basis. Here we call E the endmemers matrix and G the abundance matrix. In 3.1 we first give a short review of the applications of HSI in planetary science, as well as existing challenges for processing planetary images. We provide the definitions and short overview of the latest approaches for spectral unmixing and mention the importance of anomaly detection for planetary hyperspectral imagery. In 3.2 we give a thorough explanation of the model. First in 3.2.1 we will provide the formal definition for exact simplex volume maximization procedure and then in 3.2.2 we derive the iterative procedure of anomaly detection. In 3.2.3 we mention how the abundance values are being calculated through sparse projections onto the simplex. In 3.2.4 the two major parts of the model are brought in psudocode form in order to clarify the whole procedure. Then in 3.3 we set up a series of experiments on synthetic and real datasets. First in 3.3.1 we explain the procedure of generating synthetic datasets and the parameters which are involve in it. Then the results of applying the algorithm on these datasets are demonstrated in 3.3.2 together with their comparison with other unmixing techniques. In 3.3.3 we demonstrate the benefits of exploiting overcomplete dictionaries in estimation of the endmembers and finally in 3.3.4 we demonstrate the performance of our joint anomaly detection and spectral unmixing algorithms on a real dataset. This chapter corresponds to an article published in IEEE TGRS under the title "Joint anomaly detection and spectral unmixing for planetary hyperspectral images" ( Volume: 54, Issue: 12, Dec. 2016 ).

Introduction

Hyperspectral (HS) remotely sensed images are of prime interest in many scientific fields, since they enable us to assess a dense spectrum (generally composed of several hundreds of contiguous electromagnetic wavelengths) in each pixel, making possible the identification of various materials composing the scene [START_REF] Stein | Anomaly detection from hyperspectral imagery[END_REF]. For this reason, hyperspectral data are vastly being used in earth observation (identifying land cover, crops, ... see for example [START_REF] Bioucas-Dias | Hyperspectral remote sensing data analysis and future challenges[END_REF][START_REF] Theodore | Hyperspectral remote sensing[END_REF]) as well as planetary exploration (identification of minerals and rocks [START_REF] Ceamanos | Intercomparison and Validation of Techniques for Spectral Unmixing of Hyperspectral Images: A Planetary Case Study[END_REF]). This latter application is the one we focus in this chapter. To understand the processes that drive the formation and evolution of planets, it is crucial to characterize their surface's composition. Visible-near infrared reflectance spectroscopy has long been recognized as an extremely powerful tool to achieve this objective as it allows to determine rocks composition at both local and global scales [START_REF] John | Interpretation of visible and near-infrared diffuse reflectance spectra of pyroxenes and other rock forming minerals[END_REF]. Imaging spectroscopy therefore, has been rapidly grown in planetary exploration, and very soon dedicated instruments were sent toward other planets in our Solar System [START_REF] Bibring | Results from the ISM experiment[END_REF]. Number of pioneering instruments have followed resulting in huge amount of data available for the scientific community. From the computational point of view however, the processing of hyperspectral images is particularly challenging when considering planetary exploration, since unlike the earth, no ground truth nor field campaigns exist except few particular cases (i.e. Apollo landing sites for the Moon and rovers/landers locations for Mars). Moreover, despite continuous spectral and spatial resolution improvements of HS images, they still suffer from a considerable low spatial resolution to observe every single material. As a consequence a pixel generally contains mixed information of several elements which are of key interest when considering planetary processes, especially magmatic rocks (i.e. rocks made of the minerals olivine and pyroxenes) that are used to evidence igneous processes typical of the crust or mantle. It should also be outlined that in addition to observation conditions (e.g. illumination, shadows) and to the nature of the surface itself (e.g. roughness, composition) which create non-linearities (as multiple scattering [START_REF] Jie Chen | Nonlinear unmixing of hyperspectral data based on a linear-mixture/nonlinear-fluctuation model[END_REF] or intimate mixing [START_REF] Dobigeon | Nonlinear unmixing of hyperspectral images: Models and algorithms[END_REF]), olivine and pyroxenes can have a wide range of chemical compositions, leading for each of them to a non-linear spectral response (e.g. [START_REF] Clenet | A new systematic approach using the Modified Gaussian Model: Insight for the characterization of chemical composition of olivines, pyroxenes and olivine-pyroxene mixtures[END_REF] and references therein). Combined to this effect, rocks are themselves a combination of minerals, implying to decompose the different spectra not on a basis of pure materials but rather as a combination of mixtures of materials.

Apart from these problems related to the specificity of hyperspectral data and planetary contexts, let us also note that anomalies often occur in the acquisition process of hyperspectral images. Anomaly is referred to any pixel whose signature is considerably different from its surrounding background. Several factors give rise to anomalies : 1) Natural degradation due to time spent in space (especially regarding missions far away from the Solar System) where unrealistic responses can appear progressively and randomly ; 2) Varying natures of sensors used to cover the whole wavelength range which generate misregristrations ; 3) Correction of atmospheric effects where for some planetary objects (e.g. Mars, Titan) residues can affect spectra in specific wavelengths.

Although efficient pre-launch calibration procedures exist, they are not sufficient to perfectly clean data. All these reasons give rise to the necessity of efficient nonlinear unmixing procedures, robust to outliers, able to decompose mixed pixels into a number of pure reflectance spectra called endmembers. This problem refers to an inverse procedure to determine endmembers and estimate fractional abundances for each pixel [START_REF] Boccacci | Introduction to Inverse Problems in Imaging[END_REF]. Traditionally, anomaly detection and spectral unmixing have taken as separate subjects. A number of families of unmixing techniques which aim at recovering abundances of endmembers directly from the image (or from a separate dictionary of acquired signatures) have been presented. Despite efficient and sound existing methods (see [START_REF] Bioucas-Dias | An overview on hyperspectral unmixing: Geometrical, statistical, and sparse regression based approaches[END_REF][START_REF] Bioucas-Dias | Hyperspectral unmixing overview: Geometrical, statistical, and sparse regression-based approaches[END_REF][START_REF] Dobigeon | Nonlinear unmixing of hyperspectral images: Models and algorithms[END_REF] for an overview), most techniques still suffer from degraded performances in presence of anomalies. In such situations, endmembers are often being incorrectly estimated and consequently, associated fractional abundances are less meaningful. Therefore, together with unmixing, an accurate detection of anomalies is crucial. This constitutes the aim of this chapter.

Summary, goals of the chapter

In this introduction we have discussed the necessity of both unmixing HS images and the detection of anomalies. Also interesting techniques have cited for solving spectral unmixing. However, to the best of our knowledge, the joint estimation of anomalies and unmixing hyperspectral data has not been proposed yet 1 . In this work we introduce a technique entitled SAGA+ (Sparse And Geometry Aware +) able to jointly detect anomalies and perform spectral unmixing. This Anomaly-Detection and Spectral Unmixing procedure is unsupervised and defined in RKHS. The idea consists of increasing the volume of the manifold hull in feature space. SAGA+, in an iterative way, increases the volume of the simplex by finding the purest pixels within the dataset, while at the same time keeping track of the rate of change of the sum of projection errors (SOPE) in order to exclude anomalies from being taken as simplex vertexes.

The chapter is organized as follows: Section 3.2 introduces the Anomaly Detection and Spectral Unmixing model mentioned above and in section 3.3, evaluation of the model and comparison results using both synthetic and real datasets are presented. Finally, section 3.4 concludes this work with an overall conclusion and pointers to possible future extensions.

SAGA+ : simplex volume maximization and anomaly detection

In this section the methodological contributions are presented. The spectral unmixing relies on a matrix factorization technique denoted SAGA (Sparse And Geometry Aware) and introduced in [START_REF] Courty | Saga: sparse and geometry-aware non-negative matrix factorization through non-linear local embedding[END_REF]. In the present contribution we add an anomaly detection and removal approach to deal with aforementioned issues, yielding an approach denoted "SAGA+". SAGA is based on the geometrical concept of finding the simplex that embeds data. This simplex is being computed in a feature space associated to a kernel. As for anomalies, they are tracked through their lack of representational capacity in the dataset.

We first present the endmember extraction technique in 3.2.1 and subsection 3.2.2 is devoted to anomaly detection. Once endmembers extracted, abundances are estimated by projecting data onto the endmembers. This procedure is detailed in 3.2.3. The overall process, entitled SAGA+ is summarized in 3.2.4. Before entering into details, let us introduce some general notations:

The hyperspectral image contains D bands and N pixels

x i ∈ R D , i = 1, ..., N . Data are represented in a D × N matrix composed of N vectors X = [x T 1 , . . . , x T N ].
The unmixing is formalized as finding a representation of the form : X ≈ EG where columns of E ∈ R D× contain the endmembers and G ∈ R ×N is the matrix, consisting of projections of each pixel onto this endmember basis. In practice the interpretability of matrix G is enhanced if each column is positive and sums to 1, since in this case, each value G ij contains the abundance of endmember i in pixel j. Such a Non-negative Matrix Factorization (NMF) is the basic idea on which we rely. To estimate the endmember matrix E, we rely on Column Subspace Selection principles: As for each column of E is selected an input pixel x i . This idea suggests that at least one sample of each endmember is present in the data (since this assumption is likely to be violated, in practice we rely on overcomplete strategies, as will be explained). The following section introduces the simplex approach to estimate matrix E.

SAGA principles: Exact Simplex Volume Maximization for endmember extraction

As mentioned earlier, we rely on a Column Subspace Selection procedure to construct the simplex that embeds the volume generated by the data [START_REF] Thurau | Descriptive matrix factorization for sustainability adopting the principle of opposites[END_REF]. To this end we start from the idea that the bigger this simplex is, the more likely it will contain the other elements, and therefore the reconstruction error of any projection onto the new simplex will be small or null. The vertices of the simplex then correspond to the endmembers. In this study, we also claim that choosing an overcomplete endmembers dictionary will help in performing non-linear unmixing. This point is illustrated in Figure 3.1. In the first case (Figure 3.1(a)), the spectra are living in a perfectly linear manifold. In this case, two endmembers (E1 and E2) are sufficient and any spectrum can be described as a convex combination of the two endmembers. This is referred to as the linear mixing hypothesis. However, this ideal case is usually different from reality, as illustrated on Figure 3.1(b), where the spectra live on a non-linear manifold. If the linear assumption is nonetheless assumed, then the spectra are projected on the 1-simplex formed by the two endmembers, resulting in an error in the corresponding abundances. By choosing an appropriate additional endmember (E3, blue circle in Figure 3.1(c)), one can build a piece-wise linear approximation of the non-linear manifold. The unknown spectra are now projected on the corresponding 2-simplex. Note here that this projection is naturally sparse, i.e. for all spectra only two endmembers are required to describe perfectly the projection. However this comes at the price of having one endmember that is in itself a mixture of the first two ones. There are then two cases:

1. Either one wants only to estimate the mixture coefficients over the two "pure" endmembers E1 and E2. Assume for the example that E3 is a mixture of known coefficients (0.4 and 0.6) of E1 and E2. Then for any spectrum with decomposition g 1 = 0.2, g 2 = 0 and g 3 = 0.8, the corresponding abundance w.r.t. only E1 and E2 is now g 1 = 0.2+0.8×0.4 = 0.52 and g 2 = 0+0.8×0.6 = 0.48. More generally, in the case where λ pure endmembers are available and endmembers are chosen for the overcomplete dictionary, we suppose that we have access to a (generally sparse) linear operator A of size (λ × ), with prescribed sums of rows = 1, that can be applied on any result of unmixing w.r.t. the overcomplete dictionary to get only a decomposition over pure endmembers.

2. However, we agree that obtaining this mixing matrix might not be a trivial task even for specialists. Nevertheless in some domains, as the one explored in this work, specialists are used either to work with spectra obtained in laboratory that are already combination of pure chemical elements. The decomposition on the overcomplete basis, provided it can be interpreted, gives in itself interesting insights on the underlying geological process.

In addition of the overcomplete strategy, to deal with complex manifolds issued from nonlinearities in X, it is possible to project the data in a feature space H using a projection function φ : R D -→ H (working in the original space is equivalent to choose the identity as projection function φ). We represent this projected data as: Φ = φ(X) = [φ(x 1 ), φ(x 2 ), . . . , φ(x N )]. Under the assumption that endmembers are invariant by transformation in the feature space, the endmember extraction therefore can be written as finding an indicator matrix W ( ) ∈ {0, 1} N × where:

φ(XW ( ) ) = ΦW ( ) . (3.1)
This gives rise to a simplicial convex hull in H, noted ∆ (φ(E)), whose pre-image is a non-simplicial hull (the manifold hull) in the original space. The procedure first extracts the best 1-simplex ∆ 1 by finding two instances from Φ that results in the minimum projection error on this basis, then iteratively increases the dimension of the simplex in H. At the (p)-th iteration the algorithm selects φ(x i ) such that the volume of the ∆ p-1 -simplex is maximized, i.e.

i = arg max q V ol(∆ p-1 (ΦW p-1 ) ∪ φ(x q )), (3.2)
where the volume of ∆ p simplex in this iteration is [START_REF] Courty | Saga: sparse and geometry-aware non-negative matrix factorization through non-linear local embedding[END_REF]:

V ol(∆ p-1 (ΦW p-1 ) ∪ φ(x q )) = V ol(∆ p-1 ) × dist(φ(x q )), ΦW p-1 ) p -1 . (3.3)
In practice the optimal projection φ enabling to separate complex data is unknown and potentially of high dimension. Hopefully, thanks to the kernel theory [START_REF] Camps-Valls | Kernel methods for remote sensing data analysis[END_REF], one can rely on kernel functions k(x i , x j ) = φ(x i ), φ(x j ) that allow to rewrite the relation dist(φ(x q ), ΦW p-1 ) in a closed form:

dist(φ(x q ), ΦW (p-1) ) = 1 -(k T xq • K -1 p-1 • k xq ), (3.4)
with K -1 p-1 is the inverse of the kernel matrix of the elements of the manifold:

K -1 p-1 = (W (p-1) Φ • ΦW (p-1) ) -1 (3.5)
and k x i is a vector of length p defined as

k x i = [k(x i , x j )] x j ∈XW p .
It is worth noting that instead of inverting the kernel matrix, in the implementation, for the sake of more numerical stability we chose to apply a Cholesky factorization on the kernel matrix at each iteration and then solve for it.

As for the choice of the kernel function, working on the original space is equivalent to choose k as the usual dot product : k(x i , x j ) = x i , x j . Choosing a different function enables to perform the process in a feature space where the geometry of the manifold embedding X can be more meaningful. Among existing kernels, the standard Gaussian Radial Basis Function (RBF) is widely used since it enables a projection in an infinite dimensional space :

k(x i , x j ) = exp - x i -x j 2 2σ 2 . ( 3.6)
The overall procedure enables to extract endmembers by finding the optimal manifold hull composed of spectra in the feature space. This exact simplex volume maximization procedure is the cornerstone of endmember estimation in SAGA [START_REF] Courty | Saga: sparse and geometry-aware non-negative matrix factorization through non-linear local embedding[END_REF]. In most methods assuming a linear mixing model, the number of endmembers is exactly defined by the number of pure spectra present in the image, and cannot exceed the spectral dimension. As such, choosing a point which is already a mixture of those pure spectra would lead to a nul increase in the volume of the simplex. In the RKHS corresponding to the Gaussian Kernel, each elements of Φ are orthogonal, and therefore would lead to a positive increase of the volume. In practice, it is therefore possible to have more endmembers than the number of pure materials (and can even exceed the spectral dimension). From this point of view, the process is a dictionary learning procedure, with an overcomplete dictionary (in the sense of the original spectral space). Apart from handling implicitly the non-linearity of the mixing model, this leads to two important properties of our method:

• The spectra can be decomposed over elements that are also combination of pure elements. This can ease the interpretation of the unmixing for specialists through the obtention of more physically realistic spectra, and alleviate to some extent the purity assumption which states that at least one pixel of each material should be present in the image;

• Choosing an overcomplete endmember set can help in handling the intra-class variability of each element, as similar versions of the same material (but with different spectral signatures, as it is the case with shadows for instance), can be efficiently exploited.

In practice, this also requires some sparsity assumption in the mixing coefficients, which is guaranteed by our projection method (see Section 3.2.3).

Anomaly Detection

As mentioned above, SAGA procedure selects endmembers in a way to increase the volume of the simplex hull in the feature space. According to this strategy, in each iteration, extremal points of the dataset hull are chosen. This criterion makes the algorithm vulnerable to anomalies since they are most likely to appear in the construction of our dictionary E. This issue is illustrated in the two first rows of Figure 3.2. In this figure, all instances of the first row are generated from a uniform Dirichlet distribution within the 2-simplex identified by blue vertexes (note that these pure points are not included in the dataset). In the second row 10 anomalies have been added to the nominal data. On the right column, the corresponding endmembers extracted by SAGA are depicted (red dots) using both datasets. From this simple experience, it is visible that endmember extraction algorithm fails to identify the correct borders of the nominal data in presence of anomalies.

As soon as extreme points appear, they significantly contribute in maximizing the simplex volume and are thus selected as endmembers. This behavior is observable in all unsupervised endmember extraction procedures which work based on maximization of the geometrical volume of the simplex encompassing data. To cope with this issue, the proposed anomaly detection (AD) procedure relies on the computation of the reconstruction error (in feature space) between a data sample and its perpendicular projection onto a side of the polygon of the simplex. Although similar ideas have already been proposed (see [START_REF] Chapel | Anomaly detection with score functions based on the reconstruction error of the kernel pca[END_REF] and other AD algorithms), the process here does no attempt for identifying the nominal distribution (for example based on decomposition of the eigen-structure of dataset) but we rather exploit the iterative structure of SAGA where in step p, after selection of the p + 1 -th endmember candidate through equation (3.2), we check its validity by computing the sum-of-projection-errors (SOPE) of all instances with respect to the new simplex ∆ p+1 :

SOPE(p + 1) = N i=1 1 -(k T x i • K -1 p+1 • k x i ) N , (3.7)
where N is the number of samples and acts as the normalization constant. A perfect reconstruction gives SOPE(p + 1) = 0, whereas its value increases when the quality of the reconstruction decreases (i.e. the simplex is not a good representation of the dataset). During the iterative procedure, when a considerable drop of SOPE is perceived after adding a new vertex to the simplex, this means that previously selected prototypes (endmember candidates of previous steps) could be anomaly and are among instances significantly different from nominal data. By this procedure one can then identify and reject them. In the third row of Figure . 3.2 one can see that this AD procedure, even when data are exposed to anomalies, endmember estimation remains close to a good representation of the data manifold which is in contrast to the second row of Figure .3.2 (without any AD process). In order to be able to check this drop in each iteration we need to define a threshold τ . This point is not trivial and highly depends both on the distribution of nominal data and the distribution of anomalies. To illustrate this, Figure . 3.3 represents the evaluation of SOPE when all nominal data (green points) are represented on a unique candidate issued from a meshgrid. From this figure, one easily observes that error decreases when endmember candidates get closer to the center of the data distribution2 . The iso-contour corresponding to τ = 0.01 in SOPE value is also represented to evaluate the associated threshold that enables to remove outlier data. A discussion related to the practical selection of τ is proposed in section 3.3.

Once endmembers are identified (together with anomalies), the projection of each data point onto the new basis has to be performed. This process is discussed in the following.

Sparse Projections onto the Simplex

Let us recall that the unmixing process aims at approximating the data matrix X through X ≈ EG, where the estimation of the endmember matrix E has been presented above. As for the abundance matrix G, it corresponds to the projection of X onto a region bounded by E. This is done, for each row of G (each row being represented as G •i associated to the ith point), with a projected gradient decent :

G •i = arg min G •i ||φ(x i ) -ΦWG •i || 2 s.t.G •i ∈ ∆ λ . (3.8)
where ∆ λ is the sub-simplex included in ∆ but composed only with λ vertexes. We indeed want to impose sparsity on each column of G since in practice, a pixel is likely to be composed only with a subset of existing endmembers. Explicitly, imposing sparsity prevents from solutions where all endmembers have a (even small) contribution. To this end, a greedy selector and sparse projector (GSSP) [START_REF] Anastasios | Sparse projections onto the simplex[END_REF] method is used not only to project the data onto the simplex but also to impose a sparsity level (through a number of non null elements given as an input parameter). In practice the coordinate values of each pixel will be sorted based on their magnitude and the λ greatest ones will be projected onto a unit simplex while the rest of the values will be set to zero.

We now have all ingredients for the complete process of anomaly detection and unmixing which is being presented.

Overall process

The SAGA+ procedure is presented in Algorithm 1. In a first step we compute the manifold hull to derive endmembers together with anomaly detection (lines 4 to 10). To initialize the process, the algorithm selects a random datum x i (line 4), and computes the distances between x i and all instances in X (line 5). These distances are then sorted in descending order and will be feed to the Exact Simplex Volume Maximization and Anomaly Detection process to estimate the first endmember êi associated with a set of anomalies (line 6). This process is repeated times in order to extract the endmembers (lines 7 to 9) and to construct the manifold hull E. In each iteration the volume increase of the simplex is computed in RKHS w.r.t each x i . Following equations (3.2)-(3.5), this requires to inverse the kernel matrix K -1 which is done in practice using the Cholesky decomposition. After computation of the volume changes, the anomaly detection algorithm is called to check the trade-off between the increase of the volume of the simplex and the decrease in the sum-ofprojection-errors (SOPE). This is represented in Algorithm 2.

Once the manifold hull extracted, the projection matrix G is calculated using the Greedy Selector and Sparse Projector (GSSP) algorithm. This second step of SAGA+ follows the process described in [START_REF] Anastasios | Sparse projections onto the simplex[END_REF]. The principle of Exact Simplex Volume Maximization and Anomaly Detection algorithm is described in Algorithm 2. Here we check iteratively, on the basis of the SOPE criteria presented in equation (3.7), if successive candidates for the simplex vertexes (sorted based on their distances with the current simplex) are real endmembers or lie among outliers. SOPE is computed by taking into account the presence of already found prototypes and in case no endmembers are already present, the sum of projection error reads SOP E = [(ν • ν)/N ] -1 where ν represents the vector of the similarities between the current candidate and the rest of the dataset instances. -Sort the vector of volume changes decreasingly (see eq. (3.2)).

-Since first indexes are possible anomalies, SAGA+ to identify endmembers and outliers (see algorithm in Figure 2) 8: end for 9: return(A,E) 10:

Step 2: Abundance estimation via sparse projections

Evaluation

In this section we evaluate on synthetic and real data, SAGA and SAGA+ procedures and compare their performances with a number of state-of-the-art unmixing procedures. In particular we consider two linear unmixing methods N-Finder [START_REF] Michael | N-findr: An algorithm for fast autonomous spectral endmember determination in hyperspectral data[END_REF] and VCA [START_REF] Nascimento | Vertex component analysis: a fast algorithm to unmix hyperspectral data[END_REF] introduced above. As mentioned, these are frequently used unsupervised methods of unmixing and their endmember extraction procedures are based on the simplex volume maximization method. This makes them suitable to be compared with the current method. It should be mentioned however that these models assume linear mixing and do not detect outliers, therefore the comparison with our approach would be thought to some extent unfair. Nevertheless it should also be emphasized that they all follow the same logic of extending the encompassing simplex which makes them (to some degree) comparable to the logic behind SAGA. We also decided to compare our approach with G-SiVM [START_REF] Heylen | Non-linear spectral unmixing by geodesic simplex volume maximization. Selected Topics in Signal Processing[END_REF] which is a non-linear procedure based on simplex volume maximization through shortest-path distances in a nearest-neighbor graph. Non-linearity is performed using kernel trick which makes this approach similar to SAGA and SAGA+ in the sense that they do not rely on any specific mixing assumption. In order to evaluate and compare the efficiency of the procedures, we use different criteria to assess endmember extraction and anomaly detection processes. These criteria are introduced bellow together with experimental setups.

Experimental setups

Evaluation criteria

for synthetic data where ground truth is available, a set of criteria are used to measure our performances:

• To evaluate the adequacy of estimated endmembers Ê = [ê 1 , . . . , ê ] with ground truth endmembers E = [e 1 , . . . , e ], several possibilities are offered as Spectral Information Divergence (SID) or RBF kernel distance, each of which give a measure of similarity between two signals represented as a vector. Though efficient, these criteria are sensible to the difference of endmembers in terms of magnitude while in this application, we are more interested in comparing their shape. Indeed, because of atmospheric attenuation or variation in illumination, spectral energy can be attenuated for two samples even if they belong to the same class. Thus, we rather prefer to rely on Spectral Angle Map (SAM) which is scale invariant. The SAM reads, for spectral vectors x i = [x i1 , ..., x iD ] T and x j = [x j1 , ..., x jD ] T with positive values:

SAM(x i , x j ) = cos -1 x i • x j || x i |||| x j || (3.9)
and lives in the range [0, 1] (1 representing maximum dissimilarity and 0 indicating two identical signals). In practice in order to compare two different sets of endmembers, we need a pairing procedure to associate each estimated endmember to a unique real one. This is done as finding the best permutation of estimated endmembers that gives the minimal average SAM values when they are compared to real ones. Therefore, two sets E and Ê have a difference (E, Ê) which reads:

(E, Ê) = min η 1 i=1 SAM(e η (i) , êi ) (3.10)
where is the number of endmembers and η : {1, ..., } -→ {1, ..., } is a function corresponding to all possible permutations.

• As for the evaluation of the anomaly detection algorithm, we rely on Cohen's Kappa Statistics [START_REF] Cohen | A coefficient of agreement for nominal scales[END_REF] to measure the adequacy of the classification between outliers/endmembers. Kappa criterion, noted κ, lies between 0 (nothing is in accordance) and 1 (perfect accordance) and combines False Positive (FP), False Negative (FN), True Positive (TP) and True Negative (TN) in order to measure the inter-rater agreement of categories.

• In complement we also give an execution-time comparison among the endmember extraction procedures.

We now introduce our datasets.

Synthetic data

In order to generate synthetic datasets, we have selected 41 spectral signatures from RELAB spectral library [START_REF] Clenet | A new systematic approach using the Modified Gaussian Model: Insight for the characterization of chemical composition of olivines, pyroxenes and olivine-pyroxene mixtures[END_REF] 3 . From spectral signatures E = [e 1 , . . . , e ] (e i ∈ R D ), mixed spectra have been generated based on three different mixing assumptions: Linear Mixing Model (LMM), Bilinear Mixing Model (BMM) and Highly Concentrating Model (HCM). These mixing models are well studied in the unmixing literature, though they do not fully represent the variability of mixing processes in real world hyperspectral images. The two first models (LMM, BMM) are implemented according to the procedure in [START_REF] Jie Chen | Nonlinear unmixing of hyperspectral data based on a linear-mixture/nonlinear-fluctuation model[END_REF]. The usage of HCM model is also mentioned in some earlier works [START_REF] Schmidt | Implementation strategies for hyperspectral unmixing using bayesian source separation[END_REF]. From the set of endmembers E, the Linear Mixture Model reads :

x i = E T g i + w i (3.11)
where g i is the abundance vector related to data x i and w i is a random noise. Each element in g i lives in the range [0, 1] and their sum equals 1. The fractional abundances in LMM have been generated based on Multivariate Symmetric Dirichlet Distribution, uniformly distributed over the standard ∆ -1 simplex. More precisely for each sample x i , i = 1, . . . , N and x i ∈ R D the abundance vector g i = [g 1 , . . . , g ] is generated through a Dirichlet Distribution. This latter, parametrized by vector α = [α 1 , . . . , α ] such that α i > 0, ∀i ∈ [1, ], reads:

g i = Γ(α • ) Γ(α) j=1 e α-1 j (3.12)
where Γ(•) is the Gamma function. Dirichlet are generalization of the beta distribution in a multi-variate context. Such distributions are then very useful to simulate mixture models. From equation (3.12), using a parametrization α = [α 1 , . . . , α ], we generate mixtures of pure endmemners [e 1 , ..., e ] by computing abundance parameters g i = [g 1 , . . . , g ] as the power-product between endmembers and α. This is written in the right part of equation (3.12). The left part is only a normalization coefficient that involves the Gamma function. In practice α is set as the vector of ones: α i = 1, ∀i. As the result the vector g i contains coefficients which are uniformly distributed across the vertexes of the simplex. To illustrate this, from pure endmembers visible on the top of For Bilinear Mixing Model, the power products of reflectance is added to the linear mixing model, leading to a model of the form :

x i = E T g i LMM + -1 k=1 l=k+1 β k,l,i e k e l + w i (3.13)
where β k,l,i are new mixing coefficients and the termwise Hadamard product is defined as:

e k e l =   e 1,k • • • e D,k     e 1,l • • • e D,l   =   e 1,k e 1,l • • • e D,k e D,l   . (3.14)
From the same pure endmembers on the top of Figure 3.4(a), blue points of Figure 3.4(c) represent the same 2D slice generated under the BMM model.

As for the Highly Concentrating Model, we choose to rely on the same Dirichlet Distribution as in LMM, however with large concentration parameters α i 1 (in practice we used α i = 50, ∀i). This ensures that generated abundance coefficients g i are highly concentrated in the middle of the simplex which is formed by the endmembers. This is visible in blue points of Figure 3.4(c).

In order to generate and include a number Υ of outliers to the nominal synthetic data, we used three spectral signatures E = [e 1 , e 2 , e 3 ] without any similarities to the signatures used to generate nominal data. From known endmembers in E and our abnormal spectral signatures in E , we use again Dirichlet distribution function to generate fractional abundances for anomalies g i using three more parameters α = [α 1 , α 2 , α 3 ]. The generation function reads:

g i = Γ( +3 j=1 γ j ) +3 j=1 Γ(γ j ) +3 j=1 e γ j -1 j (3.15) whereγ = [α 1 , ..., α , α 1 , α 2 , α 3 ]
is the concatenation of parameters α and α . In practice, we used α 1 , . . . , α = 1 and α 1 , α 2 , α 3 = 50 to simulate outliers. These latter are depicted in red in Figures 3.4(b-c-d) and an example of anomaly is visible in the bottom of Figure 3.4(a). As shown on these figures, outliers are lying in a separable distance from the 2-simplex encompassing the nominal data (note that in practice outliers may not necessarily live in a similar area but we prefer to use this simple situation for validation).

In further experiments on synthetic data we use different datasets with various sizes (N = {10 3 , 104 , 10 5 , 10 6 }) and endmembers ( = {3, 5, 7, 10, 12, 15}). As for the evaluation of anomaly detection procedures, for all datasets we choose to generate Υ = 20 anomalies from three different signatures.

Real data

For the application on real data we choose in this study to focus on 4-Vesta asteroid which is the second largest asteroid in the solar system and thus is of prime scientific interest in planetary exploration. 4-Vesta is though to be a protoplanet, i.e. a building block of terrestrial planets remnant of the early formation of the Solar System [START_REF] Christopher T Russell | Dawn at Vesta: testing the protoplanetary paradigm[END_REF]. Therefore understanding its properties is of prime importance to understand the initial material that accreted to form the Earth 4 . The number "4" in the name of the asteroid is due to the fact that 4-Vesta was the fourth asteroid to be discovered. It was observed recently by the Dawn mission [START_REF] Christopher | The Dawn Mission to Vesta and Ceres[END_REF], providing high resolution hyperspectral images of the whole surface. Moreover a very large number of meteorite samples in our collection come from this asteroid [START_REF] Harry Y Mcsween | Dawn; the Vesta-HED connection; and the geologic context for eucrites, diogenites, and howardites[END_REF], providing additional constraints on the minerals mixtures possibly encountered when processing the HS dataset. As the dataset available for 4-Vesta is concerned by all the limitations described in section 3.1, the use of an improved non-linear unmixing algorithm is essential to progress from basic indexes maps (as done up to now, e.g. [START_REF] Ammannito | Vestan lithologies mapped by the visual and infrared spectrometer on Dawn[END_REF]) to detailed maps, making the link between meteorites and compositions extracted from remote sensing possible. Because of the large size of the full dataset (as mentioned in section 3.1), here we choose to focus on a small subset located in the northern hemisphere. This region was chosen because Olivine, a mineral indicative of deep magmatic processes, was first detected in this area [START_REF] Ammannito | Olivine in an unexpected location on Vestaś surface[END_REF]. This mineral is never found alone in meteorites coming from Vesta-4 [START_REF] Andrew | Diogenites as polymict breccias composed of orthopyroxenite and harzburgite[END_REF][START_REF] Combe | Composition of the northern regions of vesta analyzed by the dawn mission[END_REF], implying the use of an algorithm able to decompose spectra as a combination of mixtures of materials. Following the procedure described in [START_REF] Clenet | A deep crust-mantle boundary in the asteroid 4 Vesta[END_REF], we processed Dawn VIR [START_REF] Cristina | The VIR Spectrometer[END_REF] images from the High Altitude Mapping Orbit (HAMO) 1 and 2 and from the Low Altitude Mapping Orbit (LAMO). Raw images were downloaded from PDS in level 1B and were calibrated using ISIS3 pipeline [START_REF] Anderson | Isis cartographic tools for the Dawn Framing Camera and Visual and Infrared Spectrometer[END_REF] and photometric parameters found in the literature [START_REF] Reddy | Color and albedo heterogeneity of Vesta from Dawn[END_REF][START_REF] Cristina | Vesta's mineralogical composition as revealed by the visible and infrared spectrometer on Dawn[END_REF][START_REF] Li | Photometric properties of Vesta[END_REF]. Additional in-house routines were used to filter bad pixels and to correct the geographical misalignment between the two detectors of the VIR instrument (following the method described in [START_REF] Ammannito | Vestan lithologies mapped by the visual and infrared spectrometer on Dawn[END_REF]). Finally each image was converted from radiance to I/F and projected using the spacecraft geometry information. Such process remains standard in planetary image analysis.

After removing channels known to be noisy, the final mosaic has a size of (842 × 327) pixels, with 383 spectral bands covering the range 0.55µm to 2.47µm. Nevertheless, despite the first-order filtering of bad pixels/spectra, some inconsistency still exist in the dataset and outliers can be observed, as illustrated in Figure 3.5. In this figure the typical shape of an expected spectra is depicted together with highly corrupted spectra due to the misalignment between the visible and near-infrared detectors, either showing a sharp peak when pixels fall on two distinct lithological units or an Heaviside-step shape signature when pixels fall at a boundary between sunlit and shadow.

Parameters setting for SAGA and SAGA+

To guarantee an optimal performance, internal parameters of SAGA and SAGA+ have to be tuned according to dataset characteristics. Apart from the number of atoms/endmembers and the sparsity level λ, these parameters correspond to kernel bandwidth (σ) and for SAGA+, to the error threshold (τ ) of the anomaly detection procedure (see algorithm in Figure 1). For this reason a number of tests have been taken in order to find the best parametrization of the models. As the result of these tests, we observed empirically that the best parameter values for linear mixtures are σ = 0.15 and τ = 0.2, for bilinear mixtures : σ = 0.15 and τ = 0.04 and for highly concentrating mixtures : σ = 0.15 and τ = 0.2. This setting enabled us to detect all true anomalies (true positives) through each run of the algorithm while keeping the number of wrongly selected anomalies (false positives) minimum. These parameters kept fix in all our experiments. It should be noted that these tests are meant to find the best parameters for a model regarding the dataset under consideration and thus not to be thought as generic. For the real application presented later we also use σ = 0.15 and τ = 0.2 since they are used for LMM and HCM. As will be shown in the next section, this setting produces acceptable results. The question of an optimal parameter selection for real datasets within the unsupervised scenario however remains an open question. We now turn to experiments.

Results on synthetic data

In this part we compare the results from SAGA and SAGA+ with other unmixing procedures on the synthetic datasets, presented in section 3.3.1. We first compare SAGA with N-Finder, VCA, and G-SiVM in order to evaluate their accuracy in endmember estimation in absence of any anomaly. For this we generated various datasets of size N = 10 3 , using different numbers of signatures = 3, 5, 7, 10, 12, 15 and based on the three aforementioned mixing models. Because our approach embeds some randomness in the initialization step, we depict the average of SAM errors computed on 100 runs. They are visible in Figure 3.6 in function of different number of endmembers for each algorithm. From this figure one observes that SAGA, N-Finder and VCA have similar performances unlike G-SiVMs whose accuracy is generally diminished for the case of intimate mixtures. This first experience where no outliers are involved shows the relative good performances of all techniques to extract consistent endmembers.

In a second experience we ran exactly the same experiments, this time however, we included Υ = 20 anomalies (generated from 3 extra signatures randomly selected from the dataset as explained in previous section) in the dataset. New results are depicted in Figure 3.7. As expected the overall error increases in all the algorithms except SAGA+ which includes an anomaly detector process. These simple experiments demonstrate the ability of our approach to accurately estimate endmembers while removing outliers. It should be noticed that the 41 signatures used in this setting are not extremely different in their overall profile shape. Nevertheless the results show that the presence of AD even at the situations where the anomalies are not dramatically different (therefore have small distance w.r.t the nominal data) from the nominal signatures can effect considerably. Let us also outline that in the last dataset, outlier signatures were considerably different from the nominal data, making it an ideal condition for the anomaly detection. In order to evaluate the sensitivity of SAGA+ when anomalies are more dispersed w.r.t data, several concentration levels (issued from Dirichlet distribution as mentioned in previous section) of anomalies have been generated: the higher the concentration level, the lower the dispersion of anomalies. An example of dataset is shown in Figure 3.8 together with endmember and anomaly estimation with SAGA+. As one observes in this figure, both anomalies and endmembers are accurately detected. In another test, the anomalies are assigned on the basis of 50 executions, where in each execution, a new concentration level in the range α = 1, . . . , 50 is used to generate anomalies (see equation (3.15)). In order to quantitatively evaluate the sensitivity of AD procedure to the level of concentration/dispersion of anomalies. In Table . 3.1 we report the average of Kappa values when comparing true number of anomalies to those detected by SAGA+ over synthetic datasets of size N = 10 3 , generated with various numbers of endmembers = 3, 5, 7, 10, 12, 15 and Υ = 20 anomalies generated from three different signatures. From this table one observes that when increasing the number of endmembers, the anomaly detection algorithm is still robust since κ > 0.7 is generally interpreted as a good separation accuracy in two class classification. The way Kappa reacts in function of the concentration of outliers is visually illustrated in Figure 3.9. Let us remind that when α grows, the separation between outliers and nominal data is more clear. From this figure one observes that when α = 15, we reach κ 0.8 which is a very good accordance. In order to provide a more visible understanding of the effect of anomaly dispersion on the accuracy of the AD algorithm, in Figure 3.10 we generated different datasets using the same endmember signatures and anomaly signatures. Like other tests 1000 samples generated using = 3 endmembers and Υ = 20 anomalies generated using 2 spectral signatures. For generation of anomalies different levels of concentration have set which are reflected in the parameter α = [START_REF] Nascimento | Vertex component analysis: a fast algorithm to unmix hyperspectral data[END_REF][START_REF] Kabwama | Intra-operative hyperspectral imaging for brain tumour detection and delineation: Current progress on the helicoid project[END_REF][START_REF] Chang | A fast iterative algorithm for implementation of pixel purity index[END_REF]. The result of the AD algorithm is then evaluated according to the ratio between True Positive (TP) and False Positive (FP). It is visible that how large dispersion (low values of α) can obscure the division line between nominal data and outliers, resulting in low accuracy of anomaly detection. In the contrary high value of α isolates the anomalies and thus it would be easier for the AD algorithm to detect them.

On the application of Overcomplete Dictionary

In order to illustrate the benefits of overcomplete dictionaries in estimating reliable endmembers, we generated a dataset composed of N = 10 3 nominal data on the basis of BMM using 4 original signatures. A number of Υ = 20 anomalies have been added from 2 different signatures. When running SAGA+ we choose = 10 atoms to represent endmembers on an overcomplete dictionary and we set a sparsity level λ = 4 since mixing were generated from 4 signatures.

Results are shown in Figure 3.11 where we depicted original endmembers and outliers (Figure 3.11(a)), a 2D slice of the generated dataset (Figure 3.11(b)), the = 10 estimated endmembers and associated outliers (Figure 3.11(c)) and the resulting 4 "averaged" endmembers (Figure 3.11(d)). In this last sub-figure, as we know exactly the location of the 4 pure endmembers, each of the ten estimated ones is associated to the closest real one and weighted averages, on the basis of their abundances, are represented. One can see the good adequacy of our approach when comparing these averaged endmembers with real ones in the top of Figure 3.11(a). It is also interesting to see the effect of overcomplete dictionary on the accuracy of the abundance coefficients. For this we generated datasets in a similar way as the previous step. For the projection of dataset on the estimated endmembers by Nfinder and VCA algorithms we used the Constrained Least Squares (CLS) [START_REF] Bioucas-Dias | An overview on hyperspectral unmixing: Geometrical, statistical, and sparse regression based approaches[END_REF] method. This time we evaluated the closeness of each row of the estimated abundance matrix with its closest counterpart in the original generated abundance matrix using the mean of SAM error calculated as:

(G, Ĝ) = min η 1 i=1 SAM(g η (i) , ĝi ) (3.16)
where is the number of endmembers and η : {1, . . . , } -→ {1, . . . , } is a function corresponding to all possible permutations. The result of this comparison is represented in the table 3.2 which shows considerable improvement in the accuracy of the abundance maps in the case of SAGA+ comparing to the other methods. Finally, we report a comparison of execution times in supplementary material provided in the appendix 7.

Results on real data: 4-Vesta Asteroid

In this experiment we consider the performance of SAGA+ on 4-Vesta and try to give a high level interpretation of the results which can be of interest from the planetary observation point of view. In practice, two difficulties arise when manipulating 4-Vesta data : 1) the real number of endmembers is unknown and 2) the assumption that some pure pixels are present in the original image does not hold. In order to cope with these difficulties, as mentioned earlier, we exploit the intrinsic characteristics of overcomplete dictionary learning approach which relaxes both mentioned difficulties by selecting out a large number of endmembers (referred to as the dictionary atoms) from the image. From the geological point of view, the large number of selected endmembers/atoms might be interpreted as different compositions of materials in a scene. In practice for the case of 4-Vesta we choose the dictionary of size = 10. Because of the intrinsic characteristic of overcomplete dictionary for re-grouping the similar endmembers, this parameter is not crucial but has to be set such that the number of atoms of the dictionary is larger than the number of expected "pure" materials. The estimated endmembers are represented in the bottom of Figure 3.12 (denoted P0 to P9) while the top of this figure depicts the spatial positions of the respective estimated endmembers (in green) together with the outliers (in red). A specimen of outliers spectra signature (with a sharp peak) is represented in the bottom of Figure 3.5.

By comparing these two plots (Figure 3.12 and Figure 3.5), one can observe that the estimated endmembers signatures are generally in accordance with the expected nominal signature profiles (laboratory spectra related to the lithologies on the northern hemisphere). Although one observes some unwanted signatures (mainly contributing to the shadowing effect), it is of prime interest to note that the absence of signatures corrupted by noise is significant 5 . Indeed, endmembers related to unwanted spectra can be discarded: P1 and P3 have a low reflectance and correlate perfectly with the most dark (shadow) areas, while P2, P6 and P7 show a shift between the two detectors, with a very small amount of noise. The five other endmembers are of interest regarding mineralogical detections. P4 and P8 show characteristic signatures of pyroxenes and represent the surface background. Some outcrops, related to P0, exhibit similar signatures with stronger absorptions, in relation with impact craters (i.e. fresh rocks) and very good illumination conditions. The last two endmembers P5 and P9, detected in small spots in and around craters, show a broader absorption feature in the 1µm range. One can see the whole set of endmembers and their respective abundance maps in Figure 3.15. From these plots, it is worth noting that the spatial consistency of associated maps matches well with real physical structures in 4-Vesta, which is a very good property.

As a more physical comparison, one can also compare the spectra in Figure 3.12 with the laboratory measured spectra of Olivine and Orthopyroxene which is illustrated in Figure 3.13. As can be observed in this figure, the broader 1µm absorption seen in P9 and P5 can be correlated to an increased olivine content in the rocks, in accordance with observations made by [START_REF] Ammannito | Olivine in an unexpected location on Vestaś surface[END_REF] or [START_REF] Ruesch | Detections and geologic context of local enrichments in olivine on Vesta with VIR/Dawn data[END_REF]. From these observations, SAGA+ is thus able to extract from the hyperspectral image, spectra that are a combination of mixtures of materials, representative of distinct geological units. Additionally, using the extracted endmembers, the outcrops enriched in Olivine (i.e. spots where pixels in which P9 and P5 are needed in the unmixing) can be mapped (see Figure 3.14) and found localization are similar to those obtain with other techniques [START_REF] Ruesch | Detections and geologic context of local enrichments in olivine on Vesta with VIR/Dawn data[END_REF]. Using an anomaly detection procedure, one observes that signatures whit magnitudes higher than 1 and those with sharp peaks are not selected as endmembers.

Let us note that the shift in absorption center at 2µm is due to the chemical composition of orthopyroxene (olivine does not have absorptions in this wavelength range) and reflects the fact that the orthopyroxene used as a reference does not have the same chemical composition as the one observed on 4-Vesta's surface. Lastly, in order to represent the benefits of the anomaly detection procedure, we have applied the same method without the outlier removal process. Extracted endmembers and associated abundances are available in the supplement material provided in the appendix 7. By comparing them with Figure 3.15, we clearly observe corrupted endmembers yielding less consistent abundance maps. This experience reveals the importance of the anomaly detection process introduced here.

Conclusion

In this chapter, a joint procedure for anomaly detection and spectral unmixing has been presented for hyperspectral data. Endmembers and anomalies are extracted from the vertexes of the simplex that embeds data, this simplex being computed in a feature space to deal with non-linearities. Abundances are computed on the basis of sparse projections of spectra onto extracted endmembers. To deal both with the internal heterogeneity inside classes and with the fact that pure elements are not necessarily present in the data, overcomplete dictionaries have been used. We evaluated and compared our approach with linear and non linear unmixing methods on synthetic datasets generated from real spectra where linear and nonlinear mixing procedures have been applied. Anomalies have also been introduced in the data. These experiments have revealed the efficiency of our technique to deal with anomalies and endmember detection at the same time. Computational aspects have also been discussed and have revealed the ability of our algorithm to easily scale up to large datasets. Lastly, experiments have been performed on real data issued from the 4-Vesta Asteroid. From a geological point of view, we have shown that anomalies, endmembers and mixing abundances are consistent, demonstrating the efficiency of our approach for real applications.

Part III Optimal Transport for Unmixing

Chapter Summary

In this chapter a new supervised model for spectral unmixing based on Optimal Transport(OT) theory is presented. This theory, very well adapted to the manipulation of distributions, is exploited here on hyperspectral data viewed as distributions instead of a multivariate vector. Optimal transport consists in searching for a unique plan which transports a distribution x 1 onto x 2 within the minimum cost. The associated distance measure for this is called Wasserstein distance. In order to calculate the abundance values for each pixel in X, we project the pixel spectrum onto each endmember (also presented as distribution) in our dictionary E using OT-based projection (called Bregman projection), while an entropy-based regularization term promotes the a priori knowledge regarding the inter-endmember relationships. In this framework the information reflecting distance relationships along the spectral dimension (in the case of hyperspectral spectra it will be wavelength values) of data is being used as an initial cost function. Also the information regarding the material dependencies of the endmembers of E will be provided as another cost function. Our algorithm is also able to deal with endmember variabilities among the values of E. In 4.1 after a short recap on the notion of spectral unmixing, the necessity of having representative distance measure, being able to reveal the characteristics of data domain is mentioned which consequently follows with the introduction of the Wasserstein distance. In 4.2 Optimal transport, Wasserstein distance and the Entropy based regularization is being explained along with a closed-form solution for bregman projection. In 4.2.3 the application of two tailored cost functions within the optimization procedure is discussed and finally in 4.3 the preliminary results of applying this algorithm on a portion of 4-vesta dataset is represented.

This chapter corresponds to an article published in WHISPERS 8th workshop on Hyperspectral image and Signal Processing: Evolution in Remote Sensing, and presented in Aug 2016, Los Angeles, United States.

OT considers spectra as probability distributions. A pixel spectrum therefore will be normalized along spectral values. Normalization makes the analysis less sensitive to the total power of spectra in each pixel. This improves robustness against shadows or other large radiance changes and thus can prevent degenerate solutions.

In order to tailor a cost matrix C0 in alignment to the characteristics of the dataset, it is possible to build C0(D × D) as the square euclidean distance over the spectral values. This clearly reflects the characteristic of the spectra and the level of (dis)similarity among them.

It is possible to manually construct the cost matrix C 1 in order to incorporate the information regarding the endmembers belonging to the same material group. Two endmembers belonging to the same material share a very low cost with the corresponding material in α 0 , C1(i; j) = 0

Introduction

Unmixing is an active field in hyperspectral data analysis and consists in estimating the abundance of pure materials (named endmembers) that are part of a pixel spectrum. These pure elements are either estimated from the image (unsupervised methods) or directly given in a prior knowledge (supervised approach).

Spectra comparison. & optimal transport

In all techniques mentioned in chapter 2, a metric between spectra is required. Several criteria have been used in a hyperspectral context (spectral angle mapper and Euclidean minimum distance being the two most popular ones), however, few approaches consider spectra as probability distributions in the spectral domain. We consider in this work an original metric, designed between probability density functions to compare spectra. This metric, based on optimal transport (OT) is called Wasserstein distance, and can be designed to be mostly sensitive to shifts in the frequency domains. To exploit OT, each spectrum should be normalized (the integration along spectral values should be 1) in order to be viewed as a probability distribution. By doing this normalization we make the analysis insensitive to quantities of materials in each pixel and we rather focus on proportions. This enables for example to improve robustness against shadows or other large radiance changes and as discussed in [START_REF] Sigurdsson | Endmember constrained semi-supervised hyperspectral unmixing[END_REF], can prevent degenerate solutions.

Outline of the chapter. In the next section, we explain the basics of OT together with the entropy-based regularization term which is used for calculation of abundance coefficients. Experiments on the challenging 4-Vesta asteroid dataset are then shown and discussed in section 4.3.

Unmixing using priors and iterative Bregman projections

We first introduce optimal transport, then we present in section 4.2.2 a regularized version more efficient from a computational point of view and section 4.2.3 proposes an unmixing model based on optimal transport.

Optimal transport of spectra

Optimal Transport (OT) can be seen as the search for a plan that moves (transports) a spectrum µ 1 onto another spectrum µ 2 with a minimum cost. This cost is usually related to a metric of the embedding space. In our case, we will generally use a squared Euclidean distance L 2 2 as this cost, but it should be mentioned that other cost functions, specific to a given problem, might be used. In the relaxed formulation of Kantorovitch [START_REF] Kantorovich | On the translocation of masses[END_REF], OT seeks for an optimal coupling that can be seen as a joint probability distribution between µ 1 and µ 2 . In other words, if we define Π(µ 1 , µ 2 ) the space of probability distributions over R 2 with prescribed marginals µ 1 and µ 2 , the optimal transport is the coupling γ ∈ Π(µ 1 , µ 2 ) which minimizes the following quantity:

W c (µ 1 , µ 2 ) = inf γ∈Π(µ 1 ,µ 2 ) R 2 c(x 1 , x 2 )dγ(x 1 , x 2 ) (4.1)
where c(x 1 , x 2 ) is the cost between x 1 and x 2 (issued from distributions µ 1 and µ 2 respectively). In the discrete versions of the problem, i.e. when µ 1 and µ 2 are sampled by a sensor and expressed as vectors of R d , the previous equation reads:

W C (µ 1 , µ 2 ) = min γ∈Π(µ 1 ,µ 2 ) < γ, C > F , (4.2)
where < ., . > F is the Frobenius dot product and C ≥ 0 is a cost matrix (of size d×d), which gathers all the costs for moving from bin i to bin j, and γ is now a matrix of size d × d, with marginals defined as µ 1 and µ 2 . We note that this distance is also known as the Earth Moving Distance in the computer vision community [START_REF] Rubner | A metric for distributions with applications to image databases[END_REF]. Solving for equation (4.2) is a simple linear programming problem with equality constraints, but its dimensions scales quadratically with the number of samples in the spectra. This might generate computational problems when the dimension of distributions is high, as in hyperspectral data. For this reason, we now consider a regularized version of the problem which has the extra benefit of being faster to solve.

Entropy regularized Optimal Transport and Bregman projections

The idea is here to put a regularization term over γ that controls the smoothness of the coupling through its entropy. The entropy regularized version of the discrete optimal transport problem reads [START_REF] Cuturi | Sinkhorn distances: Lightspeed computation of optimal transport[END_REF]:

W C, (µ 1 , µ 2 ) = min γ∈Π(µ 1 ,µ 2 ) < γ, C > F -h(γ), (4.3)
where h(γ) is the entropy of γ, reading

h(γ) = - ij γ ij log γ ij = -< γ, log γ > F ,
if log is applied component-wise to γ. Denoting the Kullback Leibler divergence (KL), i.e. KL(γ, ρ)

= ij γ ij log γ ij ρ ij =< γ, log γ ρ > F
where again the log and the division are taken component wise, one can establish the following link between OT and Bregman projections: Proposition 4.2.1 [OT as a Bregman projection [119, Eq. ( 6)]] The minimizer γ of (4.3) is the solution of the following Bregman projection

γ = arg min γ∈Π(µ 1 ,µ 2 ) KL(γ, ζ), (4.4)
where ζ = exp(-C ).

Interestingly enough, the entropy regularized version of OT admits a simple algorithm for the resolution of this problem, based on the successive projections over the two marginal constraints. More specifically, solving (4.4) when one of the two marginals is unknown admits a closed form solution.

Proposition 4.2.2 [Closed form solution of OT with one unknown marginal [119, Eq. ( 7)]] For an undefined µ 2 , γ is solely constrained by one marginal µ 1 and is the solution of the following closed-form projection:

γ = diag µ 1 ζ1 ζ (4.5)
where again the division has to be understood component-wise. Consequently, by iterating over successive Bregman projections on the set of the two marginal constraints, one can derive an efficient algorithm to solve for the OT problem. We refer the reader to [START_REF] Benamou | Iterative bregman projections for regularized transportation problems[END_REF] for a more complete treatment on this subject. Finally, when several OT distances are implied in the same optimization problem, this iterative projection strategy is still applicable. We use this property in our unmixing model, described in the following.

Unmixing model

We start here by recalling our initial assumptions. We assume a linear mixture for each spectrum µ : µ = Eα where matrix E of size d × q contains q (known) endmembers of dimension d and vector α > 0 of size q × 1 has to be estimated under the sum-to-one constraint, i.e. α T 1 = 1. As such, it corresponds to discrete probability distribution function or histogram. From the image analysis, we seek to extract p abundances per pixels of distinct materials of interest. In practice the number of endmembers q exceeds the number of proportions p we aim at estimating (q ≥ p). This allows to consider several distinct versions of spectra for the same material of interest, thus accounting for endmembers variability. In addition we assume that we have a prior knowledge α 0 (p-dimensional) over the abundances.

The unmixing of a spectrum µ is the solution of the following optimization:

α = arg min α W C 0 , 0 (µ, Eα) data fitting +τ W C 1 , 1 (α, α 0 ) prior . (4.6)
As it can be seen, the problem is subdivided in a data fitting term (which tries to find the best decomposition from observations) and a regularization term (which enforces the compliance of the solution with the prior) balanced by parameter τ ∈ R + . From these two terms, regarding equation (4.3), two cost functions C 0 and C 1 are involved. Matrix C 0 is linked with a metric in the spectral dimensions of the data (and is usually related to the L 2 -norm along frequencies) while C 1 is a cost matrix of size q × p with a specific structure which contains information about the classes of endmembers. Two endmembers related to the same element will share a very low cost with the corresponding material in α 0 but a higher cost with other materials. An easy way to handcraft C 1 is to set C 1 (i, j) = 0 if endmember i is related to material j, and 1 elsewhere. Finally, two regularization parameters 0 and 1 related to the regularized version of the transport have to be set. The resolution of the optimization problem in (4.6), inspired from the weighted Wasserstein barycenter problem [START_REF] Benamou | Iterative bregman projections for regularized transportation problems[END_REF], is done using iterative Bregman projections. Algorithm 3 presents the corresponding procedure in pseudo-code.

Algorithm 3 Supervised unmixing procedure with OT

1: Input: µ 2: Parameters: α 0 ,E,τ ,C 0 ,C 1 , 0 , 1 ,maxIter 3: cpt ← 0, err ← ∞ 4: α ← 1/q 5: ζ 0 ← exp(-C0 0 ), ζ 1 ← exp(- C1
1 ) 6: while cpt <maxIter and err>threshold do 7:

ζ 0 ← diag( µ ζ01 )ζ 0 8: ζ 1 ← diag( α0 ζ11 )ζ 1 9: µ tmp ← exp(log(ζ 0 1) + τ log(Eζ 1 1)) 10: ζ 0 ← ζ 0 diag( µtmp ζ T 0 1 ) 11: ζ 1 ← ζ 1 diag( E T µtmp ζ T 1 1 ) 12:
err← ||α (cpt) -α (cpt-1) || 2 13: cpt←cpt+1 14: end while return α

Experiments and Analysis

4-Vesta properties and dataset

A mosaic of Vesta-4 asteroid is visible in Figure 4.1. The image is taken by visible and infrared mapping spectrometer (VIR) which comes with two arrays of detectors: One array covering the visible range (between .25µm and 1.05µm) with spectral sampling of 1.8nm and the other covering the infrared range (between 1.0µm and 5.0µm) with spectral sampling of 9.8nm. The total number of considered bands after removing the noisy ones are d = 383 covering the range of .55µm to 2.4µm . Although there is no exact information about the components of the crust and mantel of this protoplanet asteroid, some expert knowledge is available (see [START_REF] Ammannito | Olivine in an unexpected location on Vesta's surface[END_REF] and [START_REF] Combe | Composition of the northern regions of vesta analyzed by the dawn mission[END_REF]).

From this we know that the major surface is dominated by three rocks: Eucrite, Howardite and Diogenite. It is believed that roughly the southern hemisphere of 4-Vesta is Deogenite-rich while the northern hemisphere crust is mainly composed of Eucrite and different types of pyroxene with recently detected clusters of Olivine.

Experimental setup

We used a total number of q = 20 spectral signatures of different lithologies, found ubiquitously in the composition of the meteorites, collected on Earth. we could categorize these signatures in 3 groups: Olivine, Orthopyroxene, Eucrit. Consequently we set p = 3 in order to identify "maximally pure" endmembers which are believed to exist within the scene. As the consequence, matrix C 1 involved in equation (4.6) will be of dimension 20 × 3. As mentioned in previous section, this cost matrix contains information about the classes of endmembers (two endmembers related to the same material will share a very low cost with the corresponding material in α 0 but a higher cost with other materials). In particular we set C 1 (i, j) = 0 if endmember i is related to material group j, and 1 elsewhere.

As for the cost matrix C 0 in relation (4.6), we simply use the L 2 2 norm between each spectral value. The regularization parameters 0 and 1 have been set to 1 and 10 3 and the balance parameter is τ = 0.9 . This is set manually from data but more advanced approaches with cross-correlation could be used. For setting prior knowledge of the expected endmembers different combinations of non-negative 3-values for α 0 with non-negativity and additvity constraints were tested. By α 0 we reflect our prior expectation of the proportions of the 3 expected groups of lithologies.

Unmixing results and discussion

We made some experiments with setting different combinations of priors for each group of materials. By this we observe abundance maps with some nuances. These are represented in Figure 4.2(a). It is interesting to observe that in all cases, spatially consistent areas extracted and correspond to existing spatial patterns of rocks. This is in contrast to the unmixing results we achieved by using constrained least squares (CLS) method which is illustrated in Figure 4.2(b) b). As for the quantities of pure elements, our results are consistent with very recent studies [START_REF] Combe | Composition of the northern regions of vesta analyzed by the dawn mission[END_REF] on this part of 4-Vesta. In particular the fact that different compositions of Orthopyroxene are dominant while Olivine patches are spatially much more isolated and sparse goes well with the latest observations on the composition of the materials of this area. As the role of priors one can see the change in susceptibility/sensitiveness towards a particular material when changing the prior value over that material. This in particular is visible w.r.t the Olivine maps.

Discussion/Conclusion

In this work we have proposed a hyperspectral unmixing approach based an optimal transport. This metric, devoted to distributions, is mostly sensitive to shifts in the frequency domains and therefore differs from most usual criteria as spectral angle mapper or Euclidean distance. The problem of spectral variability and the confusion in selecting the best endmembers among the existing dictionary is addressed through an overcomplete dictionary. This results in a specific cost-function to optimize where both observation and regularization terms are based on OT. Our primarily experiments on the real 4-Vesta dataset show a good accordance of our results with the unmixing issued from the latest geological findings about the composition of lithologies on this asteroid.

In the next chapter we are going to see the application of the two proposed unmixing procedures on a hyperspectral image for urban monitoring purpose.

Part IV

Applications of Spectral Unmixing in Urban Context

Chapter Summary

The goal of this chapter is to apply and demonstrate the performance of the two spectral unmixing algorithms we introduced in the previous chapters within the context of land surface temperature (LST) estimation. An example of such temperature estimation is visible in above illustration. The global idea of the project associated with our work is to use the abundance maps generated by the unmixing procedures in order to accurately calculate the surface area of different urban components (roads, parks, building roofs with specific materials etc.). This information will be combined with emission models specific to each material to provide accurate maps of surface temperature for different regions of the city. In this chapter and for the sake of this thesis work we focus on the unmixing part of this project. Experiments have been performed on a dataset issued from the city of Madrid (Spain) using the two algorithms introduced in this thesis. The generated abundance maps have also been compared with the results of other state-of-the-art techniques. In the first section we outline the goals of the associated project (entitled CATUT, French acronym for Mapping Urban Temperatures), then we provide a description of the data, (acquisition and characteristics). In the second section we introduce the way three endmember dictionaries which are going to be used for our unmixing task are acquired and preprocessed. We also set six strategies for conducting our unmixing experiments. In the last section we provide our preliminary results. This chapter is an applicative part under development. First results are introduced here. PORT 1 covers the VIS/NIR range. Its bands have relatively low spectral resolution (28 -30 nm wide) while the spectral range is continuous from 443(nm) up to 1025(nm).

PORT 2A covers the SWIR range with an isolated band centered at 1.6(µm) with spectral resolution of 90(nm).

PORT 2, has a set of continuous, with relatively high spectral resolution (17 -18(nm) wide) laying between 1907(nm) and 2558(nm), well suited for soil/geologic studies. Information regarding the spectral range and the rate-of-change of spectral resolution across the bands and for the first three ports is visually depicted in Figure 5.3. Within this chapter, we will merely rely on the information provided by these three PORTs for the evaluation of our unmixing algorithms. All radiometric calibrations have been performed using at-lab measurements. In Figure 5.4 a complete mosaic captured through one north-south flight over the area is depicted.

SPECTRAL CHARACTERIZATION OF URBAN SUR-FACES

In order to build an urban material spectral library, through the DESIREX 2008 project, spectral signatures over different natural and artificial surfaces were acquired with the ASD FieldSpec-3 by INTA Remote Sensing Area. ASD FieldSpec-3 is a portable and general-purpose spectroradiometer which is illustrated in Figure 5.5.

The Visible/Near Infrared (VNIR) portion of the spectrum, the 350 -1050(nm) wavelength domain, is measured by a 512-channel silicon photodiode array. Each channel is geometrically positioned to receive light within a narrow 1.4(nm) bandwidth. The VNIR spectrometer has a spectral resolution (FWHM of a single emission line) of approximately 3(nm) at around 700(nm). The short-wave Infrared (SWIR), also called Near Infrared (NIR), portion of the spectrum is acquired with two scanning spectrometers. These differ from the array used in the VNIR in that they measure wavelengths sequentially, rather than simultaneously. Each spectrometer consists of a concave holographic grating and a single thermoelectrically cooled indium gallium arsenide (InGaAs) detector. The gratings are mounted in a common shaft that oscillates with a period of about 200 milliseconds (100ms/scan). Unlike the VNIR, each SWIR spectrometer has only one detector, which is exposed to different wavelengths of light as the grating oscillates. The first spectrometer (SWIR1) measures light between about 900 -1850(nm); the second (SWIR2) covers the region between about 1700 -2500(nm). The sampling interval for each SWIR region is about 2(nm), and the spectral resolution varies between 10(nm) and 12(nm), depending on the scan angle at that wavelength.

In particular Pavement, Asphalt (Parking), Red Brick, Concrete, Asphalt (Composition), Metalic Mesh, Red Tile, Asbestos, Steel, Granite, Green Grass and Bare Soil were among materials carefully measured during the campaign. Figure 5.6 represents an exemplary image for each one of interesting materials.

On the valid choice of Endmember dictionary

Considering the difference between the spectral characteristics of Desirex dataset and the ASD-based endmember dictionary, also due to the fact that for the unmixing of the dataset we are considering supervised and unsupervised procedures, one faces different possible strategies to prepare the required endmember dictionary and conduct the unmixing. In particular one can:

1. Pre-process the ASD dictionary and use it for the supervised OT-unmixing algorithm.

2. Manually extracts endmembers from the Desirex dataset and use them for the supervised OT-unmixing algorithm.

3. Using the Exact endmembers estimation part of SAGA+ to build the endmember dictionary and then use it for OT-unmixing to produce the abundance estimations.

4. Pre-process the ASD dictionary and use it for the SAGA+ algorithm (therefore skipping the endmember extraction part in SAGA+).

5. Manually extracts endmembers from the Desirex dataset and use them for SAGA+ algorithm (therefore skipping the endmember extraction part in SAGA+).

6. Taking a complete unsupervised strategy and let SAGA+ to find the endmembers from the Desirex dataset and estimate the abundances.

Considering the fact that for each of the above mentioned strategies one needs also to set a few number of model parameters (choice of the kernel in SAGA+, level of sparsity of the representations, grouping and weighting values of the endmembers in OT-unmixing etc.), an extensive amount of tests needs to be performed and comparisons to be taken in order to find the best unmixing scenario and validated endmember dictionary. This complete work is in progress and in this PhD, we will represent one unmixing result per strategy using a fix parametrization setup in order to let comparisons between performances of strategies possible. Also, in order to provide a comparison between the proposed methods with other state-of-the-art techniques we consider the performance of the sequential least squares programming (SLSQP) on the same dataset and using three endmember dictionaries presented above.

In the remaining part of this chapter we first briefly introduce the preprocessing procedure for the ASD endmembers and introduce the approach to manually generate dictionary of endmembers under the supervision of human experts in the field of geoprocessing. We also introduce a dictionary of endmembers generated by the endmember estimation part of SAGA+. Having these three endmember dictionaries, in a second part the results of the unmixing procedures (using these three dictionaries) according to the above mentioned strategies are being presented.

ASD Preprocessing

Before using the data and in-lab signatures for conducting unmixing we need to guaranty that our dataset (issued from airborne acquisitions) and the endmember dictionary (issued from lab measurements), are representing the same wavelength ranges. Their spectral resolutions need also to be comparable. When comparing Desired and ASD spectra, we first observe that the wavelength range of DESIREX is slightly superior than 2.5(µm) while the last wavelength registered for each ASD signature is 2.5(µm). This forces us to remove the last three bands of the DESIREX dataset, leaving us with 60 spectral bands.

ASD cut-off

In Figure 5.9 are depicted examples of spectra issued from DESIREX and ASD. On this figure, one can observe a large difference between their wavelength coverage across the electromagnetic field and also their spectral resolution. Therefore some ranges in ASD (out of the green area in Figure 5.9) must not be considered for the means of unmixing Desirex data. By cutting those excessive ranges off from the endmember signature profiles we will have a new endmember dictionary consisting of endmember signatures with 1185 bands, represented in Figure 5.10.

ASD downsampling

Once DESIREX and ASD endmember dictionaries are in accordance from the point of view of their wavelength range coverages, we need also to make them comparable in terms of spectral resolution in order to reduce the spectral dimension of the ASD endmembers from 1185 to 60 bands. Figure 5.11 represents the endmember signatures of Figure 5.10 after being donwnsampled. One can see that by uniformly donwnsampling the endmember profiles related to PORT1 and PORT2, the overal shape of the endmember signatures remain intact.

Even after performing the above mentioned preprocessing steps, one can see in Figure 5.11 the large level of incompatibility in terms of reflectance, between the overall signature profile of the Desirex dataset across the three PORTs and the overall signature profile of the ASD endmember dictionary. Some scaling parameters related to different sensors have been used which are different from other sensors. This disturbs the shape of all spectra along the wavelength range. Unfortunately in the present study these parameters were not available. By knowing this limitation and assuming that these parameters must be available to bring ASD dictionary and the Desirex in accordance, for the sake of comparison of our unmixing models performances we use the current dictionary of ASD endmembers (adjusted and downsampled) as one of the available sources of knowledge regarding the materials within the scene and compare their results with other possible dictionaries. Let us now introduce another endmember dictionary, prepared based on expert knowledge.

Expert Endmember Extraction

Another widely used approach towards hyperspectral unmixing is to find maximally pure pixels from a dataset and take them as endmembers. This approach can be taken manually and with the expert knowledge about the distribution of materials within the scene. Considering the fact that the ASD endmember dictionary, even after the pre-processing procedure does not look completely in accordance to to the type of materials that are common in the desirex dataset, we also used a manual approach of endmember extraction based on a human expert observation. Pixel vectors related to highly distinct materials within the scene have been extracted to be used as a dictionary of endmembers. In order to account for material variability, for each material two representatives endmembers were selected in order to build-up a dictionary of 24 signatures. This Endmember dictionary is represented in Figure 5.12. .12: Dictionary of the manually extracted endmembers. Based on expert observations, pixel vectors related to highly distinct materials within the scene have been extracted to be used as a dictionary of endmembers. In order to account for material variability, for each material two representatives endmembers were selected in order to build-up a dictionary of 24 signatures. Figure 5.13: Using the exact simplex volume maximization (E-SiVM) procedure a dictionary of 24 atoms with maximum variability among their signature profiles has been generated. One can observe even larger variability among endmember signatures while comparing this dictionary to the one extracted manually.

Endmember Extraction with E-SiVM

In order to evaluate the performance of SAGA+ in automatic estimation of pure endmembers, we use the exact simplex volume maximization (E-SiVM) procedure to generate another dictionary of endmembers. For this we chose an RBF kernel with kernel parameter σ = 2 × 10 -5 in order to get maximum variability among the signatures of 24 different pixels from a window of size 601 × 601 pixels. The choice is based on our results on various experimental tests of ESiVM on the data. Of course an automatic selection of these parameters (for example using cross validation) would be better. This is the scope of a current work. This endmember dictionary is represented in Figure 5.13. By comparing the Figures 5.12 and 5.13 one observes that the estimated endmembers with E-SiVM procedure can represent even larger diversity among the existing materials in comparison to the manually selected dictionary. This difference is more obvious when comparing the performance of ESiVM dictionary with a dictionary of randomly selected endmembers from the same dataset(that is shown in Figure5.14). In order to make a visual comparison in Figure 5.15 we provided two illustrations consisting of the superposition of the Expert dictionary with both randomly generated dictionary and the one generated by ESiVM procedure. In order to have a quantitative analysis, one could select randomly a large number of dictionaries (consisting of thousands of signatures) and compute distances between these dictionaries and endmembers (minimal distance between endmembers and all permutations of spectra in dictionaries). Unfortunately, with 24 endmembers, this would require 24! ∼ 10 23 distance computations between endmembers which is in practice impossible. Therefore only a visual inspection can be done here. As can be seen in Figure 5.15, ESiVM dictionnary seems to be more in accordance with the Expert one (more variability is captured) than random ones. It is indeed obvious that the variability in endmembers extracted by ESiVM is higher. In the next section we introduce preliminary results of our two spectral unmixing algorithms according to the six strategies mentioned above. Results provided from a state-of-the-art procedure called sequential least squares programming (SLSQP [122]) are also included in order to compare the performance of our algorithms from the point of view of computational time and the accuracy of the unmixing.

Unmixing Experiments on Desirex

Experiment conditions

Following the mentioned strategies in section 5.2, we conduct the abundance estimation tests in three phases. In the first phase (5.3.2) we run a series of unmixing tests by the OT-unmix procedure (Strategies 1, 2 and 3) while endmember dictionary is provided through one of the three mentioned ways: 1) The ASD dictionary 2) the dictionary manually extracted from the AHS Desirex dataset and 3) The estimated dictionary from AHS Desirex and using E-SiVM. For the sake of clarity from now we refer to them as ASD, Expert and the ESiVM dictionaries. As for the C1 cost matrix, we chose the 24 atoms of the dictionary to be grouped based on K-means clustering results into 5 different groups.

In the second phase (presented in section 5.3.3) we repeat the same procedure with the sequential least squares programming (SLSQP), in order to provide a comparison between two supervised scenarios. This leads to Strategies 7, 8 and 9. Finally in the third phase (presented in section 5.3.4) we chose to test the kernelbased SAGA+ while using the three dictionaries. This leads to Strategies 4, 5 and 6. For all the experiments in this phase we choose an RBF kernel with the kernel parameter σ = 2 × 10 -5 (with which SAGA+ produces mostly distinguishable abundance maps). These strategies are synthesized in table 5.2.

Experiments have been performed on a window of 601 × 601 pixels extracted from the full Desirex mosaic. This window is assumed to contain a diverse number of interesting urban materials and is illustrated in Figure 5.16.

K-means clustering for endmember grouping and estimation of the level of sparsity

For these preliminary experiments, we decided first to group the 24 endmembers on several clusters in order to evaluate our capacity to unmix the main urban structures. To be consistent, the level of sparsity λ of the SAGA+ procedure will be set to this number of clusters.

The choice of the ideal number of clusters is based on the global inertia η m written for a clustering in m classes as:

η m = m i=1 k j=1 w (i,j) ||x (i) -µ (j) || 2 2 (5.1)
where i is the sample index, j represents the cluster whose centroid is µ. Parameter w (i,j) = 1 if the sample x (i) is in cluster j and equals zero otherwise. By increasing the number of clusters m one can observe a decrease in η m . The principle is to find the optimal number of clusters where the increase in m does not change the distortion in a considerable amount. This principle is usually referred as the Elbow test. In practice we found that a grouping in m = 5 clusters was an optimal choice. The results of this test performed on two endmember dictionaries (Expert and ESiVM) are visible in Figure 5.17 where one can observe the decreasing rate of inertia (η) when the number of clusters grows. All tests performed using an intel(R) Core(TM) i7 2.4 GHz platform with 16 GB RAM.

Let us now turn to the analysis of abundance maps. In the first phase of experiments the supervised OT-unmixing algorithm is used for producing abundance maps using the three mentioned endmember dictionaries. This comprises strategies 1 to 3 represented in the first column of table 5.2. In order to parametrize the model, for these tests we defined the C0 cost matrix in alignment to the spectral characteristics of the dataset. In particular we build C0(63 × 63) as the square euclidean distance over the spectral range of the AHS Desirex dataset. This cost matrix is illustrated in Figure 5.18 has been used in strategies 1 to 3. As we saw in chapter 4 this cost function gathers all the costs for moving one spectral value to another across the signature spectra for computing the differences between two signatures µ 1 and µ 2 .

Abundance Estimation with OT-Unmixing

As for the C1 cost function which represents the relationship of the endmembers belonging to the same group (introduced in 4.2.3), we chose to use the clustering information of each endmember dictionary (presented in section 5.3.1). C1 provides a priori information regarding the groups to which different signatures of the endmember dictionary belong. Here and for each strategy, the total number of signatures within each dictionary has been associated to 5 groups, assumed based on Also in order to give a visual representation of the characteristics of the signatures in each group and comparing similar groups over the two dictionaries Expert and ESiVM, in Figure 5.20 the result of clusterings of these dictionaries are illustrated. By focusing on the intensity values of the first 21 bands (as the result of PORT1) of the endmember signatures, pairs of groups with maximum similarity among their signature envelopes are illustrated in the same colors.

By matching the information of table 5.3 and Figure 5.20, one observes that the G1 and G2 groups over the Expert dictionary embed more signatures. As for ESiVM, these groups are G2 and G4. This means that the endmembers of these groups share some common characteristics thus, most probably referring to similar materials. The same analysis goes for the pairs (G4,G3) and (G5,G5) across the two dictionaries.

As for the vector of priors a 0 , we chose to set their value as g i t where g i is the total number of endmembers belonging to the ith group and t is the total number of atoms in each dictionary. As mentioned in section 4.2.3, this prior vector defines weights on the relationship of the endmembers belonging to the same group. Having these parameter settings fix for each strategy we made three experiments on the mentioned As expected, the estimations based on the ASD dictionary are not very informative, though they exhibit interesting spatial structures. Unfortunately, because of the lack of consistency between ASD and DESIREX spectra, the quality of such abundances is limited. Regarding the other two experiments using the Expert and ESiVM dictionaries, interesting patterns are visible in each case with different degrees of emphasis on particular structures depending mainly on different ways of grouping the atoms of each dictionary. It should be notified that in general the abundance maps produced based on the ESiVM dictionary, which is the result of automatic endmember estimation of our SAGA+ procedure, seems to outperform the results of unmixing based on the Expert dictionary. It is interesting also to notice the relative similarity between the couple (G1, G2) over the Expert dictionary and (G2, G4) over the ESiVM dictionary, which indeed share similar patterns in representing the intensities.

In the absence of expert interpretation of the maps, one could conclude that our abundance maps seem to reveal consistent structures in this difficult context, which is a good starting point for further investigations. These results will be compared with state-of-the-art techniques in the next phases.

Abundance Estimation with SLSQP

In the next phase and in order to provide a comparison between the performance of OT-unmixing with a common state-of-the-art technique, we compare our technique with the sequential least squares programming (SLSQP) to perform the inversion procedure on the same dataset and using the same three available dictionaries. The results of these tests are depicted in Figure 5.22.

Comparing the results with those in Figure 5.21 one can easily observe that the abundance maps produced by OT-unmixing are of much higher spatial information. It should also be considered that the computation time of SLSQP on this data was higher than OT-unmixing by the order of 3, which means from about 1.5 hours in the case of OT-unmixing to more than 4 hours in SLSQP.

Though it is complicated to make a quantitative interpretation of these results, when comparing Figures 5.22 and 5.21, it appears clearly that spatial maps issued from our proposed technique reveal more consistent structures. This is to our opin-this point, we have computed the average pairwise difference between all groups estimated with OT-Unmixing and SLSQP techniques. More precisely, for each unmixing strategy ST i , i = {2, 3, 5, 6, 8, 9} where i refers to the strategy used (represented in table 5.2), we have computed the average difference

D i = 1 |C 2 5 | k, D(A i,k , A i, ) such that k < and D(A i,k , A i, )
is the distance between abundance map of group k and technique i and the abundance map of group and technique i. The number of pairwise combinations among the 5 groups |C 2 5 | = 5! 2!3! = 10. This measures the heterogeneity between the different groups and this is therefore an indicator of the ability of each technique to separate various kind of mixtures (the higher is this difference, the more diverse mixtures are separated).

Of course we would like to point out here that this criteria is not strong in the sense that even with reliable values, this does not state that the unmixing is correct. This only quantifies the heterogeneity between abundances and therefore this relies on the assumption that all abundance maps are spatially different (which is quite rational for 5 classes). These results are given in table 5.4 as well as a visual representation of the pairwise difference matrices related to 6 strategies (the results of those strategies over the ASD dictionary were excluded because of their poor performance) are illustrated in Figure 5.23. As one can observe, this difference is considerably higher for OT-related strategies (ST-02 and ST-03) comparing to the SLSQP-related ones (ST-08 and ST-09). This simple quantitative interpretation strengthens the fact that our approach seems to have a better accuracy in unmixing urban patterns. There is no doubt however, that a more complete comparison with ground truth data could validate this. This is the scope of a current work.

Let us now turn to the experiments using the complete SAGA+ technique (endmember extraction and abundance estimation).

Endmember Extraction and Abundance Estimation with SAGA+

In the last phase of tests, we used the sparse projection part of SAGA+ in order to produce abundance maps using the same three dictionaries.

As for this test we set the RBF kernel parameter σ = 2 × 10 -5 and the sparsity λ = 5. This choice of λ is made to make the 5-sparse maps as the result of SAGA+ comparable to the 5-groupings of the OT-Unmixing and the other method.

The 5-sparse abundance maps based on the gready selector and sparse projector algorithm are represented in Figure 5.24. Here again, the same kind of conclusion than using the OT-unmixing procedure (Figure 5.21) holds. Even if abundances are inconsistent using ASD endmembers, the spatial groups extracted using Expert and ESiVM dictionaries are consistent, visually even less noisy than using OTbased results. We also have computed the average of pairwise distances between endmembers. The values are depicted in table 5.4 and Figure 5.23 respectively. We can also conclude that the various materials of the scene appear to be more consistent here. Moreover using SAGA+, performances in terms of computational cost are lower since the same calculation took place in about 15 minutes. Ongoing works are focusing on ground truth validations to perform more complete quantitative comparisons.

Conclusion of the chapter

In this last part of this PhD thesis we provided some preliminary results on a specific challenges of spectral unmixing for surface temperature estimation in an urban context. Preliminary results were presented and some comparisons were with a state-of-the-art approach seems to promote the use of the two algorithms developed through this work. Indeed stronger validations have now to be done. This is the scope of a current campaign with ground truth data.

Conclusion and future work

Spectral unmixing has gained a considerable amount of attention during the recent years within the remote sensing community. As the aftermath of continuously improving spectral resolution of the hyperspectral sensors, more accurate measurements of materials have become feasible.

The major issue in hyperspectral unmixing is the way existing endmembers have mixed at the onset of pixel registration. The formation process of endmembers within a pixel is formally called the mixing model. Depending on the context, different mixing models are being assumed while studying remotely sensed data. Mixing models can be split in three major categories, namely linear, non-linear and intimate mixing models. Having a correct assumption of the mixing model, while analyzing hyperspectral data has a great impact on the accuracy of the analysis. Another important issue related to unmixing is the issue of endmember variability where spectral signature of the same material shows subtle differences from one sample to the other. Yet another common issue within the field of spectral analysis is the presence of anomalies which can either though as unwanted artifacts that spoil the true value of the pixels or on the opposite, to be seen as highly important but spatially scarce information hidden within the scene. In this thesis, the three aforementioned issues were mainly brought into light and by introducing two original algorithms, defined within different mathematical frameworks, solutions to these open problems has provided.

• The first contribution of this PhD thesis, brings all the three mentioned problems of mixture modeling, endmember variability and anomaly detection into focus and using the applications of kernel theory proposes a new unsupervised algorithm (SAGA+) for representation of the non-linear manifold embedding the data while through a simultaneous anomaly detection procedure makes sure that the representation of the manifold hall is not being distorted at the presence of anomalies. In fact SAGA+ can be seen as a kernel-based non-negative matrix factorization (NMF) which is coupled with an anomaly detection. It operates by coding data w.r.t the local neighbors with non-linear weights. For this the solution of NMF is divided in two steps of estimations of data points characterizing the boundaries of the manifold hull E embedding 123 data and then using this as the endmember dictionary to estimate the coefficients of each endmember contributing to the reconstruction of each pixel value of the original dataset. The locality is the aftermath of the simultaneous imposition of sparsity and convexity constraints. The set of abundance maps produced by this procedure therefore, is the convex surrogate of the original data which is achievable through a controllable level of sparsity. The matrix of abundances G, encodes the convex combinations and the columns of the endmember dictionary E characterizes a hull of dataset X in R D , as data points lying within are perfectly reconstructed. The Anomaly detection procedure is formulated in order to play the role of a supervisor through the first part of SAGA+. The idea consists of iteratively increasing the volume of the manifold hull in feature space while at the same time keeping track of the rate of change of the sum of projection errors (SOPE) in order to exclude anomalies from being taken as simplex vertexes, which is being formed based on the previously selected endmembers. In the second step of solving NMF, in order to account for the endmember variability among endmembers (as can be chosen large even to outnumber the dimensionality of the data), by manually imposing a sparsity value λ, each reconstructed point within G will be a convex combination of λ-selected atoms within the -endmember dictionary.

This is what makes this algorithm unique in addressing all the three mentioned problems within a unified framework.

• The second major contribution of this PhD focuses mainly on the issue of endmember variability and by exploiting the notion of overcomplete dictionary tries to address this problem. The main assumption is that a dictionary of endmembers, possibly larger than the effective number of pure materials is available to deal with the internal heterogeneity inside a particular material (which should correspond to only one endmember in classical approaches). This supervised algorithm (OT-unmixing) which is based on the optimal transport theory is comparable to the second step of SAGA+, as it solves an inversion problem and calculates the sparse representation of the original datum/pixels through generation of the abundance maps. For this OT-Unmixing tries to minimize a transition plan between the datum (represented as a probability distribution in the spectral domain) and the known atoms of the dictionary (also in the form of probability distributions). An original metric, designed between probability density functions is considered to compare spectra. This metric (Wasserstein distance), is also based on optimal transport and can be designed to be mostly sensitive to shifts in the frequency domains. Although this problem is solvable through linear programming with equality constraints, considering that the computational cost scales quadratically with the dimension of the spectra, a regularization term is used to help increasing the speed of the calculation. This regularization term which controls the smoothness of the transition map is based on the entropy of the map. Having the discrete version of the OT with an Entropy regularization term, it is possible to establish a link between OT and Bregman projection which by itself admits a simple closed-form solution. Therefore for each datum/pixel signature, OT-Unmixing uses an iterative Bregman projection to project the signature onto the atoms of the dictionary and calculating the Wasserstein distance in each iteration. It worths mentioning that for producing sparse abundance maps the a priori information regarding the λ groups representing endmember variabilities, also the respective weights are being provided through a separate matrix and a weight vector.

The mentioned unmixing algorithms have been used on synthesized and real datasets in different settings in order to be evaluated and compared with existing state-of-the-art procedures. In particular SAGA+ has shown its robustness to highly corrupted by anomalies planetary dataset, letting the expert to identify small patches of otherwise difficult to detect lithologies on the surface of 4-Vesta asteroid. Besides, to be the third major contribution of this thesis, both algorithms were deployed in the context of CATUT project where, SAGA+ and OT-unmixing were used in providing abundance maps of an urban scene. The utilization and comparing different vector-based spectral distance measures was among the subtasks of this thesis. In particular to compare the endmember estimation performance of SAGA+ with other endmember extraction techniques. In a current work a more complete evaluation and comparison of distance measures is considered by including other correlation weighted measures along the well-known vector-based measures. Also from the computational time point of view the algorithms are at least as good as their state-of-the-art counterparts. In the case of SAGA+ having an overall linear computational complexity, the algorithm is faster than many non-linear spectral unmixing procedures. This is an asset while processing large volumes of data like planetary datasets. In the case of OT-Unmixing depending on the maximum iterations that it takes for the minimization to converge, the results can be up to twice more expensive.

As for the future possible continuation of the current work one can mention the potential capacity of the current models to be used on multi-angular and multitemporal hyperspectral datasets. The new multi-angular sensors have been contributing to many astronomical and astrophysical projects (CRISM, OMEGA etc.) providing additional (angular) measurement dimension, therefore a better separation of the signals coming from the atmosphere and the surface become possible. This allows for more accurate characterization of the planetary materials. These hyperspectral multi-angular (HMA) observations, available in large quantities, offer new perspectives but require extensive methodological research and development because of their very nature and the complexity of natural scenes. SAGA+ seems quite a proper method to be extended for the exploitation of this angular dimension. It can simply expanded by dedicating several kernels to simultaneously unmix the pixels contributing to the same spatial position but representing different angular/temporal perspectives. This will requires us also to devise a fusion mechanism after the value of the current pixel being projected to the feature space through L kernels (where L represents the number of angular datasets). By using a segmentation strategy (Metropolis-within-Gibbs sampler for example) we can partition the data space into convex regions where separate units represented by these regions are being modeled with not all the dictionary atoms but a limited (hopefully not more than a handful) number of estimated endmembers. Also by considering the possibility of having different mixing paradigms underlying different areas of a scene, it is possible to investigate to what degree segmentation of the scene into separated units and unmixing each unit using a dedicated kernel can help increasing the overall accuracy of the unmixing procedure. In such scenarios one can think of dedicating the result of separate endmembers estimation per segmentation. Besides that by availability of datasets covering vast homogeneous areas of land surface, it will be effective to exploit the potential of domain adaptation techniques for generalizing the unmixing result of a well-studied, but limited-in-spatial-size area to a large scene while improving the computation time. Also the natures of the presented algorithms, are such that make the parallelization of the procedures feasible, thus suitable for exploiting the computation power of parallel-processing units and GPGPUs. As for the OT-unmixing it is possible to think about exploitation of other regularization techniques instead of the entropy-based regularizer presented in this work. Among other regularization methods one can mention LASSO.

Part V Appendices 7.2 SOPE maps In each plot of Figure 7.2, all points of a meshgrid of size (26 × 26) are taken as the next prototype and the SOPE value is calculated. The only difference is the number of Endmembers which have already selected to construct the representative manifold. As one can see by selecting endmember candidates from the borders of the nominal data, in each iteration the overall SOPE value drops, until the time when adding an extra simplex does not improve the SOPE significantly. In this particular case we can see that by selecting the fifth endmember the overal SOPE does not change significantly which means that the 3-simplex is almost as good as the 4-simplex in representing this dataset.

Parameter Setting

Anomalies (Υ ) 20

Dimension (D) 383

In the table 7.1 the notations represented different model parameters are gathered together with their range of values which they used for synthetic tests.

Execution Time

In table 7.2 the execution times of the tested algorithms (SAGA+, G-SiVM, Nfinder and VCA) using different endmembers and instance numbers N have been compared . From this comparison we observed that VCA outperformed any other algorithm while (G-SiVM) has a systematic larger computation time which makes it unusable for real applications. This large computation time is explained by its inherent complexity of the computation of Geodesic distances which scales quadratically with the size of data. As for SAGA+ and N-finder, only minor differences can iii Contributions to Hyperspectral Unmixing Sina Nakhostin 2017 be observed. As SAGA+ is a non-linear procedure able to capture more subtleties in the variations of signatures, this low computational time makes it a very appealing solution.

Orthogonal

  Projection Divergence (OPD) . . . . . . . . . . . . Spectral Information Divergence (SID) . . . . . . . . . . . . . Mixed SID -SAM . . . . . . . . . . . . . . . . . . . . . . . . . 2.2.2 Correlation-weighted measures . . . . . . . . . . . . . . . . . . HS measures weighted by a priori correlation . . . . . . . . . HS measures weighted by a posteriori correlation . . . . . . . 2.3 Inversion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4 Anomaly detection . . . . . . . . . . . . . . . . . . . . . . . . . . . . II Kernel methods for Unmixing 3 Joint Anomaly detection and spectral Unmixing 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1.1 Summary, goals of the chapter . . . . . . . . . . . . . . . . . . 3.2 SAGA+ : simplex volume maximization and anomaly detection . . .

3

 3

1. 6

 6 Left: Original painting. Right: Endmember map of the materials used for painting. The red corresponds with egg yolk binder, the blue is consistent with a glue binder, and the green maps the areas of Azurite in a glue binder. . . . . . . . . . . . . . . . . . . . . . . . . 1.7 Using an iterative linear spectral unmixing model the global OMEGA dataset has been analyzed. The goal was to retrieve global atlas of the main minerals contributing to the signal recorded by OMEGA on the very top surface of Mars. Top : false color composite emphasizing water ice (in blue). Bottom : false color composite with Red=Olivine, Green= Low-Calcium pyroxene, Blue = High-Calcium pyroxene. (Courtesy of ESA, OMEGA, LPGNantes) . . . . . . . . . . 1.8 Left : MUSE instrument (property of ESO) installed at Nasmyth focus of UT4. Right : An abstract illustration of image registration by MUSE across the spectral range. . . . . . . . . . . . . . . . . . . . 1.9 Hyperspectral image (Left) vs. Multipectral image (Right). Hyperspectral Images provide an unprecedented amount of spectral information which requires new approaches of image analysis in order to exploit this amount of information. (Illustration courtesy of [1]) . . . 1.10 Each pixel can be the mixing result of different endmembers. In this abstract scenario three different materials are being mixed according to the linear mixing model (LMM). . . . . . . . . . . . . . . . . . . . 1.11 Representation of the northern hemisphere of 4-Vesta asteroid together with some nominal spectral signatures and one signature which is corrupted by anomalies. . . . . . . . . . . . . . . . . . . . . . . . . 1.12 DESIREX 2008 Madrid classification map. . . . . . . . . . . . . . . .

3. 4 3 . 6 3 . 7 3 . 8 3 . 9

 436373839 Illustration of synthetic data. (a): Some endmembers (top) and the anomaly spectrum (bottom) used to generate the datasets. (b): A 2D Slice of two components of a synthetic dataset of 1010 samples generated based on LMM from the 3 nominal endmembers and 10 samples issued from the single anomaly signature. The geometrical location of the anomalies in the 2D plan makes them visually separable. (c): Similar visualization but mixtures are generated on the basis of BMM. (d): Similar visualization but mixtures are generated on the basis of HCM. . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5 4-VESTA mosaic composed of (842 × 327) pixels, with 383 spectral bands covering the range 0.55µm to 2.47µm. Examples of spectra (located on the mosaic with letters) are reported with (a) a normal spectrum, (b) discontinuity at the two detectors boundary and (c) a spectrum at the boundary of a shadow region. The gray wavelength range is where inconsistency are usually observed. . . . . . . . . . . . Endmember estimation Averages of SAM values for SAGA, G-SiVM, VCA and N-Finder in function of the number of endmembers. Each plot represents the results for three mixing models. . . . . . . . Endmember estimation Averages of SAM values for SAGA+, G-SiVM, VCA and N-Finder in function of the number of endmembers. Data contains anomalies generated based on 3 random signatures. Each plot represents the results for three mixing models. . . . . . . . Results in presence of anomalies using SAGA+. (a): Representation of the dataset; (b) original endmembers and anomalies (c) estimated endmembers and anomalies with SAGA+ . . . . . . . . . Evaluation of SAGA+: Kappa values as the function of Concentration of the anomalies. The red lines represent the threshold where the Kappa value starts to get over 0.8 which is equal to the concentration level α = 15. . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.10 Effect of Anomaly concentration: In each plot a 2D slice of a dataset generated based on the same endmember signatures and anomaly signatures is illustrated. For construction of anomalies different levels of concentration is set : α = [5, 15, 50]. . . . . . . . . . . 3.11 Results in presence of anomalies using SAGA+. (a): original signatures used for generation of dataset; (b): 2D Slice of dataset together with the estimated prototypes (c): signatures of estimated endmembers and anomalies with SAGA+ (d): the 4 weighted averaged endmembers. . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3. 14

 14 RGB composition with R=P9+P5, G=P4+P8+P0, and B=P4+P8+P0. Coloring shows the localization of outcrops enriched in olivine (appearing in red) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4. 1 5 . 2 104 5. 7

 1521047 Mosaic of northern hemisphere of Vesta-4 asteroid. The window in red (61 × 61) is used as the test bed. . . . . . . . . . . . . . . . . . . 5.1 AHS flight lines for urbane surface measurements . . . . . . . . . . . Left: Platform C-212-200 property of INTA (PATERNINA). Right: Airborn hyperspectral scanning sensor. . . . . . . . . . . . . . . . . . 5.3 Illustration of the spectral range and spectral resolution of each band for the first three ports of the AHS DESIREX dataset. One observes that the bands of DESIREX dataset are measured using different portions of the electromagnetic field and with different spectral resolutions.102 5.4 DESIREX hyperspectral images acquired through one north-south flight over Madrid in 2008. (The image is rotated by 90 • ) . . . . . . . 5.5 ASD FieldSpec3 portable, general-purpose Spectroradiometer. . . . . 5.6 Man-made materials considered for endmember dictionary construction.Spectral signatures of Endmembers which are registered through an in-situ procedure using the ASD equipment. . . . . . . . . . . . . . . 5.8 Illustration of the "wavelength range" and "wavelength sampling rate" of the in-lab measured ASD endmembers across the electomagnetic field. One can see that the bands of these endmembers have uniformly been sampled across the range .35 -2.5(µm), with the sampling rate of 1(nm) which makes them incomparable to the spectral characteristic of the Desirex dataset. . . . . . . . . . . . . . . . . . . . . . . . 5.9 Spectral signatures of a sample of 1000 pixels randomly selected from Desirex dataset versus the original ASD endmembers. . . . . . . . . . 5.10 Spectral signatures of a sample of 1000 pixels randomly selected from Desirex dataset are plotted versus the clipped ASD endmembers. . . . 5.11 Spectral signatures of a sample of 1000 pixels randomly selected from Desirex dataset versus the downsampled ASD endmembers. . . . . . .

Figure 1 . 1 :

 11 Figure 1.1: Left: Active sensor Right: Passive sensor (Courtesy of NASA)

Figure 1 . 2 :

 12 Figure 1.2: Electromagnetic wavelength: depending on their presumed utility, different spectrometers are being designed to capture specific portions of the electromagnetic field.In many applications therefore, hyperspectral data combines more than one captor.

Figure 1 . 3 :

 13 Figure 1.3: Chronological process of water shrinkage of lake Urmia. With a new composite of multi-spectral high resolution satellite observations, it is shown that the area of the Iranian lake has decreased by around 88% in the past decades, far more than previously reported ∼ 25% to 50%. (Courtesy of Kaveh Madani)

Figure 1 . 4 :

 14 Figure 1.4: Paris area mean temperature map in 2003. (Courtesy of VITO, PLANETEK)

Figure 1 . 5 :

 15 Figure 1.5: Generated maps based on two different hyperspectral images captured from the same subject with similar tissues depicted with similar colors.

  and Mapping • Examination of Underdrawings • Analysis of 3-D Objects • Recognition of Organic Materials • Forgery Detection

Figure 1 . 6 :

 16 Figure 1.6: Left: Original painting. Right: Endmember map of the materials used for painting. The red corresponds with egg yolk binder, the blue is consistent with a glue binder, and the green maps the areas of Azurite in a glue binder.

Figure 1 . 7 :

 17 Figure 1.7: Using an iterative linear spectral unmixing model the global OMEGA dataset has been analyzed. The goal was to retrieve global atlas of the main minerals contributing to the signal recorded by OMEGA on the very top surface of Mars. Top : false color composite emphasizing water ice (in blue). Bottom : false color composite with Red=Olivine, Green= Low-Calcium pyroxene, Blue = High-Calcium pyroxene. (Courtesy of ESA, OMEGA, LPGNantes)

Figure 1 . 8 :

 18 Figure 1.8: Left : MUSE instrument (property of ESO) installed at Nasmyth focus of UT4. Right : An abstract illustration of image registration by MUSE across the spectral range.

Figure 1 . 9 :

 19 Figure 1.9: Hyperspectral image (Left) vs. Multipectral image (Right). Hyperspectral Images provide an unprecedented amount of spectral information which requires new approaches of image analysis in order to exploit this amount of information. (Illustration courtesy of [1])

Figure 1 . 10 :

 110 Figure 1.10: Each pixel can be the mixing result of different endmembers. In this abstract scenario three different materials are being mixed according to the linear mixing model (LMM).

Figure 1 . 11 :

 111 Figure 1.11: Representation of the northern hemisphere of 4-Vesta asteroid together with some nominal spectral signatures and one signature which is corrupted by anomalies.

Figure 1 . 12 :

 112 Figure 1.12: DESIREX 2008 Madrid classification map.

  -the-art Chapter Summary Satellite image from a central desert of Iran (Dasht-e-Kavir ), representing a large mixture of different minerals with different levels of moisture. These areas are salt lakes and the colors show fluctuations in the amount of water present over time. The image is acquired by Landsat-1. (Courtesy of ESA)

Figure 2 . 1 :

 21 Figure 2.1:In the ideal case, where there is no effect of noise, also with the assumption of pure endmembers being present in the data, the abundance coefficients as the result of inversion part of linear unmixing will lay within a convex polygon (n-1 simplex) whose vertexes are consist of the endmembers. Illustration courtesy of Dobigeon et al.

Figure 2 . 2 :

 22 Figure 2.2: From left to right: linear mixing model(LMM), intimate mixing model (IMM) and bilinear mixing model (BMM). Courtesy of [2].

Figure 3 . 1 :

 31 Figure 3.1: Overcomplete endmember dictionary for non-linear unmixing

Figure 3 . 2 :

 32 Figure 3.2: Illustration of difficulties in presence of outliers. First Row: endmember extraction using SAGA without anomalies: in this situation extracted endmembers are a reliable representation of the manifold ; Second Row: endmember extraction using SAGA in presence of anomalies: in this situation extracted endmembers (in red) do not correspond to real borders of the data manifold ; Third Row: endmember extraction and anomaly detection In presence of anomalies using SAGA+: here extracted endmembers (in red) by are in accordance with borders of the manifold embedding data while the anomalies (in cyan) are properly selected.

Figure 3 . 3 :

 33 Figure 3.3: Evolution of SOPE (Sum Of Projection Errors) when all data in green are represented through a single endmember uniformly represented on a meshgrid. One observes a decrease in SOPE when candidates get closer to the nominal distribution.

Algorithm 1 SAGA+ 1 :Step 1 :

 111 Input: Data matrix X 2: Parameters: • Number of required endmembers • Threshold τ required for SOPE (see section 3.2.2) • Sparsity level λ for abundance matrix G • Kernel bandwidth σ 3: Output: • Endmember matrix E • Abundance matrix G • Anomaly matrix A Endmember Estimation (EE) and Anomaly Detection (AD): 4: -Select a random point x i . 5: -Compute distances w.r.t x i in RKHS (see eq. (3.4)-(3.6)). 6: -Identify the first endmember and first outliers using Exact SiVM and AD algorithm (see description in Figure 2) 7: for (p = 2, . . . , ) do -From current simplex with p-1 endmembers, compute volume changes with all instances (see eq. (3.3)).

Figure 3 .Figure 3 . 4 :

 334 Figure 3.4: Illustration of synthetic data. (a): Some endmembers (top) and the anomaly spectrum (bottom) used to generate the datasets. (b): A 2D Slice of two components of a synthetic dataset of 1010 samples generated based on LMM from the 3 nominal endmembers and 10 samples issued from the single anomaly signature. The geometrical location of the anomalies in the 2D plan makes them visually separable. (c): Similar visualization but mixtures are generated on the basis of BMM. (d): Similar visualization but mixtures are generated on the basis of HCM.

Figure 3 . 5 : 4 -

 354 Figure 3.5: 4-VESTA mosaic composed of (842 × 327) pixels, with 383 spectral bands covering the range 0.55µm to 2.47µm. Examples of spectra (located on the mosaic with letters) are reported with (a) a normal spectrum, (b) discontinuity at the two detectors boundary and (c) a spectrum at the boundary of a shadow region. The gray wavelength range is where inconsistency are usually observed.

Figure 3 . 6 :

 36 Figure 3.6: Endmember estimation Averages of SAM values for SAGA, G-SiVM, VCA and N-Finder in function of the number of endmembers. Each plot represents the results for three mixing models.

Figure 3 . 7 :

 37 Figure 3.7: Endmember estimation Averages of SAM values for SAGA+, G-SiVM, VCA and N-Finder in function of the number of endmembers. Data contains anomalies generated based on 3 random signatures. Each plot represents the results for three mixing models.

Figure 3 . 8 :

 38 Figure 3.8: Results in presence of anomalies using SAGA+. (a): Representation of the dataset; (b) original endmembers and anomalies (c) estimated endmembers and anomalies with SAGA+

Figure 3 . 9 :

 39 Figure 3.9: Evaluation of SAGA+: Kappa values as the function of Concentration of the anomalies. The red lines represent the threshold where the Kappa value starts to get over 0.8 which is equal to the concentration level α = 15.

Figure 3 . 10 :

 310 Figure 3.10: Effect of Anomaly concentration: In each plot a 2D slice of a dataset generated based on the same endmember signatures and anomaly signatures is illustrated. For construction of anomalies different levels of concentration is set : α = [5, 15, 50].

Figure 3 . 11 :

 311 Figure 3.11: Results in presence of anomalies using SAGA+. (a): original signatures used for generation of dataset; (b): 2D Slice of dataset together with the estimated prototypes (c): signatures of estimated endmembers and anomalies with SAGA+ (d): the 4 weighted averaged endmembers.

Figure 3 . 12 :

 312 Figure 3.12: Results on real 4-Vesta asteroid. Top: Vesta-4 mosaic with the spatial position of the extracted endmembers. The red points are detected anomalies and the 10 green points with labels are selected endmember prototypes. Bottom: The respective estimated endmember signatures.

Figure 3 .

 3 Figure 3.13: Top: laboratory spectra of mixtures involving olivine and orthopyroxene (spectra from the RELAB library with numbers referring respectively to the olivine and orthopyroxene contents (% olivine / % orthopyroxene). Bottom: comparison with the three extracted endmembers P0, P5 and P9. The shift in absorption center at 2µm is due to an increase in Ca and/or Fe in the pyroxene, compared to the pure orthopyroxene composition.

Figure 3 . 14 :Figure 3 . 15 :

 314315 Figure 3.14: RGB composition with R=P9+P5, G=P4+P8+P0, and B=P4+P8+P0. Coloring shows the localization of outcrops enriched in olivine (appearing in red)

Figure 4 . 1 :

 41 Figure 4.1: Mosaic of northern hemisphere of Vesta-4 asteroid. The window in red (61 × 61) is used as the test bed.

  (a) Estimated Abundance Maps with different priors α 0 . (b) Estimated Abundance Maps with Constrained Least Squares.

5 CATUT:

 5 Spectral Unmixing for Urban surface Temperature 5.1 DESIREX 2008 project Land surface temperature (LST) estimation is one of the goals of the "Dual-use European Security IR Experiment 2008" (DESIREX 2008) campaign in Madrid Spain, carried out by the European space agency (ESA) through several acquisition missions between June 23rd and 6th of July. Spectrally, geometrically and radiometrically representative datasets were generated and estimating Urban Heat Island (UHI) and Urban Thermography (UT) modeling are the two major objectives. Because of the limited spatial resolution of acquired hyperspectral data (from 10 to 30 square meters in order) estimating LST needs two steps: 1) Spectral unmixing of the pixels to recover abundances of each pure material and 2) Calculation of the area temperature based on the endmember abundances and associated emissivity models for each material. Data were acquired with an airborne hyperspectral scanner (AHS) covering the spectral regions of 445 -12450nm with varying spectral resolutions (across the spectral range) For more details regarding airborne and satellite data acquisition, one can refer to the DESIREX 2008 final report a . a https://earth.esa.int/c/document_library/get_file?folderId=21020&name= DLFE-905.pdf 5.1.1 Description of the urban area DESIREX measurements were registered over the city of Madrid in two directions north-south and west-east.

Figure 5 .

 5 [START_REF] Lu | Medical hyperspectral imaging: a review[END_REF] shows the location of these lines in green crossed mosaics.

Figure 5 . 1 :

 51 Figure 5.1: AHS flight lines for urbane surface measurements

Figure 5 . 2 : 4 FWHM

 524 Figure 5.2: Left: Platform C-212-200 property of INTA (PATERNINA). Right: Airborn hyperspectral scanning sensor.

Figure 5 . 3 :

 53 Figure 5.3: Illustration of the spectral range and spectral resolution of each band for the first three ports of the AHS DESIREX dataset. One observes that the bands of DESIREX dataset are measured using different portions of the electromagnetic field and with different spectral resolutions.

Figure 5 . 4 :

 54 Figure 5.4: DESIREX hyperspectral images acquired through one north-south flight over Madrid in 2008. (The image is rotated by 90 • )

Figure 5 . 5 :

 55 Figure 5.5: ASD FieldSpec3 portable, general-purpose Spectroradiometer.

Figure 5 . 6 :

 56 Figure 5.6: Man-made materials considered for endmember dictionary construction.

Figure 5 . 7 :

 57 Figure 5.7: Spectral signatures of Endmembers which are registered through an in-situ procedure using the ASD equipment.

Figure 5 . 8 :

 58 Figure 5.8: Illustration of the "wavelength range" and "wavelength sampling rate" of the in-lab measured ASD endmembers across the electomagnetic field. One can see that the bands of these endmembers have uniformly been sampled across the range .35 -2.5(µm), with the sampling rate of 1(nm) which makes them incomparable to the spectral characteristic of the Desirex dataset.

Figure 5 . 9 :

 59 Figure 5.9: Spectral signatures of a sample of 1000 pixels randomly selected from Desirex dataset versus the original ASD endmembers.

Figure 5 . 10 :

 510 Figure 5.10: Spectral signatures of a sample of 1000 pixels randomly selected from Desirex dataset are plotted versus the clipped ASD endmembers.

Figure 5 . 11 :

 511 Figure 5.11: Spectral signatures of a sample of 1000 pixels randomly selected from Desirex dataset versus the downsampled ASD endmembers.

Figure 5

 5 Figure5.12: Dictionary of the manually extracted endmembers. Based on expert observations, pixel vectors related to highly distinct materials within the scene have been extracted to be used as a dictionary of endmembers. In order to account for material variability, for each material two representatives endmembers were selected in order to build-up a dictionary of 24 signatures.

Figure 5 . 14 :

 514 Figure 5.14: Dictionary of 24 randomly selected signatures from 601 × 601 window of Desirex dataset.

Figure 5 . 15 :

 515 Figure 5.15: Random and ESiVM dictionaries compared to the Expert dictionary. The result of ESiVM shows more fidelity to the Expert dictionary in capturing more varieties among the endmember signatures.

Figure 5 . 16 :

 516 Figure 5.16: Window of 601 × 601 (pixels) extracted from the main mosaic on which all tests are performed.

Figure 5 .

 5 Figure 5.17: Left: Result of Elbow test on a series of k-means clusterings with increasing number of clusters performed on the endmember dictionary of manually selected atoms from the AHS Desirex dataset. Right: Result of Elbow test on a series of k-means clusterings with increasing number of clusters performed on the endmember dictionary generated by performing E-SiVM procedure on the AHS Desirex dataset.

Figure 5 . 18 :

 518 Figure 5.18: In order to set cost matrix C0 in alignment to the characteristics of the dataset, we build C0(63 × 63) as the square euclidean distance over the spectral values of the AHS Desirex dataset.

Figure 5 . 19 :

 519 Figure 5.19: C1 cost function for three dictionaries. C1 provides a priori information regarding the groups that different signatures of the endmember dictionary belong to.

Figure 5 .

 5 Figure5.19 represents these three cost functions, used as inputs of our unmixing procedure. By comparing the results of clustering on the two endmember dictionaries (Expert and ESiVM) one can see that both are grouped in a more or less the same ratio. The number of endmembers belonging to each group for each dictionary is given in table 5.3.

Figure 5 . 20 :

 520 Figure 5.20: Clustering of two dictionaries. Groups which considered to share similar semantic relationships in two dictionaries are plotted in the same color.

Figure 5 . 21 :

 521 Figure 5.21: Abundance maps produced by OT-unmixing and according to the first 3 strategies presented in the first column of table 5.2.

Figure 5 . 23 :

 523 Figure 5.23: Map distances calculated between pairs of abundance maps for 6 different strategies.

Figure 5 . 24 :

 524 Figure 5.24: Abundance maps produced by SAGA+ algorithm and according to the strategies 4,5 and 6 presented in the middle column of table 5.2.

Figure 7 . 2 :

 72 Figure 7.2: Evolution of SOPE map along the endmember extraction process. Each plot shows the map of SOPE values w.r.t the "next" endmember candidate. In each plot, we take all the cells of 26 × 26 grid one by one as the next endmember candidate and calculate SOPE. In each plot the only difference is the dimensionality of the simplex that we use to calculate the SOPE values. It is observable that the points laying outside the boundaries of the dataset always produce very high values of SOPE (almost close to 1). By setting a correct threshold for SOPE change then SAGA+ distinguishes the good candidates (points on the boundaries of the dataset manifold) from the anomalies (points outside the boundaries of the dataset manifold)

7. 5

 5 Abundance Maps for SAGA Abundance maps of 4-Vesta based on 10 estimated endmembers are presented in Figure 7.3. The selected signatures by SAGA shows that without an anomaly detection procedure this method can select corrupted signatures instead of real signatures. The produced maps therefore can be effected by the poor result of endmembers estimation.

Figure 7 . 3 :

 73 Figure 7.3: Abundance maps and respective spectral signatures extracted by SAGA (no anomaly detection). Without an anomaly detection procedure, one observes that signatures with sharp peaks are among selected endmember signatures.

  .3.2 Results on synthetic data . . . . . . . . . . . . . . . . . . . . 3.3.3 On the application of Overcomplete Dictionary . . . . . . . . 3.3.4 Results on real data: 4-Vesta Asteroid . . . . . . . . . . . . . .3.4 Endmember Extraction and Abundance Estimation with SAGA+120 5.4 Conclusion of the chapter . . . . . . . . . . . . . . . . . . . . . . . . 122 Spectral Signatures . . . . . . . . . . . . . . . . . . . . . . . . . . . . i 7.2 SOPE maps . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii 7.3 Parameter Setting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii 7.4 Execution Time . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii 7.5 Abundance Maps for SAGA . . . . . . . . . . . . . . . . . . . . . . . v Left: Active sensor Right: Passive sensor (Courtesy of NASA) . . . 1.2 Electromagnetic wavelength: depending on their presumed utility, different spectrometers are being designed to capture specific portions of the electromagnetic field. In many applications therefore, hyperspectral data combines more than one captor. . . . . . . . . . . 1.3 Chronological process of water shrinkage of lake Urmia. With a new composite of multi-spectral high resolution satellite observations, it is shown that the area of the Iranian lake has decreased by around 88% in the past decades, far more than previously reported ∼ 25% to 50%. (Courtesy of Kaveh Madani) . . . . . . . . . . . . . . . . . . 1.4 Paris area mean temperature map in 2003. (Courtesy of VITO, PLANETEK) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.5 Generated maps based on two different hyperspectral images captured from the same subject with similar tissues depicted with similar colors.

	6 Conclusion and future work V Appendices 7 Appendix : SAGA+ 7.1 List of Figures 123 127 i 1.1

.2.1 SAGA principles: Exact Simplex Volume Maximization for endmember extraction . . . . . . . . . . . . . . . . . . . . . . 3.2.2 Anomaly Detection . . . . . . . . . . . . . . . . . . . . . . . . 3.2.3 Sparse Projections onto the Simplex . . . . . . . . . . . . . . 3.2.4 Overall process . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 Evaluation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3.1 Experimental setups . . . . . . . . . . . . . . . . . . . . . . . Evaluation criteria . . . . . . . . . . . . . . . . . . . . . . . . Synthetic data . . . . . . . . . . . . . . . . . . . . . . . . . . . Real data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Parameters setting for SAGA and SAGA+ . . . . . . . . . . . 33.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

5

  . . . . . . . . . . . . . . . . . . . 5.14 Dictionary of 24 randomly selected signatures from 601 × 601 window of Desirex dataset. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.15 Random and ESiVM dictionaries compared to the Expert dictionary. The result of ESiVM shows more fidelity to the Expert dictionary in capturing more varieties among the endmember signatures. . . . . . . 5.16 Window of 601 × 601 (pixels) extracted from the main mosaic on which all tests are performed. . . . . . . . . . . . . . . . . . . . . . . .21 Abundance maps produced by OT-unmixing and according to the first 3 strategies presented in the first column of table 5.2. . . . . . . 5.22 Abundance maps produced by SLSQP algorithm and according to the last 3 strategies presented in the last column of table 5.2. . . . . . 5.23 Map distances calculated between pairs of abundance maps for 6 different strategies. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.24 Abundance maps produced by SAGA+ algorithm and according to the strategies 4,5 and 6 presented in the middle column of table 5.2. . 7.1 Example of spectra used for generating the synthetic datasets issued from the 41 spectral signatures of different Lithologies RELAB. The signatures have re-sampled in order to be comparable to the signature characteristics of 4-Vesta dataset according to the wavelength range and sampling range. . . . . . . . . . . . . . . . . . . . . . . . . . . . i 7.2 Evolution of SOPE map along the endmember extraction process. Each plot shows the map of SOPE values w.r.t the "next" endmember candidate. In each plot, we take all the cells of 26×26 grid one by one as the next endmember candidate and calculate SOPE. In each plot the only difference is the dimensionality of the simplex that we use to calculate the SOPE values. It is observable that the points laying outside the boundaries of the dataset always produce very high values of SOPE (almost close to 1). By setting a correct threshold for SOPE change then SAGA+ distinguishes the good candidates (points on List of Tables 3.1 Kappa statistics for SAGA+ in endmember estimation for the three synthetic datasets . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 3.2 Average SAM errors between the original and estimated abundance coefficients. Best results are reported in boldfont . . . . . . . . . . . . 77 5.1 AHS spectral configuration for DESIREX 2008 surveys. Due to the capability of the visible and near-to-mid infrared range of electromagnetic field in representation of the ground compositions and discrimination among different compositional units, only the signatures of the first three ports (PORT 1, PORT 2A and PORT 2) are used for the current study. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

	Spectral Analysis
	1.1 Remote sensing and its applications

5.17 Left: Result of Elbow test on a series of k-means clusterings with increasing number of clusters performed on the endmember dictionary of manually selected atoms from the AHS Desirex dataset. Right: Result of Elbow test on a series of k-means clusterings with increasing number of clusters performed on the endmember dictionary generated by performing E-SiVM procedure on the AHS Desirex dataset. . . . 5.18 In order to set cost matrix C0 in alignment to the characteristics of the dataset, we build C0(63 × 63) as the square euclidean distance over the spectral values of the AHS Desirex dataset. . . . . . . . . . . 5.19 C1 cost function for three dictionaries. C1 provides a priori information regarding the groups that different signatures of the endmember dictionary belong to. . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.20 Clustering of two dictionaries. Groups which considered to share similar semantic relationships in two dictionaries are plotted in the same color. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5the boundaries of the dataset manifold) from the anomalies (points outside the boundaries of the dataset manifold) . . . . . . . . . . . . ii 7.3 Abundance maps and respective spectral signatures extracted by SAGA (no anomaly detection). Without an anomaly detection procedure, one observes that signatures with sharp peaks are among selected endmember signatures. . . . . . . . . . . . . . . . . . v 5.2 Different strategies introduced in section 5.2 which relates different methods of unmixing considered in this chapter to the available endmember dictionaries. . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.3 Number of endmembers grouped into the 5 different clusters based on Expert and ESiVM dictionaries. . . . . . . . . . . . . . . . . . . 5.4 Statistics of the pairwise distance calculation over the 5 abundance groups produced as the result of each of 6 different strategies. The pairwise distances of each strategy over which we calculated the average value in the last column of this table are illustrated as a separate upper triangular matrix in Figure 5.23. . . . . . . . . . . . . . . . . . 7.1 Parameter setting for synthetic data generation . . . . . . . . . . . . iii 7.2 Execution Time (s) of the different algorithms for various endmembers and elements in the dataset . . . . . . . . . . . . . . . . . . . . . . . iv 1

Table 3 . 1 :

 31 Kappa statistics for SAGA+ in endmember estimation for the three synthetic datasets

		3	5	7	9	11	13	15
		LMM 0.73 0.84 0.85 0.84 0.84 0.84 0.79
	κ	BMM 0.94 0.94 0.89 0.88 0.85 0.85 0.78
		IMM 0.96 0.96 0.95 0.95 0.94 0.94 0.91

Table 3 . 2 :

 32 Average SAM errors between the original and estimated abundance coefficients.

	Best results are reported in boldfont
		3	5	7	10	12	15
	VCA	0.237	0.541	0.756	0.912	0.981	1.105
	N-finder 0.155	0.451	0.713	0.918	0.970	1.108
	G-SiVM 0.333	0.592	0.951	1.084	1.121	1.301
	SAGA+ 0.096 0.123 0.311 0.549 0.711 0.892

Table 5 . 1 :

 51 AHS spectral configuration for DESIREX 2008 surveys. Due to the capability of the visible and near-to-mid infrared range of electromagnetic field in representation of the ground compositions and discrimination among different compositional units, only the signatures of the first three ports (PORT 1, PORT 2A and PORT 2) are used for the current study.

Table 5 . 2 :

 52 Different strategies introduced in section 5.2 which relates different methods of unmixing considered in this chapter to the available endmember dictionaries.

Table 5 . 3 :

 53 The number of endmembers belonging to each group for each dictionary is given in table 5.3. Number of endmembers grouped into the 5 different clusters based on Expert and ESiVM dictionaries.

		G1 G2 G3 G4 G5
	C1 matrix over Expert	8	8	5	1	2
	C1 matrix over ESiVM 4	6	3	7	4

Table 7 . 1 :

 71 Parameter setting for synthetic data generation

		Parameter Values
	Sample Size (N)	10 3 . . . 10 6
	Endmember ( )	3 . . . 15
	Sparsity (λ)	3-5
	Kernel bandwidth (σ)	0.15
	Threshold (τ )	0.04
	Concentration LMM (α)	1
	Concentration IMM (α)

Table 7 . 2 :

 72 Execution Time (s) of the different algorithms for various endmembers and elements in the dataset

				# Samples (N )
	SAGA+		10 3	10 4	10 5	10 6
		= 3	.022	.255	3.058	31.432
		= 5 0.039	0.586	4.312	42.316
		= 7 0.044	0.462	20.716	54.562
	# ( )	= 9 0.065	0.687	10.872	53.666
		= 11 0.096	0.727	12.276	75.676
		= 13 0.085	0.833	12.012	87.943
		= 15 0.097	0.960	9.871	98.948
				# Samples (N )
	G-SiVM		10 3	10 4	10 5	10 6
		= 3 2.382 235.815 ≈ 3 × 10 3	_
		= 5 2.405			_
		= 7 2.507			_
	# ( )	= 9 2.674			_
		= 11 2.958			_
		= 13 3.174			_
		= 15 4.480			_
				# Samples (N )
	N-finder		10 3	10 4	10 5	10 6
		= 3 0.072	0.528	4.641	47.315
		= 5 0.103	0.832	7.796	86.957
		= 7 0.135	1.135	10.651	369.269
	# ( )	= 9 0.160	1.417	13.831	147.192
		= 11 0.248	2.091	20.158	225.753
		= 13 0.320	2.501	24.849	276.560
		= 15 0.364	3.010	29.752	496.775
				# Samples (N )
	VCA		10 3	10 4	10 5	10 6
		= 3 0.012	0.028	0.306	2.414
		= 5 0.015	0.027	0.257	2.501
		= 7 0.009	0.028	0.253	3.015
	# ( )	= 9 0.010	0.032	0.260	2.819
		= 11 0.010	0.032	0.302	3.184
		= 13 0.010	0.036	0.281	2.820
		= 15 0.013	0.030	0.292	3.684

Contributions to Hyperspectral Unmixing Sina Nakhostin 2017

http://muse-vlt.eu/science/publications/

Unfortunately, only after submitting this draft we became aware of a very recent work[START_REF] Altmann | Robust linear spectral unmixing using anomaly detection[END_REF], using Bayesian framework for joint anomaly detection and spectral unmixing.

In order to illustrate more into details, various maps of SOPE fields along the endmember extraction process are represented in the supplemental materials provided in the appendix 7.

These signatures together with additional information regarding them can be seen in the supplement material provided in the appendix 7.

In astrophysics, accretion is the growth of particles into a massive object by gravitationally attracting more matter, typically gaseous matter and dust, in an accretion disc.

Interested reader may find a plot of detected anomaly signatures from this mosaic in the supplement material provided in the appendix 7.

Contributions to Hyperspectral Unmixing Sina Nakhostin 2017Contributions to Hyperspectral Unmixing Sina Nakhostin 2017

Aknowledgement

The Author would like to thank Conseil Régional de Bretagne and Université de Bretagne sud for their support of the project ARED-RB OPTIMAP. Also thanks goes to ANR ASTERIX project for their support of this work. Many thanks to the authors in [START_REF] Bioucas-Dias | An overview on hyperspectral unmixing: Geometrical, statistical, and sparse regression based approaches[END_REF] for the available package and codes provided along the survey. The same appreciation goes to the authors in [START_REF] Jie Chen | Nonlinear unmixing of hyperspectral data based on a linear-mixture/nonlinear-fluctuation model[END_REF][START_REF] Nascimento | Vertex component analysis: a fast algorithm to unmix hyperspectral data[END_REF][START_REF] Heylen | Non-linear spectral unmixing by geodesic simplex volume maximization. Selected Topics in Signal Processing[END_REF] for providing the codes of their algorithms, which made the comparisons possible within this work. At the end Special thanks and gratitude goes to Nicolas Courty and Thomas Corpetti for their availability and supervision of this thesis. The positive and friendly atmosphere within the Obelix team of IRISA laboratory made the experience of accomplishing this thesis a true pleasure. Particular thanks goes to Sébastien Lefèvre for making it possible. ion a reliable behavior. In addition, the different groups extracted in Figure 5.21 are also more heterogeneous using OT-unmixing, which tends to promote a better separability of our approach. This is probably due to the non-linear property of our approach that enables to deal with the complex mixing phenomena generated by urban structures. To justify Appendix : SAGA+

Spectral Signatures

A total number of 41 spectral signatures used in order to generate synthetic datasets for this paper is represented in Figure 7.1. The signatures are related mostly to the lithologies which are found on the surface of meteorites from 4-Vesta asteroid which have hit the earth surface. Based on their compositions, the signatures were roughly categorized into 7 groups, each representing specific families of magmatic rocks. The original signatures are registered based on 461 contiguous electromagnetic bands covering the range of 0.3 -2.6µm. The signatures went through a re-sampling process in order to be approximately in comparison to the wavelength range and sampling rate of the 4-Vesta dataset. Therefore after re sampling each signature has 383 bands covering the range 0.55 -2.46µm.