# DÉMÉLANGE ROBUSTE D'IMAGES HYPERSPECTRALES EN INTERPRÉTANT LES PÔLES DE MÉLANGE COMME DES DONNÉES DIRECTIONNELLES

L. Drumetz¹, J. Chanussot², A. Iwasaki³

1: IMT Atlantique, Lab-STICC, UBL, 29238 Brest, France

2: Univ. Grenoble Alpes, CNRS, Grenoble INP, GIPSA-lab, Grenoble, France

3: Research Center for Advanced Science and Technology, The University of Tokyo, Tokyo, Japon

## Résumé

Le démélange d'images hyperspectrales vise à extraire les spectres des matériaux purs de la scène observée (pôles de mélange ou endmembers), ainsi que leurs proportions dans chaque pixel de l'image (abondances). La plupart des algorithmes s'appuie sur des techniques classiques d'extraction d'endmembers qui peuvent échouer dans des scénarios difficiles, en présence d'ombres, de données aberrantes ou encore de matériaux très corrélés. Dans cet article, nous abordons ce problème en même temps que la variabilité des matériaux en considérant qu'un endmember est une direction dans l'espace ambiant plutôt qu'un seul point. Sous ce paradigme, nous proposons un algorithme pour fournir des spectres de référence robustes. Nous montrons le potentiel de l'algorithme proposé sur un jeu de données synthétiques en utilisant des spectres réels présentant de la variabilité intra-classe, et une image d'une scène naturelle complexe.

**Mots clés :** Démélange hyperspectral, variabilité des endmembers, Modèle de mélange linéaire étendu, variété oblique

## Abstract

Hyperspectral image unmixing aims at recovering the spectra of the pure materials of the observed scene (endmembers), as well as their proportions in each pixel of the image (abundances). Most algorithms rely on classical endmember extraction techniques which may fail in difficult scenarios, i.e. when shadows, outliers or correlated materials are present. In this paper, we address this issue along with endmember variability by considering that an endmember is a direction in the feature space, rather than a single point. Under this paradigm, we design an algorithm to provide robust reference endmembers. We show the potential of the proposed algorithm on a synthetic dataset using real world spectra with variability, and a real image of a complex natural scene.

**Keywords :** Hyperspectral unmixing, endmember variability, Extended Linear Mixing Model, oblique manifold

## 1. Introduction

La télédétection hyperspectrale permet une identification automatique et précise des matériaux présents dans la scène observée, grâce à la fine résolution spectrale des images hyperspectrales (Bioucas-Dias et al., 2013). Leur résolution spatiale est cependant limitée, et le champ de vision d'un seul pixel comprend souvent plusieurs matériaux d'intérêt. Le spectre observé est alors un mélange des contributions des différents matériaux présents à cet endroit. Le démélange spectral est un problème de séparation aveugle de sources dont le but est de récupérer les spectres des matériaux purs (appelés endmembers), et d'estimer leurs proportions relatives (appelées abundances) dans chaque pixel (Keshava et Mustard, 2002). Un modèle de mélange linéaire (MML) (Bioucas-Dias et al., 2012) est généralement adopté pour ce problème. Il suppose qu'un pixel d'indice *n* parmi les *N* pixels de l'image $x_n \in \mathbb{R}^L$, où *L* est le nombre de bandes spectrales utilisées, est décomposé en :

$$x_n = \sum_{p=1}^{P} a_{pn} s_p + e_n \quad (1)$$

où *P* est le nombre de matériaux considérés, $s_p \in \mathbb{R}^L$ est la signature de l'endmember *p*, censé représenter parfaitement le matériau correspondant, $a_{pn}$ est le coefficient d'abondance de ce matériau dans le pixel *n*, et $e_n$ est un bruit additif. Les coefficients d'abondance, étant des proportions, sont soumis à la contrainte de non-négativité des abundances (CNA) et à la contrainte de somme à un des abundances (CSA), c'est-à-dire que chaque pixel doit être entièrement expliqué par des contributions positives des différents matériaux. Géométriquement, les données se trouvent dans un simplexe dont les sommets sont les endmembers.

Le MML a été largement utilisé au cours des deux dernières décennies, mais au fil des années, certaines de ses limites ont été identifiées, à savoir les non-linéarités dans les processus de mélange et la variabilité des matériaux. La première limitation est sensible dans des scénarios complexes tels que les canopées d'arbres ou les matériaux particulaires (par exemple du sable) et nécessitent des modèles de mélange plus complexes tels que des modèles de mélange bilinéaires (Heylen et al., 2014; Bigeon et al., 2014). D'un autre côté, la variabilité des
endmembers correspond simplement au fait que chaque
matériau ne peut être complètement représenté par un
seul spectre et est toujours sujet à des variations intra-
classe (Somers et al., 2011). Celles-ci peuvent être dues
à des facteurs multiples, les principaux étant les change-
ments locaux dans les conditions d'illumination (à cause
des ombres ou de la topographie), et la variabilité intrin-
sèque des matériaux, correspondant aux changements
physico-chimiques locaux dans leur composition. Cette
variabilité peut être considérée dans les domaines spa-
tial (Drumetz et al., 2017) ou temporel (Henrot et al.,
2016) quand il s'agit de séquences d'images. Ici, nous
nous concentrons sur la variabilité au sein d'une seule
image. Deux classes de méthodes pour y remédier ont
été identifiées dans (Zare et Ho, 2014) : méthodes re-
posant sur des "bundles", où les endmembers sont re-
présentés par une collection de signatures, éventuelle-
ment extraites des données (Somers et al., 2012), et des
méthodes statistiques, où des distributions statistiques
sont assignées aux endmembers pour permettre leur es-
timation dans chaque pixel (Halimi et al., 2015). Récem-
ment, des modèles explicites de variabilité ont été propo-
sés pour expliquer directement les variations possibles
des endmembers (Thouvenin et al., 2016; Veganzones
et al., 2014). Généralement, ils partagent la formulation
de l'équation (1), à ceci près que les endmembers sont
maintenant indexés par les pixels. Un de ces modèles est
le modèle de mélange linéaire étendu (MMLE), proposé
dans (Drumetz et al., 2016). Il approxime le modèle de
transfert radiatif physique de Hapke (Hapke, 2012) en une
version exploitable en démélange aveugle. Dans (Dru-
metz, 2016), la variabilité des matériaux due à l'illumi-
nation est en effet raisonnablement expliquée par des
facteurs d'échelle positifs dépendant du pixel. Dans sa
version la plus simple, le modèle s'écrit :

$$
\mathbf{x}_n = \psi_n \sum_{p=1}^{P} a_{pn} \mathbf{s}_p + \mathbf{e}_n \quad (2)
$$

où $\psi_n$ est un facteur d'échelle qui tient compte des changements de luminosité dans les signatures spectrales des pixels. Dans ce modèle, $s_p$ n'est plus l'unique signature associée au matériau $p$, mais plutôt une référence fournissant la direction d'une droite joignant l'origine et ce point, sur laquelle se trouve chaque endmember. Cette formulation a conduit à un algorithme rapide pour estimer les paramètres de ce modèle par moindres carrés non négatifs, appelé Scaled Constrained Least Squares Unmixing (SCLSU), voir (Veganzones et al., 2014) pour plus de détails. Cependant, cet algorithme ne peut pas gérer d'autres types de variabilité que des facteurs d'échelle. Tous les phénomènes de variabilité intrinsèque, qui sont beaucoup plus difficiles à modéliser physiquement en raison de la diversité de leurs causes, ne sont pas pris en compte.

Une version plus complexe de l'algorithme MMLE (Drumetz et al., 2015, 2016) a été proposée ensuite. Cette version considère différents facteurs d'échelle pour chaque matériau, mais permet également à chaque endmember de s'éloigner d'une version mise à l'échelle de la

<figure><img src="image_1.png" /><figcaption>**FIGURE 1:** Interprétation géométrique du MMLE, tel que présentée dans (Drumetz et al., 2016). Les endmembers sont en plus autorisés à s'éloigner de versions mises à l'échelle des références via un a priori gaussien.</figcaption></figure>

référence, en termes de distance euclidienne. D'un point
de vue bayésien, cela revient à mettre un a priori gaus-
sien sur chaque endmember local, centré sur la version
mise à l'échelle de la référence. De cette façon, la variabi-
lité récupérée est plus complexe qu'une simple variation
du facteur d'échelle et peut expliquer dans une certaine
mesure les effets de la variabilité intrinsèque. Ceci est fait
en résolvant le problème d'optimisation suivant :

$$
\arg \min_{\mathbf{A} \in \Delta_P, \mathbf{S}, \mathbf{\Psi}} \frac{1}{2} \sum_{n=1}^{N} \left( \|\mathbf{x}_n - \mathbf{S}_n \mathbf{a}_n\|_2^2 + \lambda_{\mathbf{S}} \|\mathbf{S}_n - \mathbf{S}_0 \mathbf{\psi}_n\|_{\mathbf{F}}^2 \right) \quad (3)
$$

où $\mathbf{S} = \{\mathbf{S}_n\}, n = 1, ..., N$, et $\mathbf{A} \in \Delta_P$ signifie que chaque vecteur d'abondance $\mathbf{a}_n \in \mathbb{R}^P$ dans chaque pixel appartient au simplexe unité à $\mathbf{P}$ sommets, c'est-à-dire satisfait les CNA et CSA. $\psi_n \in \mathbb{R}^{P \times P}$ est une matrice diagonale comportant des facteurs d'échelle dépendant du matériaux sur sa diagonale. $\mathbf{\Psi} \in \mathbb{R}^{P \times N}$ est une matrice qui rassemble tous les facteurs d'échelle pour tous les matériaux et pixels. $S_0$ est une matrice d'endmembers de référence. $||.||_F$ est la norme de Frobenius.

L'interprétation géométrique du MMLE est présentée
dans la Fig. 1. Les données se trouvent dans un cône
convexe, dont les arêtes sont les endmembers, et les
instances locales de ceux-ci (non représentés ici) sont
des points proches d'emplacements spécifiques sur ces
droites (les endmembers de référence mis à l'échelle).
Chaque pixel appartient toujours à un simplexe.

Cette formulation s'est avérée utile, mais repose de fa-
çon critique sur les endmembers S₀. Ces derniers doivent
être choisis avec soin car ils fournissent les arêtes du
cône convexe, qui conditionnent tout le processus de
démélange. L'algorithme VCA (Vertex Component Ana-
lysis) (Nascimento et Bioucas Dias, 2005) est l'un des
algorithmes d'extraction les plus connus, qui repose sur la
géométrie convexe du problème de démélange purement
linéaire. Cependant, dans les situations où des ombres ou d'autres forts effets liés à la variabilité sont présents, ces algorithmes, qui recherchent des points extrêmes dans les données, récupèrent souvent des signatures avec une très faible amplitude (et peu de signal). Ceci peut être expliqué dans le cadre du modèle conique qu'est le MMLE. Les spectres de faible amplitude sont proches de l'origine et sont dus aux ombres ; si un modèle fondé sur une géométrie de simplexe est adopté, alors ces points sont en effet extrêmes dans les données. Ils sont de piètres représentants des spectres purs des matériaux et provoquent des erreurs dans l'estimation des paramètres de démélange.

Dans cet article, nous proposons une nouvelle formulation du MMLE et un nouvel algorithme de démélange qui aborde ce problème. Nous montrons que l'utilisation d'un simple k-means avec la mesure de similarité cosinus est capable de fournir des signatures de référence initiales relativement fiables, que nous pouvons affiner dans l'algorithme de démélange, en modélisant les endmembers de référence comme des données directionnelles, c'est à dire des données qui se trouvent sur l'hypersphère unité (Mardia et Jupp, 2009). Un point sur la sphère définit de manière unique une direction dans l'espace des caractéristiques.

## 2. Modèle et algorithme proposés

### 2.1. Présentation de la méthode

D'après le MMLE, un matériau ne doit plus être résumé par un seul point, mais tous les points d'une droite joignant l'origine et une référence représentent le même matériau, avec des variations de luminosité différentes. Idéalement, un endmember devrait alors être vu comme une direction (ou une droite passant par l'origine) dans l'espace des caractéristiques. Dans cet esprit, il est plus logique pour un algorithme d'extraction d'endmembers de rechercher des directions plutôt que des points extrêmes dans les données. Trouver des directions de référence peut être simplement réalisé en utilisant l'algorithme k-means comme une étape de pré-traitement, avec la mesure de similarité cosinus :

$$d(\mathbf{x}_i, \mathbf{x}_j) = 1 - \frac{\mathbf{x}_i^\top \mathbf{x}_j}{\|\mathbf{x}_i\|_2 \|\mathbf{x}_j\|_2}. \quad (4)$$

En effet, cette mesure est insensible aux variations d'échelle des spectres. Fait intéressant, on peut montrer qu'un k-means avec la similarité de cosinus peut être relié à un mélange de distributions de Von Mises-Fisher (un analogue de la distribution gaussienne pour des données directionnelles), d'une manière similaire au fait que le k-means classique utilisant la distance euclidienne est lié à un mélange de distributions gaussiennes (Banerjee et al., 2005). Cela montre l'intérêt d'utiliser l'algorithme k-means et la similarité cosinus dans le cas d'un modèle conique.

Nous pouvons ensuite simplement utiliser les centroïdes des clusters comme endmembers de référence. Cependant, ceux-ci peuvent ne pas être parfaits et trop proches du centre du véritable cône convexe, puisque le

clustering fait des affectations uniques pour chaque pixel. Cela signifie que la moyenne des pixels appartenant à un cluster aura tendance à être biaisée vers le centre du cône en raison de l'influence des pixels mélangés.

Il est alors nécessaire de pouvoir corriger ces derniers à la volée dans l'algorithme de démélange. Une façon de le faire dans la littérature du MML est d'utiliser le volume du simplexe pour régulariser les endmembers (Li et Bioucas-Dias, 2008). Cependant, ce régularisateur est non convexe et difficile à manipuler, et une relaxation convexe plus simple est considérée dans (Berman et al., 2004) en utilisant la somme des distances euclidiennes entre paires d'endmembers possibles. Il n'est pas simple d'adapter cette dernière au modèle conique du MMLE, du fait que chaque membre peut être représenté par n'importe quel point de la droite correspondante. Comparer la distance euclidienne entre deux points quelconques sur deux droites n'est pas pertinent pour définir une distance entre ces deux droites si on leur permet d'avoir des amplitudes différentes (comme nous le confirmerons dans les expériences). L'utilisation de l'angle spectral serait une possibilité, mais cette quantité est assez difficile à gérer dans les problèmes d'optimisation.

Pour résoudre ce problème, nous proposons de modéliser les endmembers de référence comme des données directionnelles, c'est-à-dire en tant que directions dans l'espace de caractéristiques. Il y a plusieurs façons de le faire (Absil et al., 2009), la plus simple étant de contraindre les endmembers de référence à être normalisés : $||s_p||_2 = 1 \forall p$. En effet, chaque droite passant par l'origine est représentée de façon unique par un point sur l'hypersphère unité. De manière équivalente, la matrice des endmembers de référence doit avoir des colonnes de norme 1, c'est-à-dire $S_0 \in OB(L, P)$, la variété dite oblique (Absil et al., 2009). De cette façon, la distance entre les références sur la sphère influe directement sur la position des droites : plus le paramètre de régularisation associé est grand (resp. petit), plus les droites se rapprochent de la droite joignant l'origine et la moyenne des données, avec une moins bonne (resp. meilleure) fidélité aux données. Nous verrons que la position des droites de référence, ainsi que l'importance de l'a priori gaussien sur les endmembers locaux, auront également un impact sur la parcimonie des abondances.

Ainsi, nous raffinons le MMLE pour pouvoir ajuster itérativement la position des droites de référence dans l'espace des caractéristiques. Pour ce faire, nous proposons de minimiser la fonction de coût suivante :

$$ \begin{array}{ll} \underset{A,S,\Psi,S_0}{\operatorname{argmin}} & \frac{1}{2} \sum_{n=1}^{N} (\|\mathbf{x}_n - \mathbf{S}_n \mathbf{a}_n\|_2^2 + \lambda_S \| (\mathbf{S}_n - \mathbf{S}_0 \psi_n) \|_F^2) \\ \text{s.t.} & \mathbf{a}_n \in \Delta_P \quad \forall n \\ & \mathbf{S}_0 \in OB(L, P) \end{array} \quad (5) $$

où l'on note tr la trace d'une matrice, et $V = P I_P - I_P 1_P^\top$ ($I_P$, désigne un vecteur colonne composé de $P$ 1), de sorte que $\operatorname{tr}(\mathbf{S}_0 \mathbf{V} \mathbf{S}_0^\top) = \sum_{i=1}^{P-1} \sum_{j=i+1}^{P} \|\mathbf{s}_{0i} - \mathbf{s}_{0j}\|_2^2$, c'est-à-dire la somme des distances euclidiennes entre les
paires possibles d'endmembers de référence (Berman et al., 2004). $\lambda_S$ et $\lambda_{S_0}$ sont des paramètres de régularisation. Le terme $\lambda_S ||(S_n - S_0\psi_n)||_F^2$ force chaque endmember local à être proche de versions mises à l'échelle de représentants (de norme unité) des directions de références. Autrement dit, il autorise les endmembers locaux à s'éloigner du modèle défini par l'équation (2). Les facteurs d'échelle capturent la variabilité induite par les conditions d'illumination, alors que $S_n$ peut également prendre en compte les effets de variabilité intrinsèque. $\lambda_S$ est directement lié à la variance de l'a priori gaussien (prise ici égale pour tous les endmembers par simplificité). Le fait que les endmembers de référence soient normalisés présente également l'avantage de permettre de comparer facilement l'amplitude des facteurs d'échelle (et donc l'impact de la variabilité induite par l'illumination) pour différents matériaux et images. Des régularisations spatiales peuvent également être ajoutées si nécessaire, comme cela est fait dans (Drumetz et al., 2016).

## 2.2. Optimisation

Ici, nous proposons un algorithme pour obtenir un point stationnaire de la fonction de coût (5). Cette fonction objectif est difficile à minimiser pour plusieurs raisons : elle est non-convexe sur l'ensemble des variables, ce qui nécessite habituellement des méthodes de descente de coordonnées par blocs pour obtenir un minimum local. Dans ce cas, cette approche est rendue encore plus complexe car le problème n'est pas convexe par rapport à $S_0$ non plus, à cause des contraintes de norme unité, qui sont non convexes. Cependant, nous verrons que nous pouvons toujours obtenir un minimum local pour cette variable en profitant de la structure de variété riemannienne de l'ensemble de contraintes. Avant de détailler les différentes étapes de l'optimisation, nous décrivons brièvement comment nous initialisons l'algorithme. Nous utilisons d'abord l'algorithme k-means (avec la similarité cosinus) pour obtenir des centroïdes, que nous normalisons pour initialiser $S_0$. Nous initialisons $S_n$ en affectant la colonne appropriée de cette matrice au pixel courant $x_n$, en fonction de son label de clustering. Les autres colonnes sont initialisées en utilisant les centroïdes résultants. Les matrices d'abondances et de facteurs d'échelle sont initialisées en utilisant l'algorithme SCLSU avec les centroïdes comme références, ce qui est très rapide. De cette façon, nous espérons obtenir un bon minimum local malgré la complexité du problème.

L'optimisation par rapport à A est relativement aisée, la fonction objectif étant différentiable, convexe et il est facile de projeter sur l'ensemble de contraintes (simplexe unité) (Condat, 2016). Le minimum global de ce sous-problème peut alors être obtenu pixel par pixel en utilisant (par exemple) une descente de gradient projeté. Les optimisations par rapport à $S_n$ et $\psi_n$ ne posent aucun problème et possèdent des solutions analytiques (voir (Drumetz et al., 2015) pour plus de détails). L'optimisation de $S_0$ est plus difficile à cause des contraintes d'appartenance aux sphères unité, malgré le caractère différentiable de l'objectif. En utilisant le fait que l'ensemble de

contraintes a une structure de variété riemannienne pour laquelle une rétraction peut être facilement trouvée, nous effectuons une descente de gradient conjugué sur la variété oblique (Absil et al., 2009) (nous utilisons la toolbox Manopt pour MATLAB (Boumal et al., 2014)). La convergence de chaque sous-problème est garantie, mais nous ne pouvons pas prouver la convergence de l'algorithme global vers un point stationnaire de l'objectif, bien qu'en pratique la convergence soit toujours observée.

Nous arrêtons l'algorithme dès que les variations relatives entre deux itérations consécutives de tous les blocs de variables deviennent inférieures à $\epsilon = 10^{-3}$ (en norme). Nous notons que la convergence sera plus lente que pour le MMLE original avec des endmembers de référence fixes, car ces derniers sont maintenant mis à jour de manière itérative et ont un impact sur la géométrie entière du problème.

# 3. Résultats expérimentaux

Dans cette section, nous présentons les résultats obtenus sur un jeu de données synthétiques dont les matériaux présentent une variabilité spectrale réaliste, ainsi qu'un jeu de données réel complexe, avec des endmembers très corrélés et la présence d'une proportion significative de zones d'ombre.

## 3.1. Données synthétiques

### 3.1.1. Génération des données

Nous générons un jeu de données synthétique de la manière suivante. Tout d'abord, nous utilisons la vérité terrain du jeu de données bien connu de l'Université de Pavia<fnref n="1" /> pour nous fournir des spectres labellisés (203 bandes dans les domaines visible et proche infrarouge) appartenant à plusieurs classes d'intérêt, et intégrant leur variabilité spectrale. Nous considérons trois classes présentes dans cette image : la végétation, le béton et les toits métalliques. Ces classes intègrent à la fois la variabilité induite par les variations d'éclairement (les toits et les arbres ont des orientations différentes par rapport au soleil) et les sources de variabilité intrinsèque (notamment dans le béton et la végétation). Dans chaque pixel, nous choisissons les endmembers locaux comme échantillon aléatoire dans chacune de ces classes (après une normalisation pour que chaque représentant se trouve sur la sphère unité).

Des facteurs d'échelle ont été simulés en utilisant un mélange de 4 distributions gaussiennes (ajustées à partir des résultats de SCLSU sur une sous-image de des données de Pavia), ce qui reflète le fait que les facteurs d'échelle proviennent souvent de distributions multimodales (par exemple toits avec différentes orientations, ou zones ombragées).

Les abondances ont été simulées pour être relativement parcimonieuses, en utilisant une distribution de

<footnote n="1">1. http://www.ehu.eus/ccwintco/index.php?title=Hyperspectral_Remote_Sensing_Scenes\sharpPavia_University Newman</footnote>
Dirichlet telle que la densité de probabilité soit concentrée autour des arêtes et des sommets du simplexe unité (tout en comportant toujours une certaine proportion de pixels fortement mélangés).

Les données ont ensuite été générées en utilisant l'équation (2), en ajoutant un bruit blanc gaussien tel que le rapport signal sur bruit soit de 30 dB. L'image générée bénéficie alors de propriétés statistiques réalistes.

### 3.1.2. Résultats

Nous exécutons d'abord l'algorithme VCA+SCLSU pour obtenir rapidement des résultats de démélange prenant en compte la variabilité spectrale, en utilisant des endmembers obtenus via VCA. Nous montrons ci-dessous que cette approche donne de très mauvais résultats. Ensuite, nous nous concentrons sur le test de deux algorithmes avec des références obtenues via k-means : SCLSU et le MMLE de (Drumetz et al., 2015). De plus, nous désignons par MMLE-SDC (Somme des Distances au Carré) le MMLE augmenté de la régularisation convexe sur le volume, mais sans la contrainte de variété oblique. Enfin, nous comparons toutes ces méthodes à la méthode proposée, appelée RMMLE (pour une version Robuste du Modèle de Mélange Linéaire Etendu). Notez que nous ne comparons pas les résultats au démélange utilisant l'algorithme classique FCLSU (Fully Constrained Least Squares Unmixing) (Heinz et Chang, 2001), parce que cet algorithme repose sur une géométrie des données selon un simplexe, et a été mis en défaut dans de multiples scénarios où la variabilité est significative. Pour chaque algorithme, nous ajustons empiriquement les paramètres de régularisation pour obtenir les meilleures performances possibles (les valeurs choisies sont reportées dans le tableau 1). Les résultats quantitatifs sont présentés en utilisant deux mesures : l'erreur quadratique moyenne entre les abondances (aEQM) réelles et celles estimées :

$$ \frac{1}{N\sqrt{P}} \sum_{n=1}^{N} \|\hat{a}_n - a_n\|_2, \qquad (6) $$

et l'angle spectral moyen (sur l'ensemble des pixels et des matériaux) (SAM pour Spectral Angle Mapper)

$$ SAM = \frac{1}{NP} \sum_{n=1}^{N} \sum_{p=1}^{P} \text{acos} \left( \frac{\hat{s}_{pn}^{\top} s_{pn}}{\|\hat{s}_{pn}\|_{2} \|s_{pn}\|_{2}} \right) \quad (7) $$

entre les vrais endmembers de chaque pixel et ceux qui sont récupérés. Ces quantités sont rassemblées dans le tableau 1, ainsi que les temps de calcul sur un ordinateur équipé d'un processeur Intel Core i7-7820HQ CPU @ 2.90GHz, certains calculs ayant été parallélisés sur les 4 cœurs physiques de la machine.

L'approche VCA + SCLSU obtient de très mauvais résultats tant en estimation d'abondance qu'en estimation de variabilité, car en effet deux des signatures extraites sont associées à des pixels avec de petits facteurs d'échelle, et ont une très faible amplitude. Utiliser k-means à la place, avec SCLSU, conduit à de meilleurs résultats, mais loin d'être optimaux car la variabilité n'est expliquée que par des facteurs d'échelle. Le MMLE fait

<table><caption>**TABLE 1:** Résultats quantitatifs sur les données synthétiques. À l'exception de VCA + SCLSU, tous les algorithmes utilisent k-means pour obtenir la matrice de référence initiale. Les valeurs des paramètres de régularisation sont signalées le cas échéant.</caption><thead><tr><th></th><th>λ<sub>S</sub></th><th>λ<sub>S0</sub></th><th>aEQM</th><th>SAM (degrés)</th><th>Temps (s)</th></tr></thead><tbody><tr><td>VCA+SCLSU</td><td>×</td><td>×</td><td>0,2075</td><td>54,4</td><td>3</td></tr><tr><td>SCLSU</td><td>×</td><td>×</td><td>0,0654</td><td>6,32</td><td>2</td></tr><tr><td>MMLE</td><td>0,01</td><td>×</td><td>0,0642</td><td>5,62</td><td>18</td></tr><tr><td>MMLE+SDC</td><td>0,1</td><td>0,25</td><td>0,1718</td><td>10,41</td><td>88</td></tr><tr><td>RMMLE</td><td>0,1</td><td>0,5</td><td><b>0,0560</b></td><td><b>3,48</b></td><td>428</td></tr></tbody></table>

encore mieux à cause de l'a priori gaussien supplémentaire. On voit toute l'importance de la contrainte de variété oblique sur les deux dernières lignes du tableau : MMLE + SDC échoue car les termes de régularisation impliquent la comparaison de références avec des amplitudes potentiellement différentes, alors que l'introduction de la contrainte conduit aux meilleurs résultats. Nous montrons dans la Figure 2 des résultats qualitatifs en utilisant des diagrammes de dispersion des données (en utilisant les trois premières composantes principales) avec les endmembers obtenus pour les trois meilleurs algorithmes (les autres diagrammes montrent que les droites sont beaucoup trop loin du vrai cône pour être pertinentes). Des conclusions similaires peuvent être tirées de cette figure, montrant que la méthode RMMLE est capable de trouver les meilleurs endmembers dans chaque pixel.

## 3.2. Données réelles

Le jeu de données réelles que nous utilisons a été acquis en 2009 par Japan Space Systems au dessus du Tama Forest Science Garden dans la région ouest de Tokyo, avec le capteur CASI-3 (72 bandes spectrales dans les domaines visible et proche infrarouge) (Matsuki et al., 2015). La résolution spatiale est de 1m. L'image que nous utilisons est un sous-ensemble de taille 207 × 268 × 72 de toute la scène. Une représentation RVB est montrée en Figure 3 (a). Ce jeu de données a été utilisé pour la classification supervisée d'espèces d'arbres, en utilisant une vérité terrain et des données LiDAR comme une caractéristique de classification supplémentaire, puisque les différentes espèces d'arbres sont spectralement très proches les unes des autres. L'image comprend également de nombreuses zones d'ombre en raison des coufonnes des arbres, qui constituaient un obstacle important dans les études précédentes (Matsuki et al., 2015). En outre, d'autres matériaux non végétaux sont présents, tels que des toits, des routes et de la terre. Nous montrons ici que l'utilisation de k-means à la place de VCA permet de faire la distinction entre les conifères et les feuillus de manière totalement non supervisée. Quelques conifères et feuillus labellisés sont montrés dans la Figure 3 (b). Nous montrons le diagramme de dispersion des données et des pixels purs étiquetés dans la Figure 3(c).

Nous démélangeons les données en utilisant $P = 4$ matériaux avec les algorithmes mentionnés ci-dessus. Nous montrons sur la Figure 4 les diagrammes de dispersion des données et des endmembers récupérés pour
<figure><img src="image_2.png" /><figcaption>**FIGURE 2:** Dispersion des données (bleu), des vrais endmembers (noirs) et de ceux extraits (rouges) pour (a) SCLSU (b) MMLE (c) RMMLE.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>**FIGURE 3:** (a) Réprésentation RVB des données. (b) Vérité terrain pour les conifères (vert) and les feuillus (rouge). (c) Dispersion des données, et des pixels de vérité terrain (même code couleur).</figcaption></figure>

VCA + SCLSU, SCLSU, MMLE et RMMLE. Les cartes d'abondance sont représentées sur la Figure 5. Comme dans le cas des données synthétiques, les endmembers récupérés par VCA sont mauvais à cause des zones d'ombre de l'image, et les abondances correspondantes n'ont aucun sens. La plupart des données sont projetées sur la droite la plus proche du cône identifié, qui représente la végétation. L'utilisation de k-means permet en revanche de distinguer les conifères des feuillus. La végétation rase (herbe, gazon) et les ombres sont également caractérisés par des valeurs faibles et grandes des facteurs d'échelle, respectivement. Les facteurs d'échelles obtenus par l'algorithme pour chaque matériau sont représentés en Figure 6. En effet la végétation rase est spectralement très proche des deux classes d'arbres identifiées, au point que les directions correspondantes se superposent (les pixels d'herbe correspondent au prolongement des droites rouge et vertes sur la Figure 3(c). En revanche, l'herbe a tendance a mieux réfléchir la lumière, et si la forme spectrale est proche d'autres types de végétation, la réflectance est bien plus élevée, ce que l'on peut retrouver dans les facteurs d'échelle des deux classes de végétation. La faible luminosité des zones d'ombres est absorbée par les facteurs d'échelle, si bien que l'algorithme peut traiter ces pixels comme des pixels classiques. Les abondances de SCLSU et du MMLE sont

plutôt similaires, légèrement plus parcimonieuses pour le MMLE, car il est capable de mieux capturer les effets de variabilité que SCLSU (comme le montre la Figure 4 (c). L'approche RMMLE étant capable d'ajuster les références, il est possible d'obtenir des cartes d'abondance plus parcimonieuses qui correspondent de près à la vérité terrain de la Figure 3(b). Nous voyons que les droites identifiées embrassent très bien le jeu de données et sont les plus proches des pixels de vérité terrain de la Fig. 3 (c). Les valeurs des paramètres de régularisation et les temps de calcul sont donnés à titre indicatif dans le tableau 2.

<table><caption>**TABLE 2:** Résultats quantitatifs sur les données réelles. À l'exception de VCA + SCLSU, tous les algorithmes utilisent k-means pour obtenir la matrice de référence initiale. Les valeurs des paramètres de régularisation sont signalées le cas échéant.</caption><thead><tr><th></th><th>λ<sub>S</sub></th><th>λ<sub>S0</sub></th><th>Temps (s)</th></tr></thead><tbody><tr><td>VCA+SCLSU</td><td>×</td><td>×</td><td>9</td></tr><tr><td>SCLSU</td><td>×</td><td>×</td><td>11</td></tr><tr><td>MMLE</td><td>0,01</td><td>×</td><td>243</td></tr><tr><td>RMMLE</td><td>0,5</td><td>100</td><td>757</td></tr></tbody></table>
<figure><img src="image_4.png" /><figcaption>**FIGURE 4:** Dispersion des données (bleu) et des endmembers extraits (rouge) pour (a) VCA+SCLSU (b) SCLSU (c) MMLE (d) RMMLE.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**FIGURE 5:** Cartes d'abondances obtenues par les algorithmes testés.</figcaption></figure>

## 4. Conclusion

Dans cet article, nous avons proposé un nouvel algorithme pour démélanger des données hyperspectrales en tenant compte à la fois de la variabilité induite par l'illumination et de la variabilité intrinsèque des endmembers. Cet algorithme est capable d'obtenir de meilleures estimations des endmembers de référence que VCA en profitant du fait que les endmembers sont bien modélisés par des données directionnelles. Les signatures de référence sont contraintes à être sur la sphère unité puisque chaque point sur cette dernière représente complètement une droite dans l'espace de caractéristiques. Nous obtenons des estimations robustes des arêtes du cône convexe dans le cadre du modèle de mélange linéaire étendu. Les résultats sur les données synthétiques montrent la pertinence de l'approche et de très bonnes performances sont également observées sur un jeu de données réelles com-

plexe. Les travaux futurs comprendront un moyen d'estimer automatiquement les paramètres dépendant des matériaux (tels que les variances des a priori gaussiens), et des tests sur les différentes façons de représenter les endmembers en utilisant des données directionnelles.

## 5. Remerciements

Ce travail a été financé par l'Agence Nationale de la Recherche et la Direction Générale de l'Armement, par le projet ANR-DGA APHYPIS, ANR-16 ASTR-0027-01. L. Drumetz a également bénéficié d'une subvention du programme d'été de la Japanese Society For the Promotion of Science, JSPS-SP17206 et d'une bourse de mobilité post-doctorale Campus France, PRESTIGE-2016-4 0006.
<figure><img src="image_6.png" /><figcaption>**FIGURE 6:** Facteurs d'échelle obtenus par l'aglorithme RMMLE sur les données réelles. Un facteur d'échelle proche de 1 signifie que l'endmember local est proche de la référence. Dans les zones d'abondance faible pour un matériau, les facteurs d'échelle on tendance à rester proches de 1, leur valeur initiale.</figcaption></figure>

Références

Absil, P.-A., Mahony, R., Sepulchre, R., 2009. Optimization algo-
rithms on matrix manifolds. Princeton University Press, Prin-
ceton, NJ, USA.

Banerjee, A., Dhillon, I. S., Ghosh, J., Sra, S., 2005. Clustering on the unit hypersphere using von Mises-Fisher distributions. Journal of Machine Learning Research 6, 1345–1382.

Berman, M., Kiiveri, H., Lagerstrom, R., Ernst, A., Dunne, R., Huntington, J. F., 2004. ICE : A statistical approach to identifying endmembers in hyperspectral images. IEEE Transactions on Geoscience and Remote Sensing 42 (10), 2085–2095.

Bioucas-Dias, J., Plaza, A., Camps-Valls, G., Scheunders, P., Nasrabadi, N., Chanussot, J., 2013. Hyperspectral remote sensing data analysis and future challenges. IEEE Geoscience and Remote Sensing Magazine 1 (2), 6–36.

Bioucas-Dias, J., Plaza, A., Dobigeon, N., Parente, M., Du, Q., Gader, P., Chanussot, J., 2012. Hyperspectral unmixing overview: Geometrical, statistical, and sparse regression-based approaches. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing 5 (2), 354–379.

Boumal, N., Mishra, B., Absil, P.-A., Sepulchre, R., 2014. Manopt, a Matlab toolbox for optimization on manifolds. Journal of Machine Learning Research 15, 1455–1459.
URL http://www.manopt.org

Condat, L., 2016. Fast projection onto the simplex and the $L_1$
ball. Mathematical Programming 158 (1-2), 575–585.

Dobigeon, N., Tourneret, J. Y., Richard, C., Bermudez, J. C. M., McLaughlin, S., Hero, A. O., 2014. Nonlinear unmixing of hyperspectral images : Models and algorithms. IEEE Signal Processing Magazine 31 (1), 82–94.

Drumetz, L., 2016. Endmember variability in hyperspectral image unmixing. Thèse de doctorat, Université de Grenoble Alpes, France.

Drumetz, L., Henrot, S., Veganzones, M. A., Chanussot, J., Jutten, C., 2015. Blind hyperspectral unmixing using an extended linear mixing model to address spectral variability. Dans : IEEE Workshop on Hyperspectral Image and Signal Processing : Evolution in Remote Sensing (WHISPERS 2015). Tokyo, Japan, pp. 1–4.

Drumetz, L., Tochon, G., Veganzones, M., Chanussot, J., Jutten, C., 2017. Improved local spectral unmixing of hyperspectral data using an algorithmic regularization path for collaborative sparse regression. Dans : Acoustics, Speech and Signal Processing (ICASSP), 2017 IEEE International Conference on. La Nouvelle Orléans, LA, USA, pp. 6190–6194.

Drumetz, L., Veganzones, M. A., Henrot, S., Phlypo, R., Chanussot, J., Jutten, C., 2016. Blind hyperspectral unmixing using an extended linear mixing model to address spectral variability. IEEE Transactions on Image Processing 25 (8), 3890–3905.

Halimi, A., Dobigeon, N., Tourneret, J. Y., 2015. Unsupervised unmixing of hyperspectral images accounting for endmember variability. IEEE Transactions on Image Processing 24 (12), 4904–4917.

Hapke, B., 2012. Theory of reflectance and emittance spectroscopy. Cambridge University Press, Cambridge, UK.

Heinz, D., Chang, C.-I., 2001. Fully constrained least squares linear spectral mixture analysis method for material quantification in hyperspectral imagery. IEEE Transactions on Geoscience and Remote Sensing 39 (3), 529–545.

Henrot, S., Chanussot, J., Jutten, C., 2016. Dynamical spectral unmixing of multitemporal hyperspectral images. IEEE Transactions on Image Processing 25 (7), 3219–3232.

Heylen, R., Parente, M., Gader, P., 2014. A review of nonlinear hyperspectral unmixing methods. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing 7 (6), 1844–1868.

Keshava, N., Mustard, J. F., 2002. Spectral unmixing. IEEE Signal Processing Magazine 19 (1), 44–57.

Li, J., Bioucas-Dias, J., 2008. Minimum volume simplex analysis : A fast algorithm to unmix hyperspectral data. Dans : Geoscience and Remote Sensing Symposium, 2008. IGARSS 2008. IEEE International. Boston, MA, USA.

Mardia, K. V., Jupp, P. E., 2009. Directional statistics. Vol. 494. John Wiley & Sons, Hoboken, NJ, USA.

Matsuki, T., Yokoya, N., Iwasaki, A., 2015. Hyperspectral tree species classification of japanese complex mixed forest with the aid of LiDAR data. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing 8 (5), 2177–2187.

Nascimento, J., Bioucas Dias, J., 2005. Vertex component analysis : a fast algorithm to unmix hyperspectral data. IEEE Transactions on Geoscience and Remote Sensing 43 (4), 898–910.

Somers, B., Asner, G., Tits, L., Coppin, P., 2011. Endmember variability in spectral mixture analysis : A review. Remote Sensing of Environment 115 (7), 1603 – 1616.

Somers, B., Zortea, M., Plaza, A., Asner, G. P., 2012. Automated extraction of image-based endmember bundles for improved spectral unmixing. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing 5 (2), 396–408.

Thouvenin, P.-A., Dobigeon, N., Tourneret, J.-Y., 2016. Hyper-spectral unmixing with spectral variability using a perturbed linear mixing model. IEEE Transactions on Signal Processing 64 (2), 525–538.

Veganzones, M., Drumetz, L., Marrero, R., Tochon, G., Dalla Mura, M., Plaza, A., Bioucas-Dias, J., Chanussot, J., 2014. A new extended linear mixing model to address spectral variability. Dans : Proc. IEEE Workshop on Hyperspectral Image and Signal Processing : Evolution in Remote Sensing (WHISPERS). Lausanne, Suisse.

Zare, A., Ho, K., 2014. Endmember variability in hyperspectral analysis : Addressing spectral variability during spectral un-mixing. IEEE Signal Processing Magazine 31 (1), 95–104.