# VISUALISATION DES CHANGEMENTS SUR SÉRIES TEMPORELLES
RADAR : MÉTHODE REACTIV ÉVALUÉE À L'ÉCHELLE MONDIALE SOUS
GOOGLE EARTH ENGINE

Elise Colin-Koeniguer¹, Jean-Marie Nicolas², Béatrice Pinel-Puysségur³, Jean-Michel Lagrange³, Fabrice
Janez¹

1. Onera, DTIS, Chemin de la Hunière, 91123 Palaiseau

2. LTCI, Télécom ParisTech, Université Paris-Saclay, 75013, Paris

3. CEA, DAM, DIF, 91297 Arpajon Cedex, France

## Résumé

Cet article présente une méthode de visualisation d'une pile temporelle d'images SAR, appelée REACTIV, qui permet de faire ressortir en couleur les zones ayant subi des changements sur la période de temps observée. Cette méthode a été testée à grande échelle grâce à la plateforme Google Earth Engine. Elle est fondée sur l'espace de couleur TSV (pour Teinte Saturation Valeur ou, en anglais, HSV pour Hue Saturation Value) et n'exploite que des estimations dans le domaine temporel, sans aucune estimation spatiale. La saturation des couleurs est codée par le coefficient de variation temporel, dont plusieurs propriétés statistiques sont explicitées. Les limites de l'utilisation de la plateforme Google Earth Engine sont évaluées, et plusieurs cas d'applications sont proposés : agriculture, dynamique urbaine, trafic maritime.

**Mots clés :** Piles temporelles, SAR, coefficient de variation temporel, visualisation, détection de changements.

## Abstract

This article presents a method for visualizing a temporal stack of SAR images, that highlights in color the areas where changes occurred over the observed period. The method has been tested worldwide using the Google Earth Engine platform. It is based on the color model HSV (for Hue-Saturation-Value) and exploits the temporal coefficient of variation for color saturation, without any spatial average. The statistical properties of this parameter are explained. The Google Earth Engine platform is evaluated, and several applications of the visualization are proposed : agriculture, urban expansion, maritime traffic.

**Keywords :** Time-series, SAR, temporal coefficient of variation, visualization, change detection

## 1. Introduction

Depuis 2014, suite au lancement du programme d'ac- quisition Copernicus de l'ESA avec la mise en accès libre des produits, les traitements de télédétection et la re- cherche de nouvelles applications potentielles connaissent un succès sans précédent. Parmi les potentialités offertes par ce changement de contexte, on trouve notamment l'exploitation de l'axe temporel : là où autrefois l'obten- tion d'une série temporelle d'images sur un même site était rare et difficile, aujourd'hui, de nombreuses plate- formes de téléchargement permettent de constituer de telles piles très rapidement et partout dans le monde. Parallèlement à l'obtention de ces données, les solutions permettant de les traiter avec de grosses puissances de calcul se développent.

Les images SAR ont des propriétés temporelles re- marquables : comme les images ne souffrent pas de la couverture nuageuse, et que le radar est un système ac- tif, le signal est stable entre deux acquisitions en condi- tions interférométriques. Aussi, parmi les intérêts de dis- poser de longues séries temporelles, on trouve notam- ment le fait d'en améliorer la qualité, par simple filtrage

temporel, ou par des méthodes plus évoluées, mais aussi la détection de changements ou plus largement, le suivi d'activité. En contrepartie, les images SAR, de par leur nature cohérente, sont sujettes au phénomène de spe- ckle.

Lorsque l'on parle de détection de changements, on sous-entend un principe de détection binaire des zones modifiées entre deux dates. Avant même la décision de détection de changements, disposer d'une méthode de visualisation des images mettant en relief ces informa- tions est utile pour un opérateur. Classiquement, la vi- sualisation des changements passe par l'utilisation de canaux colorés complémentaires (vert et magenta, rouge et cyan ou jaune et bleu). Les deux canaux codent les in- tensités radar respectives des images du couple. De ce fait, la présence de comportements différents entre deux images se traduit par une couleur, tandis que les pixels d'amplitude stable apparaissent en niveaux de gris. Cette méthode peut être extrapolée au cas de trois images seulement, puisque l'espace RGB est à trois dimensions. En revanche, pour N>3 images, la visualisation des chan- gements n'est plus aussi immédiate. Il est possible de passer par des animations, comme dans (Koeniguer, 2017)
dont l'approche permet de différencier les changements de type permanents de ceux qui sont des changements ponctuels. Cependant, ces produits restent des objets visdeos, qui n'adressent pas la problématique de réduction temporelle permettant d'aboutir à un seul produit résuivant l'ensemble des changements survenus. Dans ce papier, nous cherchons à trouver une telle représentation basée sur une composition colorée et qui exploite au mieux l'axe temporel pour le suivi d'activité.

Plusieurs travaux de la littérature proposent des compositions colorées qui utilisent l'information de cohérence interférométrique d'une pile temporelle. Parmi eux, Boccardo et al. (2015) proposent une telle approche appelée MTC (Multi Temporal Change) pour la visualisation des zones d'extension de lave. Les canaux rouge et vert codent les amplitudes des images du couple interférométrique et le canal bleu le niveau de cohérence interférométrique. Cette composition colorée n'adresse pas réellement la réduction de la dimension temporelle, les auteurs préconisant simplement de choisir le couple minimisant à la fois la ligne de base temporelle et géométrique. De la même manière, Wegmuller et al. (2015) montrent plusieurs compositions colorées RGB où les canaux codés sont soit une amplitude, soit une cohérence calculée sur une paire interférométrique Sentinel-1 à 12, 24 et 48 jours de lignes de base temporelles. D'autres compositions colorées utilisent simplement en entrée les images de trois cohérences interférométriques. Schulz et al. (2010) proposent une représentation appelée CovAmCoh, là encore permettant de visualiser des changements sur une paire d'images en conditions interférométriques. Le canal rouge est codé par le coefficient de variation de l'image SAR, s'agissant ici du coefficient de variation spatial et non temporel. Ce paramètre donne alors une information sur l'hétérogénéité locale. Comme deux images sont considérées, il est possible de faire la moyenne des deux cartes de coefficients de variation calculés localement. Le canal vert code l'amplitude, et le bleu, la cohérence interférométrique. Toutes ces compositions ne font qu'illustrer la pertinence de l'information contenue dans le niveau de cohérence interférométrique.

Les premiers travaux qui adressent réellement l'utilisation de la dimension temporelle sont présentés dans (Amitrano et al., 2017). Ils proposent, dans un premier produit appelé de niveau 1α, de séparer les intensités radar liées à la saison sèche codées en bleu, de celles liées à la saison humide, codées en vert. Enfin, le canal rouge code la cohérence interférométrique, ce qui permet de faire ressortir dans cette couleur les cibles anthropiques. Toujours dans (Amitrano et al., 2017), les auteurs proposent un produit RGB de visualisation appelé de niveau 1β, «supérieur» au sens où il fait intervenir dans le codage un paramètre statistique temporel qui est celui de la variance des intensités radar temporelles. Le rouge code cette variance, le canal vert code l'intensité radar moyenne, et enfin le canal bleu est une combinaison de la cohérence interférométrique moyenne et d'un coefficient proportionnel à l'intensité radar. Ce dernier produit 1β utilise donc la notion de variance temporelle.

Plus récemment, l'ONERA a proposé dans (Koeniguer et al., 2018) une représentation colorée de piles temporelles de données SAR. Cette représentation est destinée à faire apparaître en couleur les pixels qui sont le siège de changements potentiels, et à laisser en niveaux de gris les pixels inchangés. Cette approche se distingue des travaux précédents sur au moins trois critères :

- Un paramètre temporel crucial utilisé est le coefficient de variation temporel, qui possède des propriétés statistiques remarquables. Aucune estimation spatiale n'est réalisée dans la composition.

- La composition colorée passe par l'utilisation de l'espace de représentation couleur HSV, avant une conversion dans l'espace RGB classique.

- Enfin, en entrée, elle n'utilise que les amplitudes, et peut donc être testée sur des données géoréférencées de type Sentinel 1 GRD.

Dans cet article, l'ambition est double. Il s'agit tout d'abord d'expliciter les propriétés statistiques temporelles du critère exploité dans notre approche et qui expliquent en partie son efficacité. Il s'agit ensuite d'analyser les résultats à l'échelle de la planète, suite à l'implémentation de notre méthode sur la plateforme Google Earth Engine, dénommée par la suite par son acronyme GEE. Ce travail d'implémentation est l'occasion de faire un retour d'expérience sur cette plateforme.

L'article s'articule de la manière suivante. Dans la section 2, nous présentons le principe de la visualisation proposée. En section 3, nous adressons le cadre théorique statistique permettant de modéliser le coefficient de variation utile pour la visualisation. Ces résultats permettent à la fois de mieux comprendre l'efficacité de la méthode et de proposer son paramétrage. En section 4, les limites et pistes d'améliorations potentielles à apporter à la plateforme GEE sont abordées. En section 5, plusieurs types de résultats sont analysés, dont découlent des applications potentielles différentes. Enfin, en section 6, une conclusion nous permet de synthétiser les principaux résultats obtenus et d'aborder les perspectives d'évolutions.

## 2. La méthode de visualisation REACTIV

### 2.1. Principe général

Le principe de la méthode est de passer dans l'espace colorimétrique HSV, où la teinte du canal H est représentative du temps, le canal de saturation S sert à mettre en relief les couleurs là où le changement est significant, et la valeur V de luminosité correspond à une notion d'intensité radar classique. Pour cette raison, nous appelons notre méthode de visualisation REACTIV, pour *Rapid and EAsy Change detection in radar Time-series by Variation coefficient*.

En pratique, une difficulté concerne le choix des valeurs de teinte (canal H) et de luminosité (canal V) en raison de la diversité des changements qu'il est possible
de rencontrer dans une série temporelle. En effet, pour N dates, différents évènements peuvent expliquer un changement, à savoir : une apparition ponctuelle d'un objet comme un véhicule, une apparition ou disparition permanente comme la construction ou la destruction d'un bâtiment, ou des événements plus complexes comme des variations saisonnières telles que rencontrées sur des parcelles végétales. Etant donnée la diversité de ces situations, il n'existe pas une solution de codage unique mais plusieurs solutions devant s'adapter aux types d'événements que l'on souhaite mettre en avant.

Actuellement, le choix retenu pour la luminosité (canal V) est l'amplitude radar maximale atteinte dans le profil temporel. Son avantage est de maximiser les chances de voir l'événement quel qu'il soit. En revanche, il ne permet pas de s'affranchir du phénomène de specke. Une alternative consisterait à utiliser la moyenne des amplitudes radar des images de la pile, en vue de bénéficier en même temps d'un filtrage de speckle et donc d'aboutir à une visualisation plus "cosmétique". Cependant, un tel choix a tendance à faire disparaître rapidement les évènements ponctuels, et de fait, ne doit être envisagé que pour le cas d'évènements de type saisonnier ou rémenent.

Pour la teinte (canal H), le choix a été fait d'associer une couleur à la date pour laquelle l'intensité radar maximale est atteinte. Ce choix conviendra particulièrement pour un événement ponctuel, par exemple en cas de présence d'un bateau. Dans le cas de l'apparition d'un objet, la couleur sera représentative de la date à laquelle le signal sera maximal ; ce ne sera pas nécessairement la date d'apparition de l'objet mais une des dates à laquelle il est particulièrement fort. De même, dans le cas de la destruction d'un objet, la date appartiendra à l'une des dates auxquelles l'objet était visible. Une alternative serait d'utiliser un algorithme de détection de ruptures dans une série temporelle afin de localiser plus précisément la date de changement. Cet algorithme est aujourd'hui à l'étude, mais sera nécessairement plus coûteux en temps de calcul que la simple recherche du maximum dans le profil.

Enfin, la saturation (canal S) est codée par le coefficient de variation temporel dont l'intérêt est de réagir à la présence de changements. Ainsi, la saturation de la couleur sera d'autant plus forte que le signal temporel attaché à un pixel aura une évolution temporelle spécifique.

## 2.2. Réglages des dynamiques

Afin de mettre au point le codage couleur, les paramètres réels utilisés pour la visualisation doivent être définis entre 0 et 1. Ici, nous explicitons les choix réalisés sur les paramètres de teinte et de luminosité. Le choix relatif au coefficient de variation s'appuie sur l'analyse des propriétés statistiques de ce paramètre ; il sera donc explicité dans la section suivante.

### 2.2.1. La Teinte H

La composition colorée affecte une couleur à une date précise. Les piles constituées ne sont pas forcément acquises avec un échantillonnage temporel régulier : dès lors qu'une acquisition manque dans la pile, cela provoque un saut dans les dates d'acquisitions. Aussi, notre choix consiste à paramétrer la date linéairement par rapport à la plage d'acquisition choisie. Si $t_1$ et $t_2$ sont les dates de début et de fin de la pile à visualiser, alors nous calculons, pour chaque image acquise à la date $t$, la fraction temporelle $f_t$ :

$$f_t = \frac{t - t_1}{t_2 - t_1} \qquad (1)$$

Ce paramètre $f_t$ compris entre 0 et 1 sert directement à coder tel quel le canal H de l'image. Il faut simplement noter, dans cette composition, que la première et la dernière date ont des couleurs très proches car la palette de couleur HSV est cyclique. D'autres alternatives seraient possibles : on peut restreindre le choix des couleurs à l'intérieur d'un secteur angulaire de la palette HSV, ou encore affecter une couleur fixe à chaque saison, dans le cas où la série couvre plusieurs années, afin de visualiser les effets saisonniers.

### 2.2.2. La luminosité V

La luminosité est codée par le canal d'amplitude radar maximale. S'il n'y a aucun changement, et que l'ergodicité temporelle est respectée, alors cette amplitude suit une loi classique de Rayleigh-Nakagami. Aussi, les choix de dynamique effectués peuvent se faire de manière classique en imagerie radar, en seuillant l'amplitude radar à une valeur maximale de dynamique. Pour les besoins de visualisation, nous avons converti les valeurs stockées en dB par la plateforme GEE en niveaux linéaires d'amplitudes. Ensuite, nous avons seuillé la dynamique d'amplitude à 1 (soit 0 dB, seuil imposé par GEE), et enfin, appliqué un facteur de puissance choisi manuellement à 0,8 pour la polarisation VH, et 1,4 pour la polarisation VV, dans le but de rendre plus similaires les contrastes présents dans ces deux polarisations, qui sont naturellement différents.

## 3. Dynamique du coefficient de variation : aspects théoriques

L'originalité de la composition REACTIV repose principalement sur l'emploi du coefficient de variation temporel pour calculer le paramètre de saturation. C'est ce paramètre qui fera ou non que notre œil interprète la couleur comme liée à un changement. Cette section s'intéresse aux propriétés statistiques de ce coefficient.

### 3.1. Coefficient de variation d'une zone de speckle

Dans les images radar, il est communément admis que l'amplitude d'une loi de speckle sans texture suit une loi de Rayleigh-Nakagami (Goodman, 1976) :

$$RN[\mu, L](u) = \frac{2}{\mu} \frac{\sqrt{L}}{\Gamma(L)} \left( \frac{\sqrt{L}u}{\mu} \right)^{2L-1} e^{-(\frac{\sqrt{L}}{\mu} u)^2}, (2)$$
où $\mu$ est un paramètre d'échelle et $L$ est un paramètre de forme.

Généralement, les paramètres de cette loi sont estimés spatialement sur une zone homogène. Dans notre approche, nous nous intéressons aux statistiques temporelles. Nous faisons alors l'hypothèse que, pour un pixel appartenant à une zone de speckle de loi Rayleigh Nakagami, et n'ayant subi aucun changement, les différentes réalisations d'amplitudes au cours du temps de ce pixel suivent également une loi de Rayleigh Nakagami. Cela revient en quelque sorte à transformer l'hypothèse de stationnarité spatiale en une hypothèse de stationnarité temporelle. Dans ce cas, on parlera dans la suite de speckle stable.

De la loi de Rayleigh Nakagami, il est possible de dériver l'expression des moments empiriques. On a notamment (Nicolas, 2006) :

$$m_1 = \mu \frac{\Gamma(L + \frac{1}{2})}{\sqrt{L\Gamma(L)}} \text{ et } m_2 = \mu^2, \quad (3)$$

ce qui permet de trouver l'expression suivante pour le coefficient de variation, défini comme le rapport de l'écart-type et de la moyenne :

$$\gamma = \frac{\sigma}{\mu} = \frac{\sqrt{m_2 - m_1^2}}{m_1} = \sqrt{\frac{\Gamma(L)\Gamma(L+1)}{\Gamma(L+1/2)^2}} - 1. \quad (4)$$

Cette expression fait apparaître une première propriété intéressante : le coefficient de variation aura la même valeur pour toutes les zones de speckle stables, et ce, quelle que soit l'amplitude moyenne de ce speckle.

Dans le cas particulier $L = 1$, on trouve $\gamma = 0,522723$. L'application d'un multi-look sur une image SAR a pour effet de modifier le paramètre $L$. Les données Sentinel-1 GRD sont délivrées avec un nombre de vues équivalent affiché de $L = 4,9$, calculé pour la valeur théorique de milieu de fauchée et en milieu d'orbite (Sentinel 1 Team, 2013). Cette valeur de $L$ donne pour le coefficient de variation associé au speckle $\gamma = 0,2286$.

Afin de savoir comment l'estimation de ce paramètre se comporte, nous avons simulé une pile d'images de speckle suivant une loi de Rayleigh Nakagami, à l'aide du logiciel SARLAB (Nicolas, 2017), qui passe par l'expression de la fonction de répartition de la loi de Nakagami sous la forme de fonction de Meijer (Nicolas et Tupin, 2016). Une pile de $N$ images de 1024 x 1024 pixels peut être ainsi simulée, et il alors est possible de calculer le coefficient de variation selon l'axe temporel, en fonction du nombre $N$ d'images choisies. On obtient ainsi une image de coefficient de variation que l'on traite comme n'importe quelle autre image : en particulier on peut, en fonction du nombre $N$ d'images générées selon l'axe temporel, calculer la valeur moyenne $m_\gamma$ du coefficient de variation.

La figure 1 montre le résultat de cette recherche pour $N$ variant de 5 à 100, et pour les valeurs $L = 1$ (speckle pleinement développé) et $L = 4,9$. Le point que l'on peut relever est que la valeur moyenne $m_\gamma$ semble converger assez vite vers les valeurs théoriques 0,522723 et

<figure><img src="image_1.png" /><figcaption>FIGURE 1: Statistiques simulées de la valeur moyenne du coefficient de variation sur un speckle "stable", en fonction de $N$, nombre d'images dans la pile.</figcaption></figure>

0,2286. L'estimateur semble donc à première vue sans biais.

Nous nous intéressons maintenant à la variance de cet estimateur. Kendall et Stuart (1977) proposent de calculer la variance d'une fonction $g(m_1, m_2)$ en effectuant un développement limité au premier ordre de la fonction autour des valeurs $m_{0,1}$ et $m_{0,2}$. Cette méthode mène à trouver la formule suivante, proposée dans (Nicolas, 2006) :

$$var(\gamma) = \frac{1}{4N} \frac{4m_2^3 - m_2^2 m_1^2 + m_1^2 m_4 - 4m_1 m_2 m_3}{m_1^4 (m_2 - m_1^2)}, \quad (5)$$

où $N$ est le nombre d'images considérées. Connaissant les moments d'ordre 1 à 4 de la loi de Nakagami, il est possible d'écrire la variance de cet estimateur en fonction de $L$ :

$$var(\gamma) = \frac{1}{4N} \frac{L v^4 (4L^2 v^2 - 4Lu^2 - u^2)}{u^4 (L v^2 - u^2)}, \quad (6)$$

avec $u = \Gamma(L + \frac{1}{2})$ et $v = \Gamma(L)$.

Dans le cas particulier $L = 1$, on trouve pour la variance la formule $var(\gamma) = 0,137881/N$, et donc pour l'écart-type $\sigma_\gamma = 0,3713/\sqrt{N}$. Pour la valeur $L = 4,9$ du nombre de vues équivalent affiché des données Sentinel-1 GRD, on obtient $var(\gamma) = 0,0216/N$ et $\sigma_\gamma = 0,1616/\sqrt{N}$.

Cette formule montre une seconde propriété intéressante sur ce paramètre : l'écart-type du coefficient de variation décroît avec le nombre d'images $N$ en $\frac{1}{\sqrt{N}}$.

A nouveau, ces formules ont été complétées par des simulations suivant le même schéma que précédemment. On peut calculer la variance du coefficient de variation empirique d'autant plus raisonnablement que le nombre de pixels est suffisamment grand pour faire une étude statistique raisonnable (on a un effet à peu près $10^6$ valeurs).

Les résultats de simulation de la figure 2 montrent une bonne adéquation entre les valeurs théoriques de l'écart-type et les valeurs empiriques calculées sur nos images simulées. Enfin, plusieurs simulations réalisées montrent que le coefficient de variation, calculé sur des amplitudes suivant une loi de Rayleigh Nakagami, semble suivre à nouveau une loi de Rayleigh Nakagami.
<figure><img src="image_2.png" /><figcaption>**FIGURE 2:** Ecart-types empiriques (simulés) et théoriques du coefficient de variation sur un speckle "stable", en fonction de N, nombre d'images dans la pile.</figcaption></figure>

## 3.2. Coefficient de variation d'un profil de speckle avec rupture temporelle

Maintenant que nous avons analysé des propriétés du coefficient de variation pour des zones de speckle pur, nous proposons de simuler la manière dont ce coefficient se comporte sur une zone de speckle qui comporte une rupture temporelle.

Cette simulation est réalisée de la manière suivante. On se fixe un nombre *N* d'échantillons temporels, et on génère alors $10^6$ profils de taille *N*, suivant une loi de Rayleigh Nakagami de paramètres $\mu = 0,3$ (-11 dB) et $L = 4,9$. On génère également le même nombre de profils pour lesquels nous remplaçons la dernière valeur du profil par une amplitude fixe, en considérant des valeurs d'amplitudes échelonnées entre $m_1 + 2,5$ dB et $m_1 + 20$ dB, où $m_1$ est le moment d'ordre 1 de notre distribution. Ces différentes amplitudes sont représentatives d'une rupture temporelle de plus en plus forte qu'on pourrait avoir par exemple avec l'apparition d'une cible ponctuelle. Ces réalisations permettent alors de calculer un coefficient de variation unique, à la fois pour du speckle pur et pour une rupture temporelle. Cette estimation est faite pour tous les *N* variant de 2 à 200 sur la figure 3. Nous avons également affiché la courbe de valeur $\gamma + \sigma_{\gamma}$, pour avoir une idée de la marge statistique de distance entre le speckle pur et la rupture temporelle. On constate que le coefficient de variation du speckle pur converge en fonction de *N* vers la valeur théorique attendue (0,22), et pour les autres cibles, converge plus ou moins vite vers des valeurs qui sont fonction du rapport entre l'amplitude de la cible et celle du speckle.

La position du maximum entre le coefficient de variation de la population avec rupture et celle sans rupture, dépend à la fois de *N* et de les rapports d'amplitudes entre cible et speckle. Un résultat remarquable est que ces courbes ne dépendent pas de la valeur moyenne du speckle : **à nombre d'images fixé, en présence d'une rupture temporelle, le coefficient de variation ne dépend que du saut de rupture entre le speckle et la cible et du paramètre L.** Il s'écarte rapidement de la valeur du coefficient de variation d'un speckle stable.

<figure><img src="image_3.png" /><figcaption>**FIGURE 3:** Statistiques simulées du coefficient de variation sur un profil de type speckle "stable" + rupture temporelle ponctuelle.</figcaption></figure>

Ce résultat a été confirmé expérimentalement en reproduisant les simulations pour d'autres valeurs du paramètre $\mu$ : les courbes résultantes de la figure 3 restent totalement inchangées.

## 3.3. Choix de normalisation

La normalisation du paramètre de saturation à partir du coefficient de variation s'impose pour deux raisons. La première est liée à la nécessité d'obtenir des valeurs de saturation comprises entre 0 et 1, et qui s'adaptent à notre perception visuelle des couleurs. La seconde est liée au fait que le nombre d'acquisitions sur une scène varie selon le point considéré du globe terrestre, et que nous souhaitons nous affranchir de la dépendance du coefficient de variation au nombre d'images acquises.

Nous supposons que les pixels considérés pour décrire notre scène appartiennent soit à une zone de speckle stable, soit à une zone sujette aux changements. La distribution du coefficient de variation obtenu sur l'ensemble de nos pixels va donc être bimodale, comme représentée sur la figure 4 avec un premier mode qui correspond au speckle stable, dont les valeurs théoriques de moyenne $m_{\gamma}$ et d'écart-type $\sigma_{\gamma}$ ont été exprimées en fonction de *L* et *N* dans la section 3.1. Le comportement des pixels sujets au changement est bien moins décrit : nous savons simplement que les valeurs moyennes du coefficient augmentent rapidement avec l'amplitude du changement. Le but de notre normalisation est de distinguer les pixels de speckle stable et ceux de changement par la saturation de la couleur, en ramenant à des valeurs moyennes de saturation faibles tous les pixels appartenant à la première catégorie, et à des valeurs de saturation élevées et bien étalées pour ceux appartenant à la seconde.

Pour mettre en place la normalisation, nous agissons sur les paramètres connus de moyenne $m_{\gamma}$ et écart-type $\sigma_{\gamma}$ du coefficient de variation du speckle stable, pour ramener le premier mode du coefficient de variation autour d'une valeur moyenne de 0,25, et d'un écart-type de 0,1. Ces deux valeurs ont été choisies empiriquement : l'œil est peu sensible aux couleurs de saturation autour de
<figure><img src="image_4.png" /><figcaption>FIGURE 4: Normalisation empirique du coefficient de satura-
tion.</figcaption></figure>

0,25 et, en pratique, un écart-type de 0,1 pour ce premier
mode donne une visualisation satisfaisante. Pour simpli-
fier la mise en pratique de cette normalisation, nous pro-
cédons comme pour une loi normale. Même si la loi sta-
tistique du coefficient de variation d'un speckle stable est
supposée être du type Rayleigh Nakagami définie sur $\mathbb{R}^+$
et non une loi normale définie sur $\mathbb{R}$, cette approximation
semble en pratique peu pénalisante pour $L$ assez grand,
et nous permet d'avoir une forme de normalisation empi-
rique au premier ordre qui ne dépend que des nombres
$L$ et $N$. Ce choix nous conduit à l'expression suivante :

$$
\gamma \leftarrow 0,25 + 0,1 \times \left( \frac{\gamma - m_{\gamma}}{\sigma_{\gamma}} \right). \qquad (7)
$$

4. Portage du code sur la plateforme
Google Earth Engine

4.1. Présentation de la plateforme

Google Earth Engine est une plateforme en ligne dé-
diée aux traitements de données numériques d'observa-
tion terrestre (Gorelick et al., 2017). Son avantage est
de stocker les pétaoctets de données collectées par les
satellites et accessibles en source ouverte, notamment
ceux de la constellation Copernicus, ce qui évite leur té-
léchargement. Elle permet ensuite d'implémenter et de
tester des traitements au travers d'une interface de dé-
veloppement (API en javascript) et de visualiser les ré-
sultats à l'échelle de la planète, en bénéficiant de toute
la puissance de calcul Google.

Cependant, comme pour tout choix de plateforme, il
convient de prendre en compte les aspects négatifs de
celle-ci. L'interface de développement est contrainte par
le langage choisi et les fonctions mises à disposition en
javascript. Enfin, même si la plateforme GEE est encore
au stade d'outil de recherche et de fait est amenée à
connaître des évolutions, certains choix techniques ac-
tuels, notamment sur la préparation des données, (éta-
lonnage, interpolation, coregistration) peuvent avoir des
conséquences sur les données elles-même.

La méthode REACTIV a été implémentée sur la plata-
teforme pour être appliquée aux données Sentinel 1, afin
de démontrer l'intérêt de la visualisation sur toute la pla-
nète, et concomitamment de mieux cerner les limites de
la plateforme. Notre démarche a consisté à comparer sur
une même zone, les résultats obtenus sur la plateforme
avec ceux obtenus avec les jeux de données télécharg-
gés et traités localement. La zone d'étude qui a été rete-
nue est Saclay en Ile de France. Le téléchargement des
piles Sentinel 1 GRD s'est fait à partir de la plateforme
de l'Agence Spatiale Européenne. Nous avons alors pré-
traité les données par le logiciel SNAP, puis affiné le re-
calage des images avec l'algorithme GeFolki (Plyer et al.,
2015), et enfin appliqué la méthode REACTIV implémen-
tée en langage Matlab.

4.2. Résultats obtenus en local et sur la plateforme

La pile constituée porte sur 57 images de passes
descendantes, acquises entre le 29-01-2016 et le 04-04-
2017, avec une incidence moyenne constante. La figure
5 montre à la fois la représentation obtenue sur la plate-
teforme, et celle obtenue localement sur données téléchar-
gées, pour la polarisation VH.

Visuellement, les représentations ne semblent pas
différer. Elles laissent apparaître les parcelles cultivées
avec des couleurs fortement saturées, ainsi que locale-
ment, des zones de construction dont nous avons pu vé-
rifier qu'elles correspondaient effectivement à des chan-
tiers mis en place à cette période.

4.3. La qualité des dynamiques statistiques

Aujourd'hui, sur la plateforme, seules les données
Sentinel-1 géoréférencées de type GRD sont disponibles.
Les données SLC de Sentinel-1 ne peuvent pas être ac-
tuellement intégrées, car Earth Engine ne prend pas en
charge les images dont les valeurs sont complexes, en
raison de son incapacité à les calculer pendant la dé-
composition en pyramide des images sans perdre les in-
formations de phases. Chaque scène GRD incluse dans
GEE a été pré-traitée avec SNAP en procédant à la sup-
pression du bruit thermique, à la calibration radiométrique
et à la correction de terrain réalisée avec le MNT SRTM
30 ou ASTER DEM pour les zones supérieures à 60 de-
grés de latitude. Aucune correction radiométrique liée au
relief n'est aujourd'hui apportée. Les valeurs des images
corrigées géométriquement sont alors converties en dé-
cibels via la mise à l'échelle logarithmique.

Afin de gagner de l'espace mémoire, la plateforme
fait le choix de convertir les valeurs float32 en entiers non
signés de 2 octets (uint16), et ne conserve que les va-
leurs comprises entre le 1er et le 99ème centile des va-
leurs avant d'appliquer la quantification. Cela peut avoir
pour conséquence de modifier les statistiques d'origine
des signaux SAR.

L'étalonnage des données temporelles de la plate-
teforme et des données locales a été vérifié par la mé-
thode proposée dans (Quin et al., 2014). Elle n'a pas
apporté de gain substantiel à la stabilité de nos images,
<figure><img src="image_5.png" /><figcaption>A Gauche : Représentation obtenue dans GEE ; à droite, représentation obtenue sur données téléchargées en local</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**FIGURE 6:** Statistiques locales du coefficient de variation (CV) sur une zone de speckle dite "stable".</figcaption></figure>

et nous avons donc conservé l'étalonnage radiométrique basé sur les métadonnées.

Nous avons analysé les distributions statistiques du coefficient de variation sur différentes zones homogènes sélectionnées dans l'image. Notons que dans ce cas, ce qui nous permet de regarder la distribution de ce paramètre sont les échantillons spatiaux, alors même que le paramètre n'a été calculé que de manière temporelle.

L'histogramme du coefficient de variation estimé sur les données locales est représenté sur la figure 6, pour une zone qui semble appartenir à une zone de speckle "stable" (zone boisée). Dans ce cas, la distribution semble être bien modélisée de manière empirique par une loi de Rayleigh Nakagami, dont l'espérance et l'écart-type calculés par maximum de vraisemblance sont égaux respectivement à $m_γ = 0,21648$ et $σ_γ = 0,0214$, valeurs proches des valeurs théoriques trouvées en section 3.1, respectivement 0,2286 et 0,0212. Un test d'hypothèses avec la loi de Rayleigh Nakagami de ces paramètres conclut à l'adéquation de cette loi.

Les mêmes estimations menées sur GEE donnent comme moyenne 0,1992 et comme écart-type 0,0200, ce qui fournit des estimations légèrement biaisées par

rapport à celles des données traitées localement. L'analyse statistique a été reproduite pour trois zones de speckle soumis à changements, choisies dans des parcelles cultivées de différentes couleurs de la figure 5. Pour ces trois champs, les valeurs moyennes trouvées pour le coefficient de variation évoluent entre 0,3 et 0,4, et les écarts-types entre 0,03 et 0,05. Selon les tests d'hypothèses mis en place, les distributions obtenues ne suivent plus la loi de Rayleigh Nakagami.

En résumé : localement, l'analyse statistique montre que le coefficient de variation d'une zone sans changement suit bien une loi de Rayleigh Nakagami dont les paramètres sont proches des valeurs théoriques trouvées. Sur la plateforme, les résultats d'estimation sont légèrement biaisés, même si au premier ordre, les erreurs sur l'estimation du coefficient de variation sont négligeables : ce sont principalement les statistiques portant sur les dif-fuseurs permanents, dont les valeurs se trouvent dans la queue de distribution, qui seront affectées. Il faut donc garder à l'esprit que le choix actuel de compression de dynamique de la plateforme peut affecter des résultats quantitatifs basés sur les statistiques, même si sur le plan qualitatif, la visualisation n'est pas affectée.

## 4.4. La qualité du recalage

<figure><img src="image_7.png" /><figcaption>**FIGURE 7:** Zoom sur une représentation effectuée sur les données de la plateforme à gauche, et des données téléchargées puis recalées à droite.</figcaption></figure>

Que ce soit sur la plateforme ou en local, les données sont géoréférencées à l'aide de SNAP, ce qui as-
sure leur mise en correspondance. Cependant, la précision du résultat obtenue peut encore être améliorée. Sur les données téléchargées localement, une méthode de recalage subpixelique, capable des précisions requises pour le mode interférométrique, a été appliquée (Plyer et al., 2015). Elle amène à trouver des décalages d'une image à l'autre non négligeables, pouvant aller jusqu'à 2 pixels, avec un écart-type de distance trouvé de 0,4 pixels. Ces décalages résiduels peuvent avoir une influence non négligeable sur la saturation et la teinte obtenue dans la composition, notamment autour des éléments saillants : un décalage sur une cible fixe peut se traduire par un changement sur les bords de cette même cible. Parallèlement, la date de l'événement ponctuel trouvé peut également être modifiée. Ces effets sont illustrés sur la figure 7, sur un zoom effectué sur le site de Saclay, où plusieurs diffuseurs sont entourés d'un halo coloré.

## 5. Résultats de visualisation à travers le monde

Tout en gardant en tête les limites de la plateforme GEE, principalement liées aux imperfections de géoréférencement, celle-ci est intéressante pour tester notre approche à n'importe quel endroit de la planète et analyser très rapidement les résultats obtenus. Elle a ainsi permis de montrer l'intérêt de notre approche sur différents scénarios de détection de changements : le suivi d'activité agricole, le suivi de l'expansion urbaine et enfin la détection de bateaux.

### 5.1. Les pratiques agricoles

L'un des premiers constats relatif à la visualisation colorée est la saturation marquée des couleurs sur toutes les parcelles agricoles. Ceci conforte la grande sensibilité des données radar aux dynamiques présentes sur les champs, mises en lumières par la composition colorée REACTIV. Ce type de composition colorée peut ainsi contribuer à la classification des parcelles agricoles : une analyse rapide nous permet de voir que les deux compositions obtenues pour la polarisation VH et pour la polarisation VV apportent des renseignements complémentaires aux données multispectrales traditionnelles.

Par ailleurs, le coefficient de variation temporel pour des parcelles agricoles prend des valeurs bien plus élevées que pour des parcelles de forêts pour lesquelles les effets saisonniers sont bien moindres. Ce constat peut être mis à profit pour distinguer parmi les parcelles végétalisées, celles de type agricole des zones boisées.

### 5.2. La croissance du tissu urbain

Outre la présence de parcelles agricoles cultivées, la composition colorée permet de mettre en avant les constructions de bâtiments. Ceci a déjà été validé sur le site de Saclay. Un autre exemple est donné sur la figure 8 qui permet de visualiser l'évolution d'un quartier de Pékin (Beijing) entre février 2016 et février 2018. Les nombreux spots colorés apparaissant dans l'image témoignent des modifications fréquentes de l'occupation

des sols. A grande échelle, il devient possible de cartographier l'ensemble des zones urbaines qui subissent des développements à un rythme rapide.

De la même manière, il est envisageable d'étudier la destruction de structures. Sur la figure 9, le résultat de la représentation REACTIV appliquée à la ville de Palmyre en Syrie en est un exemple. Elle fait apparaître un grand nombre de modifications, principalement en vert et magenta : ces deux couleurs correspondent précisément aux dates de reprises de la ville par l'Etat Islamique, accompagnées de destructions importantes, en mai 2015 et décembre 2016.

### 5.3. Le trafic maritime

La visualisation REACTIV étant particulièrement sensible aux objets qui apparaissent de manière ponctuelle, elle permet de repérer de manière très rapide l'ensemble des bateaux de taille suffisante, dans toutes les zones portuaires et côtières. Sur la figure 10, on visualise notamment les produits obtenus autour de la zone portuaire de Singapour, pour laquelle le transport de commerce maritime est particulièrement important dans le monde. Aussi, la visualisation temporelle permet de mettre en avant à la fois les accumulations de bateaux sur les zones d'attentes, mais aussi les routes maritimes.

La composition colorée regardée à grande échelle sur la figure 11 fait apparaître clairement ces routes maritimes.

<figure><img src="image_8.png" /><figcaption>**FIGURE 8:** Est de Pékin (Beijing), Chine, sur deux ans de 2016 à 2018.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**FIGURE 9:** Composition REACTIV sur le site de Palmyre, entre février 2014 et février 2017.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**FIGURE 10:** Composition REACTIV sur Singapour de 2015 à fin 2017 à gauche.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**FIGURE 11:** Composition REACTIV sur la mer Méditerranée et Océan Atlantique, permettant de visualiser les routes maritimes.</figcaption></figure>

Il est également possible d'aller au delà d'une simple visualisation en proposant une carte binaire de présence des bateaux. Pour cela, on se propose de mettre en place un simple test statistique. Si l'on approxime la distribution du coefficient de variation comme une loi normale, alors on peut fixer un seuil de détection à la valeur 0,2286 + z₁-α.0, 1616/√N où N est le nombre d'images dans la pile temporelle et z₁-α est le quantile 1-α d'une loi Normale centrée réduite. Le fait de passer par une détection temporelle uniquement et non spatiale comme dans (Santamaria et al., 2017), fait de ce simple test statistique un outil particulièrement rapide par rapport à ces homologues.

Ainsi, sous GEE, la carte de détection moyennée spatialement sur la figure 12 a été obtenue pour l'ensemble de la période d'acquisition de Sentinel 1 (de fin 2014 à fin 2017) en moins d'une minute, permettant d'obtenir une carte pouvant s'interpréter directement en densité de trafic sur toute cette période.

## 6. Conclusions

Dans cet article, nous avons proposé une composition colorée REACTIV, basée sur un codage de l'information dans le plan HSV, afin d'optimiser l'analyse de grandes séries temporelles SAR et notamment d'adresser la problématique de détection de changements.

<figure><img src="image_12.png" /><figcaption>**FIGURE 12:** Carte de densité de bateaux calculée à partir du coefficient de variation temporel sous GEE.</figcaption></figure>

L'efficacité de cette approche repose sur l'utilisation du coefficient de variation temporel dont les propriétés statistiques ont été étudiées dans le cadre de réalisations d'une loi de Nakagami. Il en est ressorti que, pour une zone de speckle, ce coefficient a une valeur moyenne constante qui ne dépend que de L, le nombre de vues équivalent, et que sa variance décroît avec le nombre d'images. Une rupture temporelle dans un profil de speckle va se traduire par une augmentation de la valeur moyenne de ce coefficient. Grâce à l'approche REACTIV, cette propriété se traduit dans l'image par une couleur vive, significative de la présence d'un changement.

La méthode a été implémentée sous la plateforme Google Earth Engine, dont les limites principales sont la compression de dynamique et la qualité perfectible du recalage par défaut des données GRD, dont il faudrait mesurer l'impact. Enfin, à l'échelle mondiale, nous avons présenté des produits de visualisation pour différentes applications qui ont permis de montrer l'efficacité de la méthode.

Ces résultats prometteurs nous conduisent aujourd'hui à envisager de passer de la visualisation à la détection proprement dite des changements. Un autre défi sera de gérer le contexte multi-événements, c'est-à-dire, à disposer d'une capacité à détecter plusieurs changements dans le temps en un même endroit.

## Distribution du code

Le code REACTIV est accessible via la plateforme de partage GitHub à l'adresse https://earthengine.googlesource.com/users/ekoeniguer/reactiv. Les paramètres sont à renseigner en tout début de code : la polarisation, le type de passe (ascendante ou descendante), ainsi que la date de début et la date de fin de la plage temporelle observée. Si vous utilisez ce code dans vos projets, merci de citer cet article !

## Références

* Amitrano, D., Di Martino, G., Iodice, A., Riccio, D., Ruello, G., 2017. RGB SAR product exploiting multitemporal : general processing and applications. Dans : Analysis of Multitemporal Remote Sensing Images (MultiTemp), 2017 9th International Workshop on the. IEEE, pp. 1-4.
Boccardo, P., Gentile, V., Tonolo, F., Grandoni, D., Vassileva, M., 2015. Multitemporal SAR coherence analysis : Lava flow monitoring case study. Dans : Geoscience and Remote Sensing Symposium (IGARSS), 2015 IEEE International. IEEE, pp. 2699–2702.

Goodman, J., 1976. Some fundamental properties of speckle. JOSA 66 (11), 1145–1150.

Gorelick, N., Hancher, M., Dixon, M., Ilyushchenko, S., Thau, D., Moore, R., 2017. Google Earth Engine : Planetary-scale geospatial analysis for everyone. Remote Sensing of Environment 202, 18 – 27, big Remotely Sensed Data : tools, applications and experiences.

URL http://www.sciencedirect.com/science/article/pii/S0034425717302900

Kendall, M., Stuart, A., 1977. The advanced theory of statistics. Vol. 1 : Distribution theory. London : Griffin, 1977, 4th ed.

Koeniguer, E., 2017. Paris changes.

URL https://youtu.be/DJJmQltFB8o

Koeniguer, E., Boulch, A., Trouvé, P., Janez, F., 2018. Colored visualization of multitemporal SAR data for change detection : issues and methods. Dans : Proceedings of EUSAR 2018 : 12th European Conference on Synthetic Aperture Radar.

Nicolas, J.-M., 2006. Application de la transformée de Mellin : étude des lois statistiques de l'imagerie cohérente. Rapport technique, Telecom ParisTech.

URL https://perso.telecom-paristech.fr/nicolas/PUBLIS/nicolas2006D010.pdf

Nicolas, J.-M., 2017. Généralités sur la plateforme SARLAB.

URL https://perso.telecom-paristech.fr/nicolas/SARLABBETA/

Nicolas, J.-M., Tupin, F., 2016. Statistical models for SAR amplitude data : a unified vision through Mellin transform and Meijer functions. Dans : Signal Processing Conference (EUSIPCO), 2016 24th European. IEEE, pp. 518–522.

Plyer, A., Colin-Koeniguer, E., Weissgerber, F., 2015. A new coregistration algorithm for recent applications on urban SAR images. IEEE Geoscience and Remote Sensing Letters 12 (11), 2198–2202.

Quin, G., Pinel-Puyssegur, B., Nicolas, J.-M., Loreaux, P., 2014. MIMOSA : An automatic change detection method for SAR time series. IEEE Transactions on Geoscience and Remote Sensing 52 (9), 5349–5363.

Santamaria, C., Alvarez, M., Greidanus, H., Syrris, V., Soille, P., Argentieri, P., 2017. Mass Processing of Sentinel-1 Images for Maritime Surveillance. Remote Sensing 9 (7).

URL http://www.mdpi.com/2072-4292/9/7/678

Schulz, K., Cadario, E., Boldt, M., Thiele, A., 2010. Improving high-resolution repeat pass SAR image interpretation by the CoVAmCoh method. Dans : Proceedings of EUSAR 2010 : 8th European Conference on Synthetic Aperture Radar. pp. 1–4.

Sentinel 1 Team, 2013. Sentinel 1 User Handbook.

URL https://sentinel.esa.int/documents/247904/685163/Sentinel-1_User_Handbook

Wegmuller, U., Santoro, M., Werner, C., Cartus, O., 2015. On The Estimation and Interpretation of Sentinel-1 TOPS InSAR Coherence. Dans : Proceedings of FRINGE'15 : Advances in the Science and Applications of SAR Interferometry and Sentinel-1 InSAR Workshop, Frascati, Italy.