Incorporation of InSAR and LandSat ETM+ with CGPS and SVM to assess the subsidence of Mexico City

This study presents an enhanced analysis of the subsidence rates and their effects on Mexico City. As a result of excess water withdrawal, Mexico City is experiencing subsidence. To analyze this subsidence of Mexico City, we integrated and analyzed Interferometric Synthetic Aperture Radar (InSAR), Continuous Global Positioning Systems (CGPS), and optical remote sensing data. This study utilized 52 ENVISAT-ASAR, nine GPS stations, and one Landsat ETM+ image from the Mexico City area to gather a better understanding of the subsidence rates and its effects on Mexico City’s community The InSAR data covers the period between March 2002 until June 2010, and the GPS data span the period from 1998 until 2012. We find that the maximum of 352-mm change in Line Of Sight (LOS) direction supports the outcomes from previous studies. This study shows that the maximum rate of subsidence Mexico is 352 mm/yr. 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 population density. We used SVM to compare Persistent Scatterer Interferometry (PSI) subsidence rates with the buildings’ distribution densities. This study improves the existing method by incorporating 52 ENVISAT images, ETM+, SVM classification, and 9 CGPs. This integrated study shows that the fastest subsidence zone (i.e. areas greater than 100 mm/yr), which falls into the above-mentioned temporal baseline, occurs in the areas of high and moderate building distribution density.


Introduction
Interferometric Synthetic Aperture Radar (InSAR) data have been available for geological and geomorphological analysis since the launching of ERS 1 in 1992. Persistent Scatterer Interferometry (PSI) has been employed for more than 15 years to monitor the surface of the Earth. In the last two decades, PSI and similar techniques have been proposed as well (Usai and Hanssen 1997;Hanssen 2001;Ferretti et al. 2001;Lanari et al. 2004;Salvi et al. 2004;Poreh et al. 2017). Developed by researchers at the Politecnico di Milano (POLIMI), the procedure of the PSI technique is known as Permanent/Persistent Scatterers Interferometry (Ferretti et al. 1999;Wang et al. 2010). Prior to monitoring similar terrain subsidence, several researchers used geodetic methods such as InSAR, PSI, and Continuous Global Positioning System (CGPS) (Berardino et al. 2002;Mora et al. 2003;Crosetto et al. 2005;Hooper 2008;Crosetto et al. 2008;Pepe et al. 2011;Perissin and Wang 2012;Navarro-Sanchez and Lopez-Sanchez 2013). A large number of radar images and techniques were considered for the estimation of historical changes on the Earth's surface. For example, similar methods, such as the Small Baseline Subset (SBAS), have emerged with the effectiveness of the PSI approach (Murillo and Manuel 1995;Rudolph et al. 2006;Gourmelen et al. 2007;Dai et al. 2016). The subsidence in Mexico City (Fig. 1), which is due to the extraction of groundwater ,began in the 1840s. Later, this phenomenon became extensive in the 1930s and 1950s (Carrillo 1947;Ortega-Guerrero et al. 1999;López-Quiroz et al. 2009;Rodell et al. 2009;Lopez-Quiroz et al. 2009;Ayazi et al. 2010;Ortiz-Zamora and Ortega-Guerrero 2010;Osmanoglu et al. 2011;Yan et al. 2012;Chaussard et al. 2014;Suárez et al. 2018;. Mexico City is built on highly compressible clays (Cuevas 2004). The mechanism of the subsidence in Mexico City lacks enough natural water recharge (i.e. no extracted water replacement); consequently, further compaction of the clay layers is expected (Carrillo 1947;Cuevas 2004;Osmanoglu et al. 2011;). The subsidence is not associated with the aquifer. Instead, this subsidence may exist because of the compaction of the overlying geological layers Suárez et al. 2018;. http://www.glcf.umd.edu (Cabral-Cano et al. 2008) In Mexico City, subsidence occurs when irrecoverable aquifer systems are heavily pumped. This operation results in the loss of aquifer storage and damage to engineered structures. Since the late 1950s, this subsidence has accelerated to a remarkable extent, and considerable correlated structural damage has been reported in the Mexico City area (Cabral-Cano et al. 2008;Osmanoglu et al. 2011;Yan et al. 2012). More than 20 million inhabitants in the metropolitan area are facing this extraordinary land subsidence hazard (http://www.unesco.org). For example, in Mexico City Metropolitan Cathedral (which took 250 years to build), one side is settled nearly 2.44 meters deeper than the other side. Furthermore, the cathedral is leaning to the left side (http://www.unesco.org). Cabral-Cano et al. (2008) studied Mexico City's subsidence based on the InSAR methodology.
By using ERS and ENVISAT-ASAR satellite data in the temporal baseline of 1996-2003, the authors pointed out that the subsidence rate in Mexico City has reached 370 mm/yr. It is controlled by compaction of the quaternary lacustrine clays and silts.  studied the land subsidence in major cities in the whole country including Mexico City. Suárez et al. (2018)   https://pubs.usgs.gov/of/1997/ofr-97-470/OF97-470L/graphic/data.htm  used GRACE and InSAR data sets to remotely assess groundwater storage loss in the Mexico City area. Using the SBAS-InSAR algorithm to reveal areas subject to ground motion related to groundwater overexploitation, they noted that GRACE satellite data sets do not entirely detect the significant groundwater losses. As a solution, these data sets should be combined with other high-resolution satellite imageries. Chaussard et al. (2014)  authors ran predictive simulations and field data to predict the deformation rates in the Mexico City area from 1984-1989. Ortega et al. (1993 used hydraulic data from a network of monitoring wells, geotechnical data from core samples, and historical information for 1984-1989. The authors ran a mathematical model to predict future subsidence under the current pumping rates. Strozzi and Wegmüller (1999) used ERS satellite data sets in the temporal baseline of [1995][1996][1997] to monitor subsidence in Mexico City. They found 400 mm/yr subsidence in the eastern part of Mexico City's community. They concluded that because of water extractions, almost nine meters of subsidence has occurred in the Mexico City area over the last century. Yan et al. (2012) compared PSI (with Gamma-IPTA chain methodology) and SBAS for Mexico City's subsidence mm/yr subsidence with InSAR. They also showed that the deformation is almost linear over time baselines. A PSI study based on ENVISAT-ASAR data was carried out in the eastern part of Mexico City during a short temporal baseline (i.e. 2004-2006) by Osmanoglu et al. (2011, and . They used only 23 images for PSI data analysis, which is the minimum number of sufficient imageries for PSI analysis. Finally, Sowter et al. (2016) used Sentinel for temporal baselines of 2014-2015. The SBAS (ISBAS) technique was utilized to measure the deformation rate of 240 mm/yr along the LOS direction, equivalent to over 400 mm/yr vertical rates (Zebker et al. 1997;).
However, despite the controlling procedure and the above-mentioned research, we used ten years of InSAR data in this study to focus on the maximum 352 mm/yr displacement rate (in LOS direction) occurring in the central and eastern parts of Mexico City. Our study covers a longer temporal baseline (between November 2002 and June 2010) and a wider area (62 × 56 km 2 ) than Page 6 of 35 the previous studies. In this study, we analyzed 52 ENVISAT-ASAR data from November 2002 until June 2010 with a coverage area that is twice the size as that used in Osmanoglu's work (Osmanoglu et al. 2011). We managed to extend the temporal baseline. As pointed out by Zebker et al. (1997), for N independent interferograms (N+1 InSAR imageries), the temporally uncorrelated noise reduces by a factor of 1/√ . Therefore, this statement supports the current study while enhancing the reliability performance of similar works (Cabral-Cano et al. 2008;Osmanoglu et al. 2011;.
In summary, the objectives and advantages of this study as compared to the existing research are (1) analyzing more InSAR imageries (longer temporal baselines), (2) working with a larger coverage area, (3) using accurate ENVISAT images in conjunction with more CGPS stations, and (4) conducting a detailed geohazard risk assessment of Mexico City based on SVM integration.

Geological and geotechnical settings
When 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 ( Fig. 1 and Fig. 2). The Spanish drained the lake over an extended time and expanded Mexico City onto the new land, where it exists today. Almost the entire city stands on layers of sand and clay with thicknesses up to 300 meters at some locations.
These soft, water-laden and loose sediments make the city uniquely vulnerable to subsidence, earthquakes, and other kinds 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 (Fig. 2). This volcanic zone is the result of the subduction of the Cocos and Rivera oceanic plates underneath the North America plate. Morphologically, the Mexico Basin includes (a) volcanic ranges, Page 7 of 35 composed of either polygenetic or monogenetic volcanoes; and (b) a series of fan-like knolls, located at the base of each volcanic range. The intercalation of pyroclastic and epiclastic deposits and (c) flat-land areas results from the accumulation of lacustrine sediments of variable thicknesses interbedded with tephra layers (Ferrari et al. 1999;Arce et al. 2013).
The lacustrine sediments of Mexico City exhibit very unusual behaviour; in some cases, the water content can surpass 500%, the plasticity index sometimes exceeds I >300%, and the compression index could be in the order of 10 . The friction index of Mexico City clayey soils is comparable in magnitude to that of sands. Furthermore, the sediments do not show strength loss, not even when the amplitude of the cyclic loading is as high as 80%. The sediments of Mexico City are a complex mixture of crystalline minerals and amorphous material with heterogeneous volcanic and lacustrine sediments (silty clay or clayey silt) (Ortega-Guerrero et al. 1993, Du et al. 2019). ground truth and calibration tool for the InSAR data. With more than six years (and eight years for one station) of CGPS and InSAR overlap, the two independent geodetical methodologies provide an updated picture of subsidence in Mexico City and its surrounding area in the first decade of the 21st century.

Data gathering
This study used the Landsat ETM+ image to classify the types of man-made structures and buildings' densities in the study area. Landsat ETM+ image with seven bands, acquired on 15/11/2005, was processed by the U.S. Geological Survey (USGS). We used DORIS for SVM analysis incorporation with InSAR/PSI processing techniques. We also applied supervised classification using ENVI 4.8 under Windows and Linux platforms for image processing.

GPS
All of the CGPS stations are inside the ENVISAT-ASAR data coverage (white frame in Fig.1).
UCHI, UGOL, UIGF, and UGAL are located either on the andesite-basalt lava or on the tephra deposits around Mexico City (Osmanoglu et al. 2011). CGPS data was provided by the University of Mexico and has been analyzed by utilizing the precise point positioning of the ITRF-2000 reference frame (Dixon et al. 2000;Altamimi et al. 2002;Hilley et al. 2004;Argus 2007). The GPS data have been converted from the ITRF coordinate system to the north, east, and vertical coordinates utilizing the national coordinate system. The comparison of PS with the GPS data is attempted, and the velocities are compared to a reference point. In the Mexico City area, the tectonics (i.e. plate movements) also plays a very important role, which has an undeniable effect on the InSAR reflectors, including PS's reference point. In other words, the selection of reliable coordinate systems is crucial for the fusion of InSAR and GPS.   Table 1 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, refer to Dixon et al. (2000) and Sella et al. (2002). The results of the CGPS data analysis are provided in Fig. 3 and Table 1. These data are presented in "absolute" coordinates (latitude, longitude, and height). All of the vertical (UP) components show negative values (subsidence) with a maximum of -275.3 ± 3.5 mm/yr. It identifies subsidence in the MRRA station or the easternmost CGPS station (Fig. 1). A minimum of -0.3 ± 2.5 mm/yr is observed in the UCHI station, which is located on the andesite-basalt lava or tephra deposits (see Fig. 1a and Fig. 2).
GPS data could be contaminated with several potential noises. Because the study area is small, orbital effects are assumed to be spatially uniform. Despite different temporal coverages 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 aligned with to the PSI points. The velocities measured by the PS technique are in the Line-Of-Sight (LOS) direction, and the GPS measurements are projected to this direction via:

Interferometric SAR and PSI
In this study, we used 52 ENVISAT-ASAR scenes to analyze the subsidence in the study area.
The study area covers a temporal baseline between November 2002 and June 2010. The 05/05/2006 acquisition has been selected as the master scene to minimize the effects of spatial and temporal baselines (Zebker et al. 1997;Scharroo and Visser 1998;Hanssen and Bamler 1999). The ENVISAT-ASAR imageries, which cover an area of 62 km × 56 km, are centered cover Mexico City's historic downtown. In the first step, we ran the crop by applying Delft Object-oriented Radar Interferometry Software (DORIS) method. Next, we oversampled the data by using a factor of two in range and azimuth to avoid any undersampling of the interferograms, especially during the resampling of the slave acquisition (Bell et al. 2002). By running this step, we prevented aliasing over the data. In the interferometry step, scenes were oversampled by a factor of two in range and azimuth to make each pixel (originally 4m × 20m) approximately 2 m × 10 m. DORIS was used to make differential interferograms of imageries (Hanssen and Bamler 1999). For reducing the orbital effect in the produced interferograms, we used precise orbit data (DOR and VOR) from the TU-Delft University that were supposed to minimize orbital errors significantly (Scharroo and Visser 1998). We appended useful orbit data to the images. We did not run the Porbits step (in DORIS software) in our data. Table 2 illustrates the baseline information (perpendicular, temporal, and Doppler baselines) for 52 ENVISAT-ASAR satellite imageries. We acquired images in the descending mode for the study area. Imageries are given in the format of YYYYMMDD (first and fifth columns). Major factors influencing the InSAR Phase measurements are (Hanssen 2001): φdef is the part connected to deformation, φDEM is the topographic phase contribution, and φorb is the orbital part error that could be minimized by 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 may not fall in the urban area. Finally, φn is the noise part of the phase, which is for strong scatterers and would be negligible. Data analysis and probable errors are explained in Hanssen and Bamler, (1999), Hanssen (2001), and Kampes (2005). We applied the TU-Delft approach for analyzing the data and determining the average subsidence rates, including the deformation time series for each pixel on the ground (Ferretti et al. 1999;Wang et al. 2010;Poreh et al. 2017;Du et al. 2019;  (2:3:5) to generate the RGB image. This study applied the SVM method on the RGB image to classify the populated buildings from the study area and to give a sense of the densities of the buildings (Fig. 6).
SVM is a well-known classification method that currently has wide applications in image processing and machine learning. This SVM statistical method is designed for recognizing different patterns, classification, and regression analysis. In the simplest form, SVM classification finds a hyper plan such as a ∀i〖y〗_i (x_i.w+b)-1 problem that could segregate two classes by the minimization of the following Lagrange equation (Hilley et al. 2004): In this study, we used three training sets to classify the patterns. We considered three different kinds of building density classes as the training sets. We selected three areas with three Regions Of Interest (ROI) to apply supervised SVM classification analysis. These three areas include (i) highly populated areas (dense building distribution), (ii) sparse building distribution with less density, and (iii) areas that are mountainous and unpopulated. After selecting these three ROI, we imposed the SVM method to classify these areas mentioned above with minimum error.

Results and discussion
The ability to make combined measurements from InSAR, GPS, and optical remote sensing imageries may be a powerful tool for studying the subsidence in Mexico City and the related risk management. The combination of CGPS and PSI methods seems straightforward. However, this process requires careful analysis since we do not have radar reflectors (PSs) in the exact location of each CGPS station. Therefore, the PSI approach is a relative method in comparison with the absolute GPS methodology. ± 3.5 mm/yr. The highest vertical subsidence in the 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 subsidence for UIGF station. The station UPEC, located farther to the west (see Fig. 1), shows vertical subsidence in rates of -82.3 ± 0.7 mm/yr with a negligible amount of seasonal variations. Some of the other

InSAR and PSI
In the InSAR and PSI techniques, the displacements are resolved in the Line-Of-Sight (LOS) direction. This study is not in the three orthogonal displacement vectors; therefore, the combination with CGPS stations must be dealt with very carefully.    Table 3 Comparisons of LOS' rates of movements for GPS and PSI  Table 3 illustrates the deformation rates of nine GPS stations inside the ENVISAT-ASAR data coverage area (see Figure 1a) and the closest PSI rates to those stations. It shows that the PSI and the CGPS rates are in agreement (see Fig. 9). Moreover, to summarize, more than 600,000 points have been selected as permanent scatterers in the study area. For these points, the authors calculated time series (2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)  The PSI deformation rates measure the time interval between the generated interferograms, while the CGPS stations measure in the constant rate of daily data. However, more than six (and eight for one station) years of overlaps can help us to compare CGPS stations with PSI data. Despite the poor quality and gaps for some CGPS stations (for instance UTEO and UJAL), data from other CGPS stations could be compared with the PSI-derived deformation rates. The accuracy of the PSI method can be measured in sub-centimetres if a sufficient number of imageries are utilized in PS analysis . We compared the outcomes of this research with nine independent CGPS stations to assess the accuracy of the PSI. We found that the is an improvement of the previous studies. Table 2 lists the LOS rates obtained from point positioning GPS and PSI analysis. For each GPS station, a search was conducted to find the closest PSI point's rates. PSI points close to the CGPS stations show similar rates of subsidence in the Mexico City metropolitan area. Fig. 9 shows the graph interpretation of Table 2. A one-toone line shows perfect agreement 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 PSI recorded displacements are from vertical movements of the terrain. This measurement justifies the perfect agreement between two independent displacement monitoring tools (PSI and CGPS).
Furthermore, recorded deformation 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 to the north, and MRRA is moving in a southwest direction. On the other hand, a comparison with PSI histogram (Fig. 8) shows that extraordinarily high rates of displacements (rates less than -20 mm/yr) are occurring within the Mexico City metropolitan area. Since the two independent methods, PSI and CGPS, with different time intervals are in agreement, the seasonal behaviors are negligible. and/or east directions (see Fig. 3), and the station MOCS also has a small amount of acceleration in the horizontal directions. However, the main displacements are in the vertical directions. 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. Next, the average of PSs displacement rates in distance of r = 1, 2, 4, 6, and 8 km were examined. In all of the tested CGPS stations, for r > 2 km, the deviation is remarkable and the stability would be poor.
Therefore, we selected the maximum distance from each CGPS station in r = 2 km. Statistics of the selected PSI points around the GPS antennas MOCS, MRAA, and MPAA are given in Table   4. This table depicts the number of selected PSs, Min, Max, Mean, and SD for each GPS station.

Risk assessments for existing subsidence in Mexico City
The comparison of Mexico City's subsidence monitoring history with the existing research and current study is depicted in Table 5. This study used larger InSAR and more CGPS temporal baselines in conjunction with the combination of other available remotely-sensed data. We used the SVM classification based on Landsat ETM+ imagery. We applied SVM to analyze the populated area in Mexico City and the comparison with the subsidence rates from PSI data as well as risk assessment. The RGB composite (2:3:5) of three available bands was generated, and then we oversampled the composite image with three ROIs. In this last stage, the SVM classifier was applied to the image to obtain the buildings' densities. Fig. 6 illustrates three classes of densely populated area in Mexico City. We overlaid the CGPS stations on the SVM classification map to compare this map with the terrain displacements maps. Comparison with the PSI deformation rates stress that the denser building zones are located on the highest  The existing subsidence due to over-pumping in the Mexico City metropolitan area has been examined in this study. Nelson (2000) showed that a maximum of nine meters of subsidence in an area as large as 225 km 2 had been observed in the metropolitan area. As mentioned previously, the main subsidence occurs because of water extraction and the consequent compaction of the alluvial sediments. In the subsiding area, the 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 the central and eastern parts of Mexico City. The area experiences rapid subsidence in regions with a subsidence rate greater than 50 mm/yr (Fig. 7). It correlates with the regions of intense groundwater extraction such as Texcoco sediments (Fig. 1) (Osmanoglu et al. 2011). The more stable area is located on the western side of Mexico City (in the mountains). As is evident from a comparison of Fig. 6 and Fig. 7, the highest deformation rates are located in the region with high densities of buildings. As the subsidence is developing continually in this region (see GPS results in Fig.3), the city's central and eastern parts are threatened. With a subsidence rate of 352 mm/yr, there will be a total of 3.5 meters of displacement in ten years. In this case, floods in rainy seasons will be of great concern in these two parts of the city. However, this may be an inaccurate estimation. Because the compaction of aquifers in response to water extraction depends on the aquifers' physical properties subsidence may slow down over time (Terzaghi 1925 (Murillo and Manuel 1995). Economically, the costs associated with the subsidence are enormous and include more than three billion USD due to the immediate collapse of 412 buildings and the severe damage of another 3,124 buildings (Murillo and Manuel 1995). Maintenance and geotechnical supports have indirectly led to the increase in flood risk, soil fractures, and other threats to human life. As these costs grow over time, it becomes increasingly important to assess the extension and magnitude of potential damage.
Meanwhile, the monitoring of subsidence must be continued with more GPS stations. As well, a new generation of InSAR satellites such as Sentinel, TerraSAR-X, and CosmoSkyMed is essential.

Conclusion
By comparing Mexico City's subsidence monitoring history ( The proper low spatial resolution CGPS stations (only nine stations) with east, north, and vertical components are used to calibrate the field observation assessment of the high spatial resolution PS displacement rates (more than 600,000 points). Fig. 9 shows a good correlation between CGPS and PSI data. PSI is a relative method by nature despite CGPS methodology, which provides absolute displacement rates. This is why we combined PSI with CGPS stations.
This study combined InSAR, CGPS, and optical remotely-sensed imageries to measure Mexico City's subsidence because of groundwater extraction and related risk assessments. Fifty-two ENVSAT-ASAR imageries and nine CGPS stations are used to study the deformations. Geodetic analyses based on InSAR and GPS methodologies give promising results for monitoring deformation rates in large areas. Long-term deformation rates based on PSI and GPS methodologies are similar to those found in previous InSAR and PSI works. The subsidence in the eastern part of the Mexico City metropolitan area shows fast and constant rates. This study concluded that the maximum subsidence rate is 352 mm/year in the LOS direction. This compactions and the permanent loss of porosity and reservoir capacity. The data we used in this study show that the mitigations have no effect on long-term compaction of the clay-rich aquifers, and the seasonal variations are small. The data also show that there is a considerable amount of riski in the metropolitan area of Mexico City. The subsidence (see Table 5) leads to the damage of buildings and infrastructures as well as economic ramifications. As the clay-rich aquifers shrink, fracturing, sinkholes, and faults develop and cause 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 subsidence is occurring in the highly populated zone of the Mexico City metropolitan area. PSI and SVM are two valuable methods to study these kinds of subsidence threats in the metropolitan areas.
Following further assessment of subsidence, mechanism of alluvial compaction, change of porosity and permeability, and geohazards, this study suggests that additional geophysical work is needed to map the exact extension of the subsidence and subsurface geology. Geodetic surveys using denser CGPS networks could estimate the precise amount of subsidence and its extensions.
Regarding InSAR and PSI, using the new and improved generation of InSAR images such as