Apport de la photogrammétrie satellite pour la
modélisation du manteau neigeux

Chaîne de production de cartes de hauteur de neige (SMASH)

 Calcul de MNE à partir d’images stéréo Pléiades 

Le calcul des MNE est effectué à l’aide du logiciel Ames Stereo Pipeline version 2.6.2 (Shean et al., 2016 ; Beyer et al., 2018), une librairie d’utilitaires pour le traitement d’images 30 Figure 2.1 – Position des sites d’études sur lesquels des cartes de HTN ont été calculées à partir d’images stéréo Pléiades traitées avec la chaîne SMASH. Toutes les cartes du manuscrit sont orientées de sorte à avoir le nord vers le haut. stéréo. Dans le cas d’une paire d’images stéréo, une des images de la paire est choisie comme référence. La position à la surface de la Terre (X,Y,Z) de chaque pixel de l’image de référence est calculée, générant un nuage de points. L’obtention de la position d’un pixel à la surface est possible d’une part grâce à l’identification du pixel homologue dans l’autre image de la paire et d’autre part grâce à la connaissance de la position et de la direction de visée du capteur lors de l’acquisition de chaque image (Figure 1.1). ASP réalise ces calculs en plusieurs étapes. Tout d’abord, l’orientation des images est ajustée à l’aide de l’identification automatique de points d’attachements dans les deux images (affineepipolar). Un ajustement des images par bundle adjust pourrait précéder cette étape. Le bundle adjust est conseillé par ASP et est de plus en plus utilisé pour le calcul de séries temporelles de MNE (Shean et al., 2020 ; Eberhard et al., 2021 ; Bhushan et al., 2021). Il permet d’assurer une meilleure cohérence entre les images d’une série, de rendre interprétable l’erreur d’intersection et d’éviter la coregistration pour une série de MNE ajustés ensemble. Le bundle adjust n’a pas été utilisé ici par simplicité et car une étape de co-registration est nécessaire en cas d’utilisation de MNE non photogrammétrique dans la série, comme dans ce cas le MNE lidar de l’ASO. La localisation initiale des images est connue grâce aux fichiers Rational Polynomial Coefficient 31 Chapitre 2. Cartographie de la hauteur de neige avec des images Pléiades Figure 2.2 – La chaîne de production de cartes de hauteur de neige à partir d’images stéréo Pléiades (SMASH). 32 (RPC) fournis avec les images. Les RPC sont les coefficients d’une fonction polynomiale qui relie les coordonnées images d’un pixel (i,j) avec ses coordonnées sur la surface de la terre (latitude,longitude, élévation). A cause de la distorsion des images par le relief, un même point de la surface correspond à un pixel différent dans chaque image de la paire. Les paires de pixels homologues (i.e. représentant la même portion de terrain) sont déterminées par de la corrélation d’image. Au cours de cette étape, des sous-portions des deux images sont comparées à l’aide d’une fonction de coût qui mesure leur similitude. Cela produit une carte de disparité dont la valeur indique la distance en pixels entre un pixel dans l’image de référence et son homologue dans l’autre image. Le nuage de points est calculé par triangulation des droites passant par les deux pixels homologues. Le nuage de points est converti en un MNE sur une grille régulière par un moyennage des points contenus dans chaque pixel de la grille. Les images panchromatique et multispectrale sont projetées sur ce MNE. Le calcul du MNE à la résolution finale est réalisé en deux itérations. La première itération produit un MNE grossier à 50 m à partir des images panchromatiques en géométrie capteur parfait. La deuxième itération produit un MNE fin à 3 m à partir des images panchromatiques projetées sur le MNE grossier. Ces itérations réduisent les artefacts (Shean et al., 2016 ; Beyer et al., 2018). ASP permet également de traiter des triplets d’images en réalisant une triangulation conjointe de deux cartes de disparité calculées à partir de deux paires du triplet. L’ordre des images du triplet importe puisque la première sera utilisée pour le calcul des deux cartes de disparité. Pour plus de détails sur le processus de restitution du relief, le lecteur pourra consulter le chapitre consacré dans Delvit et Michel (2016). De nombreuses options sont disponibles dans ASP pour adapter les calculs à la radiométrie des images (fort contraste ou non) et à la géométrie du terrain (forêt, environnement urbain, champ de neige). Seules celles utilisées dans la suite des travaux sont présentées ici. L’étape de corrélation des images qui met en relation les pixels homologues de la paire d’images est centrale et a deux variantes principales : les méthodes locales ou globales. Le corrélateur d’ASP par défaut utilise une méthode locale (local-search). Avec cette technique, la disparité est déterminée indépendamment pour chaque pixel. La disparité est mesurée 33 Chapitre 2. Cartographie de la hauteur de neige avec des images Pléiades en optimisant la fonction de coût par une exploration systématique du voisinage du pixel (block-matching). Plusieurs méthodes globales basées sur la méthode du Semi-Global Matching (SGM Hirschmuller, 2005) sont implémentées dans ASP. SGM minimise l’énergie de la carte de disparité en tenant compte de la valeur de la disparité des pixels voisins (Delvit et Michel, 2016) et en pénalisant les discontinuités. SGM est censé fournir de meilleurs résultats que le block-matching dans les zones de l’image à faible texture radiométrique (ASP user guide, https ://stereopipeline.readthedocs.io/en/latest/, dernière consultation le 1er septembre 2020). Il est également possible de résoudre des détails plus fins avec SGM que le local-search du fait de la plus petite taille des fenêtres de corrélation. Plusieurs fonctions de coût peuvent être utilisées. La fonction de coût utilisée par défault avec la méthode locale est la corrélation croisée normalisée (normalized cross correlation). Cette fonction est plus robuste aux changements de contraste entre les images de la paire que la somme des différences absolues ou que la somme au carré (Delvit et Michel, 2016). Pour les méthodes globales, les fonctions de coût dites « non-paramétriques » qui comparent les valeurs relatives des pixels sont conseillées (Zabih et Woodfill, 1994). ASP dispose pour cela des fonctions binary census transform et ternary census transform. Ces fonctions calculent un nombre binaire pour chaque pixel par comparaison de sa valeur avec la valeur des pixels de son voisinage. La comparaison avec les 8 pixels voisins produit un nombre à huit chiffres. Pour la binary census transform les chiffres valent 1 si la valeur du pixel voisin est supérieure à celle du pixel central, 0 sinon. Pour la ternary census transform, les digits peuvent prendre la valeur 00, 01, 11, selon que leur valeur est inférieure, comprise dans ou supérieure à un intervalle centré sur la valeur du pixel central. Les valeurs binaires des pixels des images de la paire sont ensuite comparées en comptant le nombre de digits différents (distance de Hamming) (Zabih et Woodfill, 1994). 

 Calcul de carte d’occupation du sol à partir des images multispectrales 

Les images multispectrales projetées sur le MNE à 3 m sont utilisées pour calculer des cartes d’occupation du sol à l’aide de la librairie Orfeo Toolbox (Grizonnet et al., 2017). L’algorithme de random tree forest est utilisé pour classer chaque pixel dans les catégories : neige, forêt, terrain stable ou eau. Les zones de forêt comprennent tout ce qui est interprété comme de la végétation haute. Le terrain stable correspond aux zones dans lesquelles la topographie est stable au cours du temps. Cela inclut les zones de terrains nus (roche, sol), de végétation basse (prairie) et les routes. Le jeu de données d’apprentissage est créé manuellement par interprétation de l’image panchromatique et de l’image multispectrale (Rouge, Vert, Bleu, PRoche Infrarouge) augmentées d’une bande NDVI ( P RI−R P RI+R ). Le jeu de données d’apprentissage est un ensemble de polygones répartis sur l’image à partir duquel un classifieur (un modèle statistique) va apprendre à classer les pixels de l’image multispectrale (+NDVI). Le classifieur est ensuite appliqué à l’ensemble de l’image. Un classifieur est calculé pour chaque date d’acquisition. Le calcul d’un unique classifieur pour l’ensemble de la série temporelle ne fonctionnerait probablement à cause de l’absence de correction atmosphérique, des différences d’illuminations et des différences de développement de la végétation. En l’absence de données de validation indépendantes les cartes d’occupation du sol sont simplement inspectées visuellement. Les cartes d’occupation du sol sont utilisées pour calculer des masques de neige et de terrain stable. Ces masques subissent une érosion morphologique pour éliminer les bordures des zones. Les bordures sont possiblement moins précisément classées du fait de la nature mixte des pixels dans ces zones. Pour la même raison, les îlots de moins de 30 pixels sont éliminés des masques. Ce choix est guidé par l’objectif d’obtenir un terrain stable de bonne qualité quitte à réduire sa surface. 

 Calcul des cartes de hauteur de neige à partir des MNE 

Les MNE Pléiades sont géo-localisés grâce à la connaissance de l’orbite et de la direction de visée du satellite au moment de l’acquisition renseignée dans des fichiers RPC. Cette géolocalisation a une précision de quelques mètres (Lebègue et al., 2012) et est susceptible de générer des différences d’élévations entre les MNE indépendantes des changements d’élévation réels. Il est donc nécessaire d’améliorer par un recalage la localisation relative des MNE hiver et été avant de les différencier. La méthode proposée par Nuth et Kääb (2011) a été implémentée dans une fonction python qui produit un MNE recalé à partir d’un MNE mobile considéré comme mal localisé et d’un MNE de référence (immobile). Le MNE été est choisi comme référence. Il n’est pas modifié lors de l’étape de recalage et le MNE hiver est recalé et interpolé sur sa grille (fonction python scipy.interpolate.RectBivariateSpline). Le vecteur de recalage est calculé par optimisation de la relation entre la différence d’élévation, l’orientation et la composante horizontale du vecteur de recalage (Nuth et Kääb, 2011). Le calcul est effectué sur les zones de terrain stable identifiées grâce à la carte d’occupation du sol Pléiades. Les pentes inférieures à 4◦ sont exclues pour le calcul du vecteur horizontal car le calcul de l’orientation est entaché d’erreurs sur les faibles pentes. Le calcul du vecteur de recalage est itératif et s’arrête lorsque la dispersion de la différence d’élévation n’est plus améliorée (i.e. réduite) ou lorsque le vecteur n’est pas modifié par plus de 0.1 m. La composante verticale du vecteur est calculée comme la moyenne de la différence d’élévation résiduelle après filtrage à 3 Normalized Median Absolute Deviation (NMAD) du résidu. La NMAD est une mesure robuste de la dispersion d’une population, X, marquée par la présence de valeurs aberrantes et est de ce fait adaptée aux résidus de différence d’élévation (Höhle et Höhle, 2009) : NMAD(X) = 1, 4826 × mediane(|X − mediane(X)|) (2.1) 

 Ressource informatique 

Les ressources informatiques requises pour le calcul d’une carte de HTN dépendent de la taille de la zone, de l’utilisation d’une paire ou d’un triplet d’images et des options de traitement. L’étape la plus longue est le calcul des MNE par ASP (fonctions ASP stereo et point2dem). Les calculs des MNE ont été effectués sur le supercalculateur HAL du CNES en utilisant 18 CPU. L’option SGM dans sa configuration optimale nécessite plus de mémoire (∼17 Go) que l’option Local-search (∼9 Go). Le temps de calcul est relativement proche pour les deux méthodes et dépend surtout de la taille de la zone (4h pour 90 km2 , 9h pour 220 km2 ) et de l’utilisation d’une paire d’image (6h) ou d’un triplet (9h). Le reste de la chaîne (recalage des MNE, calcul des cartes d’occupation du sol) prend environ 2h pour l’ensemble de la zone (∼ 300 km2 ). Des gains de temps seraient obtenus en tirant un meilleur parti des possibilités de parallélisation offertes par ASP et l’Orfeo Toolbox (OTB). Le code source de la chaîne est disponible sur https:// framagit.org/ cesardb/ SMASH. 2.2 Évaluation d’une carte de hauteur de neige Pléiades avec une carte de lidar aéroporté 

Site d’étude et données 

 Site d’étude

 Le site d’étude est une portion du bassin de la Tuolumne dans la Sierra Nevada, Californie (USA) (Figure 2.3). La rivière Tuolumne irrigue une partie de la plaine agricole californienne où sont produits les deux tiers des fruits secs et un tiers des légumes consommés aux EtatsUnis (Li et al., 2017). Elle alimente également le barrage de Hetch Hetchy, principale source d’alimentation en eau de l’agglomération de San-Fransisco et ses 7 millions d’habitants. La contribution de la fonte de la neige aux écoulements de surface est estimée à 53% du total des écoulements dans l’ensemble de l’ouest américain et à 70% dans les zones montagneuses 37 Chapitre 2. Cartographie de la hauteur de neige avec des images Pléiades (Li et al., 2017). La Californie a connu une période critique de sécheresse de 2014 à fin 2016 caractérisée notamment par des faibles accumulations de neige. En revanche, l’hiver 2016-2017 a battu des records d’accumulation, ce qui lui a valu le surnom de « snow-pocalypse »(Painter et al., 2017). Figure 2.3 – (a) Position du bassin de la Tuolumne aux Etats-Unis. (b) Des mesures de HTN sont acquises chaque hiver pendant la période de fonte sur le bassin de la Tuolumne dans le cadre du programme de l’ASO (ligne verte). Des acquisitions Pléiades (rectangle rouge) couvrent 200 km2 de cette région.

 Données de lidar aéroporté 

Le bassin de la Tuolumne est suivi dans le cadre de la campagne Airborne Snow Observatory (ASO) menée par la NASA (Figure 2.3). Une carte de hauteur de neige est mesurée par lidar aéroporté environ deux fois par mois pendant la période de fonte (Figure 2.4). Le lidar produit en premier lieu un nuage de points issus de la réflexion du rayon laser sur la surface (sol libre, neige) ou la végétation. La densité du nuage de points varie en fonction de l’angle d’intersection avec le sol et l’altitude entre 2 pnt.m−2 et plus de 10 pnt.m−2 (Jeff Deems, communication personnelle). Le nuage de points est filtré pour éliminer les réflexions sur la végétation et ainsi obtenir un MNT. Le MNT avec neige est comparé à un MNT sans 38 neige acquis le 13 octobre 2015. La carte de HTN a une résolution horizontale de 3 m pour le produit distribué gratuitement en ligne. Les 1100 km2 du bassin de la Tuolumne sont couverts à chaque vol. Le vol de l’ASO avec neige a eu lieu le 2 mai 2017. La précision verticale des HTN de l’ASO a été estimée à partir de 80 mesures par sondage manuel lors d’une campagne effectuée avant 2016 (Painter et al., 2016). Aucun biais n’a été alors détecté et la RMSE est 0,08 m (Painter et al., 2016). Une carte de hauteur de la végétation mesurée par lidar aéroporté le 13 octobre 2015 est utilisée en discussion. Le MNS hiver n’a pas pu être obtenu auprès de l’équipe qui produit les données de l’ASO malgré plusieurs demandes. En le reconstruisant à partir du MNT été et de la carte de HTN, on suppose qu’aucun post-traitement n’est appliqué à la carte de HTN de l’ASO après différenciation des MNT.