Karst landslides detection and monitoring with multiple SAR data and multi-dimensional SBAS technique in Shuicheng, Guizhou, China

Shuicheng District is a karst mountain area, located in Guizhou Province, China. Its fragile stratum and frequent underground mining activities makes it prone to landslides. Owning to its wide coverage and frequent revisit, the InSAR technology has advantages in potential landslide identification and deformation monitor. However, affected by dense vegetation and atmospheric delay, it is much difficult to get sufficient effective targets to derive the deformation in this area. Besides, deformation derived from single orbit SAR data can result in the missing identification of some potential landslides and the misinterpreting of the real kinematics process of landslides. In this study, the multi-source SAR data, atmospheric error correction by quadratic tree image segmentation method, and phase-stacking method were selected to derive the surface deformation of this area. Besides, DS-InSAR and MSBAS method were combined to derive the deformation of Pingdi landslide. First, the potential landslides in this area were identified, surface deformation result, optical remote sensing images and geomorphological features were jointly considered. Then, the landslide distribution characteristics was analyzed in terms of slope, elevation and stratum. After that, the deformation along the LOS direction was acquired using the DS-InSAR method. The MSBAS method was used to retrieve the two-dimensional deformation of Pingdi landslide. Finally, the comprehensive analysis of triggering factors and failure process were conducted according to the spatial-temporal deformation characteristics and field investigation. The results indicated that landslides in Shuicheng district were mostly located in the junction of T1 and P3 stratum and mining related. Mining activity was the main cause of the Pingdi landslide deformation, the precipitation was the driving factor of the landslide instability. The research provides an insight into the explore the unstable slope distribution characteristic and the failure process of the landslides.


Introduction
China has the largest karst area in the world, accounting for nearly 1/3 of the total area of karst in China (Wu, 1998). Because of the complex geomorphologic environment, fragile stratum environment and frequent underground mining activities, largescale landslides and severe casualties occurred frequently in karst mountain areas (Huang et al., 2011;Zhao et al., 2022;Zhong et al., 2022). For example, on 3 December 2004, a catastrophic landslide caused by mining activity occurred at the front edge of the mountain in Zongling town, China, causing 39 deaths, thirteen people injured, and five people missing (Wu et al., 2006;Wang et al., 2020); on 2 August 2017, a combination of rainfall and mining activity led to a catastrophic landslide in the Pusa Village, Nayong Town, China, resulting in nine missing and 26 deaths (Chen et al., 2020a;Chen et al., 2020b); on 8 May 2022 the Baiyan landslide triggered by underground-mining in Zhijin Town, China, causing 3 deaths (CCTV news, 2022). These landslides were all located in highaltitude areas with steep slopes, caused by mining activities, difficult to detect and forecast in advance and mostly failed suddenly with complicated formation mechanisms (Wu et al., 2006;Wang et al., 2020;Chen et al., 2020a;Chen et al., 2020b;CCTV news, 2022). Although several landslides have been managed, those undetected potential landslides are always a threat. Therefore, it is of great significance to accurately identify potential landslides and analyze the deformation characteristics induced by underground mining activities and reveal the formation mechanism of landslides for the prevention and management of geological disasters.
Previous research of karst landslides mainly focused on the occurred landslide (Zheng et al., 2016;Lin et al., 2018;Kong et al., 2020), and the research methodology was mainly field investigation, geological exploration, and physical modeling (Wu et al., 2006;Gao et al., 2008). Few focused on the large-scale potential landslide identification, distribution characteristics and deformation process of landslide, so the distribution characteristics and kinematic process of potentially unstable slopes is still unknown. Since the slope deformation occur before the landslides (McKean and Josh, 2004), surface deformation monitor can be used for potential landslide identification and provides continuous view of deformation process. To this end, the distribution characteristics, spatial extension and temporal evolution of potential landslides can be explored through the deformation of the slopes.
Interferometric synthetic aperture radar (InSAR) technology has widely used in quick identify (Xu et al., 2014;Fobert et al., 2021) and monitor landslide deformation (Achache et al., 1996;Zhao et al., 2012;Kang et al., 2017;Liu et al., 2021), because of its wide coverage, high precision and frequent revisit. In last 2 decades, lots of advanced InSAR technologies have been proposed to increase the precision and accuracy of surface deformation results. The methods used in karst landslide deformation monitor can be divided into two categories according to the type of coherent targets, persistent scatterers (PS) (Ferretti et al., 2000;Ferretti et al., 2000) and distributed scatterers (DS). The PS can be identified from long temporal series of interferometric SAR images even with baselines larger than the so-called critical baseline (Ferretti et al., 2000). The commonly used advanced PS-InSAR technologies include the Interferometric point target analysis (IPTA) (Werner et al., 2003) and the Stanford method for PS (Hooper et al., 2004). The earliest and simplest use of DS point targets for surface deformation monitoring was the small baseline subset method (Berardino et al., 2002), which weakens the effects of decorrelation by selecting the interferograms with short spatial and temporal baseline. After that, many advanced SBAS-InSAR technologies have been proposed, such as new small baseline subset (NSBAS) (López-Quiroz et al., 2009) and Multidimensional small baseline subsets (MSBAS) (Samsonov and d'Oreye, 2017). To address the issue of high accuracy monitoring of surface deformation in areas with heavy vegetation, the SqueeSAR technology was proposed in 2011 (Ferretti et al., 2011). The basic principle of SqueeSAR technology is that to increase the density of observation targets, PS targets and DS targets are combined to derive surface deformation. In recent years, many scholars have extended and improved the selection of homogeneous points and phase enhancement of DS targets based on SqueeSAR technology. They proposed some improved algorithms, such as the Eigenvalue decomposition based Maximum-likelihood-estimator of Interferometric phase (EMI) (Ansari et al., 2018) and the Fast statistically homogeneous pixel selection (FaSHPS) (Jiang et al., 2015). These advanced DS-InSAR technologies can significantly increase the density of monitoring targets in landslide deformation monitor, especially in lush mountain areas (Chen et al., 2022a).
Atmospheric error is the most important source of InSAR deformation monitoring Wang et al., 2022). The phase-stacking InSAR (Sandwell and Price, 1998) can weaken the atmospheric delay phase with random variation characteristics. And for the complex topography of the Karst mountains, the vertical stratification of atmospheric delay is very severe, which usually spatially linear or power-law related to topography (Cavalié et al., 2007;Bekaert et al., 2015) and can be expressed as a function of elevation differences. In the existing methods, a single linear model cannot provide an accurate estimate of the atmosphere in the complex topography area, but quadratic tree image segmentation method can effectively correct it very well (Kang et al., 2021;Jiang et al., 2022).
Shuicheng District, Guizhou, China, a typical karst mountain area and one of the largest coal mining areas in southwest China (Yang and Tang, 2014), was taken as an experimental area. In this study, we designed a technologic route to identify the potential unstable slopes based on Multi-source SAR datasets, which can better overcome the shadow and layover problems intrinsic in single track SAR data, and improve the reliability and accuracy of the unstable slope (Grebby et al., 2021). According to the observation conditions and InSAR error characteristics, we corrected the atmospheric error by quadratic tree image segmentation method, and the phase-stacking method to derive the surface deformation. Then, ground surface deformation information and optical remote sensing imagery were combined to map active landslides. The surface cracks, debris and slopes, and mine locations can be extracted from high resolution optical remote sensing. Whilst deformation derived with InSAR technique can quickly identify the active landslides (Fan et al., 2019;Li et al., 2019).
In addition, an improved DS InSAR and MSBAS methods were conducted to derive the two-dimensional deformation time series of the Pingdi landslide, a typical "upper-hard-and-lower-soft" stratum landslide (Chen et al., 2022b). Furthermore, the effects of trigger  Frontiers in Earth Science frontiersin.org 03 factors of slope instability were explored. Finally, the failure process of Pingdi landslide was analyzed by combining the lithology, geomorphological features and the landslide deformation.

Study area
The exposed karst mountain area  in Shuicheng District is located in Yunnan-Guizhou Plateau, Liupanshui City, Guizhou Province, with a total area of 3,054.92 km 2 ( Figure 1). The topography is high in the northwest and low in the southeast, with an elevation range of 500-2,900 m. This area belongs to humid subtropical monsoon climate, and is characterized by abundant precipitation and lush vegetation. The annual precipitation is from 1,100 to 1,300 mm, with concentration from May to October, accounting for more than 80% of the annual precipitation.
The study area has rich mineral resources, with shallow coal deposits and frequent mining activities (Zhu et al., 2022). Figure 2 is the geological structure of the study area. The Permian, Triassic and Carboniferous strata are well developed and widely distributed, while the Devonian and Jurassic strata are sporadically distributed in blocks. The Permian strata are composed of medium-thick limestone, argillaceous siltstone and coals. Among them the Upper Permian Longtan Formation (P 3 l) is the primary coal-bearing strata in southwest China (Xiong et al., 2007). The Triassic strata are composed of calcareous dolomite, medium-thick limestone, purplish-red sandstone, and mudstone.
Folded deformation and fault structures widespread in this area, and the majority of the rock formations are hard and soft rock layers overlap each other. Weathering on the rocks is severe. Under different weathering conditions, the upper part of the strata forms high steep slopes or cliffs, and the lower part forms gentle slopes. As a result, the area has many high steep slopes and cliffs of hard rocks with soft bases, which constitutes the upper-hard-andlower-soft stratigraphic structure (Chen et al., 2016). According to such geological factors, coal fields are typically found at the bottom of slopes and coal seams are extracted from within the slopes (Yao, 2020), resulting in thousands of mining-induced landslides in China's karst mountain areas (Zhu et al., 2022).
The location of the Pingdi landslide is shown in Figure 2A, which is the typical high steep slope landform. The geological map of the Pingdi landslide is shown in Figure 2C. And the geological setting along profile AA' in Figure 2B is shown in Figure 3, where we can see that the exposed strata from top to bottom are the Lower Triassic Yongning Formation (T 1 yn), the Lower Triassic Feixianguan Formation (T 1 f), and the Upper Permian Longtan Formation (P 3 l). The T 1 yn and T 1 f are exposed strata of unstable slope, composed of hard limestone and medium hardness siltstone. The lower layer P 3 l is the coal-bearing strata, composed of soft argillaceous siltstone, which forms atypical upper-hard-and-lowersoft stratigraphic structure.
With the continuous underground-mining in recent years, the hazardous areas are further expanded. It is necessary to investigate the unstable slopes and to derive its failure process with the influence of underground mining.

Datasets
In order to detect the distribution of potential landslides in Shuicheng District and to monitor the deformation time series of the Pingdi landslide, a total of 299 SAR images were collected in this study, including three SAR datasets from Sentinel-1A and ALOS/ PALSAR-2 satellites. The spatial coverage of the SAR datasets is shown in Figure 1, and the basic parameters of the SAR images are summarized in Table 1.
To mitigate the effects of temporal and spatial decorrelation, the Small Baseline Subset (SBAS) strategy (Berardino et al., 2002) was

FIGURE 3
Geological settings along profile AA' of the Pingdi landslide, the location of AA' is marked in Figure 2B.

Frontiers in Earth Science
frontiersin.org used to generate interferometric pairs. After interferogram filtering (Goldstein and Werner. 1998), phase unwrapping (Costantini, 1997), atmospheric error correction by quadratic tree image segmentation (Kang et al., 2021), and DEM error correction (Liu et al., 2018), we selected high-quality interferograms for large-area unstable slope identification and small-scale monitoring of the Pingdi landslide.
For the large-area unstable slope identification, the SRTM DEM with a resolution of 30 m was employed to remove the topographic phase. Then to implement the precise monitoring of the Pingdi landslide, the ALOS digital surface model (DSM) with a resolution of 12.5 m was applied. The interferograms were multi-looked using factors of 2 × 2 (range × azimuth) for ALOS/PALSAR-2 images and 4 × 1 (range × azimuth) for Sentinel-1A images. And Google earth images were employed to cross-validate the suspected landslides identified by the InSAR deformation map. To explore the correlation between surface deformation and precipitation, the daily precipitation data were collected from the website of global precipitation measurement (GPM).

Methods
The technical route to investigate the potential instable slopes within the study area and to derive the spatiotemporal evolution of the Pingdi landslide is shown in Figure 4.
As the atmospheric errors have significant impact on deformation results in the study area, the atmospheric error correction by quadratic tree image segmentation (Kang et al., 2021;Jiang et al., 2022) was used to weakening the effect of the atmosphere. The technical route to identify the unstable slopes is organized as the following steps: calculation of annual deformation rate by phase-stacking InSAR method, quick identification of landslide and type by optical remote sensing and

FIGURE 4
Flowchart for the active landslide identification and deformation analysis.
Frontiers in Earth Science frontiersin.org geomorphological information, summary of the distribution characteristics of landslides. The DS-InSAR method was used to suppress the decorrelation effects caused by dense vegetation and increase the density of point targets, and MSBAS was used to calculate the two-dimensional deformation results. The technical route to derive the spatiotemporal evolution of Pingdi landslides is organized as the following steps: calculation of deformation results by DS-InSAR and MSBAS methods, analysis of the triggering factors of slope instability and simulation the failure process.

The atmospheric error correction by quadratic tree image segmentation
Interferograms are strongly influenced by the topographyrelated tropospheric delay in greatly undulate terrain mountain areas, where landslides occur frequently. Segmenting the interferograms based on the height difference of DEM, and estimating the vertically stratified atmosphere can suppress the influence of tropospheric atmosphere on the deformation results (Kang et al., 2021). The quadratic polynomial fit model is as follow: μ s x, y a 0 + a 1 x + a 2 y + a 3 xy + a 4 x 2 + a 5 y 2 + a 6 h s x, y (1) where, μ s is the vertical stratified component of the tropospheric delay and the long wavelength of the atmospheric delay, h s is the elevation, a i is the fitting parameters. The tropospheric delay is separated from the interferogram according to Eq. 1. In the phasebased estimation of the tropospheric delay, the low coherence pixels are masked to weaken the decorrelation noise. Then, reciprocate the initial deformation rate (Stacking-InSAR), to mask the points with large deformation rates. Finally, the atmospheric error is refitted using Eq. 1 and separated from the interference phase.

Optimization deformation based on distributed scatterers
In order to improve the efficiency of calculation and to focus on the landslide deformation process, after the identification of landslides, we cropped the SAR images. To reduce vegetationinduced interference decorrelation, the covariance matrix is estimated robustly with homogeneous points to enhance the phase of interferograms. The SAR data satisfy the standard Gaussian distribution (Jiang et al., 2018), of which the complex coherence matrix Τ can be expressed as: T 1 γ 1,2 e jφ 1,2 ... γ 1,2 e jφ 1,N γ 2,1 e jφ 2,1 1 ... γ 2,N e jφ 2,N ... γ N,1 e jφ N,1 γ N,2 e jφ N,2 ... 1 where, e j,φ i,j is the interference phase of the i image and the j image; γ i,j is the coherence of interferogram.
Since the phase of each distributed target is composed of the backscattered signals of multiple ground objects, the different scattered signals can be separated by eigenvalue decomposition of the covariance matrix to achieve the optimized phase corresponding to the largest eigenvalue as follow: where, λ 1 represents the largest eigenvalue of covariance matrix for the generic pixel; u 1 is the corresponding eigenvector, related to the dominant scatterer within a given pixel. The correction effect of this method was shown in Supplementary Appendix B.

MSBAS InSAR
The deformation along the LOS direction is the vector sum of the deformation in the vertical, east-west and north-south directions (Fialko et al., 2001;Wright et al., 2004) as follows: Due to the limitation of the near-polar orbit of SAR satellites, the projection coefficient in the north-south direction is much smaller than those in the east-west and vertical directions. Therefore, MSBAS (Samsonov and d'Oreye, 2017) is used to calculate the east-west and vertical surface deformation as follow: −cos θ sin φA cos φA λL 0 where, θ and φ are the azimuth and incidence angle, respectively; A is the time interval of acquired SAR images; λ is the regularization factor and L is the Tikhonov regularization matrix; V E and V U are east-west and vertical surface deformation rates;φ is the LOS deformation phase.

Landslides identification in Shuicheng District
Landslide investigation allows us to quickly understand the current status of landslide distributions in an area (Hussain et al., 2021) and guides us to select typical landslides to analyze its deformation process. A total of 42 landslide clusters were identified in this area. Comparing the identified results with the historical landslide sites in the area, we found 15 landslide clusters which were known as historical landslide hazard sites and the other 27 landslide clusters were unknown before ( Figure 5). Most of them were related to mining activities, which were consistent with those results by other researchers in this area (Zhu et al., 2022). The basic information of landslides was shown in Table 2.
The different number of landslides detected from different SAR datasets are mainly because of these reasons: Since the SAR images were acquired in the side looking acquisition mode typical of SAR sensors, the steep topography of karst mountains inevitably results layover and shadowing geometrical distortions for the InSAR observation (Hanssen, 2001), and different incidence angles form different invisible areas. The slope movement perpendicular to the LOS cannot be monitored, and the descending and ascending geometries are favorable, respectively, for west and east facing slopes (Wasowski and Bovenga, 2014), which would result in large differences in the deformation results obtained from the two data . And compared to the C-band Sentinel-1A, the L-band ALOS/PALSAR-2 has better penetration of vegetation and can detect more covert landslides (Chen et al., 2022a). Therefore, both descending and ascending acquisitions are Frontiers in Earth Science frontiersin.org available the limitations related to topography (layover, shadowing) can be significantly reduced (Grebby et al., 2021), and can increase the number of identifiable landslides.

Landslides distribution characteristics of Shuicheng District
The 30 m resolution DEM is employed to calculate the slope and the elevation of the identified active unstable slopes in this area. The percentage of landslide number in each section is defined as Landslide Number Percentage (LNP), while the percentage of the total natural area of landslide in each section is defined as Natural Area Percentage (NAP) (Yao et al., 2017).
Landslides with higher slopes have more potential energy and faster sliding speeds (Guo et al., 2008). According to Figure 6A, there was no landslide on slope below 10°slope angle, and the landslides were mainly distributed in the range of 20°-35°. In this range, the number of landslides increases with the slope. In the range of 30°-35°, the LNP and NAP are 26.19% and 32.72%, respectively.
Different elevation ranges have different vegetation types and vegetation cover. The higher the elevation is, the more severe the weathering will become, resulting in a high correlation between landslides and elevation (Wang et al., 2017). Landslides in this area were mainly developed at the elevation of 1,600-1800m and the LNP and the NAP in this range were 26.19% and 30.88%, respectively. There was no landslide below 1000 m and above 2200 m ( Figure 6B).
The lithology of the strata is the material basis for the formation of landslides. Different lithology has different mechanical properties, and different combinations of lithology and slope structure have different stability (Dai and Deng, 2020). Landslides in this area were well developed in the Lower Triassic Formation (T 1 ) and Upper Permian (P 3 ). They were mostly distributed at the junction of T 1 and P 3. The LNP and NAP were 83.34% and 83.22%, respectively ( Figure 6C).  According to the above statistics, many landslides in this area were related to geological and environmental conditions, such as topography and geomorphology, stratigraphic lithology, and geological formations, which are intrinsic to the development of landslides. They were mainly distributed in the junction of T 1 and P 3 , with slope 20°-35°, and elevation between 1,600 and 1800 m. *Notes: ALOS stands for ALOS/PALSAR-2 images; S1A and S1D stand for ascending and descending Sentinel-1A SAR images, respectively.

FIGURE 6
Statistics of geomorphological elements of mining-induced active landslides in the study area; Distribution ratio of (A) slope; (B) elevation; (C) stratum.
Frontiers in Earth Science frontiersin.org

Spatiotemporal deformation of the pingdi landslide
The Pingdi landslide is located on the west bank of the Beipan River, marked as Landslide No. 38 in Figure 5D. Figure 7 shows the historical optical remote sensing images acquired on Dec.21, 2015, Apr.6, 2017, Nov. 14, 2020and Nov. 11, 2020. The cracks could be seen on the slope surface ( Figures 7E,F). The white circle indicates the location of the coal factory. As shown in Figure 7D, the slope in the red rectangular area changed obviously, and a large amount of rock debris accumulated at the position shown by the purple arc. With the continuous mining activities, the fragile rock mass will lose its balance, forming rock debris flow , causing blockage of rivers.
The LOS deformation of the three orbits were calculated independently based on the DS-InSAR method. Due to the difference of ALOS-2 and Sentinel-1A datasets in quantity and the time of access, instead of calculating the three-dimensional deformation, we calculated the two-dimensional deformation using Sentinel-1A datasets based on MSBAS. The deformation results were sufficient for the subsequent analysis, which could provide valuable information for exploring the dynamic process, analyzing the driving factors and the failure modes.
The two-dimensional deformation rate of the Pingdi landslide was shown in Figures 8C,D, respectively. The red indicated the westward displacement ( Figure 8C) or downward displacement ( Figure 8D), while the blue indicated the opposite displacement, and the green was corresponding to the stable surface. Figures 8C,D indicated that the Frontiers in Earth Science frontiersin.org greatest vertical and horizontal annual deformation rate was at the lower and part of the slope, near point S3.
To get more two-dimensional time series information about slope deformation, we have selected four points located in the front edge of the slope, the upper, middle and lower part of the slope, near the profile AA' to reflect the deformation of the slope and crossreference with the results of the profile deformation. The location of four points was shown in Figure 8. The two-dimensional time serious deformation of Point S1, S2, S3 and S4 were shown in Figure 10. In vertical direction, the downward deformation of points S2 and S3 were large, while points S1 and S4 were small. In horizontal direction, the eastward deformation of S3 and S4 were large, while point S1 and S4 were small. Of these points, the cumulative deformation of S3 was the largest, the eastward deformation was 76 cm and the downward deformation was 124 cm. After January 2021, the deformation of S2 was dramatic, especially on vertical, the deformation rate from 24 January 2021 to 3 February2,021 was 2.84 mm/d, very unstable.
As shown in Figure 9, 19 ALOS/PALSAR-2 images were employed to derive the deformation time series in the LOS direction, where negative value (red) indicated the deformation away from the satellite and the positive value (blue) indicated the deformation towards the satellite. To the southeast part of the landslide, the deformation initially occurred at the middle and lower part of the slope, then the front edge of the slope started to deform, the most intensive deformation was at the middle and lower part of the slope, around Point S3, the cumulative maximum deformation of the slope was over 130 cm. From 16 April 2019 to 29 September 2019, the northeast part of the slope did not show significant deformation, but during 29 September 2019 to 10 May 2020, the northeast part of the landslide shown slow deformation and the deformation gradually expanded, connecting to the previous main deformation area. The detail of LOS cumulative deformation of the Pingdi landslide was shown in Figure 8B, and had obvious circle-like deformation on the slope. The similar circle-like deformation also reflects in the two-dimensional annual deformation rate, which will explain subsequently.

Distribution characteristics of landslides in shuicheng district
According to statistics above, most unstable slopes wew mainly distributed in the junction of T 1 and P 3 . The reasons are as follows: First, P 3 is the important coal bearing strata in South China, with shallow coal seams, which has the most coal mining activities; Second, as shown in Figure 2A, the T 1 stratum in the study area is superimposed on the P 3 stratum. And at the junction of T 1 and P 3 , the compose of T 1 in this area is mainly T 1 yn superimposed on the T 1 l, and the P 3 is mainly composed of P 3 l; Third, rocks in this area are heavily weathered. Under different weathering conditions, the upper part of the strata forms high steep slopes or cliffs, and the lower part forms gentle slopes. The T 1 yn and T 1 f are exposed strata of slope, composed of hard limestone and medium hardness siltstone. The lower layer P 3 l is the coal-bearing strata, composed  Figure 11, and the deformation time series at points S1, S2, S3 and S4 are displayed in Figure 10; (C) and (D) are the east-west and vertical annual deformation rate maps derived from ascending and descending Sentinel-1A images, respectively.

Frontiers in Earth Science
frontiersin.org of soft argillaceous siltstone, which forms an upper-hard-and-lowersoft stratigraphic structure. With frequent mining activity under the slope, mining goat is formed at the lower part of the slope, and the soft rock cannot support the upper slope, resulting in the failure of stress balance, then the slope will become unstable (Chen et al., 2022b). Therefore, these unstable slopes are mainly distributed at the junction of T 1 and P 3 in this area.

The driving factors of karst slope instability
The driving factors of karst slope instability usually can be attributed to natural factors and human activities (Bonacci and Juračić, 2010;Ford and Williams, 2013;Gutierrez et al., 2014). Some studies have shown from a geological point of view that in karst mountain area, precipitation and mining activity are the two important driving factors of slope instability (Xiao et al., 2019;Cui et al., 2022;Emami Meybodi et al., 2022;Zhong et al., 2022). In this section, we have analyzed these two important driving factors in terms of surface deformation.
First, we added the daily precipitation data to determine if there is a correlation between precipitation and the Pingdi landslide deformation. As shown in Figure 10, the vertical deformation rate of these four points was constant or accelerated at the beginning and middle of the rainy season. And there was a sudden increase in deformation rate after heavy precipitation. In times of low precipitation, the rate of deformation became slower. The east-west deformation rate of these four points showed similar phenomenon with vertical deformation rate. In general, changes in slope deformation rate do not occur simultaneously with precipitation, but lag behind it. And there is a significant positive correlation between the magnitude of deformation rate and the magnitude of precipitation (Lin and Guo, 2001).
Then the time series deformation calculated from ALOS-2 dataset along the profile AA' and BB' in Figure 8A was extracted. The profile AA' was plotted along the landslide movement direction approximately. The deformation in section Ⅲ (Figure 11A) was like a "subsidence funnel", and the circle-like deformation mentioned above usually relates to mining activities (Modeste et al., 2021). Figure 11A showed that the deformation initially occurred in the position of caving area (section Ⅲ), which also was the position of Frontiers in Earth Science frontiersin.org

FIGURE 10
Comparison between deformation time series and daily precipitation acquired from GPM. The locations of S1, S2, S3 and S4 are shown in Figure 8. (A, B, C) and (D) are comparison between two-dimensional deformation time series, two-dimensional deformation rate and daily precipitation for S1, S2, S3 and S4, respectively.

FIGURE 11
Line-of-sight deformation time series obtained from ALOS/PALSAR-2 data along profiles AA' and BB'. (A) LOS deformation time series along profile AA'; (B) LOS deformation time series along profile BB'. The grey area is the topography profile from ALOS DSM. The location of profiles is shown in Figure 8A.
Frontiers in Earth Science frontiersin.org 13 the maximum downward deformation mainly occurred. As the deformation intensified in section Ⅲ , the upper part of the slope section Ⅱ and the front edge of the slope section Ⅰ began to deform. This deformation conforms to the general performance of slope deformation caused by mining activities (Zheng et al., 2015). Similar deformation characteristic also occurred in Jianshanying landslide , Zongling landslide (Chen et al., 2022a) and Kaiyang landslide (Chen et al., 2022b), which were all mining related. The reasons for the deformation above are as follows: The downward deformation mainly occurs above the exposed strata T 1 yn and T 1 f, which is mainly composed of hard limestone and medium hardness siltstone, as the coal seams stored in stratum P 3 l which is mainly composed of soft argillaceous siltstone were mined out, the remaining coal pillar cannot bear the pressure of the roof gradually. Since the strength of the upper strata is stronger than that in the lower layer, repeated mining activities can induce the plastic extrusion from the upper strata, resulting in the tension fracture at the back edge of the slope (Chen et al., 2022b).
Through the analysis above, precipitation and mining activity are two important factors of slope instability. Mining activity is the main cause of the slope deformation, and precipitation is the driving factor of the slope instability. The precipitation can intense mininginduced slope deformation.

The failure process of the pingdi landslide
The spatiotemporal evolution feature of the Pingdi landslide was explored above, and the driving factors such as the precipitation and mining activities were also discussed according to the deformation characteristics. This section will reveal the failure process of the Pingdi landslide and analyze the role of these factors in the process.
As shown in Figure 11B, the maximum deformation was at the middle of the profile BB', where the initial deformation occurred, then the deformation decreases from the middle to both ends of it. And the maximum deformation was at the lower and middle of the profile AA' (Figure 11A), where the initial deformation occurred. Combined the deformation characteristic of the Pingdi landslide, we could affirm that the deformation initial occurred at the middle and lower part of the slope, near the point S3, where the maximum deformation occurred, and then the deformation decreased to periphery. This phenomenon is much consistent with the characteristics of retrogressive landslide, which shows that the deformation occurs at lower part of the slope first, causing the change of stressed structure of the upper part of the slope and later the deformation (Pascal et al., 2019).
The coal seams in the Longtan Formation (P 3 l) were gradually mined out, the roof gradually subsidence, and cracks were subsequently formed. With the exploration of coal mining, ground fissures develop and expand continuously, and the roof subsidence has increased. In addition, the deformation of the minedout areas could lead to uneven forces across the slope, which would cause the formation of the cracks near the top of the slope (Zheng et al., 2015), even cause the tension fracture at the back edge of the slope. Under the continuous action of gravity and precipitation, the cracks near the top of the slope accelerated cracking, and the dangerous rock mass gradually separated with the mountain. The ground fissures or faults induced by repeated mining activities

FIGURE 12
The failure process of the Pingdi landslide.

Frontiers in Earth Science
frontiersin.org would enhance the water conductivity of the cracks and stratum, which could further erode the unstable stratum and weaken the rock strength. The rock finally disconnects from the mountain and forms rock fall or avalanche. Through the analysis above, the failure process of Pingdi landslide can be summarized as Figure 12: (A) Due to precipitation and weathering, there were some cracks at the top of the slope prior to the mining activity. (B) As mining began, deformation and new cracks began to appear on the slope, and heavy precipitation would make cracks expand (Matsuura et al., 2008). (C) With the further expansion of the mining progress, the deformation area of the slope is further expanded and the force balance is broken up. The progressive failures, such as rupture, creeping, rockfall and sliding, was generated. (D) The slope is currently at a relatively stable deformation, continual mining activities and heavy precipitation may break the force equilibrium of the slope, causing a catastrophic landslide.

Conclusion
In this study, the distribution characteristic of landslides in Shuicheng district and the failure process of a typical landslide, Pingdi landslide, were explored. On the one hand, the atmospheric error correction by quadratic tree image segmentation methods was conducted to retrieve the surface deformation of Shuicheng district. On the other hand, the DS-InSAR and MSBAS methods were conducted to retrieve the surface deformation of Pingdi landslide. Atmospheric effect on karst mountains was sever, the quadratic tree image segmentation methods can effectively weaken atmospheric effects. The DS-InSAR method can improve the precision and the density of the deformation time series results, which can guarantee the InSAR application in karst mountain areas. The MSBAS method can be successfully applied to retrieve the vertical and east-west deformation time series.
The relationship between landslides distribution and stratum in Shuicheng district was revealed, that landslides were mostly located in the junction of T 1 and P 3 stratum. And the spatiotemporal characteristics of the surface deformation in Pingd landslide prove that mining activity was the main cause of the slope deformation, and the precipitation was the driving factor of the slope instability. The research provides an insight into the explore the unstable slope distribution characteristic and the failure process of the landslides.

Data availability statement
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

Author contributions
Conceptualization, RS, CZ, and BL; methodology, RS; software, RS; validation, RS; formal analysis, RS; investigation, RS and HC; resources, LC; data curation, RS; writing-original draft preparation, RS; writing-review and editing, RS and CZ; visualization, RS; supervision, CZ; project administration, BL and CZ; funding acquisition, CZ All authors contributed to discussions and manuscript preparation.