Parler à un RADAR spatioporté : traitement et analyse des données de Sentinel-1

Les RADAR sont des outils de télémesure exceptionnellement performants : contrairement aux observations optiques, ces méthodes actives de mesure ne sont pas dépendantes de l’illumination de la scène ni du couvert nuageux, puisque les ondes électromagnétiques radiofréquences traversent les nuages, permettant une observation, quelles que soient les conditions météorologiques. Il s’agit d’outils indispensables aux observations à long terme des impacts naturels et artificiels sur l’environnement ou la gestion de crises tel que promu par les Nations Unies [1], justifiant de diffuser la connaissance du traitement de ces signaux et de mettre en commun nos techniques d’analyse les plus efficaces. Nous avions pris rendez-vous à l’issue de [2] pour étendre nos connaissances acquises sur un RADAR à ouverture synthétique au sol (GB-SAR) vers les RADAR spatioportés : la promesse sera tenue avec le présent document.

Cependant, l’analyse des signaux bruts {I, Q} issus de ces mesures n’est pas complètement triviale, comme le serait l’analyse d’une image multispectrale dont on se contente de combiner les bandes spectrales pour identifier la nature de la cible réfléchissant la lumière du soleil. L’accès aux données acquises par les anciennes plateformes spatioportées [4] transportant des RADAR (ERS européen [5], RADARSAT canadien, TERRA-SAR et TanDEM-X allemands) nécessitait l’acceptation d’un programme scientifique ou technique d’analyse, qui ne pouvait donc être déposé que par des laboratoires ou instituts entraînés à traiter ce type de données, garantissant que des individus en dehors de ces agences ne solliciteraient pas les jeux de données. La politique de dissémination des données RADAR a changé en 2010, lorsque l’agence spatiale européenne ESA a décidé de libérer ces informations en les distribuant sur son site Copernicus sur http://scihub.copernicus.eu/dhus/. Ainsi, [3, p. 67] nous informe que depuis 2011, « Sentinel Data Policy: [...] anyone can access acquired Sentinel data; in particular, no distinction is made between public, commercial and scientific uses, or between European or non-European users […]. It is expected that this open and free access to the data (for any purpose, within or outside Europe), will maximise the beneficial use of data for the widest range of applications. It will strengthen EO markets in Europe, in particular the downstream sector, with a view to promoting growth and job creation. ».

Il devient maintenant possible d’accéder aux données brutes acquises par les plateformes européennes, et en particulier les actuels Sentinel-1A et B (figure 1) [3]. Ces deux satellites en orbite polaire basse, à une altitude d’environ 700 km, couvrent l’intégralité de la surface de la Terre, avec un passage sur un même site tous les 12 jours. Ainsi, les données acquises depuis 2014 sont accessibles gratuitement, permettant autant une analyse historique que contemporaine de l’évolution de la topographie, sous réserve de savoir analyser les mesures, ce qui est le sujet de ce document.

Après avoir rappelé les bases de l’analyse du RADAR à antenne synthétique (SAR), et les particularités des plateformes spatioportées, nous allons considérer trois cas d’étude qui incluent la localisation de l’auteur à Besançon, où les mouvements de terrain ne sont pas d’un grand intérêt ; une zone sismique et une zone d’exploitation minière en latitude tempérée couverte par le modèle numérique de terrain acquis par la navette spatiale lors de la Shuttle Radar Topography Mission (SRTM) en 2000 ; et finalement, une région arctique aux latitudes trop au nord pour avoir été couverte par SRTM, qui s’est limitée aux latitudes comprises entre 60° (contraintes orbitales du lancement de la navette spatiale), imposant l’utilisation d’un modèle numérique de terrain dédié.

1. Principes de base d’InSAR dans le cas des RADAR spatioportés

Nous avons présenté les généralités sur le traitement SAR [6, 7 8] et la mesure interférométrique résultante (InSAR) [2] auparavant. Pour reprendre rapidement, la distance est mesurée par compression d’impulsion (corrélation) dans la direction de visée avec une résolution inverse à la bande passante B, tandis que la résolution azimutale est obtenue par compression le long de la direction de déplacement de l’antenne synthétique – dans le cas du satellite, le long de son orbite. Nous avions vu que ces deux transformées de Fourier impliquent une période d’échantillonnage constante, tant dans le temps (compression en distance) que dans l’espace (compression en azimut). Nous apprenons sur les divers sites décrivant Sentinel [9] qu’il s’agit d’émetteurs radiofréquences en bande C centrés sur 5,405 GHz et occupant environ 100 MHz de bande passante. Tout comme les sonars sous-marins (side-scan sonar [10]), un RADAR spatioporté ne vise pas à sa verticale, mais avec un angle sur le côté afin de profiter des différents décalages Doppler pour permettre la compression d’impulsion : dans le cas de Sentinel, cet angle est d’environ 30 à 40° dans le mode de mesure nommé Interferometric Wide (IW) Swath tel que décrit sur https://sentinel.esa.int/web/sentinel/technical-guides/sentinel-1-sar/products-algorithms/level-1/single-look-complex/interferometric-wide-swath. Le traitement préliminaire de compression des impulsions en distance et en azimut a déjà été effectué dans les produits de niveau 1, tel que nous explique https://sentinel.esa.int/web/sentinel/missions/sentinel-1/data-products avec « Level-1 Single Look Complex (SLC) products consist of focused SAR data geo-referenced using orbit and attitude data from the satellite and provided in zero-Doppler slant-range geometry ». Ces données sont donc prêtes pour les traitements ultérieurs de mesure de réflectivité [11], cohérence [12] et phase [13] (InSAR) puisque les grandeurs résultantes restent complexes (I+jQ). La taille d’un pixel est de l’ordre de 5 par 20 m : inutile donc d’imaginer viser son jardin avec ce type de mesure, la surface analysée couvre une centaine de m². Par contre, la résolution en mesure de déplacement est excellente (figure 3) : avec une longueur d’onde de 300/5400=5,6 cm, la phase présente une rotation de 2π aller-retour tous les 2,8 cm au-delà de laquelle une incertitude subsiste sur le nombre de déplacements multiples de cette valeur, mais un déplacement de quelques millimètres à grande échelle est aisément détectable. On notera, et ce sera la première étape du traitement des données, que la position du satellite dans l’espace doit être connue à l’instant de la mesure à mieux que la résolution recherchée, une prouesse de l’observation des trajectoires de satellites et l’identification de leurs paramètres orbitaux qui permet aussi l’analyse centimétrique de signaux GPS a posteriori : cela explique qu’il faille attendre une vingtaine de jours entre l’acquisition des données et la possibilité de les traiter.

L’orbite polaire des deux satellites est conçue pour repasser « exactement » sur la même zone tous les 12 jours. Une analyse des passages au-dessus de la France indique que les Sentinel-1 passent en moyenne tous les 3 ou 4 jours : cependant, seuls les passages correspondant à une même orientation de l’orbite peuvent être traités en mode interférométrique. Cette information est obtenue dans le nom de fichier après la date d’acquisition : si nous considérons les passages de Sentinel-1A au-dessus de Besançon en novembre 2020, nous obtenons les données de la figure 3. Laurent Ferro-Famil (IETR, Rennes) nous informe que le traitement d’InSAR différentiel peut s’effectuer entre Sentinel-1A et Sentinel-1B, réduisant la période de répétition à 6 jours – nous n’avons pas ici expérimenté avec une telle approche.

Donc, nous pourrions croire, en l’absence du code couleur, que l’analyse est possible tous les 3-4 jours environ selon la nomenclature décrite sur https://sentinel.esa.int/web/sentinel/technical-guides/sentinel-1-sar/products-algorithms/level-1-product-formatting. Il faut prendre garde au numéro d’orbite qui suit la date (rouge) : seules les séquences commençant par 05513x (orange), 17320x (bleu) ou 17240x (vert) peuvent être comparées. Leur répétition est bien de 12 jours. De même, il n’est pas possible selon la séquence ici de mélanger des données de Sentinel-1A et B dans une analyse interférométrique. Les dates de passage de Sentinel sont prédites et affichées dans Google Earth avec les cartes disponibles sur https://sentinels.copernicus.eu/web/sentinel/missions/sentinel-1/observation-scenario/acquisition-segments.

2. Chaîne de traitement Sentinel

Sentinel fournit une multitude de services, mais celui qui couvre continuellement toute la Terre est SLC (Single Look Complex) dans le mode IW (Interferometric Wide) : dans ce mode, nous obtenons les cartes des coefficients {I, Q} desquelles nous pourrons déduire les interférogrammes suivant la documentation proposée par Sentinel-1 Toolbox – TOPS Interferometry Tutorial (janv. 2020) sur http://step.esa.int/docs/tutorials/S1TBX%20TOPSAR%20Interferometry%20with%20Sentinel-1%20Tutorial_v2.pdf. On pourrait s’interroger de l’intérêt d’écrire un article sur le sujet si cette excellente documentation est déjà disponible :

1. nous allons diverger quelque peu de la séquence proposée pour exploiter pleinement les fonctionnalités de Makefile sur lequel les auteurs du logiciel de traitement ne se sont pas appuyés et ;
2. la mise en œuvre des étapes de traitement ne s’est pas faite sans mal, compte tenu des messages d’erreurs laconiques et incohérents fournis par la machine virtuelle Java.

En effet, le logiciel proposé par l’ESA pour traiter les données de Sentinel, SNAP pour SeNtinel Application Platform (figure 4), se télécharge sur http://step.esa.int/main/download/snap-download/ en version binaire ou par ses sources sur https://github.com/senbox-org. Pour une fois, on s’autorisera l’exploitation des paquets binaires, la compilation manuelle semblant relativement pénible. Nous utiliserons la version 8.0.0 qui semble avoir résolu un certain nombre d’écueils sur la gestion de la mémoire, tout en restant relativement capricieuse.

2.1 Premiers pas

Avant de nous lancer dans la mesure interférométrique, assurons-nous déjà que nous sommes capables d’acquérir un jeu de données pertinent pour une région donnée, d’en extraire les informations utiles, et de les projeter sur une carte. Pour ce faire, nous nous intéresserons à Besançon, sans intérêt autre que de permettre une vérification sur le terrain des informations fournies par le RADAR. Nous suivons pour cela le tutoriel proposé sur [14] qui permet de prendre en main l’interface graphique du logiciel SNAP (figure 4), sans rencontrer de problème majeur autre que des erreurs aléatoires de la machine virtuelle Java dans laquelle s’exécute ce logiciel. L’incapacité à conclure de façon reproductible une séquence de traitement nous a encouragés à passer en ligne de commande pour éviter les opérations répétitives de l’interface graphique, libérer les ressources associées à l’affichage des résultats intermédiaires, et de façon générale réduire la perte de temps liée à la manipulation de la souris (section 3.2).

Malgré les erreurs aléatoires de dépassement du tas de mémoire (heap) de la machine virtuelle, que nous verrons être au moins en partie dues à la définition des zones sans-valeurs des modèles numériques de terrain utilisés, nous obtenons des résultats intermédiaires cohérents lorsque les champs (I, Q) ne sont pas tous nuls (figure 5), pour conclure le traitement avec des cartes de réflectivité (figure 6) cohérentes avec le contexte géographique. On note dès à présent que SNAP ne conserve que les coefficients I et Q issus de chaque traitement (et éventuellement le coefficient de corrélation lors du traitement interférométrique), mais que les grandeurs dérivées que sont la phase et le module, qui sont calculées au vol, doivent être explicitement mentionnées lors des dernières étapes de traitement pour se conclure par une carte de ces grandeurs.

La figure 6 présente clairement un fort réflecteur d’opportunité qu’est un bâtiment d’architecture moderne présentant de nombreux murs plans, métalliques et agencés à angle droit.

Le résultat du traitement est converti en format GeoTIFF compréhensible par QGIS par gdal_translate -of GTiff xxx_split_Orb_deb_TC.data/Intensity_IW3_VV.img TC_intensity.tif bien que Simon Gascoin (CESBIO, Toulouse) nous informe que QGIS sait lire les images au format .img issues de SNAP.

On pourrait se demander où est passé le satellite et quand, information primordiale si un événement bref doit être observé (par exemple, sortir un catadioptre au moment de son passage, voir la fin de ce document). Diverses sources d’informations sont disponibles pour prédire les passages de Sentinel-1 et ses paramètres de prise de vue. A posteriori, les fichiers .dim issus du traitement de SNAP contiennent ces informations. Par exemple, S1A_IW_SLC__1SDV_20200821T173205_20200821T173232_034008_03F26D_1274_split_Orb_deb.dim nous indique que :

La prise de vue était le 21 août à 17 h 32 (ces 9 secondes d’acquisition prennent plus de 4 h de traitement !). Connaissant la date et l’horaire de passage, nous pourrons demander à Heavens Above les paramètres du passage vu depuis le site de réception. Ce point sera détaillé plus tard (section 3).

2.2 Chaîne de traitement Sentinel par Makefile

Afin de limiter les risques d’erreurs lors d’opérations manuelles répétitives, voir d’automatiser la chaîne de traitement, les étapes documentées dans [14] sont automatisées :

1. parmi les mesures multiples acquises par Sentinel lors de son passage, ne conserver que les segments contenant la zone à observer : split ;
2. appliquer les observations orbitales pour positionner avec précision le satellite dans l’espace et ainsi identifier finement le temps de vol de l’onde électromagnétique, en connaissant la position de l’émetteur et du récepteur : orbit ;
3. fusionner les impulsions successives et éliminer les artefacts induits : deburst ;
4. projection dans le référentiel géographique, en choisissant la bonne zone UTM pour une projection plane adéquate, et correction des effets topographiques sur le positionnement des points de mesure : terrain correction ;
5. traduction en un format GeoTIFF lisible par QGIS.

Toutes ces opérations sont relativement bien documentées dans les tutoriels de l’ESA [14] et nous n’avons aucun intérêt à les reproduire ici. Le traitement manuel de ces opérations est rapidement long et fastidieux, surtout lorsque les temps de calcul sont de plusieurs minutes ou dizaines de minutes, interdisant de vaquer à d’autres occupations entre deux étapes. SNAP propose bien un mode de séquencement des tâches sous forme de graphique, mais tout utilisateur UNIX verra immédiatement la logique de passer par make lorsqu’une étape de traitement dépend du résultat des étapes précédentes. SNAP propose en effet en ligne de commande toutes les étapes de traitement mises en œuvre par l’interface graphique sous forme de l’instruction gpt [15]. Le nom du produit généré par chaque étape de traitement est imposé par l’option -t qui permet donc, quel que soit le nom du fichier d’entrée, de prédire la séquence des répertoires créés. Nous proposons donc le Makefile suivant avec les arguments que nous avons déduits de gpt <commande> --help :

Les arguments de TOPSAR-Split, et notamment le numéro du burst et des indices au sein de ce burst, nécessitent en l’état actuel de notre connaissance une évaluation, au moins sur la première image d’une séquence, par l’interface graphique SNAP. Le reste est totalement automatisé, pour prendre 90 à 180 secondes de calcul sur deux bursts sur un Panasonic CF-19 muni de 12 GB de RAM, mais surtout d’un stockage non volatile sur SSD, un point crucial pour le confort de travail lorsque plusieurs GB de données sont transférés de la zone de stockage en mémoire volatile pour traitement.

Plus compliquée (nous ne ferons pas le jeu de mots de nommer cette approche complexe), l’exploitation de la phase répète ces opérations sur deux images – une maîtresse et une esclave – afin de calculer un interférogramme qui nécessite un ajustement fin de l’alignement des deux images en vue de calculer leur écart de phase. Comme le même modèle numérique de terrain sera appliqué aux deux images, toute erreur sera annulée et seuls les effets atmosphériques [16], en plus de la différence de distance dans la direction de visée, impacteront la phase. Les étapes additionnelles à suivre pour un traitement interférométrique se résument, après les parties orbit (appliquées sur les deux images, maîtresse et esclave) à :

1. positionner avec précision l’image esclave par rapport à l’image de référence : back-geocoding ;
2. améliorer cet alignement : enhanced spectral diversity ;
3. calculer l’interférogramme (complexe) : interferogram ;
4. éliminer les artefacts d’acquisition des données : deburst ;
5. retrait de la phase topographique : topographic phase removal ;
6. moyenne spatiale : multilook ;
7. filtrage dans le domaine spectral spatial : Goldstein filtering ;
8. projection dans le référentiel géographique et correction des effets topographiques sur le positionnement des points de mesure : terrain correction ;
9. traduction en un format GeoTIFF lisible par QGIS.

Tout ceci, mis en œuvre dans un Makefile étendu, permet d’obtenir :

On notera en particulier que sur un ordinateur Panasonic CF-19 muni d’un processeur quad-core et de 12 GB de RAM, nous nous sommes longuement heurtés à des erreurs aléatoires de corruption de mémoire et de dépassement de taille du tas, avant d’inclure les arguments -c 2048M -q 2 qui limitent le nombre de threads exécutés en parallèle et la taille des blocs traités. Le traitement est un peu plus long, mais les plantages arbitraires de la machine virtuelle Java ne sont désormais plus qu’un mauvais souvenir, qui nous avait bloqués dans ce traitement depuis plusieurs mois.

La seule subtilité dans ce Makefile, outre l’identification des arguments de chaque commande, est que seules les bandes complexes I et Q sont conservées, tandis que le module et la phase du complexe sont déduits de ces grandeurs, mais ne sont pas conservés. Ainsi, lors de la dernière étape de génération du modèle projeté, obtenir la phase et le module (et non pas seulement la partie réelle et imaginaire du complexe, qui n’a aucun intérêt) nécessite d’informer explicitement gpt de notre volonté de générer ces bandes. Ce résultat s’obtient en analysant le nom des couches calculées par les étapes précédentes (disponible dans le répertoire data de chaque étape de calcul) et à coups de remplacement d’expressions régulières (sed), de générer le nom des bandes Phase* et Intensity*.

On notera que les cartes des modèles numériques de terrain sont copiées en local et stockées dans ~/.snap/auxdata/dem/ : cette information peut être utile si un serveur de modèle numérique devient inaccessible et un téléchargement manuel est nécessaire.

2.3 Application : exploitation minière à ciel ouvert

Nous allons appliquer cette chaîne de traitement sur un cas a priori simple d’exploitation de charbon dans la mine à ciel ouvert en Allemagne (figure 8) de Garzweiler près de Cologne et Düsseldorf (51°N, 6,5°E). Nous nous attendons à observer des variations de topographie de plusieurs centimètres à décimètres sur un environnement contrôlé dans l’intervalle de temps relativement court qu’est le taux de répétition des passages de 12 jours. Nous avons sélectionné deux passages les 15 et 27 juin 2019 pour nous affranchir des risques d’inactivité de l’été 2020.

L’application du dernier Makefile proposé doit se faire sans problème en renseignant convenablement FILENAME1 et FILENAME2 avec les jeux de données téléchargés – dans notre cas S1A_IW_SLC__1SDV_20190615T054201_20190615T054228_027686_032005_E4D9.zip et S1A_IW_SLC__1SDV_20190627T054202_20190627T054229_027861_03253C_9240.zip pour fournir le résultat proposé dans la figure 9.

Les mines – tout comme les autres zones urbanisées – sont parfaitement visibles par des niveaux de corrélation élevés, en l’absence de végétation qui absorbe et diffracte les signaux électromagnétiques, et sont caractérisées par des rotations rapides de phase, chaque rotation de 2π étant induite par une variation de distance de propagation de l’onde de λ/2 = 2,8 cm.

Moins simples à analyser, les rotations de phase lentes à grande échelle dans des zones portant de fortes corrélations ne sont pas faciles à interpréter. Afin d’analyser une hypothèse d’artefact atmosphérique [16], les images MODIS obtenues sur https://wvs.earthdata.nasa.gov/ pour les mêmes dates sont proposées : la seconde acquisition s’est faite sans couvert nuageux, mais la première présentait une couverture nuageuse dense qui aurait pu ralentir localement l’onde, compte tenu de la permittivité plus élevée de l’air humide (hypothèse qui mériterait une analyse plus poussée pour être validée [17, 18]). Nous n’avons pas réussi à identifier l’heure de prise de vue de MODIS Terra, tandis que MODIS Aqua présente pour la même date une couverture nuageuse très différente.

2.4 Déplacement induit par un tremblement de terre

Le tremblement de terre survenu le 7 janvier 2020 à Puerto Rico est décrit sur le site de la NASA sur https://www.nasa.gov/feature/jpl/puerto-rico-quake-damage-visible-from-space et va nous servir de référence pour valider notre compréhension de la chaîne de traitement sur un cas concret.

Bien entendu, ces mesures de déplacement, ici encore présentées sous forme de phase enroulée de 2π toutes les demi-longueurs d’onde, ne sont pas les ondes sismiques elles-mêmes, mais les déplacements dans la direction de visée du satellite résultant de leur propagation à la date du séisme. On notera que les populations vivant en zone sismique peuvent établir la distance à l’épicentre comme nous avons l’habitude de le faire avec les impacts de foudre. Tandis que le flash lumineux de la foudre arrive « instantanément » (0,3 km/µs), le son correspondant au tonnerre se traîne à 0,33 km/s (noter que l’unité de temps est passée de µs à s) et la distance à l’impact, en km, se déduit de la différence de retard entre les signaux lumineux et sonores divisée par 3. De la même façon, pour un tremblement de terre, l’onde longitudinale P se déplace à environ 5 km/s dans la roche tandis que l’onde de Rayleigh, la plus destructrice, se propage à 3 km/s et, entre les deux, l’onde transverse S, plus subtile à détecter, de l’ordre de 4 km/s. Un ordre de grandeur à l’épicentre en km est donc l’intervalle de temps entre les trois secousses en seconde.

2.5 Cas d’une zone arctique : modèle numérique de terrain dédié

Le dernier cas que nous proposons d’analyser – et qui en pratique nous intéresse le plus dans nos activités de recherche – est une zone arctique (figure 11). Ce cas est intéressant, car il n’est pas couvert par le modèle numérique de terrain issu de la navette spatiale SRTM et nécessite donc d’utiliser un modèle numérique de terrain dédié. Les modèles numériques de terrain globaux issus de traitements optiques (p. ex. GDEM) sont très mauvais dans les régions arctiques, où la couverture nuageuse est souvent importante et l’uniformité du manteau neigeux interdit un traitement efficace par corrélation des images ASTER [19] ou SPOT5 [20] sur les glaciers. Nous exploiterons donc un modèle numérique de terrain avec une résolution de 5 m disponible sur le site de l’Institut polaire norvégien https://geodata.npolar.no/ et en particulier [21]. La difficulté que nous avons rencontrée dans ce cas est de renseigner la valeur d’absence de données – par exemple au-dessus de la mer – information sans laquelle SNAP donne une erreur inextricable de dépassement de mémoire. Ainsi, partant d’un modèle numérique de terrain au format GeoTIFF dont les valeurs absentes sont probablement renseignées comme valant -∞, nous définissons ces valeurs comme nulles par gdal_calc.py -A DEM.tif --calc="A" --outfile=result.tif --NoDataValue=0 et informons gpt que la valeur nulle signifie une absence par -PexternalDEMNoDataValue=0.0. Cette solution résout une bonne partie des dépassements de mémoire de la machine virtuelle Java.

Il est intéressant de noter que l’interaction de l’onde électromagnétique micro-onde avec le milieu complexe qu’est le manteau neigeux et la glace fournit des informations non triviales sur la couverture neigeuse et son évolution au cours de la saison estivale. En effet, [11] nous enseigne que la ligne d’équilibre peut être déduite uniquement d’une analyse de réflectivité, sans entrer dans les subtilités de l’interférométrie. Une tentative de reproduction des résultats présentés dans cet article pour le contexte qui nous intéresse est fournie en figure 12, qui fournit bien une tendance correspondant au recul du manteau neigeux avec la saison de fonte, en accord avec une disparition de la neige au niveau de la mer à partir du 17 mai 2020 telle que documentée sur https://www.yr.no/en/statistics/graph/1-2837778/Norway/Svalbard/Svalbard/Ny-%C3%85lesund?q=2020-05. L’analyse de cette évolution n’est pas simple, entre eau liquide qui sature le manteau neigeux et augmente la réflectivité de l’interface avec l’air, rugosité de surface et changement d’interface entre neige et glace.

Un catadioptre a été déployé (figure 13), probablement pour une analyse des signaux rétrodiffusés par les satellites ERS au début des années 2000, près de Ny-Ålesund et pourra éventuellement servir de référence pour compenser les variations de propriétés atmosphériques lors d’analyses interférométriques à long terme. Cette cible coopérative n’est située qu’à 6 km environ de la zone qui nous intéresse (ROI) et devrait fournir une référence pertinente (voir section suivante 3).

3. Cible artificielle coopérative

Maintenant que nous sommes convaincus de traiter les données RADAR de Sentinel-1 et de les projeter dans un référentiel géographique cohérent, serions-nous capables de parler au satellite en rétrodiffusant la puissance émise vers le satellite ? Nous avons vu [2] qu’une géométrie appropriée est le catadioptre conçu pour renvoyer vers la source monostatique toute puissance rayonnée, et ce, quelle que soit la direction d’incidence. Le coin de cube doit cependant être suffisamment grand pour que la puissance qu’il réfléchit soit significative devant la puissance réfléchie par les obstacles présents dans le pixel de 5 m x 20 m illuminé par Sentinel-1. Alors que les zones planes ou végétalisées ne renvoient que peu de puissance, nous avons vu que les zones urbanisées peuvent renvoyer une puissance significative, surtout si le hasard de l’architecture présente des structures à angle droit. Les contraintes géométriques sur la réalisation du catadioptre sont drastiques pour réfléchir une fraction significative de la puissance incidente vers le récepteur situé à plus de 700 km dans l’espace, et une étude détaillée [22] montre que l’angle entre les faces ne doit différer de la valeur idéale de 90° que de quelques degrés et la planéité de quelques fractions de millimètres. Il aura fallu toute la dextérité de soudure du second auteur pour atteindre un résultat fonctionnel après une recherche de solution alliant rigidité, réflectivité de l’onde électromagnétique et poids sous forme des panneaux de polymères couverts d’aluminium Dibond (50 euros/m²). La dimension du catadioptre doit être estimée pour présenter une surface réflectrice équivalente importante devant les cibles non coopératives présentes dans le pixel contenant le coin réflecteur. Les deux géométries classiques sont le trièdre carré et le trièdre triangulaire : bien que la première géométrie présente un coefficient de réflexion nettement plus important, la seconde est presque toujours utilisée en pratique, probablement pour sa moins grande prise au vent, l’aisance de la protéger en tendant une toile entre ses coins, et une moins grande sensibilité du coefficient de réflexion à l’élévation. Dans nos conditions contrôlées, nous opterons pour le trièdre carré, il sera toujours possible de couper les faces carrées du cube en triangles isocèles plus tard.

La littérature [23] nous informe que nous devons viser un catadioptre de dimensions métriques. Concrètement, les trois cas qui nous intéressent sont :

• la sphère de rayon R présente une section RADAR RCS = πR² m² soit, en prenant R = 2 m pour modéliser une voiture, une valeur de l’ordre de 13 m² que nous devrons noyer dans le signal que nous rétrodiffuserons volontairement ;
• un catadioptre de côté triangulaire de longueur L illuminé par une onde de longueur d’onde λ présente une section RADAR RCS = (4/3)πL⁴/λ² soit pour un réflecteur de 1 m de long illuminé à 5,4 GHz une valeur de 1400 m² ;
• un catadioptre de côté carré de longueur L illuminé par une onde de longueur d’onde λ présente une section RADAR RCS = 12πL⁴/λ² soit pour un réflecteur de 1 m de long illuminé à 5,4 GHz une valeur de 12200 m² ou dix fois plus importante que le catadioptre triangulaire.

On notera que ces grandeurs de RCS évoluent comme la puissance quatrième de la longueur du côté L et chutent donc d’un facteur 16 en réduisant la taille du catadioptre d’un facteur 2.

Un tel dispositif est réalisé afin de valider notre capacité à détecter un signal dominé par le catadioptre qui devra servir à terme de mesure ponctuelle de déplacement par une analyse de la phase dans le pixel contenant ce signal dominé par la position du réflecteur uniquement, et non des autres cibles non coopératives (véhicules par exemple) qui seraient susceptibles de perturber la mesure.

Par ailleurs, la littérature nous informe que le maximum de réflectivité d’un catadioptre est atteint pour un angle de 35° par rapport à la plaque qui forme sa base. Connaissant l’élévation du satellite à son apogée, nous devrons donc pencher le coin réflecteur afin que de maximiser le coefficient de réflexion. Comment savoir a priori où et quand passera la satellite ? L’ESA publie des cartes au format KML (facilement exploitables dans Google Earth) avec les dates et horaires de passage ainsi que les zones couvertes. Ainsi, https://sentinel.esa.int/web/sentinel/missions/sentinel-1/observation-scenario/acquisition-segments fournit les zones observées lors de chaque passage. Mais cela ne nous dit pas où se trouve le satellite lors de son passage. Comme d’habitude, Heavens Above fournit la solution avec la position du satellite lors de son survol, en particulier pour savoir s’il observe la cible depuis l’est ou l’ouest, ce que la littérature nomme un passage ascendant ou descendant. Par ailleurs, Heavens Above nous informe de l’élévation au périgée : dans l’exemple de la Fig. 14, cette élévation est 47° pour une observation depuis Besançon : il faut donc pencher le catadioptre afin de l’élever de 47-35 = 12° par rapport à l’horizontale.

Une fois le catadioptre réalisé (figure 15), il est positionné sur un site où il ne devrait pas être dérangé et supposé stable, pour valider notre capacité à détecter sa présence sur les images Sentinel-1.

Afin de nous convaincre de la capacité à obtenir un signal significatif du catadioptre, nous avons traité pour les 4 passages possibles de chaque plateforme Sentinel-1 – deux ascendantes et deux descendantes – l’intensité des pixels autour de la zone où le signal du catadioptre doit se trouver. Pour ce faire, seul le premier Makefile présenté suffit, puisque nous n’effectuons pas (encore) un traitement interférométrique, mais uniquement une mesure de puissance rétrodiffusée. La figure 16 présente l’évolution de la carte de réflectivité en fonction des dates d’acquisition, pour 4 passages pour lesquels le catadioptre n’est pas encore déployé (haut), et un passage où le catadioptre est en place (bas, gauche). Le contexte géographique est présenté en bas, avec la flèche rouge indiquant l’emplacement prévu du catadioptre. L’évolution du maximum de réflectivité dans le carré jaune indiquant la position du catadioptre est résumée en figure 17 : une nette augmentation pour les passages des satellites illuminant le catadioptre (passages ascendants) est observée, alors qu’aucune évolution pour les passages n’illuminant pas la cible coopérative (orbites descendantes) n’est visible, car passant du mauvais côté (illumination depuis l’est alors que le catadioptre est orienté vers l’ouest).

Nous attribuons les deux zones de réflectivité moyenne (figure 16) qui polluent quelque peu nos résultats par les escaliers métalliques d’issue de secours le long du bâtiment se comportant comme réflecteurs RADAR : de tels artefacts ne devraient pas être visibles en dehors d’un environnement urbain. L’absence de réflectivité de l’escalier aussi présent sur le bâtiment du milieu n’est pas expliquée.

Conclusion

Suite à une étude qui nous avait permis d’appréhender les bases du traitement SAR et InSAR des RADAR au sol (GB-SAR), nous nous sommes efforcés d’appliquer ces connaissances aux données de RADAR spatioportés mises à disposition librement par l’ESA.

Partant de la bibliothèque SNAP de traitements fournie par l’ESA, nous avons automatisé la séquence dont chaque étape dépend du succès de la précédente par un Makefile. Le résultat, sous forme d’images géoréférencées projetées contenant corrélation, phase (figure 18) et réflectivité, s’intègre efficacement dans un logiciel de gestion d’informations spatialisées tel que QGIS pour une mise en contexte géographique des déplacements fins analysés par interférométrie radiofréquence. En complément de cette analyse des réflecteurs d’opportunité, nous avons déployé une cible coopérative – un catadioptre – dont le signal a été identifié de façon univoque après son installation. Ces résultats offrent des perspectives de suivi long terme de déplacements subcentimétriques liés à des évolutions géomorphologiques naturelles ou artificielles avec une résolution spatiale de l’ordre du décamètre. L’extension de la réalisation du catadioptre vers une solution démontable et transportable, résistante aux conditions climatiques arctiques, est en cours d’évaluation.

Références & Notes

[1] UN-SPIDER, « United Nations Platform for Space-based Information for Disaster Management and Emergency Response » et leur explication sur l’utilisation de Sentinel-1 : https://un-spider.org/advisory-support/recommended-practices/mudslides-flood-sentinel-1/step-by-step

[2] J. FRIEDT, W. FENG, « Résolution azimutale d’un RADAR à bruit : analyse et réalisation d’un RADAR à synthèse d’ouverture (SAR) par radio logicielle », GNU/Linux Magazine no242, novembre 2020 : https://connect.ed-diamond.com/GNU-Linux-Magazine/GLMF-242/Resolution-azimutale-d-un-RADAR-a-bruit-analyse-et-realisation-d-un-RADAR-a-synthese-d-ouverture-SAR-par-radio-logicielle

[3] ESA, « Sentinel-1 – ESA’s Radar Observatory Mission for GMES Operational Services », SP-1322/1, p.20, mars 2012 : https://sentinel.esa.int/documents/247904/349449/S1_SP-1322_1.pdf

[4] Nous empruntons ce qualificatif à Radars aéroportés et spatioportés (2015) à https://www.rncan.gc.ca/cartes-outils-publications/imagerie-satellitaire-photos-aer/tutoriels-sur-la-teledetection/teledetection-par-hyperfrequence/radars-aeroportes-et-spatioportes/9398

[5] ESA, « ENVISAT and ERS missions – Data access guide », 2011 : https://earth.esa.int/documents/10174/13019/Envisat_ERS_Data_Access_Guide.pdf

[6] J. FRIEDT, W. FENG, « Mesure fine de déplacement par RADAR interférométrique à synthèse d’ouverture (InSAR) par radio logicielle », GNU/Linux Magazine no244, janvier 2021 : https://connect.ed-diamond.com/GNU-Linux-Magazine/GLMF-244/Mesure-fine-de-deplacement-par-RADAR-interferometrique-a-synthese-d-ouverture-InSAR-par-radio-logicielle

[7] J.A. JOHANNESSEN, F. COLLARD, « SAR Instrument Principles and Processing », 2012 : https://earth.esa.int/eogateway/documents/20142/0/01_Wednesday_OCT2013_SAR_principles.pdf

[8] H. LEE, J. MOON, « Analysis of a Bistatic Ground-Based Synthetic Aperture Radar System and Indoor Experiments », MDPI Remote Sensing 13(1), 63, 2020 : https://www.mdpi.com/2072-4292/13/1/63/htm

[9] Sentinel dans la presse : https://www.bbc.com/news/science-environment-55022567 et https://www.bbc.com/news/science-environment-54926152

[10] D.W. HAWKINS, « Synthetic Aperture Imaging Algorithms:with application to wide bandwidth sonar », PhD, University of Canterbury, Christchurch, New Zealand, 1996 : https://ir.canterbury.ac.nz/handle/10092/1082

[11] S.H. WINSVOLD et al., « Using SAR satellite data time series for regional glacier mapping », The Cryosphere 12, 867–890, 2018.

[12] T. TAMM et al., « Relating Sentinel-1 Interferometric Coherence to Mowing Events on Grasslands », MDPI Remote Sensing 8 802, 2016.

[13] N.J SCHNEEVOIGT et al., « Glacier displacement on Comfortlessbreen, Svalbard, using 2-pass differential SAR interferometry (DInSAR) with a digital elevation model », Polar Record 48 (244), 17–25, 2012.

[14] A. BRAUN et L. VECI, « Sentinel-1 Toolbox – TOPS Interferometry Tutorial », janvier 2020 : http://step.esa.int/docs/tutorials/S1TBX%20TOPSAR%20Interferometry%20with%20Sentinel-1%20Tutorial_v2.pdf

[15] L. VECI, « SNAP Command Line Tutorial », ESA, 2020 : http://step.esa.int/docs/tutorials/SNAP_CommandLine_Tutorial.pdf

[16] M.P DOIN et al. « Corrections of stratified tropospheric delays in SAR interferometry: Validation with global atmospheric models », Journal of Applied Geophysics 69(1) 35–50, 2009 : https://www.sciencedirect.com/science/article/abs/pii/S0926985109000603?via

[17] I.P. KOVÁCS et al., « How to avoid false interpretations of Sentinel-1A TOPSAR interferometric data in landslide mapping? A case study: Recent landslides in Transdanubia, Hungary », Natural Hazards 96 (2) 693–712, 2019.

[18] P. MARINKOVIC et al., « InSAR quality control: Analysis of five years of corner reflector time series », Proc. Fringe 2007 (ESA SP-649) 26, 2007.

[19] T. TOUTIN, « ASTER DEMs for geomatic and geoscientific applications: a review », International Journal of Remote Sensing 29 (7) 1855–1875, 2008 : https://www.tandfonline.com/doi/full/10.1080/01431160701408477

[20] J. KORONA, É. BERTHIER et al., « SPIRIT. SPOT 5 stereoscopic survey of polar Ice: reference images and topographies during the fourth International Polar Year (2007–2009) », ISPRS Journal of Photogrammetry and Remote Sensing 64 204-–212, 2009.

[21] Norwegian Polar Institute, « Terrengmodell Svalbard (S0 Terrengmodell [Data set] », Norwegian Polar Institute, 2014 : https://doi.org/10.21334/npolar.2014.dce53a47

[22] M.C. GARTHWAITE et al., « The Design of Radar Corner Reflectors for the Australian Geophysical Observing System – A single design suitable for InSAR deformation monitoring and SAR calibration at multiple microwave frequency bands », Record 2015/03 avec la correction de l’équation erronée en p.27 dans R. SARATPULAPA et al, « Corner reflectors pattern study with & without orthogonality errors », 13th International Conference on Electromagnetic Interference and Compatibility (INCEMIC), pp.196–201, 2015.

[23] M. JAUVIN, et al., « Integration of Corner Reflectors for the Monitoring of Mountain Glacier Areas with Sentinel-1 Time Series », MDPI Remote Sensing 11 (8), 988, 2019.

[24] Y. QIN et al., « The design and experiments on corner reflectors for urban ground deformation monitoring in Hong Kong », Int. J. of Antennas and Propagation, 2013.