2.1 Study area
The city of Tianjin (located at 116°43′-118°4′E, 38°34′-40°15′N) lies in the northeastern part of the North China Plain. It is bordered by the Bohai Sea to the southeast and the Yanshan Mountains to the north. To the northwest, it neighbors China’s capital, Beijing, and Hebei Province. The city boasts an extensive transportation network, including waterways, railways, aviation, and highways (Fig. 1a, 1b).Tianjin experiences a warm temperate continental monsoon climate, characterized by distinct seasons. Summers are hot and rainy, while winters are cold and dry, with an average annual temperature of 11–12°C. The city receives an average annual precipitation of 550–600 mm, with a pronounced peak during the summer months, particularly from June to September, accounting for approximately 75% of the total annual rainfall(Hu et al. 2013). Notably, areas with severe ground subsidence in Tianjin are concentrated in the central districts (Dongli, Xiqing, Jinnan, and Beichen), as well as in the Binhai New Area, Wuqing, and Jinghai districts.
The Tianjin Jinnan district, covering an area of approximately 434.5 square kilometers, encompasses six districts at the municipal level and seven sub-districts at the village level. The area has a permanent population of around 924,300. The urbanization rate, which represents the proportion of the population residing in urban areas, is approximately 80%(Hu et al. 2013). Tianjin Jinnan lies in a plain region, bordered to the north by the Tianjin Binhai New Area. Its terrain is relatively low, with the lowest point at an elevation of approximately − 8 meters and the highest point reaching around 27 meters(Lixin et al. 2011). The Tianjin Metro Lines 1 and 6 traverse this area. Economically, Tianjin Jinnan is a crucial development zone within Tianjin. It focuses on the development of modern urban agriculture, circular economy, and ecotourism. However, due to its geographical location and geological characteristics, the area faces specific environmental and geological challenges, such as soil erosion and subsidence, which must be carefully considered in regional planning and construction.
2.2 Data sources
This study leveraged high-resolution radar imagery data from the Sentinel-1A satellite, acquired between January 9, 2016, and July 31, 2023. The data covered a time series of deformations in the study area, with a resolution of 5 meters by 20 meters. The imagery was obtained in an ascending orbit, using VV polarization and IW imaging mode(alaska.edu). Notably, Sentinel-1A data differs from conventional remote sensing data due to its inverted SAR coordinates and geographic coordinates. For digital elevation model (DEM) data, we utilized the 90-meter resolution STRM (Shuttle Radar Topography Mission) model provided by NASA(dwtkns.com). Additionally, we incorporated precise orbit determination (POD) data from the European Space Agency (ESA), which boasts a positioning accuracy of less than 5 centimeters(esa.int). To comprehensively analyze the relationship between ground subsidence and urbanization, we also gathered rainfall data from NOAA(noaa.gov), Tianjin road network data, historical optical images from Google Earth, Sentinel-2A imagery(copernicus.eu), and on-site photography data(sohu.com). This integrated approach allows us to explore the intricate interplay between ground deformation and urban development in the study region.
3.Methods
3.1 SBAS-InSAR Technique for Ground Subsidence Detection
The Small Baseline Subset (SBAS) technique was a multi-master interferometric synthetic aperture radar (InSAR) method(Liu and Zhang 2023). It leveraged shorter temporal and spatial baselines to acquire a greater number of interferograms(Du et al. 2021). This approach helped maintain consistency in the interferometric phase information. Unlike traditional Persistent Scatterer InSAR (PS-InSAR), SBAS overcame issues related to low coherence caused by a single dominant scatterer, thereby reducing the SAR image requirements. SBAS-InSAR found extensive application in ground subsidence detection. It was based on the following equation(Cigna et al. 2021):
$$\begin{array}{c}\frac{N+1}{2}\le M\le \frac{N\left(N+1\right)}{2}\end{array}$$
1
$$\delta {\phi }_{j}\left(x,\gamma \right)=\phi \left({t}_{B},x,r\right)-\phi \left({t}_{A},x,r\right)\approx {\varDelta \phi }_{disp}+{\varDelta \phi }_{topo}+{\varDelta \phi }_{orb}+{\varDelta \phi }_{atm}+{\varDelta \phi }_{noise}$$
2
Equation (1) represented the range of M differential interferograms generated from N + 1 ordered Synthetic Aperture Radar (SAR) images acquired over time at \({t}_{0},{t}_{1},\dots ,{t}_{N}\) within the same region. Additionally, Eq. (2) described the composition of the interferometric phase in pixel (x, r) for interferogram j (generated from two images at tB and tA). Here, x and r represented azimuth and distance coordinates, respectively. The phase \({\varDelta {\phi }}_{\text{d}\text{i}\text{s}\text{p}}\) arose from distance variations along the Line-of-Sight (LOS) between the target and the radar. Furthermore, the phases \({\varDelta \phi }_{topo}\), \({\varDelta \phi }_{orb}\), \({\varDelta \phi }_{atm}\), and \({\varDelta \phi }_{noise}\) originated from residual topography, satellite orbit errors, atmospheric effects (particularly tropospheric delay), and other noise sources. By removing the residual components from the interferometric phase using SBAS-InSAR, we obtained the deformation phase \({\varDelta \phi }_{disp}\). Subsequently, a system of M equations with N unknowns can be organized into the following matrix:
$$\begin{array}{c}A\phi =\delta \phi \end{array}$$
3
Given a coefficient matrix A of dimensions M×N,\(\phi ={\left[{\phi }\left({\text{t}}_{1}\right),{\phi }\left({\text{t}}_{2}\right),\dots ,{\phi }\left({\text{t}}_{\text{N}}\right)\right]}^{\text{T}}\) represented the unknown deformation phase values as an N×1 vector, and \({\delta }{\phi }={\left[{\delta }{\phi }\left({\text{t}}_{1}\right),{\delta }{\phi }\left({\text{t}}_{2}\right),\dots ,{\delta }{\phi }\left({\text{t}}_{\text{N}}\right)\right]}^{\text{T}}\) denoted the unwrapped phase values as another N×1 vector, we can rewrote Eq. (3) as follows to obtain the deformation velocity:
$$\begin{array}{c}Bv=\delta \phi \end{array}$$
4
Given a coefficient matrix B of dimensions M×N, where v provided the average phase velocity as follows:
$$\begin{array}{c}v={\left[{\text{v}}_{1},{\text{v}}_{2},{\dots ,\text{v}}_{\text{N}}\right]}^{\text{T}}={\left[\frac{{{\phi }}_{1}}{{\text{t}}_{1}-{\text{t}}_{0}},\frac{{{\phi }}_{2}-{{\phi }}_{1}}{{\text{t}}_{2}-{\text{t}}_{1}},\dots ,\frac{{{\phi }}_{\text{N}}-{{\phi }}_{\text{N}-1}}{{\text{t}}_{\text{N}}{-\text{t}}_{\text{N}-1}}\right]}^{\text{T}}\end{array}$$
5
we can apply the least squares (LS) or singular value decomposition (SVD) algorithms to equations (4) and (5) to compute the deformation velocity(Mora et al. 2003). By performing calculations over specific time intervals, cumulative deformation can be generated.
3.2 SBAS-InSAR data processing process
In this investigation, we employed the SBAS approach within the InSAR framework to analyze ground deformation(Fig. 2). The process began with the acquisition of SAR Single Look Complex (SLC) imagery, followed by a careful selection of a small baseline subset to optimize interferometric analysis. Interferogram generation, a critical phase, utilized an external Digital Elevation Model (DEM) for data smoothing, filtering, phase unwrapping, satellite trajectory refinement, and phase discrepancy correction. Ground Control Points (GCPs) were strategically chosen to enhance the accuracy of slope phase effect elimination and inversion robustness. The methodology progressed through two pivotal inversion stages. The initial inversion targeted atmospheric phase influence reduction, deriving average deformation rates and DEM correction coefficients. The subsequent inversion refined these metrics, improving deformation rate and DEM coefficient accuracy.
After inversion, we transformed the cumulative deformation and rate data from SAR to geographic coordinates, aiding geospatial analysis application. Refinement of the SBAS-InSAR method involved generating interferogram pairs with stringent temporal and spatial baselines to reduce noise. We applied Goldstein’s adaptive filtering algorithm for phase filtering and set a multi-looking ratio of 4:1. Phase unwrapping used the minimum cost flow method with a coherence threshold of 0.3. Orbit refinement and phase removal were conducted using 34 high-coherence GCPs, avoiding areas of significant deformation, to precisely estimate and correct residual phases and slopes. The two-step inversion process first calculated displacement rates and residual topography, then computed time-series displacement, applying atmospheric filtering tailored to initial deformation rates for atmospheric phase removal. Finally, we geocoded the inversion outcomes, converting phase data into elevation information aligned with the study’s goals. This involved using radar LOS data and the DEM to generate LOS and vertical displacement maps in the GCS-WGS-84 coordinate system, offering an exhaustive ground deformation analysis.
3.3 PS-InSAR Technique for Ground Subsidence Detection
PS-InSAR technology, first proposed in 2000, addresses limitations in conventional InSAR caused by factors such as atmospheric delay, orbital errors, and terrain effects, which significantly restricted its applicability(Ferretti et al. 2000). The fundamental principle of PS-InSAR involved identifying Persistent Scatterer (PS) points that maintained stable coherence across multiple SAR images with long temporal and spatial baselines. These PS points predominantly occurred on artificial structures and exposed rock surfaces. Leveraging the strong coherence of PS points, PS-InSAR achieved true surface deformation information at millimeter precision, effectively mitigating the impact of temporal and spatial inconsistencies on interferometric phase quality. Urban surface monitoring benefited from PS-InSAR due to the abundance of buildings in cities. Most PS points were selected from fixed building corners (permanent scatterers) during technical processing. With the ability to detect displacements at the millimeter level and derive surface deformation rates from multi-temporal data, PS-InSAR proved valuable for urban management and decision-making(Zhang et al. 2023b). In summary, PS-InSAR technology provided a powerful tool for monitoring urban surface deformation, enabling precise measurements and efficient analysis of subsidence rates over time(Zhang et al. 2023c).
In the context of PS-InSAR technology, the phase resulting from deformation can be separated using the following equation based on imaging geometry:
$$\begin{array}{c}\varPhi ={\phi }_{deformation}+{\phi }_{topograpℎic}+{\phi }_{orbit}+{\phi }_{atmospℎeric}+{\phi }_{noise}\end{array}$$
6
The interferometric phase of two SAR images consisted of several components: differential deformation phase due to tilt, residual topographic phase, orbit inaccuracies, atmospheric delay phase, and thermal noise from the radar system. After removing other components, we obtained the deformation phase. Subsequently, we calculated the surface deformation rate along the LOS. Considering that the horizontal deformation in the Back Projection (BP) method can be neglected, we estimated the ground subsidence rate using the following formula(Zhou et al. 2019):
$$\begin{array}{c}v=\frac{{V}_{LOS}}{cos\theta }\end{array}$$
7
Here, we defined the following terms: V represented the ground subsidence rate (vertical deformation rate), \({V}_{LOS}\) denoted the deformation rate in the LOS direction, and θ signified the incidence angle of the radar sensor.
3.4 PS-InSAR data processing process
The study employed PS-InSAR, a powerful technique for monitoring surface deformation, by analyzing changes in surface points over time. To ensure high precision and accuracy, specific guidelines were followed. A robust dataset of 193 Sentinel-1A images, spanning January 2016 to July 2023, was selected, with at least 20 scenes and continuous time intervals of approximately one month between adjacent scenes. The PS-InSAR implementation involved five key stages: 1) network construction by linking 193 Sentinel-1A images covering Tianjin Jinnan district from 2016 to 2023, generating 668 interferometric pairs with December 19, 2019, as the master date; 2) interferometric processing, including georeferencing, interferogram generation, phase unwrapping, and amplitude deviation index computation; 3) initial model inversion to obtain displacement rates, residual topography, and multi-temporal coherence coefficients; 4) second model inversion to estimate and remove atmospheric phase components; and 5) geocoding to derive intensity data, surface deformation rates, deformation precision, and elevation precision from January 2016 to July 2023. This rigorous approach facilitated accurate and continuous ground deformation monitoring using PS-InSAR technology.
3.5 Machine Learning Models for Ground Subsidence Prediction
The study employed five machine learning methods to construct predictive models for ground subsidence rates. The Support Vector Machine (SVM) method, grounded in statistical theory, identified the optimal hyperplane that separated the training dataset while maximizing geometric distance(Dagès et al. 2023). The Gradient Boosting Decision Tree (GBDT) algorithm iteratively adjusted the weights of base learners to minimize the loss function, combining decision trees and gradient boosting techniques (Zhou et al. 2019). The Random Forest (RF) method aggregated predictions from multiple base learners in an ensemble approach(Zhang et al. 2023a). The Extremely Randomized Trees (ERT) algorithm extended from RF by introducing greater randomness in feature selection and node splitting thresholds(Shen et al. 2023a). Lastly, the Long Short-Term Memory (LSTM) model, a deep learning approach based on recurrent neural networks, effectively captured long-term dependencies in sequential data through memory units and gating mechanisms(Chen et al. 2021; Ma et al. 2023). These models were rigorously validated through cross-validation techniques.
After iterative parameter tuning, the LSTM model was determined to be the optimal approach for this application. The LSTM model utilized 509,201 InSAR-detected pixels to predict cumulative subsidence data from December 19, 2019, to July 20, 2025. A test set of 56 predictions from December 19, 2019, to July 31, 2023, validated the accuracy of the model, while predictions from August 24, 2023, to July 20, 2025, provided the final ground subsidence forecast for the Jinnan District. The LSTM approach was further employed to establish ground subsidence prediction models for four selected regions with significant deformation rates: building areas (P1, P2), along highways (P3), and along subway lines (P4). The model's time step was set to 5, utilizing cumulative subsidence data from 1 to 5 monitoring periods to predict the 6th period, and so on. A total of 191 cumulative subsidence data points from December 19, 2019, to March 17, 2026, were predicted for these four regions. The 111 predictions from December 19, 2019, to July 31, 2023, served as the test set for accuracy validation, while the remaining period (July 31, 2023, to March 17, 2026) represented future cumulative subsidence predictions. Predictions for specific dates during both the rainy and dry seasons, spanning from July 31, 2023, to March 17, 2026, were showcased.
3.6 Comparative Evaluation of SBAS-InSAR and PS-InSAR Techniques in Monitoring Ground Subsidence
Ground subsidence data in the study area was monitored from 2016 to 2023 using the SBAS-InSAR technique. Due to the unavailability of contemporaneous horizontal measurement data for comparison, the PS-InSAR method was employed to validate the subsidence values obtained from SBAS-InSAR reference points. The error assessment involved determining the coefficient of determination (R2) and applying a linear fitting function to the average deformation rates derived from both SBAS-InSAR and PS-InSAR techniques.
$$\begin{array}{c}{R}^{2}=1-\frac{{SS}_{err}}{{SS}_{tot}}\\ =1-\frac{{\sum }_{i=1}^{m}{\left({y}_{i}-\widehat{{y}_{i}}\right)}^{2}}{{\sum }_{i=1}^{m}{\left({y}_{i}-\stackrel{-}{{y}_{i}}\right)}^{2}}\sqrt{\frac{{\sum }_{i=1}^{N}{\left({X}_{i}-\stackrel{-}{X}\right)}^{2}}{N}}\end{array}$$
8
Where, the coefficient of determination, denoted as R2, quantifies the proportion of the variance in the dependent variable that is predictable from the independent variables. It is calculated by comparing the sum of squared errors, or residuals \({SS}_{err}\), which is the sum of the squared differences between the observed values and the predicted values\({\left({y}_{i}-\widehat{{y}_{i}}\right)}^{2}\) for each observation i. To the total sum of squares \({SS}_{tot}\), which represents the total variance in the observed data from their mean, calculated as the sum of squared differences between each observed value \({X}_{i}\) and the mean value \({\left({X}_{i}-\stackrel{-}{X}\right)}^{2}\) for each data point i. Here, m denotes the number of observations, and N represents the number of data points.
To quantitatively evaluate the spatiotemporal predictive performance of each model, we employed the following metrics: coefficient of determination (R2), mean absolute error (MAE), root mean square error (RMSE), and mean absolute percentage error (MAPE).
$$\begin{array}{c}MAE=\frac{1}{m}\sum _{i=1}^{m}\left|{y}_{i}-\widehat{{y}_{i}}\right|\end{array}$$
9
Where m is the total number of data points, \({y}_{i}\) is the actual observed value, and \(\widehat{{y}_{i}}\) is the predicted value. This equation calculates the Mean Absolute Error (MAE), which is the average of the absolute differences between the predicted and observed values.
$$\begin{array}{c}RMSE=\frac{1}{m}{\sum }_{i=1}^{m}{\left({y}_{i}-\widehat{{y}_{i}}\right)}^{2}\end{array}$$
10
Where m is the total number of data points, \({y}_{i}\) is the actual observed value, and \(\widehat{{y}_{i}}\) is the predicted value. This equation calculates the Root Mean Square Error (RMSE), which is the square root of the average of the squared differences between the predicted and observed values.
$$\begin{array}{c}MAPE=\frac{1}{m}{\sum }_{i=1}^{m}|\frac{{y}_{i}-\widehat{{y}_{i}}}{{y}_{i}}|\end{array}$$
11
Where m is the total number of data points, \({y}_{i}\) is the actual observed value, and \(\widehat{{y}_{i}}\) is the predicted value. This equation calculates the Mean Absolute Percentage Error (MAPE), which is the average of the absolute differences between the predicted and observed values, expressed as a percentage of the observed values.