Reconstruction of spatially continuous time-series land subsidence based on PS-InSAR and improved MLS-SVR in Beijing Plain area

ABSTRACT

Beijing has undergone severe settlement in recent years. Persistent Scatterers Interferometric Synthetic Aperture Radar (PS-InSAR) technique has been widely used to derive time-series land deformation. However, existing studies have faced two challenges: (1) the nonlinear characteristics of time-series subsidence has not been fully investigated; (2) since PS points are normally distributed in urban areas with high building density, measurement gaps usually exist in nonurban areas. To address the challenges, we presented a new method to reconstruct spatially continuous time-series deformation. First, PS-InSAR was used to retrieve the deformation based on 135 scenes of Envisat ASAR and Radarsat-2 images from 2003 to 2020. Polynomial Curve Fitting (PCF) was then used to model nonlinear time-series deformation for the PS points. In the PS measurement gaps, Iterative Self-Organizing Data Analysis Technique (ISODATA) and Multi-output Least Squares Support Vector Regression (MLS-SVR) were used to estimate the PCF coefficients and then time-series deformation considering 40 features including thickness of the compressible layers, annual groundwater level, etc. The major results showed that (1) compared to linear, quadratic, and quartic models, cubic polynomial model generated better fit for the time-series deformation (R² ≈0.99), suggesting obvious nonlinear temporal pattern of deformation; (2) the time-series deformation over measurement gaps reconstructed by ISODATA and MLS-SVR had satisfactory accuracy (R² = 0.92, MAPE < 15%) and yielded higher accuracy (R² = 0.947) than IDW (R² = 0.687) and Ordinary Kriging (R² = 0.688) interpolation methods. The reconstructed results maintain the nonlinear characteristics and ensure the high spatial resolution (120 m) of time-series deformation. Among the 40 predictor variables, ground water level datasets are the most influential predictors of time-series deformation.

KEYWORDS:

• Land subsidence
• PS-InSAR
• multi-output support vector regression
• nonlinear deformation
• time series
• spatial interpolation

1. Introduction

Land subsidence, a destructive geologic disaster caused by human activities or natural factors, can lead to many problems, including building foundation settlement, underground pipeline destruction, and aggravation of flood disaster (Garg et al. 2022; Herrera-Garcia et al. 2021; Karanam et al. 2021). Beijing is one of the most important cities for settlement prevention and control in China (Gong et al. 2018). In the 1950s, extensive groundwater extraction for industrial development and the subsidence centers gradually formed in eastern Beijing (Gao et al. 2019). By 2018, 12.05% of the Beijing Plain had an average annual subsidence of more than 30 mm (Shi et al. 2020).

Interferometric Synthetic Aperture Radar (InSAR) technique has the capability to provide ground deformation across a vast geographical area, as opposed to traditional ways such as levelling and the Global Positioning System (GPS) (Ferretti, Prati, and Rocca 2000). Particularly, Persistent Scatterers Interferometric Synthetic Aperture Radar (PS-InSAR) technique, which derives deformation over point-like stable radar targets (called persistent scatterers), can reach millimeter to centimeter accuracy and has been successfully implemented in the research of the spatiotemporal evolution characteristics of land subsidence (Ng et al. 2012; Zhu et al. 2015; Chen et al. 2017b). Most of these studies investigated annual deformation rate or cumulative deformation and focused on the linear characteristics of time series subsidence. The nonlinear characteristics of long time series subsidence have been less analyzed in the existing studies (Chen et al. 2021; Park and Hong 2021).

On the other hand, the persistent scatterers detected by PS-InSAR techniques are those with high intensity and stable radar backscattering. PS-InSAR can hardly acquire deformation measurements over areas with low backscattering energy and spatiotemporal incoherence, such as vegetation with great seasonal change, bare soil, and flat roads, resulting in measurement gaps in these areas (Gao et al. 2019). Garg et al. (2022) and Karanam et al. (2021), which highlights the damaging effects of subsidence on infrastructure, including roads and aggravation of flood disaster in areas with sparse PS points. In order to increase measurement density, distributed scatterer (DS) interferometry such as SBAS (Berardino et al. 2002; Liu et al. 2021) and SqueeSAR (Ferretti et al. 2011) techniques have been proposed. However, these methods cannot derive spatially continuous measurements (Hooper et al. 2008; Bui et al. 2021; Cao, Lee, and Jung 2016). Therefore, the practical applicability of these technologies for monitoring land deformation over a vast region are currently restricted. Spatial interpolation methods have been frequently employed in recent decades to predict deformation and rebuild spatially continuous land subsidence measurement over a large area without PS measurements. Chen et al. (2017a) and Duan et al. (2020) employed Ordinary Kriging and Inverse Distance Weighted (IDW) to interpolate PS points, respectively, to obtain settlement along Beijing Subway Line 6. Ikuemonisan et al. (2020) used a geostatistical approach to investigate the variability of land subsidence in Lagos, Nigeria. These interpolation approaches are predicted on the idea that known and unknown data have similar statistical or geometrical structures and perform best when the measurements are evenly distributed (Shen et al. 2015). However, PS points were not evenly distributed in most cases, and most interpolation techniques do not take into account the geohydrology and lithology background, which are important influences on land subsidence.

Machine learning has been applied to analyze the risk of land subsidence and the contribution of various deriving factors. Rahmati et al. (2019) used two machine learning algorithms, Maximum Entropy and Genetic Algorithm Rule-set Production, to study primary mechanisms contributing to the risk in Kashmar, Iran. Zhou et al. (2019) and Chen et al. (2020) used the Gradient Lifting Decision Tree and Random Forest algorithms to quantify the contributions of factors to land subsidence in Beijing, respectively. The potential of machine learning approach in reconstruction of spatially continuous land subsidence has not been investigated. In addition, traditional machine learning applied in remote sensing generally considers single-output problem. Multi-output machine learning models the relationships between the outputs and the relevant features, ensuring a more accurate depiction of the real-world situations (Xu et al. 2019).

The objectives of this study were to (1) model nonlinear deformation time series based on multi-platform SAR datasets, and (2) reconstruct spatially continuous time-series deformation by filling PS-InSAR measurement gaps in Beijing Plain area using integrated Iterative Self-Organizing Data Analysis Technique (ISODATA) and Multi-output Least Squares Support Vector Regression (MLS-SVR) algorithm. We expect that the proposed method will be capable of accurately reconstructing subsidence measurements with high spatial and temporal resolutions.

2. Study site and data set

2.1. Study site

Beijing Plain (115°43′E-117°19′E, 39°26′N-40°29′N) is in the southeast of Beijing city (Figure 1). Beijing Plain is formed by Juma River, Yongding River, Chaobai River, Wenyu River, and Beiyun and Jiyun River carrying a large amount of sediment and gravel deposited on the basement (Gao et al. 2019). From the top to the lower part of alluvial fan, the grain size of the aquifer changes from coarse to fine. Climatically, Beijing Plain has a typical semi-humid and semi-arid continental climate, high temperature and concentrated precipitation in summer months, with annual average temperature of 11–13°C and rainfall of 585 mm (http://swj.beijing.gov.cn/zwgk/szygb/).

Figure 1. The location of the study area. (a) Location of Beijing; (b) Geography of Beijing Plain; (c) Location of SAR datasets.

The Beijing Plain has observed the development of several active high-angle faults in (Figure 1). Nankou-Sunhe fault has the properties of stick-slip and creep-slip and is the only Northwest-Southeast trending fault. Liangxiang-Shunyi, Nanyuan-Tongzhou, and Huangzhuang-Gaoliying faults are Southwest-Northeast trending faults (Chen et al. 2016). Relevant studies have demonstrated that the location of faults determines the distribution of subsidence zones. And the linear boundaries of the land subsidence correspond to the several active fault traces in Beijing plain area (Hu et al. 2019; Lyu et al. 2020; Guo et al. 2020).

In 1969, the first subway line in China was officially opened in Beijing. Beijing had 23 subway lines totaling 699.3 kilometers and 405 stations by the end of 2019 (Figure 1). The excavation of foundation pit and underground engineering during subway construction will lead to soil movement and land subsidence (Chen et al. 2017a). The cyclic dynamic load generated during subway operation also has a certain influence on land subsidence (Wang et al. 2020).

2.2. Data set

2.2.1. SAR data

A total of 135 scenes of Envisat ASAR and Radarsat-2 were used in this study to obtain the land subsidence monitoring results. Envisat is a solar-synchronous polar orbiting satellite launched by the European Space Agency (ESA) on 1 March 2002. The Envisat ASAR sensor operates in the C-band with a wavelength of 5.6 cm (https://earth.esa.int/eogateway/instruments/asar). Radarsat-2 satellite, a collaborative project between the Canadian Space Agency (CSA) and MacDonald Dettwiler Associates (MDA), was launched on 14 December 2007 and operates in the C-band (https://earth.esa.int/eogateway/missions/radarsat). Fifty-five Envisat ASAR images in Image mode with resolution of 30 m (18/06/2003–19/09/2010) were captured. Thirty-seven Radarsat-2 images in Standard mode with medium resolution of 30 m (22/11/2010–21/10/2016) were captured. Forty-three Radarsat-2 images in Extra-Fine mode with high resolution of 5 m (25/01/2017–10/01/2020) were captured. Main information of the Envisat ASAR and Radarsat-2 images are shown in Table 1.

Table 1. Main information of the Envisat ASAR and Radarsat-2.

SAR Sensor	Envisat ASAR	Radarsat-2	Radarsat-2
Orbit direction	Descending	Descending	Descending
Polarization	VV	VV	VV
Band	C	C	C
Wavelength (cm)	5.6	5.6	5.6
Incidence Angle	22.8	33.9	22.56
Image mode	Image	Standard	Extra-Fine
No. images	55	37	43
Date range	18/06/2003–19/09/2010	22/11/2010–21/10/2016	25/01/2017–10/01/2020

2.2.2. Ancillary datasets

Ancillary datasets include the quaternary thicknesses of compressible deposits within each of the three compressible layers (depth of bottom <100 m, 100 m ~300 m and >300 m, respectively), total thickness of the compressible layers (Figure S1), annual groundwater level from 2003 to 2019 (Figure S2), annual Normalized Difference Building Index (NDBI) from 2003 to 2019 (Zha, Gao, and Ni 2003) (Figure S3), and vector layers of faults and subways (Figure S4). NDBI represents the density of buildup areas and was derived from Landsat images using Google Earth Engine platform (Dai, Guldmann, and Hu 2018). The vector lines of faults and subways were converted to density of faults and subways, respectively, representing the spatial distribution of geological structures and dynamic loads, respectively (Shi et al. 2020). All 40 datasets were then converted to raster layers and utilized as feature variables.

3. Methodology

Figure 2 outlines the methodology. First, the PS-InSAR was applied to calculate the land deformation measurements on PS points using SAR datasets (Section 3.1). Second, Polynomial Curve Fitting (PCF) method was applied for each PS points to build nonlinear deformation time-series model, and four parameters of the polynomial curve were obtained (Section 3.2). Third, Fishnet partitioning was used to identify the grids with measurement gaps where time-series deformation needs to be reconstructed (Section 3.3). Finally, an improved MLS-SVR method, which incorporates ISODATA, was utilized to reconstruct time-series deformation within the grids with measurement gaps, and the results were validated (Section 3.4).

Figure 2. Flowchart of methodology.

3.1. Generation of time series settlement using PS-InSAR

PS-InSAR technique in GAMMA software was utilized to obtain surface deformation based on Envisat ASAR images (2003–2010), Radarsat-2 images (2010–2016), and Radarsat-2 images (2017–2020). The selection of master image was based on a function of temporal and perpendicular baselines and Doppler centroid frequency, which roughly minimized the spatial and temporal de-correlation effects (Foroughnia et al. 2019; Kampes 2005). Envisat ASAR image acquired on 23 January 2008, Radarsat-2 image (Standard mode) acquired on 25 November 2014, and Radarsat-2 (Extra-Fine mode) image acquired on 24 August were selected as master images for the three time periods, respectively. Then, employing Shuttle Radar Topography Mission Digital Elevation Model (SRTM DEM) and precise orbital data, respectively, a series of differential interferograms were generated and flat earth effect were removed. Afterward, persistent scatterer candidates (PSCs) were selected by using an Amplitude Dispersion Index (ADI) threshold. According to the experience and spatial distribution of PSC points, the ADI thresholds of the three datasets were 0.6, 0.6, and 0.55, respectively. And the interferometric phase of flattened differential interferogram were divided into different components (Yang et al. 2018). The vertical deformation was then generated by converting the deformation in the Line of Sight (LOS) direction (Gao et al. 2018).

Next, according to the existing research and experience, the Nearest Neighbour Search method was adopted to merge time-series deformation in each of the three periods into a single time series (Lyu et al. 2020; Shi et al. 2020). We obtained the nearest PS point generated from Radarsat-2 (Standard mode) and Radarsat-2 (Extra-Fine mode) datasets for each PS point produced from Envisat ASAR. As there are several-month time span between two consecutive time periods (e.g. 19 September 2010 to 22 November 2010, and 21 October 2016 to 25/01/2017), we first calculated the deformation between the end of the previous period and the start of the following period, and then generated the time-series deformation of the next period (Eq. 1(1) $D_{end1-start2}=\frac{V_1+V_2}{2}\ast \frac{Da{y}_{end1-start2}}{365}$(1) , 2(2) $D_2^{{ }^{\prime }}=D_2+D_{end1-start2}$(2) ).

(1) $D_{end1-start2}=\frac{V_1+V_2}{2}\ast \frac{Da{y}_{end1-start2}}{365}$(1)

(2) $D_2^{{ }^{\prime }}=D_2+D_{end1-start2}$(2)

where $V_1$ and $V_2$ denote the annual deformation rate during the first period (e.g. 18/06/2003–19/09/2010) and the second period (e.g. 22/11/2010–21/10/2016), respectively. $Da{y}_{end1-start2}$ denotes the days from the end of the first period to the start of the second period. denotes the deformation from the end of the first period (e.g. 19/09/2010) to the start of the second period (e.g. 22/11/2010). ${\mathit{D}}_2$ denotes the time-series deformation ($D_{21},D_{22},\dots ,D_{2i},\dots ,D_{2n}$) in the second period derived from PS-InSAR, and ${\mathit{D}}_2^{\prime }$ denotes the estimated time-series deformation. This fusion method assumes that the land subsidence of several-month time span presents linear characteristics, and the effect on the nonlinear characteristics of long time series can deformation be negligible (Lyu et al. 2020; Shi et al. 2020).

3.2. Polynomial curve fitting for time-series deformation

Our previous research revealed the nonlinear trend of deformation development in Beijing due to the overexploitation of groundwater and the implementation of South-to-North Water Diversion Project (Lyu et al. 2020, 2020). PCF was adopted to model the time-series cumulative deformation for each PS point. Every PS point has a set of paired measurements ($t_1,D_1$),$\hspace{0.17em}\left(t_2,D_2\right)$, … , ($t_{135},D_{135}$). Here, $t_i$ was calculated as $t_i=\frac{Da{y}_i}{365}$, where $Da{y}_i$ represents the days between the date of ith image and 1 January 2003, and $D_i$ represents the cumulative deformation of the ith image acquisition date. The nth order polynomial equation is represented as:

(3) $D\left(t\right)=a_0{t}^n+a_1{t}^{n-1}+\dots +a_n$(3)

We computed the coefficients of the polynomial fitting function using the least square method. In this study, simple linear model ($n=1$), quadratic ($n=2$), cubic ($n=3$), and quartic ($n=4$) polynomial models were used to fit time-series deformation, and their performances were compared to determine the best suited model.

3.3. Fishnet partitioning and locating grids to be reconstructed

The PS points detected by PS-InSAR technique are mainly found in man-made things like building facades and corners. In areas covered by bare soil, grass, forest, or crops, few PS points can be detected, resulting in measurement gaps. In this step, we aimed to identify the areas with measurement gaps where time-series subsidence need reconstruction in the next step. First, we partitioned the study area into a fishnet of square grids with size of 120 m × 120 m. The number of PS points within each grid was counted and assigned to the grid (Figure 3). The grids with no PS points were considered as NoData grid. For these grids, time-series deformation need to be reconstructed. We named these grids as grids to be reconstructed (Figure 2).

Figure 3. Flowchart of identification of NoData grids. The yellow grids have no PS points and time series deformation need to be reconstructed.

3.4. Reconstruction of time series deformation using ISODATA and MLS-SVR

In this step, we determined the reconstruction method for the NoData grids. Previous studies have shown that the influencing factors such as thickness of compressible layers, groundwater depth, and static and dynamic loads may affect the time-series development of land subsidence. Here, we developed a machine learning method by integrating ISODATA and MLS-SVR to estimate the time-series deformation.

Many studies have proved that the prediction within relatively homogeneous areas can improve the model performance (Cao et al., 2003; Duan et al. 2011; Lu and Chang 2014). ISODATA is an unsupervised classification method. It employs the maximum-likelihood decision rule to determine the class mean that is uniformly distributed throughout the data space and then repeatedly clusters the remaining pixels using the minimum distance methodology. This process of machine learning is repeated until the maximum number of iterations is reached, or the change in the percentage of pixels in each class is less than the selected pixel change threshold. In this study, we utilized ISODATA to cluster the 40 feature layers into five classes. The 40 feature layers including quaternary thicknesses of compressible deposits within each of the three compressible layers (depth of bottom <100 m, 100 m ~300 m and >300 m), total thickness of the compressible layers, annual groundwater level from 2003 to 2019, Normalized Difference Building Index (NDBI) from 2003 to 2019, density distribution of faults, and density of subways.

For each class, we trained MLS-SVR model to learn the correlation between the polynomial coefficients of time-series deformation and the feature variables, and the model was then used to predict polynomial coefficients for the NoData grids. In this model, each point corresponds to a unique set of parameters, and each variable is a separate input variable. MLS-SVR is an extension of the single output Least Square SVR (LS-SVR) and can effectively solve the problem of traditional single output SVR that disregards the relationship between different outputs (Xu et al. 2013). Given a set of samples ${\left\{\left(x_i,y_i\right)\right\}}_{i=1}^l$ with $x_i\in R^d$ and$y_i\in R^m$, the aims of the multi-output SVR regression is to predict an output vector$y_i\in R^m$ from an input vector$x_i\in R^d$. In this study, the input predictor variable $x_i$ represents each of the 40 feature variables, and the outputs are the polynomial coefficients.

According to the principle of structural risk minimization, MLS-SVR needs to solve the following problems:

(4) $min\frac{1}{2}{\sum }_{j=1}^k||W_j|{|}^2+\frac{c}{2}{\sum }_{i=1}^l{\sum }_{j=1}^k{e}_{i,j}^2+c_0{\sum }_{i=1}^l{\eta }_i\text{breaks}.t.y_{i,j}=W_j^T\phi \left(X_i\right)+b_j+e_{i,j},{\eta }_i={\sum }_{j=1}^k\left|e_{i,j}\right|$(4)

where $c_0$ is the penalty coefficient of samples, $\eta$ is the error of overall sample curve fitting, $c$ is the penalty coefficient of the single output fitting error, $e$ is the output error of the sample, and $b$ is the deviation and $W$ is the weight.

Then Lagrangian function and Karush-Kuhn-Tucker (KKT) condition were used to solve the above formula, and finally, the multi-output results are predicted. Kernel function can convert lower-dimensional space into a higher-dimensional space. Here, the Radial basis function (RBF) was employed .

To establish the MLS-SVR model, 80% of the grids with PS points were chosen as training samples. For these grids, the feature variables were normalized, and the polynomial coefficients of the PS points within each grid were averaged and normalized. The remained 20% grids were used as test samples.

4. Results

4.1. PS-Insar results and validation

There were 931,543 PS points, 1,079,011 PS points, and 1,070,332 PS points generated from Envisat ASAR, Radarsat-2 (Standard mode), and Radarsat-2 (Extra-Fine mode) datasets, and the maximum deformation rate were −136.54 mm/year (2003–2010), −161.03 mm/year (2011–2016), and −134.50 mm/year (2017–2020), respectively (Figure 4a–c). Radarsat-2 (Standard mode) and Radarsat-2 (Extra-Fine mode) PS pixels that were within 10 m of Envisat ASAR PS pixels were judged to be in the same place. The fused deformation rate ranged from −138.53 mm/year to 6.20 mm/year during 2003–2020 (Figure 4d). The Land subsidence centers were found in Beijing Plain area’s north, northeast, and east. Among these subsidence centers, Dongbalizhuang-Dajiaoting (DD) and Jinzhan (JZ) subsidence centers have longer subsidence development history, and the settlement funnel center of Changping (CP) and Shunyi (SY) subsidence centers have gradually developed in recent years. These centers have gradually expanded and connected. The DD subsidence center had suffered more severe ground subsidence, with a greatest rate of 138.53 mm/year.

Figure 4. Deformation rate derived from (a) Envisat ASAR (2003–2010), (b) Radarsat-2 (2010–2016) and (c) Radarsat-2 (2017–2020); (d) time series fusion during 2003–2020.

For validation, the deformation rate and cumulative time-series deformation were compared to levelling data obtained from 2006–2013 and 2015–2016 (Figure 4). The root-mean-square-error (RMSE) of PS-InSAR deformation rates were within 8 mm/year, and the time-series cumulative deformation showed good agreement (Figure 5). From the trend of time-series leveling and InSAR measurements, it is evident that the land subsidence developed faster in 2010–2014 than that during 2003–2010, and then developed more slowly after 2014.

Figure 5. Comparison of deformation derived from PS-InSAR and levelling observations. (a) Scatter plot of deformation derived from PS-InSAR and levelling observations during 2006–2013 and 2015–2016. (b, c) Time series deformation during 2003–2020.

4.2. Fitted curves of time series deformation

According to the National Ground Settlement Prevention Plan (NGSP 2011–2020) promulgated by Chinese government in 2010, land subsidence should be controlled in areas with settlement rate larger than 15 mm/year. Time-series land subsidence is reconstructed for the subsidence area, which we designated as having a subsidence rate more than 15 mm/year (Figure 4d). The performances of simple linear model, quadratic, cubic, and quartic polynomial models were assessed and compared using R² and RMSE of the deformation time series. Figure 6 showed that the linear model resulted in the lowest accuracy, with the minimum of R² of 0.61, and RMSE ranging from 4.33 mm to 161.55 mm. R² of quadratic polynomials were also relatively low, ranging from 0.73 to 0.99, and RMSE ranging from 2.29 mm to 100.03 mm. Cubic and Quartic polynomials yielded a better fitting performance, with R² ranging from 0.86 to 0.99, and RMSE lower than 60.28 mm and 30.95 mm, respectively. On the whole, the R² of 75% PS point cubic polynomial and quartic polynomial fitting is greater than 0.99, and the RMSE of 75% PS point cubic polynomial and quartic polynomial fitting is less than 14.03 mm and 9.65 mm, respectively, showing a small difference in fitting accuracy. In this study, we selected the cubic polynomial model: $\mathrm{D}=a_0{t}^3+a_1{t}^2+a_{2}t+a_3\hspace{0.17em}$to fit the deformation time series for each PS points, and four correlation coefficients ($a_0\mathrm{t}\mathrm{o}{a}_3$) were characterized the time-series deformation curve. Figure 7 demonstrates the cubic polynomial-fitted curve for two random PS points (Point A and Point B), with R² values of 0.996 and 0.997, respectively.

Figure 6. (a) R² and (b) RMSE of fitted cumulative deformation derived from linear, quadratic, cubic and quartic fitting.

Figure 7. Time-series cumulative deformation modelled by cubic polynomial fitting curve for (a) Point A and (b) Point B.

The spatial distribution of the four coefficients ($a_0$ to $a_3$) of cubic polynomial curve has a high correlation with the spatial pattern of subsidence rate (Figure 8). The values of the coefficient of cubic term $a_0$ and quadratic term $a_1$ can better reflect whether the nonlinear characteristics are significant. It is obvious that the coefficients $a_0$ and $a_1$ are close to 0 around the edge of the subsidence area, while in the subsidence centers, the absolute values of these two coefficients are relatively large, indicating that the nonlinear characteristics of land subsidence are more obvious at the subsidence center than that around the edges. The linear term $a_2$ represents the linear trend of the time-series deformation, thus the spatial distribution of $a_2$ is basically consistent with the subsidence rate. The above results confirm that it is reasonable to characterize time-series deformation by polynomial curves as the nonlinear features of the settlement in this study area are obvious, especially in the subsidence centers.

Figure 8. Spatial distribution of four correlation coefficients (a) $a_0$, (b) $\hspace{0.17em}{a}_1$, (c) $a_2$ and (d) $a_3$.

4.3. Identification of grids to be reconstructed

This study fused the time-series settlement sequence based on the monitoring results of Envisat ASAR (30 m resolution). Therefore, when determining the proper size of Fishnet grid, 30 m was taken as the minimum grid size, and 2 times expansion was carried out to discuss the number and proportion of PS points in grids of different sizes. Table 2 shows that when the size of grid varies from 60 m and 120 m, the proportion of grids with PS points is similar, which is about 37.88%. Compared to 60 m grid, 120 m grid was used for space segmentation, which requires less grids processing and has higher computation efficiency. It is worth noting that when the size of the grid is 960 m, 97.05% of the grid has PS points, and it is not until the size of the grid is 3840 m that the grid can fully cover all PS points. This indicates that the spatial distance of the region without measurements in subsidence area is greater than 1920 m.

Table 2. Points at different scale grids.

Grid size (m)	Number of grids	Grids Containing PS Points	Proportion (%)
30	1,306,362	89,831	6.87
60	451,373	170,975	37.88
120	113,826	43,109	37.88
240	28,940	17,057	58.94
480	7,463	6,125	82.07
960	1,983	1848	93.19
1920	513	536	97.50
3840	145	145	100

Figure 9 demonstrates that the PS points are unevenly distributed over the study area, with high density on buildings and low density in area covered by vegetation, water, or roads. In particular, in the subsidence center regions (Region I-IV, DD, JZ, SY and CP subsidence centers), a lot of grids had no PS measurements, and the time-series deformation need to be estimated to fill the data gaps. The subsidence area was partitioned into a total of 113,826 grids with size of 120 m × 120 m using Fishnet analysis tool. A total of 43109 (37.87%) grids have PS points, and the maximum number of PS points per grid is 38. There are 70,717 grids without PS points, accounting for 62.13% of the subsidence area, and the time-series deformation in these grids were further reconstructed.

Figure 9. Spatial distribution of the grids with and without PS measurements in subsidence area.

4.4. Reconstructed time series deformation and validation

4.4.1. Evaluation of estimated polynomial parameters

Five clusters from the 40 input feature datasets were selected using ISODATA (Figure 10a). Type 1 is mostly located around the south of Chaoyang subsidence center, northwest and south of the subsidence area. The regions of the subsidence area identified as Type 2, Type 3, and Type 4 are located, respectively, in the northwest, east, and south. Using the normalized feature values as independent variables and normalized cubic polynomial parameters ($a_{n0}$, $a_{n1}$, $a_{n2}$, $a_{n3}$) as dependent variables, MLS-SVR model were trained and used to predict polynomial parameters for other grids. The model outputs ($a_{n0}$, $a_{n1}$, $a_{n2}$, $a_{n3}$) were then used to calculate ($a_0$, $a_1$, $a_2$, $a_3$) and generate time-series deformation with Eq (3)(3) $D\left(t\right)=a_0{t}^n+a_1{t}^{n-1}+\dots +a_n$(3) . R² and RMSE of the estimated $a_{n0}$, $a_{n1}$, $a_{n2}$, $a_{n3}$ and time-series surface deformation were calculated based on test samples for model evaluation (Figure 10b,c). In order to verify the necessity of ISODATA clustering, we also operated MLS-SVR modeling without ISODATA.

Figure 10. (a) Spatial distribution of five clusters clustered by ISODATA; (b) R² and (c) RMSE of the estimated $a_{n0}$, $a_{n1}$, $a_{n2}$, and $a_{n3}$.

Figure 10b shows that the average R²s of the predicted normalized polynomial coefficients for each cluster (ranging from 0.71 to 0.80) were significantly higher, and the RMSEs for each cluster were obviously lower than the case without ISODATA clustering, indicating that the prediction accuracy of MLS-SVR was improved when the region was partitioned with ISODATA method. Especially in Type 1, the R² can reach 0.87 and the RMSE was lower than 0.054. The RMSEs for Type 3, 5 were lower than 0.033 and 0.041. The reason can be explained that the area of type 3, 5 are mostly found in CP and SY subsidence areas with low subsidence rate, resulting in a narrow range of deformation and relatively small RMSE for settlement reconstruction. However, if the deformation is reconstructed without ISODATA partition, the heterogeneity of the characteristics of the whole region led to the low accuracy with low R² and high RMSE of model reconstruction. The prediction accuracies of MLS-SVR were considerably improved when the features were clustered with ISODATA method. Therefore, we called the MLS-SVR combined with ISODATA as “improved MLS-SVR” algorithm.

4.4.2. Evaluation of reconstructed time-series deformation

The deformation reconstructed by the improved MLS-SVR model was compared with that detected by PS-InSAR to further its accuracy (Figure 11). Figure 10a demonstrates the average RMSE and average mean absolute percentage error (MAPE) of the cumulative deformation estimated from the improved MLS-SVR for the test grids with different deformation rates. In the edge of subsidence area (−15~−40 mm/year), the average RMSE is around 22.2 mm. In the region with deformation rate greater than −120 mm/year, the average RMSE reaches 102.68 mm. The average MAPE of reconstructed deformation ranges from 11.5% to 13.9%. Figure 11b shows the scatterplots of the time-series deformation estimated from MLS-SVR model against those measured by PS-InSAR from 2003 to 2020. The MLS-SVR estimation shows high consistency with PS-InSAR measurements (R² >0.90). For the two example test samples, the overall trend of time-series deformation reconstructed by MLS-SVR is the same as that of PS-InSAR time series. The slight difference in the deformation trend after 2017 May be caused by the impact of the South-to-North Water Diversion Project on the groundwater level, which has a certain impact on the nonlinear temporal subsidence predicted by the model (Figure 11c,d). Two leveling measurements from 2013 to 2018 also proved that the reconstructed time-series deformation have high accuracy (Figure 11e,f).

Figure 11. (a) RMSE and MAPE of reconstructed deformation of test samples with different deformation rate. (b) Scatterplots of cumulative deformation the improved MLS-SVR against that derived from PS-InSAR from 2003–2020. (c) and (d) Cumulative deformation time series reconstructed by MLS-SVR and derived from PS-InSAR of two test samples. (e) and (f) Cumulative deformation time series reconstructed by MLS-SVR and measured from levelling measurements.

4.4.3. Comparison with IDW and Kriging method

Previous research generally adopted simple spatial interpolation methods, such as IDW and Kriging, to obtain surface deformation over measurement gaps (Chen et al., 2017; Duan et al. 2020). To further verify the superiority of MLS-SVR over spatial interpolation in areas with sparse observations, we compared the cumulative deformations reconstructed by improved MLS-SVR and those by IDW and Ordinary Kriging in region I-IV when all PS points were removed from each region (Figure 12). There were 5,411 PS points, 4,879 PS points, 2,104 PS points, and 3,026 PS points removed from region I-IV, respectively. And the blue grid of region I-IV in Figure 9 shows the spatial distribution of these removed PS points. It is found that the improved MLS-SVR model can well reflect the spatial variations of settlement on different ground objects compared to IDW and Ordinary Kriging. Taking the removed PS points as validation data, the time-series deformations from the improved MLS-SVR have higher agreement with PS-InSAR than those from both spatial interpolation methods (Figure 12 m-o), with R² of 0.947, compared to IDW (R² = 0.687) and Ordinary Kriging (R² = 0.688). The central points P1-P4 of the region I-IV were selected, and the time-series deformation was reconstructed by IDW; Ordinary Kriging and MLS-SVR at this location was compared with PS-InSAR measurements (Figure 13). It is obvious that the time-series deformation from the improved MLS-SVR is more consistent with the PS-InSAR measurements, while the results of IDW and Ordinary Kriging reconstruction gradually become invalid with the extension of time.

Figure 12. Comparison of cumulative deformation reconstructed by IDW Ordinary Kriging and the improved MLS-SVR after removing PS points in (a-c) Region I, (d-f) Region II, (g-i) Region III and (j-l) Region VI. Scatterplots of reconstructed cumulative deformation from (m) IDW, (n) Ordinary Kriging and (o) Improved MLS-SVR against PS-InSAR measurements.

Figure 13. Comparison the time-series deformation between measured from PS-InSAR and reconstructed from IDW, Ordinary Kriging, and Improved MLS-SVR on point (a-d) P1-P4.

4.4.4. Spatially continuous time-series deformation in Beijing

The trained MLS-SVR models were used to predict time-series deformation on the center points of NoData grids (Figure 14b), and the predicted deformation on these points merge with PS points (Figure 14a) to reconstruct deformation over the entire area (Figure 14c). The cumulative deformation ranges from −57 mm to −2269 mm from 2003 to 2020 (Figure 14c) and has good spatial continuity. The reconstructed deformation in small regions (Region I-IV, Figure 14d–o) have a high spatial resolution, and the distance between the points is no more than 120 m (Figure 14c–o). Especially, in the southwest built-up area of Region III, it can be found that the reconstructed settlement is greater than that of the surrounding non-built-up area.

Figure 14. Cumulative deformation of (a) PS points derived from PS-InSAR, (b) NoData grids derived from MLS-SVR, and (c) the subsidence area merged by (a) and (b); (d, g, j, m) are PS points in Region I-IV; (e, h, k, n) are NoData grids in in region I-IV; (f, i, l, o) are the merged results in Region I-IV.

Two small areas, A (Figure 14g–i) and B (Figure 14m–o), are selected to show the details of the reconstructed settlement. Area A and Area B are each 3 km². Sunhe Fault passes through Area A, and Changping subway passes through Area B. There are only a few PS points along the fault and the ubway, resulting in the missing of subsidence information (Figure 15a–d). After the time-series deformation reconstructed by MLS-SVR, the spatial resolution of deformation can reach 120 m (Figure 15b–e). Figure 15c–f shows the profile analysis results of the fault and subway. On both sides of the fault and subway, it is clear that the land subsidence changes significantly, indicating that the fault and subway have a segmental influence on the spatial pattern of settlement. It also proves that the deformation reconstructed by improved MLS-SVR can reflect spatial variability.

Figure 15. Cumulative deformation of (a, d) PS points derived from PS-InSAR, (b,e) reconstruction final result in Area a and Area B. (c, f) Profile analysis in regions a and B.

5. Discussion

5.1. Influence of feature variables on model performance

Control variable method was adopted to examine how feature variables affect model performances. We categorized the 40 feature variables into five types, i.e., (1) thickness of compressible layers (first, second, and third compressible layers), (2) groundwater level (including annual groundwater level from 2003 to 2019), (3) NDBI (including annual NDBI from 2003 to 2019), (4) density of fault, and (5) density of subways.

Figure 16 shows the change in R² and RMSE of deformation estimated by MLS-SVR if one type of variables was removed from the predictor variables or if only one type of variables was used. It can be clearly seen that model performances declined if any type of the variables were removed (R² decreased and RMSE increased), while the removal of groundwater level layers led to the greatest decrease in R² and increase in RMSE. If only one type of feature variables were considered, ground water level variables resulted in the smallest decrease in R² (0.0473) and smallest increase in RMSE (22.62 mm). This indicates that ground water level datasets are the most influential factors for the spatial prediction of the ime-series deformation. The thicknesses of compressible layers are slightly less important. In Beijing, the main factor causing settlement is excessive groundwater consumption (Qin et al. 2018; Zhang et al. 2014), and the distribution of groundwater funnel center is basically same as that of subsidence center in space (Li et al. 2017). Land subsidence magnitude can be influenced by the thickness of the compressible soil layer. In Beijing plain area, the compressible thickness and groundwater level contributed more than 60% to land subsidence (Zhou et al. 2019). Due to the internal and external causes of land subsidence, the model is more sensitive to groundwater level and compressible soil layer thickness.

Figure 16. The change of R² and RMSE of cumulative deformation after removing one feature type (NoC- No compressib layer; NoG-No Groundwater; NoN-No NDBI; NoF-No Fault; NoS-No subway;) or with only one feature type (OnlyC-Only compressib layer; OnlyG-Only Groundwater; OnlyN-Only NDBI; OnlyF-Only Fault; OnlyS-Only subway).

The time-series NDBI has smaller impact on the accuracy of the model. The average decrease in R² is 0.0067 and 0.6193, the mean RMSE increase are 4.64 mm and 202.98 mm under removing and only retaining NDBI, respectively. In the past few decades, the urban area of Beijing has been continuously expanding and the building area has been rapidly increasing. The rapid construction of buildings in metropolitan area may result in settlement (Li et al. 2020; Yang et al. 2018).

Removing the fault and subway feature vectors has the least impact on the accuracy of the model, and the average decreases in R² are 0.0052 and 0.0054, the increases in RMSEs are 3.37 mm and 3.41 mm, respectively. Meanwhile, only retaining fault or subway feature vectors shows poorest model prediction accuracy, and the R² decreases are 0.5845 and 0.5694, the RMSE increases reach 195.42 mm and 192.08 mm, respectively. Chen et al. (2016) found that geological faults have a segmented influence on the distribution of land subsidence. Hu et al. (2019) analyzed the subsidence boundaries and fault traces and showed that the subsidence patterns are spatially controlled by geological faults. The geographical distribution of ground subsidence on both sides of the fault and the subway is significantly different, as shown in Figure 15. During the construction period, the subway disturbed the stratum and increased the inhomogeneity of settlement (Chen et al., 2017). During the operation period, the heterogeneity of settlement rate also increased in the area along the subway. However, the affected area of subway and fault on land subsidence is limited. The impact on land subsidence is significant in regions where faults and subways are densely dispersed, whereas the influence on other areas is slight.

5.2. The advantages of the improved MLS-SVR

In section 4, it is obvious that the improved MLS-SVR performs better than traditional interpolation methods (IDW and Kriging). IDW presupposes that the similarity and relevance between neighbors are related to their distance (Yasrebi et al. 2009). Similar to IDW, Ordinary Kriging weights the nearby observations to get a forecast of the position of the missing data (Abushandi and Abualkishik 2020). It is clear that the quantity, density, and spital structure of points clearly impact the performances of IDW and Ordinary Kriging. The high accuracy of the improved MLS-SVR approach can be attributed to the high correlation of deformation and the secondary variables such as ground water, compressible thickness, etc., and these variables provide additional information to improve the prediction accuracy. In addition, ISODATA helped to create relatively homogeneous areas with similar distribution of ground water, thickness of compressible thickness, and other factors. And in the interior homogeneous areas, it is easier to understand how land subsidence and contributing variables relate to one another.

One of the innovations of this study was that we combined the temporal-domain deformation estimation and the spatial-domain prediction with multi-output machine learning algorithm. In temporal domain, the framework considers nonlinear characteristics of deformation time series using PCF model. In spatial domain, the framework considers the influence of secondary variables on subsidence. Unlike usual formulation of machine learning that a single value is predicted for each sample, multi-output machine learning algorithms support predicting cross-related multiple variables simultaneously (Xu et al. 2013). In our case, the MLS-SVR model also grabs the connections between the dependent variables in addition to learning the association between the influencing factors and the nonlinear deformation curves.

It must be stated that the use of this model may be limited when the evolution patterns of deformation are more complex and diverse (such as short-term seasonal characteristics).

6. Conclusion

In this study, we proposed a methodological framework for reconstructing spatially continuous time-series surface deformation by combining PS-InSAR, PCF, ISODATA, and MLS-SVR algorithms. We successfully captured the nonlinear characteristics of deformation and generated time-series deformation with 120 m spatial resolution in Beijing Plain area. The major findings are as follows:

1. From 2003-2020, the maximal deformation rate was −138.53 mm/year in Beijing Plain. The PS-InSAR measurements from multi-platform SAR datasets agree well with leveling observations, with R² >0.95 and RMSE <7.3 mm/year.

2. The nonlinear characteristics of time-series deformation are obvious in subsidence area (rate >15 mm/year) in Beijing Plain. Cubic polynomial curve fitting is more accurate in modeling time series deformation than linear, quadratic and quartic models (R² ≈0.99).

3. The Fishnet partitioning shows that 62.13% of the 120 m grids have no PS points. Based on 40 feature variables, the improved MLS-SVR algorithm with combination of ISODATA was used to estimate multiple cubic polynomial coefficients simultaneously, and then used to predict time-series deformation in these grids. The results showed that ISODATA can considerably improve the performance of MLS-SVR model. The reconstructed time-series deformation had good accuracy, with R² reaching 0.92 and MAPE lower than 15%. Compared to IDW and Ordinary Kriging interpolation method, our method yielded higher prediction accuracy, and can better reflect the spatial variations of deformation of ground objects.

4. Annual groundwater levels are the most influential variables for spatial prediction of time-series deformation, followed by the compressible soil layer thickness.

In summary, our research provides a new applicable method for time-series deformation reconstruction and effectively improves the spatiotemporal resolution of subsidence information in Beijing area.

Acknowledgments

The work is supported by National Natural Science Foundation of China (Grant No. 42271487, 42071396, 41930109), Beijing Science and Technology Plan (Grant No. Z201100006720001), Beijing Natural Science Foundation (Grant No. 5172002), and Beijing Outstanding Young Scientist Program (Grant No. BJJWZYJH01201910028032).

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

Upon reasonable request, the corresponding author, Professor Xiaojuan Li, Ph.D., can provide the data that support study’s conclusions.

Supplementary material

Supplemental data for this article can be accessed online at https://doi.org/10.1080/15481603.2023.2230689.

Additional information

Funding

The work was supported by the National Natural Science Foundation of China [42271487]; Beijing Outstanding Young Scientist Program [BJJWZYJH01201910028032]; Beijing Science and Technology Plan [Z201100006720001]; National Natural Science Foundation of China [41930109]; Beijing Natural Science Foundation [5172002]; National Natural Science Foundation of China [42071396]