Study on Land Subsidence Simulation Based on a Back-Propagation Neural Network Combined with the Sparrow Search Algorithm

Abstract

Rapid simulation of land subsidence can provide an effective means of facilitating its management and control. This paper innovatively introduced a back-propagation (BP) neural network subsidence simulation model which was optimized by the sparrow search algorithm (SSA), hereinafter referred to as the SSA-BP model, to simulate land subsidence in Tongzhou District, Beijing. First, based on the time series interferometric synthetic aperture radar (InSAR) monitoring, different technologies such as spatial analysis, Google Earth Engine (GEE), and machine learning were utilized together with a variety of multi-source spatial data, including groundwater level, compressible layer thickness, Visible Infrared Imager Radiometer Suite (VIIRS) nighttime light images, and the OpenStreetMap (OSM) road distribution. Furthermore, we determined the optimal grid scale for land subsidence research and worked out a multifactor-driven SSA-BP land subsidence simulation model for which sensitivity analysis was performed with great care. Main conclusions: (1) From November 2010 to January 2020, the average annual surface displacement rate in Beijing’s subcentre, Tongzhou, ranged from −133.9 to +3.9 mm/year. (2) Our SSA-BP land subsidence simulation model allows for a relatively high overall accuracy. The average root mean square error (RMSE) is 5.00 mm/year, the average mean absolute error (MAE) is 3.80 mm/year, and the average coefficient of determination (R²) is 0.83. (3) Sensitivity analysis shows that the SSA-BP model driven by multiple factors has a far better simulation performance than the model driven by any single weighting factor, and displays the highest sensitivity to the groundwater level factor among all the weighting factors. In terms of subdividing the study area, our SSA-BP land subsidence model runs with multifunctional zones displayed a higher accuracy than that without. This paper made use of a machine learning model to simulate land subsidence in Beijing’s Tongzhou District and explored the sensitivity of different factors to land subsidence, which is helpful for its scientific prevention.

1. Introduction

Land subsidence is a geological phenomenon that causes the surface elevation to depress within a certain area [1,2] as a consequence of natural activities and human practices. At present, more than 200 countries and regions around the world are experiencing land subsidence, which is especially observed in central Mexico, northern Italy, the California area of the United States, central Indonesia, and the BTH (Beijing-Tianjin-Hebei) area in China [3,4,5,6,7,8,9]. For this reason, efforts have been made by scholars worldwide to investigate this global problem. In the Chaussard et al. [10] study of land subsidence that occurred in the central Mexico region, results showed that groundwater pumpage, needed for urban development and agricultural production in particular, was recognized as the principal cause. Da Lio and Tosi [11] studied land subsidence in the Friuli Venezia Giulia coastal plain of Italy, used ERS1/2 and Envisat ASAR interferometric datasets, and identified the variability in the spatial gradient of land subsidence. Hu et al. [12] explored land subsidence in Los Angeles from 2003 to 2017, used small baseline subset interferometric synthetic aperture radar (SBAS-InSAR) technology and determined that the subsidence in this area was mainly caused by the overexploitation of groundwater. Taftazani et al. [13] studied the relationship between groundwater pumping and land subsidence in Jakarta, Indonesia, providing a reference for groundwater exploitation and piped water supply planning in this area.

Land subsidence occurrence in China is predominantly observed within the North China Plain, the Yangtze River Delta, the Fen-Wei Basin and the Pearl River Delta, and presents an array of diverse forms among different regions [2]. In particular, land subsidence in the Beijing Plain region has been recorded for more than half a century [14]. Due to the rapid expansion of urban clusters and other complex factors, together with the world’s largest groundwater depression bowl that occurred and developed in this area, the Beijing Plain has become one of the places which suffer the most severe land subsidence in the world. In the past 20 years, scientists have introduced InSAR measurements, benchmark monitoring and other cutting-edge techniques to perform multidisciplinary research on land subsidence in this area. Foroughnia et al. [15] explored the application of persistent scatterer interferometric synthetic aperture radar (PS-InSAR) in a study of land subsidence in Tehran, Iran, and proposed a large baseline iterative method based on PS-InSAR for the analysis of datasets. Cigna and Tapete [16] used differential interferometric SAR (InSAR) and small baseline subset (SBAS) methods to provide an SAR survey for the Aguascalientes Valley to learn the structurally controlled land subsidence which was caused by groundwater withdrawals. The two researchers reviewed the development process of SBAS-InSAR over the past 30 years and summarized SBAS-InSAR applications in the field of land subsidence and other geological disciplines. Recent studies proved that up to 2020, areas of accumulated subsidence greater than 100 mm cover more than 4000 km² in the Beijing Plain, with the maximum accumulated subsidence exceeding 2 m [17,18].

Tongzhou District is experiencing some of the worst cases of land subsidence in all of Beijing. Field measurements revealed that in 2016, the maximum subsidence rate in Tongzhou was 159.6 mm/year (Taihu, 2012). Land subsidence over 50 mm a year was observed in an area of 169.69 km². To mitigate Beijing’s “big city disease”, the municipal government adjusted the urban spatial pattern in September 2012 by building a new subcentre in Tongzhou District [19]. Subsequently, the numbers of people and buildings per square kilometre in this area increased dramatically, further aggravating the uneven distribution of subsidence in the study area. A series of geological disasters can be induced by severe land subsidence, such as ground collapse, ground fissure formation, and seawater intrusion, forming a disaster chain [20,21,22,23]. Furthermore, regionally uneven subsidence has a certain impact on urban residential safety, causing problems such as the rupture of underground pipelines and cracks in buildings [24,25].

Land subsidence simulating and modelling include physical models, data-driven models and geospatial-artificial intelligence (GEO-AI) models [26,27]. Physical models require many basic geological and hydrogeological parameters to accurately describe complex strata and achieve high simulating precision. However, it is not easy to obtain this kind of prior knowledge. Meanwhile, data density and calculation efficiency might be relatively low, making it even more difficult to apply physical models in practice. In contrast, data-driven models (such as the grey model (GM) and modified GM) are based on time series data and data mining, and short-term simulations are then output. Although these models provide good results, they have difficulty providing long-term simulations for prevention and control [28].

Land subsides due to various factors (groundwater pumping, dynamic loads, etc.) [26,29]; AI methods can achieve good results when applied to nonlinear problems [30,31]. Therefore, the use of machine learning to simulate land subsidence has become a frontier research direction in recent years. Taravatrooy et al. [32] fused the k-means method, soft computing models and other approaches. The research further demonstrated that a fusion method was more accurate than the traditional methods. Shi et al. [33] studied the nonlinear relationships between ground subsidence in the Beijing Plain and the corresponding influential factors using linear regression and principal component analysis. XGBoost was also used to predict ground subsidence with good accuracy. Arabameri et al. [34] used artificial neuronal networks (ANN), bagging and ANN-bagging to derive land subsidence in the Semnan Plain area in Semnan Province, Iran. The results showed that ANN-bagging obtained the highest accuracy, followed by the ANN and bagging. Additionally, all models generally indicated that land subsidence was most likely to occur in the central portion of the study area. Rafiei Sardooi et al. [35] studied the land subsidence on the Rafsanjan Plain in southern Iran. A comparison showed that the accuracy of the studied machine learning method was remarkably higher than that of a statistical model in land subsidence prediction.

The Sparrow Search Algorithm (SSA) is a population-based optimization algorithm which was proposed based on the foraging and anti-predatory behaviours of sparrow populations, and built upon existing population intelligence algorithms, such as Grey Wolf Optimization (GWO) and Particle Swarm Optimization (PSO). It presents certain advantages in terms of velocity, stability, and convergence accuracy. The back-propagation neural network (BPNN) algorithm optimized by SSA-BP offers a strong capacity for searching and a fast convergence speed, therefore it can compensate for the weak points of the traditional BPNN [36]. Thus, the SSA-BP algorithm is innovatively introduced into land subsidence simulation in this paper. Moreover, the density of persistent scatterers (PSs) varies with the category of urban functional areas. Therefore, based on long-time-series InSAR monitoring, a novel approach that combines geological conditions and urban functional areas is proposed on the basis of spatial analysis and mathematical statistical techniques, which are both used to determine the optimal grid size for subsidence research. Furthermore, multiple platforms (such as the Google Earth Engine (GEE)) and multisource spatial data (such as groundwater level, compressible layer thickness, nighttime light, and OpenStreetMap (OSM) data) are employed in the paper as the inputs of the SSA-BP model to simulate land subsidence. Finally, a comprehensive analysis to test the sensitivity of the SSA-BP subsidence simulation model from multiple perspectives is provided in this paper, with weight factors, urban functional areas, and geological conditions being considered.

2. Materials

2.1. Study Site

Tongzhou District spans an area of 906.28 km² in the southeastern Beijing Plain. The regional coordinates range from 39°36′ to 40°02′ north latitude and 116°32′ to 116°56′ east longitude. The annual average temperature here is 11.3 °C [19]. Flat and slightly undulating, local terrain inclines gently from northwest to southeast [14].

Fractures in this area include the Nanyuan-Tongxian fracture, Yan Jiao fracture, and Zhang Jiawan fracture (Figure 1). Considering the distribution of fractures, the geological tectonic units encompassed within this area include the Huangcun Bulge (ΙΙ_2 ΙΙΙ_7 ΙV_16), Mizidian New Depression (ΙΙ_2 ΙΙΙ_8 ΙV_18), Shunyi Depression (ΙΙ_2 ΙΙΙ_6 ΙV_13), and Niubaotun-Dasungezhuang Depression (ΙΙ_2 ΙΙΙ_7 ΙV_17). Quaternary sediments are distributed mainly in the middle and lower portions of the alluvial-pluvial fans, with interlayered sand, gravel and clay. The paleotopography controlled the thickness of the Quaternary sediment [17], which is less than 300 m in Taihu-Yongshun and gradually increases to the east/southeast where the maximum thickness exceeds 500 m. Types of groundwater mainly contain poring water in unconsolidated Quaternary sediment, phreatic water, shallow confined water (shallower than 100 m), and deep confined water (deeper than 100 m) (

Table 1. Information on Quaternary aquifer layers in Tongzhou District, Beijing.

Aquifer	Main Lithology	Bottom Depth
First aquifer (Ι) (phreatic)	Fine sand dominates, with less medium sand and coarse sand, followed by gravel.	40–50 m
Second aquifer (ΙΙ) (first confined aquifer)	Fine sand dominates, with less medium sand and coarse sand, followed by gravel.	80–100 m
Third aquifer (ΙΙΙ) (second confined aquifer)	Fine sand dominates, with less medium sand and silty sand, followed by medium-coarse sand.	150–180 m
Fourth aquifer (ΙV) (third confined aquifer)	Fine sand dominates, with less medium sand and silty sand, followed by medium-coarse sand.	300 m

). Tongzhou’s groundwater can be vertically categorised into four aquifers, and groundwater that is mainly pumped and utilized is no more than 300 m [19].

2.2. Materials

In this study, remote sensing, geological, hydrogeological and levelling data were used (

Table 2. Summary of the datasets.

Datasets		Data Source	Date Range
Remote sensing	RADARSAT-2	Canada Centre for Mapping and Earth Observation	22 November 2010–21 October 2016 25 January 2017– 10 January 2020
	Nighttime light VIIRS	GEE platform	2012 to 2020 on April 1st of each year
	SRTM DEM	Consortium for Spatial Information
Geological and hydrogeological	Groundwater flow field information	Beijing Institute of Hydrogeology and Engineering Geology	2010 to 2019
	Thickness distribution of the compressible layer
	Division of regional geological tectonic units
OSM data *		OSM	1 January 2012–1 January 2020
Beijing urban planning data		Beijing Municipal Commission of Planning and Natural Resources	2016 to 2035

* ©OpenStreetMap contributors.

). The remote sensing data includes radar data, Visible Infrared Imager Radiometer Suite (VIIRS) nighttime light data and Shuttle Radar Topography Mission Digital Elevation Model (SRTM DEM) data. Geological and hydrogeological data involved groundwater flow field information, the thickness distribution of the compressible layer, and the division of regional geological tectonic units. Other data such as open-source OSM road data and Beijing’s urban planning data were also considered.

2.2.1. Radar Image Data

In this paper, 81 scenes of RADARSAT-2 data (Figure 2) from 2 November 2010 to 10 January 2020 were collected to quantify the surface displacement throughout the study area by using time series InSAR technology and GAMMA software.

Table 3. SAR data parameters.

Parameter	RADARSAT-2 ①	RADARSAT-2 ②
Polarization	VV	VV
Orbit direction	De-orbiting VV	De-orbiting VV
Band (wavelength)	C-band (5.6 cm)	C-band (5.6 cm)
Repeat observation period	24 days	24 days
Spatial resolution	30 m	5 m
Track height	798 km	798 km
Incidence angle	20–45°	29.1°
Number of images	56	25
Date range	22 November 2010–21 October 2016	25 January 2017–10 January 2020

contains the specific parameters of the radar imagery.

2.2.2. Basic Geological and Hydrogeological Data

The geological tectonic units in the study area include the Huangcun Bulge, Mizidian New Depression, Shunyi Depression and Niubaotun-Dasungezhuang Depression. The accuracy of the compressible layer thickness reaches 10 m (Figure 3). The accuracy of the groundwater level information reaches 5 m, which was acquired from four layers of groundwater level contours (Figure 4).

2.2.3. Nighttime Light Data

Nighttime lights can monitor the development of urban areas and provide an important perspective which reflects the level of urbanization [37]. As land subsidence also exhibits a strong response to urban construction [38], a static load to some extent, VIIRS nighttime light data can be integrated to simulate the occurrence of land subsidence. The VIIRS dataset, archived from April 2012 to the present, is currently one of the most widely used remotely sensed nighttime light datasets. The data provide a 463 m spatial resolution for VIIRS images, and the VIIRS sensor has a strong detecting capability which enable it to detect subtle reflected light from the Earth’s surface. Using the GEE platform, VIIRS nighttime light data were collected on 1 April of each year from 2012 to 2020 (Figure 5).

2.2.4. OSM Data

In this research, road vector layers out of OSM (http://download.geofabrik.de.html, accessed on 3 March 2022) from 1 January 2012 to 1 January 2020 were employed to describe the effect of dynamic loads on land subsidence (Figure 6).

3. Methods

As noted in previous studies [33,39,40,41], land subsidence has a complex formation mechanism and nonlinear characteristics; therefore, an SSA-BP land subsidence simulation model driven by multiple factors was developed in this study. Based on time series of RADARSAT-2 satellite images in combination with basic geological and hydrological data, the optimal grid scale of the study area was determined through spatial and statistical analyses [2]. Factors with different weights, such as the groundwater level, compressible layer thickness, VIIRS nighttime light distribution, and OSM road distribution in Tongzhou District, Beijing, were used to simulate land subsidence and analyse the sensitivity of various factors in the SSA-BP model (Figure 7).

3.1. InSAR

Ferretti, Prati and Rocca [39] first proposed the PS-InSAR technique, which compensates for the weak points of differential interferometry. The PS-InSAR technology takes surface deformation information by tracking points with stable scattering characteristics (e.g., permanent scatterers (PSs), such as buildings, roads, and exposed rocks) during image processing, thereby effectively avoiding the incoherence caused by long temporal baselines. However, PS point selection with PS-InSAR technology is not ideal in vegetation-covered areas (e.g., nonurban areas), resulting in the inability to obtain complete surface deformation information in some areas.

Berardino et al. [42] proposed small-baseline interferometry (SBAS) in 2002 and 2004. This technology which integrates several small-baseline interferometry datasets has improved the temporal resolution of surface deformation monitoring, effectively reduced the impacts of temporal and spatial incoherence, thus providing a new technical support for monitoring land subsidence [43].

In this paper, RADARSAT-2① (2011–2016) and RADARSAT-2② (2017–2020) images covering the study area were selected. With GAMMA software, the SAR data from these two time-series were used to obtain the temporal line-of-sight (LOS) information of ground deformation. From November 2010 to October 2016, the PS-InSAR method was applied to a single master image (acquired on 25 November 2014) to calculate the surface deformation in the study area. The registration accuracy reached the required 1/8th of a pixel level. The space-time baseline of RADARSAT-2① is shown in Figure 8.

SBAS-InSAR technology was employed to take surface deformation information within the study area from January 2017 to January 2020. The SAR data collected after 2017 were 25 extra-fine images, which were too few in number to support PS-InSAR. SBAS-InSAR was thus chosen to acquire more coherent points. This provided a solution to the issue of limited point availability for monitoring land subsidence. The spatial baseline threshold was 500 m, and the temporal baseline threshold was 100 days (4 revisit cycles). Based on these thresholds, 45 small-baseline interference pairs were established. The amplitude deviation threshold was set to 0.6. The small-baseline interference pairs for RADARSAT-2② are shown in Figure 9.

Acting according to the method described by Guo [44], time series fusion of the deformation results of RADARSAT-2① and RADARSAT-2② SAR images was performed in the following steps. First, the nearest point method was applied with RADARSAT-2① as reference data to find the corresponding points in the image monitoring deformation results of RADARSAT-2②, and the cumulative subsidence from the RADARSAT-2① reference baseline from 22 November 2010 to 21 October 2016 was used as the initial baseline for the deformation sequence in the next time period (RADARSAT-2②, 2017–2020) for time series fusion [45].

Finally, the time-series ground deformation information from January 2010 to January 2020 in the study area was obtained. Along the radar LOS direction, PS-InSAR was used to measure the deformation component. Therefore, according to Deng et al. [46], the InSAR-derived surface displacement LOS direction was converted into the vertical direction.

3.2. Regional Division and Determination of the Optimal Grid Scale

First, based on the administrative divisions, groundwater levels and the groundwater aquifer boundaries of Tongzhou District, overlay analysis was used to determine the extent of the study area in ArcGIS software. Second, after taking into consideration the following three factors, namely the fracture distribution, the division of geological tectonic units, and the detailed regulatory plans for Tongzhou New City, the study area was divided into grids, each of which has varied resolutions, and land subsidence in different grids was separately simulated. Figure 10 presents the specific process, and Figure 11 indicates the specific types of multifunctional areas identified in Tongzhou District.

Finally, a grid division experiment was implemented for each subregion. Fourteen subregions were further divided into multiple grid scales using the Fishnet tool. The smallest grid scale was 30 m (in line with the accuracy of the RADARSAT imagery), and the scale was gradually increased to 960 m. Finally, the optimal grid size of each subregion was determined by means of a statistical method.

As shown in Figure 12, the proportion of PSs for each zone in the optimal grid exceeds 80%.

3.3. Acquisition of Weighting Factors

In this paper, the compressible layer thickness, groundwater level, nighttime light distribution and road network were used as the inputs of the SSA-BP subsidence model. The selection of weighting factors for these variables was based on previous research, and the compressible layer thickness and groundwater level were chosen because they are the two factors which produce the most significant impact on land subsidence [26]. In addition, with respect to the impact of dynamic/static loads on land subsidence [47,48], nighttime light distributions may indicate the degree of urbanization, providing indirect information on static load factors such as city size, population, and building density. The road network, such as traffic flows and rail transportation, on the other hand, makes up dynamic load factors to a certain extent. Therefore, these two factors were also incorporated in the model [49]. Below is a description of methods this paper used to work out different weighting factors for each grid.

3.3.1. Compressible Layer Thickness

With the ‘features to raster’ tool in ArcGIS, definable compressible layer thickness contours were converted into raster data, and the compressible layer thickness data were acquired.

3.3.2. Groundwater Level

Steps to acquire the groundwater level in each grid were as follows. First, an irregular triangulation (TIN)-based mesh of groundwater contours was created. The TIN mesh was then transformed into raster data. Next, ArcGIS software was used to change the grid data that have a resolution of 30 m into a point layer. Finally, the intersection tabulation tool was used to turn the average groundwater level in each grid [44,50].

3.3.3. Nighttime Lights

The nighttime light data in this paper were processed in the following way. First, though the GEE platform, the nighttime light data in the study area were acquired as a raster image with a spatially resampled resolution of 30 m. Then, the cumulative value of light pixels in each grid of the study area was obtained through spatial analysis and a statistical method.

3.3.4. Road Distribution

Weights of different roads were assigned on account of two facts: the classification of road data in OSM, and the corresponding average speeds of travel on Chinese roads, as shown in

Table 4. Road weights.

Road	Speed (km/h)
Highway	50
Main road	40
Bikeway	15
Pavement	5
Residential road	5
Subway	50
Railway	40

.

3.4. SSA-BP

3.4.1. BP Neural Network

BPNNs provide a variety of advantages in terms of simple structure and nonlinear fitting ability and are widely used to assess atmospheric haze, traffic flows and other factors. However, slow convergence speed, varied prediction and training abilities are some of the disadvantages.

3.4.2. SSA

Xue and Shen [36] proposed the sparrow search algorithm, a new swarm intelligence-based optimization approach, which permits strong optimization ability, high search accuracy and a fast convergence speed. Compared to other swarm intelligence algorithms, the SSA algorithm outperforms the PSO algorithm in terms of convergence capability and surpasses the Ant Colony Optimization (ACO) algorithm in search speed. The SSA has been innovatively deployed to optimize a BPNN and simulate land subsidence in Tongzhou District.

The SSA procedure is as follows. ① Initialize the population, the number of iterations, and the initial proportions of the predator and participants. ② Calculate the fitness values, and then sort them. ③ Update the predator position by using Formula (3). ④ Update the scrounger position by using Formula (4). ⑤ Update the producer position by using Formula (5). ⑥ Calculate the fitness values and update the positions of the sparrows. ⑦ Exit when the stop condition is met and output the result; if the condition is not met, repeat Steps ②–⑥.

The main formulas are as follows.

The positions of sparrows are:

$$\mathrm{X}=\left[\begin{array}{cccc}{x}_{1,1}& x_{1,2}& \dots & x_{1,r}\\ x_{2,1}& x_{2,2}& \dots & x_{2,r}\\ ⋮& ⋮& ⋮& ⋮\\ x_{m,1}& x_{m,2}& \dots & x_{m,r}\end{array}\right]$$

(1)

where m is the number of sparrow and r represents the variable dimension to be optimized.

All sparrows’ fitness values are denoted by:

$$F_x=\left[\begin{array}{c}\begin{array}{c}f\left(\left[x_{1,1}{x}_{1,2}\dots x_{1,r}\right]\right)\\ f\left(\left[x_{2,1}{x}_{2,2}\dots x_{2,r}\right]\right)\end{array} \\ \begin{array}{c}⋮\\ f\left(\left[x_{m,1}{x}_{m,2}\dots x_{m,r}\right]\right)\end{array} \end{array}\right]$$

(2)

where the fitness value of each individual is shown in each row.

According to the rules, the producer’s location is renewed during every iteration:

$$X_{i,j}^{k+1}=\left\{\begin{array}{c}{X}_{i,j}^k\times \exp \left(-\frac{i}{\alpha \times ite{r}_{max}}\right) if R_2<ST\\ X_{i,j}^k+Q\times L if R_2\ge ST\end{array}\right.$$

(3)

where k is the present iteration, $X_{i,j}^k$ is the value of the jth dimension of the ith sparrow in the kth iteration, iter_max is the largest number of iterations, and α is a random number (α ∈ (0, 1]). The alarm value R² (R² ∈ [0, 1]) and safety threshold ST (ST ∈ [0.5, 1.0]) must also be set. The random number Q obeys a normal distribution. In matrix L (1 × d), every element is 1. The scrounger position is updated with the following formula:

$$X_{i,j}^{k+1}=\left\{\begin{array}{c}Q\times \exp \left(\frac{X_{worst}^k-X_{i,j}^k}{i^2}\right) if i>n/2\\ X_p^{k+1}+\left|X_{i,j}^k-X_p^{k+1}\right|\times A^{+}\times L otherwise\end{array}\right.$$

(4)

where X_P represents the producer’s optimum position and X_worst is the present worst location. In matrix A (1 × d), every element is 1 or −1, and $A^{+}=A^T{\left(A{A}^T\right)}^{-1}$.

Assuming that 10% to 20% of sparrows are aware of danger, the mathematical model is as follows:

$$X_{i,j}^{k+1}=\left\{\begin{array}{c}{X}_{best}^k+\beta \left|X_{i,j}^k-X_{best}^k\right| if f_i>f_g\\ X_{i,j}^k+\mathrm{K}\left(\frac{\left|X_{i,j}^k-X_{best}^k\right|}{\left(f_{i-}{f}_w\right)+\epsilon }\right) if f_i=f_g\end{array}\right.$$

(5)

where X_best denotes the current best location, β is used for step size control, K is a random number, and K ∈ [−1, 1]. Here, f_i is the present fitness value of the selected sparrow. Additionally, f_g is the best fitness value at present, and f_w is the worst fitness value at present. The constant ε is a minimum value used to avoid the singularity caused by dividing by zero.

Algorithm 1 shows the specific parameters used in the simulation of land subsidence.

Algorithm 1 The framework of the SSA

Input:

G: the maximum number of iterations——20

PD: the number of producers——14

SD: the number of sparrows that perceive danger——4

R2: the alarm value——0.6

m: the number of sparrows——20

The sample points in the training set and test set account for 80% and 20% of all sample points, respectively.

3.4.3. Parameter Determination for the SSA-BP Model

To verify the reliability of the SSA-BP model, experiments were performed for all 14 subregions (hereinafter referred to as zones), one at a time.

Eighty percent of the research samples were randomly categorized into training samples, and twenty percent were used as testing samples. The RMSE, MAE and R² were selected as the evaluation indicators to verify the SSA-BP model. Through many experiments, the approximate ranges for the model parameters were determined.

Table 5. Experiments with different parameter combinations in the SSA-BP model.

n	PD:SD	R2	RMSE	MAE	R2
20	7:2	0.6	0.86	0.65	0.92
40	7:2	0.6	1.08	0.88	0.85
60	7:2	0.6	0.94	0.72	0.89
20	7:4	0.6	0.96	0.72	0.91
20	5:2	0.6	0.95	0.73	0.89
20	7:2	0.4	0.96	0.76	0.90
20	7:2	0.8	1.00	0.73	0.89

shows the results for zone F. When n = 20, PD:SD = 7:2, and R₂ = 0.6, the model yielded the highest accuracy.

4. Results

In view of the PS-InSAR monitoring results, this paper employed the ‘spatial join’ function in ArcGIS software to acquire the time series subsidence information for each grid from 2012 to 2019. The surface displacement information and weighting factors in zones A–N were then substituted into the SSA-BP model to simulate land subsidence.

4.1. InSAR Results and Validation

According to the InSAR results (Figure 13), a total of 128,379 PSs were attained, making the distribution density arrive at 154 points/km². The vertical displacement rate in Tongzhou District, Beijing, was between −133.9 mm/year and +3.9 mm/year. The average was −19.6 mm/year.

Throughout the study area, the area with an annual displacement rate of less than −100 mm/year was 2.6 km², the area with a displacement rate between −100 mm/year and −75 mm/year was 38.6 km², the area with a displacement rate between −75 mm/year and −50 mm/year was 79.8 km², the area with a displacement rate from −50 mm/year to −25 mm/year was 71.9 km², and the area with a rate over −25 mm/year was 510.2 km², which accounted for 0.4%, 5.9%, 12.2%, 11.0% and 70.4% of the study area, respectively. The specific ratios are shown in

Table 6. Area ratios of the InSAR results.

InSAR Result (mm/year)	Area Ratio (%)	Area (km2)
< −100	0.4	2.6
−75 to −100	5.9	38.6
−50 to −75	12.2	79.8
−25 to −50	11.0	71.9
> −25	70.4	510.2

.

Data from 11 levelling benchmarks were used for validation by means of comparing the InSAR value of the PS point that is closest with each corresponding benchmark. The InSAR results and levelling benchmarks were found to be highly consistent. From 2015 to 2016, the minimum absolute error was 0.89 mm/year, the maximum absolute error was −12.3 mm/year, the average error was −3.4 mm/year, and the RMSE was 7.3 mm/year. The correlation between InSAR measurements and levelling data was greater than 0.95 (Figure 14).

4.2. Accuracy Analysis of the SSA-BP Subsidence Simulation Model

Figure 15 and Figure 16 show the accuracy of the model’s performance in each zone. Generally speaking, the accuracy of the SSA-BP land subsidence simulation model in each zone was relatively high. The average RMSE of the 14 zones was 5.00 mm/year, the average MAE was 3.80 mm/year, and the average R² was 0.83. The RMSE among the 14 zones ranged between 0.32 mm/year and 13.87 mm/year, the MAE ranged between 0.24 mm/year and 10.61 mm/year, and R² ranged between 0.49 and 0.94. The model presented a sound overall accuracy and fitting degree, which confirm that the SSA-BP model is feasible in simulating land subsidence.

Among those in all the regions, the simulated values in zones B, C, F, I and J were only slightly different from the actual values. The RMSEs in zones B, C, F, I and J were 1.06 mm/year, 1.36 mm/year, 0.86 mm/year, 0.96 mm/year and 1.00 mm/year, respectively; the corresponding MAEs were 0.82 mm/year, 1.10 mm/year, 0.65 mm/year, 0.72 mm/year and 0.72 mm/year; and the related R² values were 0.87, 0.84, 0.92, 0.89, and 0.80. These values reflect a high degree of fit. Moreover, information on the groundwater level in this area was abundant, so the data meet the input requirements of a GEO-AI model. Therefore, the model achieved stable simulation accuracy and a high fitting degree in these subregions.

The simulation results in zones A, K, M and N featured small errors compared with the observation values, with RMSEs of 1.40 mm/year, 1.32 mm/year, 1.12 mm/year and 0.32 mm/year, respectively, and corresponding MAEs of 1.11 mm/year, 0.96 mm/year, 0.73 mm/year and 0.24 mm/year; the fitting degree, however, was poor, with R² values of 0.65, 0.65, 0.65 and 0.49, respectively.

The simulation results in zones D, E, G, H, and L were relatively far from the true values. The RMSEs in zones D, E, G, H, and L were 9.60 mm/year, 10.99 mm/year, 12.46 mm/year, 13.66 mm/year and 13.87 mm/year, respectively, and the corresponding MAEs were 7.76 mm/year, 7.93 mm/year, 9.36 mm/year, 10.41 mm/year, and 10.61 mm/year. However, the fitting degrees in zones D, E, G, H and L were good, with R² values of 0.96, 0.94, 0.85, 0.84, and 0.88, respectively.

The results of the SSA-BP land subsidence simulation model in the study area yielded small errors and a high degree of fitting. The overall accuracy was excellent, making it reliable for the simulation of land subsidence.

5. Discussion

5.1. SSA-BP Subsidence Simulation Model Results

In the fourteen subregions, the model achieved high accuracy and a high fitting degree in zones B, C, F, I and J, whereas it achieved good accuracy but a poor fitting degree in zones A, K, M and N, and poor accuracy but a high fitting degree in zones D, E, G, H and L.

Good results were obtained in zones B, C, F, I and J, most likely because these zones exhibited only small changes in the land subsidence rate and subsidence gradient [51].

The reasons for the results in zones A, K, M and N may be as follows: From the perspective of the subsidence gradient, zones A, M and N correspond to the periphery of the urban subcentre in Tongzhou District, with a low level of urbanization, a low groundwater pumping intensity, a slow rate of land subsidence, and a low average subsidence rate from 2012 to 2019. Although zone K is located in the urban subcentre of Tongzhou District, it has developed slowly, and the annual average subsidence rate is slow. Therefore, the RMSEs and MAEs in zones A, K, M and N were small.

For zones D, E, G, H, and L, the results may be attributed to the following factors: (1) In Figure 16, the accuracy of the SSA-BP simulation model is shown, where the sedimentation rates vary greatly among zones D, E, G, H and L because they are located along the edges of the Dongbalizhuang-Dajiaoting funnel, where the sedimentation gradient varies to a large extent. Therefore, errors in this zone are comparatively high. (2) Machine learning models are data driven and thus tend to be susceptible to the sample size [52]. Generally speaking, if the input sample size is small, the accuracy of the model is low. In this study, there were only 200 samples in zone E, so the accuracy of the SSA-BP model in this area was not as good as expected. Moreover, the original PSs in zone D were unevenly distributed, with more points in the north and fewer points in the south. A sample distribution such as this one is unfavourable for accurate land subsidence modelling.

5.2. Sensitivity Analysis of the SSA-BP Land Subsidence Simulation Model

To explore the applicability of the SSA-BP land subsidence simulation model, multiple sensitivity analyses were performed [53] in this study. Three major conditions that might affect the model’s sensitivity were considered as follows: (1) the influences of different weighting factors on the model’s accuracy; (2) the influences of multifunctional zones of Beijing’s subcentre; and (3) the effect of the zonal division on the simulation accuracy of the model.

5.2.1. Sensitivity of the Model Accuracy to Different Weighting Factors

Different weighting factors were introduced into the SSA-BP model to assess their impact on model accuracy (Figure 17). The accuracy and fitting degree of the SSA-BP land subsidence simulation model driven by multiple factors suggested a better performance than those when the model was driven by a single factor. For the SSA-BP land subsidence simulation model driven by multiple factors, the RMSE was 0.32–13.87 mm/year, the MAE was 0.24–10.61 mm/year, and the R² was 0.49–0.94. For the SSA-BP subsidence model driven only by the groundwater level, the RMSE was 0.33–18.51 mm/year, the MAE was 0.23–14.15 mm/year, and the R² was 0.41–0.93. For the SSA-BP model driven only by the compressible layer thickness, the RMSE was 0.4–33.78 mm/year, the MAE was 0.34–41.48 mm/year, and the R² was −0.04–0.37. For the SSA-BP model driven only by the road distribution, the RMSE was 0.43–39.78 mm/year, the MAE was 0.34–33.74 mm/year, and the R² was −0.02–0.33. For the SSA-BP model driven only by nighttime light information, the RMSE was 0.32–39.09 mm/year, the MAE was 0.24–39.23 mm/year, and the R² was 0.03–0.49.

These results show that the sensitivity of the model was the highest when the groundwater level was used as the sole weighting factor, followed by the compressible layer thickness, when it was used as the one and only weighting factor, and then the road distribution as the lone weighting factor. The sensitivity was the lowest when nighttime light information was used as the single weighting factor. The possible reasons underlying these findings were postulated as follows: (1) According to previous studies [26], the main factors affecting land subsidence are the groundwater level and compressible layer thickness, while the weights of factors related to urban dynamic and static loads are relatively small at only 13.3–25.6%. (2) The water level can be identified an external factor that changes every year, but the compressible layer thickness is an internal geophysical factor that remains unchanged.

In summary, a land subsidence simulation model driven by multiple weighting factors is strongly recommended. When the amount of data is small or data are difficult to obtain, it is highly suggested that groundwater level data need to be prioritized.

5.2.2. Sensitivity to Multifunctional Zones in the City’s Subcentre on the Simulation Accuracy of the Model

This paper assessed the sensitivity of the SSA-BP land subsidence simulation model to different multifunctional zoning schemes in Beijing’s subcentre. As shown in Figure 18, when the factor of district planning was taken into consideration, the average RMSE and MAE of the 14 functional zones were 6.10 mm/year and 4.64 mm/year, respectively, and the R² was 0.86. When district planning was not considered, the model accuracy decreased (RMSE = 20.91 mm/year, MAE = 14.59 mm/year, and R² = 0.81).

It was assumed that the density of PSs and groundwater pumpage could be closely related to the types of urban functional areas. In this paper, a priori knowledge of multifunctional zones was innovatively introduced to the inputs of the model. After this operation, the RMSE decreased by 14.81 mm/year, the MAE decreased by 9.95 mm/year, and the R² increased by 0.05. Accordingly, the accuracy and fitting degree of the model after adding the multifunctional zone factor performed better than those of the model without the factor. Therefore, in future studies of land subsidence, it will be beneficial to add the factor of functional zones for simulation.

We assessed the sensitivity of the SSA-BP land subsidence simulation model to the division of structural units. As shown in Figure 19, when considering the division of structural units, the average RMSE and MAE of the 14 zones were 5.00 mm/year and 3.80 mm/year, respectively, and an R² of 0.83. When the division of structural units was not considered, the model yielded an RMSE of 9.27 mm/year, an MAE of 6.94 mm/year, and an R² of 0.83.

The demarcation of structural units was based on both the thickness of the compressible layer and the fracture zone distribution, which are internal factors that cause land subsidence. Geological tectonic structure, as the basic factor for the occurrence and development of land subsidence, displays clear regional characteristics. In the study area, geological tectonic units include the Huangcun Bulge, Mizidian New Depression, Shunyi Depression, and Niubaotun-Dasungezhuang Depression. When delineating our study area, we took Beijing’s geological tectonic units into consideration, which has effectively improved the simulation accuracy of our model. By considering these geological tectonic units and their distributions, the RMSE decreased by 4.27 mm/year, and the MAE decreased by 3.14 mm/year. Our model that included tectonic units worked better in terms of accuracy than that without considering this feature. To sum up, geological tectonic units of a study area is recommended as a valuable feature for further division.

6. Conclusions

In this paper, we established a novel SSA-BP model based on various weighting factors in an effort to simulate land subsidence. First, time series results of land subsidence were obtained with PS-InSAR technology. Second, the study area was categorized into 14 zones by considering features such as geological conditions and urban functions. We then determined the optimal grid size for our land subsidence study in Beijing’s Tongzhou District. Third, the SSA-BP land subsidence model was built on the basis of the sparrow search algorithm, combining with multisource spatial data (groundwater level, compressible layer thickness, nighttime light, and OSM data) and using new technologies (GIS, GEE and GEO-AI). Then, the accuracy of the model was verified, and its sensitivity was explored from multiple perspectives, including driving factors, urban functional zones, and structural units. The main conclusions of our research are as follows:

(1) The average annual land subsidence rate in Beijing’s Tongzhou District from November 2010 to January 2020 ranged from −133.9 mm/year to +3.9 mm/year, and the annual average subsidence rate was −19.6 mm/year.

(2) According to the SSA-BP subsidence simulation results, the average RMSE in the 14 zones was 5.00 mm/year, the average MAE was 3.80 mm/year, and the average R² was 0.83. The overall accuracy and fitting capability of the model were impressive, and further confirmed the feasibility of using the SSA-BP model for land subsidence simulation.

(3) Our SSA-BP model which simulated land subsidence for study areas revealed that different driving factors produced varied impacts on land subsidence. The model driven by multiple input factors allowed a higher accuracy than that driven by any single factor. Among the many factors considered, groundwater level plays the most important role. Furthermore, the SSA-BP subsidence model that runs with multifunctional zones displayed a higher accuracy than that without.

In this paper, we give a detailed explanation of our applied research that aims to provide an innovative solution to land subsidence simulation with an AI model. From the result, the SSA-BP model can rapidly simulate land subsidence. Such ability is conducive to the exploration and prevention of land subsidence. In future research, we will try more forefront AI models to upgrade the accuracy of our land subsidence prediction models.