InSAR and Landsat ETM + incorporating with CGPS and SVM to determine subsidence rates and effects on Mexico City

This study presents an analysis of subsidence rates and their effects on Mexico City. Mexico City is well known for its subsidence as a result of excess water withdrawal for many years. This study focuses on this problem utilizing the integration of Interferometric Synthetic Aperture Radar (InSAR), Continuous Global Positioning Systems (CGPS), and optical remote sensing data. Fifty-two ENVISAT-ASAR, nine GPS stations, and one Landsat ETM+ image from Mexico City area have been analyzed to prepare a better understanding of the subsidence rates and its effects on Mexico City’s commune. This study has utilized InSAR methods. It includes differential interferometry and Persistent Scatter Interferometry (PSI) to monitor the existing subsidence in the Mexico City area. The InSAR data covers the temporal baseline between 2002 until June 2010, and the GPS data include temporal baseline from 1998 until 2012. Maximum of 352 mm annually change in Line Of Sight (LOS) direction is in agreement with the previous geodetic studies. InSAR data have been compared with CGPS data at the same time interval. The finding of this study reveals a high amount of correlation (up to 0.98) between two independent geodetic methods. We also implemented the Support Vector Machine (SVM) analysis method based on Landsat ETM+ image to classify Mexico City’s populated density area. This method performed comparing the subsidence rates with populated area buildings. This integrated study shows that the fastest subsidence zone (i.e., areas greater than 100 mm/yr) in the over mentioned temporal baseline occurs in the high and sparsely populated areas.

spatial and temporal baselines [3,1]. The ENVISAT-ASAR imageries cover an area, something like 62 km×56 km and are centred over the Mexico City historic downtown. In DORIS method, we ran the crop step first. Then we oversampled the data by using a factor of two in range and azimuth to avoid under-sampling of the interferograms, especially during resampling of the slave acquisition [4]. By running this step, we prevent what is known as aliasing over the data. In the interferometry step, scenes are oversampled by a factor of two in range and azimuth that makes each pixel (originally 4m×20m) approximately 2m×10m.
In contrast to conventional InSAR, Persisting Scatterer Interferometry (PSI) utilizes the identification and exploitation of individual radar backscatters, or permanent scatterers, that are reflector objects smaller than the resolution pixel cell and remain coherent over a long time interval. Mexico City's ( Figure 1) subsidence due to extraction of groundwater has been studied in several papers before [5][6][7][8][9][10][11]. Mexico City is built on the highly compressible clays. The mechanism of the subsidence in Mexico City is lack of enough natural water recharge (no extracted water replacement), and consequently, compaction of the clay minerals in the fluvial sediments [12][13][14][15]. Much of the subsidence is occurring in the heavily pumped aquifer systems which are irrecoverable, and results are loss of aquifer storage and damage to engineered structures. Since the late 1950s, this subsidence accelerated remarkably, and several correlated structural damages are reported in the Mexico City area [16]. More than 20 million (http://www.unesco.org) inhabitants in the metropolitan area are in the risk of this extraordinary subsidence. For example, in Mexico City's cathedral (which took 250 years to build), one side is settled nearly 2.44 meters deeper than the other side, and the cathedral is leaning to the left side [17]. Figure1: a) Shaded relief map of Mexico City overlaid by topography map and installed CGPS stations in the metropolitan area. White Rectangle shows the ENVISAT-ASAR data coverage area in descending mode. In the right, two cross-sections (map A-A' and B-B') form SRTM topography data that have been depicted on the shaded relief map. b) 3D view of Mexico City area with Landsat ETM+ imagery form 25/11/2005. c) 3D view of the study area (topography from SRTM 3 arcsec data [18]. Geodetic methods like Interferometric Synthetic Aperture Radar (InSAR), Persistent scatterer interferometry (PSI), and Continuous Global Positioning System (CGPS) have been used prior to monitoring similar terrain subsidence. InSAR data analysis provides a new high-resolution methodology for detecting and measurement of aquifer response to pumping and recharges. InSAR uses the phase change between time separated radar imageries. Different phase contributions in the InSAR equations could be dealt with Digital Elevation Model (DEM), precise orbit data, Atmospheric Phase Screen (APS), and deformation models. Persistent Scatter Interferometry (PSI), used as a great tool for urban area displacement and civil engineering monitoring tool during the last decade. InSAR and PSI methodologies are used to measure terrain displacements in rates of few millimeters with satellites like ERSs, ENVISAT-ASAR, ALOS1 and 2, RADARSAT, and much better with high-resolution imageries like TerraSAR-X, and CosmoSkyMed satellites. Nine CGPS stations which were installed by university of Mexico in Mexico City commune since 1998 (some of them are after 2005) have been analyzed. In this study, we combine the InSAR and CGPS methodologies to get better understand of Mexico City's subsidence with geodetic based methods.
Mexico City's subsidence based on InSAR methodology carried out by Cabral-Cano et al. [16] with ERS and ENVISAT-ASAR satellites data in the temporal baseline from 1996 until 2003. They pointed out the subsidence rate in Mexico City's is 300 up to 370 mm/yr. It is primarily controlled by compaction of the Quaternary lacustrine clays and silts, partly. Despite this controlling procedure (until 2003), still a maximum amount of 352 mm/yr displacement rate (in LOS direction) have been observed in the center and eastern part of Mexico City commune.
Castellazzi et al. [11] used GRACE and InSAR data sets to assess groundwater storage loss in the Mexico area remotely. They used SBAS-InSAR algorithm to reveal areas subject to ground motion related to groundwater over-exploitation. They have noted that GRACE satellite datasets do not entirely detect the significant groundwater losses. However, it should be combined with other high resolution satellite imageries.
Lopez-Quiroz et al. [8] analyzed 38 ENVISAT images acquired between 2002 and 2007 and reached out to 400 mm/yr subsidence with InSAR. They also showed that the deformation is almost linear over time baselines.
Strozzi et al. [19] used ERS satellite datasets to monitor Mexico City's subsidence. They found 400 mm/yr subsidence in the eastern part of Mexico City's commune. They concluded that, because of water extractions, almost nine meters subsidence happened in the Mexico City area for the last century. Persistent scatterer interferometry study based on ENVISAT-ASAR data was carried out in the eastern part of Mexico City in time interval 2004 until 2006 by Osmanoglu et al. [20]. As is event, this temporal baseline is so short, and we manage to extend this temporal baseline to make a better understanding of the subsidence functionality on a more extended time series. They also used only 23 images for PSI data analysis, which is the minimum number of sufficient imageries for PSI analysis. Yan et al. [21], compared PSI (with Gamma-IPTA chain methodology) and Small Baseline Subset (SBAS) for Mexico City's subsidence and explained the strength and weakness of each method. In our study, we analyzed fifty-two ENVISAT-ASAR data between 2002 until June 2010 with area coverage as twice as Osmanoglu's work. As pointed out by Zebker et al. [22], for N independent interferograms (N+1 InSAR imageries), temporally uncorrelated noise reduces by a factor of 1/ N. Therefore, this statement confirms the current study, and it is more reliable than Osmanoglu's work [20]. For PSI analysis, TU-Delft approach was applied to get the average subsidence rate and deformation time series for each pixel on the ground. CGPS stations in the area cover SAR's data and the region nearby. It employed to compare InSAR data results with the CGPS data for calibrations and ground truth controls.
Highest subsidence in Mexico City is in order of 352 mm/yr (in LOS direction) and is happening in the city center and eastern part of Mexico City commune. The velocities from PSI methodologies are in agreement with the CGPS results with a correlation rate up to 0.98, showing a perfect validation for the applied PSI method. Support vector machine (SVM) classification based on Landsat ETM+ imagery was applied to analyze Mexico City's populated area and comparison with the subsidence rates from PSI data and risk assessment. The results show that the area with a high population is at risk of the high amount of subsidence in the eastern and central part of Mexico City. This study shows, how to combine active and passive remote sensing data to analyze deformation phenomenon and risk assessment in similar regions.
Nevertheless, the objectives and advantages of this study as compared to the existing research are (1) analyzing more InSAR imageries (longer temporal baselines), (2) covering with a larger coverage area, (3) using accurate ENVISAT images in conjunction with more CGPS stations, and (4) a detailed geohazard risk assessment of the Mexico City-based on SVM integration. This study highlights (i) The Max rate of subsidence on Mexico City is 352 mm/yr. (ii) Improved method and utilized 52 ENVISAT images, ETM+, 9 CGPS, SVM classification. (iii) High amount of correlation (up to 0.98) between two independent geodetic methods. (iv) SVM finds the fastest subsidence zone with greater than 100 mm/yr belong to the high and sparsely populated areas, and (v) temporal baseline occurs in the high and sparsely populated areas.

Tectonic and geological settings
When the Spanish invaders conquered North America in 1521, they built Mexico City over the ruins of the Aztec civilization capital of Tenochtitlan. The old Aztec city was an island in Lake Texcoco ( Figure 1). The Spanish drained the lake over extended time and expanded Mexico City onto the new land which we know it nowadays. Almost the entire city stands on layers of sand and clay which thickness up to 300 meters in some locations. These soft, water-laden and lose sediments make the city uniquely vulnerable to subsidence, earthquakes, and another kind of geohazards. The 100 km long and 80 km wide, NE-SW oriented Mexico Basin is located in the eastern sector of the Mexican Volcanic Belt (Figure 1). This volcanic zone is the result of the subduction of the Cocos and Rivera oceanic plates underneath the North America plate. The Mexico basin has formed and evolved through complex tectonic processes such as NNW and NNE-Cañón de Lobos trending reverse fault. This tectonic process has affected the Cretaceous basement and the NW-SE Mixhuca normal fault. It displaced Oligocene-Miocene volcanic strata. The tectonic process has also affected the NE-SW Tenochtitlan fault system, and it has displaced Plio-Pleistocene rocks. Finally, E-W normal faults have affected the most recent volcanic rocks, paleosols, and lacustrine sediments [23]. Morphologically, the Mexico Basin includes (a) volcanic ranges, and it composes of either polygenetic or monogenetic volcanoes, (b) a series of knolls (fan-like), and it is located at the base of each volcanic range. It is intercalation of pyroclastic and epiclastic deposits, and (c) Flat-land areas which have been resulted from the accumulation of lacustrine sediments of variable thicknesses and interbedded with tephra layers [24]. Figure 1 shows the study area at Mexico City with CGPS station locations and InSAR data coverage area. The white rectangle indicates the InSAR data coverage of ENVISAT-ASAR in descending mode. The CGPS stations are also depicted with the blue signs and red horizontal band). Two profiles A-A' and B-B' are plotted and are showing the inferred geological settings in the Mexico City area based on the previous studies [25]. Three-dimensional representation of the study area based on the Landsat ETM+ and SRTM 3 arcsec (Figure 1c) also are illustrated in Figure 1b.

Data gathering
We used complementary datasets from nine GPS stations (daily data), thousands of observations of Permanent Scatterers (PS) form ENVISAT-ASAR satellite and optical remotely sensed imageries. They provide with updated estimates of subsidence in the Mexico City and study area including related risk assessments. This deformation data of Mexico City area from GPS and InSAR acquired between 1998-2012 and 2002-June 2010, respectively. For some GPS stations, we do not have data from 1998, and some of them have gaps and bad qualities (see Figure 2). These nine GPS stations were installed and maintained by the University of Mexico. We used CGPS satiations as ground truth and calibration tool for the InSAR data. With more than eight years of CGPS and InSAR overlap, the two independent geodetical methodologies have provided an updated picture of subsidence in the Mexico City and its surrounding area in the first decade of the 21st Century.
This study used Landsat ETM+ image to classify the type of human-made structures and buildings densities in the Mexico City study area. Landsat ETM+ image with seven bands was acquired on 15/11/2005 and processed by the U.S. Geological Survey (USGS). For InSAR/PSI processing, we used Doris/ADORE, and for SVM analysis, supervised classification approach was applied by using ENVI 4.8, under Windows and Linux platforms.

GPS
Nine continues GPS stations (CGPS) inside or near the study area are employed in the Mexico City metropolitan area. All of them are inside the ENVISAT-ASAR data coverage (white rectangle in Figure 1). UCHI, UGOL, UIGF, and UGAL are located either on the Andesite-Basalt lava or Tephra deposits, around Mexico City [20]. CGPS data was provided by University of Mexico, and have been analyzed utilizing the precise point positioning of the ITRF-2000 reference frame [26][27][28]. For each CGPS station, daily time series and estimated errors were plotted ( Figure 2). This study calculated velocity and uncertainty for each CGPS station by linear regression analysis. For details of the procedures and estimation of uncertainties for GPS, measurements see Dixon et al. and Sella et al. [29][30]. CGPS data analysis results are given in Table 1. These data have been presented in "absolute" coordinates (latitude, longitude, and height). As is evident, all of UP components are showing negative values (subsidence) with a maximum of -275.3±3.5 mm/yr which is happening in the station MRRA (the most eastern located CGPS station-see Figure 1) and a minimum of -0.3±2.5 mm/yr which is happening in the station UCHI (this station is located on the Andesite-Basalt lava or Tephra deposits-see Figure 1a and Figure 2). GPS data could be contaminated with several potential noises. Because the study area is small, orbital effects are assumed spatially uniform. Despite different temporal coverage of CGPS and InSAR data, we have managed to compare these two independent geodetic methods. For comparison of GPS and PSI data, GPS observed data should be tied up to the PSI points. GPS data have been projected to the Line-Of-Sight (LOS) direction based on the following formula: Where γ is the satellite's orbit angle relative to the true north direction and θ shows the incidence angle of the radar wave [31]. In this formula, each component is known as a directional cosine. For example, LOS's vertical component has a directional cosine of cos (θ°). It means that if there were only vertical movement, the LOS velocity should be divided by cos (θ°) to represent the exact amount of the vertical component of motion. It should be considered in the LOS interpretations. Regarding this setting, we can infer LOS data as a vertical displacement with an error rate of almost 8% (cos (23°)= 0.92) for ENVISAT-ASAR satellites. Despite a relatively small amount of horizontal deformation rates in the Mexico City area (based on CGPS data) in comparison with a vertical one, we consider PSI velocities as a combination of vertical and horizontal velocities (based on equation 1). It will help us to get more accurate results.

Interferometric SAR and PSI
Persistent Scatterer Interferometry (PSI) has been employed for more than 15 years to monitor the surface of the earth. PSI and similar techniques have been proposed in the last decades [32][33][34][35][36][37], [25,31]. PSI technique developed is by researchers at the Politecnico di Milano (POLIMI), and they named the procedure, Permanent/Persistent Scatterers Interferometry [25,32]. Soon after the effectiveness of the PSI method, some other similar methods have emerged. The Small Baseline Subset (SBAS) is one of the examples [38][39]. The notable difference between these two methods is PSI analysis and acquisitions. Whereas, in the SBAS, some of them are not analyzed because their spatial baseline is too high. PSI analyses are less sensitive to geometric and temporal decorrelation compared to SBAS methods [38]. However, the number of generated interferograms in SBAS techniques is greater than in a single master approach (like PSI). The phase unwrapping procedure for SBAS and PSI is also different. In the PSI analysis, the interferograms are unwrapped temporally and spatially. SBAS has a drawback of disconnected clusters of interfergrams. However, the most important advantage of the SBAS is that it allows measuring displacements. This measurement is done for highly stable point-like scatterers and distributed scatterers (DS). In another word, we can say the areas with intermediate coherence. Therefore, several endeavours have been reported aiming to develop techniques able to combine advantages of both PSI and SBAS [40][41]. For example, the baselines minimiztion and the use of all radar images utilizing SBAS method were proposed in by Perissin et al. [42]. Mora et al. [43] reported another similar technique to SBAS for the earth surface's change monitoring. Hooper [44] and Hooper et al. [45] reported a geophysical InSAR time series approach and a stepwise linear deformation technique with least square adjustment had been reported in [46][47]. 48. In 2003, Werner et al. [48], and in 2008, Crosetto et al. [47] reported the satble poit networks and Interferometric Point Target Analysis (IPTA). Hooper [44] used multiple image pixels within a certain radius to estimate spatially correlated parameters (e.g., deformation rates, atmospheric signal delay. Ferretti et al. [34] used the SqueeSAR algorithm, capable of simultaneous analysis of PSI (i.e., PS) and DS. With a combination of PS and DS in SqueeSAR algorithm, the rural areas where the coherency is lower can be mapped. A similar algorithm, with polarimetric based radar data, was heuristically proposed by Navarro-Sanchez and Lopez-Sanchez [49]. The authors describe the different PSI implementations and main features in this article. For all, we considered a large stack of radar images for estimation of historical changes of the Earth surface's, with the employment of proper modelling techniques. The time series deformation of the scatterers and scatterers elevation are the output of PSI algorithms. PSI technique remains radar's pixels coherent during the time With PSI. When we use a large stack of SAR images (usually more than 20 SAR images), the atmospheric errors such as the Atmospheric Phase Screen (APS) can be estimated with sufficient accuracy. The proper phase correction can also be implemented to remove them. In the PSI methodologies, a single master image with specific criteria is selected (from N+1 given images), and N differential interferograms w.r.t. the master image is generated. Then, with different algorithms, the Permanent Scatterers Candidates (PSCs) are selected. Final Permanent Scatterers (PSs) can be generated by refinements of the selected PSCs, and by using PSPs [50]. When the time series of the historical records of the Earth surface's height is changing, then the height of each PS concerning a reference point for each PS point is measurable. PSI method shows promising results in urban areas, where it can achieve an average of 100PSs/〖km〗^2(points densities) with low resolution sensors like ERS1/2 and ENVISAT-ASAR, and an average of a couple of thousands PSs/〖km〗^2 with high resolution sensors like TerraSAR-X and Cosmo-SkyMed. In the PSI method, the rural/vegetated areas might not be explored properly. The absence of stable scatterers in such areas is the main reason. Another drawback of the PSI is the need for a minimum amount of images for making appropriate phases unwrapping steps, which could severely influence the degree of correctness of the selected PSC. InSAR time series method, PSI and SBAS are relative technique. It means that we can measure the calculated time series for PS points when we consider a reference point. That is assumed to be without any movements. Nonetheless, continuous GPS or levelling could resolve this problem properly [51]. The other limitation is mainly due to the observation geometry of the satellite systems. InSAR/PSI deformation rates are only measured along with the satellite LOS direction. Therefore, the value we obtain for the deformation is just the projection of the deformation vector onto the SAR look direction. Figure 2 depicts the rapid subsidence which we have observed at some stations. The red line is the best fit for the data and shows the average rate of deformation on each of CGPS stations. The idea is looking for persistent/permanent, RADAR reflectors, such as fixed dihedral structures, such as building or other similar structures. After the selection of PSC1, these strong scatterers are further filtered to detect and remove the atmospheric phase contribution, commonly referred to as the atmospheric phase screen (APS). By applying refinement of the potential scatterers for several times, the time series and height of reflectors are final products of the PSI method.
Major factors influencing the InSAR Phase measurements have been depicted in Equation 2 [31]: Table 2 illustrates the baseline information (perpendicular, temporal, and Doppler baselines) for fifty-two ENVISAT-ASAR satellite imageries that were acquired in descending mode in the Mexico City area. Imageries are given in the format of YYYYMMDD (first and fifth columns). Note: B.Temp is the Temporal baseline (time difference) between master and slave acquisitions, B.Perp is the Perpendicular baseline between orbits of the master and slave scenes, and B.Dopp is the difference between the Doppler centroid frequencies of master and slave scenes.
φdef is the part that we are interested in, φDEM is the topographic phase contribution, and φorb is the orbital part error that could be minimized with using precise orbital data. φatm is the atmospheric phase screen, and φscat is the change in the scattering attributes of the scatterers during the time, that we believe is not much in the urban area. Finally, φn is the noise part of the phase, which for strong scatterers would be negligible. Data analysis and probable errors explain in details in some books and articles (for instance, see [31,50]).
The georeferencing accuracy of standard ENVISAT-ASAR images, concerning reference ellipsoid is about 12 m in azimuth and 60 m in range direction [51]. For having as high as possible accuracy in geolocation, SRTM 3 arcsec data have been used to improve the accuracy to the scale of 10 and 15 m in azimuth and range direction, respectively. Figure 3 shows ENVISAT-ASAR data configuration we used in this study. It depicts the temporal and perpendicular baselines of the ENVISAT-ASAR data for the study area in Mexico City. Scene 20060505 has been selected as a master image to reduce the orbital error and maximize the coherence. In Figure 4, some of the interferograms are given. Interferograms are not geo-referenced (in radar coordinates). Each color cycle (or fringe) shows a movement in the rate of 2.83 cm/yr in LOS direction.   Table 2). Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 7 September 2019 doi:10.20944/preprints201909.0082.v1

Building density classification
Landsat ETM+ image (with seven bands) dated 15/11/2005 collected from University of Maryland web page [18]. It is used to classify the extension of populated buildings in the Mexico City area. The RGB image is produced with the composition of bands (2:3:5). The SVM method applied to the generated RGB image to classify the populated buildings of the study area to give a sense of buildings' densities.
SVM is a well-known classification method that nowadays applies widely in image processing and machine learning. This SVM statistical method is designed for recognizing different patterns, classification and regression analysis. In the simplest form, the SVM classification is finding a hyper plan like: ∀i〖y〗_i (x_i.w+b)-1 problem that could segregate two classes by minimization of following Lagrange equation [53]: The coefficients known as Lagrange multipliers and s are either +1 or -1. And the solution is just differentiation of the Lagrange equation and solving the equations to find Lagrange multipliers with different methodologies.
In this study, we used three training sets to classify the patterns. We supposed that they would be the main three kinds of buildings density classes in the Mexico City area. We selected three areas with three Regions Of Interest (ROI) to apply supervised SVM classification analysis. These three areas include (i) highly populated area (densely building distribution), (ii) sparsely distributed building with less density, and (iii) area that are not populated or they are the mountainous area. After the selection of these three ROI, we imposed the SVM method to classify these areas mentioned above with minimum error. Figure 5 shows the result of the SVM method with the linear kernel of size 5×5, for the Mexico City area. As it depicts in this figure, the red area represents highly populated, green illustrates less populated, and the blue color shows none populated or mountainous area. Blue areas roughly are correlated with the mountains around Mexico City (Fig.  1). Later on, we examined the PSI data with this image to find out the risk of subsidence in the Mexico City area with the combination of InSAR, GPS, and Optical Landsat ETM+ remotely sensed imagery.
The RGB composite (2:3:5) of three available bands was generated and, then we oversampled the composite image with three ROI. In this last stage, the SVM classifier was applied to the image to get the buildings densities. Three main populated buildings are shown in Figure 5. We overlaid CGPS stations on the SVM classification map to compare this map with the terrain displacements maps.

Results and Discussion
The ability to make combined measurements from InSAR, GPS, and optical remote sensing imageries may be a powerful tool for studying Mexico City's subsidence and related risk management. The combination of CGPS and PSI methods seems straightforward; however, it requires careful analysis, since we do not have radar reflectors (PSs) in the exact location of each CGPS stations, and also PSI approach is a relative method in comparison with absolute GPS methodology.

CGPS
As pointed out by Osmanoglu et al. [20] and Cabral-Cano et al. [16], the eastern and central part of the Mexico City (i.e., the most populated area in the Mexico City commune) underwent to a high amount of subsidence. CGPS results of Mexico City's displacements are depicted in Table 1 and Figure 2 show extraordinary subsiding in the central and eastern part of Mexico City commune. All CGPS stations are showing negative values (subsidence) in UP (vertical) components ( Figure 2). This subsidence contains the regional tectonic movements and local subsiding because of the water withdraws. Natural compaction of basins is causing slow subsidence in rates of few millimetres per year in the area like Mexico City. While in comparison, subsidence based on groundwater extractions has rapid rates of tens centimetres per year [54,16].
For almost all of the CGPS stations in the Mexico City study area, seasonal variations in UP components are negligible. The GPS data shows high agreements with the previous study in this area as compared to previous studies [16,20].
Permanent GPS station UIGF (see Figure 2)  2005 with high temporal resolutions. These CGPS stations located in the high subsidence region. They record the vertical subsidence in the range of -164.7±2.8 and -275.3±3.5 mm/yr. The highest vertical subsidence in the Mexico City study area belongs to the station MRRA with a rate of -275.3±3.5 mm/yr. We recorded -2.2 ±2.7 mm/yr vertical subsidences for UIGF station. They locate in the zone out of the effective subsidence. The station UPEC, which is located further to the west (see Figure 1) shows vertical subsidence in rates of -82.3 ±0.7 mm/yr with a negligible amount of seasonal variations. Some of the other CGPS stations are showing bad data distributions and hard to rely on them. For instance, in the stations, UTEO, UJAL, and UCHI lots of gaps and missing data exist.

InSAR and PSI
In the InSAR and PSI techniques, the displacements are resolved in the Line-Of-Sight (LOS) direction. It is not in the three orthogonal displacement vectors; therefore, the combination with CGPS stations must be dealt very carefully. PSI rates in the Line-Of-Sight (LOS) direction for Mexico City's subsidence during 2002-2010 is illustrated in Figure 6. Most of the available PSs time series does not address the issue of seasonal variability, leading to this fact that the amount of aquifer recharge is meager. Rate of the PSI is increasing eastward, showing a massive amount of movement towards the center of the basin where the clay-rich sediment package is thickest (see profile A-A' in Figure 1a).
This study used the phase unwrapping method to solve the phase history of a single-pixel called periodogram [20]. As it is pointed in the previous studies (for instance, [8,20]), Mexico City's subsidence is almost linear in time. Therefore, the best model seems to be the linear model. Note that in PSI data analysis philosophy, non-linear displacements between neighbouring scatterers are negligible. Nonetheless, different deformation models and statistical hypothesis tests have been applied in unwrapping steps to minimize the phase ambiguity errors. Reference height for PSI analysis is fixed to 2240 meters to reduce the topographic phase errors efficiently. PSC selections method is based on thresholding of dispersion index, and the selection grid size is fixed to 300 meters. Input data is as big as one T bytes and the required time to analysis the datasets based on TU-DELFT approach was 15 days in a PC with four cores and 48 Gigabytes of RAMs.  Table 3 illustrates 9 GPS stations inside the ENVISAT-ASAR data coverage area (see Figure 1a), and the closest InSAR/PSI rates to those stations. It shows that the PSI and the CGPS rates are in agreements (see Figure 8). Moreover, in a nutshell, more than 600,000 points have been selected as permanent scatterer in the Mexico City study area. For these points, authors calculated time series (2002-2010) and height of radar scatterers. The PSs density in the study area is 240 PSs/km2. Figure   Preprints  7 shows the histogram of LOS deformation rates, which is strongly skewed to negative values. The J shape distribution of the deformation rates reveals the subsidence related to the groundwater extraction is the dominant process in the Mexico City study area (in comparison with tectonic processes). The other factors, like tectonic movements of the plates-which are mostly in the north direction have a small effect on the existing PSI rates. The most common deformation rates are between -20 and +2 mm/yr, which suggest the existing of mixture igneous rock and alluvial sediments dominants. In other words, the existing displacements because of the water extractions are restricted to the fluvial sediments, and the parts with igneous rocks, are more stable. A comparison with previous InSAR and PSI works shows a high amount of agreement with the shape and rate of the subsidence [16,20]. This study covered much more temporal time -from 2002 until June 2010-and much more study area coverage than previous studies (62×56 km2) (For example compare with Osmanoglu et al., [20]). The subsidence rate is classified in to six classes: a) Less than -250 mm/year (red color), b) -250 up to -200, c) -200 up to -150, d) -150 up to -100, e) -100 up to -50, f) -50 up to 0, g) 0 up to +100 (blue color). Each class has a length of 50 mm/year interval (histogram in Figure 7 has bin size of 5 mm/yr). Nine GPS stations location also overlaid on Figure 6. For each pixel, time series of displacements in temporal baseline (2002-June 2010) has been plotted, which would be used later on to compare GPS with PSI data for calibration and validation. The histogram strongly skewed to negative values and has J shape, leading to this fact that almost all of the subsidence is happening because of the groundwater extractions in the Mexico City area.

Comparison of PSI with CGPS stations and validation
The InSAR/PSI deformation rates are measured in the time interval between the generated interferograms, but the CGPS stations are in the constant rate of daily data. However, more than eight years of overlaps can help us to compare CGPS stations with PSI data. In spite of poor qualities and gaps for some CGPS stations (for instance UTEO, and UJAL), the others could be compared with PSI derived deformations rates.
The accuracy of the PSI methodology can be sub-centimetres if a sufficient number of imageries are utilized in PS analysis [34]. To assess the accuracy of the PSI in the study area, this research has compared the results with nine independent CGPS stations which are available in the study area. Table 2 lists the LOS rates obtained form point positioning GPS, and PSI analysis. For each GPS station, a search has been conducted to find the closest PSI point's rate. PSI points close to the CGPS stations show a similar rate of subsidence in the Mexico City's metropolitan area. Figure 8 shows the graph interpretation of Table 2. A one to one line shows perfect agreements of correlation between GPS and PSI points, leading to acceptable PSI analysis in the Mexico City study area.
Because of the low incidence angle (~ 23°) of the ENVISAT-ASAR imageries, most of the InSAR/PSI recorded displacements are from vertical movements of the terrain. They justify the perfect agreement between two independent displacement monitoring tools (InSAR/PSI and CGPS).
Besides, recorded deformations rates from independent CGPS stations of MOCS, MPAA, MRRA, and UPEC confirm that the horizontal movements are small and do not follow a preferred direction. For instance, MPAA and UPEC are horizontally moving in the northwest direction. MOCS is moving horizontally almost in north and MRRA is moving in a southwest direction. On the other hand, comparison with PSI histogram (Figure 7) shows that extraordinary high rates of displacements (rates less than -20 mm/yr) that are occurring within the Mexico City metropolitan area. Since the two independent methods, InSAR/PSI and CGPS, with different time intervals are in agreements, and the seasonal bias would be negligible.

LOS GPS and average PSI comparison
For each CGPS station, a subset of PSs located nearby of the CGPS stations has been selected to test the stability of CGPS stations with its surroundings. The comparison between CGPS and PSI displacement results shows in Figure 9. For this propose, this study selected CGPS stations MOCS, MPAA, and MRRA, which were installed after 2004 and are located in the high zone of subsidence. The MPAA and MRRA are showing an obvious acceleration in the north and/or east directions (see Figure 2), and the station MOCS also has a small amount of acceleration in horizontal directions. However, the main displacements are in the vertical direction. In order to tie the PSI results to CGPS data, this study projected the north, east, and vertical directions to the LOS direction (via equation 1) with directional cosine. Then the average of PSs displacement rates in distance of r =1, 2, 4, 6, and 8 km have been examined. In all of the tested CGPS stations, for r > 2 km deviation is remarkable and stability would be poor. Therefore, we stick to the maximum distance from each CGPS station in r =2 km. For three CGPS stations, CGPS data was projected on the LOS direction by Equation 1 and averaged PSI rate in subsets for r = 1, 2 km was calculated. As is evident, for MOCS, and MRRA, average PSs rates after 2 kilometres are deviating from CGPS displacement rates. MPPA stays more or less stable to the distance changes.

Preprints
For station MPAA, which is the most southern CGPS station from three stations above a perfect coloration has been observed. It is leading to this fact that until a maximum of 2 km distance from MPAA, this CGPS station has valid results. A small amount of deviation from this CGPS station after 2009 (June) has been observed. For CGPS stations MOCS and MRAA correlations are good until 2007, and from then, for both stations in r=1, 2 km small deviations have been observed.
The quick results from this analysis are (a) the slop for both CGPS and average PSI are similar. It is leading to the same amount of subsidence rate with two independent methods. (b) The effective and valid area for each CGPS station is approximately 2 km. In other words, we can trust the GPS and PSI correlation for a circle with a radius of maximum 2 km.

Risk assessments for existing subsidence in Mexico City
In Mexico City area the soil deformations are classified into (1) sudden subsidence and (2) slow subsidence. However, slow subsidence usually makes enormous economic and human-related disasters. Cities built on unconsolidated clays, silts, peats, or sands are in the danger of sudden or slow subsidence. Extreme groundwater extractions, flooding, tsunami, and an earthquake made morphological settings in Mexico basin in a dangerous situation of subsidence. Most often, the buildings and streets add weights to the region and intensify the soil's stress even more. Meanwhile, finding the buildings extensions and their relationship with ongoing subsidence is crucial.
The existing subsidence in the Mexico City metropolitan area because of the over-pumping has been studied in this study. Maximum of nine meters of subsidence in an area as big as 225 km 2 has been observed in Mexico City's metropolitan area [55]. As we mentioned earlier, the main subsidence is happening because of the water extraction, and consequently compaction of the alluvial sediments. In the subsiding area, intergranular pressure of aquifers decreases. The depletion of water at depth is the cause of the observed subsidence. The highest subsidence rate is located in Mexico City's central and eastern part. The area experiences rapid subsidence (regions with subsidence rate of greater than 50 mm/yr in Figure 6). It correlates with the regions of intense groundwater extraction (Texcoco sediments in Figure 1). The more stable area is located to the western side of Mexico City (to the mountains). As is evident from the comparison of Figure 5 and Figure 6, the highest deformation rates are located in the region that Mexico City's buildings show high densities. As the subsidence is developing continually (see GPS results in Figure 2) in this region, the city center and eastern part of the city is in a real treat. For highest subsidence in the rate of 352 mm/yr, in ten years, it would reach in a total of 3.5 meters of displacement, which in rainy sessions, floods would make real problems. However, it might be an inaccurate estimation. Because compaction of aquifers in response to water extraction, depend on physical properties of aquifers (see for instance [56]). It varies during the time (slowing down of subsidence with time). With more than 2500 m of sediments [57] in Mexico basin, estimation of stopping point of subsidence is not easy. We need more subsurface data for more accurate estimations analysis. The weaken subsoil, which mostly made up elasto-plastic clays minerals, has the capacity for instantaneous and high compressions [56]. Meanwhile, results of exceeding water withdrawal have made Mexico City area susceptible to any geohazards treats. For instance, in 1985 an 8.1 earthquake (18.2° N, 102.7° W) with the epicentre of 350 km away from Mexico City at least 40,000 people died. More than three billion USD damaged, 412 buildings collapsed completely, and another 3,124 were severely damaged [5]. In this earthquake, the central part of Mexico City suffered from an average vertical displacement of 30 cm [5]. Economically, the costs associated with the subsidence are enormous. Because of maintenance and geotechnical supports it has indirectly given the increase in flood risk, soil fractures, and other threats to human life. As these costs grow in time, it becomes increasingly important to assess the extension and magnitude of the likely happening damages. Meanwhile, it continues monitoring of subsidence with more GPS stations, and new generation InSAR satellites like Sentinel, TerraSAR-X, and CosmoSkyMed is essential.

Conclusions
Comparison of Mexico City subsidence monitoring history (Table 4) v.s existing research and the current study, this study contributed (a) using of bigger InSAR temporal baselines, (b) more and bigger CGPS temporal baselines, and (c) combination with other available remotely sensed data utilizing SVM. Coupling GPS, InSAR, and optical remotely sensed analysis has the advantage of accurately monitoring of subsidence of Mexico City's area, by exploiting the strengths and minimizing the weakness of each technique.
The proper low spatial resolution CGPS stations (only nine stations) with east, north, and up components are used to calibrate and ground truth assessment of the high spatial resolution PS displacement rates (more than 600,000 points InSAR/PSI analysis results show almost perfect agreement with the GPS data at R2 in the order of 0.98 in most of the CGPS stations. . The fast subsidence rates in the Mexico City metropolitan area is a result of a massive amount of well pumping. It has the consequence of clay-rich aquifers compactions and permanent loss of porosity and reservoir capacity forever. The data we used in this study showed that the mitigations have no effect on long term compaction of the clay-rich aquifers and the seasonal variations are small. It also shows a considerable amount of risk in the metropolitan of Mexico City area. The subsidence (see Table 4) effects on the damage of buildings and infrastructures, something that economically must be taken care of. As the clay-rich aquifers shrink, fracturing, sinkholes, and faults are developing and have direct effects on the engineered structures.
This study used SVM classification for classifying the populated area and finding its correlation with the high PSI deformation rates. The maximum amount of the subsidence is happening in the highly populated zone of the Mexico City metropolitan area. PSI and SVM are two valuable methods to study this kind of subsidence threats in the metropolitan areas.
Further assessing of subsidence, mechanism of alluvial compaction, change of porosity and permeability, and geohazards, this study suggests additional geophysical work to map the exact extension of the subsidence and subsurface geology. Geodetic surveys using denser CGPS networks could help to estimate the precise amount of subsidence and its extensions. In InSAR and PSI sense, using the better and new generation InSAR images like TerraSAR-X and CosmoSkyMed imageries could provide some more details for the existing subsidence. Landsat ETM+ imageries download from Global Land Cover Facility (GLCF) web page (http://www.glcf.umd.edu). We also thank Enrique Cabral Cano from University of Mexico City for providing with GPS data. In addition, we appreciate the research unit of the Southwest Jiaotong University for funding support.

Conflicts of Interest:
The authors declare no conflict of interest.
Data Availability: All data are available upon request.