Cartographier l'occupation du sol avec Sentinel-2 : une approche en machine learning

Contexte

Les données satellitaires Sentinel-2 fournissent une observation multispectrale régulière de la surface terrestre à différentes résolutions spatiales. Ces données peuvent notamment être utilisées pour caractériser la végétation, les surfaces agricoles, les zones artificialisées ou encore les surfaces aquatiques.

Parallèlement, l'ESA WorldCover fournit une cartographie globale de l'occupation du sol à partir de données satellitaires et constitue ici une source de labels permettant d'entraîner un modèle d'apprentissage supervisé.

Ce projet explore donc la possibilité d'apprendre la relation entre les signatures spectrales observées par Sentinel-2 et les classes d'occupation du sol proposées par WorldCover.

Sur la tuile utilisée dans cette étude, les données ESA World Cover présentent les classes suivantes :

Etape 1 : Pre-processing des données

Les données Land Cover qui serviront de référence (labels) pour cette étude sont récupérées sur le site de l’ESA : viewer.esa-worldcover.org. La tuile correspondant à la région de Toulouse (France, Occitanie) est téléchargée depuis le portail. Ce produit fait référence pour l’année 2021.

Visualisation de ESA WORLD COVER

Pour les données d’entraînement, on se tourne vers le portail Copernicus qui permet de télécharger les données Sentinel-2 L2A pour une zone et temps donnés. Par souci de cohérence avec les données ESA World Cover, on récupère une date en 2021, tout en cherchant à optimiser la couverture nuageuse. C’est la date du 15 avril 2021 qui présente la meilleure pertinence, avec un Cloud Cover très faible (2.8).

Visualisation des données Sentinel en RGB

Les données Sentinel-2 L2A comprennent l’ensemble des bandes de 01 à 12 (ainsi que AOT, SCL, TCI, WVP) à des résolutions spatiales allant de 10m, 20m ou 60m selon les bandes considérées. Les données à 10m étant limitées à quelques bandes et les données ESA World Cover se trouvant à une résolution spatiale de 10m, le choix est fait de travailler à 20m, en appliquant un rééchantillonnage sur les données ESA pour les aligner parfaitement sur la grille Sentinel-2.

Ce rééchantillonnage est réalisé en parcourant le raster en fenêtre (grille) de 4 pixels (20m de résolution spatiale) et en sélectionnant la valeur majoritaire. Un score de purity est calculé simultanément afin d’identifier les pixels ambigus (pas de classe majoritaire) et de les écarter ensuite lors de l’entraînement des modèles.

A l’issue de cette étape, un raster de labels ESA World Cover à 20m de résolution est généré et les rasters Sentinel-2 L2A sont prêts pour commencer l’approche machine learning.

Etape 2 : Sélection des Area of Interest (AOI)

La tuile World Cover entière ne peut être utilisée comme base pour l’entraînement des modèles, d’une part à cause du trop grand nombre de pixels qui alourdiraient les temps de calcul au-delà du raisonnable, d’autres part en raison de la forcément très inégale répartition des pixels qui pourrait générer des biais d’entraînement.

Le choix est donc fait de travailler sur des zones de 20x20km (soit 1 000 000 de pixels), ce qui représente une taille à la fois exploitable et suffisamment riche en information. Néanmoins, là où une approche standard en machine learning consiste à appliquer une séparation randomisée des données en un jeu d’entraînement, un jeu de validation et un jeu de test, la problématique est légèrement différente dans le domaine de l’apprentissage géospatial. En effet, les risques de biais liés à l'autocorrélation spatiale des données nécessitent une totale indépendance géographique entre les données d’entraînement, de validation et de test.

On cherche donc à identifier un triplet d’AOI (train, validation, test) en optimisant leurs caractéristiques. A l’œil, les zones autour des communes petites à moyennes semblent présenter une bonne diversité de classes, mais on laisse le soin à un script Python de confirmer et sélectionner les meilleures zones. Pour cela, on parcourt le raster en grille avec calcul de plusieurs paramètres : Entropie de Shannon, Pixels valides (purity > 0.75), nombre de classes et équilibre. Une distance minimale de 20km est imposée entre chaque zone pour garantir l'indépendance spatiale.

A l’issue de cette étape, le meilleur triplet est déterminé ainsi :

Localisation des zones Train, Validation & Test

Un script Python est ensuite appliqué pour extraire les données ESA ainsi que Sentinel-2 pour ces 3 zones, au format .parquet bien adapté au traitement de ce volume de données.

Etape 3 : Tests exploratoires

L’état de l’art et les connaissances métier dans le domaine de la télédétection / Machine Learning géospatial nous poussent à sélectionner 9 bandes de Sentinel-2 pour une première approche exploratoire :

Un script python est utilisé pour lire les .parquet générés précédemment. On y filtre les pixels sur un critère de Purity > 0.75 (ce qui élimine environ 10% des pixels) et on sélectionne un seuil maximal de 10 000 critères par classe. On configure ensuite un ensemble de 3 modèles (Scikit-Learn) :

Les critères retenus pour l’évaluation sont l’accuracy, la balanced accuracy et surtout le F1-score (moyenne des scores F1 de chaque classe) qui permet de prendre en compte la précision et le recall et s’avère plus adapté aux problèmes présentant des classes déséquilibrées.

A l’issue de ce test avec un paramétrage standard des modèles, les résultats obtenus sur l’ensemble de validation sont les suivants :

Le DummyClassifier fournit une référence volontairement naïve, avec un macro-F1 de seulement 0,068. Les modèles basés sur des arbres améliorent très nettement ce résultat : le Random Forest atteint un macro-F1 de 0,496.

Cette première expérience constitue donc un premier élément en faveur de l'hypothèse de départ : les bandes spectrales Sentinel-2 contiennent suffisamment d'information pour apprendre une discrimination non triviale entre les classes étudiées.

Etape 4 : Feature Engineering et Ablation Study

L’étape suivante consiste à optimiser les résultats obtenus en influant sur les features d’entrée (prédicteurs). On calcule un ensemble d’indices spectraux classiques de télédétection, que l’on ajoute à notre feature set de prédicteurs :

Ces nouvelles features sont ensuite évaluées en entraînant les modèles précédents sur un grand ensemble de FeatureSets :

Par définition, les indices type NDVI sont basés sur des informations déjà contenues dans les bandes Sentinel-2. Le fait de les tester ensemble et séparément sert à vérifier s’ils apportent tout de même une information discriminante par rapport aux bandes et si oui, quels indices sont les plus pertinents.

L’ajout des indices spectraux semble au final apporter peu d’amélioration mesurable par rapport aux seules 9 bandes dans le cadre de notre Random Forest, ce qui s’explique sans doute par le fait que le modèle est capable de retirer l’essentiel de l’information contenue dans ces indices depuis les bandes brutes. Pour la suite de l’exploration, on retient néanmoins 3 feature sets :

On cherche ensuite à intégrer la dimension spatiale du problème. Pour intégrer cette information, des statistiques locales sont calculées sur des fenêtres carrées centrées sur chaque pixel :

Deux tailles de voisinage sont testées : 3×3 pixels (60 m × 60 m) et 5×5 pixels (100 m × 100 m), calculées sur B04, B8A et B11.

L'ajout de statistiques spatiales produit une amélioration plus nette que celle observée avec les indices spectraux. Les datasets retenus pour la suite sont :

Cette étape nous aura permis d’évaluer la pertinence et le poids des principaux prédicteurs utilisables et dérivables à partir des données Sentinel-2. Le travail de comparaison et d’ablation permet d’éliminer certains indices potentiellement redondants et de déterminer 3 ensembles de features de référence à tester avec différents modèles de classification.

Etape 5 : Comparaison avec des modèles de Gradient Boosting

Il est maintenant temps de procéder à la comparaison de performance entre différents modèles. Si Random Forest est une référence reconnue dans la littérature, d’autres modèles de Classifier – souvent plus récents – peuvent constituer des alternatives de choix. Au regard de l’état de l’art sur le sujet, 3 modèles de Gradient Boosting apparaissent comme des candidats sérieux : XGB Classifier, Light GBM et le GradientBoosting Classifier de scikit-learn

Les modèles sont tout d’abord entraînés avec des paramètres standards et leurs résultats sont évalués selon leur score macro-F1, leur accuracy, balanced accuracy mais aussi temps de traitement.

# Chargement du script Python en cours...

En sortie, nous obtenons :

Deux configurations apparaissent comme les meilleures candidates pour la phase d'optimisation : LightGBM avec l'ensemble complet de variables, et XGBoost avec les variables spectrales et spatiales. GradientBoosting ScikitLearn offre aussi une performance solide à priori, mais son temps d’entraînement très largement supérieur pour des résultats légèrement en deça en font un moins bon candidat pour un tuning fin.

Etape 6 : Tuning des modèles finalistes

Les deux modèles sélectionnés font ensuite l’objet d’une recherche des paramètres optimaux via RandomizedSearchCV.

# Chargement du script Python en cours...

Dans ce cadre, une cross-validation stratifiée en 3 folds est utilisée sur le jeu d’entraînement, avec le macro F-1 comme métrique à optimiser pour accorder un poids similaire à chacune des classes. Cette approche nécessaire à la recherche des meilleurs hyperparamètres reste cependant potentiellement sensible à l’autocorrélation spatiale et ne remplace pas la validation spatialement indépendante mise en place précédemment.

Ainsi le tableau ci-dessus démontre que les performances très élevées obtenues en cross-validation 3 folds ne sont pas reproduits sur le dataset de validation, ce qui confirme en outre l’importance de conserver un jeu de validation géographiquement indépendant.

Etape 7 : Analyse des erreurs

Après le tuning, les modèles LightGBM et XGBoost sont évalués sur le jeu de validation spatialement indépendant. L’analyse ne se limite pas à l’accuracy globale : les métriques sont étudiées classe par classe afin de mettre en évidence les erreurs de classification et l’impact du déséquilibre entre les classes.

Les deux modèles présentent des performances très proches sur les principales classes. Le Tree cover, le Grassland et le Cropland obtiennent des F1 compris entre 0,79 et 0,88, tandis que la classe Permanent water est également bien reconnue.

Une difficulté importante apparaît cependant pour la classe Built-up. Son rappel élevé (~0,94) indique que la majorité des pixels bâtis sont correctement détectés, mais sa précision très faible (~0,28) montre que de nombreux pixels appartenant à d'autres classes sont également prédits comme bâtis. Le modèle a donc tendance à sur-prédire cette classe.

Les classes Bare / sparse vegetation et Herbaceous wetland restent particulièrement difficiles à modéliser. Cette observation doit toutefois être interprétée avec prudence : elles disposent respectivement de seulement 24 et 80 observations dans le jeu de validation. Les métriques associées sont donc beaucoup moins stables que celles des classes majoritaires.

×

Les matrices de confusion permettent de préciser la nature de ces erreurs. Les principales confusions concernent les classes Tree cover, Grassland et Cropland. Ces classes correspondent toutes à des couvertures végétales dont les signatures spectrales peuvent être proches, notamment lorsqu'elles sont observées à une résolution de 20 m.

La classe Herbaceous wetland est également difficile à distinguer, notamment vis-à-vis des surfaces en eau. Cette difficulté peut s'expliquer par la proximité spectrale entre certaines zones humides et les pixels contenant de l'eau ou de la végétation humide.

Les différences entre les deux modèles restent faibles. XGBoost obtient une légère avance sur l'accuracy et le macro-F1, tandis que les deux modèles présentent une balanced accuracy très proche. Aucun modèle ne se distingue donc nettement sur l'ensemble des critères.

Etape 8 : Généralisation sur la zone de test

Nous procédons maintenant à l’application des modèles LGBM et XGBoost sur la zone de test, qui contient uniquement des pixels totalement indépendants de tout entraînement et validation réalisé jusqu’à maintenant. Les zones de train et validations sont fusionnées en un seul dataset qui va servir à entraîner les modèles avec les paramètres optimaux déterminés lors des étapes précédentes.

Visualisation des résultats

Les performances quantitatives sont maintenant complétées par une analyse spatiale des prédictions. L'outil ci-dessous permet de visualiser les prédictions de LightGBM et XGBoost au format raster et de les comparer à la référence ESA World Cover. Il propose aussi d'afficher les rasters d'erreurs de classification pour les 2 modèles. La présence de pixels noirs (no data) sur les prédictions s'explique par les pixels exclus sur le cirtère de purity (pixels dont le label n'est pas considré comme assez fiable suite au rééchantillonnage à 20m).

Les deux modèles présentent une organisation spatiale globalement similaire, ce qui est cohérent avec leurs performances très proches sur la zone de test. Une tendance déjà observée dans l'analyse des erreurs apparaît également visuellement : les deux modèles semblent sur-prédire les surfaces de type Built-up.

Gauche : ESA World Cover (Référence) Prédiction LightGBM Prédiction XGBoost Erreurs LightGBM Erreurs XGBoost

Droite : ESA World Cover (Référence) Prédiction LightGBM Prédiction XGBoost Erreurs LightGBM Erreurs XGBoost

--- ---

Légende Occupation du Sol (ESA WorldCover) :

Tree Cover

Grassland

Cropland

Built-up

Bare/sparse vegetation

Permanent water

Herbaceous wetland

No data (purity filter)

Légende Matrice d'Erreur :

Prédiction Correcte

Erreur de Classification

No data (purity filter)