Spatiotemporal Characterization of Land Subsidence in Guandu (China) Revealed by Multisensor InSAR Observations

Excessive groundwater exploitation has brought about severe ground subsidence in Guandu (China), threatening the stability of urban infrastructure. Mapping of the spatiotemporal variations of ground deformation is urgently needed for disaster prevention and mitigation. In this study, multisensor interferometric synthetic aperture radar (InSAR) observations were applied to Guandu to derive the time series deformation from 2007 to 2019. The annual deformation velocity revealed three severe subsiding regions in Guandu. Based on the ascending and descending Sentinel-1 images with overlapping temporal and spatial coverage, two-dimensional vertical and horizontal east–west deformation was calculated and indicated that the deformation in Guandu was dominated by vertical direction. After connecting the multisensor results, long-term ground deformation spanning from January 9, 2007, to September 1, 2019, was produced and showed that the north subsiding region experienced fast followed by slow subsidence, whereas the south subsiding region experienced slow followed by fast subsidence. This difference was due to the changes of groundwater pumping centers and rates. The cumulative maximum subsidence reached 400mm during the period of 2007–2019. The similar variations in temporal domain between the change of groundwater level and ground deformation suggested that groundwater exploitation accounted for the severe subsidence in Guandu. Our results may provide scientific evidence regarding the sound management of groundwater exploitation to mitigate potential damage to infrastructure and the environment.


Introduction
Guandu is located at a longitude of 102°38 ′ -103°03 ′ east and a latitude of 24°54 ′ -25°17 ′ north [1]. Situated in the mideastern Yunnan-Guizhou plateau in Southwest China, Guandu is a typical lacustrine sediment basin with altitude from 1886 to 2731 m [2]. Since the launch of China's reform program in 1978, the population and economics of Guandu have experienced rapid development. Accordingly, the demands of water for domestic and industrial consumption increased rapidly over the past forty years. In this context, phenomenon of groundwater exploitation was serious in Guandu. The related document showed that groundwater level of Guandu had declined about 40.3 m between 1984 and 1998 (InSAR) has demonstrated its potential for high-density spatial mapping of ground deformation associated with earthquakes [10], volcanoes [11], and other geologic processes [12]. Recent advanced InSAR techniques have improved our understanding of the process of ground deformation, such as the deep learning approach [13], artificial intelligence technique [14], optimal phase unwrapping algorithm [15], signal retrieval for decorrelating targets [16], four-dimensional filtering approach [17], improved synthetic aperture radar (SAR) image coregistration algorithm [18], atmospheric delay correction method [19], and coherent point selection algorithm [20]. Additionally, an increasing amount of SAR satellites have provided a large set of multisensor SAR images, which have been employed to reconstruct the spatiotemporal evolution of ground deformation [21]. These advanced techniques have facilitated land subsidence monitoring in many urban areas, such as in Beijing [22], Xian [23], Shenzhen [24], Taiyuan [25], Mexico [26], Italy [27], and Turkey [28]. Regarding Guandu, existing research focused on the whole city and there was no information available for spatiotemporal characterization analysis of ground deformation so far [7][8][9].
To better understand the spatiotemporal evolution of ground deformation in Guandu, the multisensor InSAR observation was utilized to obtain the long-term and twodimensional time series deformation. For this, 20 ALOS-1 images, 40 COSMO-SkyMed (CSK) images, 24 ALOS-2 images, and 91 Sentinel-1A images were collected in this study. The annual deformation velocity was estimated to characterize the spatial pattern of subsidence in Guandu. Subsequently, two-dimensional and long-term time series deformation was obtained to analyze the temporal evolution of subsidence. Finally, the correlation between subsidence and groundwater changes was discussed.

Study Area and Datasets
2.1. History and Geologic Setting of Study Area. Guandu, situated in the eastern part of Kunming (China) (Figure 1(a)), has one of the mildest climates in China, characterized by short, cool, dry winters with mild days and crisp nights and long, warm, humid summers [1]. The annual average precipitation is about 1,002 mm, of which 86%-90% concentrates in the summer and autumn (May to October). The annual average evaporation is about 1,900 mm, which is greater than the annual average precipitation [3]. Groundwater is the main water supply in Guandu, where the annual average groundwater exploitation is more than 20,000 m 3 , among which 90% is used as domestic water consumption and 10% is used as industrial water consumption [8]. Due to the pollution of shallow pore water, exploitation of groundwater is mainly from the karst water in the study area. Figure 1(b) shows the stratigraphic profile along the line PP′ in Figure 1(a), where clay, sand, and dolomite are widespread in the study area. It is observed from Figure 1(b) that the contact layer between soil and bedrock is not the water storage layer, indicating that the soil aquifer is directly connected to the bedrock aquifer. In this context, the groundwater level of soil aquifer rapidly declines with the decline of groundwater level of bedrock aquifer. This decline causes the water-release compression and soil consolidation, resulting in the ground subsidence. Thus, groundwater exploitation of karst water leads to the several subsidence over the study area. Early studies showed that groundwater level of Guandu had declined about 40.3 m between 1984 and 1998, which caused the cumulative subsidence of 250 mm [8], as shown in Figure 1(c).
2.2. SAR Datasets. 175 SAR images in total, including 20 ALOS-1 images, 40 COSMO-SkyMed (CSK) images, 24 ALOS-2 images, and 91 Sentinel-1A images, were collected to derive the deformation time series over the study area. Table 1 and Figure 1(a) show the parameters and coverage of these SAR images, respectively. The collected multisensor SAR datasets with different imaging parameters, e.g., azimuth and incidence angle, spatial and temporal resolution, orbit direction, and wavelength, allow us to comprehensively analyze the observed surface deformation [29]. Meanwhile, two-dimensional ground deformation was calculated by integrating of the multisensor SAR datasets with overlapping temporal and spatial coverage [30]. After setting the spatial and temporal baseline thresholds, a total of 624 interferograms was generated from the collected multisensor SAR images. Based on the generated interferograms, the ground deformation between January 2007 and September 2019 was estimated over the study area. Light detection and ranging (LIDAR) digital elevation model (DEM), which had a spatial resolution of 3 m and centimeter-level height precision, was acquired as an external DEM to remove the topographic phase from the differential interferograms.

Methodology
Using the collected 175 multisensor SAR images, ground deformation was estimated by multitemporal SAR interferometry technique. The interferometric pairs were generated by setting small temporal and spatial baselines, which limited the noise effects and preserved the temporal and spatial coherence characteristics of the interferograms. As a result, totally 624 interferometric pairs were produced, among which 59 pairs were from ALOS-1 images, 88 pairs were from COSMO-SkyMed images, 71 and 15 pairs were, respectively, from ALOS-2 ascending and descending images, and 91 and 300 pairs were, respectively, from Sentinel-1 ascending and descending images. The topographic and orbital phases were simulated and subsequently removed from the differential interferometric phase by using the collected LIDAR DEM and SAR orbital data. Then, the process of adaptive spectral filtering with a window size of 32 was executed to suppress the interferometric noise [31]. After that, phase unwrapping using the Minimum Cost Flow (MCF) method was carried out to retrieve the absolute phase [32]. Finally, the atmospheric delay errors were reduced through the cascade of a lowpass filtering implemented in the two-dimensional spatial domain, followed by a temporal high-pass filtering [33]. After these processes, the unwrapped interferograms from multisensor datasets were resampled to a common study area (red rectangle in Figure1(a)) and used to estimate 3.1. Vertical Deformation Velocities. The vertical deformation velocities of multisensor SAR datasets with different imaging parameters were calculated to characterize the spatial pattern of ground deformation over the study area. Using the produced unwrapped interferometric pairs, the vertical deformation velocities were estimated by [34] where Ph rate was the deformation velocities, N was the number of interferograms involved in the estimation, Δt j was the time interval between master and slave pair, and φ j was the unwrapped interferograms. Meanwhile, the standard deviation of Ph rate was estimated by [35]

Journal of Sensors
where l was the radar wavelength and VarðPh rate Þ was the standard deviation of Ph rate . In this case, six vertical deformation velocities from ALOS-1, COSMO-SkyMed, ascending and descending ALOS-2, and ascending and descending Sentinel-1 datasets were produced over the study area, which allowed us to compare the observed deformation.

Two-Dimensional Ground Deformation.
In areas with overlapping spatiotemporal coverage, it is possible to map the two-dimensional vertical and horizontal east-west ground deformation through combing of multisensor SAR datasets. In this case, the Sentinel-1 datasets in ascending and descending orbits were used to derive the twodimensional ground deformation. It is worth noting that the ALOS-2 datasets were not involved in the twodimensional deformation estimation due to the limited acquisitions in descending orbit. The multidimensional small baseline subset (MSBAS) method was adopted to calculate   Journal of Sensors the two-dimensional ground deformation [36]: where matrix A consisted of time intervals between consecutive SAR acquisition, θ was the azimuth angle, ϕ was the inci-dence angles; V E and V U represented unknown horizontal east-west and vertical velocities that were to be determined, b Φ represented the unwrapped, geocoded, and resampled interferometric phase, λ was a regularization parameters, and L was a zero-, first-, or second-order difference operator. The unknown parameters V E and V U for each pixel were solved by applying the singular value decomposition (SVD), and deformation time series were reconstructed from the computed deformation rates by numerical integration. In this case, the first order regularization with λ equal to 0.

Vertical Deformation Velocities with Different SAR
Datasets. Figure 2 shows the vertical deformation velocities, which were derived from ALOS-1, ALOS-2, and Sentinel-1 datasets with ascending orbit and COSMO-SkyMed, ALOS-2, and Sentinel-1 datasets with descending orbit. It is worth noting that all measurement values were relative to a common reference point (red cross in Figure 2) that was considered to be stable over the entire time period. Positive and negative values represented uplift and subsidence. Three prominent deformation regions were observed from Figure 2. The first deformation region named as Guanshang was located at the center of deformation maps and presented the shape of approximate pear. The second deformation region with the shape of approximate oval was named as Yangfeng and located at the lower center of deformation maps, which was considered as the extension of Guanshang deformation region. The third deformation region was located at the lower left of deformation maps and named as Harbor. For these three deformation regions, it was found that they showed the different deformation magnitudes with respect to the different datasets. The Guanshang deformation region displayed the severe subsidence in ALOS-1 and COSMO-SkyMed datasets, while it was not prominent in ALOS-2 and Sentinel-1 datasets. As for the Yangfeng deformation region, the severe subsidence was observed in COSMO-SkyMed, ALOS-2, and Sentinel-1 datasets, while it was not prominent in ALOS-1 datasets. Compared to the Guanshang deformation region, the Harbor deformation region presented the opposite deformation magnitude: severe subsidence appeared in ALOS-2 and Sentinel-1 datasets, while insignificant deformation was observed in ALOS-1 and COSMO-SkyMed datasets. This phenomenon was due to the different imaging time regarding to the different SAR datasets. As shown in the introduction, the subsidence in Guandu was highly correlated with groundwater extraction, meaning that the degree of groundwater extraction with respect to the different times may have led to the different deformation magnitudes. The relevant information was investigated, and the degree of groundwater extraction was found to be extremely intensive in Guanshang in early acquired ALOS-1 and COSMO-SkyMed datasets. Therefore, the Guanshang deformation region showed the significant signal in ALOS-1 and COSMO-SkyMed datasets, while it was not significant in ALOS-2 and Sentinel-1 datasets. The Yangfeng and Harbor deformation regions experienced the 8 Journal of Sensors slow followed by fast subsidence in temporal domain and therefore showed the significant signal in later acquired ALOS-2 and Sentinel-1 datasets. Further inspection indicated that there were subtle differences in deformation shapes between descending and ascending orbits, particularly for the Yangfeng deformation region, where the deformation shape presented in descending orbit was elongated in ascending orbit. Through analyzing the vertical deformation velocities from different SAR datasets, it was summarized that three prominent subsiding regions were observed over the study area.

Two-Dimensional Vertical and Horizontal Deformation.
The difference in deformation shape between descending and ascending orbits indicated the existence of horizontal deformation. In order to investigate this deformation, the deformation velocities and time series deformation maps along the vertical and horizontal east-west directions were derived from the ascending and descending Sentinel-1 images, which were acquired from May 23, 2015, to September 16, 2016. The positive and negative values in vertical direction, respectively, represent uplift and subsidence, while positive and negative values in east-west direction, respectively, represent eastward and westward movement. The vertical and eastwest deformation velocities are presented in Figure 3. Compared to Figure 2, the Yangfeng and Harbor deformation regions were retained in Figure 3, while the deformation in Guanshang region was not obvious. As analyzed in Section 9 Journal of Sensors 4.1, this phenomenon was due to the temporal variations of deformation. The statistics showed that the maximum subsidence and eastward movement rates were 38 mm/year and 11 mm/year for Yangfeng region, respectively. For Harbor deformation region, the maximum subsidence rate was up to 25 mm/year, and the maximum eastward movement rate was 12 mm/year. Compared with the Yangfeng deformation region, the phenomenon of decorrelation was relatively severe  Journal of Sensors at Harbor deformation region. This decorrelation was due to the ongoing constructions according to the filed investigation. The cause of observed deformation in Figure 3 was attributed to the groundwater exploitation considering the history and geologic setting of study area. The pioneering researches suggested that this aquifer deformation occurred in both the vertical and horizontal directions in response to fluid withdrawal from confined aquifer system [37][38][39]. In the vertical direction, the decline of groundwater level led to the water-release compression and soil consolidation with the result of ground subsidence. In the horizontal direction, the eastward movement was related with a fault located at the east of the study area (as shown in Figure 1(b)), which impeded the horizontal propagation of fluid-pressure changes. The magnitude of aquifer deformation was a function of the distance from the center of pumping, the pumping rate, the time since the initiation of pumping, the hydraulic diffusivity of the aquifer, and the presence of heterogeneous barriers or boundaries that may impede flow and/or displacement [39]. Unfortunately, it was difficult to evaluate the observed vertical and horizontal deformation because of the lack of relevant supporting data. However, the standard deviations of two-dimensional deformation velocities were below 3 mm for most regions, suggesting the reliability of our observed deformation.
To characterize the spatiotemporal variation of deformation over the study area, the deformation rates along the profile line AA ′ and time series deformation at point P1, both of which are marked in Figure 3, are displayed in Figure 4. It was observed in Figure 4(a) that the vertical deformation presented the shape of a funnel, while east-west deformation presents the shape of approximate parabola. As for the temporal variations in Figure 4

Long-Term Vertical Deformation Time Series.
To investigate the temporal variation of ground deformation over the study area, the collected multisensor SAR datasets were combined to produce the long-term vertical deformation time series. During this procedure, the line-of-slight (LOS) deformation was directly converted to the vertical direction since it was dominated by the vertical deformation in this case, which has been proved in the last section. The final time series deformation was from January 9, 2007, to September 1, 2019, as shown in Figure 5. For saving the page space, we only showed the results in September or October from 2007 to 2019, as shown in Figure 6. The prominent subsiding region was detected, where the cumulative maximum subsidence reached 400 mm. Further observations suggested that the subsidence experienced two different stages of development. In the first stage spanning from 2007 to 2014, the rapid subsidence was observed at the north subsiding region, while the subsidence was relative slow for the south subsiding region. However, this phenomenon was changed in the second stage, which spans from 2014 to 2019. The accelerated subsidence appeared at the south subsiding region, while the decelerated subsidence was observed at the north subsid-ing region. We think this phenomenon may be related with the changes of groundwater pumping centers and rates. Figure 7 displays the time series deformation along the line BB ′ marked in Figure 6. It was found that the subsiding region was composed of two different subsiding funnels, where the cumulative maximum subsidence was 330 mm and 370 mm, respectively. Meanwhile, two subsiding funnels presented the different developing processes: the north region showed the fast followed by slow subsidence, while the south region showed the slow followed by fast subsidence.

Correlation between Ground Deformation and Groundwater Exploration
Groundwater is the main water supply in Guandu, where the annual average groundwater exploitation is more than 20,000 m 3 . To investigate the changes of groundwater, the groundwater observations at Guanshang and Yangfeng stations were collected in this study, as shown by the blue rectangles in Figure 8. It is observed in Figure 8(a) that there is a rapid decline between 2008 and 2010 followed by a slow decline between 2010 and 2014 for Guanshang groundwater level. It is different with Guanshang that Yangfeng station (Figure 8(b)) presents a gradual decline between 2008 and 2012 followed by a rapid decline between 2012 and 2014. The cumulative declines of groundwater at Guanshang and Yangfeng stations between 2008 and 2014 are, respectively, 12.2 and 13.8 m. Considering the presented stratigraphic structure in Figure 1(b), the subsidence easily occurred in the presence of groundwater exploitation, which caused the water-release compression and soil consolidation. The red circles in Figure 8 present the similar variations in temporal domain, particularly for Yangfeng station. However, there is a subtle inconsistent variation between the change of groundwater level and ground deformation from 2008 to 2011 for Guanshang station; the groundwater level illustrates the accelerated decline while ground deformation shows the steady decline. This inconsistency may have been caused by the inaccurate groundwater observation in 2010 since Guanshang groundwater gauging station was temporarily damaged during the building construction between 2010 and 2011 [7]. Based on the analyses of ground deformation and groundwater changes, groundwater exploitation accounts for the severe subsidence in Guandu.

Conclusion
With the rapid development of population and economics in Guandu, ground subsidence has caused substantial damage to homes, roads, canals, pipelines, and other types of infrastructure. In this context, time series deformation of Guandu was derived from multisensor InSAR observations to characterize and monitor its spatiotemporal variations. This observation may provide scientific evidence regarding the sound management of groundwater exploitation. Based on this study, the following conclusions were made.
(1) Severe subsiding regions were detected through multisensor InSAR observations. Using 20 ALOS-1 11 Journal of Sensors images, 40 COSMO-SkyMed (CSK) images, 24 ALOS-2 images, and 91 Sentinel-1A images, the annual deformation velocity and time series deformation maps of Guandu were retrieved through multisensor InSAR processing. The results showed that three severe subsiding regions were observed, where the cumulative subsidence reached 400 mm during the period of 2007-2019. The deformation of Guandu was dominated by the vertical direction through analyzing the two-dimensional time series deformation (2) Groundwater exploitation accounted for the subsidence in Guandu. Comparison between the change of groundwater level and ground deformation indicated that they presented the similar variations in temporal domain, suggesting that groundwater exploitation caused the severe subsidence in Guandu

Data Availability
Sentinel-1A/B images were from the European Space Agency (ESA) project, and these data were downloaded from the Sentinel-1 Scientific Data Hub. ALOS-2 images were provided by the Japan Aerospace Exploration Agency (JAXA). Several figures were prepared using Generic Mapping Tools software.

Conflicts of Interest
The authors declare no conflicts of interest.