Spatio-temporal evolution of land subsidence and susceptibility zonation of associated ground fissures in the urban area of loess plateau: a case in Xianyang city, China

Abstract

For over 50 years, the Xianyang city on the Chinese loess plateau, has been deeply affected by land subsidence and its associated ground fissures. In this study, 67 Sentinel-1A images from 2015 to 2022 were analysed to obtain the spatio-temporal evolution of land subsidence. Subsidence centers identified in Yunyang town and Luqiao town exhibit annual subsidence rates of −31 mm/a and −26.7 mm/a, respectively, culminating in total subsidence of −258 mm in Yunyang town and −139 mm in Luqiao town. Moreover, time series analysis revealed that both towns exhibit pronounced seasonal subsidence patterns. It is believed that human cultivation and groundwater overuse, combined with the effects of active faults, have led to uneven land subsidence. Then, the gradient of land subsidence and its potential relation to ground fissure hazards were evaluated. The results highlight a strong spatial correlation between high subsidence gradient areas and dense ground fissure distribution. Therefore, we proposed for the first time to introduce the land subsidence gradient factors for the susceptibility mapping of ground fissures using the artificial neural network algorithm. The susceptibility zonation results show that integrating the gradient factor can improve the rationality of the susceptibility evaluation. The above understanding provides valuable support for disaster prevention and mitigation in urban areas affected by land subsidence and ground fissures in loess regions.

Keywords:

• InSAR
• land subsidence
• susceptibility
• ground fissure
• ANN

1. Introduction

Land subsidence and ground fissures are two typical slow-onset geological hazards that impact urban development. These hazards are prevalent in numerous countries and regions across all over the world, such as China, Mexico and the United States (Conway 2016; Cigna and Tapete 2021b; Zhou, Yao, et al. 2022). They not only damage engineering structures, buildings, transportation facilities, and land resources, but also lead to severe environmental issues (Peng et al. 2018; Zhang, Qian, et al. 2023). Particularly, regions such as the Fenwei Basin and the North China Plain in China have long been subjected to the hazards of land subsidence and ground fissures, resulting in cumulative economic losses (Jia et al. 2020; Bagheri-Gavkosh et al. 2021). The government has shown a growing interest in the susceptibility zonation of regional land subsidence and its associated ground fissure hazards (Ye et al. 2018). Accurate ground fissure hazard mapping and the identification of essential factors in hazard development are vital for hazard reduction. However, evaluating ground fissure susceptibility is challenging due to its complexity, multidisciplinary aspects, and the uncertainties in modelling. Therefore, it is crucial to develop efficient methods for accurately identifying key factors of ground fissure occurrences and for precisely mapping their susceptibility zones.

Unlike traditional point-based monitoring methods such as GPS and leveling, InSAR technology enables rapid monitoring over a large regional scale (Gao et al. 2016; Zhou et al. 2016; Guo et al. 2020; Bai et al. 2022). Moreover, with the help of time-series InSAR techniques, the spatiotemporal evolution features of surface deformation in the study area can be easily extracted, making it an essential research method in the field of ground subsidence and ground fissure monitoring (Yang et al. 2014; Cigna and Tapete 2021a). Compared to regional scale land subsidence hazards, the accurate identification of ground fissures, as a linear hazard, remains challenging (Aldaajani et al. 2022; Zhan et al. 2022). The traditional susceptibility evaluation of ground fissures primarily relies on geological factors, such as geology, topography, or underlying faults. However, there is a clear correlation between land subsidence and ground fissure (Zhang et al. 2009; Hu et al. 2019), indicating a multi-factor coupling effect in the occurrence of ground fissures. Zang et al. (2021) attempted to introduce land subsidence factors into the zonation of ground fissure susceptibility. Cigna and Tapete (2021a) further incorporated angular distortion indices into the evaluation of ground fissure hazards, improving the accuracy of risk assessment for ground fissures. Therefore, the introduction of land subsidence factors into the evaluation of ground fissure hazards can enhance the capability of hazard zonation to a certain extent (Castellazzi et al. 2017). However, in comparison to spatially uniform subsidence, the harm caused by differential subsidence is more severe (Fernández-Torres et al. 2022). High horizontal stress concentrated in narrow areas may lead to the formation of ground fissures, and the presence of high subsidence gradient helps in their identification (Cabral-Cano et al. 2008; Wang et al. 2019). Therefore, the susceptibility zonation indicators currently used for ground fissures do not fully cover all relevant factors, and their sensitivity is not high enough to adequately reflect ground fissure hazards.

Many studies have shown that both knowledge-driven and data-driven models are effective for zoning areas susceptible to ground fissures (Cao et al. 2022; Wang et al. 2022; Lyu et al. 2023). Knowledge-driven models, such as the analytic hierarchy process and expert scoring, heavily rely on expert knowledge and may be influenced by subjectivity during the evaluation process. Data-driven methods involve analyzing the relationship between acquired data and ground fissure hazards, enabling predictive analysis. Commonly used data-driven methods include probabilistic methods, deep learning, machine learning (Nadiri et al. 2018; Merghadi et al. 2020). In recent years, machine learning models have emerged as powerful tools for multi-factor coupling analysis (Moayedi et al. 2019; Nefeslioglu et al. 2021). Arabameri et al. (2020) employed an artificial neural network (ANN) model for the evaluation and prediction of land subsidence hazards in Iran. Choubin et al. (2019) conducted a study on ground fissure susceptibility zoning in the Yazd-Ardakan Plain using various machine learning techniques. While these studies often focus on comparing the accuracy of different models, they tend to overlook the crucial impact of reasonable factor selection on model precision. Therefore, exploring effective combinations of factors to enhance the accuracy of ground fissure susceptibility zoning is imperative.

In this paper, the Small Baseline Subset (SBAS)-InSAR technique is employed to research land subsidence and ground fissures in the plain area of Xianyang city on the Loess Plateau. The long-time series deformation characteristics of the study area from 2015 to 2022 were extracted. The accuracy of the InSAR results was validated using the results of previous studies (InSAR and leveling data) and field surveys. By integrating geographical information factors and overlaying the subsidence rate and subsidence gradient derived from InSAR, susceptibility zoning for ground fissures was generated. The spatial relationship between ground fissures and geological factors such as land subsidence and subsidence gradient was analyzed, and the primary causes of ground fissure formation were explored. This study primarily aims to assess the benefits of utilizing Sentinel-1A data to gather information on subsidence gradient for evaluating potential ground fissure risks. Such an approach offers crucial perspectives for mitigating and managing ground fissure disasters.

2. Study area

2.1. Geologic and hydrogeologic settings

Xianyang city, located in the main development center of the Guanzhong basin, China, exhibits distinct geomorphic and structural conditions (Figure 1). In terms of geomorphology, Xianyang city can be divided into two parts: the north characterized by the gully region of the Loess Plateau, and the south being the Weihe alluvial plain. The southern plain features a variety of landforms, including loess tablelands, alluvial plains, and groups of alluvial fans. Numerous active faults are predominantly distributed in the southern plain of Xianyang city, include Kouzhen-Guanshan fault, Fufeng-Sanyuan-Pucheng fault, Yidian-Qianxian-Meiyuan fault and Weihe fault. After extensive engineering geological investigations, 101 ground fissures were identified in the plain area of Xianyang city. The lengths of these ground fissures range between 20 and 6017 m, and the total length of all ground fissures is 65510 m (Figure 1). Ground fissures are mainly distributed near the main faults, as shown in Figure 1a. Unlike the ground fissures found in the Xi’an city, those in Xianyang area are scattered with a relatively short extension lengthy (Peng et al. 2020).

Figure 1. Study area of Xianyang city in the Weihe basin: (a) Ground fissure and active fault in Weihe basin and (b) location of the study area, ground water monitoring well and the coverage of Sentinel-1A data.

Since the Cenozoic era, the basin has been deposited with extremely thick loose sediments under the influence of topographic structure. The main sedimentary lithologies include loess, clay and sandy soil (Figure 2). The 300 m depth of plain area can be divided into phreatic aquifer, shallow confined aquifer and deep confined aquifer. Additionally, the plain area of Xianyang city is characterized by a dense distribution of canal irrigation systems, with the major ones being the Baoji gorge irrigation district and the Jinghui canal irrigation district (Figure 2a). Xianyang city plays a crucial role as a major production base for grains and vegetables (Wang et al. 2021).

Figure 2. (a) The distribution of irrigation district in the Xianyang plain and (b) schematic map of aquifer and stratigraphic distribution in the study area (Gao et al. 2022).

2.2. Groundwater exploitation

Xianyang city has gradually expanded in the past few decades, and its population has steadily risen. In 2022, Xianyang city has a permanent population of 4.213 million, with an urbanization rate of 57.3% (Figure 3a). As reported in the water resources bulletin of Shaanxi province in 2015, there were still 92 water supply wells in urban areas, more than 300 wells for industrial and mining enterprises, and more than 2000 wells in rural areas. Within the entire water supply system, groundwater plays a significant role, constituting as much as 73.6% of the total supply (Figure 3b). With the full closure of enterprises’ own wells in urban conditions and the implementation of policies such as diverting water from the Shitou river to the Weihe river, the phenomenon of groundwater overexploitation in Xianyang city has been alleviated to some extent. However, groundwater remains the primary source of water supply in Xianyang city.

Figure 3. (a) The evolution of urbanization rate in Xianyang city and four major districts (counties) (b) proportion of water resources supply in Xianyang city.

Due to the long term over exploitation of groundwater in Xianyang city, its groundwater level has been showing a continuous downward trend. Environmental problems such as ground fissures, building fissures, groundwater quality pollution and land subsidence have emerged (Figure 4), which have seriously affected social and economic development and the production and life of residents.

Figure 4. (a) Farmland damage caused by shuanghuaishu ground fissure, (b) road damage caused by Jinghe fault, (c) farmland damage caused by Kouzhen-Guanshan fault, (d) wall cracks caused by ground fissures in Laonan village, Fuping county, (e) farmland damage caused by Wuxing village ground fissures, Qian county and (f) house cracks caused by ground fissures in Tingzi village, Fuping county.

3. Materials and methods

The comprehensive analysis presented in this paper is segmented into three sections (Figure 5). First, the ground deformation characteristics of Xianyang city from 2015 to 2022 were measured using the SBAS-InSAR method. The accuracy of the InSAR result was verified using leveling data and field investigations results. Subsequently, the subsidence gradient was obtained through gradient analysis. Lastly, the land subsidence rate and subsidence gradient were integrated with multisource data and input into the ANN model to obtain the susceptibility zoning results for ground fissures.

Figure 5. The flowchart of the data processing.

3.1. Ground deformation detection based on SBAS method

SBAS-InSAR is a commonly used time series analysis method for obtaining ground deformation characteristics, as proposed by Berardino et al. (2002) and Riccardo et al. (2004). This method adopts the strategy of obtaining interference pairs from multiple primary images, which improves the time sampling rate of deformation measurement and expands the spatial coverage of the study area. It’s noteworthy that this method has been extensively applied in surface deformation monitoring (Du et al. 2018).

For the acquired N + 1 SAR images, a reasonable spatio-temporal baseline threshold was set, and differential interferometry processing was applied to obtain M interferograms. These generated differential interferograms underwent filtering and phase unwrapping processes. The unwrapped phase of any pixel (x, r) of the j-th differential interferogram is expressed as: (1) $$\Delta {\phi }_j(x,r)=\phi (t_B,x,r)-\phi (t_A,x,r)\approx \frac{4\pi }{\lambda }[d(t_B,x,r)-d(t_A,x,r)]+\Delta {\phi }_{\mathit{\text{topo}}}^j(x,r)+\Delta {\phi }_{\mathit{\text{atm}}}^j(x,r)+\Delta {\phi }_{\mathit{\text{noise}}}^j(x,r)$$(1) where φ(t_A, x, r) and φ(t_B, x, r) represent the phase values of SAR images at t_A and t_B, respectively; d(t_B, x, r) and d(t_A, x, r) are the cumulative deformation in the line of sight direction relative to d(t₀, x, r) at moments t_B and t_A; λ is the wavelength of the radar; Δφ^j_topo(x, r) is the residual topographic phase; Δφ^j_atm(x, r) is the atmosphere delay phase and Δφ^j_noise(x, r) is the noise error phase.

In order to obtain the deformation data, after removing the above phase error, the formula (1) can be expressed as: (2) $$\Delta {\phi }_j(x,r)=\phi (t_B,x,r)-\phi (t_A,x,r)\approx \frac{4\pi }{\lambda }[d(t_B,x,r)-d(t_A,x,r)]$$(2)

Assuming that the deformation rate between different interferograms is V_k, then the cumulative deformation between any period t_A − t_B can be expressed as: (3) $$\Delta {\phi }_j(x,r)=\frac{4\pi }{\lambda }\sum _{k=t_A}^{t_B}(t_{k+1}-t_k)v_{k,k+1}$$(3)

This equation can expresse in matrix form as: (4) $$\mathit{\text{Bv}}=\Delta \phi$$(4)

B is a matrix of M × N. Each row of B corresponds to an interferogram and each column corresponds to a scene SAR image. If M ≥ N, the B can be obtained using the least squares method. Otherwise, matrix B is obtained by the singular value decomposition (SVD) method. Subsequently, the deformation rate at the acquisition time of each image can be calculated.

This study used a total of 67 ascending Sentinel-1A images to analyze the ground deformation characteristics of Xianyang city. Among them, the Sentinel-1A image acquired in July 2019 serves as the primary image for the dataset. The time and spatial baselines set for this study are 120 days and 150 m, respectively. Then, a total of 366 interferograms are obtained (Figure 6). For this research, the Goldstein filtering technique was applied to reduce random noise in the interferogram, and the phase unwrapping was carried out using the minimum cost flow approach. Finally, the deformation results in the radar coordinate system were transformed into the geographic coordinate system through geocoding.

Figure 6. Acquisition dates and interferometric pairs of satellite radar imagery. The green triangles represent Sentinel-1A acquisitions and the yellow triangle represents the primary image. Y-axis shows the length of perpendicular baselines.

3.2. InSAR results validation

Due to the lack of measured data that overlaps with the monitoring period of InSAR technology, InSAR monitoring results from Qu et al. (2022) from 2015 to 2019 were employed to validate the accuracy of the InSAR results. According to the incidence angle of the satellite, the average annual deformation rate monitoring results were transformed into the vertical direction. To address the inconsistency in spatial resolution between the two results, the corresponding InSAR results are obtained by calculating the annual deformation rate within a 100 m radius of the points of Qu et al. (2022).

Field investigations were conducted in selected typical surface deformation areas. The purpose of these investigations was to examine the correlation between the observed deformations and the actual ground conditions, thereby providing additional evidence for the accuracy of the deformation results.

3.3. Inferring gradient of land subsidence

In surface deformation or subsidence monitoring, it’s common to utilize indicators such as the subsidence rate or cumulative subsidence (Zhou, Lan, et al. 2022). These indicators help describe the velocity or the cumulative impacts of deformation over time. Despite their usefulness, it’s important to note that in engineering applications and disaster prevention, spatially uniform subsidence often results in less severe consequences compared to differential subsidence (Zhou et al. 2018; Novelo-Casanova et al. 2021). Therefore, it becomes essential to further study the differential subsidence. In this study, the concept of subsidence gradient is utilized to quantify the magnitude of differential subsidence. And the subsidence gradient (β) is defined as the ratio of the change of cumulative subsidence ΔU between two points to the distance l between them: (5) $$\beta =\frac{\Delta U}{l}$$(5)

Specially, l represents the original spacing between the identified coherent points, or in this context, the size of the grid employed for resampling. To ensure uniformity in the resolution of factors assessing susceptibility, l has been set to 60 m in this study. Then, we normalize the value of β to the range of 0 to 1. For clearer differentiation, the subsidence gradient is categorized into various levels using natural break classification.

Localized differential subsidence induces strain and bending of the surface layer, which ultimately leads to tension cracks and fractures. Considering that ground fissures predominantly occur in areas where subsidence patterns are spatially variable, the subsidence gradient is analyzed to map areas prone to ground fissures.

3.4. Ground fissure susceptibility mapping

Natural and human factors are the primary controlling factors for the occurrence and development of ground fissures. Among them, natural factors, such as geological characteristics, geomorphology and fault distribution have been widely recognized as the main disaster-causing factor (Choubin et al. 2019). Meanwhile, changes in groundwater caused by human activities are the key reason that triggers the rupture and expansion of ground fractures (Gambolati and Teatini 2015; Ye et al. 2018). However, it’s worth noting that there is still no consensus on which influencing factors should be chosen as the evaluation criteria for mapping the susceptibility of ground fissures. Considering the development characteristics, genetic mechanism and data availability of ground fissures in the Xianyang city, four types of auxiliary data are utilized to train the machine learning model for predicting the susceptibility of ground fissures. These data include geological, geomorphology, surface deformation and groundwater change (Table 1). The obtained data are outlined as follows.

Table 1. The SAR and ground fissure influencing factors data used in this study.

Type	Data	Source	Spatial resolution	Temporal coverage
Geological	Lithology	Shaanxi Institute of Geo-Environment Monitoring	Polygon	–
	Fault	Institute of Geology, China Earthquake Administration	Polygon	–
Geomorphology	Geomorphology	Shaanxi Institute of Geo-Environment Monitoring	Polygon	–
Surface deformation	Sentinel-1A Track 84	European Space Agency	5 m × 20 m	20150620–20221228
Groundwater change	Change of groundwater	Shaanxi Institute of Geo-Environment Monitoring	–	2017–2021

As a linear geological disaster, the spatial distribution of ground fissure disaster is significantly different from that of landslide and collapse (Zhu et al. 2018). According to the distribution of ground fissures in Xianyang city, more than 80% of the ground fissures are less than 500 m in length. To avoid the redundancy of samples, the disaster points are extracted on the ground fissure line at equal intervals (100 m). Subsequently, the inventory map was processed using random sampling. 70% of the data points were chosen as training data to develop the susceptibility model for ground fissures, while the remaining 30% served as test data for model validation (Hakim et al. 2023).

Factors influencing the occurrence and development of ground fissures need to be quantified and categorized before being used to evaluate their spatial relationship with the susceptibility distribution of ground fissures. The Frequency Ratio (FR) method is a commonly used approach for assessing the susceptibility of geological hazards and is widely applied for quantifying factors (Mohammady et al. 2019). The FR method evaluates the impact of a specific geological factor on the formation of geological hazards by comparing its frequency distribution in known hazard occurrence areas and non-occurrence areas. In this study, the FR values are determined by comparing the ratio of the ground fissure length within the region to the total ground fissure length, and the ratio of the area where the ground fissures occur to the total area. ANN is an information processing system modeled after the structure and function of neurons (Kulsoom et al. 2023). It is composed of numerous neurons, capable of processing and learning from vast amounts of input data. Each neuron receives input signals, processes them using an activation function, and produces an output signal, which is then passed on to the next layer. Through iterative training, selecting appropriate activation functions, and adjusting connection weights to minimize the root mean square error of the outputs, model accuracy can be improved (Duwal et al. 2023). In this study, four combinations of factors were input into a well-trained ANN model: four geological factors (ANN), four geological factors with land subsidence (ANN-velocity), four geological factors with deformation gradient (ANN-gradient) and four geological factors with deformation gradient and land subsidence (ANN-grad. + vel.) to generate the susceptibility map of ground fissures (Table 2).

Table 2. Four combinations of factors in ANN method.

	ANN	ANN-velocity	ANN-gradient	ANN-grad. + vel.
Geomorphology	√	√	√	√
Lithology	√	√	√	√
Distance to fault	√	√	√	√
Change of groundwater	√	√	√	√
Land subsidence velocity		√		√
Subsidence gradient			√	√

The susceptibility of ground fissures in the study area is determined based on prediction results and categorized into five levels. The higher the level, the greater the probability of new ground fissure occurrences (Mohebbi Tafreshi et al. 2020). The evaluation model for ground fissure susceptibility is validated using both the disaster point FR method and the receiver operating characteristic (ROC) curve. The closer the area under the curve (AUC) value is to 1, the more accurate the assessment of geological hazard susceptibility (Mehrnoor et al. 2023).

4. Results

4.1. Spatial and temporal characteristics of land subsidence

4.1.1. Spatial pattern of land subsidence

To track the urban land subsidence development patterns occurred across the Xianyang city, the average subsidence rate obtained by SBAS-InSAR from June 2015 to December 2022 projected along the line of sight (LOS) is shown in Figure 7. The subsidence distribution in the study area of Xianyang city is uneven. The study area exhibits severely subsidence in the northeast and southeast, relatively stable deformation in the central part, and ground rebound characteristics in the eastern part. Two main subsidence areas were observed: one is the Yunyang town and the other is the Luqiao town subsidence area (rectangle b in Figure 7). The most significant deformation occurred in Yunyang town, and it experienced an annual average subsidence rate ranging from −31 to 18 mm/year. In Luqiao town, the rate varied between −26.7 and 4.7 mm/year. As for the cumulative subsidence from 2015 to 2022, Yunyang town recorded a subsidence of −258 mm, whereas Luqiao town had a maximum subsidence of 139 mm. Another subsidence center exists in the Shi village, southern part of Jingyang county, with the maximum annual average subsidence rate being 25.5 mm/year and a maximum cumulative deformation of −200 mm over 8 years (see rectangle c in Figure 7). The western part of Weicheng district, the southern region of Qindu district, and the central-southern area of Xingping city exhibit more pronounced surface deformation, with annual average deformation rates of −18.2 mm/year, −13.9 mm/year, and −10.7 mm/year, respectively. Their corresponding cumulative deformation amounts is 114.3 mm, −93.6 and 79.4 mm, respectively. Conversely, Qian county experiences significant ground rebound, with the maximum annual average rebound rate being 10.3 mm/year and a maximum cumulative rebound amount of 93 mm over 8 years. As shown in Figure 7, the average ground deformation in other areas of Xianyang city is −10 mm/year to 10 mm/year, and overall remains stable.

Figure 7. (a) Mean velocity map of city subsidence in xianyang city. (Qd: Qindu district; Wc: Weicheng district). A negative value in the average velocity map signifies that the surface is in motion away from the sensor. The positions of the main deformation areas are outlined by two black dashed boxes, (b) Mean velocity map of subsidence center in Luqiao-Yunyang town, and (c) mean velocity map of subsidence center in Shi village-Chongwen town. The position of the triangles indicates the locations of the field survey points. The position of the square indicates the locations of the accuracy assessment points.

4.1.2. Temporal evolution of land subsidence

Utilizing InSAR-derived deformation results from 2015 to 2022, the study obtains the long-term trends of land subsidence in Xianyang city. Two prominent feature points (P1 and P2) with significant deformation were selected in each study area to create cumulative deformation time series graphs (Figure 8). For P1, the cumulative subsidence was approximately 126.1 mm, showing two distinct stages of deformation: rapid subsidence (−21.7 mm/year) from June 2015 to December 2020, followed by noticeable seasonal variations from December 2020 to December 2022. P2 experienced a cumulative subsidence of about 171.2 mm, reflecting a long-term deformation rate of −24 mm/year in Yunyang town, slightly higher than the surface deformation rate in Luqiao town (Figure 8). However, a decreasing trend in deformation rate was observed in Yunyang town starting from 2021.

Figure 8. Cumulative deformation in the subsidence centers of Xianyang city. The gray shading represents the decreasing trend observed in the data, while the purple shading indicates the increasing or rebound trend observed.

While the deformation magnitude of the two feature points was different, they exhibited highly similar subsidence trends over the past seven years. In addition, the land subsidence shows a clear seasonal deformation trend, which will be discussed further in Section 5.1.

4.2. InSAR accuracy assessment

4.2.1. Accuracy assessment

The InSAR results were validated by comparing them with those presented in Qu et al. (2022). The locations of the points used for this comparison are illustrated in Figure 7a. Given that the reference study spans from 2015 to 2019, the deformation rates for the same period were used to facilitate the validation. The average difference between the two sets of results is 1.9 mm/a, with a standard deviation of 1.6 mm/a (Figure 9). This indicates a relatively good agreement between the two datasets.

Figure 9. Comparison of the InSAR measurements with the result of Qu et al. (2022). The location of points 1–8 shows in Figure 7a.

4.2.2. Field accuracy assessment

Field investigations were carried out in typical regions with significant surface deformation to validate the current deformation patterns in the study area. The investigation revealed evident ground fissures caused by differential subsidence in the subsidence centers of Luqiao town and Yunyang town. The observed fissures mainly appeared on residential buildings with widths ranging from 0.5 to 2 cm (see Figure 10). Apart from the subsidence centers, notable ground fissures were also observed in two areas with abnormal ground deformation within the main urban area of Xianyang city. Numerous ground fissures were observed in the Taiqu parking lot and nearby shops in Bianfang village, Weicheng district, Xianyang city (see Figure 10c, d). The fissures predominantly aligned in a north-south orientation. And the longest fissure measured up to 65 m with a width of 0.1–1 cm. Along the road in Weicheng district, Xianyang city, a ground fissure with a width of 1–5 cm was observed, and the road surface also exhibited evident deformation and unevenness (Figure 10e, 10f). Given the close correspondence between the InSAR results and field investigations, it can be concluded that the SBAS-InSAR results in this study are reliable and accurate.

Figure 10. Buildings and road damage caused by uneven subsidence. The locations of the field survey points are shown in Figure 7, where triangles are used to represent these points.

4.3. Subsidence gradient

As described in Section 3.3, the subsidence gradient distribution map was constructed based on the deformation results (Figure 11). The main features of the gradient map were manually outlined to estimate potential ground fissure-prone areas associated with land subsidence. The larger spatial gradient is mainly concentrated within and at the edges of the subsidence funnel in the northern part of the study area.

Figure 11. (a) Overview of subsidence gradient, (b) detailed view of the subsidence gradient in subsidence centers. (c)–(e) On-site investigation of predicted ground fissure areas.

For assessing the correlation between subsidence gradient and ground fissures, the subsidence gradient was divided into five categories using the natural break method. The lengths and proportions of both existing and the predicted ground fissures within each category were determined (Table 3). The results indicate that approximately 81.11% of the predicted ground fissures were situated in areas of moderate to high gradient categories.

Table 3. Correlation of subsidence gradient class intervals and cumulative length of ground fissures in Xianyang city.

Category	Gradient	Existing ground fissure (km)	Proportion	Predicted ground fissures (km)	Proportion	Cumulative ground fissure (km)
Very high	>0.25	2.65	4.03	7.10	52.21	9.75
High	0.12-0.25	6.86	10.44	2.18	16.03	9.04
Intermediate	0.06-0.12	14.51	22.08	1.75	12.87	16.26
Low	0.02-0.06	20.68	31.47	1.31	9.63	21.99
Very low	<0.02	21.01	31.97	1.26	9.26	22.27
Total		65.71		13.6		79.31

4.4. Ground fissures susceptibility map

4.4.1. Factors influencing ground fissure susceptibility and data processing

Ground fissure susceptibility is affected by geological, geomorphological, surface deformation, and groundwater change characteristics. These four categories of data are employed for training the machine learning model in the prediction of ground fissure susceptibility. For this study, lithology and fault data were utilized to illustrate the lithology and fault development conditions in Xianyang city. The distance to the fault factor was derived using the Euclidean distance tool in ArcGIS. Geographical aspects, specifically the boundary zones between different landforms, have been observed to influence ground fissure development significantly. Land subsidence distribution map was created through InSAR technology, from which a gradient distribution map was subsequently extracted utilizing the gradient calculation function in ArcGIS. Groundwater change data from 2017 to 2021 for the Xianyang city was acquired from monitoring sources and interpolated using the inverse distance weighting (IDW) method.

All the gathered data were projected to the CGCS2000 108 E projection and resampled to a 60 × 60 m resolution, in alignment with the InSAR deformation map output (Figure 12).

Figure 12. Influential factors of ground fissures: (a) geomorphology, (b) lithology, (c) distance to fault, (d) groundwater change, (e) land subsidence and (f) subsidence gradient.

4.4.2. Ground fissures susceptibility map

The four ground fissures susceptibility maps derived from ANN are shown in the Figure 13. The results indicate that high susceptibility zones are mainly distributed near the Kouzhen-Guanshan fault and the Weihe fault. In the ANN model, the very high susceptibility zone accounted for approximately 0.04% of the total study area, while the high susceptibility zone covered around 9.64%. Compared to the ANN model, the ANN-velocity model has slightly larger area (4.5%) for very high susceptibility zone, whereas the high susceptibility zone encompassed about 10.55%. In the ANN-gradient model, the very high susceptibility zone occupied around 3.90% of the total study area, with the high susceptibility zone covering approximately 23.40%. In the ANN-grad. + vel. model, the very high susceptibility zone occupied around 4.03% of the total study area, with the high susceptibility zone covering approximately 21.54%. The above differences indicate that the ANN-gradient model and ANN-grad. + vel. model exhibit a greater proportion of very high and high susceptibility zones compared to the other two models.

Figure 13. Ground fissure susceptibility maps of the Xianyang city: (a) ANN, (b) ANN-velocity, (c) ANN-gradient and (d) ANN-grad. + vel.

5. Discussion

5.1. Nonlinear deformation behaviour and seasonal effects

The cumulative deformation map (Figure 8) reveals distinct seasonal subsidence and rebound patterns in both Yunyang town and Luqiao town subsidence centers. To effectively illustrate these seasonal deformation trends, the deformation results from the summer and winter stages were graphed based on the time series data (Figure 14a). Both Yunyang town and Luqiao town undergo accelerated subsidence from April to October each year, which subsequently transitions into a period of gradual rebound from March to April of the succeeding year. Prior to October 2018, the rebound deformation exhibited a relatively low magnitude, featuring a subsidence-to-rebound ratio of approximately 2:1. However, following October 2018, an observable increase in rebound deformation was noted, resulting in a subsidence-to-rebound ratio that approximates 1:1.

Figure 14. Map of seasonal deformation at subsidence centers in Xianyang city: (a) the deformation of the subsidence center at different stages, (b) the optical images in Luqiao town, (c) the optical images in Yunyang town and (d) the area of land subsidence centers, where 1 represent the time periods for 10/2015–04/2016, 2 represent the time periods for 04/2016–10/2016 and so on.

Analysis of Google historical imagery features since 2010 revealed that both areas were primarily used for agricultural cultivation, and there were no significant land use changes in the past decade (Figure 14b,c). Field investigations confirmed that Luqiao town and Yunyang town served as major vegetable cultivation areas in Shaanxi province, heavily reliant on groundwater for production and domestic use. According to local villagers, irrigation activities usually start around April each year and end in October. This is in agreement with what is shown by the deformation time series in Figure 14a.

Considering the spatial distribution of subsidence centers, the subsidence center in Luqiao town in 2022 remained relatively stable compared to 2015. And the deformation location did not change significantly. This consistency aligns with the findings derived from optical remote sensing imagery (Figure 14b). However, in Yunyang town, significant changes have been observed in the past eight year. These changes were marked by a gradual expansion of the ground deformation area (Figure 14c,d). Additionally, the subsidence center shifted from the southern to the northern sections, establishing the present pattern by 2019. Due to the lack of local groundwater monitoring data, this part will be studied in subsequent work.

5.2. Relationship between surface faults associated to land subsidence and the subsidence gradient

By plotting profiles of section A-A’ and B-B’ across the typical subsidence centers and faults, the relationship between subsidence gradient and faults was analyzed (Figure 15). The positions of the profile lines are shown in Figure 7. The A-A’ profile crosses the Luqiao-Yunyang subsidence center, traversing the Yidian-Qianxian-Meiyuan fault and the Fufeng-Sanyuan-Pucheng fault from north to south. From Figure 15a, it can be observed that the maximum differential subsidence at the Fufeng-Sanyuan-Pucheng fault zone is 107.2 mm, with a maximum gradient of 123.8 mm/km. Additionally, Figure 15a identifies a new crack with a differential subsidence of 82.3 mm and a maximum gradient of 78.4 mm/km. Field investigations confirmed that this newly discovered crack is located in Nanjie village, Luqiao town, and several housing cracks were observed (Figure 10i,j).

Figure 15. Distinctive land subsidence attributes and the gradient of subsidence in the horizontal direction on either side of the fault or ground fissure: (a) section A-A’ and (b) section B-B’. I represent the Kouzhen-Guanshan fault, III represent the Fufeng-Sanyuan-Pucheng fault, Ⅴ represent the Weihe fault, II and IV represent the ground fissures. The positions of the profile lines are shown in Figure 7.

The B-B’ profile runs from north to south, intersecting the Kouzhen-Guanshan fault, the Fufeng-Sanyuan-Pucheng fault, and the Weihe fault. Within the Kouzhen-Guanshan fault zone, the highest differential deformation reaches 75.9 mm, accompanied by a maximum gradient of 84.1 mm/km. Similarly, in the Fufeng-Sanyuan-Pucheng fault zone, the differential subsidence records 82.6 mm, with a gradient of 71.6 mm/km. Regarding the Weihe fault zone, the maximum differential subsidence measures 37.8 mm, featuring a maximum gradient of 40.6 mm/km. Furthermore, Figure 15b identifies two distinct differential subsidence induced ground fissures. The first ground fissure (II) is located in Beishuifeng village, Yunyang town, with a detected differential subsidence of 86.8 mm and a maximum gradient of 59.7 mm/km. This fissure has caused house cracks ranging from approximately 0.5 to 1 cm in width (Figure 10g,h). Near Xianyang Airport, numerous N-S oriented ground fissures (IV) were found in a parking lot, with a detected differential subsidence of 37.8 mm and a maximum gradient of 40.6 mm/km (Figure 10c,d).

The above analysis reveals that in areas traversed by active faults, the land subsidence profiles exhibit clear transitions or discontinuities, indicating significant uneven subsidence on both sides of the faults. The presence of faults has a pronounced limiting effect on land subsidence. Simultaneously, areas with higher subsidence gradient also exhibit apparent crack formation, indicating that differential subsidence is the primary cause of ground deformation (Figueroa-Miranda et al. 2020; Lei et al. 2021; Zhu et al. 2022).

5.3. Susceptibility model performance evaluation

Assessing model accuracy is a crucial aspect of susceptibility evaluation (Yao et al. 2023). And, the susceptibility results were validated using the disaster point FR method and ROC curves (Rahmati, Golkarian, et al. 2019). First, we examined the distribution of ground fissure disaster points in each susceptibility zone and calculated the FR value (Figure 16). From the FR results, it can be observed that when using four factors for evaluation, 76.4% of the ground fissures in the study area fell into the high and very high susceptibility zones. By introducing the land subsidence factor, the proportion of the total ground fissure length in the high and very high zones rose to 77.5%. When the assessment was adjusted to consider the gradient factor instead of land subsidence factor, it further increased to 82.26%. The incorporation of both the gradient and land subsidence factors resulted in an increase in the proportion of total ground fissure length in the high and very high zones to 83.54%.

Figure 16. Susceptibility zones distribution for the four susceptibility models. For clarity, the area in each category is displayed as “area/100”, and the FR value in each category is displayed as “FR value* 5.”

After incorporating the subsidence gradient factor, the very high susceptibility zone covered 3.90% of the total area but accounted for 50.98% of the ground fissures, yielding a high FR value of 13.12. Compared to the evaluation results without deformation factors and with deformation rate factors, the FR value increased by 4.04 and 1.77, respectively. When considering both gradient and ground deformation factors, 53.51% of ground fissures fall within the zone of very high susceptibility, with a frequency ratio of 12.49. Consequently, the application of the gradient factor notably raises the FR value of areas classified as very high and high susceptibility. Meanwhile, this approach ensures that areas experiencing high-gradient deformations are accurately categorized into high-susceptibility zones, thereby appropriately increasing their recognized risk level.

To evaluate the precision of the susceptibility maps for ground fissures generated using ANN, we employed ROC and AUC, with results displayed in Figure 17. The ROC-AUC values indicate that all four models demonstrate high accuracy in predicting ground fissure susceptibility. The ANN-gradient model and ANN-grad. + vel model have the highest AUC value (0.88), followed by the ANN-velocity (0.86), with the ANN model having the lowest AUC value (0.85). In terms of susceptibility zoning results, the introduction of deformation characteristic factors, especially the subsidence gradient, increased the proportion of high and very high susceptibility zones, enhancing the identification ability for disasters. The newly added high and very high-risk zones mainly corresponded to areas with high gradient, which were not identified as high susceptibility areas based on the previous factor combinations.

Figure 17. The AUC results of four susceptibility evaluation models.

Therefore, based on ROC curve analysis and the results of newly identified ground fissures, the introduction of deformation factors enhanced the model’s ability to capture details and improved its performance, enabling better prediction of ground fissure susceptibility and yielding more reasonable results.

The evaluation results indicated that regions most susceptible to ground fissures predominantly lie in the northwest and central parts of the study zone, especially within the alluvial plains area (Rahmati, Falah, et al. 2019). Subsidence gradient and faulting were crucial in the development of ground fissures. Additionally, the interface between the alluvial plains and geomorphic units notably influenced the areas of high susceptibility to ground fissures. The use of subsidence gradient maps obtained through InSAR can provide reliable support for ground fissure mapping and hazard assessment (Awasthi et al. 2022; Hu et al. 2022). However, it’s vital to acknowledge that a pronounced gradient of vertical subsidence doesn’t always equate to noticeable surface cracking. A comprehensive assessment should also consider the geological conditions and land use at the site. Furthermore, due to data limitation in the study, only one-dimensional deformation information has been extracted for the identification and susceptibility assessment of ground fissures. Considering the significant role of horizontal deformation in fissure development and characterization, adopting a multi-track combined monitoring approach is advisable for future studies (Zhang, Cheng, et al. 2023). Additionally, integrating data from other sources, such as GPS and field surveys, could enhance the accuracy and reliability of ground fissure detection and characterization.

6. Conclusions

Utilizing SBAS-InSAR technology, 67 Sentinel-1A images captured between 2015 and 2022 were processed. The deformation rate map and time series map of the plain area of Xianyang city were obtained. The susceptibility of ground fissure disasters in the Xianyang city was evaluated using four ANN methods, taking into account various geological data sources, land subsidence, and subsidence gradient data. The primary conclusions drawn from this study are outlined as follows:

1. Two main subsidence centers were identified in the north-east of Xianyang city, namely Luqiao town and Yunyang town subsidence center. The annual average subsidence rate in Yunyang town is -31 mm/year, while in Luqiao town, it is -26.7 mm/year. These results align well with the displacement data from previous studies (InSAR and leveling data) in Xianyang city.

2. The time series analysis of the subsidence centers in Yunyang town and Luqiao town reveals distinct seasonal subsidence and rebound patterns. The subsidence is particularly noticeable from April to October each year, followed by a gradual rebound until March to April of the following year.

3. The subsidence gradient has been identified as a favorable method for detecting ground fissures. Areas with high gradient correspond well with regions where ground fissures are densely developed. Additionally, the surface subsidence gradient is an important influencing factor in the assessment of ground fissure susceptibility.

4. This study proposed a new approach that integrated the land subsidence gradient factors for the susceptibility mapping of ground fissures using the artificial neural network algorithm. The susceptibility zonation results from Xianyang city validate the precision and viability of this technique.

Disclosure statement

No potential conflict of interest was reported by the authors.

Data availability statement

The data that support the findings of this study are available from the corresponding author, up-on reasonable request.

Additional information

Funding

This research was financially supported by the Opening Fund of Key Laboratory of Earth Fissures Geological Disaster, Ministry of Natural Resources (EFGD2021-05-01), the National Natural Science Foundation of China (Grant No. 41877250), the Young Talent Fund of Xi’an Association for Science and Technology (Grant No. 959202313094), the Megacity (Xi’an) Slow-onset Geohazard Early Warning Project (Grant No. 202307), the Opening Fund of State Key Laboratory of Geohazard Prevention and Geoenvironment Protection (Chengdu University of Technology) (Grant No. SKLGP2022K006) and the Fundamental Research Funds for the Central Universities, CHD (Grant Nos. 300102263501 and 300102263401).