Student Member, IEEE Clément Rambour
email: clement.rambour@telecom-paristech.fr

Loïc Denis
email: loic.denis@univ-st-etienne.fr

Senior Member, IEEE Florence Tupin
email: florence.tupin@telecom-paristech.fr

Hélène Oriot
email: helene.oriot@onera.fr

Introducing spatial regularization in SAR tomography reconstruction

Keywords: tomographic SAR inversion, dense urban areas, TerraSAR-X, Compressed Sensing, inverse problems, 3-D reconstruction I

come

scatterers at different 3-D locations and thus reconstruct the 3-D reflectivity of the scene. The 3-D representation of cities has many applications such as crisis management, movement monitoring or change detection.

The spatial resolution of modern sensors such as TerraSAR-X or COSMO-SKYMED can be better than one meter, which is much finer than the average height of observed buildings in city areas. At these resolutions, urban surfaces (fac ¸ades, roofs) can be recovered. However, as the dense urban scenes contain a high number of back-scattering objects, scatterers from different structures are projected within the same radar resolution-cell. The actual tridimensional location of these scatterers can be reconstructed by combining images from repeated passes with slightly different angles of observation, by SAR tomography. SAR tomography can be seen as an extension of the SAR focusing to 3 dimensions, the separation of scatterers projected within a given radar cell being performed thanks to a synthetic aperture formed by the different trajectories of the SAR antenna in the tomographic stack [START_REF] Reigber | First demonstration of airborne SAR tomography using multibaseline L-band data[END_REF], [START_REF] Moreira | A tutorial on synthetic aperture radar[END_REF].

In the ideal case of equi-spaced trajectories, focusing in the direction orthogonal to the line of sight can be simply performed by application of the inverse discrete Fourier transform. 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, hence 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 [START_REF] Stoica | Introduction to Spectral Analysis[END_REF], [START_REF] Gini | Multibaseline cross-track SAR interferometry: a signal processing perspective[END_REF], leading to satisfying results on homogeneous areas dominated by volumic scattering such as forests or landscape areas [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]. These estimators rely on the analysis of the covariance matrix. The estimation of this matrix is not a trivial task because it involves some form of spatial averaging [START_REF] Deledalle | NL-SAR: A Unified Nonlocal Framework for Resolution-Preserving (Pol)(In)SAR Denoising[END_REF] and becomes very challenging in dense urban areas due to the high spatial heterogeneity [START_REF] Deledalle | Exploiting Patch Similarity for SAR Image Processing: The nonlocal paradigm[END_REF], [START_REF] Rambour | Similarity Criterion for SAR Tomography over Dense Urban Area[END_REF]. Rather than exploiting the covariance matrix, some works on very high resolution data use directly the time domain backprojection of the signal [START_REF] Tebaldini | Imaging the internal structure of an alpine glacier via l-band airborne sar tomography[END_REF]. When backscattering elements are discrete rather than continuous, the well known Compressed Sensing (CS) approach can produce super-resolved tomographic reconstruction by performing the backprojection of the data under sparsity constraints [START_REF] Budillon | Three-Dimensional SAR Focusing From Multipass Signals Using Compressive Sampling[END_REF], [START_REF] Zhu | Tomographic SAR Inversion by L 1 -Norm Regularization ; The Compressive Sensing Approach[END_REF]. Since this approach does not require the estimation of a covariance matrix, it involves no spatial filtering and is more effective at reconstructing multiple closeby point scatterers.

CS appears to be well suited to the inversion of stacks of tomographic SAR images obtained over dense urban areas. It can however be criticized for its lack of spatial regularity: pixels are processed independently while many urban structures are organized in lines or planes. Including a spatial regularization can be expected to improve the robustness to noise, reduce the ambiguities and better preserve the signal coming from the buildings. Moreover, the sparse estimation step of the signal is generally hard to tune. Due to the high dynamic of SAR images, it often appears that either weak parts of the signal are suppressed during the estimation, or outliers are detected. To refine the results, the volume reconstructed with CS is generally post-processed by estimating the number of scatterers in order to select only the most significant points, see [START_REF]Super-Resolution Power and Robustness of Compressive Sensing for Spectral Estimation With Application to Spaceborne Tomographic SAR[END_REF], [START_REF] Fornaro | Compressive sensing and generalized likelihood ratio test in SAR tomography[END_REF]. The result is then a set of discrete points defined by their 3-D localization and complex reflectivity. As the distribution of the reflectivity in dense urban configurations is hard to model, this approach often fails to select points with low reflectivities. Close to CS, the recent FAST-SUP-GLRT detector [START_REF] Budillon | A fast support detector for superresolution localization of multiple scatterers in sar tomography[END_REF] 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 doesn't take into account the geometry of the scene, the structure of the data can be more accurately represented than with conventional CS.

A recent strategy to enhance the discrete reconstructions obtained with CS consists of post-processing the SAR tomographic (post-processed) points clouds. It is mainly based on shape detection and surface fitting using priors on the urban geometry. It has been shown to recover some of the information lost during point selection and even to fill missing parts of the scene [START_REF] Zhu | Facade Reconstruction Using Multiview Spaceborne TomoSAR Point Clouds[END_REF] [START_REF] Ley | Regularization and Completion of TomoSAR Point Clouds in a Projected Height Map Domain[END_REF]. Those approaches give satisfying geometrical 3-D modeling of urban environments but lack from a direct connection to the data. In this paper, we consider an alternative approach that enforces spatial regularity during the tomographic inversion step. Since the output of our method is a sparse volume, a discrete set of points can be extracted and any pipeline already designed to process tomographic point clouds can possibly be applied as a post-processing.

In this paper, we present a new algorithm that performs the tomographic 3-D inversion of a stack of SAR images combining sparse priors as in CS and geometrical ones using a ground coordinate back-projection operator as in [START_REF] Tebaldini | Imaging the internal structure of an alpine glacier via l-band airborne sar tomography[END_REF]. To derive our algorithm, we account for two inherent characteristics of SAR imaging: the geometrical distortions induced by the geometry of acquisition and the use of Single Look Complex (SLC) data. First, we present in details the concept of SAR tomography and the used notations (section II). Then, we address the problem of using priors corresponding to the natural properties of the urban environment. The third section is dedicated to solving the tomographic inverse problem, which is a nonconvex and high-dimensional optimization problem (section III). In the fourth section, we illustrate that our algorithm performs well both on synthetic data with controlled sensor parameters, scene composition and geometry and on real data (section IV). The quality of our results is evaluated using the accuracy and completeness defined in [START_REF] D'hondt | Impact of non-local filtering on 3D reconstruction from tomographic SAR data[END_REF] [START_REF] Hondt | Nonlocal Filtering Applied to 3-D Reconstruction of Tomographic SAR Data[END_REF]. We use a stack of 40 TerraSAR-X images of the city of Paris to test our algorithm on real data. The results are evaluated on some buildings using a 3-D modeling of the scene given by the French Geographic Institute (IGN) as ground truth.

II. SAR TOMOGRAPHY

A. Forward model for SAR tomography in radar geometry

A SAR tomographic stack consists of N SAR Single Look Complex (SLC) images acquired in interferometric configuration. Each SAR image of the stack corresponds to a different 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. To simplify the notations, in the following this master image will correspond to the first image of the stack. The geometry of a simple 3-D scene and of the sensor trajectories is illustrated in Fig. 1.

The complex value v n (x, r) obtained at pixel coordinates (x, r) after SAR synthesis of the n-th image corresponds the convolution of the complex 3-D reflectivity 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 n (x, r) = f (x -x , r -ρ n;y,z )u(x, y, z). exp - 4jπ λ ρ n;y,z + jϕ atmo dx dydz + (1)
To better identify real-valued and complex-valued variables in the description of our algorithm (section III), we underline all complex-valued variables throughout the paper. Also the x and azimuth axis coincide to avoid complicated notations without loss of generality. Here x and r stand for the discrete coordinates of the focused data, ρ n;y,z is the distance between the n th sensor and a scatterer at a position (y, z) for a given x. ϕ atmo is the phase shift corresponding to the Atmospheric Phase Screen (APS). f corresponds to the PSF and depends of the sensors and the processing of the data. is a white additive Gaussian noise.

In this paper, we will use two classical approximations in SAR tomography. Firstly, the PSF is substituted with a boxcar function that may have a near zero width depending on the knowledge of the data. Then, to be in the same framework as most of the other state of the art tomographic approaches, we use the interferometric images with respect to the master image as the input tomographic stack. The phase difference is then approximated as a linear function of the elevation of the scatterers [START_REF] Fornaro | Three-dimensional focusing with multipass SAR data[END_REF]. To remove the APS contribution we use the method detailed by Ferretti and Al. [START_REF] Ferretti | Permanent scatterers in SAR interferometry[END_REF]. Without assuming movement of the scatterers, the equation (1) becomes the sum of the complex amplitudes scattered by all objects that project within the resolution cell [START_REF] Budillon | Three-Dimensional SAR Focusing From Multipass Signals Using Compressive Sampling[END_REF], [START_REF] Zhu | Tomographic SAR Inversion by L 1 -Norm Regularization ; The Compressive Sensing Approach[END_REF], [START_REF] Fornaro | Three-dimensional focusing with multipass SAR data[END_REF]

v n (x, r) = (y,z)∈∆r u(x, y, z)exp -jξ n h(y, z) dydz + (2)
The r-th radar resolution cell is defined by ∆ r = {(y, z) | rδ r /2 ≤ ρ 1;y,z ≤ r + δ r /2 and z ≤ z max (y)} with δ r the step Fig. 1. Configuration of the tomographic acquisition of a 3-D scene. The scene is naturally described by the coordinate system (x, y, z) whereas the data are described in the sensor coordinate system (x, r, h).

in range and z max (y) the maximum height for a point at (x, y) location to be illuminated by the radar. The chosen approximations allow us to use classical algorithms to compensate the APS and to reduce the size of the operators we will then define to perform the inversion. Here as we are also interested in a comparison with other SAR tomographic estimators we chose to keep this classical interferometric framework. Moreover, using directly the equation (1) can also be done but the phase calibration may be harder.

We denote by x, y and z the coordinates of a reflector in 3-D space and by x, r and h its coordinates in radar space. In equation ( 2), ∆ r is the extension of the radar resolution cell along the elevation direction h, u(x, y, z) is the complex reflectivity of the scatterer located at position (x, y, z), and h(y, z) its elevation. The parameter ξ n = 4πbn λr is the spatial frequency associated to the sampling of the scene for each baseline. b n is the baseline n and λ the radar wavelength. Random variable models the complex additive white noise due to thermal noise.

Consider a given radar resolution cell. We will drop in the following the indices identifying the spatial location of that cell. We form the vector u = u 1 • • • u D T that collects the complex reflectivities of all scatterers seen within the same radar resolution cell and

h = h 1 • • • h D T their elevation.
After discretization of (2), we can express v ∈ C N , the collection of SLC values in the tomographic stack at the pixel of interest, as a linear combination of the complex signal backscattered by each of the D scatterers:

v = A(h)u + n , (3)
where A(h) ∈ C N ×D is the steering matrix. The d-th column [A(h)] d of A(h) 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 (4)

B. Spectral Analysis in SAR Tomography

Many strategies exist to invert (3) in order to estimate the heights h and reflectivities u. Most of them come from the Direction Of Arrival (DOA) literature and have been successfully used for SAR tomography on homogenous scenes such as forest or ice. Non-parametric methods such as beamforming or Capon filters are fast to implement and offer a good global representation of the reconstructed scene [START_REF] Stoica | Introduction to Spectral Analysis[END_REF], [START_REF] Gini | Multibaseline cross-track SAR interferometry: a signal processing perspective[END_REF], [START_REF] Capon | High-resolution frequency-wavenumber spectrum analysis[END_REF], [START_REF] Porfiri | Building profile reconstruction using TerraSAR-X data time-series and tomographic techniques[END_REF]. Parametric methods achieve super resolution by exploiting priors on the estimated data. In many cases, a limited number of point-like scatterers is assumed. When considering a given number of such scatterers, the Maximum Likelihood (ML) estimator can be applied to estimate the location and reflectivity of each scatterer. Exact computation of the ML estimator requires a prohibitive computation time when more than 2 or 3 scatterers are considered, because of the combinatorial nature of the optimization problem. Suboptimal methods are therefore generally used specially when dealing with big images. MUSIC [START_REF] Schmidt | Multiple emitter location and signal parameter estimation[END_REF] is one of the most popular parametric estimators used to estimate the position of a finite number of sources. The localization of the D sources is obtained by finding the steering vectors with minimal norm in the eigen-subspace of R associated to the (N -D) smaller eigenvalues. Since this approach is adapted, by construction, to separate discrete scatterers, it is not well suited to the reconstruction of continuous reflectivity profiles. In urban areas, it offers fair performances when the reflectivity profiles can be accurately represented by sparse distributions. The pseudo-spectrum given by the inverse of the projection is also sometime used to analyze even continuous tomographic profiles [25] [26]. Very similar to MUSIC, but derived as an approximation of ML, Weighting Subspace Fitting (WSF) [START_REF] Viberg | Detection and estimation in sensor arrays using weighted subspace fitting[END_REF] methods are designed to achieve an even better resolution than MUSIC. By using both the signal and the noise subspace, it has been proven that WSF can be used to separate the distribution of the leaves in a forest and an artificial target underneath [START_REF] Huang | Under foliage object imaging using SAR tomography and polarimetric spectral estimators[END_REF]. All those methods require to estimate the covariance matrix at each pixel, which can be a very hard task in heterogeneous regions such as dense urban areas. Since it takes at least N 2 pixels to get a full-rank sample covariance matrix, the more images in the tomographic stack, the strongest the spatial filtering necessary to estimate the covariance matrices (hence, the strongest the range and azimuth resolution loss).

C. Compressed Sensing in SAR Tomography

The use of Compressed Sensing (CS) for SAR tomography is recent. It is well suited to the urban environment [START_REF] Reale | Tomographic Imaging and Monitoring of Buildings With Very High Resolution SAR Data[END_REF] as it does not require to estimate R. The estimation of the reflectivities profile 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 0 ) u (5)
As the vector of heights h is one of the unknown, the elevation axis is sampled into D 0 bins forming the array h 0 . The solution to the combinatorial problem (5) can be approximated using the classical convex relaxation of the 0 pseudo-norm into an 1 norm:

û = arg min u 1 2 ||A(h 0 )u -v|| 2 2 + µ 1 ||u|| 1 (6)
where µ 1 is a Lagrange multiplier. Despite its very good performances on urban areas, this approach has some drawbacks. Indeed, the matrix A(h 0 ) does not in general guarantee to satisfy either the Restrictive Isometry Property (RIP) or an incoherence condition [START_REF] Zhu | Super-Resolution Power and Robustness of Compressive Sensing for Spectral Estimation With Application to Spaceborne Tomographic SAR[END_REF] required to obtain reconstruction guarantees when using the relaxed problem [START_REF] Tebaldini | Multibaseline Polarimetric SAR Tomography of a Boreal Forest at P-and L-Bands[END_REF]. When choosing a large number of samples D 0 in the height direction in order to achieve a super-resolution, these conditions are largely violated. RIP essentially indicates that any random sub-matrix of A(h 0 ) composed from D random columns of A(h 0 ) is nearly orthogonal which preserves the energy of any D-sparse vector u. The incoherence insures that a highly localized signal gets spread over all the columns of A(h 0 ) which makes the reconstruction robust to random under-sampling. The consequences of the violation of RIP and incoherency are the creation of outliers and/or multiple detections of the same scatterer in adjacent positions.

In order to refine the results obtained via a sparse reconstruction, one can include spatial priors that describe the urban geometry. As illustrated in figure 1, when performing a tomographic inversion in the radar geometry, natural geometrical properties (e.g., predominance of vertical and horizontal structures and of straight angles) are not straightforward to express. Rather than using the radar frame (r, h), we suggest in this paper to perform the tomographic inversion in the ground frame (x, y, z) in order to include structural information in the form of priors (i.e., regularization terms).

III. PROPOSED APPROACH: TOMOGRAPHIC INVERSION WITH SPATIAL REGULARIZATIONS

A. Forward model for SAR tomography in ground geometry

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 reflectivities in 3-D space to the complex amplitudes in the tomographic stack of SAR images, see Fig 2 . Since we chose to align the x axis and the azimuth, an element of Φ is defined as:

Φ i,j =      exp(-jϕ) if x i = x j and r i -δr 2 < ρ 1;yj ,zj < r i + δr 2 , 0 otherwise. (7)
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;yj ,zj 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. 3. Matrix Φ is sparse: only a few entries are different from zero, so that products of the form Φu can be computed efficiently.

The observed SAR tomographic stack v can be modeled by the following (complex-valued) linear model:

v = Φ u + ( 8
)
where stands for the noise. This corresponds to a generalization of the tomographic direct model ( 3) 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.

B. Spatial regulariations for tomographic inversion

It is often desired to reconstruct volumes with a discretization in heights that is finer than the resolution given by the synthetic aperture in the height direction (i.e., superresolution). The inversion of equation ( 8) 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 specific regularization can be enforced on the phase of our unknown complex reflectivities u. We define the regularization R as a function of the modulus of u only. The reconstruction of the volume of complex reflectivities û is thus obtained by solving an optimization problem of the general form:

û = arg min u 1 2 ||Φ u -v|| 2 2 + R(|u|) (9)
To design the regularization function R, we need to select a function that favors volumes of reflectivities |u| that are often present in urban environments. Many different 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 filled 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 (10)
where the matrices D x , D y and D z stand for the finite differences 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 finite difference operators enforce a spatial smoothness. As µ 1 controls this behavior, it plays a crucial role in the reconstruction.

In urban environment, 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 first 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 supposed to show straight walls and then a natural direction of elongation along the z axis. 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, the conjoint minimization of the 1 norm of the reflectivity and the square 2 norm of its gradient tends to favor sparse and continuous reflectivity profile. The obtained reconstructions present then less outliers to the benefit of structural elements.

When reconstructing images of urban areas, we found that it was impossible to find 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 Unified Nonlocal Framework for Resolution-Preserving (Pol)(In)SAR Denoising[END_REF] or simply using the average intensity depending of the configuration of the acquisitions. We define 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 verifies r -δr 2 < ρ 1;yj ,zj < r + δr 2 ). The equation ( 10) is then modified 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 (11)
where 1 is the vector of size N x .N y .N z with each entry equal to 1.

C. Optimization algorithm for the 3-D reconstruction

The minimization problem ( 9) 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. 4: 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 fixed value of u 1 ) corresponds to a set of complex values with equal modulus (the white circle drawn on Fig. 4).

To solve the minimization problem ( 9), 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: w ∈ R +(Nx.Ny.Nz) ). Problem ( 9) is formally equivalent to the following constrained problem:

û = arg min u 1 2 Φ u -v 2 2 + R(w) (12)
s.t. u = f |f | = w
To solve this constrained optimization problem, we apply the method of multipliers [START_REF] Bertsekas | Constrained optimization and Lagrange multiplier methods[END_REF], i.e., we look for a saddle point of the augmented Lagrangian L defined by:

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) (13)
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). We follow the hierarchical approach 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 (14) = β 1 • |u -d 1 | + β 2 • (w + d 2 ) β 1 + β 2 + exp j • arg u -d 1 (15)
where . + is the identity on [0, +∞[ and is the constant null function on ] -∞, 0[.

Proof.

The second term in ( 14) only depends on the modulus of f . The phase of f * is then driven by the first 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 (16)
ρ 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 [START_REF] Budillon | A fast support detector for superresolution localization of multiple scatterers in sar tomography[END_REF].

This leads to the following algorithm 1 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 reflectivities) 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 * (û,
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
1 compared to the well-known alternating directions method of multipliers (ADMM) [START_REF] Boyd | Distributed optimization and statistical learning via the alternating direction method of multipliers[END_REF], 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 differentiable (since w ≥ 0, w 1 = 1 T w, which is differentiable; 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 efficient 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 2 . 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 sufficient since the algorithm is warm-restarted. In the definition 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 when applying the chain rule since ∂C ∂[f * ]i = 0 for all i. 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, the convergence is assumed when we met a local minimum. In practice the algorithm can then be stopped after a fixed number of iterations (around 60) or when the distance between two iterations is 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.

IV. EXPERIMENTS

A. Evaluation Protocol

In this part we compare 3-D scenes reconstructed by five different approaches: classical beamforming [START_REF] Stoica | Introduction to Spectral Analysis[END_REF], Capon beamforming [START_REF] Capon | High-resolution frequency-wavenumber spectrum analysis[END_REF], MUSIC [START_REF] Schmidt | Multiple emitter location and signal parameter estimation[END_REF], CS [START_REF] Zhu | Tomographic SAR Inversion by L 1 -Norm Regularization ; The Compressive Sensing Approach[END_REF] [11] and our algorithm. Evaluation of different SAR tomographic estimators can be hard as they originate from various fields and are thus not 2 https://github.com/emmt/OptimPackLegacy suited to the same applications. Moreover, they do not always estimate the same quantities: classical and Capon beamforming give the estimated power of the scatterers, CS and our 3-D inversion estimate the scatterers reflectivity, and MUSIC only provides us with their position. Here, we will only evaluate the quality of the localization, as all methods can be used to estimate the 3D location of scatterers and this represents a rich information in urban areas. It can be compared to a ground truth (digital elevation model).

We use the accuracy and completeness criteria introduced in [START_REF] Hondt | Nonlocal Filtering Applied to 3-D Reconstruction of Tomographic SAR Data[END_REF] to evaluate the tomographic SAR reconstructions. Those two metrics are defined for the evaluation of point clouds reconstructions. Tomographic inversions must thus be converted into a discrete representation in order to compute these performance criteria.

1) 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 ||p j -p k || 2 (17)
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. Accuracy indicates whether reconstructed points are correctly located.

2) Completeness: The completeness corresponds to the mean distance from each point of the ground truth to the points in P:

C( P, P) = 1 N p Np k=1 min j ||p j -p k || 2 , (18)
with N p the number of points in the ground truth. Completeness indicates whether the ground truth is well represented by the set of points in the reconstruction. These two metrics are illustrated in Fig. 5 where green dots represent estimated points pj and blue squares the ground truth points p k . Accuracy and completeness provide 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)).

The accuracy and completeness can be jointly used to evaluate the quality of a discrete reconstruction. Beforehand, 3-D points must be extracted from the volume of voxels obtained by the tomographic reconstruction algorithms. We performed this extraction by finding the local maxima in each radar resolution cell. This produces a first point-based representation of the data. Then, to reduce the sensitivity to noise and side-lobes, we applied a threshold on the obtained point cloud. Starting from a reconstructed volume û, we obtain the collection of 3-D points Pû,t :

Pû,t = p j = (x j , y j , z j ) ∈ R 3 , such that |û xj ,yj ,zj | ≥ |û xj ,νy,νz | and |û xj ,yj ,zj | > t (19)
where ν y and ν z are the coordinates of the direct neighboring voxels restricted to the same radar resolution cell and t is a threshold. In the following experimentss, when we visually compare reconstructions obtained by different algorithms, we always select the threshold t as the value t * that achieves the best trade-off between accuracy and completeness:

t * = arg min t A( Pû,t , P) 2 + C( Pû,t , P) 2 . ( 20
)
When evaluating different tomographic estimators, we are generally interested in different regimes (accurate reconstruction of the strongest scatterers, dense reconstruction of most scatterers). To capture these different cases, rather than considering a single reconstruction obtained with threshold t * , we also represent the accuracy score A( Pû,t , P) 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). Many works have considered linear dependencies of the phase of permanent scatterers on various parameters such as their height but also their speed or the temperature [START_REF] Budillon | Extension of a Fast GLRT Algorithm to 5D SAR Tomography of Urban Areas[END_REF] [START_REF] Reale | Tomographic Imaging and Monitoring of Buildings With Very High Resolution SAR Data[END_REF]. If the reconstruction of surfaces is the aim, completeness should be favored, even if this degrades the accuracy.

B. Simulated Data

In this section, we compare our 3-D tomographic reconstruction algorithm to state-of-the-art methods, using simulated data. We simulate two stack of images corresponding to two different configurations 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. 17 (referred to in the following as TSX). As in [START_REF] Zhu | Tomographic SAR Inversion by L 1 -Norm Regularization ; The Compressive Sensing Approach[END_REF], the decorrelation effect is introduced by adding a Gaussian white noise and the phase noise is modeled by adding a uniformly distributed random phase between -π and π. Fig. 6 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 reflectivity of the scatterers is set constant on all the scene. The additive noise level is of -0.6 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 bigger for each cell resulting in clusters of neighboring scatterers being projected in the same radar cell. The signal is thus corrupted with a speckle effect induced by the coherent addition of the backscatterered signals. Due to the relatively strong additive noise and the geometrical decorelation the average coherance 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 and the resulting images have a bigger size in range.

C. Influence of the regularization parameters

To show the influence of each regularization parameter, we performed different reconstructions of the Reg scene for various sets of regularization values (in particular, with one parameter chosen so as to illustrate the effect of either underregularization or over-regularization). In these expreriments, the phase noise is not present to show clearly the structural influence of the spatial smoothing. The results are presented in figures 7, 8, 9 and 10. On this well-sampled scene with a good resolution, the parameter with the largest influence is µ 1 , associated to the sparsity constraint. Fig. 7 illustrates that when µ 1 is too large, there are some holes in the reconstruction, while a value of µ 1 that is too small leaves side-lobes and outliers. The effect of over-smoothing (values of µ x , µ y or µ z too large) is the extension of structures in the direction of the spatial smoothing. This is visible in particular with structures whose orientation differs from that of the smoothing, see in Fig. 9 the widening of the wall due to excessive smoothing in the horizontal direction. The same effect appears on the rooftop in Fig. 10. Insufficient smoothing translates into residual fluctuations (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. 9 and 10. 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 flexible and well-suited both to the reconstruction of urban areas where scatterers are sparsely distributed over 2-D surfaces and to vegetated areas where scatterers are spread in volumes in a smooth 3-D distribution. The downside of this flexibility 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 offers the best performance, as measured by the Minimum Accuracy / Completeness Tradeoff (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-off. We advice to perform the one dimensional optimization, by importance: µ 1 , then µ z , then µ x and µ y . While tuning the sparsity parameter, the other smoothing parameters should be put to 0.

Other derivative free methods can be used to fix all the parameters at once. We compared the described alternating minimization strategy with Nelder-Mead's simplex method by initializing with different sets of parameters and found similar values for the lowest error. The differences for µ 1 is of 21% and below 14% for the smoothing parameters. As being more precise, this approach is highly computational and needs to be reinitialize several times due to the non-convexity of the error with respect to the parameters. In the absence of ground truth, a simple numerical simulation using the same geometrical configuration and sensor noise power can be generated in order to automatically tune the parameters. If a simulated scene isn't 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 different sets of parameters. These curves are drawn for our two simulation cases Reg and TSX respectively in Fig. 11 and12. We observe that the best accuracy /completeness trade-off (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], [START_REF] He | Alternating Direction Method with Self-Adaptive Penalty Parameters for Monotone Variational Inequalities[END_REF]. Here, we fixed those parameters to 10 after having tested different values.

D. Comparaison of different tomographic estimators

The figures 13 and 14 give the 3-D reconstructions obtained using several state-of-the-art tomographic reconstruction methods and our spatially regularized approach: (a) classical beamforming, (b) Capon beamforming, (c) the parametric estimators MUSIC, (d) Compressed Sensing, and (e) our algorithm. The corresponding curves of accuracy as a function of completeness are presented in figure 15 for the two scenes.

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. In the TSX case where the theoretical resolution is worst, few of the wall is retrieved leading to a poor completeness score. A closer inspection of the TSX reconstructions indicates that 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 sup-press 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 effects on the estimation of the position of the scatterers.

In the TSX scene, the best accuracy / completeness trade-off from the MUSIC pseudo-spectra produces a sparse scene with no outliers. However the bias present in the estimation of the covariance matrix strongly deforms the wall and introduces All the spectral estimators presented (classical beamforming, Capon beamforming and MUSIC) suffer from the averaging step needed to estimate the covariance matrix, which produces visibly distorted walls and/or rooftop in the TSX reconstruction. An elongation of the rooftop and a reduction of the ground is also visible in the Reg case. However the naive covariance estimation step makes those estimators efficient for uniform areas which is illustrated by the generally accurate ground reconstruction.

As expected, CS is one of the top estimator 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. As the method uses a sparse prior, some information may be lost but thanks to the geometrical constraints, most of the rooftop and ground are filled. Moreover, in the TSX case the wall and rooftop are better located around their true position and the reconstruction is more robust to deformation introduced by the proximity of the scatterers.

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 modified with increasing perturbation values. The results are shown in the table I. Up to an isotropic perturbation of 50, the proposed method still shows a better Minimum Accuracy-Completeness Tradeoff (MACT) than the other tested methods. An other concern is whether the parameters need to be adapted to the scatterers power. Real-data experiments conducted over a very heterogeneous scene indicates that the reconstruction doesn't present artifacts due to excessive smoothing in some areas. To confirm 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 different experiments correspond increasing ranges. The results of the experiments are shown in table II. 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. azimuth resolution 0.87 m. The observed scene is presented in Fig. 16 3 . The spatial and temporal baselines are shown in Fig. 17. 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 :

δ z = δ h sin(θ) = λr sin(θ) 2∆b = 6.99m (21)
with an incidence angle θ = 0.6 rad and wavelength λ = 0.0311 m. Sparse reconstructions in urban SAR tomography have been shown to significantly improve this resolution [START_REF] Zhu | Very High Resolution Spaceborne SAR Tomography in Urban Environment[END_REF].

To evaluate the different tomographic estimators we use a rough ground truth of the scene. The ground truth surface is shown in Fig. 18 4 . 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. 18.

It is noteworthy that MUSIC performs particularly well on areas containing a lot of ground surface. In the first test, the sub-scene is composed of a mixture of tall buildings, 3 These data have been provided by the DLR in the framework of the project LAN1746. 4 Ground truth on Paris complimentary provided by the IGN, France. small structures near the ground level and smooth flat areas (streets and docks). As the diffuse signal back-scattered by flat surfaces is much weaker than the one reflected 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. 18, 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 again be due to its ability to retrieve more points on the ground and the small gap between the two methods needs to be taken with caution as the used ground truth is coarse.

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 that in all the possible ways of selecting the reconstructed voxels, our tomographic 3-D inversion algorithm performs a better reconstruction than the compared estimators. For this scene the parameters are µ 1 = 12, µ z = 300, µ y = 100 and µ x = 100 and were computed from the TSX simulations.

The 3-D plain representation of the scene obtained by our approach is presented in Fig. 21. 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 Fig. 22, the same image is shown except that the near zero intensity voxels are made fully transparent and the ground truth is superimposed. The results obtained on the same scene with conventional CS are shown in Fig. 19 and Fig. 20. The colormap are the same for both representations and are chosen to have the dark blue color associated with intensities below 10 -4 times the maximum voxel amplitude. This visually set to zero dark voxels that still have a non-zero intensity as the algorithm may not have fully converged. An other way to have a sparser representation of the modulus of the signal is to display the image w at convergence of the algorithm. In the 3-D representation, most of the isolated outliers are suppressed thanks to the spatial smoothing and the buildings where most of the high reflectivity voxels are massed are retrieved. The smoothing denoising however comes with the cost of slightly blurred scatterers. We can see that most of the very bright buildings present in the SAR images in Fig. 16 are well reconstructed. The two towers that were mostly missing in the 2D intensity image have a point-like representation but are fairly visible in this 3-D rendering.

By using the georeferencing of the voxels, we can project the corresponding point cloud into Google Earth c . The results are presented in Fig. 23. 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 identified parts of the structures on the dock and a large amount of buildings hidden behind the Mirabeau tower.

V. CONCLUSION

In this paper a new tomographic 3-D inversion method has been proposed. This method is shown to be effective to inverse a SAR tomographic stack in a dense urban area. The chosen spatial regularization strategy takes advantage of both the natural sparsity of targets in urban environments and the spatial smoothness (i.e., spatial proximity of scatterers). Our approach generalizes compressed sensing techniques and includes, in an adaptive way, a spatial smoothing.

To be able to compare our method with other estimators used in SAR tomography we extended the protocol of evaluation based on accuracy and completeness proposed in [START_REF] Hondt | Nonlocal Filtering Applied to 3-D Reconstruction of Tomographic SAR Data[END_REF]. The tests done on numerical simulations and on a stack of TerraSAR-X images in a urban area showed that our method leads to better results than existing approaches in terms of accuracy / completeness compromise. This 3D inversion algorithm opens a door to more sophisticated approaches taking into account different priors describing the data. Indeed, the main structure of the tomographic 3D inversion would remain unchanged if other regularization terms were considered within a variable splitting approach. The research of appropriate priors to scenes with a volumetric density of reflectivity such as ice or forest will be the subject of further works.

Fig. 2 .Fig. 3 .

 23 Fig. 2. The operator Φ performs a projection of voxels from the 3-D scene in ground geometry on the SAR tomographic stack.

Fig. 4 .

 4 Fig.4. Illustration of the non-convexity of spatial regularizations expressed on the modulus of the complex reflectivities.

Fig. 5 .

 5 Fig. 5. 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.

Fig. 6 .

 6 Fig.6. Two simulations of a simple 3-D scene: first row, considering a regular sampling of trajectories along the h axis; second row: using the trajectories of TerraSAR-X satellites. Sub-figures (a), (b), (d) and (e) show two views of the scenes while sub-figures (c) and (f) correspond to the average of the simulated SAR images. The sub-figures in the top row correspond to the Reg case and the ones in the bottom correspond to the TSX case. The Reg case consists of 20 images acquired with uniformly distributed sensors along the elevation axis. In the TSX case the geometry of the scene is given by the true TerraSAR-X satellites position for each of the 40 acquisitions. In both case the SNR is 1.7dB. In the TSX experiment, more scatterers are projected in the same resolution cell which produces a important speckle effect.

Fig. 7 .Fig. 8 .Fig. 9 .Fig. 10 .

 78910 Fig.7. Influence 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.

Fig. 11 .Fig. 12 .

 1112 Fig. 11. Study of the influence of each parameter on the metrics for the Reg scene (cf. Fig. 6 (a)). The plots in (a), (b), (c) and (d) correspond respectively to different values of µ 1 , µx, µy and µz with the three other parameters being fixed.

  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

Fig. 13 .

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

Fig. 14 .

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

Fig. 15 .

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

Fig. 16 .Fig. 17 .

 1617 Fig. 16. (a) Image of the temporal mean of all the iamges in the stack and (b) the corresponding optic image obtained from Geoportail (IGN).The building back-scattering a strong signal in the center af the image is the Mirabeau Tower. In the bottom left of the image, the structure presenting a strong periodic back-scattering signal is the Ministry of Foreign Affair.

Fig. 18 .Fig. 19 .Fig. 20 .Fig. 21 .Fig. 22 .

 1819202122 Fig. 18. Scores accuracy versus completeness based on the ground truth in (a) for the 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.

Fig. 23 .

 23 Fig. 23. Visualization of the point cloud obtained from the reconstructed cube with projected in Google Hearth c .

TABLE I EVOLUTION

 I OF THE MINIMUM ACCURACY-COMPLETENESS TRADEOFF (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

TABLE II EVOLUTION

 II OF THE MACT VALUE WHEN THE SCATTERER AMPLITUDE IS FIXED UNIFORMLY FOR EACH AZIMUTH. THE EXPERIENCE IS REPEATED FOR INCREASING AMPLITUDE INTERVAL SIZE.

VI. ACKNOWLEDGEMENTS

This project has been funded by ANR (the French National Research Agency) and DGA (Direction Générale de l'Armement) under ALYS project ANR-15-ASTR-0002 and by the RTRA Digiteo Labs.