Quantitative Assessment and Impact Analysis of Land Surface Deformation in Wuxi Based on PS-InSAR and GARCH Model

1. Introduction

Ground deformation refers to the fluctuating trend of surface morphology within a certain spatiotemporal range, including subsidence, uplift, and horizontal movement. Among them, subsidence and uplift have significant impacts on human life and the natural environment [1,2,3]. Ground subsidence is a geological phenomenon in which the altitude of the Earth’s surface continuously decreases over a certain period under the influence of anthropogenic or natural factors [4]. In contrast, ground uplift refers to the phenomenon of upwarping or elevation of a certain area on the Earth’s surface due to factors such as crustal movement or changes in groundwater level [3]. Whether subsidence or uplift, both can lead to reduced surface stability and trigger severe geological disasters. For instance, Texas’s Katy area witnessed an annual subsidence rate of 14 mm, forming a distinct circular pattern that heightens flood and seawater intrusion risks. [5]. In Indonesia’s Bandung Basin, subsidence impacted several towns and industrial areas within the central basin, with rates as high as 187 mm/year. This led to extensive damage to buildings and infrastructure, increasing safety risks [6]. This ongoing geological movement continues to afflict neighboring regions with natural disasters, including earthquakes and landslides [7].

Subsidence can significantly diminish the water retention ability of subterranean aquifers, causing ground collapse, damaging structures and infrastructures, and potentially triggering natural calamities like floods [8]. Conversely, uplift not only alters the flow direction of water sources, affecting water flow and velocity, but also leads to land desiccation and vegetation reduction [9]. Moreover, ground deformation can change the stress and strain state of the Earth’s crust, exacerbating surface instability trends and even triggering earthquakes [10]. According to projections by the United Nations Hydrological Organization, by 2040, ground deformation will impact over 1.6 billion people globally [11]. Therefore, it is imperative to conduct dynamic and high-precision monitoring of surface elevation changes.

Ground deformation monitoring traditionally relies on techniques like geodetic leveling and GNSS measurements [12,13]. GNSS technology offers high precision, robust real-time performance, wide coverage, and widespread application in fields such as three-dimensional terrain reconstruction, aircraft navigation, and mechanical automation guidance [14]. However, the observation spatial density of GNSS methods is relatively low, which affects the investigation and detection of subsidence funnels. Although traditional geodetic leveling techniques have high accuracy [15], the labor cost of surveying is high, and the measurement range is limited, which are significant drawbacks [16,17]. In the late 1990s, Interferometric Synthetic Aperture Radar (InSAR) technology was developed based on Synthetic Aperture Radar (SAR) and is a new type of space-to-earth observation technology [18]. InSAR utilizes two or more SAR images acquired at different times and calculates the slight deformation of the ground surface by comparing the phase differences of radar echoes at the same location [19,20]. This technology provides high-precision three-dimensional ground deformation data and has the advantages of wide coverage, immunity to climate influence, and real-time monitoring [21]. However, InSAR technology has certain requirements for interferometric conditions and is subject to various error sources such as instrumental noise, orbital errors, and atmospheric delays, resulting in limited observation accuracy. In the early 21st century, time-series InSAR (TS-InSAR) technology was further developed based on InSAR technology [22]. TS-InSAR mainly obtains high-precision time-varying signals by modeling multiple interferograms and weakening the influence of errors in traditional methods. TS-InSAR technology is mainly divided into two categories: the first category is represented by Persistent Scatterer InSAR (PS-InSAR), which identifies permanent scatterers with high reflection intensity and good coherence over a long period in SAR images, calculates the deformation of scatterers using relevant techniques, and ultimately obtains ground deformation information of the study area; the second category is the multi-master image time-series InSAR technology represented by Small Baseline Subset InSAR (SBAS-InSAR). Additionally, there are also the Coherence Target (CT) method and the Interferometric Point Target Analysis (IPTA) method derived from the above two methods [23]. In recent years, TS-InSAR technology has been widely applied. For example, Zhou et al. used SBAS-InSAR technology to analyze the spatial and temporal changes of land subsidence in Wuhan from 2017 to 2021 and found that land subsidence in Wuhan was uneven and concentrated in a few areas [24]. Meer Muhammad Sajjad et al. used PS-InSAR technology to assess the impact of land subsidence caused by groundwater depletion and aquifer degradation in Lahore, Pakistan, and proposed a theory of sustainable urban development based on this [25].

Wuxi City, situated in the Su-Xi-Chang subsidence zone and the Yangtze River Delta alluvial plain [26], experiences significant surface deformation, necessitating surveys in the area [27]. In recent years, numerous scholars have investigated this region. Yang et al., for instance, employed PS-InSAR technology and three-dimensional continuous medium simulation to analyze ground subsidence and ground fissure activities in the Su-Xi-Chang area of China, confirming that prohibiting groundwater extraction is an effective measure to prevent and control these disasters [28]. Similarly, Wang et al. utilized SBAS-InSAR technology to monitor ground subsidence in Wuxi City from 2007 to 2011, concluding that the subsidence might be related to large-scale land consolidation projects [29].

A literature review reveals that current research focuses primarily on the Su-Xi-Chang area and the Yangtze River Delta region, with fewer studies conducting long time series and large-scale investigations in the Wuxi area [30]. To examine surface deformation trends and internal mechanisms in the Wuxi area in recent years, this study employs PS-InSAR technology to process 100 Sentinel-1A SAR images from November 2015 to June 2023. The analysis incorporates multiple perspectives, including groundwater, precipitation, and soil geology. Furthermore, this study innovatively applies the GARCH model to analyze subsidence displacement in Wuxi City for the first time and proposes a novel “amplitude factor” index to evaluate surface fluctuations.

2. Study Area and Datasets

2.1. Study Area

Wuxi City is situated in the southeast of Jiangsu Province, with geographical coordinates of 31°22′N—32°02′N and 119°59′E—120°38′E. It is a central city in the Yangtze River Delta economic belt of China, with a population of over 7.46 million. The study area of this paper includes the main urban districts of Wuxi City (Xinwu District, Xishan District, Huishan District, Binhu District, and Liangxi District) and Jiangyin City (Figure 1). The study area borders Suzhou to the east, adjoins Taihu Lake to the south, and is adjacent to Zhejiang Province; it borders Changzhou to the west, relies on the Yangtze River to the north, and is separated from Jingjiang City of Taizhou City by the river [31]. The total area of the study area is 2,574.71 square kilometers, with plains as the main landform, supplemented by low mountains and remnant hills. The overall topography exhibits a distribution pattern of high in the southwest and north and low in the central part, with an average elevation of approximately 8 meters. Wuxi City belongs to the northern subtropical monsoon climate, located between the Yangtze River and Taihu Lake, with a humid climate and rain in hot season. According to the records of the Jiangsu Provincial Bureau of Statistics, the average annual precipitation in Wuxi City is 1,142 mm, concentrated in the plum rain season from June to August. There are approximately 5,635 rivers within the territory of Wuxi City, with a total area of 939.61 square kilometers, accounting for 20.3% of the city’s total area. These rivers, together with the Yangtze River, the Beijing-Hangzhou Grand Canal, and Taihu Lake, form a vast water network [32]. Within the study area, almost all regions are covered by Quaternary loose sediments, mainly composed of silty clay and sub-sand. These two soil layers are interbedded, forming the soft soil geology of the Wuxi area, which has a significant impact on the surface deformation in the region [28,30].

2.2. Data Sets

This study utilizes 100 Sentinel-1A images provided by the European Space Agency (ESA), spanning from November 2015 to June 2023. These images cover the main urban areas of Wuxi City and Jiangyin City, with the absolute difference of acquisition time intervals not exceeding 2 months. The images were acquired using the Interferometric Wide Swath (IW) mode and VV polarization. Detailed information about the SAR data is presented in Table 1. Additionally, this study employed the 30-meter Shuttle Radar Topography Mission (SRTM) Digital Elevation Model (DEM) from the National Aeronautics and Space Administration (NASA) to remove topographic phase, and the Precise Orbit Data (POD) released by ESA to refine orbital information, thereby improving the accuracy of georeference and baseline estimation.

This study collected various data related to surface changes in the Wuxi area. Groundwater data were obtained from the Jiangsu Bureau of Geological Survey. We have selected six evenly distributed monitoring wells and obtained water level data corresponding to these wells to study the relationship between surface deformation and pore confined water. The China Meteorological Administration provided precipitation data, and the Chinese Academy of Geological Sciences provided geological data. This study utilized geological borehole data to determine the depth of Quaternary sediments in the Wuxi area.

3. Methodology

3.1. PS-InSAR Technique

Ferretti et al. found that during the InSAR data processing, certain ground pixels (such as buildings and exposed rocks) can maintain high coherence in the images over a longer time scale between two SAR satellite acquisitions [33]. If these discrete and sparse PS points with high coherence can be selected from the SAR images, and interferometric phase modeling can be performed to solve for surface deformation information, deformation time series, and residual elevation, the influence of incoherent pixels in SAR images can be effectively avoided, overcoming the temporal and spatial decorrelation effects of traditional differential interferometric synthetic aperture radar (DInSAR) techniques, thus obtaining reliable deformation information with high accuracy. Based on this theory, Ferretti et al. first proposed the concept of “Permanent Scatterer (PS)” and the PS interferometry technique. This technique utilizes multiple image data, extracts PS points based on the amplitude or coherence coefficient method, and then establishes a parametric equation that includes average deformation rate, atmospheric delay phase, orbital residual error, and elevation error [34]. Through iterative calculations, the time-series information of surface deformation in the study area can be obtained, and the solution accuracy of this method can reach the millimeter level [35].

The processing flow of PS-InSAR is as follows: From N complex SAR images of the same scene, one image is selected as the master image, and the rest are slave images; the slave images are georeferenced with the master image, obtaining N-1 georeferenced SAR images; PS points are selected based on a set amplitude dispersion index threshold, and M interferometric pairs are generated [23]. The interferometric phase of each pixel in the interferometric pairs consists of multiple components, including:

$$\phi ={\phi }_0+{\phi }_{topo}+{\phi }_{defo}+{\phi }_{atm}+{\phi }_{noise}$$

(1)

where $\phi$ is the total interferometric phase, ${\phi }_{defo}$ is the deformation phase, ${\phi }_0$ is the phase on the reference ellipsoid, ${\phi }_{topo}$ is the topographic phase, ${\phi }_{atm}$ is the phase caused by changes in atmospheric conditions during the two radar acquisitions, and ${\phi }_{noise}$ is the random noise phase [22].

Using the known digital elevation model (DEM) data, differential interferometric processing is performed on the M interferometric pairs to obtain M differential interferograms. By establishing a phase difference model between adjacent PS points, the surface subsidence information is calculated. The phase difference between adjacent PS points in the i-th differential interferogram can be expressed as [36]:

$${∆\phi }_i=\frac{4\pi }{\lambda }{T}_i∆v+\frac{4\pi B_i^{\perp }∆\epsilon }{\lambda R\sin \theta }+∆{\phi }_i^{res}$$

(2)

where, $∆v$ is the deformation velocity difference between two adjacent PS points; $∆\epsilon$ is the phase difference caused by the elevation error; $B_i^{\perp }$ and $T_i$ are the spatial baseline and temporal baseline of the i-th interferogram, respectively; $\lambda$ is the radar wavelength; R is the distance from the sensor to the target; $\theta$ is the radar incidence angle; $∆{\phi }_i^{res}$ is the residual phase, including atmospheric delay phase, nonlinear deformation phase, and noise phase.

Ferretti et al.’s research showed that under the condition that the absolute value of the residual phase is less than $\pi$, the deformation velocity difference $∆v$ and the elevation error $∆\epsilon$ can be estimated by maximizing the phase coherence coefficient model of M interferograms, with the objective function as follows [37]:

$$\gamma =\left|\frac{1}{M}{\sum }_{i=1}^M\left(\cos ∆wi+j·\sin ∆wi\right)\right|=maximum$$

(3)

where, $\gamma$ is the model coherence coefficient; $j$ is the imaginary unit, i.e., $j=\sqrt{-1}$; $∆wi$ represents the difference between the observed phase and the fitted phase of the i-th interferogram, namely:

$$∆wi=∆{\phi }_i-\frac{4\pi }{\lambda ·R·\sin \vartheta }·B_i^{\perp }·∆\epsilon -\frac{4\pi }{\lambda }·T_i·∆v$$

(4)

Although the objective function is highly nonlinear and the observed phase exhibits wrapping phenomena, by solving for the maximum value of the model coherence coefficient, the unknowns $∆w$ and $∆v$ for each baseline can be determined within a predefined solution space, avoiding the unwrapping process. Using $\gamma$ as weights and applying the weighted least squares method, the time-series deformation velocity $v$ and elevation error $\epsilon$ of PS points on a sparse grid are solved from known reference points. By removing the linear deformation component of PS points, the nonlinear deformation component, atmospheric delay, and decorrelation noise are obtained. Based on the spatial and temporal frequency characteristics of each component, filtering techniques are employed to decompose the signals. Among them, the atmospheric delay appears as a low-frequency signal in the spatial domain and a high-frequency signal in the temporal domain, exhibiting strong instability. By applying a high-pass filter to the residual phase in the temporal domain and a low-pass filter in the spatial domain, the atmospheric delay phase of each image relative to the master image is obtained [36].

3.2. PS-InSAR Data Processing

Considering the large amount of SAR data, long time span, and wide geographical coverage used in this study, we divided the 100 images covering the Wuxi area into two groups for processing. One group is used for detailed year-by-year analysis, while the other group is used for long-term trend analysis, with the aim of better understanding the surface deformation patterns in the Wuxi area. We divided the period from November 2015 to June 2023 into three time intervals, and Figure 2 shows the temporal baseline distribution and time span of these three intervals. The specific processing flow is shown in Figure 3, which mainly includes five parts: data collection and preprocessing, differential interferometry calculation, first-step inversion, second-step inversion, and geocoding.

1. The data collection and processing part mainly includes downloading 100 Single Look Complex (SLC) images from the Alaska Satellite Facility (ASF) website and 100 precise orbit data from the European Space Agency (ESA) that are georeferenced with the images. Furthermore, we chose the VV polarization mode with stronger penetration, rather than the VH polarization mode. Subsequently, ENVI/IDL technology was used to perform batch preprocessing on the downloaded images, including master image selection, master-slave image georeference, and digital elevation model (DEM) phase simulation, thereby obtaining interferometric pairs with spatial baseline, temporal baseline, and Doppler centroid frequency baseline coherence coefficients.

2. Differential interferometry calculation mainly includes interferometric phase calculation, permanent scatterer candidate (PSC) point selection, flat earth and topographic phase removal, and differential interferogram calculation. Using the interferograms obtained in the first step, target points with good coherence are selected, and by performing spatial-temporal characteristic analysis on each phase component in the differential phase, two-dimensional regression and objective function optimization methods are used to solve for the linear components. The linear deformation components of the permanent scatterer (PS) points are removed from the original differential interferometric phase to obtain the residual phase components.

3. The first-step inversion uses the first inversion model to calculate residual height and displacement velocity, which is used to process complex interferograms. This method determines a certain number of “Persistent Scatterers (PS)” and then processes the pixels around each PS point. PS points need to satisfy two conditions: first, they must be stable; and second, they can be detected by the Synthetic Aperture Radar (SAR) antenna through proper orientation [37].

4. The second-step inversion is the final step of the inversion process, using the linear model results from the previous step to estimate the atmospheric phase components, thereby removing the atmospheric delay phase ${\phi }_{atm}$ and noise phase ${\phi }_{noise}$. Finally, the model of the final displacement velocity is fitted, and the displacement for each date is extrapolated.

5. The results of PS-InSAR are geocoded and output in both vector (shapefile) and raster formats. In this study, the Product Temporal Coherence Threshold was set to 0.85, which not only retained most of the effective feature points, preventing insufficient points for analysis, but also eliminated a large number of low-coherence feature points to ensure data accuracy.

3.3. Amplitude Factor Design

In mathematics and statistics, volatility is a measure of the magnitude of variation in a random variable, describing the degree of dispersion or instability of the random variable [38]. The Generalized Autoregressive Conditional Heteroskedasticity (GARCH) model is considered one of the most mature models for studying volatility due to its ability to effectively characterize features such as volatility clustering, asymmetry, and long memory [39]. The basic equations of the GARCH model are as follows:

$$\left\{\begin{array}{c}{\mu }_t={\sigma }_t{\epsilon }_t\\ {\sigma }_t^2={\alpha }_0+\sum _{i=1}^q{\alpha }_i{u}_{t-i}^2+\sum _{j=1}^p{\gamma }_j{\sigma }_{t-j}^2\end{array}\right.$$

(5)

where ${\sigma }_t^2$ is the current period volatility, ${\epsilon }_t$ is an independent and identically distributed white noise with zero mean and unit variance, and ${\mu }_t$ is the current period residual.

In the GARCH evaluation model, ${\sigma }_t^2$ represents the magnitude of current volatility, and the mean of the conditional heteroscedasticity represents the magnitude of historical volatility. If the mean of the conditional heteroscedasticity is high, it indicates that the data has high volatility. Moreover, the standard deviation of the conditional heteroscedasticity can reflect the stability of the data, the kurtosis can reflect the degree of peakedness, and the median can reflect the central tendency of the data [38].

To enhance the precision of our research findings, we divided the study area in ArcGIS into 0.005° × 0.005° grid cells, creating 14,637 evaluation units. We then input time-series displacement data for each grid cell into the GARCH model to compute key statistics - mean, standard deviation, kurtosis, and median of conditional heteroscedasticity. These values were then used in ARCGIS for an initial classification based on the quartile method, ensuring an equal number of evaluation units in each segment. We established a scoring system to assess surface volatility in the area (refer to Table 2). Each index was given a maximum score of 7 points, with scores and weights assigned according to the respective intervals of each index.

4. Results

4.1. Annual Results of Rebound Deformation in Wuxi Area

After processing 100 SAR images using the method introduced in Section 3.2, we obtained the surface deformation velocity field of the Wuxi area from November 2015 to June 2023. This study takes one year as a research time interval, allowing an error of plus or minus two months. Figure 4 shows the surface deformation velocity map of the Wuxi area obtained using the PS-InSAR technique, where positive values indicate uplift and negative values indicate subsidence. The velocity map is overlaid on a satellite image of the Wuxi area obtained from the Mapbox platform.

As shown in Figure 4, the Wuxi area does not exhibit a consistent subsidence or uplift trend in the year-by-year analysis, but there are significant fluctuations, mainly concentrated in regions BH, XS, and LX. Region BH shows a trend of large changes initially, followed by gradual stabilization, and then large changes again, with poor stability of surface deformation in recent years. Region XS is located between regions XW and JY, with its eastern part mainly consisting of villages and towns, showing good stability; while the western part is an urbanized area, with more evident surface deformation compared to the eastern part. Region LX is located in the central part of the main urban area of Wuxi City, and as time progresses, the surface fluctuations become more evident.

Region BH is located in the Binhu District of Wuxi City. Since the Jiangsu Provincial People’s Congress proposed the “Decision on Prohibiting Groundwater Extraction in Suzhou, Wuxi, and Changzhou Areas within a Time Limit” [40] in 2000, Wuxi City has comprehensively banned groundwater extraction since 2005. After a long period of development, the groundwater in the area has been fully replenished. Therefore, compared with previous studies [30], the surface subsidence occurring within the study area of this paper has been significantly alleviated. However, instability still exists in the area. The Lingshan Grand Buddha Scenic Area on its west side experienced a certain degree of subsidence in 2016, followed by significant uplift in 2017 and 2018; while after 2020, it experienced subsidence of varying degrees again. The Hongsha Bay area on its southeast side also experienced significant fluctuations in 2016-2017, followed by a smoothing trend, but in the first half of 2023, it experienced relatively evident surface subsidence again.

Region XW is located in the Xinwu District on the southeast side of Wuxi City. From the time-series subsidence rate map, it can be seen that the annual average (uplift/subsidence) rate in the area is between -5.19 and 1.93 mm/year, with little change, but there are evident fluctuations. According to the “Adjustment Plan for the Overall Land Use Planning of Xinwu District, Wuxi City (2006-2020)” [41] issued by the People’s Government of Xinwu District, Wuxi City has strictly restricted the urban and industrial construction land in Xinwu District. From 2014 to 2020, the urban construction land increased by only 2.23%, so there has been no significant subsidence in the area.

Region JY, situated within Jiangyin City under Wuxi City’s jurisdiction, witnessed significant developments following the 2012 endorsement of the “Jiangyin City Master Plan (2011-2030)” by the local government [42]. This plan aimed to foster rapid urbanization and enhance central urban area construction, while stringently safeguarding ecological assets like rivers, essential farmland, and riverside ecological wetlands, as well as conserving resources such as land and groundwater. These strategies have propelled Jiangyin City’s swift progress in recent years, simultaneously moderating surface deformation fluctuations. From 2015 to 2023, the southern townships and villages saw no notable surface subsidence. However, the northern riverside area, constituting Jiangyin’s primary urban zone, experienced the most pronounced deformation fluctuations in our study. Its proximity to the Yangtze River significantly influences surface changes through pore confined water levels, a topic we delve into as a focal point in Section 5.

4.2. Long Time Span Results of Surface Deformation in Wuxi Area

The surface deformation velocities of the Wuxi area from 2015-2018, 2018-2020, and 2020-2023 obtained through the PS-InSAR method are shown in Figure 5. Positive values in the figure indicate ground uplift, while negative values indicate land subsidence.

The period from December 2015 to June 2018 is a set of SAR images containing 16 scenes. Figure 5a shows the surface deformation velocity map along the LOS direction for this period. From 2015 to 2018, a total of 5,652,151 PS points were obtained, with deformation velocities ranging from -23.97mm/year to 29.09mm/year and an average of 6.52mm/year. 90.3% of the surface deformation velocities in the study area are between -6.02mm/year and 11.01mm/year. During this period, the surface deformation in the Wuxi area generally shows an uplifting trend, with some areas remaining unchanged. The Huishan District in western Wuxi exhibits a significant large-scale uplift, with an average annual uplift rate between -2.47mm/year and 18.78mm/year. Additionally, several key deformation zones with uplift rates exceeding 12mm/year appear along the Yangtze River in the northern Jiangyin City. Compared with the results of the previous segmented study, the subsidence zone in the south in 2016 partially overlaps with the uplift zone in the south in 2017, and the uplift zone in the north in 2016 partially overlaps with the subsidence zone in the north in 2017. Therefore, in these two overlapping areas, the surface deformation has undergone a certain neutralization effect. Furthermore, in the first half of 2018, the overall area experienced a significant uplift, causing the overall deformation effect to match the aforementioned image results.

The period from February 2018 to December 2020 is a set of SAR images containing 18 scenes. Figure 5b shows the surface deformation velocity map along the LOS direction for this period. From 2018 to 2020, a total of 2,544,221 PS points were obtained, with deformation velocities ranging from -40.03mm/year to 22.90mm/year and an average of 0.29mm/year. Compared to the period from 2015 to 2018, the overall subsidence trend is relatively gentle, with no large-scale obvious uplift or subsidence. By comparing the three annual average velocity maps from 2018 to 2020 in Section 4.1, it can be seen that although the Binhu District in southwestern Wuxi and the Liangxi District near the central area experienced large-scale subsidence in 2020, they both experienced slow uplift in 2018 and 2019; the subsidence and uplift have a neutralization effect, resulting in a stable change trend in the long time series. Moreover, the Xinwu District in southeastern Wuxi, the Xishan District in the central area, and the junction area with Jiangyin City, although experiencing significant uplift in certain specific years, show a stable trend in the remaining years; this implies that the uplift effect is weakened, resulting in a stable trend in the overall deformation.

The period from June 2020 to June 2023 is a set of SAR images containing 19 scenes. Figure 5c shows the surface deformation velocity map along the LOS direction for this period. From 2020 to 2023, a total of 3,536,105 PS points were obtained, with 88.53% of the PS points having deformation velocities between -25.98mm/year and 30.01mm/year. Unlike the significant uplift trend from April 2022 to June 2023, the surface deformation during this study period is relatively small. The main reason is that during the period from June 2020 to March 2022, some areas experienced continuous small-scale subsidence, thus weakening the subsequent uplift trend.

This section conducts mathematical statistics on the deformation information obtained by the PS-InSAR technique for two time periods. To further understand the time series changes of the deformation velocity, this section separately counts the surface subsidence and uplift information, and the statistical results are shown in Figure 6.

According to the subsidence statistics chart, from 2015 to 2023, the subsidence rate in Wuxi shows a decreasing trend. From 2015 to 2018, PS points with subsidence rates between 0-4mm/year account for 62.16% of all subsidence points, and this proportion reaches 92.97% in 2018-2020 and exceeds 95% in the last study period. However, for PS points with subsidence rates above 4mm/year, the number in the first time period far exceeds the latter two time periods, and PS points with subsidence rates exceeding 8mm/year rarely appear after 2018. From the uplift statistics data, the uplift in the Wuxi area shows a fluctuating trend of increasing first and then decreasing. In 2018-2020, PS points with uplift rates between 0-4mm/year are fewer than the two time periods before and after, while PS points with uplift rates above 4mm/year are more than the two time periods before and after. This also indirectly reflects the periodic change trend of surface deformation in the Wuxi area.

In this study, we obtained the surface deformation of the Wuxi area from 2015 to 2023. To study the characteristics of early surface deformation and compare them with our results, we summarized the surface subsidence studies conducted in the Wuxi area in the past and presented them in Table 3.

From 2004 to 2011, Wuxi City formed two major subsidence areas: Jiangyin and Xishan. Lu et al. analyzed the surface subsidence in Wuxi City from 2004 to 2017 and found that from 2004 to 2011, large subsidence funnels appeared in the southern part of Jiangyin City and the Huishan District bordering Jiangyin City, with a maximum subsidence rate exceeding 20mm/year [30]. After 2011, InSAR results indicated that the subsidence areas in Jiangyin City and Huishan District tended to flatten, and the overall subsidence expanded to the southeastern part of Wuxi [27]. After 2015, the subsidence in the junction area of Xishan District and Jiangyin City in Wuxi City showed a slowing trend, and relatively obvious subsidence areas appeared in Xinwu District and Binhu District. Although previous studies [43] used different data and methods to monitor surface deformation, they were able to identify similar subsidence/uplift areas and surface deformation rates. Compared with the research results of previous studies, the spatiotemporal distribution characteristics of the subsidence areas detected in this paper are basically consistent with the research results of previous studies [28].

4.3. Accuracy Verification of Rebound Results

Using 100 Sentinel-1 SAR images covering the Wuxi area, surface deformation information from November 2015 to June 2023 in the Wuxi area was obtained through the PS-InSAR technique. After three data processing iterations, more than 117 thousand PS points were obtained. We calculated the amplitude dispersion index (i.e., the ratio of the mean to the standard deviation) for all PS points as a precision quality factor, thus effectively assessing the internal precision quality of the obtained surface deformation information.

In Figure 7a, 94.38% of the total points have a precision quality better than 7.0mm/year. In Figure 7b, 96.89% of the total points have a precision quality better than 7.0mm/year. In Figure 7c, 94.41% of the total points have a precision quality better than 7.0mm/year.

The above results indicate that using Sentinel-1 ascending data with IW polarization mode to obtain surface information in Wuxi, combined with PS-InSAR technology for processing and analysis of the area, has very high precision and reliability. In addition, we employ cross-validation methods to verify accuracy, with specific details provided in the supplementary material.

5. Discussion

5.1. Time Series Analysis of Characteristic Points

To better analyze the spatial distribution characteristics of surface rebound changes in Wuxi City, this section randomly selected six different characteristic points distributed across Wuxi City for detailed analysis. The method for selecting characteristic points was as follows: first, the location of the characteristic points was determined, then a circle with a radius of 500m centered on these points was drawn, designating the characteristic area. The temporal cumulative deformation of the characteristic points is shown in Figure 8a, with their specific distribution locations illustrated in Figure 8b.

PS1 is located near Hongsha Bay in the Binhu District, close to the Taihu Lake area. This area exhibited continuous small-scale fluctuation changes, with the maximum cumulative uplift being 8.55mm and the maximum cumulative subsidence being -16.62mm. Due to changes in the water level of Taihu Lake, which lead to variations in the confined water level in this area, water level changes are the main factor causing these regional variations.

PS2 is located in the Shuofang Airport area of the Xinwu District. This area showed a trend of slight subsidence followed by slow uplift. The area has many new technology and industrial parks, the construction of which increases the compressive stress on the surface, causing subsidence; this is followed by slow rebound due to substantial groundwater recharge.

PS3 is located near the Liangxi District People’s Government. Prior to 2021, this point experienced wide-ranging fluctuation changes, but following the issuance of the “Fourteenth Five-Year Plan for National Economic and Social Development of Liangxi District, Wuxi City and the Vision for 2035”[45], Liangxi District strengthened its ecological civilization system construction, strictly controlled urban development boundaries, and scientifically managed the relationship between urban development and land resource utilization; therefore, in subsequent periods, this area showed a stable trend.

PS4 is located near the Huishan District railway station. There were no significant fluctuation changes here, showing a trend of slow uplift.

PS5 is located near the Xibei Town People’s Government in Xishan District. This area has many factories, such as vehicle parts, railway parts, and plastic factories, with a significant proportion of industrial water use compared to other areas; and since 2018, the region has focused on creating major industrial projects and extensive construction, leading to poor surface stability in this area.

PS6 is located near Xiangshan Village, at the junction of the southern part of Jiangyin City and Xishan District. The figure shows that this area has weak fluctuation, with an overall trend toward stability. After April 2022, this region experienced a significant uplift, related to increased precipitation and groundwater changes, which will be discussed in detail in the following text.

Overall, the fluctuation changes in the Wuxi area are uneven. This situation is mainly related to groundwater, which will be discussed in focus in the next section.

5.2. Changes of Surface Volatility and Groundwater Level

In the groundwater layer, the extraction of groundwater leads to a decrease in the water level, which inevitably reduces the pore water pressure in the aquifer itself and the relatively impermeable layers above and below it. According to the effective stress principle, the total stress in the soil caused by the overburden load is jointly borne by the water in the pores and the soil particle skeleton. The portion borne by water is called pore water pressure $P_W$, which cannot cause soil layer compression and is thus also called neutral pressure; the portion borne by the soil particle skeleton can directly cause upper layer compression and is thus called effective stress $P_S$; the sum of the two equals the total stress [46]. Assuming that the internal stress of the soil layer remains unchanged during the water extraction process, the decrease in pore water pressure inevitably leads to an equal increase in the effective stress in the soil, which in turn causes a reduction in pore volume and compression of the soil layer. Under lower stress conditions, sand layers exhibit low compressibility and mainly demonstrate elastic, reversible deformation, while clay layers exhibit much higher compressibility and primarily undergo inelastic, permanent deformation. Therefore, under conditions of lower effective stress increase, the consolidation of cohesive soil layers plays a dominant role in land subsidence, while the expansion and rebound of sand layers are crucial during the water level recovery process. The schematic diagram of effective stress is shown in Figure 9, where the soil particles above the groundwater level are not displayed.

This section selects the groundwater level data from six local confined aquifer groundwater monitoring wells from January 2018 to December 2021, using monthly intervals to form time series of groundwater level changes, with the well locations distributed as shown in Figure 10G. The InSAR elevation time series are derived from all PS pixels within a 200 m radius around each observation well, and the LOS direction displacement data (red line) and water level changes (blue line) are normalized and comprehensively compared, ultimately plotted in Figure 10A–F.

The stratigraphic structure in the Wuxi area is relatively complex, with thick sand and gravel-bearing sand layers, while the clay layers are thin and unevenly distributed. According to the above figures, it can be observed that all areas near the wells exhibit a slow rebound trend, which is closely related to the expansion and rebound of the sand layers. Furthermore, the changes in confined water levels have a high correlation with the InSAR-derived surface elevations. We employ the GRA method to calculate the relationship between elevation and water level, obtaining GRA values of 0.52, 0.71, 0.61, 0.75, 0.57, and 0.65 for wells (A), (B), (C), (D), (E), and (F), respectively. Well (A) is located in the Shuofang Economic and Trade Service Center, Binhu District, Wuxi City, while well (E) is situated in Huangtang Village, Jiangyin City, one of the county-level cities in Wuxi. From January 2018 to December 2022, the cumulative displacement components in the LOS direction for wells (A) and (E) are 14.29 mm and 12.59 mm, respectively, while the groundwater level changes are -0.82 mm and 10.85 mm, respectively. The groundwater level changes in well (E) are consistent with the LOS direction displacement changes, while well (A) shows significant differences. Moreover, the displacement change trends of wells (A) and (E) have low correlations with water level changes. The reasons for the differences are as follows: (1) the local soil and geological conditions of the two wells are different, reflecting varying impacts of confined water on surface deformation; (2) the influence of groundwater level changes on surface deformation has a certain lag, and for some areas, the impact is not immediately apparent but has a certain delay; (3) the groundwater observation accuracy of well (E) is insufficient, with some abnormal values; (4) surface subsidence is affected by multiple factors, with groundwater being only one of them, and the combined effects of other aspects lead to suboptimal fitting results. Well (B) is located in Yangjian Town, Xishan District, while well (D) is situated in the Wenlin Community Service Center, Jiangyin City, both located in the border area between Xishan District and Jiangyin City. From the Figure 10, it can be seen that the surface deformation and water level change trends of the two wells are generally similar; when the groundwater level rises, the surface undergoes uplift and rebound, and when the groundwater level decreases, the ground subsides. Their gray relational degree indices both exceed 0.7, indicating that changes in confined aquifer water levels are the main cause of surface deformation in this subsidence area. Well (C) is located in Qianzhou Street, Huishan District, while well (F) is situated in Lingang Street, Jiangyin City. From the Figure 10, it can be observed that these two areas experienced significant uplift in 2018, with well (C) uplifting by 24.60 mm and well (F) by 23.45 mm. The groundwater level exhibited a corresponding increase, causing changes in surface floating effect. The reasons for the water level rise may be related to artificial groundwater recharge and natural precipitation, which will be discussed in detail in Section 5.3.

Overall, surface subsidence deformation is closely related to confined aquifer water. In recent years, with the government’s attention and legal protection, artificial groundwater recharge has effectively suppressed surface subsidence deformation in the Wuxi area. However, reasonable limits should also be set for artificial groundwater recharge to avoid excessive recharge leading to excessive surface rebound. In the future, ground-based SAR imaging radar devices can be deployed around the water wells, combined with geodetic surveying, to conduct real-time monitoring of subsidence around the wells, making the surface more stable.

5.3. Impact of Precipitation on the Surface Flotation Effect

According to the “Notice on Carrying Out Systematic and Comprehensive Demonstration Work for the Construction of Sponge Cities” issued by the Ministry of Finance, Ministry of Housing and Urban-Rural Development, and Ministry of Water Resources, Wuxi City was designated as one of the first batch of sponge city construction demonstration cities in China in 2021 [47]. As one of the key elements of the sponge body, soil is directly affected by precipitation. Previous studies have shown that rainwater infiltration increases soil moisture content and reduces soil porosity and bulk density. This implies that when precipitation increases, the soil behaves like a water-saturated sponge, undergoing uplift and elastic rebound while indirectly recharging groundwater.

Wuxi is located in the humid subtropical monsoon climate zone of North Asia, bordering the Yangtze River to the north and Taihu Lake to the south, characterized by abundant precipitation and concurrent rain and heat. This study collected monthly average precipitation data from January 2018 to December 2021 provided by the official website of the China Meteorological Administration (https://weather.cma.cn/) and divided the Wuxi area into five latitudinal zones from north to south (N1:31°51′36″-31°59′24″; N2:31°43′12″-31°51′36″; N3:31°35′24″-31°43′12″; N4:31°27′36″-31°35′24″; N5:31°19′48″-31°27′36″), with the longitude ranging from E119°41′13″ to E120°39′58″.

As shown in Figure 11, the precipitation amounts in the five latitudinal zones are very similar, with more precipitation in the south and less in the north, consistent with the meteorological conditions of a humid south and arid north in China. Therefore, this study averages the monthly precipitation of these five latitudinal zones and compares them with the four characteristic regions in Figure 12b to analyze the impact of precipitation on surface floating effect.

From Figure 12a, it can be observed that from January 2018 to early 2020, precipitation is concentrated from May to September each year, with less precipitation in other months; after 2020, precipitation is concentrated from June to August each year, with the precipitation during this period far exceeding the flood season precipitation of the previous period, while in other dry seasons, precipitation is reduced compared to the previous period. This indicates that the precipitation distribution characteristics from 2020 to 2022 are concentrated and high in volume. Overall, due to the increased precipitation during the flood season from 2020 to 2022, the annual precipitation after 2020 is greater than that before 2020. P1 is located near Shuofang Airport in Xinwu District, and P4 is located in the Hongsha Bay area of Binhu District. The common feature of both is that during seasons with increased precipitation, the surface undergoes uplift, while during seasons with less precipitation, a subsidence trend occurs. This indicates that precipitation and surface elevation exhibit a partial linear correlation, reflecting the surface flotation effectcaused by precipitation infiltration. However, during a few periods, when precipitation increases, the surface does not exhibit a corresponding response, indicating that surface subsidence is related to multiple factors, with precipitation being a secondary factor at this time. P2 is located in the Qianzhou Street area of Huishan District, near well (C) in Section 5.2; P3 is located in Lingang Street, near well (F) in Section 5.2. We find that although both groundwater and surface levels in these two regions exhibited a certain degree of increase in 2018, the precipitation in 2018 was relatively low compared to other years. This directly reflects that the uplift occurring in P2 and P3 in 2018 has a smaller relationship with precipitation. We speculate that artificial groundwater recharge measures directly increased the groundwater level, enhancing the neutral stress and resulting in significant surface uplift; additionally, the soil and geological structure of the Wuxi area may also have a certain impact on this deformation trend, which will be specifically explained in Section 5.4.

In summary, precipitation affects the surface flotation effectthrough rainwater infiltration, increasing soil pore volume and leading to soil layer expansion and surface rebound. However, for certain periods, the impact of precipitation on the surface is not significant, as surface elevation is influenced by a combination of multiple factors. Furthermore, this section also analyzes the correlation coefficients between precipitation and surface changes, which are 0.67 (P1), 0.72 (P2), 0.71 (P3), and 0.70 (P4), confirming that the influence of precipitation is relatively limited.

5.4. The Impact of Soil Stratigraphy and Quaternary Sedimentary Regions

Wuxi, situated in the heart of the Yangtze River Delta, is a vital city that links the Yangtze River and Taihu Lake. The region is characterized by plains interspersed with low mountains and hills [26]. To better study the impact of soil layer distribution on surface fluctuation, we obtained data for three typical boreholes in the Wuxi area from the China Geological Science Institute, focusing on soil layers within a depth of 100 meters. ${\mathrm{B}\mathrm{H}}_a$ is located near Wanda Plaza in Jiangyin City, close to P3, reaching a final depth of 99.2m; ${\mathrm{B}\mathrm{H}}_b$ is near the intersection of Shudong River and Outer Ring South Road in Jiangyin City, close to P2, with a final depth of 94.8m; ${\mathrm{B}\mathrm{H}}_c$ is in the eastern part of Shuofang Airport in the Xinwu District of Wuxi City, with a final depth of 92m. The specific distribution of the boreholes and the internal distribution of rock and soil are shown in Figure 13.

As shown in Figure 13, the stratigraphic structure of Wuxi is quite complex, with sub-clay, silt, clay in the upper part, sub-clay and sub-sandy soil in the lower part, which together form the soft soil layer of the Wuxi area. During the Quaternary sedimentation process, affected by sea-land alternation, changes in the Yangtze River water level, and river course changes, the lithology of the region is complex, characterized by multi-layered cohesive soil layers and alternating loose sand layers, which is the direct reason for the interlayering of sub-clay and silt [28]. Additionally, the stratigraphic changes in Wuxi are highly uneven in the horizontal direction. Wu et al. through geological drilling and controlled source audio-frequency magnetotellurics, found that the resistivity of the Quaternary loose layer in this region varies greatly in the horizontal direction, reflecting the differences in the physical properties of the strata, and having a direct impact on the surface deformation [48].

The groundwater level in Wuxi is usually between -1 and -30 meters, with the soil mainly consisting of silt and sub-clay in the soft soil layer. The sub-clay has poor permeability and a dense structure, making it difficult for water to flow, while the silt has a looser structure and good permeability. Therefore, groundwater mainly affects surface deformation through the silt layer. In the early stages, due to excessive groundwater extraction, the water level dropped significantly, leading to severe water loss in some silt layers. At the same time, rapid industrial development increased the surface load pressure, compressing the silt layer, reducing the pores, and making the structure denser, directly causing large-scale subsidence. In recent years, due to Wuxi City’s implementation of a groundwater extraction ban, the silt layer’s moisture has been greatly replenished, leading to elastic changes in the surface. However, this kind of rebound change has certain limits. Whether it’s cohesive soil or sandy soil, due to long-term compression by load, irreversible plastic deformation occurs, and due to pore closure, it’s difficult for water to re-enter the soil layers. In other words, even with substantial groundwater recharge, the surface is unlikely to return to its previous height. Combining analysis with Section 5.3, P2 and P3, the significant uplift in these areas in 2018 was due to groundwater recharge and water level rise. However, after 2019, the uplift significantly slowed, and even subsidence occurred, because the soil layer’s elastic deformation had recovered, and most of the remainder was plastic deformation. Furthermore, the pores cannot continue to expand, and water cannot enter, thus displaying the trend shown in Figure 12 due to the combined effect of these factors.

Overall, the soil layers in the Wuxi area are primarily composed of sub-sandy soil and silty clay. As these soil layers are relatively soft, they have a strong correlation with surface fluctuation changes. Moreover, overly low groundwater levels in the soil layers can lead to plastic deformation, causing lasting and irreversible harm to the area; thus, we should strengthen real-time monitoring of soil layer deformation and groundwater levels.

5.5. Volatility Evaluation—Amplitude Factor

Continuous ground subsidence or uplift increases its instability. When the deformation amplitude reaches a certain level, it will cause incalculable losses to surrounding buildings and road facilities; therefore, it is essential to monitor and assess surface fluctuations. Traditional regional assessment methods require selecting factors such as fracture characteristics, rock and soil properties, and human activities as evaluation factors based on the overview of the study area. Due to the large research area, high survey difficulty, and multiple data sources, data missing problems are prone to occur in actual operations. Moreover, common surface stability assessment methods only propose evaluation areas at a relatively large level, and their accuracy is insufficient to guide specific engineering construction and earthquake disaster reduction [49]. Therefore, this study designed a volatility evaluation index based on the GARCH model—the amplitude factor. Among them, conditional heteroscedasticity is an important indicator of the model. Its mean, standard deviation, kurtosis, and median respectively reflect the magnitude, stability, kurtosis degree, and concentration trend of the data volatility. Using the method in Section 3.3, the surface stability evaluation results of the Wuxi area are shown in Figure 14a.

It can be seen from Figure14a that the northern riverside area of Jiangyin City is the region with poor stability in this study, with more than 64.29% of the amplitude factors greater than 5.4. In addition, the southeastern junction of Xishan District and Xinwu District, as well as the Taihu Lake coastal area of Binhu District, all exhibit unstable change characteristics. These deformation characteristics are consistent with the velocity field change trend analyzed in Section 4.1, verifying the accuracy of the model. According to the scoring system in Table 2 of Section 3.3, we obtained the minimum value of the amplitude factor in the Wuxi area as 1 and the maximum value as 7.

To perform a qualitative assessment of surface stability, this study employs the equal interval classification technique to categorize the values of the amplitude factor, as illustrated in Figure 14b. It’s evident that most areas in Wuxi are categorized as sub-stable or relatively stable, suggesting that the prohibition on groundwater extraction since 2005 has yielded notable outcomes. The northern riverside zone of Jiangyin, predominantly consisting of alternating layers of silt and clay, is highly susceptible to fluctuations in the Yangtze River’s water level, leading to compromised stability. Overall, the Wuxi surface stability evaluation model demonstrates significant precision, and its findings offer valuable guidance for urban planning and construction. This contributes to ensuring urban safety and preserving the successes of major projects in this new era.

6. Conclusion

This study utilized 100 Sentinel-1A images from November 2015 to June 2023 to investigate the spatiotemporal characteristics of surface deformation in the Wuxi area using the PS-InSAR technique. To better understand the causes of volatility changes, a comprehensive analysis of groundwater, precipitation, soil geology, and other aspects was conducted. Additionally, an amplitude factor was developed to assess the surface stability of the Wuxi area.

The main findings are as follows:

1) The internal accuracy of PS-InSAR was verified using the precision quality factor evaluation method, revealing that most errors are within 7mm, with almost no errors exceeding 10mm, indicating high accuracy.

2) Surface deformation in the Wuxi area are highly uneven, with significant spatial and temporal variations. The northern Jiangyin City area experienced strong surface deformation from 2015 to 2018, which slowed down slightly but showed significant uplift after 2021. The southern Binhu District and Xinwu District experienced substantial subsidence in 2016, followed by minor fluctuations, with the Binhu District experiencing some subsidence again from 2022 to 2023. Over time, the subsidence amplitude of all PS points has gradually decreased, with more than 90% of the characteristic points having a subsidence rate of 0-4mm/year after 2018, while the uplift amplitude fluctuates.

3) Multiple natural factors influence surface volatility in the Wuxi area. This study focused on analyzing the combined effects of groundwater, precipitation, and soil geology. Data from six groundwater monitoring wells showed a strong correlation between groundwater level changes and surface deformation. For wells C and F, the analysis revealed that surface changes are not only related to precipitation, groundwater abstraction bans, and recharge but also to the area’s soft soil geology. Sandy soil and clay are prone to water loss, shrinkage, and deformation under load pressure, further increasing surface volatility.

4) The GARCH model was used to analyze the time-series subsidence displacement in Wuxi City in recent years. By comprehensively analyzing the conditional heteroscedasticity’s mean, standard deviation, kurtosis, and median, an innovative “amplitude factor” surface volatility evaluation index was proposed. The surface stability of the Wuxi area was qualitatively displayed using a hierarchical grading method based on the magnitude of the amplitude factor. The study found that sub-stable areas are primarily concentrated in northern Jiangyin City and the southeastern junction of Xinwu District and Xishan District, consistent with the research results presented in the previous text.

These findings provide insights into the surface deformation characteristics of Wuxi following groundwater abstraction bans and recharge, offering a basis for future construction and safety control by the local government. Continuous monitoring and analysis of surface volatility changes in the Wuxi area are essential for the city’s future safety and development.