Large-Scale Surface Deformation Monitoring Using SBAS-InSAR and Intelligent Prediction in Typical Cities of Yangtze River Delta

Abstract

Large-scale short-term monitoring and prediction of surface deformation are of great significance for the prevention and control of geohazards in rapidly urbanizing developing cities. Most studies focus on individual cities, but it would be more meaningful to address large urban agglomerations and consider the relevance of the regions within them. In addition, the commonly used linear fitting prediction methods cannot accurately capture the dynamic mechanisms of deformation. In this study, we proposed an automatic PS extraction method (named PS-SBAS-InSAR) that improves SBAS-InSAR to extract surface deformation and an Informer-based short-term surface deformation prediction method for case studies in 16 typical cities of the Yangtze River Delta (YRD). The results show that PS-SBAS-InSAR successfully extracted accurate surface deformation sequences of the YRD. During the period from January 2019 to January 2021, the YRD experienced a slight deformation with an average deformation rate within [−4, 4] mm/year. Geographically neighboring cities may have associated deformation distributions and similar deformation trends, as indicated by average deformation rate maps and landscape metrics. Both types of deformation (i.e., subsidence/uplift) tend to occur simultaneously, with specific areas of subsidence/uplift occurring in close proximity to areas of concentrated deformation. The Informer model effectively captured the time-series variation in surface deformation, suggesting a slowdown of deformation over the next two months (February 2021–March 2021). Our work contributes to a better understanding of changes and trends in large-scale surface deformation and provides useful methods for monitoring and predicting surface deformation in coastal areas.

1. Introduction

Surface deformation is a slowly changing environmental geohazard caused by both natural conditions and human activities [1,2]. More than 150 cities and regions worldwide are experiencing or have experienced surface deformation, with about 80% occurring in riverbeds, floodplains, or deltas [3,4,5]. Ground subsidence in deltas, when combined with sea level rise, heavy rainfall (e.g., monsoons), or storms, could worsen the risk of storm surges and tidal inundation, causing great economic losses [5,6,7]. Conventional engineering surveys are generally used to monitor the deformation of individual projects or buildings, but they are time-consuming and labor-intensive for large-scale monitoring such as monitoring the Yangtze River Delta (YRD). Therefore, it is necessary to integrate satellite remote sensing and deep learning to monitor and predict the surface deformation in large-scale areas.

Interferometric synthetic aperture radar (InSAR) extracts ground elevation and deformation information by acquiring phase information from coherent images [8,9]. In comparison to other methods (e.g., leveling and real-time kinematic equipment), InSAR technology offers significant advantages, such as a broad monitoring range, low operating costs, and a high spatial resolution [7,10,11]. The time-series InSAR (TS-InSAR) technique is an advancement of the traditional InSAR technique that uses more observations and an improved signal-to-noise ratio to reconstruct the deformation history of a deformation region [12]. Common TS-InSAR methods include permanent scatterer InSAR, small baseline subset InSAR (SBAS-InSAR), and temporally coherent point InSAR for deformation monitoring in different scenarios [13,14].

SBAS-InSAR technology attenuates the decoherence caused by excessively long spatiotemporal baselines by defining optimal baseline thresholds and segmenting a whole SAR dataset into several smaller datasets [14,15]. The SBAS-InSAR method imposes stringent criteria for identifying ground control points (GCPs), requiring a high degree of consistency across all GCPs. However, the geography of a coastal/reclaimed area with dense vegetation and numerous lakes presents challenges in finding stable and highly coherent GCPs through manual visual interpretation. In general, GCPs can be selected [16] from relevant targets that maintain backscattering characteristics for a long period (i.e., permanent scatterers, PSs). In 2022, Feng et al. [17] developed a new method (auto-PS-GCP-Thin) to automatically extract robust GCPs from candidate PSs, enhancing the efficiency and accuracy of GCP selection. We thus utilized this method to automatically select GCPs and applied the GCPs to SBAS to develop a new deformation detection method (PS-SBAS-InSAR), aiming to solve the insufficient and unreliable samples in coastal lowlands and improve subsidence monitoring accuracy.

Long-term surface deformation is not only affected by urban construction but also by crustal movement and other human activities, and long-term prediction has rather high uncertainties, so short-term prediction has commonly been used to support the analysis of urban construction and development. Current short-term prediction methods for surface deformation can generally be divided into two categories. The first category estimates the deformation value by simulating the physical evolution process [18,19], which is limited by the difficulty of acquiring parameters and constructing models. The second category applies mathematical and statistical methods to predict ground deformation. This is performed by extracting the intrinsic relationships and development patterns from extensive historical monitoring datasets [20,21]. Currently, machine learning and deep learning methods (e.g., neural networks) can effectively address the highly complex nonlinear relationships involved in surface deformation prediction.

The Transformer model is a neural network that relies entirely on an attention mechanism to characterize the global dependencies between the inputs and outputs. It has a great capability to capture short- and long-term dependencies and interactions [22]. The Informer model is an improvement on the Transformer model developed for time-series prediction tasks [23]. The Informer model proposes a multi-headed ProbSparse self-attention layer based on the standard encoder–decoder structure. This improves model efficiency and reduces the number of model parameters that can be applied to power load and traffic flow forecasting [23]. Nevertheless, the potential of the Informer model for time-series trend prediction of surface deformation has not been elaborated.

The YRD is a region dominated by impact plains [24,25], and studies have found that its geological structure is not sufficiently stable or geologically strong [10,26,27]. Since cities in the YRD are part of the same hydrogeologic unit and the groundwater aquifer systems are interconnected, surface deformation in one area can have a significant impact on neighboring areas. Surface deformation in the YRD varies with the urbanization process and has been influenced by industrialization and the prevention and control policies in different periods, showing stage-specific characteristics [28,29]. Although there have been some studies on the surface deformation of individual cities in the YRD [12,30], surface deformation and its regional correlation at large scales in the urban agglomeration have not been sufficiently elaborated.

Large-scale short-term monitoring and prediction of surface deformation can provide timely and important information that can help to identify potential problems early and reduce hazard risk, ensuring the sustainable development of urban areas. This study addresses three questions: (1) Can the automatically extracted PS points be applied to SBAS to extract urban surface deformation? (2) Is there a regional correlation of surface deformation in a large-scale urban agglomeration? (3) Is it possible to use deep learning methods for the short-term prediction of surface deformation? To answer these questions, we selected 16 typical cities in the YRD, used the PS-SBAS-InSAR technology to monitor the surface deformation of each city from January 2019 to January 2021, measured the spatial pattern of the deformation using six landscape metrics at the class level, and then applied the Informer model to predict the deformation trend for the next 2 months. Our study should provide methodological references for monitoring and predicting surface deformation in similar coastal regions.

2. Materials and Methods

Figure 1 illustrates the main technical processes: (a) data pre-processing, including remote sensing image clipping of the study area as well as the atmospheric correction and radiometric calibration of the SAR images using ENVI tools to ensure that the amplitudes of all images correctly reflect the scattering from ground targets; (b) PS-SBAS InSAR processing, including using the auto-PS-GCP-Thin method to intelligently identify highly coherent GCPs and applying the GCPs to the SBAS process to calculate the spatiotemporal variation in regional deformation; and (c) using Informer models to predict short-term future surface subsidence and uplift in 16 typical cities.

2.1. The PS-SBAS-InSAR Technology

Conventional InSAR technology mainly performs coherent processing of radar signals and uses the relationship between phase/spatial distance differences in two observations to extract topographic and surface deformation information [8]. The interference phase (${\varphi }_{int}$) is given by

$${\varphi }_{def}={\varphi }_1-{\varphi }_2=-\frac{2\pi }{\lambda }\left(R_1-R_2\right)=-\frac{2\pi }{\lambda }\mathsf{\Delta }R$$

(1)

where ${\varphi }_1$ and ${\varphi }_2$ are the phases of the echo signals from the two observations of the same target; $\lambda$ is the radar wavelength; R₁ and R₂ are the distances from the two antennas to the observation target; and ΔR is the difference between the distances.

The Differential InSAR (D-InSAR) technique requires two SAR images of the same area in different time periods and removes the constant phase from the two observed phases with differential processing to derive the surface deformation information. The acquisition of the deformation phase (${\varphi }_{def}$) using the D-InSAR technique can be defined as [14]

$${\varphi }_{def}=\varphi -{\varphi }_{flat}-{\varphi }_{dem}-{\varphi }_{atm}-{\varphi }_{noise}$$

(2)

where ${\varphi }_{flat}$ is the flat-earth phase (i.e., reference surface phase); ${\varphi }_{dem}$ is the terrain phase; ${\varphi }_{atm}$ is the atmospheric phase; and ${\varphi }_{noise}$ is the noise phase. Usually, the D-InSAR technique extracts deformation information that is not only affected by the constant phase error. There may also be spatiotemporal decoherence phenomena, resulting in insufficient extracted deformation information [1,10].

PS-SBAS-InSAR is divided into two modules: a PS module and an SBAS module. The PS module considers temporal stability and spatial continuity to select a suitable number of PSs and then identifies the GCPs by automatically evaluating the coherence, average amplitude, and amplitude deviation index (ADI) of the PSs [17]. The evaluation of the coherence coefficient, average amplitude, and ADI are defined as

$$\left\{\begin{array}{c}standard.\gamma :\left\{\begin{array}{c}{\gamma }_{pix}=\frac{1}{K}{\displaystyle \sum }_1^K\gamma \\ s.t.{\gamma }_{pix}>{\gamma }_{pix-thd}\end{array}\right.\\ standard.A:\left\{\begin{array}{c}{A}_{img-thd=min\left\{{\displaystyle \sum }_1^K{A}_{img}\right\}}\\ s.t.A_{pix}>A_{img-thd}\end{array}\right.\\ standard.ADI:\left\{\begin{array}{c}ADI=\frac{A_{pix-\sigma }}{A_{pix}}\\ s.t.ADI<AD{I}_{thd}\end{array}\right.\end{array}\right.$$

(3)

where ${\gamma }_{pix}$ is the average coherence of pixels at an identical location across K time-series images; ${\gamma }_{pix-thd}$ is a pre-defined coherence threshold and is defined here as 0.4 [31]; $A_{img}$ is the average amplitude of an SAR image; $A_{img-thd}$ is the minimum value of the average amplitude in the K time series and can use $A_{img-thd}$ as the threshold to eliminate the points with low amplitudes but high coherence; $A_{pix}$ is the average amplitude of a pixel at the same position in K time-series images; $A_{pix-\sigma }$ is the standard deviation ($\sigma$) of the amplitude of a pixel in K time-series images; and $AD{I}_{thd}$ is the pre-defined ADI threshold value, which is set to 0.15 here [32]. After these three screening steps, the highly consistent and strongly stable GCP points were derived. Then, the SBAS method was applied to calculate the time-series deformation and the annual average deformation rate.

The SBAS module is used to filter and collect combinations with higher quality from interferograms, and then the singular value decomposition (SVD) method is used for phase unwrapping according to the relationship between phase and acquisition time, calculating the surface deformation sequence for the whole period [33]. According to the interferometric combination condition, a free combination of N + 1 SAR images of the same region can produce M interferograms, where M satisfies the following [14]:

$$\frac{N+1}{2}\le M\le N\left(\frac{N+1}{2}\right)$$

(4)

Assuming that the SAR images at moments t_A and t_B (t_A < t_B) generate the j-th interferogram, after removing the common quantities in the two observed phases using the previously acquired GCPs (Equation (2)), the differential interference phase at (x, r) in azimuthal and distance coordinates can be expressed as [14]

$$\delta {\varphi }_j\left(x,r\right)=\varphi \left(t_B,x,r\right)-\varphi \left(t_A,x,r\right)\approx \frac{4\pi }{\lambda }\left[d\left(t_B,x,r\right)-d\left(t_A,x,r\right)\right]$$

(5)

where $\varphi \left(t_A,x,r\right)$ and $\varphi \left(t_B,x,r\right)$ are the phase of the pixel point (x, r) at moments t_A and t_B, respectively, and $d\left(t_A,x,r\right)$ and $d\left(t_B,x,r\right)$ are the displacement information in the radar line-of-sight direction (RLOS) relative to the reference moment $d\left(t_0,x,r\right)$ = 0 at moments t_A and t_B, respectively. In this study, we followed the literature [6,26] to define the positive and negative displacements as uplift and subsidence, respectively. In addition, the two deformation types were assumed to be vertical, and horizontal deformation was neglected.

The specific steps for PS-SBAS-InSAR processing (Figure 2) include (1) defining spatiotemporal baseline thresholds to filter image pairs for matching; (2) differential interference processing to generate interferograms; (3) conducting phase wrapping of the interferograms with flattening using the Delaunay 3D method; (4) constructing the observational equations based on the GCPs; (5) using the SVD method to solve the equations, then calculating the annual average surface deformation rate; and (6) separating the modeling residuals and deriving the time-series surface deformation.

Many studies [33,34] have used intra-accuracy analysis to validate a method’s performance and results when field monitoring data are not available. The estimated residuals of the average surface deformation rate are calculated, and smaller residuals imply better model performance.

2.2. The Landscape Metrics

Landscape metrics are typically used to quantitatively reflect the structural composition and characteristics of spatial patterns and can be calculated using Fragstats software [34]. Usually, landscape metrics can be divided into path, class, and landscape. We selected six class-level landscape metrics (

Table 1. Remote sensing datasets of typical cities in the YRD.

Metric	Abbreviation	Category	Description
Largest patch index	LPI	Area-Edge	Indicates the percentage of the largest patches of a given type over the entire landscape area.
Edge density	ED	Area-Edge	Refers to the length of the edge between patches of landscape elements.
Perimeter–area fractal dimension index	PAFRAC	Shape	Reflects the complexity of the patch shape.
Landscape shape index	LSI	Aggregation	Reflects the complexity of the landscape shape.
Landscape division index	DIVISION	Aggregation	Represents the degree of separation between patches of a certain class.
Aggregation index	AI	Aggregation	Indicates connectivity between patches of the same class.

) to further analyze the spatial pattern of the deformation fields [35,36,37].

2.3. The Informer Model

For univariate time-series forecasting, it is usually necessary to find a suitable function to establish the relationship between observations and time covariates, which can be expressed as [38]

$$\left\{\begin{array}{c}{v}_t~f\left(Y_t\right)\\ y_t\triangleq \left[d_{t-1}\circ x_t\right]\\ Y_t={\left[y_1,\dots ,y_t\right]}^T\end{array}\right.$$

(6)

where $v_t$ (t = 1, 2, …) is the observed value at time step t, which is the cumulative deformation rate in this study; $x_t$ (t = 1, 2, …) is a vector of time-based covariates that are assumed to be known throughout a time period, such as a day of the week or an hour of the day; and $y_t$ is an augmentation matrix where $\left[\cdot \circ \cdot \right]$ is the concatenation.

In the traditional attentive layer, the multi-headed self-attentive sub-layer simultaneously transforms Y into query matrix $Q=Y{W}^Q$, key matrix $K=Y{W}^K$, and value matrix $V=Y{W}^V$, where $W^Q,W^K$, and $W^K$ are the weight matrices. The self-attention (A) is defined based on the tuple inputs $A\left(Q,K,V\right)$, and the i-th query’s attention is defined as [22]

$$\left\{\begin{array}{c}A\left(Q,K,V\right)=Softmax\left(\frac{Q{K}^T}{\sqrt{d}}\right)V\\ A\left(q_i,K,V\right)={\displaystyle \sum }_{j}p\left(\left.k_j\right|q_i\right)v_j\end{array}\right.$$

(7)

where $q_i$, $k_i$, and $v_i$ are the i-th rows of the Q, K, and V matrices, respectively, and $p\left(\left.k_j\right|q_i\right)$ is the distribution of keys given the condition of query i. For multi-head self-attention, each head takes the same optimization strategy [22]. However, an element in the sequence will generally only have high similarity/association with a few elements, and the involvement of each element in the computation leads to substantial computational complexity. The Informer model proposes a ProbSparse self-attention from a probabilistic view, which applies a criterion to measure the sparsity of self-attention. This criterion measures the difference between the attention probability distribution (p) and the uniform distribution (q) using Kullback–Leibler (KL) divergence and screens out the most important queries. The measure of self-attention sparsity can be given as [23]

$$\left\{\begin{array}{c}Attentionprobability:p\left(\left.k_j\right|q_i\right)=\frac{\exp \left(\frac{q_i{k}_j{}^T}{\sqrt{d}}\right)}{{{\displaystyle \sum }}_l\exp \left(\frac{q_i{k}_l{}^T}{\sqrt{d}}\right)}\hfill \\ Uniformprobability:q\left(\left.k_j\right|q_i\right)=\frac{1}{L_k}\hfill \\ KL\left(q\left|\right|p\right)={\displaystyle \sum }_{j}qln\left(\frac{q}{p}\right)\hfill \end{array}\right.$$

(8)

If $p\left(\left.k_j\right|q_i\right)$ is close to $q\left(\left.k_j\right|q_i\right)$, it indicates that query i cannot select the important key and will be defined as ‘lazy’; conversely, it will be defined as ‘active’. These active queries contribute more to self-attention and need to be selected to participate in dot product operations. The similarity of $p\left(\left.k_j\right|q_i\right)$ and $q\left(\left.k_j\right|q_i\right)$ can be measured using KL divergence. The larger the value of query sparsity, the greater the difference between the probability distribution and the uniform distribution and the more active the query is. Therefore, the self-attentiveness is further optimized in Informer as [23]

$$A\left(Q,K,V\right)=Softmax\left(\frac{\overline{Q}{K}^T}{\sqrt{d}}\right)V$$

(9)

The evaluation of sparsity followed by dot product operation can greatly reduce the computational complexity and thus improve the efficiency of timing prediction. In addition, the Informer model proposes self-attention distilling for down-sampling operations to reduce the number of dimensions and network parameters, and it proposes a generative style decoder for all predictions in one step. These optimizations make the Informer model more suitable for the simulation and prediction of long time series. Here, the Informer model is only applied to achieve the prediction of future surface deformation rates, and more details about this model can be found in a previous study [23].

In addition, we applied four common metrics to measure the ability of the Informer model to capture deformation time-series information, including the Mean Square Error (MSE), Root-Mean-Square Error (RMSE), Mean Absolute Error (MAE), and Pearson product-moment correlation coefficient (PCC).

3. Study Area and Datasets

The YRD (Figure 3) is located at the intersection of coastal and riverine economic zones, encompassing a land area of 358,000 km². By the end of 2019, the regional population was about 227 million, with more than 60 percent being urban residents.

The urban spatial pattern of the YRD can be described as “one core, five circles, and four belts” [39]. “One core” refers to the core city of Shanghai; the “five circles” include Nanjing, Hangzhou, Hefei, Su-Xi-Chang, and the Ningbo metropolitan area; and the “four belts” refer to the four economic belts through which the major economic corridor runs right through these five cities (Figure 3). The YRD has high urbanization rates and vulnerable geological environments, resulting in a region that is experiencing surface deformation [7,10]. The surface deformation is strongly correlated with the urbanization rate and intensity. Therefore, we selected 16 important cities across these five metropolitan areas (excluding Zhoushan) as representatives. City centers are usually situated on flat terrain and densely built up, so high-quality interference images can be acquired readily. Therefore, we selected the central urban parts of these cities as our study areas and identified the boundaries of each urban center according to local land-use planning departments. The ascending SAR data for each city were collected from the ESA’s Sentinel-1A single-observation composite data for a one-month interval from January 2019 to January 2021, totaling 25 images (search.asf.alaska.edu (accessed on 13 August 2023)).

Table 2. Remote sensing datasets of typical cities in the YRD.

Category	City	Abbreviation	Sentinel-1A Footprint	Angle of Incidence
Core city	Shanghai	SH	171-96	39.032°
Nanjing metropolitan area	Nanjing	NJ	69-99	39.037°
	Zhenjiang	ZJ	69-99	39.037°
	Yangzhou	YZ	69-99	39.037°
Hangzhou metropolitan area	Hangzhou	HZ	69-94	39.083°
	Jiaxing	JX	69-94	39.083°
	Huzhou	HuZ	69-94	39.083°
	Shaoxing	SX	69-94	39.083°
Hefei metropolitan area	Hefei	HF	142-101	39.007°
	Wuhu	WH	142-96	39.072°
	Maanshan	Ma	69-99	39.037°
Su-Xi-Chang metropolitan area	Suzhou	SZ	69-99	39.037°
	Wuxi	WX	69-99	39.037°
	Changzhou	CZ	69-99	39.037°
Ningbo metropolitan area	Ningbo	NB	171-91	38.867°
	Taizhou	TZ	171-91	38.867°

For the “Sentinel-1A footprint’ item, for example, for 171-96, 171 indicates the orbital number and 96 indicates the map width number.

shows the image information for the 16 cities, while

Table 3. The acquisition dates for different orbital data.

Image No.	Acquisition Date
	Orbital No. 69	Orbital No. 142	Orbital No. 171
1	17 January 2019	22 January 2019	24 January 2019
2	22 February 2019	27 February 2019	17 February 2019
3	18 March 2019	23 March 2019	13 March 2019
4	11 April 2019	16 April 2019	6 April 2019
5	29 May 2019	22 May 2019	24 May 2019
6	22 June 2019	27 June 2019	17 June 2019
7	16 July 2019	21 July 2019	11 July 2019
8	21 August 2019	26 August 2019	16 August 2019
9	26 September 2019	19 September 2019	21 September 2019
10	20 October 2019	25 October 2019	15 October 2019
11	13 November 2019	18 November 2019	8 November 2019
12	19 December 2019	24 December 2019	14 December 2019
13	12 January 2020	17 January 2020	7 January 2020
14	17 February 2020	22 February 2020	12 February 2020
15	24 March 2020	29 March 2020	19 March 2020
16	17 April 2020	22 April 2020	12 April 2020
17	23 May 2020	28 May 2020	18 May 2020
18	16 June 2020	21 June 2020	11 June 2020
19	22 July 2020	27 July 2020	17 July 2020
20	15 August 2020	20 August 2020	10 August 2020
21	20 September 2020	25 September 2020	15 September 2020
22	14 October 2020	19 October 2020	9 October 2020
23	7 November 2020	12 November 2020	2 November 2020
24	13 December 2020	18 December 2020	8 December 2020
25	18 January 2021	23 January 2021	13 January 2021

shows each orbit imaged for no more than 10 days in the same period (one month). In addition, we applied the DEM elevation data to eliminate the topographic phases during the PS-SBAS-InSAR processing, where the 90 m resolution STRM DEM was acquired from the Geospatial Data Cloud (www.gscloud.cn (accessed on 13 August 2023)).

4. Results

4.1. Deformation Results and Assessment

A core concept of SBAS is to extract surface deformation information using only interferometric pairs with short spatiotemporal baselines. We pre-processed the SAR images by defining thresholds of 45% for spatial baselines and 200 days for temporal baselines to segment datasets into subsets. Figure 4 shows the connection pattern of the interferometric image pairs in SH, where the yellow dots indicated the super-master image, while the green dots were the rest of the image. The connection pattern was uniform and dense, satisfying the required number of image pairs.

We then applied the PS-SBAS-InSAR technique to calculate the annual average deformation rates of 16 typical cities over the two years (January 2019–January 2021). Figure 5 shows that there was a correlation among the spatial patterns of deformation fields in geographically neighboring cities. During the monitoring period, the geographically contiguous zones in the northwest corner of the WX city and the southwest corner of the CZ city showed a more pronounced uplift, and a similar trend was observed in the HZ and SX cities. In addition, we categorized the annual mean surface deformation rate into 12 classes at 2 mm/year intervals and calculated the percentage of pixels in each grade (Figure 5b). Most of the pixels were distributed within [−4, 4] mm/year, suggesting that the YRD typically experienced slight deformation during the observational period (January 2019–January 2021).

Table 4. RLOS annual average deformation rate and its estimated residuals (January 2019–January 2021).

Category	City	Average Deformation Rate (mm/year)			Estimated Residuals (mm/year)
		Low	High	Mean	Mean	Standard Deviation
Core city	SH	−14.988	6.004	−1.817	2.103	1.922
Nanjing metropolitan area	NJ	−59.050	9.016	−1.345	2.486	1.139
	ZJ	−29.431	21.386	1.359	2.843	1.344
	YZ	−24.027	15.314	−0.034	2.635	1.380
Hangzhou metropolitan area	HZ	−40.213	12.898	−3.937	1.961	0.916
	JX	−20.627	9.832	−2.511	2.379	0.969
	HuZ	−52.280	18.304	−1.824	2.516	1.072
	SX	−28.226	9.599	−0.697	2.647	1.611
Hefei metropolitan area	HF	−28.017	17.880	−0.152	2.955	1.350
	WH	−51.116	67.897	−0.956	3.274	1.168
	Ma	−34.247	6.921	−2.972	3.343	1.571
Su-Xi-Chang metropolitan area	SZ	−16.831	9.391	−0.312	2.203	1.230
	WX	−29.147	11.381	−0.282	2.337	1.737
	CZ	−16.259	13.532	3.832	2.417	1.223
Ningbo metropolitan area	NB	−67.455	37.083	−3.247	3.514	1.909
	TZ	−73.640	34.860	−3.551	3.355	1.577

presents the average deformation rates with estimated residuals for the period January 2019–January 2021 in the 16 study areas. The standard deviation of the estimated residuals for all cities was <2 mm/year, demonstrating the validity of the Informer model. The mean of the annual average deformation rates of all pixels is a more important indicator of the deformation to characterize the entire region. The mean values showed an overall uplift in the ZJ and CZ cities during the observation period and a subsidence in other cities. The HZ, NB, and TZ cities exhibited typical significant subsidence trends with average deformation rates of <−3 mm/year.

The respective spatial visualization (Figure 6) reveals that most of the study areas displayed uniformity and minor local changes, such as in the YZ and HF cities. In contrast, the average deformation rate in the HZ city had a more distinct spatial pattern, gradually transitioning from an uplift in the northwestern region to subsidence in the southeastern region. Similarly, the JX city displayed a gradual shift from subsidence to uplift from west to east. Furthermore, there were clear uplift areas to the northwest of the WX city and southeast of the CZ city.

We classified local areas of each city into three categories based on average annual deformation rates: subsidence (<−2 mm/year), uplift (>2 mm/year), and terrain stability zones (others). We then selected six class-level landscape metrics to delineate the spatial patterns of the surface subsidence and uplift classes in terms of area/edge, shape, and aggregation.

Table 5. Landscape metrics (class-level) for RLOS annual average deformation rate.

City	Subsidence						Uplift
	Area-Edge		Shape	Aggregation			Area-Edge		Shape	Aggregation
	LPI (%)	ED (m/ha)	PAFRAC	LSI	DIVISION	AI (%)	LPI (%)	ED (m/ha)	PAFRAC	LSI	DIVISION	AI (%)
SH	9.137	25.745	1.263	66.69	0.990	93.741	0.026	1.519	1.255	30.397	1.000	78.069
NJ	0.765	11.281	1.261	61.849	1.000	88.521	0.346	4.196	1.267	40.356	1.000	86.967
ZJ	0.046	3.654	1.196	27.98	1.000	77.849	1.801	18.892	1.265	42.553	0.999	90.016
YZ	0.098	9.165	1.223	57.061	1.000	79.217	0.414	9.563	1.26	51.796	1.000	83.596
HZ	6.633	28.834	1.274	64.266	0.993	91.756	0.448	9.344	1.268	43.079	1.000	88.641
JX	5.784	22.396	1.273	34.192	0.996	90.848	0.386	5.856	1.227	20.558	1.000	87.533
HuZ	0.521	19.259	1.256	47.76	1.000	85.650	0.124	3.241	1.202	24.258	1.000	78.181
SX	2.331	14.676	1.265	36.628	0.999	88.549	1.136	11.861	1.259	33.508	1.000	88.190
HF	0.259	15.051	1.277	57.557	1.000	81.654	0.261	12.367	1.271	50.56	1.000	82.807
WH	1.062	11.967	1.247	67.633	1.000	88.582	0.156	9.659	1.228	71.2	1.000	84.280
Ma	21.353	51.903	1.295	37.896	0.954	90.950	0.093	1.224	1.199	9.266	1.000	78.197
SZ	0.772	11.434	1.245	44.454	1.000	85.050	0.163	7.334	1.261	37.496	1.000	83.475
WX	0.316	15.628	1.25	57.279	1.000	86.054	2.505	10.121	1.252	39.172	0.999	90.016
CZ	0.068	1.518	1.195	24.082	1.000	79.840	11.851	26.357	1.258	60.167	0.986	92.637
NB	6.103	16.099	1.279	92.438	0.996	92.047	0.219	8.330	1.24	86.289	1.000	86.608
TZ	3.757	33.883	1.3	86.533	0.998	89.427	0.083	13.720	1.262	78.294	1.000	78.599

shows that the subsidence class in the SH city had an LPI of 9.1%, a PAFRAC of 1.3, and an AI of 93.7%, indicating the existence of subsidence areas with high aggregation. The LPIs in the NJ, YZ, HuZ, and SZ cities did not exceed 1%, while the DIVISIONs all reached 1. These results suggest that the subsidence and uplift were relatively steady in these four cities, and the deformation was sporadically distributed across the region. The ZJ, WX, and CZ cities had more dispersed subsidence areas and more concentrated uplift areas. The CZ city had a wide range of uplift areas, with an LPI of 11.9% in the uplift class. In addition, the HZ, JX, and NB cities showed more concentrated subsidence (LPI > 5.0) and more dispersed uplift (LPI < 0.5). These metrics indicate that similar deformation patterns may exist among geographically neighboring cities.

4.2. Deformation Time-Series Analysis

Time-series analysis was required to further explore the spatial and temporal distribution characteristics of surface deformation in the YRD. Using Jan 2019 as the starting point, we calculated the displacements against Jan 2019 for the remaining 24 months. Negative and positive displacements denote subsidence and uplift, respectively. We selected SH (a coastal city), HF (an inland city), and CZ (a typical uplift city) as examples to show the time-series surface deformation every three months (Figure 7, Figure 8 and Figure 9). The maps show the changes in the spatial distribution of deformation, while each curve represents an uplift or subsidence process at a different sampling point against January 2019.

Figure 7 shows the spatial distribution of surface deformation in the SH city from January 2019 to January 2021. Overall, the SH city was observed to have minor deformation in local areas, including subsidence in the peripheral regions and uplift in the central region. Subsidence in the SH city showed an initial rise until July 2021, reaching a peak of about 10 mm, followed by a gradual decline until the end of the observation period. In contrast, there were more fluctuations in the uplift in the SH city, with two typical inflection points in January 2021 and July 2021, but it still showed an increase at the end of the observation period.

Figure 8 illustrates the changes in surface deformation over time in the HF city, with a typical subsidence funnel in the southeastern part. The subsidence in this funnel increased over time and stabilized by September 2021. The uplift in the HF city was characterized by alternating subsidence and uplift, as in the southern part. The uplift curves also show that two typical inflection points occurred in January 2021 and July 2021. These suggest that the deformation pattern was not temporally constant and may have shifted between subsidence and uplift.

The surface deformation in the CZ city from January 2019 to January 2021 showed high regularity (Figure 9), with a general upward trend that continued to rise slowly, peaking in September 2020 (>15 mm). Over time, the distribution of the uplifted areas shifted from dispersed to clustered, with a trend of a gradual southward shift in location. At the monitoring start time, the subsidence area was distributed in the southernmost part of central CZ city and then gradually moved northward. The subsidence first increased then decreased, reaching about −5 mm at the end of the observation.

4.3. Future Deformation Prediction

As the Informer model requires long time-series data, we expanded the data from 25 samples (with an interval of 1 month) to 2194 (with an interval of 8 h) using cubic spline interpolation while randomly selecting 10% of the samples to add Gaussian noise (error level 0.1). We used 80% of the new samples for training and the remaining 20% for testing. The Informer model was assigned an input sequence length of 120, a label length of 60, and a predicted future sequence length of 180. By default, the training epoch was set to 6, with a batch size of 32 and a learning rate of 0.0001. We then used the Informer models to predict future surface deformations for the next 2 months (180 time steps).

The surface deformation of the 16 selected typical cities in the YRD varied temporally but was expected to stabilize in the future (Figure 10). Some cities, such as YZ and WH, even showed a tendency for subsidence accumulation approaching zero. Geographically neighboring cities tended to have similar deformation trends. For example, the subsidence in the SZ, WX, and CZ cities typically showed an increasing pattern and then decreased over time, with a final subsidence of about −5 mm. These phenomena demonstrate the regional correlation of surface subsidence throughout the YRD. In addition, a few cities, including HZ, JX, and SZ, displayed more evident inflection points and thus required further attention.

Figure 11 shows the changes in uplift for typical cities in the YRD over the next 2 months, as predicted by the Informer model. Throughout the observational period, there were increased levels of uplift in the SH and NB cities, which were more volatile, with multiple extremes and inflection points. In contrast to the subsidence trends, most cities will still experience uplift in the future, which may reach about 5 mm by March 2021. Furthermore, the deformation patterns of the SH, NB, and TZ cities displayed great fluctuations compared to other cities, likely due to their proximity to the East China Sea, which is subject to more tidal action and sedimentation [3,40].

Table 6. Assessment of Informer model for RLOS time-series deformation (January 2019–January 2021).

City	Subsidence			Uplift
	MSE (mm/year)	RMSE (mm/year)	MAE (mm/year)	PCCs	MSE (mm/year)	RMSE (mm/year)	MAE (mm/year)	PCCs
SH	5.298	2.302	1.547	0.963	1.167	1.080	0.899	0.980
NJ	1.265	1.125	1.048	0.953	0.796	0.892	0.745	0.986
ZJ	1.034	1.017	0.929	0.998	1.123	1.060	1.000	0.933
YZ	0.890	0.944	0.647	0.978	0.185	0.430	0.374	0.992
HZ	3.127	1.768	1.610	0.870	0.089	0.299	0.241	0.946
JX	5.989	2.447	2.226	0.782	0.358	0.598	0.511	0.906
HuZ	4.081	2.020	1.968	0.969	2.298	1.516	1.377	0.964
SX	1.510	1.229	1.079	0.298	0.524	0.724	0.541	0.987
HF	2.599	1.612	1.383	0.994	0.148	0.384	0.294	0.956
WH	1.769	1.330	1.181	0.754	3.225	1.796	1.666	0.690
Ma	1.233	1.111	0.924	0.978	0.533	0.730	0.566	0.928
SZ	0.384	0.620	0.504	0.954	0.608	0.780	0.697	0.819
WX	0.142	0.377	0.305	0.955	1.134	1.065	0.998	0.596
CZ	0.528	0.727	0.570	0.587	1.957	1.399	1.217	0.783
NB	1.846	1.359	1.144	0.943	4.485	2.118	1.512	0.651
TZ	1.000	1.000	0.837	0.941	0.945	0.972	0.869	0.996

shows that the RMSE of the time-series fit in most cities did not exceed 2 mm/year and that all the PCCs were greater than 0.9. This indicates that the Informer model can accurately capture the historical subsidence changes.

5. Discussion

Large-scale short-term monitoring and prediction of ground deformation can effectively prevent and mitigate geohazards, ensuring sustainable urban development. Rapid economic development has led to the integration of towns and cities; therefore, geohazards are not usually based on administrative districts. Therefore, our study focused on the correlation of deformations among regions in the rapidly urbanizing and developing YRD. We monitored the surface deformation of 16 typical cities in the YRD from January 2019 to January 2021 using PS-SBAS-InSAR technology and successfully captured the spatial pattern and temporal evolution of deformation in each city. We used time-series analysis and landscape metrics to further delineate the spatiotemporal distribution and evolution of deformation in the YRD. Then, we reprocessed the time-series InSAR results and applied the Informer models to derive the uplift/subsidence intensity and trend for each city for the next 2 months (February 2021–March 2021).

The PS-SBAS-InSAR technique can provide complete and precise surface deformation sequences, and the standard deviations of the estimated residuals are basically < 2 mm/year for the representative cities of the YRD. Our findings are consistent with the conclusions in the literature [3,12] on urban surface deformation in the YRD region. In 2017, Han et al. [3] reported a maximum subsidence rate of −30.2 mm/year for the typical area of Shaoxing in 2013–2015, while the maximum subsidence rate in our study was −28.2 mm/year, both of which were very close to each other. These results prove that our approach is correct and effective for short-term monitoring and forecasting on a large scale [41].

All 16 representative cities experienced minor subsidence and uplift from January 2019 to January 2021, with mean deformation rates in the range of [−4, 4] mm/year. The deformation intensities within the monitoring period varied across space, but geographically neighboring cities had relatively similar deformation distributions (e.g., the WX and CZ cities). The distribution of the average deformation rate (Figure 5 and Figure 6) suggests that if one type of deformation (subsidence/uplift) is spatially concentrated, another type of deformation (uplift/subsidence) may be spread out in its vicinity. Landscape metrics subsequently quantified this phenomenon, which led to the recommendation that if local deformation is identified, it is necessary to broaden the monitoring area and make a unified decision about a large-scale area.

The trends of time-series deformation in 16 cities could not be fully expressed by the same specific mathematical expression [32]. However, cities within the same metropolitan area had highly similar temporal variations, such as the ZJ and YZ cities in the Nanjing metropolitan area and the HZ and JX cities in the HZ metropolitan area. Time series of deformation in many cities had one or several distinct inflection/extreme points, suggesting that the pattern of deformation is not static and may shift between subsidence and uplift. Interestingly, most of these points occurred near March, in spring, and July, in summer. In the Yangtze River Delta region, the early spring is cold and dry, and the summer is hot and rainy with frequent typhoons, so the change in soil water content is obviously seasonal. Studies have shown that changes in soil moisture content would lead to decreases in the cohesion and friction angle of compacted soil, thereby leading to ground deformation [42,43]. This further demonstrates the necessity for surface deformation monitoring on a large scale, especially in areas sharing the same hydrogeological unit. Therefore, we propose actively taking advantage of policies to achieve coordination and cooperation in the development, utilization, and management of groundwater resources in the YRD and jointly managing the regional surface deformation.

The Informer model provides a powerful tool for predicting future deformation patterns [23]. It is capable of capturing the variability in time-series deformation in the YRD. The predicted results of subsidence and uplift indicate that the 16 selected cities will have relatively stable deformation trends in the short term and will not show large fluctuations. This cumulative time-series slowdown was also evidenced in the early literature [44]. Nevertheless, consistent and precise monitoring of these regions is still necessary to assess potential risks and formulate appropriate measures for hazard prevention.

Additionally, the Informer model used in this study offers greater flexibility and adaptability than conventional time-series prediction models [23]. Several time-series models for InSAR deformation monitoring, such as recurrent neural networks, are thought to be sensitive to fixed input sequence lengths and require manual adjustment and preprocessing [45]. The Informer model features a highly configurable architecture that can process time-series data of varying lengths, frequencies, and characteristics. We compared the Informer model and the Long Short-Term Memory (LSTM) network in predicting the sinking of the city of Shanghai over the next two months (Figure 12). The Informer model had a smaller (by 1.5) RMSE than the LSTM model, indicating a better performance of the former in capturing the variations in the time-series data. The two methods are both neural network models, and they need to update their network layers in real time to capture temporal features when making predictions, whereas to capture the variation patterns of small samples requires identifying global signatures. LSTM is a model that captures local features and is susceptible to gradients, resulting in a possible failure to reflect the hidden patterns of small samples [45]. The Informer model introduces a global focus mechanism that allows the model to establish periodic dependencies between different time steps of the input time series. This helps the Informer model to be suitable for short-term predictions of surface deformation with smaller samples.

In the future, we will consider acquiring local groundwater data and engineering construction data to quantitatively analyze their correlations with deformation, optimize prediction models, and offer recommendations for managing surface deformation on a large scale.

6. Conclusions

In this study, we acquired 25 Sentinel-1 images from January 2019 to January 2021. We used the PS-SBAS InSAR technique to capture the surface deformation of 16 representative cities in the YRD and then applied the Informer model to make short-term predictions of future deformation trends. The results show that from January 2019 to January 2021, the YRD experienced slight deformation, with an average deformation rate within [−4, 4] mm/year. The Informer model successfully captured the time-series change in deformation, and its predicted deformation trend will slow down in the next 2 months. Nevertheless, we still recommend continuous monitoring of surface deformations in the YRD to support sustainable urban development.

We found that surface deformation in the YRD is regionally correlated. At a small scale (within an urban area), a concentrated distribution of deformation (subsidence/uplift) may be accompanied by some scattered distribution of another kind of deformation (uplift/subsidence). At a large scale (within an urban agglomeration), spatially neighboring cities may have a correlated deformation distribution and similar deformation time-series trends. Therefore, we suggest that there is a critical need to establish a comprehensive, large-scale monitoring network for ground deformation dynamics within urban conglomerates (e.g., the YRD). In addition, it is essential to take advantage of policies to jointly prevent and manage geological hazards.

Our methods contribute to resolving the problems of insufficient monitoring of large-scale deformation in coastal areas and attaining short-term high-precision time-series predictions with small samples. Future work will expand the scope and duration of monitoring while identifying the drivers of surface deformation, such as infrastructure development, groundwater extraction, and land reclamation. The drivers will then be added to the predictive model to optimize model performance and inform the prevention and management of large-scale surface deformation.