Y Yan

M P Doin
email: doin@geologie.ens.fr

P López-Quiroz

F Tupin

B Fruneau
email: fruneau@univ-mlv.fr

V Pinel

E Trouvé

TOPICS IN APPLIED EARTH OBSERVATIONS AND REMOTE SENSING 1 Mexico City subsidence measured by InSAR time series: Joint analysis using PS and SBAS approaches

In multi-temporal InSAR processing, both the Permanent Scatterer (PS) and Small BAseline Subset (SBAS) approaches are optimized to obtain ground displacement rates with a nominal accuracy of millimeters per year. In this paper, we investigate how applying both approaches to Mexico City subsidence validates the InSAR time series results and brings complementary information to the subsidence pattern. We apply the PS approach (Gamma-IPTA chain) and an ad-hoc SBAS approach on 38 ENVISAT images from November 2002 to March 2007 to map the Mexico City subsidence. The subsidence rate maps obtained by both approaches are compared quantitatively and analyzed at different steps of the PS processing. The inter-comparison is done separately for low-pass (LP) and high-pass (HP) filtered difference maps to take the complementarity of both approaches at different scales into account. The inter-comparison shows that the differential subsidence map obtained by the SBAS approach describes the local features associated with urban constructions and infrastructures, while the PS approach quantitatively characterizes the motion of individual targets. The latter information, once related to the type of building

foundations, should be essential to quantify the relative importance of surface loads, surface drying and drying due to aquifer over-exploitation, in subsoil compaction.

Index Terms

InSAR time series, subsidence, SBAS, PS, joint analysis I. INTRODUCTION Differential interferometry, using the phase difference between two radar images taken at two different dates combined with a digital elevation model (DEM), provides a measurement of the ground displacement with a pluri-centimetre accuracy [START_REF] Massonnet | Radar interferometric mapping of deformation in the year after the Landers earthquake[END_REF] [START_REF] Zebker | On the derivation of coseismic displacement-fields using differential radar interferometry -The Landers earthquake[END_REF]. This technique has been successfully applied to the monitoring of landslides, earthquake deformations, volcanic activities and urban subsidence. Its limitations result from DEM errors, atmospheric propagation delays of the radar wave and decorrelation due to the increase of the temporal and spatial baseline between satellite passes. To overcome these difficulties and produce long time series of ground motion, the Permanent Scatterer (PS) and the Small BAseline Subset (SBAS) approaches have been developed.

Permanent Scatterer Interferometry (PSI) [START_REF] Ferretti | Nonlinear Subsidence Rate Estimation Using Permanent Scatterers in Differential SAR Interferometry[END_REF] [4] [START_REF] Colesanti | Monitoring landslides and tectonic motions with the PS technique[END_REF] distinguishes itself from other SAR interferometric processing by the use of a single master image to generate a stack of differential interferograms without limitations in temporal or spatial baselines. PS candidates, which a priori carry reliable phase information across the interferogram stack, are selected based on their backscattering properties. On these points, the PS approach adopts essentially a model-based, temporal unwrapping strategy. Accordingly, a priori information on the displacement is necessary, from which a deformation model can be established. In general, the average linear displacement velocity and the DEM error are considered as the two major parameters of 2D linear regression of the wrapped phases. As temporal unwrapping is performed on local phase differences, the PS approach includes schemes to integrate in space relative displacement rates and DEM errors.

The SBAS approach [START_REF] Berardino | A new algorithm for surface deformation monitoring based on small baseline differential sar interferograms[END_REF] [START_REF] Schmidt | Time-dependent land uplift and subsidence in the santa clara valley, california, from a large interferometric synthetic aperture radar data set[END_REF] [8] [START_REF] Lanari | A Small-Baseline Approach for Investigating Deformations on Full-Resolution Differential SAR Interferograms[END_REF] increases the spatial coverage over which one extracts reliable phase delay time series, especially outside urban areas, by taking the speckle properties of most targets in SAR images into account. To maximize coherence, interferograms are computed only for image pairs separated by small temporal and spatial baselines. Interferograms form a redundant network linking between images in the temporal and spatial baseline space. Decorrelation noise in the interferograms is partly removed by range filtering of the non-overlapping part of the spectrum and by applying a spatial filter, thus reducing [START_REF]Interferometric point target analysis for deformation mapping[END_REF] DRAFT the interferogram spatial resolution. Interferograms are then spatially unwrapped. The inversion of the whole set of interferograms by Singular Value Decomposition (SVD) provides phase delay time series.

These two approaches are applied here to a data set centered on the Mexico City basin, a large, flat, endoreic basin (at 2240 m elevation) surrounded by volcanic chains. Its subsoil is classically divided into three main geotechnical zones: Foothills Zone (I), Transition Zone (II) and Lake Zone (III) (Figure 1).

The foothills subsoil consists of heterogeneous volcanoclastic deposits and lava. In contrast, the lake zone subsoil corresponds to highly compressible lacustrine clays. In between, the transition zone is mainly defined by sand and gravel alluvial deposits, intercalated with volcanic materials and clay lenses. The Mexico City subsidence, which is characterized by a wide spatial extent, very large rate (reaching 38 cm/yr), and extreme subsidence gradients (variation of up to 15 cm/yr in 200 m), is due to drying and compaction of the low permeability clay layers, driven by over-exploitation of the underlying aquifers.

Local differential subsidence gradients threaten the integrity of the structure and infrastructure, whereas global subsidence produces tilts in the drainage network and in water reservoirs, and changes the flood patterns during the rainy season. Accurate characterization of Mexico City subsidence rate at different scales is therefore of great importance.

[Fig. 1 about here.] However, the application of multi-temporal interferometry in this area is very difficult. Coherence is lost due to vegetation cover, soil occupation changes and the relief of the volcanic chains outside the city make the application of multi-temporal interferometry very challenging. Previous InSAR studies of Mexico City subsidence include the computation of ERS, JERS, and ENVISAT interferograms [START_REF] Strozzi | Land subsidence in Mexico City mapped by ERS differential SAR interferometry[END_REF] [START_REF] Strozzi | Jers SAR interferometry for land subsidence monitoring[END_REF] [12] [START_REF] Cabral-Cano | Space geodetic imaging of rapid ground subsidence in Mexico City[END_REF], and a PS time series analysis based on 23 ENVISAT images from January 2004 to July 2006, focused on the western part of Mexico D.F. [START_REF] Osmanoglu | Mexico City subsidence observed with persistent scatterer InSAR[END_REF]. In this paper, we focus on the measurement of Mexico City subsidence using a data set of 38 ENVISAT images from November 2002 to March 2007. InSAR time series analysis is performed on this data set using both a modified SBAS approach dedicated to Mexico City subsidence specificity [START_REF] Cavalié | Ground motion measurement in the lake Mead area, Nevada, by differential synthetic aperture radar interferometry time series analysis: Probing the lithosphere rheological structure[END_REF] [16] and the PS approach of Gamma-IPTA [17] [18]. The intercomparison between these two approaches is realized in order firstly to validate subsidence maps and their given uncertainties, and secondly to benefit from their complementarity for characterization of subsidence rate at different scales. The available ground truth data of a few continuous GPS measurements [13] [14] and numerous levelling data spread across the basin do not reach the mm/yr accuracy necessary to validate InSAR time series analysis. In section II, we describe the application of the PS Gamma IPTA chain to the uncertainties associated with temporal unwrapping, particularly in the presence of nonlinear motion are presented in section III. In section IV, we describe briefly a modified SBAS approach [START_REF] López-Quiroz | Time series analysis of Mexico City subsidence constrained by radar interferometry[END_REF] to solve the specific problem of spatial unwrapping on Mexico City. In section V, we present a quantitative comparison between the PS and SBAS approaches and investigate the complementary information provided by both approaches to the characterization of Mexico City subsidence. We describe here the application of the Gamma-IPTA chain to the ENVISAT Mexico City data set (Figure 2). First, we produce a stack of SLC (Single Look Complex) radar images coregistered to the same master image, chosen in the middle of the temporal and spatial baseline space. PS candidates are then identified based on the temporal variability and the intensity of the backscattered echo. The mean intensity (ā) and the standard deviation (σ a ) of each pixel, and the overall mean intensity ( Ā) of all the pixels are calculated. A pixel is considered as a PS candidate if it satisfies simultaneously two following empirical criteria:

     D a = σa ā 0.58 ā 1.2 Ā (1)
where D a is the intensity dispersion index.

13.2×10 4 pixels are selected, the vast majority in the flat basin covered by the city and a few on human made structures on the lower parts of the volcano flanks surrounding the basin. The wrapped differential phase of PS candidates on the multi-temporal data stack is then computed using an interpolated SRTM DEM (Figure 1) [START_REF] Farr | Shuttle Radar Topography Mission produces a wealth of data[END_REF] and the DELFT precise orbit data [START_REF] Scharroo | Precise orbit determination and gravity field improvement for the ERS satellites[END_REF]. The unwrapped differential phase delay φ k for a point target in interferogram k can be expressed as:

φ k = 4π λ B k ⊥ Rsinθ h + 4π λ T k υ + φ k atmo + φ k orb + φ k nl + φ k noise ( 2
)
where the first term corresponds to the DEM error, h, the second term corresponds to the linear displacement rate, υ, φ k atmo is the atmospheric phase delay, φ k orb is the residual orbital error phase, φ k nl corresponds to the nonlinear ground motion and φ k noise is related to the decorrelation noise, which is assumed to be low on PS candidates. The phase of the topographic error varies linearly with the perpendicular baseline

B k
⊥ , with a proportionality coefficient depending on the wavelength of the radar carrier signal, λ, the August 24, 2012 DRAFT radar to ground distance, R, and the incidence angle, θ. The linear displacement phase is written as a function of the temporal baseline, T k . In order to extract the displacement time series from the wrapped differential phase, ϕ k , the PS software starts with a local unwrapping in the (T k , B k ⊥ ) space of the phase difference δϕ k between nearby PS candidates. This leads to increments in DEM error and displacement rate (δυ, δh), that are then integrated through space to get a global (υ, h) solution. The phase residuum, ϕ k res , is also obtained after spatial integration assuming that the local δϕ k res differences remain included in the [-π, π] interval. As the phase residuum includes the residual orbital trend, atmospheric delay, and nonlinear motion, their progressive filtering leads to a more accurate solution of displacement time series by iterations. The Gamma IPTA chain proposes seven iterative steps, but only the first three steps were successively applied to our data set. These are described below:

1) First temporal unwrapping. An initial estimation of the displacement rate (Figure 3 (a)) and height correction (Figure 3 (b)) is obtained on PS candidates through a 2D linear regression on the wrapped differential phase as a function of the perpendicular baseline and of the time interval of the temporal series. During the first step, only 34 interferograms whose perpendicular baseline is smaller than 800 m are used. Each interferogram is divided into 1 × 1 km 2 wide patches. The regression is performed on the phase of each PS in a given patch with respect to the patch reference phase, the phase difference at short distance including only small contributions from ϕ orb , ϕ atmo and possibly ϕ nl . The phase standard deviation of the regression, σ, is given as an assessment of the quality of the modeled (δυ, δh) and reflects the point decorrelation noise, plus other terms not taken into account in the regression. We keep 12.4 × 10 4 points with σ below the threshold of 1.1 radians. This threshold is chosen on the basis of synthetic tests described in section III. Beginning with the global reference point, the patch reference points are unwrapped by connecting each point to others in a propagating way. The local velocity and height difference, together with residual phase difference, are also propagated across the interferogram and referenced to the global reference point.

The velocity map presents the same pattern and amplitude as that derived from levelling [START_REF] Dgcoh | Fichas de nivelación de banco de nivel en la[END_REF].

[Fig. 3

about here.]

2) In the second step, the velocity and DEM error are re-estimated for the unwrapped phases of the 25 interferograms identified as being correctly unwrapped in the first step. Note that in this second regression without unwrapping, single patch processing is applied to avoid patch errors. About 8% of pixels with large phase standard deviation are masked.

3) In the third step, temporal unwrapping (as in step 1) is again performed only to refine the solution [START_REF]Interferometric point target analysis for deformation mapping[END_REF] DRAFT obtained in step 2, therefore allowing only a limited variation of 0.5 cm/yr for the velocity and of 3 m for the DEM error. For accuracy, 37 interferograms are used. After this step, 30 residual interferograms are qualitatively considered as correctly unwrapped based on the visualization of the residue.

The later steps proposed in the Gamma-IPTA chain have not been applied with success to our data set. In step 4, the application of spatial unwrapping does not appear to really improve the 7 residual interferograms, which are incorrectly unwrapped at step 3, although the visual appreciation of unwrapping quality on a sparse and noisy set of data points is quite subjective. The next steps, 5 and 6, include successive corrections of the residual orbit ramp and the atmospheric phase screen associated with changes in water vapour stratification. The first is done by optimizing the baselines such that the deviations between modeled and differential phases are minimal in the least squares sense. The second is obtained by linear regression between phase and elevation for each residual interferogram. However, there is only a very small number of PS outside the deformation area, and their locations are dissymetric: they are almost entirely located on the western side of Mexico City on the gentle slopes of the mountains. Estimation of residual ramp and stratified atmospheric contributions on stable ground is thus not possible here.

Furthermore, as the residual orbital contribution may be mapped onto the estimated velocity field and DEM error map (Figure 3), removing the velocity and DEM error cannot help the residual orbital ramp correction. As residual orbital ramps still remain in the phase time series, it appears inappropriate to apply the last correction step, which consists of smoothing out the turbulent atmospheric phase screen by a weighted temporal filter. Although the Gamma-IPTA chain could not be applied in its entirety, the subsidence map obtained will be compared and discussed with respect to that obtained using the SBAS processing.

The Gamma-IPTA chain was applied using different parameters before selecting the run that gives the best results, evaluated from both visualization of the velocity field, residual interferograms, and the comparison with the SBAS subsidence rate field (see below). The first problem encountered derived from aliasing of the phase values sampled at acquisition dates during the temporal unwrapping step. As temporal sampling is δt =35 days or a multiple of 35 days, the phase can be exactly equally adjusted with 4π λ v 0 × t or 4π λ (v 0 ± λ 2δt ) × t. There is thus an ambiguity in the retrieved velocity that amounts to 29.3 cm/yr. In order to reduce the impact of aliasing, we chose a global reference point in an area subsiding at a velocity of 20 cm/yr and limited the interval of adjustment for velocities in the first temporal unwrapping step to ±8 cm/yr. The second problem was to set the width of patches small enough to reduce the relative phase [START_REF]Interferometric point target analysis for deformation mapping[END_REF] DRAFT noise between reference points of neighbouring patches while limiting the number of patches to shorten the integration path. A patch width of 1 km appeared to yield good results. Finally, the threshold on the phase standard deviation was set to the optimum value in order to avoid erroneous estimated velocities (that could propagate through the reference point network), while keeping the maximum number of PS candidates. Synthetic tests (section III) show that a threshold at 1.1 radians is appropriate.

III. ANALYSIS OF THE TEMPORAL UNWRAPPING STEP USING SIMULATED DATA

A particularity of Mexico City subsidence is the very high subsidence rates, reaching 35 cm/yr in the line of sight (LOS) direction, and the high subsidence gradients across lithological boundaries between lacustrine deposits and volcanoes flanks. Furthermore, López-Quiroz et al. (2009) showed that subsidence is almost perfectly linear in time. These two characteristics make temporal unwrapping a priori more robust than spatial unwrapping to reconstruct the spatio-temporal subsidence evolution. However, temporal unwrapping may fail due to phase noise (decorrelation, atmospheric phase screen and residual orbital trend) and deviation from the assumed linear model. In this section, we use synthetic phase data sets to test the effect of nonlinear motion and noise. We construct synthetic time series of phase differences between two points δϕ k , wrap them, and then adjust the synthetic wrapped phases by a model, δΦ k , expressed as:

δΦ k = 4π λ (δυT k + δh B k ⊥ Rsinθ ) (3)
where δΦ k is the modeled phase difference between two points, δυ is the incremental ground velocity between two points and δh is the contribution associated with relative DEM error between two points.

Both δυ and δh are estimated by maximizing the norm of complex coherence γ:

γ = 1 N N k=1 e i(δϕ k -δΦ k ) ( 4
)
where N is the number of SAR data (N=38).

The simulation is performed for synthetic data sets without noise and with an added Gaussian noise.

In each case, we have three assumptions for the temporal deformation behaviour: (a) linear deformation (b) linear deformation with added acceleration (c) linear deformation with added periodic deformation.

The dates, T k , and the perpendicular baselines, B k ⊥ , are given by the 38 acquisitions of Mexico City. The wrapped phase differences are constructed from an input velocity (10 cm/yr) and an input DEM error (5 m). The tested values of the acceleration and the period of the periodic deformation are based on the results of the SBAS approach.

August 24, 2012 DRAFT (a) As expected, in the case of linear deformation and without added noise, the retrieved coherence peak reaches 1.0 at 10 cm/yr and 5 m. As soon as Gaussian noise is added and the noise level increases, the maximum coherence decreases and the coherence of secondary peaks increases. However, as long as σ remains lower than 1.1 radians, the retrieved velocity and DEM error remain at 10+/-0.1 cm/yr and 5+/-1 m with a few exceptions (Figures 4 (a) and (c)).

(b) In the case with an acceleration of 0.5 cm/yr 2 , the input velocity is the average velocity (10 cm/yr) during the acquisition period. Without noise, we observe that the highest coherence peak is wide and divided into three. The maximum coherence occurs at 9.6 cm/yr and 5.6 m with a coherence value equal to 0.64. Two secondary peaks are seated at 8.7 cm/yr, 4.5 m and 11 cm/yr, 3 m respectively.

When different Gaussian noise levels are taken into account, the main peak remains, with a few exceptions, at 9.6+/-0.1 cm/yr and 5.6+/-1 m as long as σ remains smaller than 1.1 radians (Figures 4 (b) and (d)). However, at the chosen velocity of 9.6 cm/yr, one unwrapping error always occurs during the temporal unwrapping procedure (see the example in Figure 5).

(c) In the case of inter-annual deformation (represented by a sinusoid of period 2 years), without noise, three coherence peaks (at 10 cm/yr, 8.5 cm/yr and 11.5 cm/yr) appear instead of one. As long as the sinusoid amplitude remains lower than 0.6 cm, the main peak remains at 10 cm/yr. However, when the sinusoid amplitude exceeds 0.6 cm, the coherence becomes maximum at 8.5 cm/yr with a value larger than 0.6. Applying a threshold on σ at 1.1 radians does not exclude all points with inter-annual deformation larger than 0.6 cm, resulting in unwrapping errors on points considered as reliable.

[Fig. 4 Although, according to Figure 4, PS candidates with σ value lower than 1.0 radians carry reliable phase information. However, we took 1.1 radians as the σ threshold for PS candidate elimination during the PS processing present in the previous section. This threshold was chosen on the basis of a compromise between the PS point quality and the spatial coverage.

To quantify the effect of the Gaussian noise level on the percentage of selected PS and the percentage of errors (retrieved wrong velocity and DEM error values), we perform statistics on 1000 synthetic data sets at each noise level without and with acceleration. Results are displayed in Figure 6. Without acceleration, as long as the noise level remains lower than 0.4 cm, all simulated PS candidates have a phase standard 6), we conclude that the chance of retrieving a wrong velocity value from selected PS should be low for the case of linear deformation. However, when acceleration is included, the number of selected PS candidates decreases drastically even for low noise levels (0.2 -0.5 cm). For noise levels of 0.4 cm or above, 3% of the PS candidates have both a phase standard deviation lower than 1.1 radians and a velocity value outside the 9.3 ∼ 10 cm/yr range (9.6 cm/yr is the velocity value retrieved without added noise). Out of these, most have σ value larger than 1 radians. For the case of nonlinear motion, velocity errors are not randomly located but present clusters, as shown in Figures 4 (a) and (b).

We conclude from the synthetic tests that the probability of obtaining erroneous velocities on selected PS candidates is very low. However, a large proportion of PS targets presenting relative nonlinear motion (of amplitude above 0.6 cm) will be discarded based on their σ values obtained during the temporal unwrapping step. Some may be kept, but present one unwrapping error with a high probability.

[Fig. 6 about here.]

IV. SBAS PROCESSING

The ENVISAT Mexico City data set was also processed with a modified SBAS approach dedicated to Mexico City subsidence specificity. This approach, designed specifically to solve the unwrapping problem in areas of large and relatively stable deformation, is explained in detail in [START_REF] López-Quiroz | Time series analysis of Mexico City subsidence constrained by radar interferometry[END_REF]. Here, we just outline the main characteristics. Firstly, 72 differential interferograms are constructed with perpendicular baselines less than 500 m and temporal baseline less than 9 months. These small baselines allow a maximization of the spatial coherence and a reduction of the signal associated with DEM errors. Filtering, unwrapping, correction of stratified atmospheric delays and residual orbital trend, and inversion of interferograms are performed in three successive iterations.

In the first iteration, the raw differential interferograms (with a multi-looking of 1:5 in range and azimuth direction respectively) are adaptively filtered using the Goldstein filter that damps the spectrum in the Fourier domain [START_REF] Goldstein | Radar interferogram filtering for geophysical applications[END_REF] and are unwrapped spatially by the branch-cut algorithm of ROI-PAC [START_REF] Rosen | Updated Repeat Orbit Interferometry Package released[END_REF].

Some could not be unwrapped correctly due to high deformation fringe rates. Nevertheless, the stratified August 24, 2012 DRAFT atmospheric delay and the residual orbital contribution are jointly estimated by linear adjustment to the differential phase outside the deformation area, i.e., outside the flat portion of the basin. Then, a SVD inversion allows an examination of the closure of the redundant interferometric network and a quantification and identification of unwrapping errors.

The next two iterations repeat the previously described procedure, but with a "guide" to the unwrapping step. [START_REF] López-Quiroz | Time series analysis of Mexico City subsidence constrained by radar interferometry[END_REF] pick the 5 best interferograms with high signal to noise ratio and no phase unwrapping error, and stack them to represent an average deformation rate (the so called deformation model hereafter).

An adaptative filter is applied to the stack to reduce noise. This deformation "model" is then scaled by least squares adjustment to each interferogram unwrapped in the previous iteration. Residual orbital ramp and stratified atmospheric delays are also estimated from the unwrapped interferograms obtained in the previous iteration. All these terms are removed from the raw differential wrapped interferograms. The resulting residual interferograms present a limited number of fringes, including turbulent atmospheric patterns, the deformation that does not follow the "model", and noise. After spatial filtering, they are unwrapped by SNAPHU [START_REF] Chen | Phase unwrapping for large SAR interferograms: Statistical segmentation and generalized network models[END_REF]. A new inversion allows an examination of closure errors and a quantification of unwrapping errors. The latter decreases significantly. This unwrapping "guide" is repeated twice to refine the "model" scaling, the estimation of stratified atmospheric delay and the residual orbital contribution in the first step. A further mask is applied to interferograms in which the residual phase exceeds 4 radians. [START_REF] López-Quiroz | Time series analysis of Mexico City subsidence constrained by radar interferometry[END_REF] verify that the interferometric system mis-closure drops at each iteration.

Finally, a set of 71 unwrapped interferograms is obtained, successfully corrected from stratified atmospheric delay and orbital residual contribution, and with a phase reference equal to zero outside the flat basin, masked in noisy areas. After inversion, the phase delay time series in general show a remarkably linear subsidence through time. However, in some areas, a non negligible acceleration or deceleration occurs. Furthermore, for some pixels, the matrix for inversion has a rank deficiency, i.e. at least one critical link in the interferometric network is missing. In these cases, the acquisition data set is split into two or more independent image groups. Since a non negligible acceleration or deceleration is observed in some areas, one additional constraint, stating that the phase varies as a quadratic polynomial in time and linearly with perpendicular baseline, is then added to the design matrix with a sufficiently small weight so that it only fixes the offsets between phase delay time series of independent image groups.

The average ground motion rate is shown in Figure 7 (a), and a zoom on the root mean square (RMS) mis-closure maps of the inconsistencies ϕ RM S in the interferometric network is shown in Figure 7 (b). The uncertainty of the subsidence velocity is equal to 0.7 mm/yr, for pixels without unwrapping error (i.e. for ϕ RM S values lower than 0.35 radians) and for which the interferometric set is complete.

August 24, 2012 DRAFT This standard deviation value is computed from the linear regression of phase versus time, assuming an independent Gaussian distribution for phase errors. Slightly larger uncertainties are expected if fewer images or interferograms are available for a given point. Note also, that due to the slight adaptative filtering applied to 1×5 looks interferograms, the solution can be considered as regularized in space, and thus cannot quantitatively provide the subsidence velocity of an individual target that could subside differently from its neighbors.

[Fig. 7 about here.]

V. JOINT ANALYSIS OF PS AND SBAS APPROACHES

The PS and SBAS approaches are applied to the same data set, which allows a detailed and quantitative comparison between subsidence velocity maps built from both approaches. The first aim is to validate the results obtained by each of the two approaches and to discuss whether the uncertainty claimed is compatible with their inter-comparison. The second aim is to highlight their complementarity in measuring

Mexico City subsidence at different scales.

The displacement rate estimated by the first three steps of the PS approach is compared to the average velocity obtained by the SBAS approach. Note that, because of inaccuracies in the relative positioning of PS and SBAS velocity maps in ground geometry, we need coregister them by applying a subpixel shift that maximizes their cross correlation. This step is important as the lateral variation in subsidence velocity is extremely large in Mexico City. The SBAS velocity at PS locations is calculated by bilinear interpolation using the available grid points around the PS point (4 grid points in most cases, less in a few cases). Finally, to analyze in detail the dispersion of the PS-SBAS results, we use three additional parameters, the PS phase standard deviation, σ, the map of nonlinear deformation derived from the SBAS analysis, ϕ nl , and the SBAS RMS mis-closure, ϕ RM S . Out of the 12.3 × 10 4 PSs, for which we also have a SBAS value at PS step 1, 3.3% have a σ between 1.0 and 1.1 radians (selected PS with possible unreliable phase value), 6.6% are in areas with nonlinear deformation greater than 1.3 radians, and 1.1% have a ϕ RM S greater than 0.35 radians (with a possible unwrapping error from SBAS). Unreliable results might thus be expected in these three cases.

A. Global comparison

The PS velocity map at step 1 and the average SBAS velocity map at the same scale are shown respectively in Figure 3 (a) and in Figure 7 (a). The velocity map obtained by the SBAS analysis extends August 24, 2012 DRAFT beyond the limit of urban areas and of the valley, wider than the area covered by the PS approach. The same velocity pattern is observed, emphasizing the general agreement between both approaches. A strong difference appears only on a small number of PS patches at PS step 1. To quantify the difference, the velocity estimated at the PS step 1 is plotted in Figure 8 Let us now analyze quantitatively the distribution of the remaining points presenting large differences, at step 3, knowing that large unwrapping errors have been eliminated. Apart from the aliasing problem, only 143 points (∽ 0.14%) present differences larger than 2.5 cm/yr. These may represent isolated points unseen by the SBAS approach due to spatial regularization, whose subsidence rate strongly differs from the neighbouring pixels. Some of these points may be explained by large σ values (10% for σ larger than 1.0 radian) and/or large ϕ nl (24% with ϕ nl > 1.3 radians) and/or large ϕ RM S (19% with ϕ RM S > 0.35 radians), these proportions being clearly larger than in the 'normal' PS population. Finally, only 85 points (less than 1/1000) remain, for which one may suspect that they carry specific subsidence information as isolated points. Most of these points have larger PS subsidence velocity than that measured by the SBAS approach, which could represent isolated buildings subsiding faster due to their weight. Possible strong differential uplift with respect to the surrounding pavement, as could affect isolated buildings rigidly anchored by piles in the hard layer ("Capa dura"), does not appear in the PS-SBAS difference maps, although SBAS interferograms were spatially filtered.

[Fig. 8 about here.]

Now we focus on the 'normal' PS population which displays moderate PS-SBAS differences, lower than 2 cm/yr. After the removal of a velocity ramp due to residual orbital trends in the PS interferograms (see next subsection), the distributions of PS-SBAS velocity difference are similar at PS steps 1, 2 and 3, August 24, 2012 DRAFT with a standard deviation of 1.7, 1.6 and 1.6 mm/yr respectively. Similar cross-comparison values have been obtained in a different area by [START_REF] Shanker | Comparison of Persistent Scatterers and Small Baseline Time-Series InSAR Results: A Case Study of the San Francisco Bay Area[END_REF]. To further analyze the moderate PS-SBAS differences, it is interesting to map them both with some spatial smoothing (low-pass maps) and at a local scale (highpass maps). The PS approach can retrieve punctual deformation with accuracy, but may not propagate the solution well spatially, whereas the SBAS approach, due to slight adaptative spatial filtering, may not recover punctual deformation. The difference maps are analyzed in the next subsections after the application of a low-pass and a high-pass filter.

B. Comparison of low-pass subsidence maps

The low-pass (LP) filtered difference maps are obtained by averaging the PS-SBAS velocity difference in 450 m wide sliding windows (Figure 9). Distributions of LP difference maps are provided in Figure 10.

For PS step 1, the LP difference map (Figure 9 and 3, the spatial integration of temporal unwrapping results thus seems to produce small errors of the order of 2 mm/yr, mainly outside the city. Features common to the three PS-SBAS difference maps can be considered robust. In particular, we note that the computed average subsidence velocity in the main subsident area is about 1 to 2 mm/yr lower for the PS approach than for the SBAS approach. However, the filtering used here to compute the LP maps is based on a simple moving average that may produce offsets if isolated points present large SBAS-PS differences.

[Fig. 9

C. Comparison of high-pass subsidence maps

High-pass (HP) SBAS-PS difference maps are obtained by removing the LP maps described above simply from the raw difference maps. The HP difference map, displayed for PS step 3 in Figure 11,

shows that, in stable areas with null or low subsidence rates, the variability is small. By contrast, in the areas where the subsidence is large, as in Mexico City centre, the variability is significantly larger. The histograms of HP difference maps are represented in Figure 10 for the first three PS steps. They show very similar distributions for the three steps, with standard deviations (σ HP ) of 1.1 mm/yr. The standard deviation decreases further (0.8 mm/yr) within areas with little subsidence (less than 2 cm/yr), well within the uncertainties given by both the PS and SBAS approaches, valid if no unwrapping error occurs. On the contrary, in areas with subsidence rates larger than 2 cm/yr, the standard deviation reaches 2.4 mm/yr (Figure 10 (d)). These statistics were performed excluding points with σ > 1.0 radian, ϕ nl > 1.3 radians and ϕ RM S > 0.35 radians. Therefore, we can explain part of this variability by a slightly heterogeneous deformation behaviour depending on the characteristic of roads and buildings. As mentioned above, the SBAS approach provides a spatially filtered measurement that does not take into account the possible point-like displacement and building height of individual targets.

In conclusion, we find that the amplitude of HP PS-SBAS deviations is surprinsingly small, given the urban nature of the area of interest (elevated railways or roads, skyscrapers, etc...), the variability of underground basement (none, pilars, floating shafts, ...), and the variability of the time elapsed since construction (Aztec pre-consolidation, hispanic times, recent constructions, none). The good agreement between both approaches gives us confidence in the interpretation of local features using not only the PS measurements, but also the SBAS motion map.

[Fig. 11 about here.]

D. Effect of Nonlinear motion

In the literature, nonlinear deformation is retrieved using the PS approach by applying specific filters [START_REF] Ferretti | Nonlinear Subsidence Rate Estimation Using Permanent Scatterers in Differential SAR Interferometry[END_REF] [26] [START_REF] Liu | Estimating Spatiotemporal Ground Deformation With Improved Persistent Scatterer Radar Interferometry[END_REF]. However, our simulated tests of the PS unwrapping step, presented in section III, indicate that, due to the assumption of a linear deformation model, the nonlinear contribution will increase the apparent noise estimation, the σ value, thus increasing the chances of the pixel being rejected (Figure 6).

Furthermore, in the case of non negligeable nonlinear deformation (i.e., larger than ± 0.6 cm), either in form of interannual fluctuations or of acceleration/deceleration, the retrieved velocity differs from the prescribed average velocity during the study period, and at least one unwrapping error occurs (Figure 5), August 24, 2012 DRAFT with the effect of linearizing the phase delay evolution with time.

Note that the assumption of all PS approaches is that nonlinear motion between neighbouring points is small enough to limit these effects. The success of retrieving nonlinear motion should thus depend on its lateral gradients and on the distance between PS candidates and their reference point. To test the ability of the present PS study, with 1 km wide patches, to recover nonlinear motion, we compare the amplitude of nonlinear motion derived using the SBAS approach, ϕ nl and the PS phase standard deviation, σ, obtained during the first temporal unwrapping step. A zoom on maps of σ and ϕ nl (Figure 12) shows that the phase standard deviation (σ) is large where the nonlinear residue (ϕ nl ) is also large. A few PS patches also present large phase standard deviation but with little nonlinear motion: this might be due to phase noise on the patch reference points. Statistics on σ values as a function of ϕ nl are displayed in

E. Combination of PS and SBAS results

The spatial regularization performed in the SBAS approach results in a far broader coverage of the Mexico Valley and of the volcanoes flanks, outside the main urban areas, than that obtained with the PS approach. It allows residual orbital trends and stratified atmospheric effect to be separated from the deformation signal and thus provides a relatively stable reference against which subsidence is measured.

Even within the urban areas, the target density from which subsidence information is retrieved is far larger with the SBAS approach than using the PS approach. This provides a continuous description of the subsidence field. The PS approach could not be applied everywhere due to the limited PS point density, especially outside the urban areas. The failure of the correction of residual orbital trends and stratified atmospheric effect limits the accuracy of the PS derived subsidence field at large scales and of the referencing to non deforming areas. However, in urban areas, the PS approach captures the punctual height and deformation of single object, relevant for monitoring possible differential motion between August 24, 2012 DRAFT structures and pavements at a local scale, while taking into account the point target height contribution to the phase. We thus consider local measurements of subsidence by the PS approach as particularly accurate.

In order to combine the PS and SBAS results, we assume that the continuous subsidence field obtained with the SBAS approach can serve as a reference for the general subsidence features against which relative PS local subsidence rates can be studied. The LP SBAS filtered subsidence map is computed by adjusting quadratic polynomial surfaces to sliding windows, iteratively removing outliers. We thus expect this surface to adjust to the most populated smoothly varying subsidence pattern. The HP residual SBAS subsidence map is displayed in Figure 14 for the NW and center parts of the studied area. High HP velocities contour the two volcanoes outcropping in the city, the Peñon de los Banos (close to the airport) and the Peñon de los Marques (close to cierro de la Estrella). They are also found preferentially in areas with clear nonlinear motion. The large amplitude of HP features in the vegetated areas in the East of the Figure is not well constrained due to strong decorrelation noise. Besides these features, the HP SBAS subsidence map shows interesting signals, some being clearly aligned with the roads and block geometry, while others appear to reflect single structures.

[Fig. 14 about here.]

The HP PS subsidence map is obtained after correction of the residual orbital ramp and removal of the LP SBAS map. The discontinuous nature of PS measurement does not allow linear features apparent in the HP SBAS map to be distinguished. We display in Figure 15 a comparison between SBAS and PS HP subsidence rates in a small area centered on a linear feature observed in the SBAS HP map. Other PS results in the same area were previously described by [START_REF] Osmanoglu | Mexico City subsidence observed with persistent scatterer InSAR[END_REF]. The linear feature appears to follow an aerial portion on subway line 4 between stations Morelos and Talisman. Large relative uplift rates characterize the subway stations Morelos (1.7 and 2.1 cm/yr for SBAS and PS, respectively) and Consulado (2.0 and 2.1 cm/yr for SBAS and PS, respectively), but smaller uplift rates are also apparent at Canal del Norte, Bondojito, and Talisman stations. The agreement between PS and SBAS maps is striking, but not perfect. Isolated velocity anomalies appear correctly retrieved by the SBAS approach, but attenuated.

The PS sampling seems biased towards structures presenting relative uplift or subsidence, with an underrepresentation of the pixels with 'background' velocity values: The latter may have on average darker backscattering properties. To conclude, the PS-SBAS HP comparison allows a qualitative validation of the HP SBAS map, which is useful for interpretation of local subsidence pattern, while local PS velocities are better constrained and localized, but subject to sparse and uneven sampling. In this paper, the results of the two main multi-temporal InSAR processing approaches are obtained to measure Mexico City subsidence, characterized by very high subsidence rates and large deformation gradients. Moreover, their common features and discrepancies are analyzed and discussed, in order to validate the results from both approaches, verify their claimed nominal uncertainties and investigate the complementarity of both approaches for characterization of subsidence at different scales.

Our first observation is the very good agreement found between the used PS approach and the SBAS approach, with differences statistically close to the nominal mm/yr uncertainty. The LP difference maps between the SBAS and PS velocity fields show an agreement reaching 1.2 mm/yr. The HP filtered difference maps also yield an excellent agreement outside the subsident area reaching 0.8 mm/yr. More surprisingly, within subsiding areas, the local differences amount to 2.4 mm/yr, a value interestingly small given that part of this difference may come from coregistration misfits and approximations in the interpolation of SBAS values on PS points. The validation provided here is credible because both approaches adopt completely different data processing strategies: their possible error sources, that are not taken into account in the nominal uncertainty calculation, can be considered as independent.

The claimed norminal uncertainties of both approaches, in the order of mm/yr, have been verified

through the inter-comparison in this application. This nominal uncertainty is given from the standard deviation of the phase versus time, once unwrapped, after applying various corrections, and assuming independent random phase noise around the assumed linear trend. In this study, we identify the main error source as coming from the unwrapping step. Because of large subsidence gradients and temporal decorrelation, spatial unwrapping of small baseline interferograms is particularly difficult in the SBAS approach. The problem was solved by using a regional subsidence field to 'guide' unwrapping and by applying a LP filter to wrapped interferograms. However, unwrapping may still fail, if and where the regional subsidence field is not accurate enough, or if and where there exists large local differential subsidence. On the other hand, the PS temporal uwrapping strategy is a priori well suited to study Mexico City subsidence, which evolves mostly linearly in time. However, it is hindered by the large spatial extent of the subsident area, the small but non negligible nonlinear motion in certain parts of the city, and by the low PS density outside the urban area.

Both approaches appear to complement each other: the SBAS approach gives a continuous description of the subsidence rate patterns and temporal behaviour without linearity assumption, while the PS August 24, 2012 DRAFT approach allows a quantitative discussion of the subsidence disparity at small scales between various man-made structures. The combination of PS and SBAS subsidence measurements shows that, despite the very large subsidence rate (up to 38 cm/yr) affecting the Mexico City basin, the local differential motions are quite limited, reaching a few cm/yr of excess subsidence or uplift, on individual targets or along linear structures (e.g., roads, subways, tramways). A detailed comparison between differential subsidence and both the history of settlement (from pre-hispanic to recent constructions) and the type of foundations should be useful to quantify where and why compaction occurs. The limited amplitude of local differential motions suggests that the main factor controlling compaction is the deep consolidation by progressive depressurization of the clay layers above the aquifer, surface consolidation associated with the urban load being only secondary. Consequently, the combination of these two approaches [START_REF] Hooper | A multi-temporal InSAR method incorporating both persistent scatterer and small baseline appoaches[END_REF] [START_REF] Liu | Estimating Spatiotemporal Ground Deformation With Improved Persistent Scatterer Radar Interferometry[END_REF] appears promising in multi-temporal InSAR processing, from which a more detailed measurement and a wide range of applications can be expected in the future.

ACKNOWLEDGMENT

This work was supported by the EFIDIR project (http://www.efidir.fr) granted by the French National Agency (ANR) (ANR-07-MDCO-004). The authors wish to thank the ANR for their support. We would like to thank the European Space Agency (ESA) for supplying the ENVISAT images through the category -0.

  Finally, conclusions and open problems are addressed in section VI. II. PS PROCESSING [Fig. 2 about here.]

  about here.] [Fig. 5 about here.]

  (a) as a function of the SBAS average velocity. This plot first shows the very good agreement, with points aligned on the line y = x. However, points located away from this line with parallel trends mark large unwrapping errors at PS step 1. A few points with a PS-SBAS difference of about ± 29 cm/yr correspond to the aliasing problem. The dispersion of the PS-SBAS velocity difference is not correlated with the PS phase standard deviation, σ, (Figure8 (b)), because all points with σ above 1.1 radians have been discarded (as shown in section III). The same plots are shown in Figures8 (c) and (d) for PS step 2, for which the velocity is estimated by linear regression using only well unwrapped interferograms. A few points with large PS-SBAS differences have σ values larger than 1 radian. All linear trends parallel to the line y = x have disappeared, except for the aliasing trend. This step is therefore very efficient to mask areas with large unwrapping errors. Aliasing errors cannot be detected and remain in all PS steps.

  (a)) displays a NW-SE ramp of ±5 mm/yr, which arises from residual orbital trends in the interferograms used in the PS approach. The residual ramp is estimated and removed from the LP velocity difference map. The detrended LP map is shown in Figure 9 (b) (σ LP = 1.2 mm/yr, Figure 10 (a)), with a narrower color scale that emphasizes the differences. A small PS unwrapping error affects the NW part of the map, bounded by patch boundaries, in the form of a staircase. At PS step 2 (Figure 9 (c)), the LP map after ramp correction appears only marginally better constrained and without patch errors (σ LP = 1.1 mm/yr, Figure 10 (b)). At PS step 3, temporal unwrapping reintroduces a few patch errors, visible in the LP map in the NW corner in form of a small SE-NW trending band bounded by staircase limits (Figure 9 (d), σ LP = 1.1 mm/yr, Figure 10 (c)). At steps 1

Figure 13 .

 13 Figure 13. At low ϕ nl values, the phase noise peaks at 0.55 radians, which must thus be typical of the phase noise for PS candidates selected on the basis of their low D a values and without the contribution of nonlinear motion. As soon as ϕ nl reaches 1.5 radians, one obtains statistically larger σ values, around 0.8 radians. For a nonlinear motion with an amplitude larger than ∽ 3 radians (i.e., ∽ 1.4 cm), most PS σ values are above 1.1 radians. We conclude that the present PS analysis fails to retrieve the nonlinear subsidence temporal behaviour in the Mexico Basin, where it is the strongest.

Fig. 1 .Fig. 2 .Fig. 3 .

 123 Fig. 1. Shaded SRTM elevation relief around Mexico City Valley. The ENVISAT frame is contoured by the white rectangle.The limits of the geotechnical zones are displayed in black (Zone I: bedrock, Zone II: Alluvial deposits, Zone III: lacustrine deposits)[START_REF]Mexico Building Code. Normas Tècnica complementarias para diseño y construcción de cimentaciones Federal District Government[END_REF].

Fig. 4 .Fig. 5 .

 45 Fig. 4. (a) (b): Distribution of retrieved DEM error and velocity that maximize coherence without and with acceleration (0.5 cm/yr 2 ). (c) (d): Residual phase standard deviation as a function of the retrieved velocity without and with acceleration. Different colours represent different Gaussian noise levels added to the synthetic phase series. Each point corresponds to one of the 1000 random tests.

Fig. 6 .

 6 Fig. 6. (a) Percentage of PS candidates with a phase standard deviation, σ, below the threshold of 1.1 radians as a function of the Gaussian noise level. (b) Percentage of PS candidates presenting both a phase standard deviation, σ, below the threshold of 1.1 radians and a wrong velocity value, as a function of Gaussian noise level. Synthetic tests (1000 random tests per noise level) are performed without (solid lines) and with acceleration (dashed lines). The filled triangle in panel (a) corresponds to the noise distribution of the PS candidates in areas without nonlinear motion.

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

 789 Fig. 7. (a) Average LOS ground displacement rate estimated by the SBAS approach. A colour cycle represents 15 cm/yr. To be compared with Figure 3 (a). (b) Zoom on the RMS mis-closure map of the SBAS approach (ϕRSM ), in the area located by the black rectangle in (a). Some unwrapping errors may have occured in orange areas with ϕRMS > 0.35 radians.

Fig. 10 .Fig. 11 .Fig. 12 .

 101112 Fig. 10. (a, b, c) : Normalized histograms of PS-SBAS differences of HP (continuous lines) or LP (dashed lines) filtered maps, at PS steps 1, 2, 3 respectively. (d) Normalized histogram at PS step 3 of HP filtered difference map, separating points located outside (with v > 2 cm/yr) or inside (with v < 2 cm/yr) the main subsiding area.

Fig. 13 .

 13 Fig.[START_REF] Cabral-Cano | Space geodetic imaging of rapid ground subsidence in Mexico City[END_REF]. PS phase standard deviation (σ) as a function of the amplitude of the nonlinear motion (ϕ nl ) computed by the SBAS approach. Each PS is shown by a dot. The background color map represents the PS density with a log10 scale and is evaluated by Gaussian statistics.

Fig. 14 .Fig. 15 .

 1415 Fig.14. HP filtered subsidence velocity map obtained using the SBAS approach, surperimposed on the radar amplitude map. The figure is drawn in radar geometry. The black rectangle delimits the area depicted in Figure12. Note that the color scale saturates towards large positive (excess subsidence) or negative (relative uplift) velocities.

  lower than 1.1 radians and the retrieved velocity is always correct. When the noise level reaches 0.7 cm and above, only few PS candidates have σ value lower than 1.1 radians, but 80% of them have wrong retrieved velocity values. The percentage of errors (PS candidates that are selected but present wrong velocity values) is about 3% at a noise level of 0.6 cm or above. Because PS candidates have statistically low D a values and thus low phase noise (see gray filled triangle in Figure

August 24, 2012 DRAFT deviation, σ,

  1 project no. 3979. Some figures were prepared with the GMT software by P. Wessel and W. H. F. Smith. Shaded SRTM elevation relief around Mexico City Valley. The ENVISAT frame is contoured by the white rectangle. The limits of the geotechnical zones are displayed in black (Zone I: bedrock, Zone II: Alluvial deposits, Zone III: lacustrine deposits)[START_REF]Mexico Building Code. Normas Tècnica complementarias para diseño y construcción de cimentaciones Federal District Government[END_REF]. . . . . . . . . . . . 23 2 Perpendicular baseline of the 37 computed interferograms relative to the master date (2004-12-31), versus acquisition time. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 3 (a) LOS linear deformation velocity (a colour cycle represents 15 cm/yr) and (b) Elevation correction (a colour cycle represents 20 m) estimated at PS step 1. They are superimposed on the radar backscatter amplitude map. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 4 (a) (b): Distribution of retrieved DEM error and velocity that maximize coherence without and with acceleration (0.5 cm/yr 2 ). (c) (d): Residual phase standard deviation as a function of the retrieved velocity without and with acceleration. Different colours represent different Gaussian noise levels added to the synthetic phase series. Each point corresponds to one of the 1000 random tests. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 5 (a) Example of a constructed phase time series before (in black) and reconstructed after (in red) temporal unwrapping. The best coherence (γ = 0.58) during unwrapping is obtained with a velocity of 9.54 cm/yr. Choosing this velocity produces an unwrapping error around May, 2006. (b) Residual obtained after linear adjustment to the input phase time series (in black) or to the unwrapped output phase time series (in red). The synthetic phase time series is built with a velocity of 10 cm/yr, an acceleration of 0.5 cm/yr 2 and a Gaussian noise of 0.3 cm. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 6 (a) Percentage of PS candidates with a phase standard deviation, σ, below the threshold of 1.1 radians as a function of the Gaussian noise level. (b) Percentage of PS candidates presenting both a phase standard deviation, σ, below the threshold of 1.1 radians and a wrong velocity value, as a function of Gaussian noise level. Synthetic tests (1000 random tests per noise level) are performed without (solid lines) and with acceleration (dashed lines). The filled triangle in panel (a) corresponds to the noise distribution of the PS candidates in areas without nonlinear motion. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 7 (a) Average LOS ground displacement rate estimated by the SBAS approach. A colour cycle represents 15 cm/yr. To be compared with Figure 3 (a). (b) Zoom on the RMS mis-closure map of the SBAS approach (ϕ RSM ), in the area located by the black rectangle in (a). Some unwrapping errors may have occured in orange areas with ϕ RM S > 0.35 radians. . . . . . . 29 8 (a) (c) Velocity estimated by the PS approach as a function of SBAS average velocity. (b) (d) Difference between velocities estimated by PS and SBAS approaches as a function of PS phase standard deviation. (a) (b) at PS step 1. (c) (d) at PS step 2. . . . . . . . . . . . . 30 9 LP difference maps between PS and SBAS velocity field. (a) at PS step 1, σ LP =3.3 mm/yr. (b) at PS step 1, with removal of the residual orbit ramp, σ LP =1.2 mm/yr. (c) at PS step 2, with removal of the residual orbit ramp, σ LP =1.1 mm/yr. (d) at PS step 3, with removal of the residual orbit ramp, σ LP =1.1 mm/yr. A colour cycle represents 1 cm/yr in (a) and 0.5 cm/yr in (b-d). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Zoom on the PS phase standard deviation defined during temporal unwrapping at PS step 1. (b) Zoom on nonlinear deformation detected by the SBAS approach. (a) and (b) correspond to the same area. The areas with the highest nonlinear deformation (red to yellow in (b)) correspond to areas with high σ values (in green in (a)). . . . . . . . . . . . . . . . . . . . 34 13 PS phase standard deviation (σ) as a function of the amplitude of the nonlinear motion (ϕ nl ) computed by the SBAS approach. Each PS is shown by a dot. The background color map represents the PS density with a log10 scale and is evaluated by Gaussian statistics. . . . . 35 14 HP filtered subsidence velocity map obtained using the SBAS approach, surperimposed on the radar amplitude map. The figure is drawn in radar geometry. The black rectangle delimits the area depicted in Figure 12. Note that the color scale saturates towards large positive (excess subsidence) or negative (relative uplift) velocities. . . . . . . . . . . . . . . 36 15 Comparison of HP filtered subsidence velocity maps obtained using the SBAS (a) and PS (b) approaches, surperimposed on the radar backscatter amplitude map. The subway station locations are slightly offset towards the West in order not to mask the velocity field. This figure can be compared to Figure 11 in [14]. . . . . . . . . . . . . . . . . . . . . . . . . . . 37

	August 24, 2012 1 10 (a, b, c) : Normalized histograms of PS-SBAS differences of HP (continuous lines) or LP LIST OF FIGURES . . 31 DRAFT (dashed lines) filtered maps, at PS steps 1, 2, 3 respectively. (d) Normalized histogram at PS step 3 of HP filtered difference map, separating points located outside (with v > 2 cm/yr) or inside (with v < 2 cm/yr) the main subsiding area. . . . . . . . . . . . . . . . . . . . . . 32 11 HP difference map August 24, 2012 DRAFT 12 (a) August 24, 2012 DRAFT

between PS and SBAS velocity field at PS step 3, σ HP =1.1 mm/yr. A colour cycle represents 1 cm/yr. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33

August 24, 2012 DRAFT


Guillaume Pizaine

Raphaël Prevost

Elsa D Angelini

Isabelle Bloch

Sherif Makram-Ebeid

Raphael Prevost

Segmentation-free and multiscale-free extraction of medial information using gradient vector flow -Application to vascular structures

Keywords: gradient vector flow, diffusion, medialness, skeleton, shape analysis

come

INTRODUCTION

Gradient Vector Flow (GVF) has first been introduced as an external force field for active contours and active surfaces in [START_REF] Xu | Snakes, Shapes, and Gradient Vector Flow[END_REF] [START_REF] Xu | Snakes, Shapes, and Gradient Vector Flow[END_REF]. The GVF of an image is the vector field obtained by diffusing image gradients in homogeneous regions while keeping strong gradients untouched. The diffusion process spreads edge information into uniform regions and acts as a long range force (see Fig. 1). Consequently, it also introduces more robustness against initialization and speeds up convergence.

Formally, the GVF of an image I over a domain Ω is defined as the global minimizer V (Xu et al., 2000 [START_REF] Xu | Global Optimality of Gradient Vector Flow[END_REF]) of the following energy functional E:

E = Ω g(x) ∇V 2 (x) + h(x)|V (x) -∇I(x)| 2 dx , (1)
where g : Ω → R and h : Ω → R are spatially-varying weighting functions and ∇V is the vector norm for tensors given by √ ∇V .∇V . The first term is a regularization term that controls the diffusion over the whole image domain. The second term is a data attachment term which ensures that V is close to the image gradient at strong edges. This is the General Gradient Vector Field (GGVF) devised by [START_REF] Xu | Generalized Gradient Vector Flow External Forces for Active Contours[END_REF] [START_REF] Xu | Generalized Gradient Vector Flow External Forces for Active Contours[END_REF], which comes down to the original formulation of the GVF [START_REF] Xu | Snakes, Shapes, and Gradient Vector Flow[END_REF] [START_REF] Xu | Snakes, Shapes, and Gradient Vector Flow[END_REF]) if g is constant and h(x) = |∇I(x)| 2 . The most widely used functions are g(x) = e -|∇I(x)| 2 /K 2 , K ∈ R and h(x) = 1 -g(x), and will be used in this paper too. Since both formulations yield similar results, we will use the term GVF for both in the remaining of the paper.

The first variation of the functional E yields the following Euler-Lagrange equation1 :

g(x)∆vi(x) -h(x)(vi(x) -∇I(x)) = 0 , ( 2
)
where vi is the i-th component of the vector field and ∆ is the Laplacian operator. The GVF is then the steady state of Eq. 2.

Recently, GVF has become popular in the field of medial information extraction. Many ways of using it have been proposed since it can be viewed as an improved gradient vector field to compute various features. For instance, Bauer et al., 2009 [4] propose to recover the centerlines of airways by computing the Hessian matrix from the GVF. Then, they determine the cross-sectional planes of the tubular structures and compute a tube-likeliness map from flux measures in those planes, based on the GVF, again. Flux measures were also used in Engel et al., 2008 [5] for medial features detection. Previous works also exhibit GVF-based medialness map derived from observations. Among them, the tube-likeliness from Bauer et al., 2009 [4] has already been mentioned. In Yu et al., 2004 [6], the authors propose to build a skeleton strength map from the GVF norm for grayscale image segmentation. Finally, the GVF has also been used to extract skeletons from binary shapes. In this context, the GVF is used in Hassouna et al., 2007 [7] in a front propagation setting to design a speed function allowing faster propagation at the center of structures.

Although the GVF has already been used to extract medial information, few works have proposed approaches to recover scale information. Unlike multiscale filters, which retain the maximum response over several scales, the GVF diffuses information without keeping track of the scale. Although one benefits from this by freeing oneself from scale constraints (e.g. Hessian matrices can be computed on a 3x3 neighborhood only), scale information is still paramount for skeletons or medialness maps. Knowing the centerlines, the method in Bauer et al., 2009 [4] goes back to the airways wall by tracking the GVF back to the edges in the image, which is quite time-consuming. In Engel et al., 2008 [5], the authors recover the size of the structures as the radius yielding a maximal circular (or spherical) flux. It seems in contradiction to the multiscale-free approach of the GVF.

In this paper, we propose a simple, segmentation-free and multiscale-free algorithm to extract medial information from images, based on the GVF. Since our approach heavily relies on one-dimensional analysis of the GVF (line by line in different directions), Sect. 2 gives a thorough review of the analytic solution to the one-dimensional case. Section 3 details the algorithm, especially how scale information is recovered. Finally, we discuss parameters and show results on 2D angiographies in Sect. 4.

ANALYTICAL SOLUTION FOR THE ONE-DIMENSIONAL CASE

Equation 2 is a diffusion-reaction equation whose analytical solution is not obvious without further assumptions on h and g (as defined in Sect. 1). For a better understanding, we will analyze the one-dimensional case. We consider edges as ramps which lead to plateau-like patterns in the original gradient (Fig. 2). The equation is decomposed and can be solved onto subdomains {Ω k } 0≤k≤N where g k and h k , the restrictions of g and h to Ω k , are constant. In the following developments, f k will denote the restriction of a function f to Ω k . Two cases arise. If Ω k is a homogeneous region, ∇I k = 0 so g k (x) = 1 and h k (x) = 0. Equation 2 is then the one-dimensional heat equation ∂ 2 V k ∂x 2 = 0, so the solution is a linear function:

V k (x) = m k x + p k , m k , p k ∈ R . (3)
If Ω k is a region where the gradient is non-zero, then ∇I k is constant (due to the ramp model) and so are g k and h k . Equation 2 has then the form

∂ 2 V k ∂x 2 -a 2 (V k - ∂I ∂x ) = 0 , a 2 = 1 -g k g k , 0 < g k ≤ 1 . (4)
Solutions to this second order linear equation with constant coefficients are of the form

V k (x) = c (1)
k e ax + c (2)
k e -ax + b(x), where c

(1) k , c (2) k
∈ R and b is a particular solution. Since ∇I is constant over Ω k , it satisfies the equation. Finally, the solutions on such subdomains are of the form:

V k (x) = c (1) k e ax + c (2) k e -ax + ∇I(x) . (5)
The parameters m k , p k , c

(1) k and c

(2) k

for each subdomain Ω k are given by the Dirichlet boundary condition V = 0 on ∂Ω, the C 0 and the C 1 properties of the global solution V at boundaries between the N subdomains. This yields the following linear system (in the same order): where xi denotes the point limiting Ωi-1 and Ωi, and 0 < k < N .

p1 = 0 mN xN + pN = 0 m k-1 x k + p k-1 = c (1) k e ax k + c (2) k e -ax k + v(x k ) m k+1 x k+1 + p k+1 = ac (1) k e ax k+1 -ac (2) k e -ax k+1 , (6)
If there are M plateau-like patterns, this yields a linear system of 4M + 2 equations. A numerical solution and the corresponding analytical solution, computed from a two ramps gradient, are illustrated in Fig. 2. In practice, subdomains Ω k where ∇I = 0 tend towards ∅, which means that the GVF can be approximated by a piecewiselinear function. Although this is a mere approximation, we will use this property to derive our scale measure.

DETECTION OF MEDIAL POINTS AND THEIR CORRESPONDING SCALE

The GVF energy functional in Eq. 1 contains a diffusion term which is equivalent to a multiscale analysis, from a scale-space point of view. The method proposed here is driven by two ideas. First, scale information should be available directly from the GVF, without any further multiscale analysis. Second, since all the work has been done by the GVF, recovering scales should not use overcomplicated analysis schemes of the solution.

In contrast-enhanced images, vascular structures are considered as homogeneous regions surrounded by strong gradients. In those regions, the GVF matches gradients having opposite directions, in some sense. This interpretation still holds in the one-dimensional case: thanks to the separability property of the GVF, one can consider working on the projections of the solution V along each dimension instead of working on the gradient vector field itself. It means that analyzing the d-th component V d of V along the d-th dimension only is relevant. In this outlook, the separability of the GVF and results from Sect. 2 are exploited both to detect medial points and to estimate the radius of structures.

Detection of medial points

Matching gradients having opposite directions comes down to matching projections along each dimension d having opposite signs (see Fig. 2). According to Sect. 2, the GVF may be approximated by a linear function and vanishes between those two gradients. To ensure that zero-crossings happen in the center of structures, both corresponding gradients must have exactly the same magnitude. This is why we choose to diffuse the normalized image gradient. In practice, the Point Spread Function (PSF) of the acquisition sys-tem interfere with the linearity of the solution inside homogeneous regions so that the slope of the solution V is weaker near edges. Along a given dimension d, medial points are thus detected as maxima of dV d dx d , which can still be emphasized by taking the normalized solution Ṽ . Responses are summed over all dimensions to obtain the final measure for medial points:

M = div( Ṽ ) = d d Ṽd dx d . ( 7
)

Estimation of the radius of the structures

Following the remarks formulated in the previous paragraph concerning the linear approximation, the slope of V d is inversely proportional to the radius of the structures. Let r d,k be the size of the structures along dimension d, delimited by two matching gradients V d (x k ) and V d (x k+1 ) at positions x k and x k+1 . The slope m k can be recovered where V d vanishes and the radius can be estimated as:

r d,k = V d (x k ) -V d (x k+1 ) 2m k . ( 8
)
Knowing the positions x k is not obvious. This is why previous works usually resort to an exhaustive search through multiscale analysis. On the contrary, since we are able to detect structures of interest thanks to zero-crossings, we have all the necessary information to approximate V d with a piecewise-linear function. We are only interested in the positions where two linear functions intersect, thus the approximation does not have to be accurate (see Fig. 3). A position x k responsible for linear regions Ω k and Ω k+1 with corresponding zero-crossings c k and c k+1 is thus recovered as:

x k = m k+1 c k+1 -m k c k m k+1 -m k . (9)
The actual radius r k can now be computed with simple geometrical considerations. For example, for 2D images, the radius is:

r k = r1 sin arccos r1 r 2 1 + r 2 2 , (10)
where r1 and r2 are the radii estimated along each direction.

EVALUATION OF THE ESTIMATED SCALES AND APPLICATION TO VASCULAR STRUCTURES

Equation 2 can be solved with various explicit, implicit or semiimplicit schemes. We implemented the common explicit scheme for simplicity (see Boukerroui, 2009 [8] for more efficient explicit and implicit schemes). In particular, unconditionnally stable explicit schemes exist (the Alternating Direction Explicit scheme, for example). In practice, the straightforward explicit scheme is still widely used and is very useful for investigation. We recall this scheme (Xu et al., 1998 [3]):

V n+1 i = (1 -h∆t)V n i + g∆t ∆x (V n i-1 + V n i+1 -2V n i ) + h∇I∆t , (11)
where ∆x is the spatial resolution.

Validation on synthetical vessel templates

As mentioned in Sect. 3.1, the PSF of the acquisition system and partial volume effects impact the estimation of the vessels radius.

To study their influence, we apply our algorithm to vessel templates with various radii and PSF. Vessels are modeled by convolved barlike cross-sections with radii r 0 ranging from 1 to 25 pixels, and the scale of the convolution σP SF is 0.5, 1 and 2 pixels (we approximate the PSF by a Gaussian distribution).

The relative error of the estimation with respect to the ground truth err(r) = |r-r 0 | r 0 is illustrated in Fig. 4. The algorithm introduced in Sect. 3.2 is represented by blue lines. We compare it with two other approaches. The first one, represented by red lines, is the radius evaluated by taking r = min (|x k -c k |, |x k+1 -c k+1 |). The second one, represented by the green lines, correspond to the distance from c k to the closest local maximum of V . Finally, the evaluation was performed on noise-free profiles (K = 5) in the first column, and on profiles with a 10% random additive Gaussian noise (K = 15 to compensate for the presence noise) in the second column.

It is clear that our algorithm performs better for all PSF values and is globally more robust to noise. When r 0 ≤ σP SF with no noise, the estimation is clearly unreliable but usable since the error is still less than one pixel. For r 0 > σP SF , the error is less than 3% for noise-free profiles, and remains low (around 10%) in the presence of an additive Gaussian noise . However, for radii smaller than the PSF, zero-crossings of the GVF may disappear and thus our algorithm fails to recover the structure, which corresponds to the very high errors in Fig. 4.

Skeleton extraction of vascular structures

Our algorithm was also tested to extract the skeleton of vascular structures in 2D angiographies. The medialness map M from Eq.7 and the radii are computed from the 2D GVF of the image. Seed points are selected as directional maxima of M and those lying in regions with low local contrast are discarded. Finally, centerlines are extracted as the ridges of M going through seed points. The centerlines and a segmentation reconstructed from both types of information are illustrated in Fig. 5. Most vessels are correctly recovered, with accurate radii (they are slightly overestimated in the case of very small vessels, as one should expect from Sect. 4.1).

CONCLUSION

We presented a new segmentation-free method to extract scale information of vascular structures from the GVF of an image, without any additional multiscale analysis. We demonstrated that, through fast and effective one-dimensional analysis of the GVF, we are able to devise a method which is both accurate and robust to noise. The result can serve as an input for deformable model-based algorithms, to further refine the segmentation. The current bottleneck of our approach lies in the computation of the GVF which is highly timeconsuming, as any processes involving diffusion. Efforts will be put on efficient schemes to solve this variational problem. In the future, we believe that our approach will prove to be a good alternative to multiscale analysis.

Fig. 1 :

 1 Fig. 1: Original image and its normalized GGVF.

Fig. 2 :

 2 Fig. 2: (a) Original signal and (b) the analytical solution of the GVF equation for K = 3, K = 15 and K = 30 (where K is the parameter of function g). The dotted line represents the original normalized gradient, the analytical solution is plotted in plain red, and the numeric solution is in plain blue. Both solutions overlap almost completely. The zero-crossings are preserved for all values of K but the positions of the maxima of the solution are clearly impacted.

Fig. 3 :

 3 Fig. 3: Solution to the GVF (in blue) for a one-dimensional extracted from Fig. 1 and its corresponding piecewise linear reconstruction (in red).

Fig. 4 :

 4 Fig. 4: Relative error err of the estimated radius for radii ranging from 1 to 25 pixels and a Gaussian PSF with (a-b) σP SF = 0.5, (c-d) σP SF = 1, (e-f) σP SF = 2. The first column shows the result for profiles with no noise, while a 10% random Gaussian noise has been added to vessel templates in the second column (see the text for further details).

Fig. 5 :

 5 Fig. 5: Two examples of centerlines extracted from the medialness map M and their corresponding vessel segmentation, on 2D X-ray angiographies.

As stated inXu et al., 1999 [3], the calculus of variations yields a third term ∇g(x), ∇v i in the corresponding Euler-Lagrange equation, which does not change the result much in practice.