# OBSERVER LA TERRE EN 3D AVEC PLEIADES

Jean-Marc Delvit, Céline L'Helguen
CNES, Toulouse, France, Jean-Marc.Delvit@cnes.fr

## Résumé

La connaissance de modèles numériques d'élévation est un élément essentiel pour un grand nombre d'applications de télédétection, notamment à très haute résolution. Cette information d'élévation peut être déduite de couples d'images stéréoscopiques en combinant des techniques photogrammétriques et des techniques de type corrélation. L'amélioration de la résolution des système spatiaux d'observation de la Terre, comme Pléiades (pas d'échantillonnage au sol de 70cm) ainsi que les capacités croissantes d'acquisitions en mode stéréoscopique et multi-stéréoscopique permettent de générer de manière plus systématique des modèles numériques de terrain prenant également en compte le sursol fins (bâtis, véhicules, végétation), avec une précision de restitution de l'ordre du mètre.

Une méthode générique de restitution de modèles en « vrai » 3D fortement parallélisable est ici utilisée autorisant le traitement simultané de 2 à N images sur une même zone. Cette méthode est applicable directement à tout type de capteurs aériens ou spatiaux, comme Pléiades, Worldview2&3, SPOT6&7 ou Skybox.

**Mots-clés** : photogrammétrie, 3D, imagerie spatiale, Pléiades, corrélation, filtrage, nuage de points 3D

## Abstract

The knowledge of ground elevation is essential in most remote sensing applications especially for very high resolution images. This ground elevation information can be retrieved from a pair of stereoscopic images, by correlation methods. The improving resolution of Earth observation systems like Pléiades (70 cm GSD) and their increasing stereoscopic capabilities open up new horizons for automatic Digital Elevation Model generation and allow to consider buildings reconstruction with an accuracy better than one meter RMS. A method to compute 3D cloud points from Pléiades multi stereoscopic views is presented in this paper.

This method is generic and can be implemented on several processors. This method could be applied to Pléiades satellite, Worldview2&3, SPOT6&7 and Skybox.

**Keywords:** photogrammetry, 3D, space imagery, Pléiades, correlation, filtering, 3D cloud points

## 1. Introduction

Suite aux lancements des deux satellites agiles Pléiades 1A et 1B, respectivement en décembre 2011 et décembre 2012, de nombreuses acquisitions ont pu être réalisées en stéréoscopie et même en mode vidéo (Gleyzes et al., 2012, 2013). Une des exploitations possibles de ce type d'acquisition est la reconstruction du Modèle Numérique de Surface ou d'Elévation (MNS ou MNE) à une résolution fine.

Les méthodes mises en œuvre jusqu'à présent au CNES, nécessitent la mise en géométrie épipolaire des acquisitions ; les couples sont traités deux à deux et le résultat obtenu est une carte d'élévation, c'est-à-dire un résultat en 2,5D.

La nouvelle méthode présentée ici abandonne la géométrie épipolaire et génère un MNE directement sous la forme d'un nuage de points 3D. Les principales évolutions sont les suivantes :

* Travail en géométrie capteur ;

* Intégration d'une étape d'affinage dans l'algorithme ;

* Possibilité de traiter plus de deux images en même temps pour la génération du MNE ;

* Résultat sous forme de nuage de points 3D avec possibilité d'avoir des mesures sur les façades des bâtiments.

Le papier ci-après détaille l'algorithme mis en œuvre avant de présenter quelques exemples d'utilisation.

## 2. Méthode utilisée

### 2.1 Détails sur la corrélation

La méthode de corrélation utilisée repose sur la mesure de la ressemblance radiométrique entre les voisinages des pixels des deux images fournies en entrée de traitement. Pour calculer cette ressemblance, on définit une fenêtre de corrélation (ou vignette) qui correspond au voisinage du pixel sur lequel la fenêtre est centrée. La ressemblance est alors calculée entre la vignette de l'image de référence et la vignette de l'image secondaire. La mesure est généralement de la corrélation linéaire mais d'autres mesures peuvent être utilisées (corrélation linéaire : spatiale ou fréquentielle ; mesure statistique : algorithme du CRA (Cluster Reward Algorithm), distance à l'indépendance, divergence de Kullback, K-divergence de Lin, information mutuelle) (Giros et al., 2004 ; Delon et al., 2007).

Alors que la vignette de l'image de référence reste centrée sur le pixel $(i,j)_{REF}$ dont on cherche le correspondant dans l'image secondaire, la vignette de l'image secondaire est translatée au sein d'une seconde fenêtre, nommée fenêtre de recherche, centrée sur le pixel $(i + \text{décalage en ligne initial estimé}, j)$
+ décalage en colonne initial estimé) SEC (le décalage initial estimé entre les deux images pouvant être nul). Les coefficients de corrélation sont calculés dans l'ensemble de cette fenêtre d'exploration et la position à laquelle le coefficient est le plus élevé est considérée comme le correspondant du pixel (i,j) REF de l'image de référence. Cette position n'est pas nécessairement entière.

La corrélation est utilisée pour calculer un ensemble de points homologues pour l'affinage, ou pour calculer le modèle 3D. Toutes les images traitées deviennent tour à tour référence.

Les calculs de corrélation sont réalisés sur une grille de corrélation dont le pas est paramétré. Si les résultats du calcul de corrélation sont utilisés pour déterminer des points homologues en entrée d'affinage géométrique, le pas peut être choisi lâche à condition de fournir un nombre suffisamment grand de points appariés par rapport au nombre d'inconnues à déterminer lors de l'affinage des modèles. Si les résultats du calcul de corrélation sont utilisés pour déterminer les directions de visée, un pas dense (proche de 1) permettra d'obtenir un MNT ponctuel le plus détaillé possible.

## 2.2 Correction des modèles géométriques

L'objectif de l'affinage géométrique est de recaler par spatiotriangulation les modèles géométriques des différentes prises de vue traitées à partir des points mis en correspondance par corrélation.

En entrée on dispose des mesures (images et terrain) appariées. Ces mesures sont entachées d'une certaine erreur modélisée par une incertitude $σ_I$ pour les coordonnées image, $σ_{XY}$ pour les coordonnées planimétriques (latitude et longitude) et $σ_Z$ pour l'altitude.

<figure><img src="image_1.png" /><figcaption>Figure 1 : Principe général de la spatiotriangulation.</figcaption></figure>

En plus de ces mesures, les modèles géométriques possèdent des inconnues dont la valeur est contrainte dans un certain domaine. L'objectif est alors de calculer les valeurs des inconnues des modèles géométriques qui permettent de minimiser les écarts entre les mesures fournies en entrée de traitement (sorties de la corrélation) et les estimations de la position des points obtenues par les modèles géométriques corrigés, tout en respectant les incertitudes indiquées sur les mesures (cf. fig 1.). Cet affinage des modèles géométriques prend en compte tous les étalonnages fins réalisés lors de la recette en vol Pléiades : biais d'alignement (DeLussy et al., 2012), cartographie du plan focal (Greslou et al., 2012; Kubik et al., 2012; Lebegue et al., 2012). Cet affinage est réalisé par moindres carrés, d'autres techniques de minimisation sous contraintes peuvent également être utilisées.

## 2.3 Algorithme de calcul d'intersection de visées

L'algorithme général de calcul d'intersection de visées est décrit ci-dessous (Delvit et al., 2007).

Soit $\Delta$ une droite passant par un sommet S et de vecteur directeur $\vec{V}$.

Soit M un point quelconque et $\delta$ la distance de M à la droite $\Delta$. On a :

$$ 
\begin{aligned}
\delta^2 &= \| \overrightarrow{SM} \wedge \overrightarrow{V} \|^2 = \overrightarrow{SM}^2 - (\overrightarrow{SM} \cdot \overrightarrow{V})^2 \\
d(\delta^2) &= 2 \cdot (\overrightarrow{SM} \cdot d\overrightarrow{M} - \overrightarrow{SM} \cdot \overrightarrow{V} \cdot d\overrightarrow{M} \cdot \overrightarrow{V}) = 2 \overrightarrow{dM} \cdot (\overrightarrow{SM} - \overrightarrow{SM} \cdot \overrightarrow{V} \cdot \overrightarrow{V}) \\
d(\delta^2) &= 2 \overrightarrow{dM} \cdot (\overrightarrow{SM} - VV' \overrightarrow{SM}) = 2 \overrightarrow{dM} \cdot [(Id - VV') \overrightarrow{SM}] \\
d(\delta^2) &= 2 \overrightarrow{dM} \cdot [(Id - VV') [M - S]]
\end{aligned}
$$

Pour un ensemble de n droites de visées $\Delta_i$ avec $(S_i, \vec{V}_i)$, on recherche le point M qui minimise la somme des carrés des distances :

$$ \sum_i \delta_i^2 = \sum_i \| S_i \cdot \overrightarrow{M} \wedge \overrightarrow{V}_i \|^2 = \sum_i \left[ \overrightarrow{S_i M}^2 - (\overrightarrow{S_i M} \cdot \overrightarrow{V}_i)^2 \right] $$

En différenciant, on obtient :

$$ 
\begin{aligned}
d\left(\sum_{i} \delta_{i}^{2}\right) &= 2 \overrightarrow{d M} \cdot\left[\sum_{i}\left(\mathrm{Id}-V_{i} V_{i}^{\prime}\right)\left[M-S_{i}\right]\right] \\
d\left(\sum_{i} \delta_{i}^{2}\right) &= 2 \overrightarrow{d M} \cdot\left[\sum_{i}\left(\mathrm{Id}-V_{i} V_{i}^{\prime}\right) M-\sum_{i}\left(\mathrm{Id}-V_{i} V_{i}^{\prime}\right) S_{i}\right]
\end{aligned}
$$

La solution est donc :

$$ M = \left[ \sum_{i} (Id - V_{i} V_{i}^{\prime}) \right]^{-1} \sum_{i} [(Id - V_{i} V_{i}^{\prime}) S_{i}] $$

L'intérêt de cette méthode est de plusieurs ordres : pour chaque nouvelle droite, on constitue chaque terme $Id - V_i V_i'$ puis $(Id - V_i V_i') S_i$ qu'on cumule à la volée ;

une seule inversion d'une matrice 3x3 ;
pour éliminer une visée jugée trop lointaine, on soustrait simplement sa contribution dans $\sum_i (Id - V_i V_i')$ et

$$ \sum_{i} [(Id - V_{i} V_{i}^{\prime}) S_{i}] $$

Dans le cas particulier de la génération de MNT ponctuel par intersection de visées, on calcule dans un premier temps les intersections de visées entre deux images puis dans un second temps on calcule ces intersections de visées de façon inter-couple.

## 2.4 Filtrage de points aberrants

Un certain nombre de filtres peuvent être utilisés sur le nuage de points obtenu par intersection de visées. Il s'agit de supprimer les points statistiquement aberrants. Deux catégories se distinguent :

* des filtres globaux qui considèrent l'ensemble du nuage de points :
o filtrage global en Z : une altitude minimale en dessous de laquelle les points sont supprimés est définie;

o filtrage global en Z à K-sigma : un coefficient multiplicatif K est défini. Les points dont l'altitude est située en dehors de l'intervalle « moyenne ± K*écart-type » sont supprimés. La moyenne et l'écart-type sont calculés sur l'altitude de l'ensemble des points du nuage;

* des filtres locaux qui ne considèrent qu'une partie de l'espace occupé par le nuage de points :

o filtrage local en Z à K-sigma : un coefficient multiplicatif K et un pas sont définis. Le sol est divisé en N vignettes carrées de côté « pas » (sans recouvrement). Pour chaque vignette, la moyenne et l'écart-type sont calculés sur l'ensemble des points dont les coordonnées au sol sont comprises dans cette vignette et les points dont l'altitude est située en dehors de l'intervalle « moyenne ± K*écart-type » sont supprimés ;

<figure><img src="image_2.png" /></figure>

o filtrage local par sphère : un rayon de sphère et un nombre de points minimal qui doit être situés dans cette sphère sont définis. La sphère est centrée sur le point dont on veut déterminer la validité (supprimé ou conservé).

<figure><img src="image_3.png" /></figure>

o D'autres types de filtres peuvent être appliqués sur les nuages de points en sortie de l'algorithme. Proche du filtrage local par sphère on peut citer le filtrage local par cylindre ou par parallélépipède rectangle. La librairie PCL (Point Cloud Library) (Rusu et al., 2011) propose ainsi un certain nombre de filtres qui peuvent être pratiques à utiliser sur le nuage brut ou sur lequel persistent certaines erreurs.

### 3. Exemples d'utilisation

Plusieurs tests ont été menés sur des produits Pléiades 1A et Pléiades 1B. Les images sont soit acquises en stéréoscopie (deux images acquises successivement sur un même site : Gizeh, Everest, Monument Valley),

soit sous la forme d'une tristéréoscopie (trois images acquises successivement sur un même site : Port-au-Prince), soit sous la forme d'une vidéo (jusqu'à une vingtaine d'images acquises successivement sur un même site).

L'objectif de ces tests est de montrer dans un premier temps les résultats qu'il est possible d'obtenir sur de larges zones sans urbanisation puis sur des zones plus urbaines. L'intérêt de la reconstruction en 3 dimensions est ainsi mis en évidence avec la présence de points sur les façades des bâtiments. Pour le cas test sur Port-au-Prince (Pléiades) des références altimétriques sont disponibles. L'écart à ces références est donc calculé. Pour le test de l'Everest (Fig. 2) les deux images utilisées en entrée de traitement sont les suivantes. Le couple stéréoscopique Pléiades à faible b/h (angle stéréoscopique) inférieure à 0,2 radian est composé de deux images visuellement très semblables, ce qui est favorable aux algorithmes de corrélation. Le faible b/h a le grand avantage de minimiser les occlusions. Afin de profiter de tel angle stéréoscopique (b/h), il est nécessaire de disposer de modèle géométrique précis.

<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /></figure>

**Figure 2** : Images fournies en entrée de traitement pour le test sur l'Everest – le couple stéréoscopique Pléiades à faible b/h est composé de deux images visuellement très semblables, ce qui est favorable aux algorithmes de corrélation.
<figure><img src="image_6.png" /><figcaption>**Figure 3 :** Résultats de corrélation et du filtrage aller-retour du test sur l'Everest.</figcaption></figure>

Les étapes de corrélation et de filtrage associé (Fig. 3) sont fondamentales pour obtenir un bon résultat final (Fig. 4). Les résultats ci-dessous montrent que le filtrage permet de rejeter une large majorité des points faux. Le nuage de points obtenu par intersection des directions de visées est sans artefact, le filtrage local permettant de supprimer les quelques points faux persistants.

<figure><img src="image_7.png" /><figcaption>**Figure 4 :** Résultats sur l'Everest.</figcaption></figure>

Le test de reconstruction 3D sur les pyramides de Gizeh permet de bien voir la présence de points sur les façades visibles des pyramides (Fig. 5). Les terrains alentours sont également bien reconstruits avec la présence de dénivelés fins.

<figure><img src="image_8.png" /><figcaption>**Figure 5 :** Résultats sur les pyramides de Gizeh.</figcaption></figure>

Les tests sur Monument Valley montrent que des points sont présents même sur des faces abruptes (Fig. 6). Par contre les zones non visibles et dans l'ombre sont vides de points. Même des cheminées rocheuses de plus petite surface sont bien reconstruites. Dans le deuxième test, la radiométrie de l'image a été associée à chaque point.

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 6 :** Résultats sur Monument Valley avec un couple stéréoscopique Pléiades très dépointé (40° en roulis).</figcaption></figure>

## 4. Résultats quantitatifs

Les tests menés sur Port-au-Prince utilisent une acquisition tri stéréo Pléiades 1A du 19 juillet 2012. Quatre références sont disponibles sous la forme de 4 zones acquises par un LIDAR. L'objectif est de comparer les résultats obtenus par intersection des directions de visée et ceux de la référence LIDAR sur chacune des 4 zones.

Les 4 zones sont reconstruites sans affinage préalable. Les 6 couples de prises de vue possibles sont utilisés pour la reconstitution 3D (1/2, 1/3, 2/3, 2/1, 3/1 et 3/2). Deux filtrages sont appliqués aux nuages de points obtenus : un filtrage à K-sigma (sigma=2 et pas=15m) suivi d'un filtrage sphère (rayon sphère=10m et nombre de points dans la sphère=5).
La comparaison à la référence est menée en utilisant un outil de calcul de distance entre deux nuages de points. Les nuages de points sont projetés selon Z pour obtenir une carte en 2.5D du résultat. Seules les 2 zones les moins bonnes sont présentées dans le détail ci-après.

La zone 1 est une zone fortement urbanisée (Fig 7) avec des bâtiments de petite taille. Cette zone est donc difficile à reconstruire avec des produits Pléiades. De plus un dénivelé assez important est visible dans la référence (Fig. 8), sur la droite de l'image. Malgré tout la reconstruction 3D obtenue est de bonne qualité. Son calage planimétrique est considéré acceptable pour cette étude. Les statistiques sont très bonnes puisqu'en moyenne l'écart (biais) en altitude est de seulement 0,26m avec un écart-type de 1,64 (Fig. 9).

<figure><img src="image_12.png" /><figcaption>**Figure 7 :** Image Pléiades 1A nadir du 19/07/2012 sur la zone 1 de Port-au-Prince, zone de type urbain dense.</figcaption></figure>

Contrairement à la zone 1, la zone 2 est une zone qui comporte très peu d'habitations (Fig. 10). Par contre le terrain est plus élevé et plus accidenté également. La reconstruction 3D est de bonne qualité. Les principales erreurs sont dues à une trop faible estimation des extrema. La référence LIDAR est décrite sur la figure 11. Ainsi les points les plus bas sont surélevés tandis

<figure><img src="image_13.png" /><figcaption>**Figure 8 :** Référence LIDAR sur la zone 1 de Port-au-Prince (en rouge hautes altitudes, en bleu altitudes basses).</figcaption></figure>

que les points les plus hauts sont abaissés. Les statistiques (Fig. 12) sont donc un peu moins bonnes

que pour la zone 1. La moyenne de la différence entre le nuage de points et la référence est de 0,40m tandis que l'écart-type est de 1,76.

<figure><img src="image_14.png" /><figcaption>**Figure 9 :** Différence selon Z (altitude) entre le nuage de points obtenu par reconstruction 3D et la référence LIDAR de la zone 1 (de bleu foncé -5m à rouge +5m).</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 10 :** Image Pléiades 1A nadir du 19/07/2012 sur la zone 2 de Port-au-Prince, zone de type naturel.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Figure 11 :** Référence LIDAR sur la zone 2 de Port-au-Prince (en rouge hautes altitudes, en bleu altitudes basses).</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>**Figure 12** : Différence selon Z (altitude) entre le nuage de points obtenu par reconstruction 3D et la référence LIDAR de la zone 2 (de bleu foncé -5m à rouge +5m).</figcaption></figure>

Les cartes de différence en altitude, entre les cartes obtenues par reconstruction 3D et les références LIDAR, ont des décalages moyens tous inférieurs au mètre et des écarts-types faibles (inférieurs à 1,8m) (Tab. 1). Ceci confirme donc la qualité des résultats obtenus à partir des images Pléiades acquises sur Port-au-Prince. Malgré tout un affinage de bonne qualité en prenant des points d'appui au sol sur les références LIDAR pourrait permettre d'améliorer un peu plus encore les statistiques des cartes de différence. Ainsi on pourrait s'affranchir des erreurs de calage en longitude et latitude et corriger une possible inclinaison par rapport au plan des références. Par contre il est essentiel que l'affinage, s'il a lieu, prenne en compte l'ensemble des zones étudiées pour conserver une cohérence entre les résultats.

<table><caption>**Tableau 1** : Statistiques des cartes de différence en altitude entre les cartes obtenues par reconstruction 3D et les références LIDAR (en mètres).</caption><thead><tr><th rowspan="2">en<br>mètre</th><th colspan="2">en X</th><th colspan="2">en Y</th><th colspan="2">en Z</th></tr><tr><th>Moy</th><th>Moy</th><th>Moy</th><th>Sigm<br>a</th><th>Moy</th><th>Sigm<br>a</th></tr></thead><tbody><tr><td>Zone 1</td><td>-0,09</td><td>0,12</td><td>0,26</td><td>1,64</td><td>-0,09</td><td>0,12</td></tr><tr><td>Zone 2</td><td>-0,28</td><td>0,09</td><td>0,39</td><td>1,5</td><td>-0,28</td><td>0,09</td></tr><tr><td>Zone 3</td><td>-0,04</td><td>0,02</td><td>0,68</td><td>1,23</td><td>-0,04</td><td>0,02</td></tr><tr><td>Zone 4</td><td>0,01</td><td>0,11</td><td>0,84</td><td>1,21</td><td>0,01</td><td>0,11</td></tr></tbody></table>

## 5. Conclusion

Au vu de la forte augmentation des acquisitions multi-stéréoscopiques spatiales, il est nécessaire de disposer d'algorithmes et outils autorisant une restitution du relief en « vraie » 3D. La méthodologie présentée reste très générique et peut s'appliquer à tout type de prises de vue (terrestre, aérien ou spatial, matriciel ou pushbroom) – en particulier pour les systèmes spatiaux tels Pléiades, Geoeye, Worldview, SPOT6&7, ou même sur les vidéos Skybox. Le MNS (Modèle Numérique de Surface 3D) pourrait devenir un produit standard au même titre qu'une simple ortho-image, ou une orthomosaïque.

Les analyses quantitatives montrent que la restitution altimétrique pour un système bien dimensionné a une précision de 2 pas au sol RMS environ. Si la modélisation géométrique du système stéréoscopique est correctement maitrisée, des B/H de l'ordre de 0,12 permettent, pour des acquisitions stéréoscopiques nadir, de minimiser le nombre d'occultations tout en ayant une bonne performance altimétrique.

## Références

De Lussy F., Greslou D., Dechoz C., Amberg V., Delvit JM., Lebègue L., Blanchet G., Fourest S., 2012. *Pléiades HR in flight geometrical calibration: Localisation and mapping of the focal plane*, The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, XXXIX-B1, 519-523.

Delon J., Rougé B., 2007: “Small Baseline Stereovision”, Journal of Mathematical Imaging and Vision, Volume 28 Issue 3.

Delvit JM., Favé P., Gachet R., 2007: “The geometric supersite of Salon de Provence”, La Revue Française de Photogrammétrie et de télédétection.

Giros A., Rougé B., Vadon H., 2004. “Appariement fin d'images stéréoscopiques et instrument dédié avec un faible coefficient stéréoscopique”, French Patent No. 0403143.

Gleyzes A., Perret L., Kubik P., 2012. *Pléiades system architecture and main performances*, The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences., Vol. XXXIX-B1, 537-542

Gleyzes A., Perret L., Cazala-Hourcade E., 2013. *Pléiades system is fully operational in orbit*, Proceedings EARSeL, 445-460, Matera, Italy.

Greslou D., de Lussy F., Delvit JM., Dechoz C., Amberg V., 2012. *Pléiades-HR innovative techniques for Geometric Image Quality Commissioning*, The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences., Vol. XXXIX-B1, 543-547.

Kubik P., Lebègue L., Fourest S., Delvit JM., de Lussy F., Greslou D., Blanchet G., 2012. *First in-flight results of Pléiades 1A innovative methods for optical calibration*, ICSO Ajaccio, France.

Lebègue L., Greslou D., de Lussy F., Fourest S., Latry C., Kubik P., Delvit JM., 2010. *Pléiades-HR image quality commissioning foreseen methods*. Proceedings IGARSS: 1675-1678, Honolulu, Hawaii

Rusu R. B., Cousins S., 2011. “3D is here: Point Cloud Library (PCL)”, IEEE International Conference on Robotics and Automation (ICRA), Shangai, China.