Equipe d'imagerie Radar

Ces travaux ont été réalisés dans la thèse de Charles Deledalle dirigée par Florence Tupin et co-encadrée par Loïc Denis. Le but était d'adapter l'algorithme des moyennes non locales (A. Buades, 2005) aux images Radar multi-modalités (images d'amplitude, interférogrammes SAR, images polarimétriques). Au final, un filtre performant a été conçu pouvant se généraliser aux images présentant des bruit non nécessairement Gaussien, aux images multidimensionnelles et en particuliers aux différents types d'images Radar.

L'algorithme des moyennes non locales débruite une image en calculant pour chaque pixel s la moyenne des valeurs des pixels similaires u_t. Cette similarité est exprimée en fonction de la distance Euclidienne entre les voisinages des pixels s et t. Cette distance Euclidienne est bien adaptée aux images présentant un bruit additif Gaussien.

Les moyennes non locales procèdent à un moyennage en tout pixels s des valeurs des pixels t pondérées en fonction de la similarité entre deux fenêtres centrées respectivement sur s et t

Pour étendre le filtre aux bruits non Gaussien, nous proposons d'estimer en chaque pixel s le meilleur paramètre sous-jacent au sens du maximum de vraisemblance pondérée par rapport aux pixels u_t :

où le poids w(s,t) définit la similarité entre les pixels s et t. Sur le principe de l'algorithme des moyennes non locales, on définit ce poids par la similarité entre deux fenêtres Δ_s et Δ_t centrées respectivement sur s et t. On suppose alors que cette similarité est directement liée à la probabilité de similarité par :

avec λ un réel positif arbitraire et h un paramètre contrôlant la puissance du filtre. La probabilité de similarité peut se décomposer en deux termes, la vraisemblance de similarité d'une part p(u_Δs, u _Δt | θ_Δs = θ_Δt), et la similarité a priori d'autre part p(θ_Δs = θ_Δt) que l'on propose de modéliser à l'aide d'une estimation précédente.

(a)

(b)

(c)

(d)

(e)

(f)

(a-c) Images bruitées respectivement avec un bruit Gaussien d'écart-type σ = 40, un bruit impulsionnel à 70% et un bruit multiplicatif de speckle. L'image (c) est une image radar de Toulouse ©DGA ©ONERA. (d-f) Images débruitées obtenues avec le filtre proposé dérivé pour chaque type de bruit.

Le filtre a été dérivé pour différents types de bruit, dont le bruit Gaussien avec écart type σ = 40, le bruit impulsionnel à 70% et le bruit multiplicatif de speckle présent en imagerie radar à ouverture synthétique. Les résultas sont donnés ci-dessus. Dans le cas du bruit Gaussien, on retrouve bien l'algorithme des moyennes non-locales auquel on a rajouté un schéma itératif raffinant l'estimation. Cette méthode nous permet aussi d'exhiber une mesure de similarité bien adaptée au cas de l'imagerie radar (mono-vue ou multi-vue) et donnée par

avec A_s et A_t deux valeurs d'amplitudes et R_s et R_t les réflectivités sous-jacentes.

Ces travaux sont réalisés dans la thèse d' Aymen Shabou dirigée par Florence Tupin. L'idée consiste à exploiter les approches Markoviennes et des algorithmes d'optimisation par coupe minimale sur graphes pour résoudre efficacement et robustement le problème de déroulement de phase en interférométrie radar multi-canal. La fonction de vraisemblance est définie à partir du modèle de bruit interférométrique et les données multi-canaux. Le model a priori est la Variation Totale.

avec E la fonction d'énergie, Φ l'image des phases mesurées (phases enroulées), h l'image des phases à estimer (phases déroulées), p et q les indices respectifs des pixels et des canaux des données (M interférogrammes utilisés), f la fonction de vraisemblance mono-canal, γ le coefficient de coherence, β le coefficient de régularisation, w la pondération dépendant de la connexité utilisée et α caculée à partir des paramètres d'acquisition du canal considéré (base et fréquence).

La première partie de ces travaux est réalisée en collaboration avec Giampaolo Ferraioli (Dipartimento per le Tecnologie- Università di Napoli - Parthenope) et est consacrée à l'étude de l'approche de déroulement de phase multi-canal en se fondant sur ce nouveau modèle d'énergie et les algorithmes de minimisation par coupe minimale proposés dans la littérature, principalement les algorithmes de minimisation exacte d'Ishikawa et non-exacte de Boykov et al. (l'α-extension). Les résultats obtenus, sur des données simulées et réelles, comparés à ceux obtenus par d'autres approches de déroulement de phase multi-canal, montrent bien l'efficacité de notre approche en terme de temps de calcul et sa robustesse par rapport aux fortes discontinuités et bruits présents dans les données.

(a)

(b)

(c)

(d)

(a) Carte de profils originale (vérité terrain), (b) Image de phases enroulées (un des interférogrammes bruités), Profils obtenus par (c) l'approche proposée utilisant l'algorithme d'Ishikawa, (d) une approche classique (algorithme des différences de phases).

La deuxième partie est réalisée en collaboration avec Jérôme Darbon (UCLA, Californie, USA) et est consacrée à la propsition de nouveaux algorithmes de minimisation non-exacte par coupe minimale présentant un compromis entre qualité du minimum atteint et complexité en mémoire nécessaire pour la construction du graphe. En pratique, ces algorithmes sont capables de converger vers le minimum global de l'énergie en utilisant nettement moins de mémoire que les approches de minimisation exacte. Ils sont principalement utiles si nous devons traiter des données de mesure de phases interférométriques de grandes tailles, très bruitées et à fortes discontinuités. En effet, dans ce cas-là, les approches de minimisation non-exacte par mouvement de partition binaire, telle que l'α-extension, convergent vers des minima de mauvaise qualité et les approches de minimisation exacte, telle que celle d'Ishikawa, nécessitent la construction de graphes très coûteux en mémoire. Les algorithmes proposés se fondent sur des mouvements de partition dits multilabels qui permettent d'explorer plus d'optima locaux que les algorithmes à base de mouvements binaires. Ils convergent ainsi vers des solutions de meilleure qualité, voir même la solution exacte au problème de minimisation en construisant des graphes de plus petites tailles.

(a)

(b)

(c)

(d)

(e)

(f)

(a) Carte de profils originale (vérité terrain), (b) Image de phases enroulées (un des interférogrammes bruités), Profils obtenus par l'approche proposée utilisant l'algorithme de minimisation (c) α-extension, (d) Extension avec 8 labels, (e) Extension avec 16 labels, (f) Ishikawa (ou Extension avec 128 labels).

Ces travaux sont réalisés dans le cadre de la thèse d'Hélène Sportouche, dirigée par Florence Tupin (Telecom ParisTech) et co-encadrée par Léonard Denise (Thalès Communications). Il s'agit d'une thèse sous contrat Cifre menée dans le cadre d'un partenariat entre le département TSI de Telecom ParisTech (Paris) et le départemnt IMINT de l'entreprise Thalès Communications (Massy).

L'objectif de ces travaux est d'exploiter conjointement des données optiques et radar à résolution métrique en vue d'une extraction et reconstruction 3D des bâtiments présents dans une scène urbaine. Nous proposons un enchainement de méthodes permettant d'obtenir, de façon semi automatique, une détection et une reconstruction du bâti. L'idée consiste à combiner les informations complémentaires issues des données optiques et radar, en utilisant des primitives dites caractéristiques du bâtiment et propres à chaque type de capteur. Nous modélisons simplement notre bâtiment par un parallépipède 3D (façades verticales et toit plat) décrit par deux sortes de paramètres: son contour planimétrique rectangulaire et sa hauteur. L'approche proposée se décompose en trois principales étapes: une première étape d'extraction des empreintes potentielles des batiments en monoscopie optique, une seconde étape de projection et de recalage des arêtes au sol détectées avec leurs primitives radar homologues (échos brillants issus des coins réflecteurs sol /mur) et une dernière étape visant simultanément à valider la présence des bâtis et à estimer leur hauteur, en vue d'une reconstruction 3D.

Etape 1 : Extraction des empreintes potentielles des batiments en monoscopie optique. Il s'agit ici d'extraire les empreintes des bâtis à partir d'une seule image optique. La plupart de nos images étant acquises avec une visée quasiment au nadir, la localisation planimétrique des empreintes extraites peut aussi bien être associée à celle des arêtes au sol qu'à celles des bords de toits. Un procédé en deux phases est mis en oeuvre afin d'obtenir d'abord une première carte identifiant globalement les zones de bâtis potentielles, et ensuite une seconde carte dite affinée, indiquant la localisation précise des contours de bâtis détectés. La première carte résulte d'un procédé hiérarchique à 2 passes, basé sur un critère d'adéquation géométrique entre une forme de référence (ici rectangulaire) et une forme test. Ce critère est testé sur l'ensemble des régions appartenant aux images constituant le Profil Morphologique Différentiel (PMD) de l'extrait optique considéré. L'utilisation d'un tel profil nous permet de tirer profit de ses propriétés de simplification des images à différentes échelles et de préservation des formes. L'affinage des contours fait appel à l'utilisation d'une transformée de Hough locale, appliquée sur l'image de gradiant binaire associée à l'extrait considéré et tenant compte de contraintes d'orthogonalité. Celle-ci est mise en oeuvre sur chacune des zones de recherche déduites de la carte d'identification globale et permet de détecter un ensemble d'angles droits susceptibles de correspondre à des coins de bâtiments. Un procédé de génération de rectangles, à partir des angles droits détectés, et basé sur l'optimisation d'une énergie sur l'image de gradient, est alors appliqué pour obtenir la carte affinée des contours dans l'optique. On présente ci-dessous les résultats obtenus après la phase 1 et 2 sur un extrait optique panchromatique, situé dans la région de Marseille (France). Après la phase 1, une identification globalement satisfaisante des zones potentielles de bâtis est obtenue. Bien sûr, la principale limitation est la forme a priori rectangulaire. Quelques fausses alarmes, essentiellement dues à la présence de parcelles rectangulaires homogènes, sont présentes. Après la phase 2, une délimitation précise des contours de bâtiments est obtenue.

(a)

(b)

(c)

(d)

Résultats de l'étape 1 (extraction des emprises des bâtiments). Extrait optique panchromatique, capteur : Quickbird, résolution : 0.6 m, région : Marseille (France). ©DigitalGlobe (a) Extrait initial (b) Résultat à l'issue de la phase 1 (identification de zones potentielles) : Les régions issues de la première passe du procédé sont respectivement représentées en rouge et en bleu pour celles relevant de la partie "ouvertures" et "fermetures" du PMD. Celles issues de la seconde passe sont représentées en vert et en jaune. (c) Résultat à l'issue de la phase 2 (affinage des contours) : Les contours obtenus sont représentés en rouge. (d) Résultat à l'issue de la phase 2 : image binaire associée à (c) avec labels des bâtis obtenus.

Etape 2 : Projection et Recalage des contours optiques détectés avec leurs primitives radar homologues. La projection des rectangles optiques dans l'image radar requiert la connaissance des paramètres capteur (systèmes d'acquisition) ainsi qu'une information de hauteur approximative. On dispose pour l'image optique Quickbird d'un modèle dit RPC (Rational Polynomial Coefficients) et pour l'image radar TerraSAR-X du véritable modèle physique. L'information de hauteur est fournie par un MNE DIMAP produit par la société SPOT. L'étape de projection est réalisée en utilisant les équations fondamentales optiques et radar, via un référentiel cartographique intermédiaire. Afin d'améliorer l'estimation du modèle joint de projection (optique / radar), l'opérateur peut saisir quelques points d'appui. On présente ci dessous les résultats de la projection obtenus en utilisant le modèle joint (avec ou sans points d'appui). L'étape de recalage vise à supprimer les légères erreurs de superposition encore présentes entre primitives homologues optiques et radar. Celles-ci sont essentiellement dûes à deux raisons : Premièrement nous avons une connaissance insuffisante du MNT. Deuxièmement, le modèle RPC optique utilisé est susceptible d'introduire de petites erreurs de localisation dans le repère cartographique intermédiaire. Ainsi des erreurs aussi bien en altimétrie (coordonnées z) et qu'en planimétrie (coordonnées x et y) peuvent se produire et sont à l'origine de décalages (translations) des structures projetées dans l'image radar. Ainsi, nous proposons une stratégie d'amélioration du recalage, visant simultanément à affiner le MNT et à obtenir un très bon appariement des primitives homologues caractéristiques au sol. Etant donné qu'après projection, le recalage obtenu est déja de bonne qualité, seules quelques petites translations locales sont envisagées. L'idée consiste à recaler, pour chaque rectangle optique, les deux segments susceptibles d'être imagés au sol dans l'image radar, avec leurs primitives radar homologues c'est à dire avec les échos brillants correspondant aux doubles réflexions des ondes radar au niveau des coins sol / murs. Un critère énergétique est défini sur l'image radar des radiométries et permet d'obtenir une superposition optimale entre les arêtes optiques au sol et les échos doubles radar. Enfin, une étape de régularisation, tenant compte notamment des variations lentes du MNT, est mise en oeuvre. Celle-ci permet de sélectionner et propager les tranlations les plus robustes. On présente ci dessous les résultats de projection et de recalage.

(a)

(b)

(c)

(d)

Résultats de l'étape 2 (projection et recalage). Extrait radar, capteur : TerraSAR-X, résolution : 1.1 m, région : Marseille (France). ©Infoterra Les segments et coins optiques projetés et recalés sont respectivement représentés en rouge et en bleu. (a) Résultat à l'issue de la projection avec modèle joint initial. (b) Résultat à l'issue de la projection avec modèle joint amélioré par points d'appui (5 points). (c) Résultat à l'issue de la projection avec modèle joint amélioré et avant le recalage (zoom sur un batiment). (d) Résultat à l'issue de la projection avec modèle joint amélioré et après le recalage (zoom sur un batiment).

Etape 3 : Validation de la présence des bâtis et estimation des hauteurs en vue d'une reconstruction 3D. Il s'agit alors de tirer profit des informations apportées par l'image radar en réduisant le taux de fausses alarmes parmi les bâtis potentiels (phase de validation) détectés dans l'optique et en déduisant conjointement une hauteur estimée pour chaque bâti validé (phase d'estimation de la hauteur du bâti). Ces phases reposent sur la présence de primitives ou régions caractéristiques du bâtiment dans l'image radar, parmi lesquelles on recense : les échos doubles sol / mur, les échos provenant des bords de toits, les zones de toits (toit avec recouvrement ou toit seul) et les zones d'ombres. On s'intéresse donc aussi bien aux primitives au sol qu'à celle en hauteur. L'idée consiste à calculer conjointement un critère statistique global et un critere radiométrique local, pour un ensemble de partitionnements tests de la scène en différentes régions caractéristiques, générées, pour chaque bâti, en fonction des hauteurs testées et des emprises au sol connues. Une optimisation sur ces critères, suivie d'un seuillage, nous permet alors de valider les bâtis (adéquation effective avec le modèle prédit) et de déduire une hauteur pertinente pour chacun d'eux. On présente ci dessous, les résultats issus de la validation, de l'estimation des hauteurs et de la reconstruction 3D, obtenus pour la scène urbaine considérée.

(a)

(b)

Résultats de l'étape 3 (reconstruction 3D obtenue suite à la validation des bâtis et à l'estimation des hauteurs). Résultats issus de l'optimisation des critères. (a) Carte des hauteurs estimées (en niveaux de gris) pour les bâtis validés. (b) Reconstruction 3D de la scène.