1. Introduction
The mainstream of the Yangtze River flows through the entire territory of Jiangsu, China. The geological conditions in the Nan**g area along the river are characterized by loose geology, poor bearing capacity, and shallow groundwater levels. Due to long-term river water accumulation, the area is prone to collapse and sedimentation, resulting in land subsidence. Additionally, with the rapid development of the economy in the Yangtze River Delta region, the demand for engineering construction is increasingly expanding. Various engineering activities, such as groundwater extraction, subway engineering, and large ground building construction, have accelerated the occurrence of regional ground subsidence [
1,
2]. Under the joint action of the above natural and human factors, the geological environment system along the Yangtze River Delta region has been damaged, leading to urban land subsidence disasters, which, in turn, cause serious hazards such as soil collapse, road facility damage, and building wall cracking, posing a certain threat to the security of residents. Therefore, it is crucial to strengthen the monitoring and analysis of surface deformation within a certain range along the Yangtze River using advanced techniques.
Traditional approaches for ground subsidence monitoring mainly include precision leveling, trigonometric leveling, and GNSS measurements [
3]. Precision leveling requires less cost, a simple construction process, and can ensure measurement accuracy. However, it is limited in the number and scope of leveling points and has a long monitoring cycle, resulting in low temporal and spatial resolution [
4]. Trigonometric leveling, as an indirect height measurement method, features minimal terrain limitations and fast measurement speed, but its accuracy is easily affected by atmospheric refraction and is difficult to calibrate [
5,
6]. While GNSS boasts all-weather operation, is fully automatic, and has high accuracy, the high monitoring cost makes it challenging to conduct large-scale ground subsidence monitoring [
7,
8]. Traditional manual measurement methods face various limitations such as measurement range, manpower consumption, timeliness, and cost. Synthetic aperture radar interferometry (InSAR), as a form of active microwave remote sensing technology enabling continuous, all-weather monitoring, can provide low-cost solutions for acquiring wide-ranging and high-precision deformation information [
9,
10,
11]. Therefore, it serves as an important tool for large-scale ground deformation monitoring.
Conventional InSAR technology obtains DEM by calculating the phase difference between two SAR images [
12], while D-InSAR technology utilizes DEM data for interferometric processing and filters out other phase information from the interferometric image, retaining only the deformation phase resulting from displacement [
13,
14]. This enables the extraction of precise ground deformation information at millimeter-level accuracy. However, the monitoring accuracy of D-InSAR is significantly hindered by various factors, including temporal and spatial decorrelation, atmospheric effects, and satellite orbit errors. Further, this method cannot monitor the development process of deformation [
15,
16]. To overcome the limitations and shortcomings of conventional InSAR techniques, researchers have proposed the Persistent Scatterer InSAR (PS-InSAR) [
17,
18] and Small Baseline Subsets InSAR (SBAS-InSAR) [
19,
20] approaches. These multi-temporal InSAR techniques effectively mitigate the impact of atmospheric effects and spatio-temporal decorrelation on the accuracy of surface deformation monitoring. Consequently, they have found widespread applications in various fields such as land subsidence monitoring, geological hazard assessment, and building deformation analysis [
21,
22,
23,
24].
This research acquired 42 Sentinel-1A images along the Yangtze River in Nan**g and applied PS-InSAR and SBAS-InSAR technologies to monitor and analyze the land subsidence in the area. The accuracy and reliability of the two technologies in land subsidence monitoring were evaluated through cross-validation of the continuous deformation field obtained from both techniques and comparison with external leveling data. Additionally, the causes of land subsidence in this area were analyzed from three perspectives: underground transportation engineering, building loads, and geological structure.
2. Research Area and Experimental Data
2.1. Research Area
Nan**g, situated in the lower reaches of the Yangtze River, serves as the capital of Jiangsu Province and holds strategic importance as a major intersection of the Belt and Road Initiative and the Yangtze River Economic Belt. The city is traversed by the Yangtze River, flowing from southwest to east, with its urban area encompassing significant water bodies such as the Qinhuai River, Qinhuai New River, and Xuanwu Lake. The region along both banks of the Yangtze River in Nan**g is profoundly influenced by the impact and sedimentation of river water, resulting in the formation of a broad floodplain. The geological structure in this area is characterized by soft soil, rendering it susceptible to collapse and sedimentation, and the bearing capacity is poor.
In terms of tectonic units, Nan**g is in the northern margin of the Yangtze ancient land, and the basement is mainly light metamorphic schist and metamorphic volcanic rocks. For the area along the Yangtze River in Nan**g, the overall geological conditions are not very ideal, the main characteristics are loose geology, poor bearing capacity, and shallowly buried groundwater levels. In recent years, Nan**g has experienced rapid urbanization, marked by the extensive construction of underground railways, intensive implementation of ground construction projects, and excessive exploitation of groundwater. Consequently, some areas along the Yangtze River in Nan**g have witnessed land subsidence, posing significant challenges not only to the economic development of the Yangtze River region, but also to the safety and well-being of urban residents. In light of these circumstances, this research focused on the land subsidence along the Yangtze River, especially the 5 km buffer zone on both sides of the Yangtze River in Nan**g as a typical area for land subsidence monitoring and analysis, as shown in
Figure 1. A–F were the positions of six leveling points. The findings aim to provide valuable technical insights for urban management and protection.
2.2. Experimental Data
The Sentinel-1A satellite, operating in a solar synchronous orbit, was launched in April 2014. It is the inaugural environmental monitoring satellite deployed under the European Space Agency’s Copernicus program and is equipped with a C-band SAR. Notably, it boasts high revisit frequency, wide coverage, and multiple polarimetric and working modes, which can meet the needs of time-series InSAR processing and analysis. For this study, approximately one image was selected every two months, culminating in a total of 42 Sentinel-1A images in Level-1 SLC format, which covered the areas along the Yangtze River in Nan**g from 8 April 2015 to 14 November 2021. The images were acquired in IW mode, with VV polarization, an incidence angle of 35.9°, and a spatial resolution of 5 m × 20 m (azimuthal direction × range direction). The Frame of the used data were 99 and 100, and the Path was 69.
Furthermore, this experiment utilized precision orbit data of the Sentinel-1A satellite, which was provided by the European Space Agency, to refine the orbit information. The reference digital elevation model (DEM) data utilized in this study were sourced from the SRTM elevation data with a resolution of 30 m. To validate the accuracy of the land subsidence information extracted by two methods, on-site leveling data spanning from 2015 to 2020 were employed.
Figure 1 illustrates the spatial distribution of the six leveling points within the study area.
3. Technical Principles
3.1. PS-InSAR Technology
The Persistent Scatterer InSAR (PS-InSAR) technology, initially introduced by Ferretti in 2000, is designed to analyze objects within the study area that exhibit stable scattering characteristics over prolonged periods [
17,
18]. This involves identifying a set of ground target points with high coherence and stable reflection characteristics in image data, termed PS points, to extract phase changes. By doing this, this approach effectively avoids the problems of temporal and spatial decorrelation and reduces the influence of atmospheric effects on deformation information. As a result, PS-InSAR enables the acquisition of high-precision ground deformation information [
25].
The fundamental principle of PS-InSAR technology is to select one SAR image as the master image, followed by the registration of N additional SAR images with the master image to generate N interference pairs. Subsequently, the differential processing of these pairs yields N interferograms, from which PS points with high coherence and stable reflection characteristics are extracted. Furthermore, the interference factors in phase information extraction are eliminated using orbital parameters, external DEM, and an established model. Ultimately, high-precision ground deformation information can be obtained. The interference phase of the PS point can be expressed as:
Here, is the interference phase, is the deformation phase caused by the displacement of the target object during the satellite revisit cycle; refers to the flat phase caused by interference on a flat surface caused by radar oblique range imaging, which can be eliminated with the assistance of precise orbital parameters; is the terrain phase caused by surface undulation, which can be eliminated with the assistance of DEM; is the atmospheric delay phase caused by the influence of the ionosphere and troposphere on radar signals, resulting in signal propagation delay; and is the phase error caused by radar thermal noise and inherent interference noise, which can be reduced by filtering.
With the assistance of precise orbit parameters and an external DEM, the interference image is flattened, and the terrain phase is eliminated. To express the relationship between the terrain residual increment and the interference phase of PS points, a model is established. Iterative calculations are performed using a linear model to estimate the surface deformation rate and the residual increment of terrain errors. Subsequently, after removing the residual deformation phase and terrain phase, further iterations are conducted to estimate the atmospheric delay phase error and noise error. This iterative process leads to the ultimate retrieval of the deformation rate and the residual elevation.
In this experiment, the selection of the master image was based on the principle of maximizing correlation among interference pairs. Ultimately, the SAR image obtained on 2 September 2019 was chosen as the master image. It was combined with the remaining 41 SAR images to form 41 interference pairs. As depicted in
Figure 2, the spatial baseline exhibited a distribution primarily ranging −150~150 m, with the longest spatial baseline measuring 132.44 m and the shortest at −8.68 m. Subsequently, external DEM data were utilized to facilitate image registration and mitigate the impact of flat ground effects. The amplitude discrete index was employed to preliminarily select high-coherence target points. To obtain the initial deformation rate and residual elevation of the PS points, a linear inversion model was applied. Furthermore, high-pass filtering in space and low-pass filtering in time were employed to estimate the atmospheric delay phase. Geocoding of the PS points was accomplished using the reference DEM, with a coherence coefficient threshold of 0.7. After removing PS points with low coherence, a total of 976,871 effective PS points were obtained, enabling subsequent temporal and spatial analysis.
Due to the large coverage area of the acquired Sentinel-1A image, the processing was time-consuming. As a result, 42 images were cropped using geocoded vector data, and subsequent experimental processing was carried out based on these cropped images.
3.2. SBAS-InSAR Technology
The Small Baseline Subset InSAR (SBAS-InSAR) technique was initially proposed by Berardino et al. in 2002 [
19]. The fundamental concept of SBAS-InSAR technology involves the selection of M SAR images for flexible combination, ensuring that each image forms at least one interference pair with another image. By establishing thresholds for temporal and spatial baselines to exclude image data that do not meet the criteria, the adverse effects of spatiotemporal incoherence on monitoring results are mitigated. Subsequently, following the small baseline combination principle, differential interference processing is performed on all image pairs to derive N differential interferograms. The flat-earth phase and terrain phase in the interference phase are eliminated by removing the flat-earth effect and leveraging external DEM data. The interferometric phase of the interferogram generated by two SAR images can be expressed as the following formula,
Here, and are the interferometric phases obtained at the imaging time of two SAR images, respectively. is the residual terrain phase error, and is the atmospheric delay phase. is the noise phase. In general, SBAS-InSAR technology operates under the assumption that surface deformation follows a linear pattern, and hence, a linear model was constructed to characterize the surface deformation. The residual terrain phase present in the interferometric phase was estimated and removed using the singular value decomposition method. Additionally, the atmospheric phase and nonlinear deformation phase were eliminated through spatial domain low-pass filtering and time domain high-pass filtering. Subsequently, the least-squares method was employed to calculate the average phase change rate and surface deformation for each given time.
In this experiment, the spatial baseline threshold was set to 48%, and the time baseline threshold was 180 days. The SAR image acquired on 7 September 2018 was selected as the master image. The generated image pairs and baselines are illustrated in
Figure 3. Subsequently, the interference processing was assisted by external DEM data. The Goldstein filtering method and the Delaunay MCF unwrap** method were employed [
26]. After evaluating the noise area in the unwrap** output and considering the phase information of the research area, the coherence threshold was finally set at 0.2. Following this, orbit refinement and re-phasing were conducted. The residual terrain error and atmospheric delay phase were estimated using two inversions, and finally, geocoding was performed.