Keywords: Projection and prior

Je tiens tout d'abord à remercier les rapporteurs qui m'ont permis d'améliorer ce manuscrit et de préparer ma soutenance. Je remercie plus généralement l'ensemble des membres du jury qui ont assisté à ma soutenance et tout particulièrement ceux venant de loin.

Je remercie les membres du DEMR qui m'ont accueilli au début de ma thèse et n'ont jamais été avares de conseils ou propositions. Parmi eux Frédéric et Abigael, j'ai apprécié que vous m'intégriez dans vos (multiples) pauses café ainsi que vos coups de main quand j'ai pu en avoir besoin. Hélène, l'intérêt que tu portais à mes travaux ainsi que tes encouragements ont été très motivants tout au long de ma thèse.

La liste des gens à citer à Telecom est longueet je crains d'en oublier certains. Bien sûr, Florence, je te suis très reconnaissant de ton investissement. Ta rigueur et ton exigence m'ont beaucoup appris. Mais ce ne serait pas te faire justice que de citer ces seules qualités car j'ai aussi apprécié les valeurs que tu portes (du moins celle que j'ai pu te voir défendre notamment en conseil de département). Jean-Marie, sans toi je ne sais pas sur quoi j'aurais pu travailler ! Ta facilité d'accès et les échanges que nous avons pu avoir sur de nombreux sujets m'ont été plus qu'utiles. Enn, Pietro, Jean-Marc, Alasdair, Samy et Said, c'était un plaisir de pouvoir plaisanter au détour d'un couloir ou autour d'une bière. Chez les doctorants, je suis très reconnaissant à ceux qui ont su si bien m'accueillir. Flora, Sylvain (2D), Hélène et Guillaume notamment, sans vous ça n'aurait dénitivement pas été pareil ! Sonia, Paul et Antoine H., parce que vous n'étiez pas toujours à Telecom, nous n'avons pas pu autant nous voir mais je garde des bons souvenirs de nos échanges. Chenguang et Weiying, je garderai en mémoire notre partage de bureau et les moments que nous y avons partagé ! Alessio, on a fait notre thèse en parallèle et je n'aurais pas pu demander meilleur compagnon de route. Peut être qu'un jour je vais vraiment m'inscrire à la boxe avec toi qui sait. Je sais que le ambeau s'est transmis car Nicolas, Antoine P, Nicolas, Xiangli et Arthur, je suis très content que vous soyez dans les parages ! Je vous souhaite d'aller au bout de votre thèse et d'être ers de vous. Enn, Bastien si j'ai pu te prendre pour une nuisance intra-bureau, je n'aurais pas pu plus me tromper ! Ta bonne humeur et ton ouverture d'esprit font de toi le co-bureau parfait. Sylvain (3D), Adrien, Malik, Jeremy et les autres, vous étiez bizarres avec votre dimension en plus mais quand on arrivait à se croiser ça faisait plaisir. Laurent, ça a toujours été un plaisir de partager avec toi que ce soit sur des questions scientiques ou générales. J'ai aussi trouvé agréable que, bien que tu connaisses beaucoup de monde, tu restes très inclusif en conférence. C'est bien sûr très apprécié quand on est jeune chercheur. Yue, ta disponibilité et ta bonne humeur ont également été des gros plus dans cette collaboration. Loïc, je sais que je te dois beaucoup dans cette thèse. Tu réunis des qualités qui font de toi un excellent encadrant : clarté, ecacité, disponibilité et gentillesse. Enn, tu es toujours prêts à t'intéresser à une idée quand bien même elle pourrait être bancale ce III qui est très stimulant et encourageant. Rémy, en plus d'être un super post-doc, tu étais bienveillant et de précieux conseil. J'espère pouvoir travailler avec toi à l'avenir. Mélissa, merci pour ton soutien et ton amour sans lesquels, il aurait été bien moins simple d'aller au bout de ces trois ans. Mes parents, merci de croire en moi, de m'encourager et de vous intéresser à ce que je fais. La Mifa, vous êtes des champions, merci à vous aussi de votre soutien sans faille. Les anciens de prépas qui sont entrés dans Ph.D. game, c'est un peu grâce à vous que je me suis lancé la dedans. Donc, spéciale dédicace comme on dit ! Les Pandor-Margulis, merci de m'inclure autant et de votre amitié.

Les anciens collocs, j'ai fait tout ça pour vous impressionner, j'espère que ça a marché.

Enn, merci Florence Besnard, Salimatou Yansane et Marie-Laure Chauveaux pour votre aide et votre patience (notamment vis à vis de mon temps de réponse).

Introduction

Le suivi de la terre par images aériennes ou spatiales est depuis longtemps un sujet d'intérêt dans une grande variété de disciplines allant de l'urbanisme à la séismologie en passant par l'étude de la biodiversité. La possibilité d'observer certaines scènes directement en trois dimensions apporte une information qui peut être très utile à la compréhension de certains milieux. Ainsi, de même qu'observer des organes en 3-D plutôt que de simples coupes peut être déterminant en imagerie médicale, pouvoir exploiter diérents angles de vues de scènes complexes telles que des zones densément peuplées ou certains reliefs naturels trouve des intérêts dans de nombreuses applications.

Ainsi l'étude des activités humaines ou le suivi de déformation du sol ou des bâtiments peuvent grandement bénécier de ces techniques.

De nombreuses modalités existent pour obtenir une visualisation 3-D d'une scène.

Utiliser plusieurs images du même objet prise sous un point de vue diérent est une approche classique. De bons résultats sont d'ailleurs obtenus de cette façon à partir d'images naturelles. Ainsi des algorithmes permettant d'obtenir des modèles 3-D de zones urbaines avec une résolution centimetrique sont maintenant couramment utilisés. Une technique d'imagerie qui a fait ses preuves depuis de nombreuses années est l'imagerie par RAdio Detection and Ranging (RADAR) à synthèse d'ouverture ou SAR.

Bien que la résolution atteignable est généralement en deçà de ce qui peut être obtenu en optique, cette technique d'acquisition présente d'autres avantages. Ainsi, le RADAR étant actif, des images peuvent être prises indépendamment de la luminosité ambiante.

De plus les ondes émises ne sont pas stoppées par la couverture nuageuse ce qui peut se révéler utile lorsque la zone imagée est soumise à des intempéries ou pour observer un volcan actif. L'imagerie SAR peut donc être un instrument de prédilection pour observer des zones à l'épreuve des éléments.

Les images SAR correspondent à l'addition cohérente des ondes rééchies par les diérents réecteurs présents au sol. Chaque pixel possède une information de phase qui peut être liée à la distribution volumique des réecteurs qui se projettent dans ce pixel. La tomographie SAR exploite la diérence de phase entre plusieurs images recalées pour extraire la localisation des objets au sol. Cette technique a été utilisée avec succès pour reconstruire des volume de végétation [START_REF] Reigber | First demonstration of airborne SAR tomography using multibaseline L-band data[END_REF] [START_REF] Huang | Under foliage object imaging using SAR tomography and polarimetric spectral estimators[END_REF] ou de glace [START_REF] Tebaldini | Imaging the internal structure of an alpine glacier via l-band airborne sar tomography[END_REF]. Des paysages plus complexes comme des XIV CONTENTS zones urbaines denses ont également pu être reconstruits avec cette technique. Toutefois l'intensité très variable des signaux rééchis par les bâtiments ainsi que la quasi omniprésence de phénomènes de layovers (diérents réecteurs se projetant dans la même case) rend l'utilisation de la tomographie SAR plus compliquée.

Contributions Si les zones urbaines densément construites peuvent être diciles à analyser, elles présentent toutefois la propriété d'être très structurées. Ainsi, il est raisonnable de s'attendre à ce que des murs soient droits, des toits et des rues plats et horizontaux. Cette structuration peut être essentielle pour améliorer les reconstructions tomographiques. Deux approches peuvent être envisagées pour en tirer partie :

Analyser les redondances dans les images SAR et se servir des pixels les plus simmilaires pour améliorer l'estimation de signal reprojeté.

Inverser le signal SAR en s'assurant que la reconstruction réponde à des contraintes structurelles.

Le premier point a été exploité pour la première contribution présentée dans ces travaux :

Contribution (1): Critère de similarité pour large pile d'images SAR L'algorithme Non Local -SAR (NL-SAR) est une méthode ecace pour débruiter des images SAR ainsi que des piles d'images interférométriques et/ou paramétriques. La matrice de covariance du signal est fondamentale pour la plupart des estimateurs utilisés en tomographie SAR . Une meilleure estimation de celle-ce pourrait donc aboutir à une meilleur reconstruction tomographique.

Lorsque la dimension devient trop importante, les performances de NL-SAR ont cependant tendancent à diminuer. Ici nous proposons un nouveau critère robuste à l'augmentation de la dimension en modiant le critère de similarité entre pixel utilisé par l'algorithme.

Ce nouveau critère est eectivement plus robuste à l'augmentation de la dimension et permet une sélection prenant explicitement en compte les réecteurs les plus puissants pour chaque pixel. Cette approche reste toutefois limité car les échantillons similaires peuvent s'avérer soit trop peu nombreux soit contraints à un ensemble depoints réparti de façon isométrique autour du pixel d'intérêt. Les résultats tomographique ne sont donc pas sensiblement améliorés.

Le second axe de recherche envisagé peut sembler plus prometteur pour les paysages densément construits. En eet, l'approche Compressive Sensing (CS) [START_REF] Budillon | Three-Dimensional SAR Focusing From Multipass Signals Using Compressive Sampling[END_REF] (Zhu and Bamler, 2010a) permettant d'inverser chaque pixel sous contrainte de parcimonie s'est imposée ces dernières années comme l'un des meilleurs estimateurs pour la reconstruction de milieu urbain par tomographie SAR. CS s'inscrit donc com-CONTENTS XV plètement dans la deuxième approche. De plus, comme il ne dépend pas de la matrice de covariance du signal il est adapté à des milieux hétérogènes.

Contribution (2): Estimation tomographique sous contraintes exprimées en géométrie sol L'approche CS classique en tomographie SAR n'utilise pas d'autres a priori que la parcimonie. Nous proposons dans ces travaux d'aller plus loin en introduisant une régularisation géométrique du volume reconstruit. Pour cela, l'estimation est faite directement sur un cube i.e. un ensemble de voxels et la fonction objectif à minimiser fait intervenir la norme quadratique du gradient dans les directions horizontales et verticale.

Cette extension naturelle de l'approche CS classique permet de favoriser des distributions de réectivité qui soient structurées et de réduire les erreurs de reconstruction.

Le cadre proposé dans ce nouvel algorithme pour l'introduction de contrainte spatiales est relativement simple et peut être repris pour ajouter d'autres a priori.

L'analyse de reconstruction 3-D par tomographie SAR n'est pas toujours facile.

Ainsi, même après avoir projeté la scène dans le repère sol, la forte dynamique des voxels, la présence de lobes secondaires et/ou d'erreurs rend l'interprétation dicile.

Pouvoir situer la surface urbaine i.e. l'ensemble des objets (murs, toits, routes, etc) ayant rééchi l'onde incidente n'est donc pas exactement une tâche triviale. Exploiting multiple views of the same object to obtain a 3-D interpretation has been a research topic in many elds and with quantities of applications. When the observed object is, in fact, the earth, many dierent imaging techniques can be consider depending on the soil, the vegetation, the relief, the luminosity... 3-D reconstruction obtained by exploiting multiple optical images can be computed with a great precision providing centimeter resolution. Complex landscapes such as, for instance, dense urban areas can then be processed to retrieve the 3-D distribution of the scene. The 3-D rendering of urban areas has many various applications such as city management, architecture, crisis management, building and ground deformation monitoring or demography expansion. SAR images are an other well known technique to observe the earth from the sky or space. The resolution of SAR images is generally worse than what can be obtained with optical sensors although modern sensors can provide images of few centimeters resolution. The SAR imaging technique however provides other advantages. The sensor being active, images can be taken even by night or through clouds. Depending on the operating bandwidth, the wave emitted by the RADAR can penetrate some environments (vegetation, ice) and provide information about their structure or the presence of underneath objects [START_REF] Reigber | First demonstration of airborne SAR tomography using multibaseline L-band data[END_REF] [START_REF] Huang | Under foliage object imaging using SAR tomography and polarimetric spectral estimators[END_REF] [START_REF] Tebaldini | Imaging the internal structure of an alpine glacier via l-band airborne sar tomography[END_REF]. Finally, as the obtained images are composed of the coherent summation of back-scattered signals, the pixels have a phase information that can be linked to their 3-D distribution.

Exploiting phase dierences between SAR images to obtain the height of the scatterers is an important research eld that backs to the 1980. SAR interferometry consists in estimating the height of the main back-scattering object for each pixel exploiting the phase dierence between a pair of images. As more data with higher resolution became available, using more than two images became possible with the advantages to avoid fringes unwrapping, drastically lowering ambiguities or even suppressing them.

Of course with more accurate information, not only one but several scatterers height could be evaluated. Depending on the imaged scene even a 3-D reectivity distribution could be observed. SAR tomography is the technique consisting in evaluating the scatterers reectivity in the 3-D space using a stack of co-registered SAR images. Many results have been shown on uniform areas with a volume distribution of scatterers such as ice or forest [START_REF] Reigber | First demonstration of airborne SAR tomography using multibaseline L-band data[END_REF] [START_REF] Lombardini | Adaptive spectral estimation for multibaseline sar tomography with airborne l-band data[END_REF] [START_REF] Guillaso | Scatterer characterisation using polarimetric sar tomography[END_REF] [START_REF] Huang | Under foliage object imaging using SAR tomography and polarimetric spectral estimators[END_REF] [START_REF] Tebaldini | Multibaseline Polarimetric SAR Tomography of a Boreal Forest at P-and L-Bands[END_REF]. Results on urban areas also present accurate scatterers localization. Dense urban scene are however hard to process due to the very heterogeneous back-scattering mechanisms and the importance of layover phenomena. However, urban areas are very structured and geometrical priors can be used to enhance the tomographic reconstruction. In the past years algorithms exploiting the sparsity of the signal such as MUSIC [START_REF] Schmidt | Multiple emitter location and signal parameter estimation[END_REF] or the well-known (Zhu and Bamler, 2010a) [START_REF] Budillon | Three-Dimensional SAR Focusing From Multipass Signals Using Compressive Sampling[END_REF] have led to promising results with many scatterers being retrieved. Nevertheless, the redundancy of the buildings and geometrical shapes were never used to perform the tomographic inversion. The subject of this PhD is to explore dierent approaches to exploit this structural information for the SAR tomographic reconstruction of urban areas.

Contributions

At least two strategies can be tried to exploit the particular geometry of urban areas: Analyze the redundancies in the SAR images to select pixels presenting the same behavior.

Perform the tomographic inversion under some priors on the reconstructed result.

The rst approach is motivated by the good performances of non-local denoising algorithms such as NL-SAR [START_REF] Deledalle | NL-SAR: A Unied Nonlocal Framework for Resolution-Preserving (Pol)(In)SAR Denoising[END_REF]. By exploiting similarities between pixels and patches in the images, these approaches manage to obtain a good estimation of the covariance matrix. This parameter is central for most of the estimators used in SAR tomography

Contribution (1): Similarity criterion for large stack of SAR images

We present a similarity criterion robust to the increase of SAR images in the tomographic stack. This new criterion is based on a decomposition of the signal in a deterministic, stable part and a stochastic one. This criterion is tested over dierent congurations of dense urban areas.

The second approach may be actually more suited to complex landscapes such as dense urban areas. Indeed the well-known (Zhu and Bamler, 2010a) [START_REF] Budillon | Three-Dimensional SAR Focusing From Multipass Signals Using Compressive Sampling[END_REF] approaches that estimate the reectivity for each pixel as the sparsest projected 1.2. CONTRIBUTIONS 3 vector lead to correct results for urban areas. This approach may be consider as one of the top state-of-the-art method for urban landscape and may be the best suited for dense concentration of buildings. Contrary to most of the other methods used in SAR tomography, Compressive Sensing (CS) doesn't need any estimation of the covariance matrix and thus is not impacted by inaccurate parameter estimation. The main issue to add geometrical priors to the reconstruction is to deal with the distortions induced by the active side-looking sensor. Once it is done, dierent priors and optimization strategies can be consider.

Contribution (2): Ground coordinate based geometrical priors for SAR tomography

We present a new algorithm allowing to perform the tomographic reconstruction of a scene under structural priors. Its spirit is very close to CS as it also takes benet from the sparsity. However, we go one step further by constructing the projection operator as a mapping from ground to radar coordinates. This allows us to express geometrical priors in the scene frame and favor smooth vertical and horizontal distribution of scatterers.

Having a framework in which the estimated reectivity is described in its natural coordinate system opens the door to many applications. From the 3-D reectivity one may want to retrieve the actual urban surface on which are located the scatterers. This information can moreover be used to rene the quality of the tomographic inversion.

Contribution (3): Graph-Cut based urban surface segmentation from SAR tomographic reconstruction

We present a segmentation algorithm to extract the urban surface from 3-D tomographic reconstruction expressed in ground geometry. The priors on the surface are more sophisticated than previously and allow to select shapes corresponding to the following considerations:

The surface is located near high intensity voxels.

The buildings are described by smooth, uniform structures. The urban surface is then expected to correspond to a Total Variation (TV) minimization along the vertical and horizontal directions.

The scatterers on buildings are supposed to be opaque for the electromagnetic wave. The surface is then expected to go through each ray from the sensor to the scene exactly once.

To avoid distortions in the vertical structures, it is impossible for the surface to intersect multiple time a vertical line The cost function associated to these hypotheses is highly non-convex. However by discretizing the problem, it can be formulated as minimum cut / maximum ow search on a well constructed graph. The method proves to be ecient on real data and performs particularly well when associated to CS like tomographic reconstructions.

Contribution (4): Alternated 3-D reconstruction and surface segmentation

The estimated urban surface can be used to tune more accurately the parameters of the estimators used to perform the inversion. This strategy is used to design a new algorithm that performs iteratively the reconstruction in ground coordinates and segment the urban surface from it.

Organization of the manuscript

This manuscript is divided in two parts. First basic concepts of SAR imaging and stateof-the-art methods for SAR tomography are presented. In chapter 2, the SAR images acquisition is summarized. In chapter 3, the dierent statistical models used in SAR tomography are presented. Finally, the chapter 4 details dierent algorithms used in SAR tomography.

The second part is centered on the dierent developed methods during this PhD.

Chapter 5 is a small introduction summarizing the presented research axes in this part.

Chapter 7 is focused on covariance matrix estimation. Chapter 8 presents the proposed tomographic reconstruction method integrating 3-D geometrical priors (submitted in (Rambour et al., 2018a)). Finally, chapter 9 corresponds details the graph-cut based urban surface segmentation and the AlteRnatEd 3-D REconstruction and Surface Segmentation (REDRESS) algorithm (submitted in (Rambour et al., 2018b)).

Part I SAR tomography

Chapter 2 Synthetic Aperture Radar

Introduction

This chapter is an introduction to the SAR techniques. First, the bases of the SAR images acquisition are described and the principle of the aperture synthesis are given.

Then, simple multibaseline phase model and foundation of SAR tomography are described.

Acquisition

SAR is a coherent, microwave imaging technique that produces an image of the spatially distributed complex reectivity of the scene with a resolution depending on the RADAR characteristics. As many other active imaging techniques, the principle is to transmit an electromagnetic wave in a direction called the Line of Sight (LOS) to an object of interest. The time from the transmission to the reception of the wave gives an indication on the distance between the object and the transmitter. In SAR imagery, the sensor is embedded on a moving platform (a plane or a satellite) and illuminates the ground according to an incidence angle θ and a squint angle α. The geometry of acquisition for a single antenna is shown in Fig. 2.1.

Range resolution: After its transmission, the wave reaches a back-scattering object at a given time delay, it is then received after propagating back to the sensor. The time between the transmission and the reception is then:

∆t = 2d c (2.1)
with d the distance from the RADAR to the ground and c the speed of light. The RADAR transmits short pulses towards the scene that then back-scatters part of that incoming signal depending on its composition and geometry. To be able to separate two scatterers in the range direction, their echoes must be temporally separated by more than the duration τ between two pulses. The slant range resolution is then dened as:

δ r 0 = cτ 2 (2.2)
To obtain a good resolution in range one should then transmit close pulses. However, to obtain a good Signal-to-noise ratio (SNR), the pulses should be long to transmit more energy. To still achieve a good resolution, a frequency modulated pulse is transmitted.

With a chirp of central frequency f 0 and frequency excursion [f 0 -Bc 2 , f 0 + Bc 2 ], the resolution is, after signal processing (see section 2.3):

δ r = c 2B c (2.3)
B c here is the chirp bandwidth dened as B c = KT with T the transmission duration and K the rate of frequency change or chirpyness. Therefore, a high bandwidth leads to a well resolved signal. Due to the incidence angle of the sensor, the ground range resolution is (under the assumption that the squint angle α is negligible and the ground at): Azimuth resolution: The angular spread given by a conventional RADAR is related to the size of the antenna L and the wavelength λ by:

δ gr = c 2B c sin(θ) (2.4)
γ = λ L , (2.5)
The resolution in azimuth for a scene located at a distance R from the antenna is then:

δ az 0 = γR = λR L , (2.6)
For SAR system in space, the order of magnitude of R can be 10 5 m and for X-band systems λ is a few centimeters which leads to an order of magnitude of the numerator in (2.6) of 10 3 m to achieve a metric resolution. Of course, such a dimension for an antenna mounted on a satellite is not realistic. To improve this low resolution, the SAR technique consists of coherently combine the echoes of the scatterers as they are seen from successive azimuth positions of the RADAR. This creates a synthetic antenna of length L sa as illustrated in Fig. 2.2. The new angular spread is:

γ sa = λ 2L sa , (2.7)
where the factor 2 in the denominator comes from the return trip of the transmitted wave. The length of the synthetic antenna correspond to the displacement of the sensor for which a given point back-scatters a signal to the antenna. L sa is then given by the size of the antenna footprint L sa = δ az 0 . This gives the achievable resolution in azimuth:

δ az = γ sa R = L 2 (2.8)
It is worth to notice that the resolution in azimuth does not depend on the distance between the RADAR and the ground. Moreover the smaller the antenna the larger its footprint and then the better the resolution.

Basics of SAR raw data processing

Range compression: Let's consider that the emitted signal consists in a pulse dened as:

∀t ∈ [- T 2 , T 2 ], s e (t) = exp 2jπ(f 0 t + K 2 t 2 )
(2.9)

To better identify real-valued and complex-valued variables which will be useful in the next chapters, all complex-valued variables are underline. The received echo from a scatterer located at a distance R to the sensor is then

s r = s e t - 2R c (2.10)
To exploit the phase modulation and achieve the resolution stated in (2.3) a matched lter is applied to the signal. This matched lter is nothing else than the transmitted pulse itself and one can verify that:

s * e (-t) s e (t) sinc Bc (t) exp 2jπf 0 t (2.11) where sinc Bc (x) = sin(πBct) πBct with B c = KT . After convolution by the pulse replica, the received signal becomes then:

s r = sinc Bc t - 2R c exp -2jπ( 2R λ -f 0 t) (2.12)
After centering of the spectrum, the received signal is then:

s r = sinc Bc t - 2R c exp -4jπR λ (2.13)
The achieved resolution is given by the width of the sinc function and corresponds then to the one stated in (2.3).

Azimuth compression:

The azimuth compression step is analog to the previous one.

The phase modulation here is introduced by the motion of the sensor. From 2.3, the distance R from the antenna at position x and a point on the ground is given by:

R = R 2 0 + (x -x 0 ) 2 (2.14)
R 0 is the distance from a reference point along the track to the target. As the distance from the sensor to the target is generally higher than the synthetic antenna, a second order approximation gives:

R = R 0 + (x -x 0 ) 2 2R 0 (2.15)
The received signal from the antenna in x is then (by a change of origin):

s r = sinc Bc t - 2R c exp -4jπR 0 λ exp jπK x x 2 (2.16) with K x = 2 λR 0
From the relation (2.14), the convolution of the signal by its replica gives then:

s r = sinc Bc t - 2R c sinc Bx (x) exp -4jπR 0 λ (2.17)
The bandwidth is dened as B x = K x L sa . The azimuth resolution is then:

δ az = 1 B x = λR 0 2L sa = L 2 (2.18)
These two processes describe the basic of the range and azimuth compression for a point like signal. In practice advanced techniques have to be employed to treat the SAR images depending on the acquisition system.

SLC SAR images

The complex value v(x, r) obtained at pixel coordinates (x, r) after SAR synthesis of the SAR image corresponds to the convolution of the complex 3-D reectivity u(x, y, z) with the Point Spread Function (PSF) of the sensor [START_REF] Reigber | First demonstration of airborne SAR tomography using multibaseline L-band data[END_REF][START_REF] Fornaro | Three-dimensional focusing with multipass SAR data[END_REF]:

v(x, r) = f (x -x , r -ρ y,z )u(x, y, z). exp - 4jπ λ ρ y,z + ja n (x, y, z) dx dydz + (2.19)
To avoid complicated notations x and azimuth axis coincide. Here x and r stand for the coordinates of the focused data, ρ y,z is the distance between the sensor and a scatterer at a position (y, z) for a given x. a n (x, y, z) is the phase shift corresponding to the Atmospheric Phase Screen (APS) contribution. f corresponds to the PSF and depends of the sensors and the processing of the data. is a white additive Gaussian noise and models the thermal noise. The formula (2.19) is illustrated Fig.

v(x, r) = (y,z)∈∆r u(x, y, z)exp - 4jπ λ ρ y,z + ja(x, y, z) dxdydz + (2.20)
The r-th radar resolution cell is dened by ∆ r = {(y, z) | r -δ r /2 ≤ ρ y,z ≤ r + δ r /2 and z ≤ z max (y)} with δ r the step in range and z max (y) the maximum height for a point at (x, y) location to be illuminated by the radar. x, y and z stand for the coordinates of a reector in 3-D space and x, r and h its coordinates in RADAR space.

The SAR range imaging system is the cause of geometrical distortions in the observed images. Indeed, due to the incidence angle of the RADAR, the structured objects presenting an angle β between their normal and the LOS inferior to 0 • are stretched toward the sensor whereas the ones for which β > 0 • are compressed in the image. These eects are illustrated in Fig. 2.5. Objects presenting an angle β > 0 • are then projected backwards in the range direction as illustrated in Fig. 2.6 where a building is projected in decreasing range positions as its height increases.

As the emitted electromagnetic waves are coherent, the total backscattering amplitude is aected by a speckle eect. This phenomenon is produced by the coherent summation of all the contributions in the radar cells that may be constructive or destructive. Depending on the application and the resolution, the speckle may be considered as a noise. In this document, it is rather described as a phenomenon corresponding to complex back-scattering mechanisms.

Phase model and calibration

The topography of a scene can be obtained through the analysis of the phase dierence in a stack of N co-registered SAR images. Even if theoretically the 3-D location of the scatterers could be obtained directly from the dierent ρ n;y,z

for n ∈ {1, • • • , N },
it is generally more convenient to work with the optical paths dierences ie. the in- terferograms. This way, the height or elevation of scatterers can be seen as directly proportional to the phase. To be able to use classical phase calibration algorithms and be in the same framework as the other state-of-the-art tomographic estimators, the SAR images are used as interferometric pile through this manuscript. Starting from the set of SAR images, the interferometric SAR stack is built by removing the phase of the master image to every SAR image. For simplicity the rst image of the stack is chosen as the master one. The value of a pixel v n (x, r) from the n th is then:

v n (x, r) = (y,z)∈∆r u(x, y, z)exp -jϕ n;x,y,z dydz + (2.21)
where the interferometric phase is dened as:

ϕ n;x,y,z = - 4jπ λ (ρ n;y,z -ρ 1;y,z ) + j(a n (x, y, z) -a 1 (x, y, z)) (2.22)
The optical path dierence is generally split into two terms relative to the elevation h(y, z) and the range variation of the scatterer ∆ρ 1;y,z from the center of the RADAR cell (D'Aria et al., 2010):

ρ n;y,z -ρ 1;y,z = b n h(y, z) r + b n ∆ρ 1;y,z r tan(θ) (2.23)
The distance dierence in (2.23) is obtained using a rst order approximation whereas a second order one was used for the azimuth compression in (2.15). This means that a phase factor depending on the square of the distance variation will corrupt the estimation of u(x, y, z) and needs to be compensated via post-processing. However, here as in the vast majority of the tomographic applications, we are mostly interested in the amplitude of the reectivity and the second order term is left in the phase of the reectivity.

The analysis of interferometric data has a longer history than SAR tomography and many phase models have been proposed for various applications and the scene. Due to the high concentration of dihedral and trihedral structures occurring in urban areas, many bright point-like echoes appear on the corresponding SAR images. These points are generally very stable between the acquisitions and the corresponding pixels are well described by the PS model [START_REF] D'aria | SAR Calibration Aided by Permanent Scatterers[END_REF]. The PS interferometric phase model for a point located in (y, z) is built as a linear combination of its elevation, its potential motion along the LOS direction and linear perturbation along the azimuth and range direction:

ϕ n;x,y,z = a 0 + k (n) x x + k (n) r r + ξ n h(y, z) + t n v + e (2.24)
where a 0 is a constant phase oset for all the images, k λr is the spatial elevation frequency associated to the sampling and b n is the baseline n as illustrated in Fig. 2.7, t n the time between the rst acquisition and the current one and v is the unknown scatterer slant range velocity. This phase model makes the assumption that only one scatterer is present in each cell. This strong hypothesis is true if the pixel stability is high enough so that almost all the received echoes come from the same scatterer. The system (2.24) is composed of six unknowns and has to be estimated for the whole stack of SAR images. As there is no knowledge of the unwrapped phase or on the number of scatterers in the PS-like cells, the problem formulated in 2.24 is highly non-linear and cannot be inverted from a set of L bright stable pixels. Dierent algorithms have been proposed these last two decades to estimate the phase model coecients and address the problem of optimizing the phase model given by:

(â 0 , k(n) x , k(n) r , ĥ, tn , ê) = argmax N n L l=1 v n (x l , r l ) exp(-jϕ n;x,y,z ) (2.25)
where ϕ n;x,y,z is given in (2.24). The spaceborne data presented in this manuscript were processed using a similar approach as the method described in (D'Aria et al., 2010) and the optimization was performed alternatively on the linear APS components and the elevation. Moreover, even if this hypothesis will not always hold, the scatterers are supposed to remain still between the acquisitions. After calibration of the data, the complex value of a pixel in the SAR interferometric images is dened by: v n (x, r) = (y,z)∈∆r u n (x, y, z)exp -jξ n h(y, z) dydz +

(2.26)

SAR tomography

SAR tomography is the extension of the 2-D SAR imaging to three dimensions. As conventional 2-D SAR imaging uses a synthetic aperture in the azimuth direction, 3-D SAR imaging is performed by a synthetic aperture in the elevation direction by collecting several images from parallel tracks. Using multibaseline interferometry techniques on the well-calibrated SAR images stack, it is possible to retrieve the localization of the scatterers in the third dimension. This approach allows separating scatterers mapped in the same resolution cell, which is likely to happen on the dense urban area due to the layover phenomenon.

A SAR tomographic stack consists in N SAR SLC images perfectly co-registered.

Each SAR image of the stack corresponds to a slightly dierent trajectory of the sensor over the scene. We consider all images to have been co-registered with respect to a master image in a preprocessing step. Each image is acquired from a slightly dierent angle at each pass of the sensor. This angular diversity induces a dierent distance ρ n;y,z to each antenna thus a dierent phase shift which can be exploited to retrieve the 3-D location of the scatterers. After phase calibration and under a far-range approximation, it is generally the dierent baselines b n that are used to characterize the received echoes rather than the dierence in the wave propagation.

Through the equation of the SAR interferometric pixels (2.26), it can be seen that the value of v n (x, r) is given by the Fourier transform of the 3-D reectivity function. A common way to retrieve u(x, y, z) is then to apply a Discrete Inverse Fourier Transform (DIFT) to the vector v(x, r)

= v 1 (x, r), • • • , v N (x, r)
T collecting the SLC values for all the pixels located at the position (x, r) in the dierent images. However, doing so, the obtained estimation of the reectivity is described in the RADAR reference frame (x, r, h). The interpolation of the reectivity function is then given by:

û(x, r, h) = a(r, h) H v(x, r) (2.27)
where the vector a(x, r) is the steering vector associated to location (r, h) and dened as:

a(r, h) = exp(-jξ 1 h), • • • , exp(-jξ N h) T (2.28)
This approach only leads to poor estimation of the reectivity as the resolution δ h is given by the maximal orthogonal baseline ∆b:

δ h = λR 2∆b
(2.29)

The resolution associated to a maximal orthogonal baseline of 100m with a wavelength λ = 0.0311m would then be around 55m in a spaceborn conguration (R 10 5 m).

Chapter 3

Multibaseline signal models

This chapter presents the signal statistical models used in SAR tomography.

The rst model corresponding to distributed sources such as rugged surfaces (oor,street,rooftops) or volumes (forest,ice) is generic whereas the second one is a derivation of the rst one for strong stable scatterers.

Sensor array signal model

From now on to the end of the rst part of this manuscript when no spatial regularization is involved, the subscripts identifying the location of the cells are dropped and all the equations are relative to one radar cell. In the case where dierent pixels are involved as in the empirical covariance matrix computation, they will receive a discrete subscript but no 2-D or 3-D spatial information.

Under the hypothesis that there is a nite number D of point like scatterers in the observed cell, we dene the vector h

= h 1 • • • h D T
containing all the elevations of the dierent scatterers in the observed cell. For 1 ≤ d ≤ D, the vector

u d = u 1,d • • • u N,d
T collects the complex reectivity of one scatterer for each acquisition. After discretization of (2.21), we can express v n , the SLC value of the pixel corresponding to the n th track as the sum of the complex signals back-scattered by each of the D scatterers:

v n = d u n,d exp(-jξ n h d ) + , (3.1)
Two models can be considered depending on the behavior of the scatterers reectivity from one image to another. The rst one is general and does not make any assumption on the scatterers correlation. The second one is adapted to very bright stable scatterers or campaign where the acquisitions were close enough in time to provide coherent images with no temporal decorrelation. The dierent models and the type of scatterer distribution they are considering are illustrated in Fig. 3.1. any assumption on the scatterers correlation, the vector v collecting the signals received by all the antennas is written:

LOS

v = d u d a(h d ) + (3.2)
with being the Schur-Hadamard product (elementwise multiplication). a(h d ) is the steering vector for the elevation h d :

a(h d ) = exp(-jξ 1 h d ) • • • exp(-jξ N h d ) T (3.3)
The noise is supposed to be a stationary, independent, white Gaussian noise with power σ 2 n and covariance matrix σ 2 n I N ∈ C N , (I N being the matrix identity of size N ). The vectors u d are modeled as the product of a random vector x d times their amplitude:

R = E{vv h } = d τ d C d a(h d )a(h d ) H + σ 2 n I N (3.5)
Another form of the previous equation without the Schur-Hadamard product is:

R = d τ d L d C d L h d + σ 2 n I N (3.6)
where the matrix L d ∈ C N is the diagonal matrix associated to a(h d ).

Conditional (determinist) signal model

The second model is called conditional or deterministic by simplication even though it is not associated to a fully deterministic model since the additive noise is still present. As the vectors u d are now considered constant, the former equation (3.2) can be written as:

v = d u d a(h d ) + = A(h)u + (3.7)
A(h) ∈ C N ×D is the so-called steering matrix and its d-th column [A(h)] d corresponds to the steering vector a(h d ) associated to the elevation of the d-th scatterer:

[A(h)] d = a(h d ) = exp(-jξ 1 h d ) • • • exp(-jξ N h d ) T (3.8) u = u 1 • • • u D T
is the vector collecting the scatterers reectivity. As the only nondeterministic part of equation (3.7) is due to the additive gaussian noise, the covariance matrix is then σ 2 n I N . Many spectral estimators however are dened on the non-central

R = E{vv H } = A(h)SA(h) H + σ 2 n I N (3.9)
with S = uu H ∈ C D×D being the signal covariance matrix.

Hybrid signal model

The hybrid model is the combination of the two preceding ones corresponding to a mixture of D d fully coherent and D u decorrelated signals:

v = Du d=1 u d a(h d ) + Du+D d d=Du+1 u d a(h d ) + (3.10)
The non-central covariance matrix is then given by:

R = Du d=1 τ d C d a(h d )a(h d ) H + Du+D d d=Du+1 A(h)S A(h) + (3.11)

Phase Model

It is possible to link the conditional model to the phase model given in section 2.5 when the modulus of the scatterers is temporally stable. Under this consideration, the equation (3.2) can then be written as in the conditional model with the phase model being included in the steering vectors. As stated in the section 2.5, phase models on urban areas may include potential movement of the scatterers. In this case the signal model is then:

v n = d u d exp(-jξ n h d + η n v d ) + (3.12) And then, v = A(h, v) + (3.13)
where η n = 4πtn λ is the so-called "velocity frequency" and t n is the time interval between the acquisitions of the master and the n th image. The scatterers slant range velocities

correspond to the vector v = v 1 • • • v D T . The matrix A(h, v
) is again dene as the concatenation of the steering vectors:

[A(h, v)] d = a(h d , v d ) = exp(-j(ξ 1 h d + η 1 v d )) • • • exp(-jξ N h d + η N v d ) T (3.14)
Here, only one extension of the phase model to include the displacement of the objects is considered. Other models have been proposed including for instance the temperature (Budillon et al., 2017b) [START_REF] Weissgerber | Joint measurement of height and deformation by radar Interferometry: the example of the Eiel Tower[END_REF] or the clutter decorrelation [START_REF] Aghababaee | Multiple Scatterers Detection Based on Signal Correlation Exploitation In Urban SAR Tomography[END_REF]. As stated before, only the non-moving scatterers model is used in this manuscript. The extension described in (3.12) is conceptually easy to add but may enlarge an already high number of unknowns when working on many pixels at once.

Conclusion on urban signal models

Introduction

The aim of SAR tomography is to retrieve an estimation of the 3-D reectivity of a given scene. After co-registration of the SAR tomographic stack this corresponds to the inversion of equation (2.26):

v n (x, r) = (y,z)∈∆r u n (x, y, z)exp -jξ n h(y, z) dydz + (2.26)
If the phases are correctly calibrated and the APS removed, in the ideal case of equispaced trajectories, focusing in the direction orthogonal to the line of sight can be simply performed by application of the inverse discrete Fourier transform:

û = A( h) H v (4.1)
where the vector h stands for the sampling along the elevation direction. The matrix A( h) is then the inverse Discrete Fourier Transform (DFT) matrix. The resolution of this focusing is inversely proportional to the maximal orthogonal baseline. Increasing the number of tracks within this maximal orthogonal baseline improves the sampling in Fourier domain, hence it reduces height ambiguities. However, the vertical resolution is generally much worse compared to the resolution in azimuth and range directions. Moreover, the baselines are generally irregularly distributed which produces side-lobes higher than expected and degrades the interpretation of the reconstructed volume. Several spectral super-resolution techniques have been introduced to overcome these phenomena. This chapter details the two main categories of approaches used in SAR tomography. The rst category is composed of the spectral analysis estimators coming mainly from the Direction Of Arrival (DOA) literature, they exploit the covariance matrix of the received signal. The second category corresponds to the more recent CS approaches which achieve super-resolution without resorting to an estimation of the covariance matrix and are more suited to sparse and heterogeneous areas. In the previous section two signals models were described. Only the conditional model is considered in the following tomographic approaches as it is well suited to urban areas. Some of them are nevertheless robust to decorrelation phenomena.

Except for the unconditional Maximum Likelihood, all the presented method are used in the presented works. The methods are illustrated on simulations and on a real slice. The simulations are built following the conditional model with the scatterer having their amplitude xed for all the images. The SNR corresponding to the additive noise is of 1.4 dB. The amplitude of the scatterers is set equal to 1 and the PSF function is taken as a sinc function inducing a mixture between adjacent cell. Simulated slice presents scatterers well separated in the range direction. Real data are however hard to model accurately due to various decorrelation mechanisms. It is then useful to look at results obtained from both a simple and a real scenario.

Spectral Analysis Techniques

Several estimation algorithm from the spectral analysis eld that are used for SAR tomography are presented in this section. They allow estimating the parameters of a multibaseline SAR signal even if it is corrupted by a speckle eect.

Beamforming techniques

Beamforming is an important technique in signal array processing to estimate a DOA. This approach has been used in many dierent applications such as RADAR, SONAR, wireless telecommunications or medical echography. Beamforming techniques consist in nding a Finite Impulse Response (FIR) lter whose output is maximal when the received signal originates from a source located in a direction of interest.

Beamforming lters are the rst algorithms used in SAR images to perform the unmixing of scatterers [START_REF] Homer | High resolution 3-d sar via multi-baseline interferometry[END_REF]. Many beamforming estimators are derived in the DOA literature and the ones used in SAR tomography are designed to suppress the noise under a directional constraint. The problem is then to nd the lter f (h) that maximizes the signal to noise ratio i.e., that minimizes the noise level for a given output level:

min f E{|f (h) H | 2 } s.t f (h) H a(h) = 1 (4.2)
The solution f of the directional beamforming problem are also called Minimum Variance Distortionless Response (MVDR) lters. The resolution of this linearly constrained quadratic minimization uses the following results (P. [START_REF] Stoica | Introduction to Spectral Analysis[END_REF]

: Theorem 4.1 Let D ∈ C n be a positive denite matrix, X ∈ C n×m , G ∈ C n×k and C ∈ C m×k with k ≤ n. Then the following minimization problem: min X X H DX s.t X H G = C
has a unique solution given by: X

0 = D -1 G G H D -1 G -1 C H
Conventional Beamforming: According to equation (3.9), the term to minimize in the equation (4.2) can be rewriten using the noise covariance matrix R n = σ 2 n I N :

E{|f (h) H | 2 } = f (h) H R n f (h) = σ 2 n f H f (4.3)
The previous optimization problem then becomes:

min f f (h) H f (h) s.t f (h) H a(h) = 1 (4.4)
From the Theorem 4.1, the lter solution of the conventional beamforming problem is then:

f BF (h) = a(h) a(h) H a(h) = a(h) N (4.5)
The estimation of the power for the elevation h is given by the ltered signal:

P BF (h) = E{|f BF (h) H v| 2 } = a(h) H Ra(h) N 2 (4.6)
with R the non-centered covariance matrix of the signal R = E{vv H }. An image interpretation of the conventional beamforming for an additive white gaussian noise can be expressed through the empirical covariance matrix:

R = 1 L L l=1 v l v H l (4.7)
where L corresponds to the number of samples used for the estimation. The expression of the estimated power using this covariance matrix estimation is then:

P BF (h) = a(h) H R a(h) N 2 = 1 LN 2 L l=1 |a(h) H v l | 2 (4.8)
When using the empirical covariance matrix the power of the ltered signal can be seen as an averaging of the inverse DFT of the signal for dierent pixels. Capon Beamforming: The conventional beamforming lter is derived from a model of the noise distribution. However when dealing with SAR images stacks, specially when the dimension N becomes high, or when the echoes are produced by distributed scatterers, the unconditional model is more appropriate and decorrelation must be taken into account. When no prior on the noise is available, the beamforming lter can be designed to attenuate the power from every direction other than h under a unit gain constraint [START_REF] Capon | High-resolution frequency-wavenumber spectrum analysis[END_REF]:

min f f (h) H Rf (h) s.t f (h) H a(h) = 1 (4.9)
From the theorem 4.1, the solution is given by:

f C (h) = R -1 a(h) a(h) H R -1 a(h) (4.10)
The estimated power for the elevation h is then:

P C (h) = E{|f C (h) H v| 2 } = 1 a(h) H R -1 a(h) (4.11)
This method has empirically shown better resolution and lobe suppression than conventional beamforming (P. [START_REF] Stoica | Introduction to Spectral Analysis[END_REF] [START_REF] Gini | Multibaseline cross-track sar interferometry: a signal processing perspective[END_REF]. When the covariance matrix is poorly estimated this estimator may however presents a more hieratic behavior than conventional beamforming due to its dependency with 1/ R-1

.

MUSIC

The MUSIC estimator [START_REF] Schmidt | Multiple emitter location and signal parameter estimation[END_REF]) is designed to retrieve a nite number of DOA for a signal corrupted by an additive white noise. This estimator is easy to compute and is one of the rst to produce a sparse representation of the estimated signal. The MUSIC estimator derivation comes from the analysis of the non-centered covariance matrix subspaces.

The signal v consisting of the echoes produced by D < N scatterers following the conditional model has a covariance matrix structure given by the equation (3.9):

R = E{vv H } = A(h)SA(h) H + σ 2 n I N (3.9)
As R is Hermitian, there exists an orthogonal basis where this matrix is diagonal. The eigenvalues and eigenvectors of R are denoted by λ

N ≥ • • • ≥ λ 1 and e N , • • • , e 1 .
Under the assumption that the kernel of R is empty, the last N -D eigenvalues are The signal subspace:

E s = e 1 , • • • , e D and λ s = λ 1 , • • • , λ D (4.12)
The complementary subspace or so-called noise subspace:

E n = e N -D , • • • , e N and λ n = λ N -D , • • • , λ N = σ 2 n , • • • , σ 2 n (4.13)
Let D s and D n be the diagonal matrices built from λ s and λ n . The covariance matrix R can now be written as:

R = E s D s E H s + E n D n E H n (4.14)
As the subspaces are orthogonal, multiplying the previous equation by the matrix E n from the right side leads to

R E n = E n D n = σ 2 n E n (4.15)
From equation (3.9) the same matrix can also be described as:

R E n = σ 2 n E n = A(h)S A(h) H E n + σ 2 n E n (4.16)
Combining equations (4.15) and (4.16) and supposing that S is full rank leads to the following relation from which MUSIC is derived:

A(h) H E n = 0 (4.17)
The equation (4.17) states that the steering vectors associated to the back-scatterers are orthogonal to the noise subspace. The position of the scatterers can then be found by minimizing the projection:

ĥ = argmin h A(h) H E n E H n A(h) H (4.18)
The problem (4.18) is non-convex and depends on the unknown number of scatterers D. To avoid these diculties ĥ is generally estimated by looking at the maximums of the MUSIC Pseudo Spectrum (MUSIC-PS) dened as:

P M U (h) = 1 a(h) H E n E H n a(h) H (4.19)
For discrete scatterer distributions, the MUSIC algorithm performs better than beamforming estimators. Moreover, it is also robust to multiplicative noise for a small number of scatterers as the orthogonality between the steering vectors and the noise subspace is likely to be preserved. When the back-scattered echoes are highly correlated, the kernel of R can be non-empty which produces a degeneracy as the orthogonality relation is not exclusive to the signal steering vectors.

WSF

The WSF methods are in the same spirit as MUSIC but introduce the empirical covariance matrix subspace distribution in the estimators. At least two dierent WSF estimators can be derived. First, it can be observed that the orthogonality relation used in the MUSIC derivation implies that, if the rank of S is equal to D, then the range space of A(h) coincides with the one of E s . Then it can be stated that there is an unknown linear transformation T such that:

E s = A(h)T (4.20)
This linear relation and the orthogonal one dened in equation ( 4.18) can be used to derived two estimators based on the distance minimization between the matrix A(h) and the weighted subspace of the empirical covariance matrix R:

The Noise Subspace Fitting (NSF) estimator obtained by minimizing the following criterion:

||E H n A(h)|| 2 F W (4.21)
The Signal Subspace Fitting (SSF) estimator obtained by minimizing the following criterion:

||E s -A(h)T || 2 F W (4.22)
Where ||X|| 2 F M = tr(X M X T ) is the weighted Froebenius norm and W is an Hermitian positive semidenite weighting matrix. Consistent estimates of W are based on the empirical covariance matrix subspace distribution and allow to asymptotically reach the Cramer-Rao lower bound. They are given by [START_REF] Huang | 3D imaging for underfoliage targets using L-band Multi-Baseline PolInSAR Data and sparse estimation methods[END_REF]Viberg and Ottersten, 1991;[START_REF] Stoica | Maximum likelihood methods for direction-of-arrival estimation[END_REF]:

W SSF = D s -σ 2 n I N 2 D -1 s (4.23) W N SF = A(h) H E s W -1 SSF E H s A(h) -1 (4.24)
WSF techniques are supposed to provide high elevation resolution in SAR tomography.

However, the cost functions are non-convex and multimodal and thus hard to optimize.

A way proposed in (Viberg et al., 1991) is to choose the result given by a suboptimal minimization criterion such as MUSIC as an initialization.

ML

The ML techniques could refer to dierent approaches depending on the signal model used. For a signal following the conditional model (3.7), the negative log-likelihood is given by (up to a constant):

L(h, u) = ||v -A(h)u|| 2 2 (4.25)
When no assumption can be made on the signal covariance matrix, the negative loglikelihood becomes then:

L(h, u) = v -A(h)u H R -1 v -A(h)u = ||v -A(h)u|| 2 R -1 (4.26)
where ||x|| 2 M = x T M x is the weighted 2 norm. The minimization of L in (4.25) and (4.26) with respect to both h and u can be hard to achieve due to the phase dependency in h and the unknown number of scatterers. For given elevations, the conditional ML coincides with the Least-Square (LS) estimate of the signal. A LS approximation of the reectivity û can thus be computed after using a spatial component estimator such as MUSIC or WSF .

When multiple independent and identically distributed samples are available, the negative log-likelihood for the conditional model can be written using the empirical covariance matrix [START_REF] Stoica | Maximum likelihood methods for direction-of-arrival estimation[END_REF]:

L(h) = tr I N -A(h) A(h)A(h) H -1 A(h) R (4.27)
This function is highly nonlinear, multimodal and hard to minimize in a reasonable amount of time. Moreover, the conditional ML has been proven statistically less ecient for big number of samples than WSF techniques as it does not achieve Cramer-Rao lower bound [START_REF] Stoica | Maximum likelihood methods for direction-of-arrival estimation[END_REF].

M-RELAX

Multilook-RELAXation spectral estimator (M-RELAX) is an iterative algorithm that can be used to minimize (4.25) or (4.27) when multiple samples sharing the same scatterer elevation distribution are available [START_REF] Li | Ecient mixed-spectrum estimation with applications to target feature extraction[END_REF] [START_REF] Gini | Layover solution in multibaseline sar interferometry[END_REF]. At each step, the previously estimated sources are rened to account the newly detected one. When the number of scatterers is high, the correction step may be time consuming as it must be run multiple time until convergence. This procedure is described in the following algorithm: Initialization :

1: R ← 1 L l vv H 2: ĥ1 ← argmax h a(h) H Ra(h) 3: û1,l ← a( ĥ1 ) H v l N 4: d ← 2 5: while d < D do 6: while κ < κ do 7: k ← d 8: while k > 0 do 9: v(k) l ← vl - d i=1 i =k ûi,l a(h i ) 10: R(k) ← 1 L l v(k) l v(k)H l 11: ĥk ← argmax h a(h) H R(k) a(h) 12: ûk,l ← a( ĥk ) H v(k) l N 13: k ← k -1 14:
end while 15:

κ ← 1 L l ||v l - d i=1 û(i) l a(h (i) )||
ĥ(k) ← argmax h a H (h)v
The RELAX and M-RELAX algorithms without the correction process are forms of the CLEAN algorithm which presents no guaranty to converge to the global minimum but is much faster to optimize.

SPICE

The recent SPICE method [START_REF] Stoica | SPICE: A Sparse Covariance-Based Estimation Method for Array Processing[END_REF] is a fully non-parametric sparse algorithm based on the minimization of the covariance tting criteria:

f = ||R -1/2 R -R R-1/2 || 2 F (4.28)
The matrix R is structured as stated in the conditional model (3.9) except that sources are fully incoherent (S is diagonal) and the noise may be colored. The covariance matrix has then the form:

R = E{vv H } = D d=1 τ d a(h d )a H (h d ) +    σ 2 1 σ 2 N    (4.29)
The previous equation can be rewrite into a product of matrices:

R = EDE (4.30) with E = a(h 1 ), • • • , a(h D ), I N = e 1 , • • • , e D , e D+1 , • • • , e D+N and D =            τ 1 τ D σ 2 1 σ 2 N            =            d 1 d D d D+1 d N +D           
Under this covariance matrix structure assumption, a consistent estimate of the solution of the minimization of f is given by the following optimization problem:

min τ d 0 tr R 1 2 R -1 R 1 2 s.t. D+N d=1 w k d k (4.31) with w k = e H k R-1 e k N (4.32)
The minimization algorithm is presented in detail in [START_REF] Stoica | SPICE: A Sparse Covariance-Based Estimation Method for Array Processing[END_REF]. This algorithm achieves a good performance when the covariance matrix is correctly estimated and when the uncorrelated conditional model is respected. This method presents no parameter tuning and is relatively easy to optimize which makes it very promising. The signal model used however is even more restrictive than the conditional model.

Compressive Sensing

Except for conventional beamforming, the estimators detailed in the previous section may all achieve super resolution under some hypotheses. However, they all need an estimation of the covariance matrix R and/or the number of scatterers D. Finding D is still a hard task as the true signal distribution may be hard to model. Over uniform areas, R can be eciently estimated locally and the presented spectral analysis techniques may be used to obtain the reectivity distribution along the elevation. When the scene is heterogeneous, more complex approaches are needed to estimate R (see Chapter 7). When considering large stacks of SAR images and when the scene is very heterogeneous, the estimation of R becomes very challenging. The CS approach introduced more recently (Zhu and Bamler, 2010a) [START_REF] Budillon | Three-Dimensional SAR Focusing From Multipass Signals Using Compressive Sampling[END_REF] than most of the previous estimators is an ecient way to overcome this diculty as it uses directly the back-projection of v as stated in (4.1).

Projection and prior

The use of CS for SAR tomography is relatively recent. The estimation of the reectivity prole u along the height direction, for a given SAR resolution cell, is obtained by solving the following optimization problem:

min u ||u|| 0 s.t. v = A( h)u (4.33)
Whereas many spectral estimators are designed to retrieve a discrete set of signal parameters, the CS approach tries to retrieve the sparsest reectivity prole depending on the elevation sampling h ∈ R N h . The CS theory insures the existence of an exact solution for the problem (4.33) if the matrix A( h) satises some conditions. Chronologically, the rst one is the Restrictive Isometry Property (RIP) condition:

Denition 4.1 A matrix M satises the RIP property of order k is there exists a δ k ∈ [0, 1] such as

(1 -δ k )||s|| 2 2 ≤ ||M s|| 2 2 ≤ (1 + δ k )||s|| 2 2
where s is any vector of sparsity at most k.

This property can be understood as any subset of at most k columns of M must be as close to orthogonality as possible. A matrix satisfying the RIP of order 2k can then be seen as preserving approximately the distance between k-sparse vectors. From the RIP order of a matrix, one can compute a bound on the signal sparsity to ensure an exact reconstruction. In practice, it is however very hard to verify that a given matrix submatrices. The coherence is a more intuitive and more easily computable matrix descriptor dened by:

Let M = c 1 , • • • , c m ∈ C n×m , (4.34) µ(M ) = max 1≤i<j≤m |c H i c j | ||c i || 2 ||c j || 2
The coherence basically indicates how much the columns of M are correlated. The more they are the harder it becomes to retrieve the exact sparse signal. Again, the coherence is an indication on how orthogonal are the columns of a matrix. In the presence of additive noise, bound errors in the reconstruction with k-sparse vector can be computed from the coherence [START_REF] Ben-Haim | Coherence-Based Performance Guarantees for Estimating a Sparse Vector Under Random Noise[END_REF]. Other conditions that may be more adapted to sensing matrix consisting of an oversampled DFT and thus to SAR tomography may be found.

Going back to the SAR tomography, when a high oversampling is applied along the elevation axis, the columns of A( h) are almost fully correlated with their neighboring ones. This means that under too much oversampling, it is almost impossible to know exactly from which direction a signal is received which corresponds to a nite maximum resolution. Avoiding oversampling is also not an option as the distribution of the signal may be continuous in space.

The solution to the combinatorial problem (4.33) can be approximated using the classical convex relaxation of the 0 pseudo-norm into an 1 norm:

min u ||A( h)u -v|| 2 + µ 1 ||u|| 1 (4.35)
Dierent algorithms have been proposed to solve this problem such as Least Absolute Shrinkage and Selection Operator (LASSO) [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF], basis pursuit [START_REF] Chen | Atomic decomposition by basis pursuit[END_REF] or greedy algorithms as matching pursuit [START_REF] Mallat | Matching pursuits with time-frequency dictionaries[END_REF].

Drawbacks

The matrix A( h) is over-complete and does not guarantee to satisfy either the RIP or low coherence. The sparse reconstruction obtained through the resolution of (4.35) has nevertheless led to successful reconstructions of sparse urban scenes. However, artifacts can generally be found in those results. For instance, small spurious impulses far from the true localization of the objects or spreading of the scatterers to adjacent lines due to the oversampling. The parameter µ 1 is also generally hard to tune globally on the image because of the high dynamic of SAR images: a large value of µ 1 leads to the suppression of low intensity structures whereas a low value of it cannot allow outliers suppression.

Scatterers selection

To rene the results, the volume reconstructed with CS is generally post-processed by estimating the number of scatterers in order to select only the most signicant points.

The result is then a set of discrete points dened by their 3-D localization and complex reectivity. As the distribution of the reectivity in dense urban congurations is hard to model, this approach often fails to select points with low reectivity.

The Scale-down by 1 norm Minimization, Model selection, and Estimation Recon- struction (SL1MMER) algorithm (Zhu and Bamler, 2012a) estimates the number of non-zero points D in the cell using a Model Order Selection (MOS) technique. The MOS approaches essentially consist in nding the number of scatterers that minimize the penalized log-likelihood of the data:

D = argmin D -log p(v| ĥ(D)) + C(D) (4.36)
Under conditional model hypothesis, the last equation becomes:

D = argmin D ||v -A( ĥ(D))|| 2 2 2σ 2 n + C(D) (4.37)
where ĥ(D) is composed of the elevations of the D most powerful pixels in the cell. Dierent penalties can be used for C(D) such as the Bayesian Information Criterion (BIC) or Akaike criterion [START_REF] Burnham | Multimodel inference: Understanding aic and bic in model selection[END_REF]) [START_REF] Stoica | Model-order selection: a review of information criterion rules[END_REF].

Close to CS , the recent FAST-SUP-GLRT detector (Budillon et al., 2017a) avoids post processing selection by applying a sub-optimal statistical test taking into account the distribution of the data based on an approximated 0 norm minimization. Even if it does not take into account the geometry of the scene, the statistic of the data can be more accurately represented than with conventional CS.

As the distribution of the reectivity in dense urban congurations is hard to model, this approach may fail to select points with low reectivity.

Summary

Non-parametric estimators like conventional beamforming or Capon beamforming are easy to implement, fast to compute and adapted to continuous reectivity distributions. the assumption that the signal distribution is correctly modeled. M-RELAX can be used to maximize sequentially the likelihood but does not insure global convergence. SPICE is a recent non-parametric algorithm optimizing the reectivity of the scatterers to t the estimated covariance matrix to its model given by (3.9). Shared priors, input parameters and output representation between estimators are illustrated in the diagram also be sparse as most of the objects are punctual in ar least one direction: a wall or a roof are seen by the sensor as a surface located inside a cube. To use these priors at least two strategies may be considered: exploit the redundancies to have the best representation of the data or inverse the data and regularize the estimation.

The rst concept consists in cleverly use similar data to estimate their parameters.

This idea is at the origin of the patch based restoration methods where the neighborhood of the pixels is used to characterize them. Similar pixels, i.e. described by akin patches, are regrouped to estimate their parameters. For SAR tomography the desired parameter is generally the covariance matrix as the spectral estimators are function of it. The prior are then used in the date space to dene which and how patches are alike.

The second concept is closer to CS as no local or non-local average is used and the SLC vectors are directly inverted. The priors are then used in the estimation space to favor the reconstruction satisfying the desired behavior. With conventional CS , sparsity is the only structural prior used whereas in the following work, geometrical constraints are also proposed.

After presenting the data and the evaluation protocol, the rst strategy is addressed through chapter 7. The second approach is detailed in chapter 8 and chapter 9 for two variation of geometrical priors.

Presentation of the datasets 6.2.1 Simulations

Two stacks of images are simulated in two dierent congurations of trajectories: a regular sampling along the elevation axis h (referred to in the following as Reg ), and the actual trajectories of TerraSAR-X given in Fig. 6.5 (referred to in the following as TSX ). As in (Zhu and Bamler, 2010a), the decorrelation eect is introduced by adding a Gaussian white noise and the random scatterer phase is chosen uniform between -π and π. The reectivity distribution is then constituted of discrete point-like scatterers following the conditional signal model (3.7). The reectivity corresponding to a scatterer located at the position (x, y, z) has the following expression:

u(x, y, z) = √ τ exp(jϕ) (6.1)
with ϕ ∼ U(-π, π) strong additive noise and the geometrical decorrelation the average coherence is 0.63 for the Reg images and 0.68 for the TSX ones. Finally as the resolution is not the same in the two experiments the TSX structure is taller (20m in the Reg case against 30m in the TSX one) and the resulting images have a larger size in range.

Real data

The dataset is composed of 40 TerraSAR-X images of the south west of Paris. The corresponding temporal mean intensity image and the optical view of the scene are presented in Fig. 6.3 and 6.4. The images were acquired using the spotlight mode of the sensor in which the varying squint angle allows to improve the size of the synthetic azimuth antenna. The wavelength λ is of 3.11 cm. The images have a size 2048 × 2048 pixels with a resolution of 45 cm in range and 87 cm in azimuth. The spatial and temporal baselines are presented in Fig. 6.5. The total spatial baseline ∆b span is more than 775 m and the total temporal baseline more than 5 years with a large gap of almost two years. The theoretical resolutions in h and z are given by :

δ h = λR 2∆b = 1.75m (6.2)
δ z = δ h sin(θ) = 6.99m

Accuracy and Completeness

In (D 'Hondt et al., 2018), the authors present two metrics to compare the dierent SAR tomographic results: the accuracy and the completeness. These two errors give a complementary evaluation of an estimated a point cloud.

Accuracy For a given discrete reconstruction P, the accuracy represents the mean distance from each point in P to the ground truth P.

A( P, P) = 1 N p N p j=1 min k dist (p j -p k ) (6.4)
where pj ∈ P is the j th point of the estimated point cloud P and p k ∈ P is the k th point of the ground truth. N p is the number of points in the estimated reconstruction. The function dist is the distance used to compute the evaluation. Accuracy indicates whether reconstructed points are correctly located.

Completeness The completeness corresponds to the mean distance from each point of the ground truth to the points in P: In the presented work we chose to use the Euclidean distance for the evaluation.

The two metrics are illustrated in Fig. 6.6 where green dots represent estimated points pj and blue squares the ground truth points p k . Accuracy and completeness provide URBAN ENVIRONMENT complementary information: accuracy improves when reconstructed points are close to actual points from the ground truth, but does not indicate when points are missing (holes in the reconstruction). For instance, if we could retrieve a single point of the scene (the strongest permanent scatterer, for example) with a location very close to the ground truth, the accuracy would be excellent (A( P, P) near zero) while the completeness C( P, P) would be large, indicating that much of the scene is missing in the reconstruction. Conversely, if we loosely select points (many points for each resolution cell), we would obtain a dense volume, thus a good completeness (low completeness value C( P, P)), but erroneously selected points lying far from the true surfaces would lead to a poor accuracy score (large accuracy value A( P, P)). Pû,t = p j = (x j , y j , z j ) ∈ R 3 , such that (6.6) |û x j ,y j ,z j | ≥ |û x j ,νy,νz | and |û x j ,y j ,z j | > t (6.7)

Evaluation protocol

where ν y and ν z are the coordinates of the direct neighboring voxels restricted to the same radar resolution cell and t is a threshold.

Evaluating tomographic estimators can be done for dierent regimes (accurate reconstruction of the strongest scatterers, dense reconstruction of most scatterers). To capture these dierent cases, rather than considering a single reconstruction obtained with a giving threshold, the accuracy score A( Pû,t , P) is represented as a function of the completeness C( Pû,t , P). Applications that focus on the reconstruction of permanent scatterers will favor algorithms that achieve the highest accuracy values (even if the completeness is poor). If the reconstruction of surfaces is the aim, completeness should be favored, even if this degrades the accuracy. The evaluation protocol is summarized in the diagram shown in Fig. 6.7.

Chapter 7

Covariance Matrix Estimation

The spectral estimators used in SAR tomography are almost all based on the analysis of the covariance matrix. When the observed area is very heterogeneous with a high variety of scatterers height, material and orientation the estimation of the local covariance matrix may be very challenging. Dierent strategies have been proposed to perform this estimation. Local estimation techniques exploit only the information close to the pixel of interest whereas non-local ones are designed to exploit self similarity in the image.

Local Approaches

The estimation of the covariance matrix in SAR tomography is often done by local averaging using an isotropic kernel such as a boxcar or a gausscar lter. The value of the covariance matrix is then given by: R = I coecient of variation of the current window. When no structures appear inside the window, γ 2 I should tend toward γ 2 S and k Lee is then equal to 0. The intensity estimation is then given by the local averaging. When bright structures are present, the local variance is expected to increase resulting in a coecient tending to 1. When k is near 1, the subimage is too heterogeneous and intensity is left unchanged.

To improve the denoising, the smoothing can be combined with oriented windows allowing to preserve straight edges more eciently. The windows are composed of a rectangular mask and 8 edge-aligned oriented ones. Although continuous smooth edges can be well preserved with this method, sudden discontinuities and orientation changes produce artifacts. Moreover the limited number of orientations in the mask limits the denoising power for complex structures and textures.

Non-Local SAR algorithm

The NL-SAR algorithm [START_REF] Deledalle | NL-SAR: A Unied Nonlocal Framework for Resolution-Preserving (Pol)(In)SAR Denoising[END_REF] became in the recent years one of the top denoising algorithm for SAR images. It can be used to restore multichannel data and has shown very good results for SLC images as well as interferometric and polarimetric data. Due to its performances, this approach may be seen as very promising for SAR tomography. In the following the algorithm is summed up and results for increasing number of channels are studied.

The estimation of the pixel value can be extended to non-connected neighborhoods.

The pixels used to perform the estimation are then selected based on their similarity.

This allows to use far apart information when the local neighborhood is not sucient to perform a satisfactory denoising. This led to Lee's sigma lter [START_REF] Lee | Digital image smoothing and the sigma lter[END_REF] where the restoration of the image is done based on a similarity between the pixels.

Rather than using only a pixel wise similarity, using the similarity between the local neighborhoods (i.e. patches) of the pixels allows to use the structural information to select the pixels. This idea has been popular in various image processing algorithms that followed the seminal NL-means approach [START_REF] Buades | A review of image denoising algorithms, with a new one[END_REF]. The NL-SAR algorithm is an extension of this method for the SAR images restoration.

The similarity used in NL-SAR is derived from the fully unconditional model where the N -dimensional complex pixels follow a circular complex Gaussian distribution:

p(v|R) = 1 π N |R| exp -v H R -1 v (7.4)
where R ∈ C N ×N = E{vv H } is the complex covariance matrix, and |R| its determinant.

The empirical covariance matrix corresponds to the averaging of the L elements inside a local window:

R = 1 L L l=1 v l v H l (7.5)
L is referred as the equivalent number of looks. When N ≤ L, the empirical covariance matrix follows a complex Wishart distribution given by:

p ( R | R) = L LN | R| L-D Γ N (L)|R| L exp -Ltr(R -1 R) (7.6)
When L < N , the complex empirical covariance matrix is singular. The matrix is then said to have a degenerate distribution.

Similarity Criterion Several methods have been proposed to express a similarity criterion between pixels corrupted by a speckle eect. The one used in the NL-SAR algorithm is based on a hypothesis test on the empirical covariance matrix distribution.

For two empirical covariance matrices R1 and R2 , the test is derived as follows:

H 0 : R 1 = R 2 R 12 (7.7a) H 1 : R 1 = R 2 (7.7b)
With independent Wishart distributed empirical covariance matrices R1 and R2 , the Generalized Likelihood Ratio (GLR) corresponding to the hypothesis test in (7.7a)-(7.7b) is given by [START_REF] Conradsen | A test statistic in the complex wishart distribution and its application to change detection in polarimetric SAR data[END_REF]:

L G ( R1 , R2 ) = | R1 | L | R2 | L | 1 2 R1 + R2 | 2L (7.8)
The closer the empirical covariances matrix R1 and R2 , the larger L G is. For identical matrices, it is equal to 1 whereas it tends toward 0 for very dierent ones. Note that the GLR is only dened if the number of samples used in the empirical covariance matrices computation is larger than the number of channels N . The patches are then compared by computing the similarity between all the pixels. For two patches centered at two positions i and i , the similarity ∆ is then:

∆(i, i ) = p -log L G R(i + p), R(i + p) (7.9)
As in equation ( 7 As all the patches are compared, ∆ represents how much the structures surrounding the pixels of interests are not-alike.

Weights computation A common way to dene the weights from the (dis)similarity measure is to use an exponential kernel w(i, i ) = exp [-∆(i, i )/h], with h > 0 a ltering parameter. The NL-SAR algorithm proposes a framework that makes the computation of the weights independent of the parameters such as the search window and the patch size, the number of samples used for the estimation of the empirical covariance matrices or the number of channels. To do so the exponential kernel is composed with two other functions:

F the cumulative distribution function of ∆ under the hypothesis H 0 G -1 the reciprocal of the χ 2 distribution.

The weight for a pixel i respectively to a pixel i is then dened as:

w(i, i ) = exp - G -1 • F [∆(i, i )] h (7.10)
The χ 2 function corresponds to the distribution of the dissimilarity ∆ under a Gaussian additive noise where the exponential kernel has proven to be ecient for denoising [START_REF] Buades | A review of image denoising algorithms, with a new one[END_REF]. The NL-SAR weights computation corresponds then to a mapping to the SAR image statistics to an additive Gaussian scenario. Finally to avoid having any pixel having a higher weight than the central one and favor pixels following the same distribution, the kernel is modied as follow:

w(i, i ) =      exp -G -1 •F[∆(i,i )-c] h if i = i 1 otherwise (7.11) with c = E G -1 • F [∆(i, i )] |H 0 .
The non-local estimation of the covariance matrix is then given by:

RNL (i) = i w(i, i )v i v H i i w(i, i ) (7.12)
Bias reduction To reduce the spreading of bright structures induced by the empirical covariance matrix pre-estimation step, a bias correction step is applied to the non-local estimate. This step is similar to the Lee lter strategy and the corrected Non-Local Bias Reduced estimation is then:

RNLRB (i) = RNL (i) + α v i v H i - RNL (i) (7.13)
The parameter α computation is detailed in [START_REF] Deledalle | NL-SAR: A Unied Nonlocal Framework for Resolution-Preserving (Pol)(In)SAR Denoising[END_REF] and depends on the variance of the non-local estimation.

Unsupervised local adaptation

The NL-SAR algorithm present many parameters:

the size of the patches, the size of the search window, the number of samples used for the local estimation of the empirical covariance matrices or the scale parameter h.

However, one of the strength of the method is its unsupervised adaptation to the local context. The quality of the estimation is given in terms of variance reduction and the set of parameters is set locally to have the best reduction. As the variance decreases monotonously with the number of samples (or looks) used for the estimation, the best set of parameters is the one maximizing the number of looks. After non-local estimation and bias-reduction, the equivalent number of looks is given by [START_REF] Deledalle | NL-SAR: A Unied Nonlocal Framework for Resolution-Preserving (Pol)(In)SAR Denoising[END_REF]:

LNLRB (i) = LNL (i) (1 -α) 2 + α 2 + 2α(1-α) i w(i,i ) (7.14)
The set of parameters maximizing this quantity is then used for the estimation. Examples of denoising using the presented local and non-local algorithms are presented in Fig. 7.1 and 7.2.

Proposed pixel similarity for large tomographic stacks of urban areas

In urban areas where the amplitudes of bright pixels are signicantly higher than their surrounding background and where the SAR images are very heterogeneous, nding similar patches may be a dicult task. When the dimension of the SAR stack increases patches may then appear more and more unique within the search window. From the covariance matrix structure detailed in the MUSIC algorithm, it is possible to get an intuition justifying the reduced eciency of the Wishart GLR :

R 1 = V 1            µ 1 µ d 1 σ 2 n σ 2 n            V H 1 , R 1 = V 2            ν 2 ν d 2 σ 2 n σ 2 n            V H 2 and R 12 = V 12            ρ 2 ρ d 12 σ 2 n σ 2 n            V H 12 (7.15)
where V 1 , V 2 and V 12 are composed of the eigenvectors of the corresponding matrices, d 1 , d 2 and d 12 are the number of signicant eigenvalues; {µ n } 1≤n≤d 1 , {ν n } 1≤n≤d 1 and {ρ n } 1≤n≤d 1 , the corresponding eigenvalues. σ 2 n is the additive noise power. Now, let's look at the simple example where R1 and R2 follow very similar distributions except for the last subspace. Here

d 1 = d 2 = d 12 = d and V 1 H k V 2 k =    1 if k = d 0 otherwise and µ k = ν k = ρ k if k = d (7.16)
The last assumption states that all the eigenvectors of V 1 and V 2 are the same except for the d th ones that are orthogonal. Then the GLR under Wishart distribution for the empirical covariance matrices as dened in (7.8) is:

L G ( R1 , R2 ) = | R1 | L | R2 | L | 1 2 R1 + R2 | 2L = d-1 i=1 µ L i ν L i µ i +ν i 2 2L 2 4 σ 2 n µ d ν d 2L (7.17)
As the noise power is generally at least 2 or 3 dB lower than the scatterers power and

µ i ν i µ i +ν i 2 2 ≤ 1, L G is
L G ( R1 , R2 ) = | R1 | L | R2 | L | 1 2 R1 + R2 | 2L = d-1 i=1 µ L i ν L i µ i +ν i 2 2L 2 4L 2 4 σ 2 n µ d ν d ρ d ρ d+1 2L (7.18)
The ratio

µ d ν d ρ 2 d
is also smaller than 1 through Weyl's Inequalities (see below). Again the ratio

2 4L σ 2 n ρ 2 d+1
may be very low. The more acquisitions, the more likely the back-scattering signals are to dier and the more probable it becomes to nd dissimilar subspaces.

Property 7.1 Weyl's Inequalities Let M , N and P ∈ C N with at least two of them being Hermitian and M +N = P . Let µ n 1≤n≤N , ν n 1≤n≤N and ρ n 1≤n≤N be their respective eigenvalues in decreasing order, then for all i = 1, • • • , N :

µ i + ν N ≤ ρ i ≤ µ i + ν 1 and µ N + ν i ≤ ρ i ≤ µ 1 + ν i
If the matrices are also positive denite then:

∀i = 1, • • • , N, µ i < ρ i and ν i < ρ i
It is also noticeable that the orthogonality example is not a theoretic conguration in urban areas. Indeed, let's consider two signals with only one strong contribution.

Their covariance matrices are then:

R 1 = µa(h 1 )a H 1 + σ 2 n I N R 2 = νa(h 2 )a H 1 + σ 2 n I N (7.19)
where h 1 and h 2 stand for the elevation of the scatterer of the rst and second pixel.

The signal subspaces are orthogonal if:

N i=1 a(h 1 ) H a(h 2 ) = N i=1 exp -jξ i (h 1 -h 2 ) = 0 (7.20)
When the baselines are uniformly distributed with a sampling distance ∆b, the last condition is veried if ∆h = λR 2∆bN . For the TerraSAR-X sensors, this corresponds to a dierence in elevation of 5m or equivalently 3m variation with respect to the vertical axis. For very bright scatterers with high backscattering power, this implies that they should be seen as more dissimilar than a bright point versus pure noise. This overly discriminative power may be softened by adapting the scale parameter and/or applying a stronger smoothing when estimating the empirical covariance matrices. Moreover, when the number of channels is high, empirical covariance matrices need to be estimated with a lot of samples or need to be shrinked to ensure having non singular matrices.

Then, due to the reduction of the pre-estimation quality, it is harder to dierentiate neighboring pixels. Large scale and smoothing parameters tend also to make the bias reduction parameter closer to 1 resulting in an output image closer to the input one.

Examples of denoising in a dense urban scenario and the similarity maps corresponding to the GLR under Whishart distribution are presented for increasing numbers of et al., 2017). It has also been used successfully for small sets of data of hybrid areas (D 'Hondt. andAl., 2017) (D'Hondt et al., 2018). Few strategies have been proposed to extend its good performances to high number of channels. A recently proposed approach is to perform the denoising for each pair of interferometric images and then reconstruct the covariance matrix [START_REF] Shi | Non-Local Compressive Sensing Based SAR Tomography[END_REF]. This however does not guaranty a coherent estimation of the weights between pairs or to take into account long-term variation of the data.

Proposed PS and background based similarity

For SAR tomography, the interferometric phase is the most important information to retrieve. The similarity could then be improved in focusing more on the scatterers geometric distribution inside the SAR pixels as illustrated in Fig. 7.5. In urban areas, most of the usable information is carried by the PS that stay highly coherent between the acquisitions. A new similarity function could then be based on the conditional or hybrid model to have an explicit information of the strong scatterers in the pixels distribution. This contribution was presented in [START_REF] Rambour | Similarity Criterion for SAR Tomography over Dense Urban Area[END_REF]. Under hybrid model, the pixel v follows an non-centered N dimensional circular Gaussian distribution:

p (v | R, h, u) = 1 π N |R| exp -v -A(h)u H R -1 v -A(h)u (7.21)
To determine if two pixels v 1 and v 2 are similar or not based on their PS like scatterers, we adopt an iterative approach where we jointly extract their PS components, update the mean and covariance of their distribution and test if they are likely to share the same distribution or not. Rather than testing the equality of height and reectivity of the targets the similarity is based on their likelihood as in equation ( 7.8). The extraction of the main components of the signal is done using the M-RELAX [START_REF] Li | Ecient mixed-spectrum estimation with applications to target feature extraction[END_REF] scheme to solve the following non-linear least squares:

û, z = argmin (u,z) ||v -A(z)u|| 2 R -1 (7.22)
The proposed algorithm iteratively estimates the main signal components from the empirical covariance matrix of the two pixels v 1 and v 2 . The extracted targets are compared with the ones of the mixed pixel v 12 = 1 2 v 1 + v 2 to compute the similarity.

After the rst iteration of the algorithm, the rst components from R1 , R2 and R12 = 1 2 R1 + R2 are extracted, respectively (û

(1) 1 , ẑ(1) 1 ), (û (1) 2 , ẑ(1)
2 ) and (û

(1)
12 , ẑ(1) 12 ). The estimation of the empirical covariance matrices of v 1 and v 2 can then be updated as in the algorithm M-RELAX :

vi (1) = v i - û(1) i a(ẑ (1) i ) (7.23) R(1) i = 1 L L l=1 v(1) i l v(1) H i l for i = {1, 2}. (7.24)
Then, extracted scatterers similarity is derived from the following hypotheses:

   H (1) 0 : u (1) 1 = u (1) 2 = u
(1) 12 and z

(1) (7.25) the observed buildings are lower and closer to each others. In this example, the Wishart GLR manages to select points corresponding to buildings but the similarity image does not show a lot of contrast. The proposed similarity presents a more structured selection with a clear preference for walls. The selected points are roughly concentrated at the same height as the central pixel. Finally, the discriminative power of the proposed similarity is robust to an increasing number of images. This point is illustrated in Fig. 7.8 on the Keller Tower and 7.9 on dense smaller buildings for tomographic stacks of size 10, 20, 30 and 40.

1 = z (1) 2 = z (1) 12
H (1) 1 : u (1) 1 = u (1) 2 or z (1) 1 = z (1) 2

Limits of Non Local approaches for urban tomography

Although we have shown that the similarity used for non-local selection can be improved, this approach still presents drawbacks for SAR tomography over urban areas. Firstly, the proposed similarity function is time consuming and can then hardly be consider for a non local approach. Secondly, the current debiasing step can lead to singular matrix specially for big dimensions which may completely distort the tomographic estimations.

The max operator used in the debiasing step of NL-SAR could thus be changed to a less radical one for big dimensions. Moreover, even for dense urban areas, no clear benet can be notice using a kernel based on the similarity as illustrated in the Fig. Chapter 8

Tomographic inversion with spatial regularization

The presented work in this chapter comes mostly from (Rambour et al., 2018a). Some detailed were added in the optimization section.

In order to enforce some spatial smoothness, the tomographic inversion has to be performed globally in ground coordinates. Rather than considering the collection of measurements at a given radar pixel, from now on the notations u ∈ C Nx.Ny.Nz and v ∈ C Nx.Nr.N will respectively refer to column vectors obtained by stacking all the values in the 3-D volume, and all the values in the tomographic stack of SAR images.

N x , N y and N z are the number of voxels in each direction in ground geometry while N r is the size of SAR images in the range direction and N is the number of images in the tomographic stack. The linear operator Φ ∈ C (Nx.Nr.N )×(Nx.Ny.Nz) maps the volume of complex reectivities in 3-D space to the complex amplitudes in the tomographic stack of SAR images, see Fig 8 .1. An element of Φ is dened as:

Φ i,j =          exp(-jϕ) if x i = x j and r i -δr 2 < ρ 1;y j ,z j < r i + δr 2 , 0 otherwise. (8.1)
with ϕ = ξ i h(y j , z j ) the phase shift due to the path between voxel j with coordinates (x j , y j , z j ) and antenna i. As in the previous equations, ρ 1;y j ,z j corresponds to the distance between antenna 1 (of the master image) and the point with ground coordinates (y j , z j ). The size of a radar pixel in the range direction is noted δ r .

The construction of the matrix Φ is illustrated on Fig. 8.2 and Φ is sparse: only a few entries are dierent from zero, so that products of the form Φu can be computed eciently.

The observed SAR tomographic stack v can be modeled by the following (complex- valued) linear model: where stands for the noise. This corresponds to a generalization of the tomographic direct model (4.1) where we additionally consider the geometric transformation from ground geometry to SAR geometry and model at once the measurements for all the pixels. In order to invert this tomographic model, it is necessary to introduce some regularization terms.

v = Φ u + (8.2)

Ground base spatial priors for SAR tomography

It is often desired to reconstruct volumes with a discretization in heights that is ner than the resolution given by the synthetic aperture in the height direction (i.e., superresolution). The inversion of equation (8.2) is therefore ill-posed (more unknowns than measurements) and requires some regularization. In the following, we denote R the regularization function. Since the intrinsic phase of a scatterer is typically modeled as uniformly distributed and independent from one scatterer to another, no specic regularization can be enforced on the phase of our unknown complex reectivities u. We dene the regularization R as a function of the modulus of u only. The reconstruction of the volume of complex reectivities û is thus obtained by solving an optimization problem of the general form:

û = arg min u 1 2 ||Φ u -v|| 2 2 + R(|u|) (8.3)
To design the regularization function R, we need to select a function that favors volumes of reectivities |u| that are often present in urban environments. Many dierent such functions could be considered, we selected a function based on the two following remarks:

The 3-D scene can be represented as a sparse volume. Indeed a good reconstruction should retrieve only the illuminated part of the buildings and of the ground.

The estimated volume is then mostly lled with zero intensity voxels.

The illuminated structures are spatially smooth (continuous surfaces: frontages, rooftops, ground).

These remarks suggest the following regularization function:

∀w ∈ R Nx.Ny.Nz , R(w) = µ x 2 D x w 2 2 + µ y 2 D y w 2 2 + µ z 2 D z w 2 2 + µ 1 w 1 (8.4)
where the matrices D x , D y and D z stand for the nite dierences operators in the x, y and z directions, and parameters µ x , µ y , µ z and µ 1 weight each term. The 1 norm favors vectors w with many zeros while the terms with the nite dierence operators enforce a spatial smoothness. As µ 1 controls the sparsity of the reconstructed volume, it plays a crucial role.

In urban environments, dihedral and trihedral structures produce very strong echoes.

When only an 1 norm is minimized, it is hard to reconstruct at the same time very strong scatterers and weaker scatterers on the ground or rooftops. The role of the rst three terms is thus to favor spatial smoothness and hence preserve these scatterers whenever they are close to other scatterers. In urban areas most of the buildings are expected to show straight walls along the z direction. Of course the true orientation of the objects is a priori unknown and it is most likely that they may not always follow the

x and y directions. However, when the horizontal smoothing is small compared to the power of the scatterers, the only major cost is to introduce some thickness in vertical objects with the benet to ensure a global reduction of isolated outliers and retrieve part of at areas with a lower intensity. When reconstructing images of urban areas, we were unable to nd a value of the sparsity parameter µ 1 that would both preserve the weakest scatterers and successfully suppress side-lobes in areas with strong scatterers.

To improve the reconstructions, we introduced a spatially variant regularization based on the square-root of the estimated intensity of the master image which can be obtained using a denoising algorithm such as NL-SAR [START_REF] Deledalle | NL-SAR: A Unied Nonlocal Framework for Resolution-Preserving (Pol)(In)SAR Denoising[END_REF] or simply using the average intensity depending of the conguration of the acquisitions. We dene the diagonal matrix D 1 whose j-th diagonal entry is equal to the square-root of the estimated intensity at the corresponding azimuth and range coordinates (i.e., such that the range r veries r -δr 2 < ρ 1;y j ,z j < r + δr 2 ). The equation (8.4) is then modied into:

∀w ∈ R Nx.Ny.Nz , R(w) = µ x 2 D x w 2 2 + µ y 2 D y w 2 2 + µ z 2 D z w 2 2 + µ 1 1 T D 1 w (8.5)
where 1 is the vector of size N x .N y .N z with each entry equal to 1.

Optimization algorithm

The minimization problem (8.3) is not easy. It is indeed large scale (millions up to several billions unknowns) and non-convex. The non-convexity is illustrated in a simple case where the vector u has only two elements in Fig. 8.3: the spatial smoothness favors vectors such that the modulus of each entry is close. Since the regularization is independent on the phase, the minimum (for a xed value of u 1 ) corresponds to a set of complex values with equal modulus (the white circle drawn on Fig. 8.3).

We have to use an algorithm that can deal with the non-convexity of the problem.

To solve the minimization problem (8.3), we use a variable splitting approach in order to break down the problem into a sequence of simpler problems. We introduce two new vectors: f (complex-valued: f ∈ C (Nx.Ny.Nz) ) and w (real-valued and non-negative: w ∈ R +(Nx.Ny.Nz) ). Problem (8.3) is formally equivalent to the following constrained problem:

û = arg min u 1 2 Φ u -v 2 2 + R(w) (8.6) s.t.    u = f |f | = w
Dierent methods could be consider to nd a solution of (8.7). The Lagrangian dual ascent for instance would be one them if it was for the non-convexity introduced by the modulus. To insure convergence we then look at algorithms exploiting the augmented Lagrangian of the previous problem:

L(u, f , w, d 1 , d 2 ) = 1 2 ||Φu -v|| 2 2 + β 1 2 ||f -u + d 1 || 2 2 + β 2 2 w -|f | + d 2 2 2 + R(w) (8.7)
where d 1 ∈ C (Nx.Ny.Nz) and d 2 ∈ R (Nx.Ny.Nz) are the scaled dual variables and β 1 and β 2 are penalty parameters (relevant only to the optimization method, i.e., impacting the convergence). Alternating Direction Method of Multipliers (ADMM) is a well known augmented Lagrangian method. It backs from works in the 1970s [START_REF] Glowinski | Sur l'approximation, par éléments nis d'ordre un, et la résolution, par pénalisation-dualité d'une classe de problèmes de Dirichlet non linéaires[END_REF])(gab, 1976) and is very close to other dual decomposition algorithms such as the Lagrangian Methods of Multipliers (LMM) [START_REF] Bertsekas | Constrained optimization and Lagrange multiplier methods[END_REF], the Douglas-Rachford splittings [START_REF] Douglas | On the numerical solution of heat conduction problems in two and three space variables[END_REF]), Dykstra's alternating projections [START_REF] Boyle | A method for nding projections onto the intersection of convex sets in hilbert spaces[END_REF] or Bregman iterative algorithm for 1 minimization prob- lems [START_REF] Yin | Bregman Iterative Algorithms for _1-Minimization with Applications to Compressed Sensing[END_REF]. ADMM is based on the LMM algorithm which for the current problem consists in the following iterations:

u (k+1) , f (k+1) , w (k+1) = argmin u,f ,w L(u, f , w, d (k) 1 , d (k) 2 )
(8.8)

d (k+1) 1 = d (k) 1 + β 1 2 f (k+1) -u (k+1) (8.9) d (k+1) 2 = d (k) 2 + β 2 2 w (k+1) -|f (k+1) | (8.10)
Minimizing the augmented Lagrangian jointly with respect to u, f , and w may be hard to achieve. The ADMM solves the last problem by taking advantages of the dual decomposition and optimizes sequentially the primal variables. For the given problem, REGULARIZATION the algorithm is described by iterating the following steps:

u (k+1) = argmin u 1 2 ||Φ u -v (k) || 2 2 + β 1 2 ||f (k) -u + d (k) 1 || 2 2 (8.11) d (k+1) 1 = d (k) 1 + β 1 2 f (k) -u (k+1) (8.12) f (k+1) = argmin f β 1 2 ||f -u (k+1) + d (k+1) 1 || 2 2 + β 2 2 ||w (k) -|f | + d (k) 2 || 2 2 (8.13) d (k+1) 2 = d (k) 2 + β 2 2 w (k) -|f (k+1) | (8.14) w (k+1) = argmin w β 2 2 ||w -|f (k+1) | + d (k+1) 2 || 2 2 + R(w) (8.15) (8.16)
Dierent convergence proofs of ADMM exist in the litterature. From [START_REF] Boyd | Distributed optimization and statistical learning via the alternating direction method of multipliers[END_REF], the convergence of LMM and ADMM is guaranteed under the following two hypothesis:

The objective functions are closed, proper and convex

The unaugmented Lagrangian L 0 has a saddle point:

L 0 (u, f , w, γ 1 , γ 2 ) = 1 2 ||Φu -v|| 2 2 + R(w) + γ H 1 u -f + γ H 2 |f | -w (8.17)
with γ 1 and γ 2 the unscale dual variable.

These conditions are very general and allow to use these algorithms for many applications. An other strength of these approaches is that they converge even if the minimization steps are inexact. In our case, the convexity assumption is not satised so LMM and ADMM can only insure local convergence and the obtained solution may depend on the initialization.

To solve this constrained optimization problem, we apply a variation of the precedent methods described in [START_REF] Mourya | Augmented Lagrangian without alternating directions: Practical algorithms for inverse problems in imaging[END_REF] and minimize jointly on the variables u and w while f is substituted with its optimal value f * (u, w):

f * = arg min f β 1 2 ||f -u + d 1 || 2 2 + β 2 2 w -|f | + d 2 2 2 (8.18) = β 1 • |u -d 1 | + β 2 • (w + d 2 ) β 1 + β 2 + exp j • arg u -d 1 (8.19)
where . + is the identity on [0, +∞[ and is the constant null function on ] -∞, 0[.

Proof:

The second term in (8.18) only depends on the modulus of f . The phase of f * is then driven by the rst term and must be chosen equal to that of ud 1 so as to minimize the cost function. There remains to estimate the modulus of f * which is a solution of a 1-D quadratic problem arg min

ρ≥0 β 1 2 ||ρ -|u + d 1 ||| 2 2 + β 2 2 w -ρ + d 2 2 2 (8.20)
ρ is either given by the unconstrained solution:

ρ * = (β 1 • |u -d 1 | + β 2 • (w + d 2 ))/(β 1 + β 2 ),
or ρ = 0. The optimal value of ρ and of the phase together lead to the expression (8.19).

This leads to the following algorithm which alternates between the joint minimization with respect to variables u and w, and the update of dual variables.

Algorithm Tomographic SAR 3-D Inversion

Input: v (stack of SLC SAR images)

Output: û

(3-D cube of complex reectivities)

Initialization :

1: d1 ← 0 2: d2 ← 0 3: while not converged do 4:

{û, ŵ} ← approximate_min(v, d1 , d2 , û, ŵ) 5:
d2 ← d2 + ŵ -|f * (û, ŵ)| 6: d1 ← d1 + f * (û, ŵ) - û
7: end while 8: return û Compared to ADMM , this hierarchical approach has been shown in [START_REF] Mourya | Augmented Lagrangian without alternating directions: Practical algorithms for inverse problems in imaging[END_REF] to converge faster and to be less sensitive to the tuning of the optimization parameters β 1 and β 2 ; moreover, we improve the convergence by constraining w to be positive, such a constraint would be costly to enforce with ADMM .

After replacing f by its optimal value f * (u, w) and constraining w to be positive, the objective function is dierentiable (since w ≥ 0, w 1 = 1 T w, which is dierentiable; moreover, after substituting f by the optimal value f * in the quadratic terms, the cost function becomes smooth, see [START_REF] Mourya | Augmented Lagrangian without alternating directions: Practical algorithms for inverse problems in imaging[END_REF]). The cost function can thus be minimized using a limited memory quasi-Newton algorithm that handles positivity constraints, such as L-BFGS-B [START_REF] Zhu | Algorithm 778: L-BFGS-B: Fortran subroutines for large-scale bound-constrained optimization[END_REF], or the slightly more ecient algorithm VMLM-B [START_REF] Thiébaut | Optimization issues in blind deconvolution algorithms[END_REF] that we used, with Eric Thiébaut's freely available implementation 1 . The minimization step described in Procedure approximate_min does not need to be performed up to a high precision. A few (e.g., ten) iterations of the quasi-Newton algorithm are sucient since the algorithm is warm-restarted. In the denition of the gradient of the cost function (lines 2 to 4 of Procedure approximate_min), we used the fact that

∂C ∂[f * ] i ∂[f * ] i ∂[u] j = ∂C ∂[f * ] i ∂[f * ] i ∂[w] k =
0 for all i, j, k, even though f * depends on variables u and w, by applying the chain rule and noting that ∂C ∂[f * ] i = 0 for all i. Output: {û, ŵ}

(approximate solution) 1: C (v, d1 , d2 , û, ŵ) = 1 2 ||Φu -v|| 2 2 + R(w) + β 1 2 ||f * (u, w) -u + d1 || 2 2 + β 2 2 w - |f * (u, w)| + d2 2 2 , 2: ∇ u C = Φ H (Φu -v) + β 1 (u -f * (u, w) -d1 ) 3: ∇ w C = (µ x D T x D x + µ y D T y D y + µ z D T z D z )w + µ 1 1 +β 2 (w -f * (u, w) + d2 ) . 4: ∇C (v, d1 , d2 , û, ŵ) = ∇ u C ∇ w C ,
Call quasi-Newton minimization algorithm :

5: {û, ŵ} ← quasi_Newton C (v, d1 , d2 , û, ŵ), ∇C (v, d1 , d2 , û, ŵ), constraint: ŵ ≥ 0
The necessary condition for f * being a minimizer of the augmented Lagrangian is then respected. As the function we want to minimize is non-convex, only the convergence to a local minimum can be expected. In practice the algorithm can then be stopped after a xed number of iterations (around 60) or when the distance between two iterations falls below a given threshold.

The priors proposed in this section are designed to be simple and to lead to a function R that is easy to optimize. In the following experiments, the relevance of this method is shown by comparing to other state of the art estimators on simulated and real data. The proposed framework is very general and can easily be adapted to include other spatial regularizations expressed in ground coordinates.

Inuence of the dierent parameters

To study the inuence of each regularization parameter, we performed dierent reconstructions of the Reg scene for various combinations of regularization values (in particular, with one parameter chosen so as to illustrate the eect of either under-regularization or over-regularization). In these expreriments, the phase noise is not present to show clearly the structural inuence of the spatial smoothing. The results are presented in gures 8. 4, 8.5, 8.6 and 8.7. On this well-sampled scene with a good spatial resolution, the parameter with the largest inuence is µ 1 , associated to the sparsity constraint. Fig. 8.4 illustrates that when µ 1 is too large, there are some holes in the reconstruc-tion, while a value of µ 1 that is too small leaves side-lobes and outliers. The eect of over-smoothing (values of µ x , µ y or µ z too large) is to expand the structures in the direction of the spatial smoothing. This is visible in particular with structures whose orientation diers from that of the smoothing, see in Fig. 8.6 the widening of the wall due to excessive smoothing in the horizontal direction. The same eect appears on the rooftop in Fig. 8.7. Insucient smoothing translates into residual uctuations (i.e., large variance) that are reduced by increasing spatial smoothing, see in particular the corner between the roof and the wall, or the ground and the wall in Fig. 8.6 and 8.7.

Outliers located far from the actual surfaces can also be observed when the spatial regularization is too weak.

By combining sparsity and spatial smoothness constraints, our algorithm is very exible and applicable to the reconstruction of dierent kinds of areas. The downside of this exibility is the necessity to tune four regularization parameters. In numerical simulations, the ground-truth can be used to select the set of regularization parameters {µ x , µ y , µ z , µ 1 } that oers the best performance, as measured by the Minimum Accuracy / Completeness Trade-o (MACT) (MACT) min A( Pû,t * , P) 2 + C( Pû,t * , P) 2 .

The regularization parameters can then be tuned in order to reach the best possible trade-o. We recommend performing this tuning by order of importance: µ 1 , then µ z , then µ x and µ y . While tuning the sparsity parameter, the other smoothing parameters should be put to 0.

Derivative-free methods can also be used to set all the parameters at once. We compared the described alternating minimization strategy with Nelder-Mead's simplex method by initializing with dierent set of parameters and found similar values for the best accuracy/completeness trade-o. The dierences for µ 1 is of 21% and below 14%

for the smoothing parameters. In the absence of ground truth, a simple numerical simulation using the same geometrical conguration and SNR can be generated in order to automatically tune the parameters. If a simulated scene is not an option, the alternating minimization strategy can still be applied in an interactive fashion: the user tunes each parameter, by order of importance, in order to reach a satisfying reconstruction.

To further illustrate the behavior of our algorithm with respect to its regularization parameters, we plotted accuracy as a function of completeness for dierent sets of parameters. These curves are drawn for our two simulation cases Reg and TSX respectively in Fig. 8.8 and 8.9. We observe that the best accuracy / completeness trade-o (point of the curve closest to the origin of the axes) is reached for a unique set of parameters that can thus be found for example by binary search.

The optimization parameters β 1 and β 2 have an impact on the convergence speed.

We found that, when starting from a volume initialized at zero, using large penalty parameters β 1 and β 2 produces very quickly a sparse reconstruction while lower contrasted structures are correctly reconstructed after many more iterations. Smaller values of the parameters help to reconstruct those structures, at the cost of a slower convergence (i.e., sidelobes suppression) in the brightest areas. Penalty parameters β 1 and β 2 can be set according to methods described in [START_REF] Boyd | Distributed optimization and statistical learning via the alternating direction method of multipliers[END_REF]; In the easy case where the antennas are uniformly distributed, classical beamforming presents good performances with respect to the accuracy-completeness metrics, achieving a smaller error than Capon beamforming or MUSIC. However this method is limited by its available resolution introducing thickness in the wall. In the TSX case where the theoretical resolution is worse, very little of the wall is the retrieved, leading to a poor completeness score. A closer inspection of the TSX reconstructions indicates the presence of side-lobes (i.e., many outliers).

Both Capon beamforming and MUSIC, which rely on the covariance matrix to achieve super resolution, manage to suppress most of the outliers. However, the averaging introduced for the covariance matrix estimation produces an extension of the ground and walls. For the Reg case, this explains the bounded accuracy scores even when the completeness is poor. Capon beamforming manages to reduce the side lobes in the Reg case but not in the TSX simulation where the irregular sampling produces dramatic biasing eects on the estimation of the position of the scatterers.

In the TSX scene, the best accuracy / completeness trade-o from the MUSIC pseudo-spectra corresponds to a sparse scene with no outliers. However the bias present in the estimation of the covariance matrix strongly deforms the wall and introduces some fake discontinuities.

All the spectral estimators presented (classical beamforming, Capon beamforming and MUSIC) suer from the averaging step needed to estimate the covariance matrix, which produces noticeably distorted walls and/or rooftop in the TSX reconstruction, and an expansion of the rooftop and a corresponding reduction of the ground in the Reg case.

As expected, CS is one of the top estimators in accuracy and often manages to achieve the most accurate representations. However, the precision comes with the cost of a loss of information which explains the relatively bad scores in completeness.

The proposed 3-D inversion algorithm shows improved performances compared to the CS approach. Some outliers that were present even when imposing an 1 penalty are suppressed thanks to the spatial smoothness prior.

The tuning of the parameter in the simulated experiences has been done by exploiting the knowledge of the scene. Of course, when dealing with real data, this information is not available. Then, to test if the proposed method is able to exhibit better performances when the tuning of the parameter is not ideal, the model parameter are modied with increasing perturbation values. The results are shown in the table 8.1.

Up to an isotropic perturbation of 50, the proposed method still shows a better Minimum Accuracy-Completeness Tradeo (MACT) than the other tested methods. An other concern is whether the parameters need to be adapted to the scatterers power.

Real

Real Data

We now compare our algorithm to other SAR tomographic estimators on a stack of 40 TerraSAR-X images acquired in spotlight mode over the front de Seine in the south-west of Paris, France. The slant-range resolution is 0.45 m and the azimuth resolution 0.87 m. The observed scene is presented in Fig. 6.3 and 6.4 in 6. The spatial and temporal baselines are shown in Fig. 6.5. The total spatial baseline span ∆b is more than 775 m and the total temporal baseline more than 5 years with a large gap of almost two years.

The theoretical resolutions in h and z are given by :

δ z = δ h sin(θ) =
λr sin(θ) 2∆b = 6.99m

(8.21)

with an incidence angle θ = 0.6 rad and wavelength λ = 0.0311 m. Sparse reconstructions in urban SAR tomography have been shown to signicantly improve this resolution (Zhu and Bamler, 2010b).

To evaluate the dierent tomographic estimators we use a rough ground truth of the scene. The ground truth surface is shown in Fig. 8.15 2 . The big structures such as the skyscrapers or the tall buildings are well represented by smooth polygons. However, the vegetation and the small structures on the docks, bridges or streets are not represented.

We therefore performed the evaluation only on some areas where we found the ground truth to be accurate enough. The results of the evaluation are given in Fig. 8.15.

It is noteworthy that MUSIC performs particularly well on areas containing a lot of ground surface. In the rst test, the sub-scene is composed of a mixture of tall buildings, small structures near the ground level and smooth at areas (streets and docks). As the diuse signal back-scattered by at surfaces is much weaker than the one reected by the dihedral or trihedral structures, it is generally not taken into account by sparse representations of the scene. CS is then unable to correctly represent ground areas.

Taking more points into account results only in more outliers generally due to the sidelobes of the brightest points. The averaging step used to estimate the covariance matrix allows MUSIC to retrieve part of the ground points or to extend the signal coming from one punctual target close to the ground to neighboring pixels. Unlike the non-parametric beamforming methods, MUSIC is designed to retrieve a sparse scene which removes most of the side lobes and outliers. MUSIC is then able to outperform CS according to our evaluation method that includes scatterers on the ground. On the second tested area corresponding to the red rectangle in Fig. 8.15, the performances of MUSIC and CS in terms of accuracy are very similar. MUSIC seems to perform a little bit better than CS. However, this may be due to its ability to retrieve more points on the ground. The small dierence between the two methods and the limited accuracy of the ground truth makes it dicult to draw a clear conclusion between the two methods.

In both experiments, the presented algorithm achieves the best scores in term of accuracy and completeness. Moreover, the analysis of the metrics on real data shows 2 Ground truth on Paris complimentary provided by the IGN, France. Chapter 9

Segmentation in urban SAR tomography

The presented work in this chapter comes mostly from (Rambour et al., 2018b). Some results were added in the experimental part.

Graph-cut based surface segmentation

Starting from the tomographic reconstruction (a 3-D volume u) obtained with one of the methods described in the previous chapters (see chapters 4 and 8) we aim to recover the urban surfaces (ground, building facades, roofs). Following a typical approach in computer vision for surface reconstruction, we formulate the problem as an energy minimization problem [START_REF] Kolmogorov | Kolmogorov and Zabih's Graph Cuts Stereo Matching Algorithm[END_REF]. We seek a surface S corresponding to an elevation map: (x, y) → z = E (x, y) that both ts well the reconstructed tomographic volume and that is smooth. We rst formulate a cost function that captures these two properties, then we describe an ecient graph-based algorithm to perform the minimization of the cost function.

Denition of the cost function

The rst component of the cost function favors surfaces that are faithful to the reconstructed tomographic volume. We seek surfaces such that, when considering a given ray direction in 3-D space, the scatterer encountered along the ray falls close to the ray-surface intersection, see Fig. 9.1. The reectivity prole along the ray may display several local maxima due to residual sidelobes after the tomographic inversion.

Rather than detecting these maxima and deciding for the most meaningful maximum, we consider that a satisfying location of the surface is a location such that the reectivity prole is split into two well-balanced halves. We dene the cumulative reectivity

C -(r s ) from the antenna to the surface S and the cumulative reectivity C + (r s ) from C -(r equi )). Then, if r > r equi , the rst integral equals r r equi (C -(r s ) -C + (r s )) dr s . Conversely, the second integral in (9.3) is non-zero only if the distance r s is smaller than the distance of equilibrium r equi . It is then equal to

r equi r (C + (r s ) -C -(r s )) dr s .
D(r) is thus a function that monotonically increases with the distance |r -r equi | and that is minimal and equal to zero when r = r equi .

The second component of the cost function guarantees that the segmented surface be smooth. To prevent the surface from oscillating in order to pass through the position of equilibrium r equi for each ray, we penalize the area A (S ) of the surface. In order to favor surfaces with horizontal or vertical parts, we suggest measuring the area with respect to the 1 distance (i.e., Manhattan distance p

1 = |p x | + |p y | + |p z |).
To summarize, we suggest dening the segmentation as the surface S that is a solution to the following variational problem:

min S ray∈R D ray (r ray→S ) dR + β A (S ) , (9.4)
where S is required to be representable as an elevation map E (x, y) (formally, there exist a function E : (x, y) → E (x, y) such that S be the boundary of the epigraph of E ). To prevent from introducing too many notations, we denote 'ray' for the generic denition of a ray in an adequate parameterization (a line in 3-D space), R represents the set of all rays, r ray→S is the distance from the radar to the surface S along the direction dened by 'ray', D ray is the penalty dened by equation ( 9.3) for the direction specied by 'ray'. Finally, β is a parameter that balances the delity to the tomographic reconstruction and the spatial smoothness of the surface.

Graph-cut algorithm for minimization

The variational problem (9.4) is very challenging to solve. We show in this paragraph that, after discretization of the surface and of the set of rays, it can be transformed into a minimum cut problem on a particular graph. By computing the minimum cut using available ecient graph-cut libraries, we obtain a fast method to solve the surface segmentation problem.

The surface is represented by an elevation map E : (x, y) → z = E (x, y) (which guarantees that it is representable as an elevation map). The horizontal location (x, y) and the elevation z are discretized. To make an easier connection between the elevation map and the surface it denes, we consider the layer cake decomposition of the elevation.

With this decomposition, a discrete elevation map corresponds to a binary volume (a discrete version of the epigraph of E ) and the boundary in that volume denes the discrete surface.

We build a graph as depicted in gure 9.2, with a node to represent each voxel of the binary volume of E . Two special nodes, called the source (denoted 's') and the sink (denoted 't') are added in order to simulate a ow from the source to the sink. Nodes are connected together by directed edges with specic capacities and a ow is said to be admissible if and only if the ow along each edge is non negative and smaller or equal to the edge capacity, and there is no ow accumulation/creation at nodes (except at the source and at the sink). By the max-ow min-cut theorem, algorithms that identify the maximum admissible ow on the graph can also identify the minimum cost cut among all possible cuts in the graph 1 , see for example [START_REF] Boykov | An experimental comparison of min-cut/maxow algorithms for energy minimization in vision[END_REF]. During the graph construction, by creating edges with well-chosen capacities, we can make the cost of any cut exactly match the cost of the corresponding surface in the variational formulation (9.4).

To represent the rst term in equation (9.4), we substitute D ray with its denition in equation ( 9.3):

ray∈R

D ray (r ray→S ) dR = ray∈R r ray→S r min C - ray (r s ) -C + ray (r s ) + dr s dR + ray∈R rmax r ray→S C + ray (r s ) -C - ray (r s ) + dr s dR . (9.5)
Each of the two terms corresponds to summations over a half-space whose boundary is S : the half-space that contains the radar and the half-space with the farther ranges, respectively. We add an edge directed from the source to node i, the node that represents the 3-D position (x i , y i , z i ) and that is located at the distance r i from the radar antenna.

The capacity 2 of this edge is set to C - i (r i ) -C + i (r i ) + , where C - i and C + i are the cumulative reectivities computed along the ray directed from the radar through the 1 the cost of a cut is the sum of the capacities of all edges cut that are directed from a node in the source partition to a node in the sink partition 2 note that an edge with zero capacity can be suppressed because it carries no ow and has no contribution to the cost of the cuts of the surface is obtained by adding bi-directional edges between neighboring nodes in the x and y directions. (c) to prevent the cut from severing twice a column of nodes along the z direction, ascending edges with innity capacity are added. These edges are counted in the total cost of the cut only when they go down-stream: from the partition containing the source to the partition containing the sink. point of coordinates (x i , y i , z i ). Another directed edge is added from node i to the sink, with capacity C + i (r i ) -C - i (r i ) + . To separate the graph into two parts by a cut, some edges must be severed (unless the cut passes through the distance of equilibrium r equi ) and the sum of the capacities of those edges corresponds to a discretization of equation (9.5), see Fig. 9.3(a).

Additional edges are created to account for the regularization term β A (S ): bidirection edges between pairs of nodes that are direct neighbors in the x or y directions, with capacity β, see Fig. 9.3(b). Finally, ascending edges with innite capacity are included between neighboring nodes in the z direction. These edges are necessary to guarantee that the cut denes a surface that is representable by an elevation map, see 9.3(c). Similar edges are added in Ishikawa's graph construction that is also based on the layer-cake decomposition [START_REF] Ishikawa | Exact optimization for markov random elds with convex priors[END_REF].

In our implementation, we computed eciently the summations along the rays by resampling the reconstructed tomographic volume in ray geometry so that sums could be carried out along columns in this new geometry. For the construction of the graph and the computation of the minimum cut, we used the graph-cuts library by Boykov and Kolmogorov [START_REF] Boykov | An experimental comparison of min-cut/maxow algorithms for energy minimization in vision[END_REF].

Joint reconstruction and surface segmentation

The knowledge provided by the segmented urban surfaces can help to improve the inversion described in chapter 8. The reconstruction algorithm that can most readily be extended to include segmented surfaces is the 3-D inversion method described in equation (8.3). Under the assumption that the signal retrieved over urban areas is mainly constituted of punctual bright points, sparsity may be an ecient enough prior to obtain clean tomograms. Nonetheless, this implies that the tuning of the sparsity parameter be done locally according to the position of the scatterers. In CS for SAR tomography, the sparsity constraint is generally set locally in the range and azimuth direction but constant for each radar cell. Here we propose to use the 3-D information provided by the estimated surface to go one step further and perform a spatially varying penalization of the sparsity.

When applying CS or the 3D inversion, the sparsity parameter µ is set proportional to the level of spurious elements in the reconstruction. Generally µ is set according to the noise level (Zhu and Bamler, 2012b), but as decorrelation mechanisms and side-lobes should also be discarded, the knowledge of the sensor thermal noise may not be enough.

Many SAR tomographic algorithms propose to estimate the number of backscattering elements in order to extract the largest scatterers in each radar cell. This step cleans the estimated tomograms from residual outliers, but is also a challenging task for large multitemporal stacks in dense environments. Moreover the CS approach may then lose one of its asset with respect to MUSIC or WSF if it also needs an estimation of the number of targets.

Under the assumption that the location of the urban surface is known, the sparsity parameter µ can be spatially tuned to lead to rened tomograms. Even when the surface is roughly known, it provides information on where the reconstructed signal should be located. In the proposed iterative algorithm, µ is computed as a function of the distance to the surface in the 3-D space and the number of iterations: The rened tuning of the sparsity according to the surface allows to considerably improve the scatterers localization and main lobe reduction. In some cases, however, the segmented surface follows the lobe main extension direction and is not as localized as would be expected for a collection of point-like scatterers. In the global reconstruction of the scene, most of the artifacts due to the TV penalization are suppressed after 10 iterations. The obtained surface is then very close to the ground truth and provides the lowest error according to table 9.1.

Experiments

To validate both the generality of the segmentation method and its eciency on real A third experiment presents the evolution of the reconstructed slice as the REDRESS algorithm iterates cf. Fig. 9.8. It can be observed that the distribution of reectivities becomes much sharper after a few iterations.

To estimate the covariance matrix at each point, we used a 7 × 7 Gaussian lter.

For MUSIC and WSF , the number of scatterers is set constant and equal to 2 to avoid selecting too many outliers while allowing multiple scatterers within each radar resolution cell. For these two estimators, the reectivity is estimated by mean square minimization, to keep a physical interpretation of the tomograms. As the scene is very heterogeneous with a lot of layover, this step introduces some undesired mixing of the information in the image. The surfaces estimated from tomographic reconstructions using spectral analysis techniques present noticeable artifacts in the dense areas. Some structures are too extended, partially lling streets or the building atrium. Meanwhile, the averaging step makes the tomographic estimation smoother in homogenous areas for the fully sparse approaches MUSIC and WSF .

The CS technique performed on all the data set presents results that seem visually the closest to the ground truth. Many details can be observed in this reconstruction: most of the rooftops and visible streets are well segmented and the buildings atrium are also retrieved.

For all the previous estimators, the TV minimization produces some building elongation resulting in phantom structure in low intensity signal area. This can be seen for instance in the bottom right part for Fig. 9.6 or around the position 450 for Fig. 9.5. Google Earth © cf. g. 9.11.

To conduct a quantitative comparison of the segmentation results, we report the mean error for each estimated surface to the ground truth cf. 9.1. The TV parameter β is set, for each method, as the one minimizing this error.

Conclusion

In this chapter we introduced a graph-cut based segmentation algorithm to estimate the urban surfaces from a SAR tomographic reconstruction. The proposed approach is very general and can be used in combination with many dierent tomographic algorithms.

Conclusion and perspectives

Conclusion

The main objective of this PhD was to provide methods to use the structural information to perform the SAR tomographic estimation of the 3-D reectivity for dense urban scenes. Two dierent strategies were explored:

Exploit the redundancy in the SAR images

Regularize the estimation

The rst method was based on the non-local estimation of the covariance matrix.

This approach proved to be challenging and didn't oer guaranties to outperform basics local covariance matrix estimation approaches.

The second approach was expected to be relatively successful as it was close to the well known CS approach which may be consider as one of the top tomographic estimator for urban areas. The proposed algorithm allowed to perform the tomographic 3-D inversion of a scene under geometrical priors described in the ground coordinate system. This allowed to reduce the number of outliers and favor structured distribution of scatterers.

Tomograms and volume representation of the reectivity may be hard to analyze, specially with close buildings of comparable height. We proposed a post-processing graph-cut based segmentation to extract the urban surfaces of a tomographic reconstruction. Due to the very generic design of the method, it can used in association to any estimator providing a measure of the scene reectivity (amplitude or power). The surface is expected to be retrieved were high intensity voxels are found. The geometrical priors consist of a TV minimization under constraints avoiding the structures to be transparent for the electromagnetic wave and crenelation like distortions.

Finally, we linked the 3-D reconstruction with the surface segmentation, to design an iterative algorithm that rene the reectivity estimation based on the estimated surface.

All the proposed method were tested on a stack of TerraSAR-X images and the results were compared to a ground truth. The REDRESS algorithm lead to a average error less than 2m and preserves many details such as buildings courtyard or some roof reliefs.

  Contribution (3): Segmentation de surface urbaine par graph-cut La troisième contribution présentée dans ces travaux est un algorithme de segmentation de surface urbaine dans des cubes tomographiques. Il est attendu que la surface réponde à plusieurs contraintes : La surface est proche des voxels de forte intensité la surface est lisse et composée de zones planes selon les axes horizontaux et verticaux Chaque rayon partant du capteur vers la scène intersecte la surface exactement une fois La surface n'intersecte pas plus d'une fois chaque ligne verticale La fonction de coût correspondant à ces contraintes est non-convexe et dicile à optimiser. Toutefois ce problème après discrétisation peut s'écrire comme la recherche d'une coupe minimale dans un graphe bien construit. L'algorithme de segmentation proposé est très général au sens où n'importe quelle estimation proposant une mesure physique de la réectivité (amplitude ou intensité) peut être utilisée en entrée. Associée à une approche parcimonieuse type CS, les résultats obtenus sont particulièrement prometteurs. XVI CONTENTS L'information apportée par la surface est, de plus, utilisable pour régler les paramètres de certains algorithme. CS est en eet sensible au réglage du paramètre de parcimonie. Contribution (4): Reconstruction 3-D et segmentation de surface alternées La quatrième contribution de cette thèse est l'association naturelle des deux dernières. Nous proposons ainsi un algorithme itératif où après chaque reconstruction 3-D, la surface estimée est utlisée pour raner le réglage du paramètre de parcimonie. Ainsi, la parcimonie est calculée comme une fonction croissante de la distance à la surface. De cette façon, plus on est loin de la surface, plus il est coûteux de trouver de la matière. Peu d'itérations susent pour corriger les tomogrammes d'erreurs résiduelles et pour diminuer les artefacts générés par la minimisation de la variation totale. Conclusion Au moins quatre contributions originales sont présentées dans ces travaux. Si le premier axe de recherche n'a pas pu aboutir à une réelle amélioration des techniques d'inversion tomographique, le deuxième axe a été plus fructueux. Ainsi, les approches développées basées sur une régularisation structurelle exprimée en géométrie sol ont permis d'une part d'obtenir des résultats équivalants ou au dessus de l'état de l'art et de poser un

Figure 2 . 1

 21 Figure 2.1 Geometry of acquisition of a SAR image.

Figure 2 . 2

 22 Figure 2.2 Construction of the synthetic antenna. A point seen multiple times during the sensor motion can be localized precisely along the x axis by combining all the received echoes.

Figure 2 . 3

 23 Figure 2.3 Construction of the synthetic antenna. The phase variation induced by the motion of the sensor is exploited for the azimuth compression.

Figure 2 . 4

 24 Figure 2.4 The value of each pixel in the Single-Look Complex (SLC) SAR image corresponds to the coherent summation of the back-scattered signals inside a given radar cell of resolution.

Figure 2 . 5

 25 Figure 2.5 The gray areas cannot be seen by the sensor as it is hidden by the stair step structure. In the blue areas dierent elements are projected in the same radar resolution cells producing a layover eect.

Figure 2 . 6

 26 Figure 2.6 Distortions induced by the side-looking sensor. When β > 0 • the order is preserved (a) whereas it is inverted when β > 0 • (b).

  coecients of the linear phase perturbation across the image, ξn 2π = 2bn

Figure 2 .

 2 Figure 2.7 SAR tomographic geometry of acquisition. The scene on the ground is seen multiple times from a slightly dierent angle by the sensor at each pass. The phase dierence between the images is then exploited to retrieve the 3-D distribtution of the reectivity of the scene.

Figure 3 . 1

 31 Figure 3.1 Models depending on the structure of the back-scattering elements: distributed scatterers like vegetation or rough objects back-scatter the wave in all the directions. The signal from these scatterers obeys the unconditional model where no stability assumption can be made between the acquisitions. Dihedral or trihedral objects reect most of the wave in the sensor direction. These scatterers follow the conditional model. The signal retrieved from distributed scatterers is generally of low intensity and poorly localized compared to conditional ones.

Figure 3 . 2

 32 Figure 3.2 Example of objects back-scattering a signal following the unconditional model in a SAR image obtained from the sensor SETHI of the ONERA. This model is well suited for natural rough elements such as trees or grass. Depending on their granularity and the used bandwidth, roads may back-scatter some signal following this model.

3. 2

 2 Unconditional (distributed) signal model The rst model named unconditional or stochastic signal model corresponds to elements presenting a varying reectivity across the images. Such signals are generally reected by rough at surfaces or volumic scatterer distributions as seen in Fig. 3.2. Without

Figure 3 . 3

 33 Figure 3.3 Example of objects back-scattering a signal following the conditional model in a SAR image obtained from the sensor SETHI of the ONERA. Most of these strong bright scatterers come from man-made objects such as buildings, cars or pylons.

Figure 3 . 4

 34 Figure 3.4 Example of objects back-scattering a signal following the hybrid model in a SAR image obtained from the sensor SETHI of the ONERA.

  Urban areas are, by denition, lled with many man-made structures presenting many straight walls, metallic objects, corners, etc. The resulting SAR images present a high level of bright point-like scatterers. Many of those scatterers can be considered as stable in time and are correctly described by the conditional signal model. This model or the hybrid one are then generally used for describing back-scattering mechanisms over

  Fig 4.1 illustrates conventional Beamforming performances for scatterer localization and on real data.

Figure 4 . 1

 41 Figure 4.1 Unmixing power of the conventional beamforming. The rst 4 images present the estimation of two scatterers separated by 2, 1, 0.5 and 0.25 times the ambiguity height. The position of the scatterers is indicated by the dashed red lines and the blue curves is the estimated intensity P BF . The last row illustrates the estimated signal intensity for a simulation where the scatterers are well separated and a real urban slice. The ground truth is indicated in white.

Fig 4. 2

 2 Fig 4.2 illustrates conventional Beamforming performances for scatterer localization and on real data.

Figure 4 . 2

 42 Figure 4.2 Unmixing power of the Capon beamforming. The rst images present the estimation of two scatterers separated by 2, 1, 0.5 and 0.25 times the ambiguity height. The position of the scatterers is indicated by the dashed red lines and the blue curves is the estimated intensity P C . The last row illustrates the estimated signal intensity for a simulation where the scatterers are well separated and a real urban slice. The ground truth is indicated in white..

Fig 4 .

 4 Fig 4.3 illustrates MUSIC performances for scatterer localization and on real data.

Fig 4 .

 4 Fig 4.4 illustrates WSF performances for scatterer localization and on real data.

Figure 4 . 3 Figure 4 . 4

 4344 Figure 4.3 Unmixing power of MUSIC through its pseudo-spectrum. The rst 4 images present the estimation of two scatterers separated by 2, 1, 0.5 and 0.25 times the ambiguity height. The position of the scatterers is indicated by the dashed red lines and the blue curves is the pseudo-spectrum P M U value. The last row illustrates the estimated signal for a simulation where the scatterer are well separated and a real urban slice. The ground truth is indicated in white.

Fig

  Fig 4.5 illustrates ML performances for scatterer localization and on real data.

  Figure 4.5 Unmixing power of ML according to formula 4.27. The rst 4 images present the estimation of two scatterers separated by 2, 1, 0.5 and 0.25 times the ambiguity height. The position of the scatterers is indicated by the dashed red lines and the blue dot the estimated scatterers position and intensity. The last row illustrates the estimated signal for a simulation where the scatterer are well separated and a real urban slice.

Fig

  Fig 4.6 illustrates M-RELAX performances for scatterer localization and on real data.

  Figure 4.6 Unmixing power of M-RELAX . The rst 4 images present the estimation of two scatterers separated by 2, 1, 0.5 and 0.25 times the ambiguity height. The position of the scatterers is indicated by the dashed red lines and the blue dot the estimated scatterers position and intensity. The last row illustrates the estimated signal for a simulation where the scatterer are well separated and a real urban slice. The ground truth is indicated in white.

Fig 4 .

 4 Fig 4.7 illustrates M-RELAX performances for scatterer localization and on real data.

  Figure 4.7 Unmixing power of SPICE according to formula 4.27. The rst 4 images present the estimation of two scatterers separated by 2, 1, 0.5 and 0.25 times the ambiguity height. The position of the scatterers is indicated by the dashed red lines and the curve is the estimated intensity. The last row illustrates the estimated signal for a simulation where the scatterer are well separated and a real urban slice.

  Figure 4.8 Unmixing power of CS . The rst 4 images present the estimation of two scatterers separated by 2, 1, 0.5 and 0.25 times the ambiguity height. The position of the scatterers is indicated by the dashed red lines. The last row illustrates the estimated signal for a simulation where the scatterer are well separated and a real urban slice. The ground truth is indicated in white. For the simulation the SNR corresponding to the additive white noise is 0.3 dB.

Figure 4 . 9

 49 Figure 4.9 Diagram representing the estimators presented in this part around three axes: sparsity, MOS and covariance based.

Fig

  Fig. 4.9.

Figure 5 .Figure 5

 55 Figure 5.1 shows an example of reconstructed proles illustrating the decorrelation, distortions and outliers issues. In this example, a building facing the sensor (the Mirabeau Tower) sends back a strong signal whereas most of the other scatterers on the prole are relatively low intensity. The reconstruction obtained with conventional beamforming is then highly impacted with high side lobes propagating in the radar cells were the building is projected. When looking at the conventional beamforming spectrum as well as the MUSIC-PS , it appears the signal estimated is determined by the brightest scatterers in it. This side lobe phenomena combined with the main lobe size may also produce outliers even with sparse approaches such as MUSIC or CS .In the rst part, state-of-the-art estimators for SAR tomography have been described. As stated in the preamble of chapter 3, with these methods, the estimation of the scatterers height and reectivity is done pixel wise. 3-D SAR tomography is then seen as a concatenation of local estimations. Dense urban areas are however very structured as most of scatterers are distributed along roughly vertical or horizontal plans.The inherent diculty induced by dense city landscape can then be reduced by taking into account the geometrical behavior of the scene. Urban 3-D reconstruction should

  Figure 6.1 Two simulations of a simple 3-D scene: (a) considering a regular sampling of trajectories along the elevation axis (Reg case); (b) using the trajectories of TerraSAR-X satellites (TSX case). For each scene, two 3-D views are presented (left and middle) as well as the temporal mean intensity image (right).

Figure 6 . 2

 62 Figure 6.2 Illustration of the dient elevation samplings in the Reg case (left) and TSX case (right).

Fig. 6 .

 6 Fig.6.1 shows the theoretical distribution of the scatterers. Both scenes are composed of a ground at a constant altitude, a wall and a roof, resulting in a large layover area in the SAR images. The simulated building is higher than the estimated elevation resolution (as given by Fourier inversion) in the TSX case. The reectivity of the scatterers is set constant on all the scene. The SNR for the additive noise level is 1.4 dB and is the same for both experiments. In the Reg case, the scatterers are well separated in the azimuth and range directions. In the TSX experiment, the density of scatterers is larger in each cell, resulting in clusters of neighboring scatterers being projected in the same radar cell, as illustrated in the Fig 6.2. The distance between the scatterers is set as a third of the cell resolution. All the TSX images are thus corrupted with a speckle eect induced by the coherent addition of the back-scatterered signals. Due to the relatively

  angle θ = 0.6 rad, a wavelength λ = 0.0311 m and a distance R = 6.1510 5 m. The characteristics of the sensor and of the scene are summarized in the table 6.1. The scene is very heterogeneous with dierent back-scattering mechanisms. Most of the scene is composed of buildings producing very bright point-like echoes. The majority of the constructions are composed of either 20 m to 40 m height buildings or characteristics of the sensor and the observed scene 80 m to 110 m skyscrapers. Other very typical Parisian monuments can also be observed such as the Eiel Tower in the upper right corner or the Maison de la Radio in the left side of the image. The Seine river with many bridges crossing it ows from the top right to the bottom left corner. As the water is very smooth, almost all the signal is scattered in the specular direction which creates these typical very low intensity areas in the image. Many vegetation areas are visible in the optical image but they produce very low intensity and low coherency signal in the SAR images. Moreover most of the trees answers are projected in the same radar cell as neighboring buildings (layover phenomena) and are thus hidden by brighter scatterers. Vegetation is still visible in some places such as on the island (l'ile aux Cygnes) where no buildings are present.

Figure 6 . 3

 63 Figure 6.3 Up: temporal intensity average of the stack of SAR images, bottom: corresponding optical view of the same zone.

Figure 6 . 4

 64 Figure 6.4 Annotated SAR and optical images of the observed scene. The position of the sensors with respect to the scene is indicated for both images. The typical distortions induce by the SAR ranging acquisition system are well visible as the tall structures such as the Eiel Tower or the Mirabeau Tower are projected along the range direction.

Figure

  Figure 6.5 Temporal and spatial baselines

Figure 6 . 6

 66 Figure 6.6 Accuracy (a) et completeness (b) criteria. The accuracy computes the mean of the error for each estimated point. The completeness gives an indication on the proportion of holes in the rendering.

Figure 6 . 7

 67 Figure6.7 Evaluation protocol: each RADAR cell is described by the set of its maxima. Dierent thresholds are applied to further select the points which are then compared to the ground truth. This operation is done for all the cells of the image.

2 S 2 S

 22 are the coecients of the isotropic lter and the subscript p corresponds to a 2-D shift p ∈ [-H, H] 2 indicating the support of the lter. The number of samples is then L = (2 × H + 1) 2 . Although these lters are fast and easy to compute, they lack of spatial adaptivity and induce resolution loss in the images. Moreover, as SAR images have a very high dynamic range, bright pixels may have amplitudes several orders of magnitude larger than the background. These basic local lters then produce blurry estimation of point-like bright pixels.Other local methods like the Lee lter[START_REF] Lee | Speckle analysis and smoothing of synthetic aperture radar images[END_REF] are designed to adapt locally on the content inside the support of the function. The Lee lter is built to reduce the speckle phenomena while preserving sharp structures in the image. This lter is designed to minimize the mean square error for a linear speckle noise model. The Lee lter gives then as output for each channel n of the tomographic image i.e. each 2-D SAR image:În = I n + k Lee I n -I n (7.2)where I n = |v(n)| 2 is the intensity of the noisy pixel, I n the mean intensity within the lter window and În the estimated intensity. The coecient k Lee adjusts the local smoothing and is computed as: standing for the theoretical coecient of variation of the speckle: γ

  .1), p ∈ [-H, H] indicates a local 2-D shift for a window of size L = (2 × H + 1) 2 . The quantity ∆ is then a dissimilarity indication between patches.

Figure 7 . 1

 71 Figure 7.1 Images of buidings near Nimes-Garon airport obtained by the Sethi airborn sensor -ONERA. From left to right : noisy SLC SAR image and its denoised version using a gaussian kernel, the Lee lter or NL-SAR algorithm.

Figure 7 . 2

 72 Figure 7.2 Images of buildings in the south-west of Paris obtained by the TerraSAR-X spaceborn sensor -DLR. From left to right : noisy SLC SAR image and its denoised version using a gaussian kernel, the Lee lter or NL-SAR algorithm.

Figure 7 . 3

 73 Figure 7.3 Deterioration of the denoising power of NL-SAR for an increasing number of images. From top left to bottom right the number of channels is 2, 4, 6, 8, 10 and 12.

Figure 7

 7 Figure 7.4 Deterioration of the discriminative power of the GLR for increasing number of images. From top left to bottom right the number of channels is 2, 4, 6, 8, 10 and 12. On top the amplitude image and one interferogram. The central pixel in red in the amplitude image is compared to the rest of the window.

Figure 7 . 6 Figure 7 . 7 Figure 7 . 8

 767778 Figure 7.6 Comparison of the selective power of the proposed similarity and the Wishart GLR. On top the amplitude image and one interferogram. The central pixel indicated in red in the amplitude image is compared to the rest of the window. On the middle the similarity map obtained with Wishart GLR and the 2000 most similar pixels. On the bottom the similarity map obtained with the proposed method and the 2000 most similar pixels.

Figure 7 . 9 Figure 7 .

 797 Figure 7.9 Comparison of the selective power of the proposed similarity for increasing number of SAR images in the tomographic pile for a pixel of dense buildings. From top left to bottom right the number of images is N = 10, 20, 30, 40. The selection of the pixels is robust to the increasing size of the images stack: most of the selected pixels are located at the same height and on the same building.

REGULARIZATIONFigure 8

 8 Figure 8.1 The operator Φ performs the projection of the voxels from the 3-D scene in ground geometry on the SAR tomographic stack.

Figure 8 . 2

 82 Figure 8.2 Construction of the matrix Φ. The radar cell associated to each position of the 3-D scene is computed as well as the phase term corresponding to the distance between a point in the 3-D space and a given antenna.

Figure 8 . 3

 83 Figure 8.3 Illustration of the non-convexity of spatial regularizations expressed on the modulus of the complex reectivities.

Figure 8 Figure 8

 88 Figure 8.4 Inuence of the parsimony µ 1 parameter. In (a) the parameter is way above its optimal value resulting in a lot of holes in the structure. In (b) we use a small value of µ 1 resulting in a high number of outliers.

  Figure 8.9 Study of the inuence of each parameter on the metrics for the TSX scene (cf. Fig. 6.1 (b)). The plots in (a), (b), (c) and (d) correspond respectively to dierent values of µ 1 , µ x , µ y and µ z with the three other parameters being xed.

Figure 8

 8 Figure 8.10 Best representation on the Reg scene using classical beamforming (a), Capon beamforming (b), MUSIC (c), CS (d) and our tomographic 3D inversion (e).

  Figure 8.11 Best representation on the Reg scene using classical beamforming (a), Capon beamforming (b), MUSIC (c), CS (d) and our tomographic 3D inversion (e).

  Figure 8.13 Best representation on the TSX scene using classical beamforming (a), Capon beamforming (b), MUSIC (c), CS (d) and our tomographic 3D inversion (e).

  Figure 8.14 Accuracy vs. completeness for classical beamforming, Capon beamforming, MUSIC, Compress Sensing and our 3-D Inversion algorithm. In (a) the evaluation is done on the Reg scene and in (b) on the TSX scene cf. 6.1.

Figure 8

 8 Figure 8.15 Accuracy as a function of completeness, computed based on the ground truth displayed in (a), for the following tomographic estimators: classical beamforming, Capon beamforming, MUSIC, CS and the proposed 3D inversion. In (c) the evaluation is done over the area corresponding to the red rectangle in (b). In (d) the test area is inside the green rectangle.

Figure 8

 8 Figure 8.16 Visualization of the module of the reconstructed cube (Andriy Fedorov et al., 2012) from a stack of 40 TerraSAR-X images of Paris (cf. Fig.6.3) using the conventional CS algorithm.

Figure 8

 8 Figure 8.18 Visualization of the module of the reconstructed cube (Andriy Fedorov et al., 2012) from a stack of 40 TerraSAR-X images of Paris (cf. Fig.6.3) using the 3D inversion algorithm.

Figure 8

 8 Figure 8.20 Visualization of the point cloud obtained from the reconstructed cube and projected in Google Earth © .

Figure 9 . 1 Figure 9 . 2

 9192 Figure 9.1 We seek a surface S that, for each ray, is close the scatterer(s) found along the ray.

Figure 9 . 3

 93 Figure 9.3 The capacities of the edges are chosen so that the cost of the cut corresponds to the energy of the surface. (a) the delity to the tomographic reconstruction is enforced via edges originating from the source or leading to the sink. (b) the spatial smoothness

  µ k (p, S ) = µ 0 + b (n -1) p, S ) is the Euclidean distance from the point p = (x, y, z) T to the estimated surface, k is the current iteration and n the total number of iterations. We dene by µ(S ) ∈ R Nx.Ny.Nz the 3-D sparsity parameter map. As the surface location estimation may be subject to errors in the rst iterations, it is important to avoid over-penalizing points moderately close to the surface during the rst reconstructions. This why we use the square of the distance d multiplied b a factor smaller than 1. As the number of iterations increases, the reconstruction and thus the surface estimation should be more accurate (and better in match) which suggests an increase of the penalization of the distance from a reconstructed voxel to the surface. µ 0 + b is then the desired minimal sparsity that need to be applied to voxel not on the surface.The proposed iterative reconstruction and surface segmentation is summarized in the

2 (

 2 Figure 9.4 Observed urban areas : optical image (a), temporal mean of the corresponding SAR image (b), and the 3-D model from Google Earth used as a ground truth. The red line in (b) and white one in (c) correspond to the slice shown in Fig. 9.5

  Figure 9.5 Urban surface estimation using graph-cut segmentation of the tomograms, as described in section 3. The estimated surface corresponds to the red prole. The ground truth for the given slice is shown in green. The tomograms are obtained using Capon beamforming (a), SPICE (b), MUSIC (c), WSF (d), the 3-D inversion approach (e) and REDRESS (f ).

Figure 9

 9 Figure 9.6 Ground truth height (a), scene surface estimation using SPICE (b), MUSIC (c), WSF (d), Capon beamforming (e), the 3-D inversion (f ) and REDRESS (g). For each results, the image shows the surface colored according to its height.

  Experiments done on a set of 40 TerraSAR-X images of Paris show good results for dierent tomographic estimators(Capon beamforming, MUSIC , WSF , SPICE ,. As the 3-D inversion algorithm is designed to use 3-D priors, we also present an algorithm that alternatively reconstructs the 3-D distribution of reectivities, segments the urban surfaces from the volume of reectivities and updates the regularization so as to improve the subsequent 3-D reconstruction. While the noniterative 3-D inversion algorithm fails in some cases to reduce the main lobes of the strong scatterers, the alternating scheme achieves a much sharper estimation of the distribution of reectivities.

  Figure 9.8 Three dierent iteration steps from the alternate reconstruction algorithm. On the left column, the estimated reectivities are shown for the prole presented in Fig 9.4. On the right, the estimated surface (red) and the ground truth (green) are superimposed in addition to the estimated reectivities. Rows (a), (b) and (c) correspond respectively to the rst, third and fth iterations (last one).

  Figure 9.10 Visualization of the module of the reconstructed cube using the REDRESS algorithm (left), super-imposed with the ground truth (right) (Andriy Fedorov et al., 2012).

  Figure 9.11 Plotting of the non-zero intensity scatterers in Google Earth ©. For computational purpose 1 voxel out 10 is drawed.

  Titre: Approches structurelles pour la tomographie SAR haute résolution en milieu urbain dense Mots clés: SAR, Tomographie, Problèmes inverses, parcimonie, graph-cut Résumé: La tomographie SAR exd'un point de vue légerement diérent pour reconstruire la densité complexe de réectivité au sol. Cette technique d'imagerie s'appuyant sur l'émission et la réception d'ondes électromagnétiques cohérentes, les données analysées sont complexes et l'information spatiale manquante (selon la verticale) est codée dans la phase. De nombreuse méthodes ont pu être proposées pour retrouver cette information. L'utilisation des redondances naturelles à certains milieux n'est toutefois généralement pas exploitée pour améliorer l'estimation tomographique. Cette thèse propose d'utiliser l'information structurelle propre aux structures urbaines pour régulariser les densités de réecteurs obtenues par cette technique. Title: Structural Approaches for SAR tomography over dense urban areas Keywords: SAR, Tomography, Inverse problems, sparsity, graph-cut Abstract: SAR tomography consists in exploiting multiple images from a same area acquired from a slightly different angle to retrieve the 3-D distribution of the complex reectivity on the ground. As the transmitted waves are coherent, the desired spatial information (along the vertical axis) is coded in the phase of the pixels. Many methods have been proposed to retrieve this in-formation in the past years. However, the natural redundancies of the scene are generally not exploited to improve the tomographic estimation step. This Ph.D. presents new approaches to regularize the estimated reectivity density obtained through SAR tomography by exploiting the urban geometrical structures. Université Paris-Saclay Espace Technologique / Immeuble Discovery Route de l'Orme aux Merisiers RD 128 / 91190 Saint-Aubin, France

  expected to be smaller than 1. Depending on the ratio 2 4 σ 2

	n
	µ d ν d

, L G may be very low. Now if the d th columns of V 1 and V 2 are dierent but not orthogonal, the previous equation becomes:

Table 8 .

 8 -data experiments conducted over a very heterogeneous scene indicate that the reconstruction doesn't present artifacts due to excessive smoothing in some areas. To conrm this observation with a controlled experiment, reconstruction scores of the TSX scene with random scatterers amplitude are also tested. The amplitudes were chosen randomly for each azimuth in a given range. The dierent experiments correspond increasing ranges. The results of the experiments are shown in table 8.2. The proposed method still performs better, even when the amplitudes of the scatterers varies by 3 orders of magnitude and the regularization parameters are kept constant. 2 Evolution of the MACT value when the scatterer amplitude is xed uniformly for each azimuth. The experience is repeated for increasing amplitude interval size.

	Estimator		MACT
	Beamforming		0.96
	Capon Beamforming		0.98
	MUSIC			0.66
	CS			0.71
	3-D inversion (µ * )		0.57
	3-D inversion (µ * -50%)	0.62
	3-D inversion (µ * -100%) 0.62
	3-D inversion (µ * + 50%)	0.59
	3-D inversion (µ * + 100%) 0.65
	3-D inversion (µ * + 200%) 0.74
	Table 8.1 Evolution of the Minimum Accuracy-Completeness Tradeo (MACT) when
	the parameters are getting far from their optimal value µ * = {µ * 1 , µ * x , µ * y , µ * z }. The
	smaller the MACT the better.
	Estimator/Amplitudes [1, 1] [10 -1 , 10] [10 -2 , 10 2 ] [10 -3 , 10 3 ]
	Beamforming	0.96	1.09	1.10	1.25
	Capon Beamforming	0.98	1.35	1.38	1.58
	MUSIC	0.66	0.66	0.68	0.69
	CS	0.71	0.71	0.72	0.72
	3-D inversion	0.57	0.58	0.60	0.60

Remerciements

and the corresponding likelihood ratio :

12 , u

(1)

12 , z

(1)

2 )

(7. 26) The log of the generalized likelihood ratio test is computed by injecting the parameters values in (7.26) and takes the following form:

1 )û

(1)

2 )û

(1)

12 )û

(1)

where ||x|| 2 M = x H M x. The procedure can be repeated for all the scatterers contained in the tested pixels, e.g for the k th iteration the following hypothesis are tested:

where the covariance matrix, the altitudes and amplitudes are updated at each iteration.

As just a few strong scatterers are expected to be present in each pixel, only three targets are tested for each pixel. The similarity criterion is then given by the sum of all the log likelihood ratio : (7.28) This ensures to take into account the similarity between all the dierent contributions and avoids to have any parameter to tune. Through the algorithm iterations and the parameters updates, dierent decompositions of the signal into stable strong scatterers and residual are jointly tested. Moreover as the rst contributions are several magnitude higher than the following ones, the value of L is conditioned by the rst elements in the sum in (7. to well select the pixels containing iso-height scatterers whereas the proposed similarity detect the correct samples on the building. More important, the shadow area is clearly seen as dissimilar as the building which is not the case with Wishart GLR. In Fig. 7.7, REGULARIZATION that our method dominates the other approaches for all accuracy/completeness tradeos.

The 3-D plain representation of the scene obtained by our approach is presented in Fig. 8.18. The heat color is chosen to be proportional to the intensity of the voxels normalized by the averaged intensity of the corresponding radar cell. This normalization is used only to help the visualization as the dynamic range of SAR images is very high.

In By using the georeferencing of the voxels, we can project the corresponding point cloud into Google Earth ©. The results are presented in Fig. 8.20. This step allows us to have more details in the visualization of the scene. We can now see that some points above the Mirabeau Tower are relevant: they correspond to the structure of its rooftop.

Our reconstruction method also correctly identied parts of the structures on the dock and several buildings hidden behind the Mirabeau tower.

Part of the PhD was also dedicated to the study and phase calibration to a set of airborne SAR images obtained with the sensor SETHI from the ONERA. This work is however still in progress and the data are not usable for SAR tomography yet.

Perspectives

SAR tomography is still a relatively recent research eld. Its application to dense urban areas is a trending topic and new algorithms are proposed either to perform the tomographic inversion or to apply some post-processing regularization. The framework presented here to perform the tomographic inversion in ground coordinates allows to used many dierent regularization strategies in a plug-and-play fashion. Some of alternative to the presented priors could then be tried:

Going back to the smoothing presented in chapter 8, strategies allowing to know the local orientation of the scatterers distribution could be used to locally tune the gradient parameters. This could be done rst as an improve REDRESS algorithm.

Rather than using the quadratic norm of the gradient, minimizing the TV for the tomographic reconstruction should also be tried.

Besides trying sophisticated priors, the 3-D projection operator could be used to retrieve a SAR image associated to a given sensor track. This could be useful for several applications:

The proposed methods may introduce an increase of the number of parameters.

We discussed how to set these parameters in chapter 8 but this remain an open question. This consideration is also true for conventional CS . A tuning strategy could be to nd the set of parameters that allows to best preserve the input SAR images.

Being able to exploit a stack of SAR images to generate a new one could be used for changes detection or image restoration.

When dealing with a new set of SAR images, phase calibration can be a challenging task. An idea could be to perform iteratively the calibration and the tomographic reconstruction. At each step the linear contribution of the APS and the scatterers height could then be estimated using the information from the previous iteration to enhance the results.