Surface Deformation Monitoring in Zhengzhou City from 2014 to 2016 Using Time-Series InSAR

In recent years, with the development of urban expansion in Zhengzhou city, the underground resources, such as underground water and coal mining, have been exploited greatly, which have resulted in ground subsidence and several environmental issues. In order to study the spatial distribution and temporal changes of ground subsidence of Zhengzhou city, the Interferometric Synthetic Aperture Radar (InSAR) time series analysis technique combining persistent scatterers (PSs) and distributed scatterers (DSs) was proposed and applied. In particular, the orbit and topographic related atmospheric phase errors have been corrected by a phase ramp correction method. Furthermore, the deformation parameters of PSs and DSs are retrieved based on a layered strategy. The deformation and DEM error of PSs are first estimated using conventional PSI method. Then the deformation parameters of DSs are retrieved using an adaptive searching window based on the initial results of PSs. Experimental results show that ground deformation of the study area could be retrieved by the proposed method and the ground deformation is widespread and unevenly distributed with large differences. The deformation rate ranges from −55 to 10 mm/year, and the standard deviation of the results is about 8 mm/year. The observed InSAR results reveal that most of the subsidence areas are in the north and northeast of Zhengzhou city. Furthermore, it is found that the possible factors resulting in the ground subsidence include sediment consolidation, water exploitation, and urban expansion. The result could provide significant information to serve the land subsidence mitigation in Zhengzhou city.


Introduction
With the rapid economy development and urban expansion, the demand for underground resources, such as water, increase greatly in China recently [1][2][3][4].Due to insufficient rainfall in the North China Plain (NCP), groundwater has become the main water resources for farmland irrigation, industrial activities, and daily water supply of urban inhabitants since 1970s [5,6].With the long-term over-withdrawal of groundwater, the NCP is suffering from ground subsidence due to human activities, such as underground water, oil and coal mining exploitation, and urban expansion [7][8][9].It was found that most land subsidence in middle-east plain of the NCP were caused by extraction of deep confined groundwater [7].Peng et al. (2016) found that more than 1000 ground fissures were distributed on about 300 sites on the NCP [10].The land subsidence has become one of the major geological hazards in the NCP.Zhengzhou city experienced a rapid growth in population and socioeconomic development, resulting in a great expansion of the urban area with the urban construction land increase from 98.57 km 2 to 221.11 km 2 during last two or three decades [11].Previous studies show that large subsidence funnel have been found in Zhengzhou and the deformation areas were expanding year by year [7,12].Such serious subsidence could destroy surface urban infrastructure, building and farmland, which limited the development of underground space.Therefore, it is important to monitor the ground subsidence and understand the distribution of the subsidence in Zhengzhou city.
Conventional methods for ground deformation measurement include ground levelling and GPS techniques could provide accuracy deformation information at point scale.Those point-based methods have been established in many areas for deformation monitoring [13,14].However, these methods cannot provide a deformation map with high spatial sampling density.Compared with those single-point-measurement methods, Synthetic Aperture Radar Interferometry (InSAR) is a powerful technology that can obtain highly precise elevation and surface deformation information along the line of sight (LOS) using phase information over wide areas.InSAR technique has been applied in many fields, such as earthquake, volcano, surface subsidence, and permafrost [15][16][17].Nevertheless, the InSAR technique often suffers from spatial, temporal and atmospheric decorrelations, which limit the application of InSAR [18].
To overcome the shortages of InSAR, several advanced time-series InSAR techniques, such as Persistent Scatterers Interferometry (PSI) [19], the Stanford Method for Persistent Scatterers (StaMPS) [20], The Coherent Pixels Technique (CPT) [21], and the Small Baseline Subset (SBAS) algorithm [22], have been proposed in recent years.Those techniques could monitor wide area deformation with millimeter accuracy by exploiting persistent scatterers (PSs) (such as buildings, man-made structures) or high coherence targets within resolution cells.Those point-targets are stable scattering points, which should remain coherent in the whole observation period.Those time-series InSAR techniques have been successfully applied in ground subsidence monitoring, due to underground resources exploitation and urban exploitation [23][24][25].However, those advanced time-series InSAR methods have limitations on the suburb, due to the lack of enough measured points.In recent years, methods have been proposed by exploring distributed scatterers (DSs) [26][27][28], Quasi Persistent Scatterers (QPS) [29], Temporarily Coherent Point (TCP) [30], to overcome the limitation of insufficient coherence points in sub-urban and natural areas.Those techniques extended the measured points from persistent or coherent point to DS by using special selection and filtering techniques.In the plain area, the DSs are widely distributed, which could be used as the measurement points.
The goal of this paper is to apply time-series InSAR technique in the ground deformation monitoring of Zhengzhou city with PSs and DSs and try to find the possible factors affecting the ground motion.In particular, a polynomial method has been introduced to correct the phase ramps in the interferograms, due to the orbit errors and atmospheric delays.Furthermore, the deformation parameters of the DSs are retrieved using a region-grown method with an adaptive searching window based on the initial results of PSs. 17 C-band Radarsat-2 images acquired from February 2014 to October 2016 have been collected and processed in this study.

Study Area
Zhengzhou is selected as our study area, which is located in the central of Henan province and situated at the transitional zone among piedmont, torrential plain and the Yellow River alluvial plain (Figure 1).It belongs to the north China stratigraphic area, the NCP stratigraphic division, with the thickness of Cenozoic sedimentation ranging from 500 to 2200 m, the thickness of the quaternary sedimentation 60~80 m, southwest, northeast thick, Neogene 400~1900 m sedimentation thickness [31].Zhengzhou strata belong to the North China strata, mainly for the Quaternary coverage.Affected by the erosion of Yellow River, geomorphologic types of Zhengzhou city can be classified into three classes: Loess platform, flood plains, and sandy dunes (Figure 2).With the rapid urban expansion, the construction of infrastructure and huge of long-term large foundation pit dewatering engineering have increased the demand of water, which contributed to the reduction of shallow groundwater reserves.Groundwater use in Zhengzhou has increased by more than 50% over the last quarter century, with the largest increase over the years 1980-1995 [6].Together with the decreasing water supplement of Yellow River and other rivers, many deformation areas have been detected in Zhengzhou [7].
Remote Sens. 2018, 8, x 3 of 19 groundwater reserves.Groundwater use in Zhengzhou has increased by more than 50% over the last quarter century, with the largest increase over the years 1980-1995 [6].Together with the decreasing water supplement of Yellow River and other rivers, many deformation areas have been detected in Zhengzhou [7].groundwater reserves.Groundwater use in Zhengzhou has increased by more than 50% over the last quarter century, with the largest increase over the years 1980-1995 [6].Together with the decreasing water supplement of Yellow River and other rivers, many deformation areas have been detected in Zhengzhou [7].As shown in Figure 1, there are many mountains and hills in the west of Zhengzhou and the east is plain.The elevation decreases from west to east.The Yellow River runs from west to east through Henan province.The average annual precipitation is about 500-900 mm and 50% of the precipitation is concentrated in summer.Previous studies show that the maximum subsidence rate is over −50 mm/year in the northern part of Zhengzhou city, and the ground surface continues to subside [7,25].

Dataset
To monitor the ground deformation of the study area, 17 Radarsat-2 ascending images acquired from February 2014 to October 2016 with an incidence angle of 31.2 degrees and 'HH' polarization have been used (Table 1).The coverage of the Radarsat-2 images is presented as the yellow rectangular box in Figure 1a.The mean SAR amplitude image of the study area (red rectangular box in Figure 1a, approximately 150 × 180 km 2 ) from all the images is shown in Figure 1b.The pixel spacings of range and azimuth are 12 m and 5 m, respectively.To restrain the effect of spatial and temporal decorrelations, forty-seven interferometric pairs are generated with the perpendicular baseline less than 250 m and the temporal baseline less than 200 days.Figure 3 shows the combination of the selected interferometric pairs.The Shuttle Radar Topography Mission (SRTM) DEM data with a 30 m (1 arc-second) grid [32] is used to remove the topographic phase.As shown in Figure 1, there are many mountains and hills in the west of Zhengzhou and the east is plain.The elevation decreases from west to east.The Yellow River runs from west to east through Henan province.The average annual precipitation is about 500-900 mm and 50% of the precipitation is concentrated in summer.Previous studies show that the maximum subsidence rate is over −50 mm/year in the northern part of Zhengzhou city, and the ground surface continues to subside [7,25].

Dataset
To monitor the ground deformation of the study area, 17 Radarsat-2 ascending images acquired from February 2014 to October 2016 with an incidence angle of 31.2 degrees and 'HH' polarization have been used (Table 1).The coverage of the Radarsat-2 images is presented as the yellow rectangular box in Figure 1a.The mean SAR amplitude image of the study area (red rectangular box in Figure 1a, approximately 150 × 180 km 2 ) from all the images is shown in Figure 1b.The pixel spacings of range and azimuth are 12 m and 5 m, respectively.To restrain the effect of spatial and temporal decorrelations, forty-seven interferometric pairs are generated with the perpendicular baseline less than 250 m and the temporal baseline less than 200 days.Figure 3 shows the combination of the selected interferometric pairs.The Shuttle Radar Topography Mission (SRTM) DEM data with a 30 m (1 arc-second) grid [32] is used to remove the topographic phase.

Methodology
The objective of this paper is to retrieve the ground subsidence of Zhengzhou city and using time-series interferometry with PSs and DSs.The processing flow chart of the method is shown in Figure 4.In this section, the main steps of the proposed method in this paper includes three steps: (1) PSs and DSs selection, (2) Phase ramps correction, and (3) deformation estimation.

Methodology
The objective of this paper is to retrieve the ground subsidence of Zhengzhou city and using time-series interferometry with PSs and DSs.The processing flow chart of the method is shown in Figure 4.In this section, the main steps of the proposed method in this paper includes three steps: (1) PSs and DSs selection, (2) Phase ramps correction, and (3) deformation estimation.

PSs and DSs Selections
First, all the SAR images are co-registered to the master image (the first SAR image acquired at 20 February 2014).Then, the interferograms are generated based on the selected interferometric pairs (Figure 3).After that, the SRTM DEM of the study is converted from geographic coordinate system to radar coordinate system.The differential interferograms are obtained by removing the topographic

PSs and DSs Selections
First, all the SAR images are co-registered to the master image (the first SAR image acquired at 20 February 2014).Then, the interferograms are generated based on the selected interferometric pairs (Figure 3).After that, the SRTM DEM of the study is converted from geographic coordinate system to radar coordinate system.The differential interferograms are obtained by removing the topographic phase from the interferograms.After the differential interferograms have been generated, the deformation parameters could be retrieved on the selected measured points.In this paper, the PSs and DSs are both selected as measured points to improve the density of the measured points.PSs are those ground objects that can hold high coherence in all the interferograms.Several methods have been proposed and used for identifying the PSs, such as amplitude dispersion index and coherence stability [19,33].In this paper, the coherence method is adopted for PSs selection.
DSs are statistical homogeneous pixels (SHP) sharing the same scatter behavior and often corresponding to areas of non-cultivated land with short vegetation and bare soil.In this paper, the DSs are selected combing classified information and statistic information [28].The DSs are selection by two steps.Considering the classes of DSs, the DSs candidates are firstly selected from the classification map generated by the unsupervised Iterative Self-Organizing Data Analysis Technique algorithm (ISODATA) method.For the selected candidates, an adaptive sample selection strategy (such as the Anderson-Darling (AD) test) is applied to refine them [27].More description about the DSs selection method used in this paper can be found in Reference [28].

Correction for Phase Ramps
The first step of the time-series InSAR method is to generate the interferometric pairs with small spatial and temporal baselines to decrease the orbital errors and temporal decorrelations.For the Radarsat-2 images, the orbital error could reach several centimeters in some cases, which could cause serious errors in final deformation result, especially in mountainous areas [34].Therefore, it is necessary to remove those phase ramps.The orbital phase ramps can always be estimated by using a bilinear or biquadratic model [35][36][37].In mountainous area, the stratified troposphere could result in serious atmospheric delays in interferograms, which could introduce several centimeters errors into the final displacement result.The topographic related atmospheric delay can be removed by a linear model based on elevation data [35].Considering the above errors, a combined model comprising a biquadratic model for orbital phase errors and a linear model for the elevation dependent errors (stratified atmospheric delay and topography errors) was applied.For each pixel in interferograms, the phase ramps can be modeled by: where ϕ(x, y) ramp is the modeled phase ramps.x and y are the range and azimuth coordinate of the pixel; h is the elevation of the pixel; ε(x, y) represents the random phase error; a i represents the estimated parameters.After removing the phase ramps from interferograms, the deformation parameters can be retrieved by the time-series InSAR method, which would be described in the following section.

Deformation Estimation on PSs and DSs
After the generation of interferograms and identification of PSs and DSs, the deformation parameters of the selected measured points can be obtained through conventional PSI method [19].Even though the atmospheric delay can be reduced trough the linear model in (1), those selected point should be connected to further eliminate the effects of the atmosphere.Generally, the phase difference between two neighboring points, x and y, on each edge of the network in the kth interferogram can be expressed as the following phase model: where T k , B k ⊥ , R and θ denote the time baseline, normal baseline, slant range, and incidence angle, respectively; λ is the wavelength, and ∆v and ∆ε are the deformation rate variation and the elevation error variation of two neighboring pixels.The differential height error and deformation velocity between the two points of the edge can be generated by maximizing process [33].
where N is the number of the interferograms.∆ϕ k phase is the differential phase of the two points in the kth interferogram.The maximum value is called temporal coherence.After this maximization processing of each edge, a quality test is performed to reject links with temporal coherence lower than the threshold.In this study, the temporal coherence threshold is set to 0.7.
The deformation rate and DEM error of the PSs can be estimated through the integration step, like the CPT [29,33].Considering the physical properties of the PSs and DSs, a layered strategy is used in this study.The valid deformation results from the PS set are used as the reference data for the parameter estimation of the DS set.In our latter experiment, the valid PSs are those points with the temporal coherence value larger than 0.7.Each of the valid PSs is referred as processed object (PO).
For the deformation estimation of DSs, we only point out the major differences from the previous study [38].For each DS, referred as unprocessed object (UPO), we used an adaptive search window not a fixing search window to search the neighboring POs from the PO set through an iterative process.Figure 5 shows the process of searching neighboring POs.In this paper, the searching process will be executed iteratively with the searching window size of 2n + 1 (n = 1, 2, 3, . . .), respectively.During the searching process, when the number of the selected neighboring POs reaches the set threshold, the searching process will be finished.In this paper, the threshold is set to 20.Using the adaptive search window, more neighboring POs will be identified for the processing UPO.After selecting the neighboring POs, the ∆v and ∆ε of the links between UPO and PO can be estimated using Equation (2).Then the displacement rate and DEM error of the UPO can be obtained through the following equations: where V C and ε C are the estimated deformation rate and DEM error; M is the number of validated links; v i and ε i are the deformation rate and DEM error of the ith PO, respectively; ∆v i , ∆ε i , and ξ i denote the deformation rate, DEM error, and temporal coherence of the link between UPO and the ith PO, respectively.In the same way, all the UPOs are processed to retrieve the deformation rate and DEM error.After that, the processed and valid point is added in the PO sets.Once deformation rate and DEM error of all the selected points (PSs and DSs) have been estimated, the residue phase can be generated, which contains atmospheric phase, nonlinear phase, and noise components.The atmosphere and noise phase components are characterized by high space correlation, but show a low temporal correlation, which can be separated using spatial and temporal filters [33].After removing the atmosphere and noise phases component, the nonlinear deformation component can be obtained.Finally, by integrating the linear and residue nonlinear deformation components, a full deformation map of the study area is obtained.
There are two advantages to apply the layered strategy in jointly processing PSs and DSs.Firstly, it could control error propagation.A layered strategy is introduced in conventional PSI to deal with PSs and DSs in a multi-layer schema from high phase quality to low.Secondly, it will improve computing efficiency.Because of the large number of DSs, it is computationally expensive and difficult to retrieve deformation parameters if constructing Delaunay triangulation on PSs and DSs simultaneously.Once deformation rate and DEM error of all the selected points (PSs and DSs) have been estimated, the residue phase can be generated, which contains atmospheric phase, nonlinear phase, and noise components.The atmosphere and noise phase components are characterized by high space correlation, but show a low temporal correlation, which can be separated using spatial and temporal filters [33].After removing the atmosphere and noise phases component, the nonlinear deformation component can be obtained.Finally, by integrating the linear and residue nonlinear deformation components, a full deformation map of the study area is obtained.

Phase Ramps Correction
There are two advantages to apply the layered strategy in jointly processing PSs and DSs.Firstly, it could control error propagation.A layered strategy is introduced in conventional PSI to deal with PSs and DSs in a multi-layer schema from high phase quality to low.Secondly, it will improve computing efficiency.Because of the large number of DSs, it is computationally expensive and difficult to retrieve deformation parameters if constructing Delaunay triangulation on PSs and DSs simultaneously.

Phase Ramps Correction
Three typical interferometric pairs with obvious orbit errors and topographic error are shown in Figure 6a-c.It would produce significant artifacts in the final displacement result if those phase ramps are not removed.Figure 6d-f indicate the corrected interferograms using the phase ramps correction model (1).It can be seen that most of the phase ramps (orbit error, topographic related atmospheric error) have been removed and the deformation signal are retained and observed in the corrected interferograms.In practically, there are about 10 orbit error fringes throughout the interferometric pair 20141018-20150215 (Figure 6c) before correction, which make the deformation signal weaken.In the interferograms after correction (Figure 6f), two subsidence centers are clear to find in the central of the interferogram.All the interferograms have been corrected using the error correct model ( 1), which would control error propagation into the final deformation signal.After that, deformation parameter and DEM error can be estimated at PSs and DSs based on the corrected interferograms using the proposed method.

Time-Series InSAR Results
The pixels with the mean coherence value larger than 0.7 are selected as the PSs.DSs candidates, such as non-cultivated land with short vegetation and bare land, are selected from the classification map generated by the ISODATA method.Then the AD test with a window of 11 × 11 is applied to refine the DSs that has an average coherence higher than 0.5 and number of SHPs higher than 70.A total of 498,836, of which 127,200 are PSs, are obtained from the SAR data.Based on the method described in Section 3.3, the mean deformation velocity in the line of sight (LOS) of the study area based on PSs and DSs is shown in Figure 7. Red triangle is selected as the reference point in the Zhengzhou city, which is stable based on the previous study.The deformation rate ranges from −55 mm/year to 10 mm/year at approximately 90% of the measured points.The standard deviation of the deformation velocity is 7.8 mm/year.Most of the study areas show a stable condition, in which the deformation rate is less than −5 mm/year.It is found that obvious ground motions have been detected in Zhengzhou city during the whole observation, with the maximum value reaches up to −55 mm/year.Figure 7b shows the estimated DEM errors, which are between the 40 m in the mountainous area to −20 m in the plain area, which are consistent with the relative accuracy of the SRTM DEM [32].

Time-Series InSAR Results
The pixels with the mean coherence value larger than 0.7 are selected as the PSs.DSs candidates, such as non-cultivated land with short vegetation and bare land, are selected from the classification map generated by the ISODATA method.Then the AD test with a window of 11 × 11 is applied to refine the DSs that has an average coherence higher than 0.5 and number of SHPs higher than 70.A total of 498,836, of which 127,200 are PSs, are obtained from the SAR data.Based on the method described in Section 3.3, the mean deformation velocity in the line of sight (LOS) of the study area based on PSs and DSs is shown in Figure 7. Red triangle is selected as the reference point in the Zhengzhou city, which is stable based on the previous study.The deformation rate ranges from −55 mm/year to 10 mm/year at approximately 90% of the measured points.The standard deviation of the deformation velocity is 7.8 mm/year.Most of the study areas show a stable condition, in which the deformation rate is less than −5 mm/year.It is found that obvious ground motions have been detected in Zhengzhou city during the whole observation, with the maximum value reaches up to −55 mm/year.Figure 7b shows the estimated DEM errors, which are between the 40 m in the mountainous area to −20 m in the plain area, which are consistent with the relative accuracy of the SRTM DEM [32].  Figure 8a shows the estimated deformation rate in Zhengzhou city (see red rectangle in Figure 6a).The histogram of the estimated deformation rate in Zhengzhou city is shown in Figure 8b.We can see that most of Zhengzhou city shows a movement away from the satellite with the maximum deformation rate up to −50 mm/year.It can be seen that the subsidence rate in the north is larger than that in other areas.Comparing the deformation rate and the distribution of the fault, no correlation has been found between them.However, previous study found that underground fault developed several fissures, which could cause different surface deformation [10,39].Figure 9a presents the time-series displacement in position P1 (Figure 8), which presents a linear trend during the observation period, with a mean velocity of −78.5 mm/year.Figure 9c display the time-series displacement in position P2, with a mean deformation velocity of −11.4 mm/year.In Figure 9b, it can be seen that, point P1 is near the building construction site.Figure 8a shows the estimated deformation rate in Zhengzhou city (see red rectangle in Figure 6a).The histogram of the estimated deformation rate in Zhengzhou city is shown in Figure 8b.We can see that most of Zhengzhou city shows a movement away from the satellite with the maximum deformation rate up to −50 mm/year.It can be seen that the subsidence rate in the north is larger than that in other areas.Comparing the deformation rate and the distribution of the fault, no correlation has been found between them.However, previous study found that underground fault developed several fissures, which could cause different surface deformation [10,39].Figure 9a presents the timeseries displacement in position P1 (Figure 8), which presents a linear trend during the observation period, with a mean velocity of −78.5 mm/year.Figure 9c display the time-series displacement in position P2, with a mean deformation velocity of −11.4 mm/year.In Figure 9b, it can be seen that, point P1 is near the building construction site.Figure 8a shows the estimated deformation rate in Zhengzhou city (see red rectangle in Figure 6a).The histogram of the estimated deformation rate in Zhengzhou city is shown in Figure 8b.We can see that most of Zhengzhou city shows a movement away from the satellite with the maximum deformation rate up to −50 mm/year.It can be seen that the subsidence rate in the north is larger than that in other areas.Comparing the deformation rate and the distribution of the fault, no correlation has been found between them.However, previous study found that underground fault developed several fissures, which could cause different surface deformation [10,39].Figure 9a presents the timeseries displacement in position P1 (Figure 8), which presents a linear trend during the observation period, with a mean velocity of −78.5 mm/year.Figure 9c display the time-series displacement in position P2, with a mean deformation velocity of −11.4 mm/year.In Figure 9b, it can be seen that, point P1 is near the building construction site.

Validation
It is difficult to get the levelling data to validate our result directly.Several attempts have been made to validate the InSAR-derived results in this paper indirectly.Firstly, the estimated DEM error has been obtained using the deformation retrieval model [19].The estimated DEM error value ranges from 40 m in the mountainous area to −20 m in the plain area, which is consistent with the relative accuracy of the SRTM DEM [32]; Secondly, the derived deformation results have been compared with the previous study [39].In [39], the authors used the Interferometric point target analysis (IPTA) technique to retrieve the ground deformation of Zhengzhou city from 2012 to 2013.The precision of deformation velocity in this case study was up to 3.16 mm/a based on the comparison with the ground leveling measurement.The deformation areas are mainly distributed in the northeast of Zhengzhou city.The deformation pattern and distribution in our paper is consistent with that of in [39].
The subsidence rate ranges from 10 to −50 mm/year, which is lower than that in 2013.Finally, we have compared the time-series amplitudes and optical images of the study.It is found that the deformation areas are mainly in the urban expansion areas.To validate the estimated deformation results in details, more supporting works, such as filed investigation and leveling data, should be taken.

Discussion
Through the above result in Section 4, it is found that the Zhengzhou city shows the most serious deformation.Due to the complex of the geological environment in the NCP and rapid expansion of urban area in Zhengzhou city, the reasons for the ground deformation are complex and variable.In this section, the possible factors causing the ground deformation are analyzed.

Geology and Subsidence
In this section, the effects of geomorphological features on deformation are analyzed.Zhengzhou city is located in the alluvial plain of the Yellow river.As shown in Figure 2, most of the geomorphologic features are in near northwest direction, which is similar to that of the faults.In Figure 8a, it can be seen that most of the deformation areas are in the fluvial plains.Figure 10 shows the histogram of the estimated deformation velocity in different geomorphologic areas.We can see that the mean deformation velocity in loess platform areas are less than −20 mm/year.While the mean deformation velocity in flood plain areas are larger than −20 mm/year.Those deformation rate difference are correlated with different geomorphologic types.It is much more obvious in Figure 11, which shows the subsidence rate along profile AA' (see Figure 8a).It can be seen that subsidence rate shows change along profile AA'.The fluvial plain area undergone the most serious subsidence with the maximum value of −55 mm/year.The deformation rate ranges from −20 to −30 mm/year in flood plain.The loess tableland has the lowest deformation rate, with the deformation rate varies from −5 to −20 mm/year.This difference reflects the soil property of each geomorphologic type.Based on previous studies, it was found that the shallow groundwater was widely distributed in the fluvial plain, which was the main water source for agriculture urban domestic uses.With the over exploitation of shallow groundwater, serious subsidence has been induced due to the large soil porosity in fluvial plain areas.For deeper analyzing the deformation mechanism of different geomorphologic types, more geological information and field investigation are required.

Urban Expansion and Subsidence
Zhengzhou city has experienced rapid urban expansion over the past decades.From 1988 to 2002, the urban area of Zhengzhou city increased from 73.73 km 2 to 159 km 2 , with the average annual growth rate of ~6% [40].The speed of the urban expansion does not slowdown in recent years.The north-eastern part of Zhengzhou city has the largest urban expansion.It is because this area is near the Yellow river, in which the water resource is reach and the relief is flatten.For the convenient transportation and construction, cities are prone to be situated in flat areas, particularly in the vicinity of rivers and coasts [41].
The Zhengzhou New District is located in the east of Zhengzhou city.Looking at the subsidence rate in Figure 8a, it is found that the serious subsidence areas are mainly distributed in the northern and eastern Zhengzhou city, which is consistent with the expansion direction of Zhengzhou city.Figure 12a give the deformation rate of the northern Zhengzhou city, with the maximum value reaches up to −50 mm/year.Figure 12b shows the time-series displacement of the three points (Figure

Urban Expansion and Subsidence
Zhengzhou city has experienced rapid urban expansion over the past decades.From 1988 to 2002, the urban area of Zhengzhou city increased from 73.73 km 2 to 159 km 2 , with the average annual growth rate of ~6% [40].The speed of the urban expansion does not slowdown in recent years.The north-eastern part of Zhengzhou city has the largest urban expansion.It is because this area is near the Yellow river, in which the water resource is reach and the relief is flatten.For the convenient transportation and construction, cities are prone to be situated in flat areas, particularly in the vicinity of rivers and coasts [41].
The Zhengzhou New District is located in the east of Zhengzhou city.Looking at the subsidence rate in Figure 8a, it is found that the serious subsidence areas are mainly distributed in the northern and eastern Zhengzhou city, which is consistent with the expansion direction of Zhengzhou city.Figure 12a give the deformation rate of the northern Zhengzhou city, with the maximum value reaches up to −50 mm/year.Figure 12b shows the time-series displacement of the three points (Figure 12a).It can be found that all the selected points undergone serious ground subsidence, with the deformation of −90 mm, −30 mm, −34 mm respectively during the whole observation period.It can be found that three points share the same linear deformation trend.From the amplitude images of the study area during the whole observation period, it is found that the area of the Zhengzhou New District increased.Figure 12c,d  Figure 13 shows the google earth maps of region I and II (Figure 12) in 10 March 2014 and 15 December 2016, corresponding to the start and end period of the whole observation.Comparing Figure 13a,b, it is found that most of the region A is covered by vegetation in 10 March 2014.After two and a half years, many buildings and roads have been constructed or under-construction.There was no water area in region A before October 2012.It is the under-construction of building that contributed the surface subsidence and the formation of subsidence lake.In Figure 12a, it is shown that in the west of the region II serious ground deformation has been detected.Comparing the google maps of region II in 10 March 2014 and 15 December 2016, it is obvious that the ground surface had undergone great change and serious subsidence caused subsidence lake in region II (see Figure 13d) during the whole observation period.be found that three points share the same linear deformation trend.From the amplitude images of the study area during the whole observation period, it is found that the area of the Zhengzhou New District increased.Figure 12c,d   Figure 13 shows the google earth maps of region I and II (Figure 12) in 10 March 2014 and 15 December 2016, corresponding to the start and end period of the whole observation.Comparing Figure 13a,b, it is found that most of the region A is covered by vegetation in 10 March 2014.After two and a half years, many buildings and roads have been constructed or under-construction.There was no water area in region A before October 2012.It is the under-construction of building that contributed the surface subsidence and the formation of subsidence lake.In Figure 12a, it is shown that in the west of the region II serious ground deformation has been detected.Comparing the google maps of region II in 10 March 2014 and 15 December 2016, it is obvious that the ground surface had undergone great change and serious subsidence caused subsidence lake in region II (see Figure 13d) during the whole observation period.

Ground Water and Subsidence
Zhengzhou is located on the impact plain of the Yellow River.Loose rock ground water is the most important and most widely distributed type of groundwater in Zhengzhou city, while clastic rock crevice water and carbonate karst water have the smallest distribution range in Zhengzhou, and their water-rich properties are poor.To analyze the relationship between the ground water and subsidence, it is necessary to obtain the count of ground water exploitation during the whole obversion.The ground water type and water-rich map of Zhengzhou city is collected and shown in Figure 14 [42].Comparing Figures 2 and 14, We can see that the distribution of ground water types

Ground Water and Subsidence
Zhengzhou is located on the impact plain of the Yellow River.Loose rock ground water is the most important and most widely distributed type of groundwater in Zhengzhou city, while clastic Remote Sens. 2018, 10, 1731 14 of 17 rock crevice water and carbonate karst water have the smallest distribution range in Zhengzhou, and their water-rich properties are poor.To analyze the relationship between the ground water and subsidence, it is necessary to obtain the count of ground water exploitation during the whole obversion.The ground water type and water-rich map of Zhengzhou city is collected and shown in Figure 14 [42].Comparing Figures 2 and 14, We can see that the distribution of ground water types is similar with that of geomorphologic types.Due to the difference of geomorphologic types, the aquifer lithology and water richness show obvious heterogeneity.As shown in Figure 14, the water outflow capacity ranges from over 3000 ton/day to less than 100 ton/day.Regions I and II (Figure 12) with serious subsidence are located the areas with the water outflow capacity ranges from over 1000 ton/day.Comparing the Figures 8a and 14, it is found that the deformation velocity is consistent with the water outflow capacity.It can be concluded that the ground water exploitation is one of the main factors causing the ground deformation in Zhengzhou city.

Conclusions
In Zhengzhou city, due to the underground resources exploitation and urban expansion, serious ground deformation has been observed in many areas.To measure the deformation and understand the characteristic of the subsidence of those areas, a time-series InSAR technique has been applied in Zhengzhou city.Considering the complex natural and geological environment, the proposed approach used in this study area attempts to overcome those limitations in two ways: (1) The orbit and topography-related phase errors have been corrected by a phase ramp correction method.(2) PSs and DSs are combined for deformation estimation, which improve the density of measurement points.Furthermore, the deformation parameters of PSs and DSs are retrieved based on a layered strategy from high phase quality to low.
The proposed approach was applied to 17 Radarsat-2 images acquired from 20 February 2014 to 31 October 2016.The general conclusions are as follows based on a quantitative analysis of the experimental results.
(1) The ground displacement in our study area could be retrieved by the proposed method

Conclusions
In Zhengzhou city, due to the underground resources exploitation and urban expansion, serious ground deformation has been observed in many areas.To measure the deformation and understand the characteristic of the subsidence of those areas, a time-series InSAR technique has been applied in Zhengzhou city.Considering the complex natural and geological environment, the proposed approach used in this study area attempts to overcome those limitations in two ways: (1) The orbit and topography-related phase errors have been corrected by a phase ramp correction method.(2) PSs and DSs are combined for deformation estimation, which improve the density of measurement points.Furthermore, the deformation parameters of PSs and DSs are retrieved based on a layered strategy from high phase quality to low.

Figure 1 .
Figure 1.(a) The study area is outlined by a red rectangle, background is the google map.The Radarsat-2 image coverage is represented by the yellow rectangle; (b) the mean amplitude map of the study area in radar geometry (red rectangular box in (a)).

Figure 2 .
Figure 2. Simplified geological map of Zhengzhou city.

Figure 1 .
Figure 1.(a) The study area is outlined by a red rectangle, background is the google map.The Radarsat-2 image coverage is represented by the yellow rectangle; (b) the mean amplitude map of the study area in radar geometry (red rectangular box in (a)).

Figure 1 .
Figure 1.(a) The study area is outlined by a red rectangle, background is the google map.The Radarsat-2 image coverage is represented by the yellow rectangle; (b) the mean amplitude map of the study area in radar geometry (red rectangular box in (a)).

Figure 2 .
Figure 2. Simplified geological map of Zhengzhou city.

Figure 2 .
Figure 2. Simplified geological map of Zhengzhou city.

Figure 3 .
Figure 3. Time-position of Radarsat-2 image interferometric pairs.Blue line represents the interferometric pairs.Red point is the master Synthetic Aperture Radar (SAR) image and black points are the slave images.

Figure 3 .
Figure 3. Time-position of Radarsat-2 image interferometric pairs.Blue line represents the interferometric pairs.Red point is the master Synthetic Aperture Radar (SAR) image and black points are the slave images.

Figure 4 .
Figure 4. Processing flowchart of time-series InSAR applied in this study.

Figure 4 .
Figure 4. Processing flowchart of time-series InSAR applied in this study.

19 Figure 5 .
Figure 5.A sample of the adaptive estimation for distributed scatterers (DSs) (a) the first estimation process with the initial search window; (b) the same process with a new search window; (c) the process iteratively by increasing the search window size until the pixel connecting enough POs.

Figure 5 .
Figure 5.A sample of the adaptive estimation for distributed scatterers (DSs) (a) the first estimation process with the initial search window; (b) the same process with a new search window; (c) the process iteratively by increasing the search window size until the pixel connecting enough POs.
Remote Sens. 2018, 8, x 9 of 19 correct model(1), which would control error propagation into the final deformation signal.After that, deformation parameter and DEM error can be estimated at PSs and DSs based on the corrected interferograms using the proposed method.

Figure 7 .
Figure 7. Deformation velocity in line of sight (a) and DEM error (b) estimated by both PSs and DSs.Red triangle is the location of the reference point.Background map is the DEM map of the study area.

Figure 7 .
Figure 7. Deformation velocity in line of sight (a) and DEM error (b) estimated by both PSs and DSs.Red triangle is the location of the reference point.Background map is the DEM map of the study area.

Figure 8 .
Figure 8.(a) Deformation velocity in Zhengzhou area.(b) The histogram of the estimated deformation velocity.

Figure 9 .
Figure 9. (a,c) displacement time-series of points P1 and P2, respectively; (b,d) are the corresponding optical images.

Figure 8 .
Figure 8.(a) Deformation velocity in Zhengzhou area.(b) The histogram of the estimated deformation velocity.

Figure 8 .
Figure 8.(a) Deformation velocity in Zhengzhou area.(b) The histogram of the estimated deformation velocity.

Figure 9 .
Figure 9. (a,c) displacement time-series of points P1 and P2, respectively; (b,d) are the corresponding optical images.

Figure 9 .
Figure 9. (a,c) displacement time-series of points P1 and P2, respectively; (b,d) are the corresponding optical images.
displays the amplitudes of the northern Zhengzhou city in 20 February 2014 and 31 October 2016, respectively.It can be seen that those areas have changed greatly due to city construction, especially in regions I and II.And the maximum deformation rate reaches up to −20 mm/year during the whole observation.
displays the amplitudes of the northern Zhengzhou city in 20 February 2014 and 31 October 2016, respectively.It can be seen that those areas have changed greatly due to city construction, especially in regions I and II.And the maximum deformation rate reaches up to −20 mm/year during the whole observation.

Figure 12 .
Figure 12.(a) the deformation rate of the northern Zhengzhou city; (b) the time-series displacement of the selected points in Figure12a; The amplitude maps of northern Zhengzhou city in 20 February 2014 (c) and 31 October 2016 (d).Regions I and II changed greatly during the whole observation period.

Figure 12 .
Figure 12.(a) the deformation rate of the northern Zhengzhou city; (b) the time-series displacement of the selected points in Figure 12a; The amplitude maps of northern Zhengzhou city in 20 February 2014 (c) and 31 October 2016 (d).Regions I and II changed greatly during the whole observation period.Remote Sens. 2018, 8, x 15 of 19

Figure 13 .
Figure 13.(a,b) Google earth maps of region I and II in 10 March 2014, respectively.(c,d) Google earth maps of region I and II in 15 December 2016, respectively.

Figure 13 .
Figure 13.(a,b) Google earth maps of region I and II in 10 March 2014, respectively.(c,d) Google earth maps of region I and II in 15 December 2016, respectively.

Figure 14 .
Figure 14.Ground water type and water-rich map of Zhengzhou city.Region I and II are the urban expansion area (Figure 12).

Figure 14 .
Figure 14.Ground water type and water-rich map of Zhengzhou city.Region I and II are the urban expansion area (Figure 12).