Failure process and three-dimensional motions of mining-induced Jianshanying landslide in China observed by optical, LiDAR and SAR datasets

ABSTRACT

The occurrence of collapses and landslides due to underground mining has its unique failure mechanism, especially in the Karst mountainous regions of China. Spaceborne and airborne remote sensing observations provide rapid and effective tools for assessing surface changes and monitoring surface deformation of such landslides. In this study, we take the Jianshanying landslide, a typical mining-induced and fast-deformed landslide, as an example, and reveal the failure mechanism of such landslide by investigating the historical surface displacement. First, the complete evolution of the landslide surface was investigated from its original state to the overall sliding. The data include the satellite and Unmanned Aerial Vehicle (UAV) optical images, UAV three-dimensional (3-D) real scene models, high-resolution Light Detection and Ranging (LiDAR) DEM, and field survey. The results show that the head region entered the high deformation stage after 2019, the maximum deformation rate was 12.3 m/yr. The landslide morphology was formed after the overall slide occurred in September 2020. Then, the pre-event 3-D surface deformation after the landslide entered the high deformation stage was recovered using Interferometric Synthetic Aperture Radar (InSAR), differential DEM, and SAR/optical offset-tracking techniques. The vertical deformation was recovered around −30 m from 2019 to 2020. In particular, we solved the problem of unequal accuracy of SAR and optical offset-tracking observations in 3-D deformation inversion by employing the Helmert variance component estimation method. The maximum deformation was 6 m and 3 m within 4 months in the NS and EW directions, respectively. Finally, we revealed the failure mechanism of the Jianshanying landslide based on the disparity of horizontal and vertical deformation. That is, underground mining causes a significant subsidence of the rear part of the landslide body, resulting in different stress changes in the rear and front parts of the landslide body, which eventually led to sliding of the front part of the slope along the free surface. This work investigates and monitors the typical underground mining-induced Jianshanying landslide by using multi-sensor remote sensing approaches to trace the pre-event surface motions and to reveal its failure mechanism.

KEYWORDS:

• Optical remote sensing
• SAR
• DEM
• mining-induced landslide
• deformation monitoring
• failure mechanism

1. Introduction

Due to the comprehensive effects of fragile geological environment, climate change, and anthropogenic activities, landslide disasters have occurred frequently in the mountainous areas of southwest China, one of the largest karst landforms in the world, and thousands of people have lost their lives due to landslide disasters in recent decades. For example, on 28 August 2017, a long-runout rock flow triggered by a ridge-top rockslide in Pusa Village, Zhangjiawan Town, Guizhou Province, China, buried residential areas and caused 26 deaths and 9 missing (Fan et al. 2019). On 23 July 2019, a landslide occurred in Pingdi Village, Jichang Town, Shuicheng County, Guizhou Province, China, resulted in 42 deaths and 9 missing (He et al. 2022). Therefore, revealing landslide mechanisms and failure mode is crucial for early warning of the disasters. In this study, we selected the landslide caused by underground mining as a typical study area, which is characterized by huge changes, high and rapid deformations (Chen et al. 2021; Wang et al. 2020). We applied multi-sensor remote sensing (RS) approaches to recover the long-historical surface evolution process in this area and explore the landslide failure mode under the influence of underground mining.

To date, RS observations has been successfully applied to individual landslide investigation and deformation monitoring, mainly involving in spaceborne optical and Synthetic Aperture Radar (SAR) images, airborne optical, SAR, and Light Detection and Ranging (LiDAR)-derived data, and others (Bhuyan et al. 2023; Casagli et al. 2017; Dai et al. 2022; He et al. 2019; Scaioni et al. 2014; Spinetti et al. 2019; Zhao and Lu 2018). For optical RS from all available platforms, qualitative investigation of surface change can be directly compared with multi-temporal optical images (Hölbling et al. 2017; Plank, Twele, and Martinis 2016; Qu et al. 2021; Rossi et al. 2018) or automated change detection methods using multi-spectral images (Kyriou and Nikolakopoulos 2020). Long-term or high surface deformation in north-south and east-west directions can be recovered by optical offset-tracking based on spaceborne and Unmanned Aerial Vehicle (UAV) optical images (Lacroix, Dehecq, and Taipe 2020; Teo et al. 2023; Xiong et al. 2020), which is limited to obtaining only horizontal deformation. For SAR RS, the multi-temporal Interferometric Synthetic Aperture Radar (InSAR) techniques are widely applied in landslides and it is suitable for slow-moving landslide deformation monitoring (Intrieri et al. 2018; Kang et al. 2021; Traglia et al. 2021; Xu et al. 2020). The observation obtained by this method is in SAR imaging geometry, and the deformation is one-dimensional (1-D) in terms of line-of-sight direction.

The SAR offset-tracking method is effective for rapid and large-magnitude deformation monitoring when the surface motion exceeds the maximum detectable displacement gradient of InSAR. It can obtain the two-dimensional (2-D) deformation in range and azimuth directions, even if the surface reaches meter-level movement in a short period of time. In addition to satellite SAR data from the same platform (Cai et al. 2022; Raucoules et al. 2013; Xu et al. 2020), historical and time-series deformation can also be retrieved based on this method using data from different satellite platforms (Liu et al. 2020), and airborne platforms (Hu et al. 2020). In addition, point clouds generated by UAV photogrammetry or LiDAR can be used to obtain multi-dimensional deformation of landslides (Bernard, Lague, and Steer 2021, Wang et al., 2022; Teo et al. 2023). In particular, the 3-D deformation can be retrieved by differencing the multi-period and high-resolution point clouds. High-resolution and multi-temporal DEMs can monitor deformation in the 1-D vertical direction, which can be reconstructed by InSAR, optical RS (Pieczonka et al. 2013; Zhang et al. 2021; Zhou, Li, and Li 2017), UAV, or LiDAR methods (Burns et al. 2010; Eker, Aydın, and Hübl 2018; Kakavas and Nikolakopoulos 2021; Shugar et al. 2021). Finally, the 3-D deformation of landslides can be inferred mathematically based on multi-track deformation observations under the same platform and external constraints (Hu et al. 2020; Liu et al. 2021; Samsonov et al. 2020). Furthermore, some studies combined SAR and optical deformation monitoring from different platforms to invert surface 3-D deformation of landslide (Li et al. 2020), but the unequal accuracy of the two types of observations was not fully considered. Consequently, relying on a single method cannot provide a comprehensive and accurate investigation and monitoring of mining-induced landslides with significant surface damage. Thus, it is necessary to integrate multiple RS methods to trace their historical evolution.

Accordingly, to explore the failure mechanism of landslides caused by underground mining, we chose the Jianshanying landslide as the study area because it has been affected by underground coal mining for decades with huge surface deformation. Multi-sensor RS approaches were applied to investigate and monitor this landslide, the pre-event surface motions of which were traced, and then revealed its failure mechanism. First, we conducted a qualitative study of the surface changes through satellite and UAV optical images, UAV 3-D real scene models, high-resolution LiDAR DEM and field survey to reconstruct the complete evolution of the landslide surface from its original state to small-scale deformation, then to high deformation, and finally to overall sliding. Then, we recover the 3-D surface deformation after the landslide entered the high deformation stage using InSAR, differential DEM, and SAR/optical offset-tracking techniques. In particular, we improve the fusion method by using Helmert variance component estimation to address the unequal accuracy of SAR and optical offset-tracking observations in 3-D deformation inversion. Finally, the failure mechanism of underground-mining landslide is revealed from surface deformation.

2. Study area

The Jianshanying landslide is located in Faer Town, Guizhou Province, China, and has been induced by underground mining since 2005 (Dong et al. 2021; Yu et al. 2020; Zhao et al. 2021). The landslide boundaries before and after sliding are mapped based on optical images as shown in Figure 1. The area belongs to the western plateau of Guizhou, and the topography fluctuates greatly. The free faces of Jianshanying's unstable slope are in the northeast and south regions. The highest elevation of this area is approximately 1526 m located in the southeast, the lowest elevation is about 949 m. The length of the cliff slope area is about 1300 m, while its width is about 150 m. The average annual temperature in this area is about 15.8°C, ranging from −6.3°C to 32°C. The average annual precipitation is approximately 1027.2 mm, ranging from 757.2 mm to 1192.2 mm, concentrated from May to September, which is one of the key factors of this landslide mechanism (Chen et al. 2021).

Figure 1. The full view of Jianshanying landslide on the UAV image, where the landslide boundaries before and after large-scale slide are shown in red and green polygons, respectively. The underground mining working lanes are indicated by black dashed polygons, and the inset is the location of this landslide in China.

Furthermore, there are 13 underground coal seams in this area, six of which have been mined from 2004 to 2018 (Li et al. 2022). The underground working faces are shown in Figure 1. Therefore, obvious deformation has occurred since 2014, and the landslide began to undergo large-scale creep after 2019, with significant surface damage and many cracks with a width of more than several tens of meters, while causing surface collapse. In September 2020, an overall sliding occurred, and thereafter the landslide boundary shown in Figure 1 was formed. There were mainly sparse trees and dense shrubs in this area before landslide occurred, but after the landslide, they were mostly sparse shrubs and bare ground.

3. Datasets and methods

3.1. Datasets

To investigate the pre-event surface changes of Jianshanying landslide and monitor its surface deformation, multi-temporal satellite/UAV optical images, UAV reconstructed 3-D models, UAV/LiDAR DEM, and satellite SAR images are involved, as shown in Table 1. First, the overall pre- and post-event surface changes from 24 March 2013 to 11 November 2020 are investigated through five Google earth optical images with a pixel spacing of 1 m. Then, 2 scenes UAV 3-D real scene models are used to detect changes of key blocks after large-scale creep, which are generated by Smart Capture 3-D processing software based on UAV oblique photographs, and LiDAR DEM is used to identify cracks and landslide morphology. The multi-dimensional deformations of this landslide are quantitatively monitored by multi-source data, in which the vertical material changes are obtained by differencing three scenes UAV/LiDAR DEMs with a pixel spacing of 1 m; the cumulative east-west (EW) and north-south (NS) deformations are calculated by 13 scenes optical images from Sentinel-2 Band 8a images, we choose the images with resolution of 20 m to reduce the noise in deformation caused by surface fragmentation. Besides, 19 scenes L-band ascending ALOS/PALSAR-2 acquisitions are used to recover the azimuth (AZI) and line-of-sight (LOS) deformations with single-look pixel spacing of 2.1 and 1.4 m. The metadata of the ALOS/PALSAR-2 acquisitions can be found in Table 2. This L-band datum has a unique advantage in penetrating vegetation and acquiring high coherence compared to C-band SAR data during InSAR processing.

Table 1. The multi-source RS data used in the research of Jianshanying landslide.

	No.	Type	Source	Resolution (m)	Time spanning	Scenes	Methods	Objectives
Pre-event	1	Optical	Google earth	0.1	2013.03.24–2018.12.24	5	Visual comparison	Surface changes
	2	3-D real scene model	UAV	0.1	2019.04.19–2020.08.02	2	Visual comparison	Key block Changes
	3	DEM	UAV	1	2019.04.19–2020.08.02	2	Difference	Vertical deformation
	4	DEM	LiDAR	1	2019.06	1	Difference	Vertical deformation
	5	Optical	Sentinel-2	20	2016.07.31–2020.03.22	13	Offset-tracking	EW/NS deformation
	6	SAR	ALOS/PALSAR-2	2.1×1.4	2017.04.16–2020.05.10	19	InSAR/Offset-tracking	Azimuth/LOS deformation
Post-event	7	Optical	Google earth	0.1	2020.11.11	1	Visual comparison	Surface changes

Table 2. The metadata of the ALOS/PALSAR-2 and Sentinel-2 images.

No.	Platform	Date (YYYY.MM.DD)	Mode	No.	Platform	Date (YYYY.MM.DD)	Mode
1	ALOS-2	2017.04.16	SM1	17	ALOS-2	2019.12.22	SM1
2	ALOS-2	2017.05.28	SM1	18	ALOS-2	2020.01.19	SM1
3	ALOS-2	2017.08.20	SM1	19	ALOS-2	2020.05.10	SM1
4	ALOS-2	2017.10.29	SM1	20	Sentinel-2	2016.07.31	-
5	ALOS-2	2017.11.26	SM1	21	Sentinel-2	2016.12.08	-
6	ALOS-2	2018.02.04	SM1	22	Sentinel-2	2017.04.07	-
7	ALOS-2	2018.02.18	SM1	23	Sentinel-2	2017.11.08	-
8	ALOS-2	2018.05.13	SM1	24	Sentinel-2	2018.01.17	-
9	ALOS-2	2018.06.24	SM1	25	Sentinel-2	2018.04.02	-
10	ALOS-2	2018.08.05	SM1	26	Sentinel-2	2018.11.03	-
11	ALOS-2	2018.09.16	SM1	27	Sentinel-2	2019.02.06	-
12	ALOS-2	2019.01.20	SM1	28	Sentinel-2	2019.04.22	-
13	ALOS-2	2019.03.03	SM1	29	Sentinel-2	2019.09.29	-
14	ALOS-2	2019.04.14	SM1	30	Sentinel-2	2019.12.08	-
15	ALOS-2	2019.07.07	SM1	31	Sentinel-2	2020.01.22	-
16	ALOS-2	2019.09.29	SM1	32	Sentinel-2	2020.03.22	-

3.2. Methods

To reveal the surface changes and deformations of landslides caused by underground mining, we adopt different approaches in the investigation of low and high deformations by fusing satellite and aerial images shown in Figure 2. First, the red-green-blue (RGB) composite optical images taken by satellite and UAV are used to investigate the surface changes through visual comparison. For low deformation, phase-based InSAR is used to acquire the 1-D deformation in the LOS direction. For high deformation, changes in key blocks are investigated using UAV 3-D real scene models and LiDAR DEM. Then, the 1-D elevation change in the vertical direction is monitored by multi-temporal UAV/LiDAR DEMs. The 1-D deformation in EW and NS directions are obtained by optical offset-tracking technique. The 1-D deformation in LOS and AZI directions are obtained by amplitude-based SAR offset-tracking method. Finally, we combine the multi-dimensional deformation observations acquired by SAR and optical offset-tracking techniques to decompose the real 3-D deformation of the landslide surface.

Figure 2. The flowchart of surface changes investigation and deformation monitoring for Jianshanying landslide.

3.2.1. 1-D deformation in vertical by multi-temporal DEMs difference

Since the multi-temporal DEMs are obtained and processed by different methods, the spatial references are inconsistent. Therefore, it is necessary to match them into the same coordinate system. We determine the first DEM in April 2019 as the reference DEM image and the other two DEM images are matched to the reference image with ground control points, which are manually selected by combining the corresponding optical images. The number of control points is to be more than three. After that, the transformation formula is established by the least-squares fitting function.

(1) $\begin{array}{c}\left[\begin{array}{c}x{ }^{\mathrm{\prime }}\\ y{ }^{\mathrm{\prime }}\end{array}\right]=\left[\begin{array}{c}\mathrm{\Delta }x\\ \mathrm{\Delta }y\end{array}\right]+m\left[\begin{array}{cc}cos\alpha & -sin\alpha \\ sin\alpha & cos\alpha \end{array}\right]\left[\begin{array}{c}x\\ y\end{array}\right]\end{array}$(1)

where $x$ and $y$ are the coordinates of the columns and rows in the original DEM image to be aligned, $x{ }^{\mathrm{\prime }}$ and $y{ }^{\mathrm{\prime }}$ are the coordinates in the reference DEM image, $m$ is the scaling parameter, $\alpha$ is the rotation parameter, $\mathrm{\Delta }x$ and $\mathrm{\Delta }y$ are the translation parameters. Then, the DEM images under the uniform spatial reference can be obtained by using the nearest neighbor interpolation. We calculate the vertical deformation by subtracting the elevation values between the two DEMs. Finally, a first-order fitting polynomial is built for each differential DEM to remove the trend term.

3.2.2. 1-D deformation in LOS by InSAR

Based on SAR images, the low surface deformation of the Jianshanying landslide in LOS direction is obtained by the InSAR technique. First, interferometric pairs are obtained by setting spatial and temporal baseline thresholds. Then, the interferograms filtering is performed using the Goldstein method (Goldstein and Werner 1998), which is a very effective filtering method provided by the GAMMA software. This method can improve the fringe visibility and reduce phase noise significantly, which is beneficial for phase unwrapping. Phase unwrapping is applied using the minimum cost flow (MCF) method (Werner et al. 2002). Subsequently, the trend term and atmospheric phase are removed to obtain high-quality interferograms. Finally, the deformation rate maps are acquired by using Stacking interferograms method (Lyons and Sandwell 2003), which can significantly weaken the random errors in time.

3.2.3. 1-D deformation in EW and NS directions by optical images offset-tracking

The surface EW and NS deformations between the two optical images can be obtained by calculating the phase difference after Fourier transforms. The deformation between images can be obtained when the number of cross correlation peaks at sub-pixel accuracy (Leprince et al. 2007; Rott et al. 1998). We calculated the surface deformation using MicMac software (Rupnik, Daakir, and Deseilligny 2017), which calculates the deformation by using regularization. Therefore, it allows to obtain high accuracy deformation with small sliding windows (Rosu et al. 2015), which is critical for the analysis of small-scale landslides.

3.2.4. 1-D deformation in azimuth and LOS directions by SAR offset-tracking

The SAR offset-tracking method can also be used to recover 1-D landslide deformation in the AZI and LOS directions, based on a principle similar to the cross-correlation of optical images, which is based on amplitude information of SAR images. In the case of landslides with large deformation gradients, there are many low-coherence targets whose offset values change significantly with the window size (Bamler 2000). To solve this problem, multiple-deformation results are obtained using variable windows, and finally the best result is obtained among all offset results based on the median operation (Liu et al. 2020; Zhao, Lu, and Zhang 2013).

3.2.5. 3-D deformation decomposition by combining optical and SAR offset-tracking measurements

The 1-D deformations in the AZI and LOS, and EW and NS directions are obtained from the SAR and optical offset-tracking methods, respectively, and the 3-D deformations can be calculated jointly. For this purpose, we first resample the high-resolution SAR results to maintain consistency with the low-resolution Sentinel-2 results. Based on the geometric relationship of SAR satellites, the observation equation between AZI, LOS deformation, and 3-D surface deformation can be established (Li et al. 2020; Wang and Jónsson 2015). Then, the observation equation of 3-D deformation is built by combining the EW and NS deformations from optical offset-tracking, shown as equation (2)(2) $\begin{array}{c}\left[\begin{array}{ccc}cos\theta & sin\phi sin\theta & -cos\phi sin\theta \\ 0& cos\phi & sin\phi \\ 0& 1& 0\\ 0& 0& 1\end{array}\right]\left[\begin{array}{c}{D}_{ud}\\ D_{ns}\\ D_{ew}\end{array}\right]=\left[\begin{array}{c}{D}_{LOS}\\ D_{AZI}\\ D_{NS}\\ D_{EW}\end{array}\right]\end{array}$(2) .

(2) $\begin{array}{c}\left[\begin{array}{ccc}cos\theta & sin\phi sin\theta & -cos\phi sin\theta \\ 0& cos\phi & sin\phi \\ 0& 1& 0\\ 0& 0& 1\end{array}\right]\left[\begin{array}{c}{D}_{ud}\\ D_{ns}\\ D_{ew}\end{array}\right]=\left[\begin{array}{c}{D}_{LOS}\\ D_{AZI}\\ D_{NS}\\ D_{EW}\end{array}\right]\end{array}$(2)

where $D_{AZI}$ is the AZI deformation, $D_{LOS}$ is the LOS deformation, $D_{EW}$ and $D_{NS}$ are the optical offset-tracking deformations in the EW and NS direction, and $D_{ew}$, $D_{ns}$, $D_{ud}$ are the deformations in the east-west, north-south, and vertical directions, respectively, which are unknown parameters. The $\theta$, $\phi$ represent incidence angle and flight azimuth angle.

Moreover, to determine the weights of SAR and optical observations with unequal accuracy, we use Helmert variance component estimation method for iterative weighting processing (Hu et al. 2012; Xu et al. 2009). Specifically, the initial variances are first given, calculated from each central pixel in AZI, LOS, EW, and NS deformations with a fixed window of 3 × 3. Then, an error equation (3)(3) $\begin{array}{c}\left\{\begin{array}{c}{V}_1=A_1\hat{X}-L_1\\ V_2=A_2\hat{X}-L_2\end{array}\right.\end{array}$(3) can be rewritten from formula (2):

(3) $\begin{array}{c}\left\{\begin{array}{c}{V}_1=A_1\hat{X}-L_1\\ V_2=A_2\hat{X}-L_2\end{array}\right.\end{array}$(3)

where $L_1={\left[\begin{array}{cc}{D}_{LOS}& D_{AZI}\end{array}\right]}^T$, $L_2={\left[\begin{array}{cc}{D}_{NS}& D_{EW}\end{array}\right]}^T$, $\hat{X}={\left[\begin{array}{ccc}{D}_{ud}& D_{ns}& D_{ew}\end{array}\right]}^T$, ${\left[\begin{array}{cc}{A}_1& A_2\end{array}\right]}^T$ is coefficient matrix of equation (2)(2) $\begin{array}{c}\left[\begin{array}{ccc}cos\theta & sin\phi sin\theta & -cos\phi sin\theta \\ 0& cos\phi & sin\phi \\ 0& 1& 0\\ 0& 0& 1\end{array}\right]\left[\begin{array}{c}{D}_{ud}\\ D_{ns}\\ D_{ew}\end{array}\right]=\left[\begin{array}{c}{D}_{LOS}\\ D_{AZI}\\ D_{NS}\\ D_{EW}\end{array}\right]\end{array}$(2) . Accordingly, the unit weight variance of the two types of observations ${\sigma }_{01}^2$ and ${\sigma }_{02}^2$ after iteration are obtained (Hu et al. 2012). Then, the weights are determined again according to formula (4) based on the unit weight variance of the two types of observations after iteration.

(4) $\begin{array}{c}{P}_i^{\left(k+1\right)}=\frac{c}{{\hat{\sigma }}_{0i}^2}{P}_i^{\left(k\right)}\end{array}$(4)

where $k$ is the number of iterations, $P_i$ is the weighted value of $i$-th type of observation, and $c$ is a constant. Then, iterations are performed until the unit weight variance of the two types of observations is equal. Finally, the 3-D surface deformation is obtained.

4. Results and analyses

4.1. Historical surface change of Jianshanying landslide

4.1.1. Time series of surface changes before large-scale deformation

As recorded, the Jianshanying landslide has deformed since 2003. However, the earliest high-resolution archived optical RGB images were acquired on 24 March 2013. Accordingly, the time series of surface changes detected by Google earth images from 2013 to 2018 are shown in Figure 3(a–d). It can be seen that the debris of the Jianshanying landslide was clearly formed in 2013 and occurred mainly in three sections. Since 21 December 2015, these sections were connected into one whole. Subsequently, material kept sliding down and accumulating at the toe of the slope. Starting in 2017, localized collapse occurred in the back edge area, as shown in Figure 3c. Then, the boundary of the trailing edge did not change significantly until 24 December 2018. There was no obvious large-scale slope movement during 2017 and 2018, and the deformation was mainly focused on the formation of surface cracks and slope damage at the front edge.

Figure 3. Time series of surface changes detected by Google earth optical images. (a), (b), (c) and (d) are Google earth images, acquired in March 24, 2013, December 21, 2015, April 6, 2017 and December 24, 2018, respectively, coloured lines present main scarp boundary in different dates.

4.1.2. Comparison of large-scale landslide movement

After 2019, the Jianshanying landslide entered a stage where the trailing edge slope suffered a significant downward movement. To further study the large-scale displacement of the trailing edge area, we first made a comparison between 3-D real scene models obtained on 19 April 2019 and 2 August 2020 by UAV. Landslide surfaces on different dates are shown in Figure 4(a,b), respectively, where the trailing edge boundary was marked with a white dash line. In 2019, the front landslide body was basically aligned with the back slope, and local collapses occurred in this area, along with multiple cracks. However, compared with 2019, the front body slid obviously in 2020, and the main body slid along the slope direction, forming the head and scarp area, as shown in Figure 4b. We selected four feature points for elevation comparison, three of which were on the main body, and one was at the trailing edge of slope, as shown in Figure 4(a,b). For the three points on the slope, there was a decrease in elevation relative to 2019, but Point 4 was stable during that period.

Figure 4. The comparison of large-scale movement characteristics of Jianshanying landslide using UAV 3-D real scene models and LiDAR DEM. (a) and (b) are 3-D real scene models acquired on April 19, 2019 and August 2, 2020, respectively, the white dash line is trailing edge boundary, points 1–4 are typical surface points in two stages, where points 1–3 on the landslide body and other one on the trailing edge, the elevations are given on the point.

Figure 4c is LiDAR DEM acquired in June 2019. The solid lines were landslide boundaries and deformation boundaries of the posterior edge, and the white rectangle was an enlarged map of cracks, area I was mainly deformation region, area II was new landslide area after overall slide in September 2020. Then, we surveyed the landslide surface characteristics after 2019 by LiDAR DEM in Figure 4c. In the mainly deformed area I, there were two large-scale sliding areas and two small-scale sliding areas. The blocks slid and formed obvious deposition area, which can be clearly seen through LiDAR DEM. For area II, we cannot see any obvious signs of a slide except for the southeast area, which was a small-scale landslide. Therefore, a large amount of material at the back of the landslide body slides forward, which gave an increased load for the front edge slope and was the main factor for the eventual sliding in area II.

4.1.3. Surface characteristics after overall slide

We investigated the surface characteristics after the overall slide in September 2020 by combining Google earth images and field survey. After the landslide slid, the morphology changed significantly in the image acquired on 11 November 2020, as shown in Figure 5a, where roads and houses at the toe of landslide were distinctly destroyed and the landslide influence scope was expanded widely. Areas I and II in Figure 4c were connected as a result of the overall slide. Figure 5b shows a full view of the post-event Jianshanying landslide, where existed clear signs of rock fragmentation and sliding in the rear, and it can be seen that the intact rear blocks were sliding forward. Figure 5c shows the main source area, which was sliding along the slope as a whole, but had not yet completely fallen. It was still in a stage of high deformation and contained many cracks. In the middle and lower parts of the landslide, there was mostly rubble. Figure 5d shows large rock deposits at the toe of the landslide. Figure 5e shows the head area, where the main scarp can be clearly seen. Thus, after the overall instability occurred, the large deformation of slope was still mainly concentrated in the upper slip source area, which was still in the stage of continuous movement.

Figure 5. Surface characteristics of landslide after overall slide in September 2020. (a) Google earth image, (b)-(e) field photos, which are taken at locations shown as blue symbols in (a).

4.2. Multi-dimensional deformation of Jianshanying landslide

4.2.1. 1-D vertical deformation

The vertical deformation was obtained by differentiating the DEMs of three dates from April 2019 to August 2020, as shown in Figure 6, with the latter date subtracting the former one. The magnitude of the vertical deformation was on tens of meters. The precision is better than the meter level, which can describe the landslide characteristics with such high deformation. We first focused on DEM changes after April 2019, the main deformed area in Figure 6(a–c) can be divided into two parts, the top and bottom deformation zones of slope, which were shown in red and blue color, respectively. The results showed that the elevation at the top of the slope has decreased about −30 m and increased approximately 29 m at the bottom of the slope. In order to reveal the long-term deformation of the landslide, we extracted the time-series DEM differences of the point P, the Figure 6d shows that the landslide experienced relatively slow deformation from April 2019 to June 2019 and entered rapid deformation from June 2019 to August 2020. Therefore, we subsequently focused on the high deformation phase. We further extracted the DEM changes along profile A-A’ from the three results, as shown in Figure 6e. To reveal the causes of deformation, we focused on three typical deformation zones I, II, and III in Figure 6e. For zone I, the DEM change decreased, indicating that downward sliding along the slope direction or subsidence has occurred in this area. However, for zone II, an increase in elevation was observed in some areas between April 2019 and June 2019, indicated that the material at the trailing edge slid into this area for a short period of time, but there was no large-scale sliding in the entire area. Furthermore, for the deformation from June 2019 to August 2020, there was a decrease in elevation. This indicated that the landslide deformation in this area was mainly sliding along the slope or settlement. For zone III, the increased elevation change was mainly caused by the accumulation of trailing edge material. The surface deformation in these sections can also provide an important reference for analyzing the landslide failure mechanism.

Figure 6. Surface vertical deformation obtained by DEM difference. (a) shows the LiDAR DEM in June 2019 minus the UAV DEM in April 2019. (b) shows the UAV DEM in August 2020 minus the LiDAR DEM in June 2019. (c) shows the UAV DEM in August 2020 minus the UAV DEM in April 2019. (d) are time-series DEM differences at point P in (a) and (b). (e) are the DEM differences along profile A-A’ in (a)-(c).

4.2.2. 1-D LOS deformation

The LOS deformation rate obtained by InSAR is shown in Figure 7a, where the deformation along the line of sight was negative. The magnitude of the deformation was concentrated on the decimeter level. The precision of the results is up to millimeters in low deformation areas. It is apparent that the surface of the landslide area suffered obvious deformation, and the monitored deformation was mainly concentrated at the front and back edges of the slope. The maximum deformation rate was −0.3 m/yr and was located on the northwest side. Meanwhile, we show that the InSAR method failed to monitor the high deformation of the main landslide body, where no sufficient effective deformation points can be captured. This phenomenon occurs because the gradient of surface deformation exceeds the maximum limit of InSAR and that causes the loss of coherence (Chen et al. 2021). Therefore, only the low deformation area can be monitored with high accuracy by the InSAR method in the 1-D LOS direction.

Figure 7. The deformation rate maps comparison among InSAR and SAR offset-tracking methods. (a) the LOS deformation obtained by Stacking-InSAR. (b) and (c) the AZI and LOS deformation rate maps acquired by SAR offset-tracking.

4.2.3. 1-D deformation in AZI and LOS directions

Figure 7(b,c) showed the 1-D deformations in the AZI and LOS directions, respectively, acquired by an improved SAR offset-tracking technique. The magnitude of deformation was on the level of meter. The precision is up to 1/20–1/10 pixels (Raucoules et al. 2013). The deformation along the azimuth direction in Figure 7b was positive, and along the LOS direction in Figure 7c was positive. Firstly, in Figure 7(b,c), the high deformation field of the Jianshanying landslide can be monitored completely, especially for the area where the InSAR method failed. Secondly, the low deformation area can be monitored successfully with InSAR method in Figure 7a, especially, in the rear areas of landslide, where it can hardly be detected with SAR offset-tracking method in Figure 7c. The large deformation retrieved in Figure 7(b,c) mainly occurred in the heavily damaged area, which is on the back edge of the slope. The deformation direction was consistent with the topographic change direction. For instance, the LOS deformation showed a long strip, which was consistent with the topography of the trailing edge shown in Figure 4. The maximum deformation rate in both directions reached 12.3 m/yr, which suggests that this method can successfully obtain 1-D deformation in high deformation landslide.

4.2.4. 1-D horizontal deformation

The surface deformation in the EW and NS directions of the Jianshanying landslide from 31 July 2016 to 22 March 2020 was recovered by the optical offset-tracking method. The magnitude of deformation was also at the level of meter. The precision was around 1/10 pixel. The horizontal deformation field is generated by combining the accumulated deformation in both directions, as shown in Figure 8. The maximum deformation monitored was located on the landslide body, with values exceeding 45 m. From Figure 8, it can be seen that the horizontal deformation direction was mainly concentrated in the northeast direction, which agreed well with the main slope direction. The results verified that the optical offset-tracking method can reliably recover the high 1-D deformation of the landslide in case of gravel collapse and excessive deformation. However, it cannot monitor vertical deformation.

Figure 8. The horizontal deformation obtained by Sentinel-2 images based on optical offset-tracking method.

4.2.5. 3-D surface deformation

The 3-D surface deformation can be obtained by fusing optical and SAR offset-tracking observations. In order to obtain more accurate results, the cumulative deformation during the overlapping time of SAR and optical data were chosen for calculation. The optical images were acquired from 29 September 2019 to 22 January 2020, and the SAR data were acquired from 29 September 2019 to 19 January 2020. There was only a three-day difference between optical and SAR observations, so it was reasonable to be considered as the same duration. The 3-D deformation is shown as Figure 9, where the color represents the deformation in the vertical direction, and the positive value indicated the uplift of the surface relative to the previous surface, and the negative values were the subsidence. Meanwhile, the arrows show a horizontal deformation, whose length corresponds to the magnitude. It can be seen from Figure 9 that the vertical deformation can be divided into three parts, namely, the uplift zone of the eastern slope, the subsidence zone of the slope, and the uplift area of the trailing edge of landslide. The first part was mainly caused by the forward movement of the upper sliding slope. The second part is well consistent with the regular failure mode of the landslide, where the surface falls vertically with the sliding of the slope. In the third part, the trailing edge of the landslide and the surface settled, but the decomposition result is uplifted, and this phenomenon may be caused by the low accuracy of the offset-tracking method in the region of large material change. Meanwhile, the horizontal deformation coincides well with the slope direction, indicating that the results obtained by this method are reliable.

Figure 9. The 3-D surface deformation of Jianshanying landslide. The arrow and length indicate the direction and magnitude of horizontal deformation, respectively. The colour represents the vertical deformation.

To reveal the different characteristics of surface deformation in different parts of the landslide body, the horizontal deformation obtained from the decomposition was analyzed along the profile line, and the location of the profile was shown as line X-X’ in Figure 9. The 2-D deformation along the profile is shown in Figure 10. For the deformation in the NS direction, it was basically concentrated in the northward direction. The large deformation section was between about 200 and 700 m. The maximum deformation was approximately 6 m between 29 September 2019 and 22 January 2020. For the deformation in the EW direction, it mainly moved eastwards, with a maximum deformation close to 3 m. Similarly, in the NS direction, the largest deformation occurred between 200 and 700 m. Overlaying the underground mining range shown in Figure 1, we discovered that below the surface of 200 and 700 m was underground mined-out area, the boundaries of which were in good agreement with the abrupt change points of the 2-D horizontal deformation, as shown in Figure 10. This indicates that underground mining is the main triggering factor for this landslide deformation (Chen et al. 2021; Li et al. 2022).

Figure 10. The horizontal deformation distribution of Jianshanying landslide along profile X-X.

As for the UD deformation, its spatial deformation is relatively complex. To reveal the precision-influencing factor of the UD deformation, the UD deformation along profile X-X’ was analyzed by combining the DEM changes and two UAV 3-D real scene models between April 2019 and August 2020 (Figure 11). It can be seen that the small error areas mainly focused on sections A1 and A3. In section A1, a low deformation occurred and the surface change was around 2 m. Secondly, for section A3, the vertical deformation was about −2 m, but the surface material change was relatively small. Thirdly, the large deformation errors mainly occurred in section A2, which was caused by two factors. On the one hand, the largest DEM changes occurred, which meant that a large amount of surface material has broken up and moved down. On the other hand, although the height change in the transition zone was not significant, the material also transformed from high to low altitude. It can be seen form Figure 11(b,c) that the surface was heavily eroded, the material in this section also changed greatly. In summary, the accuracy of vertical deformation depends on the degree of surface material change. The greater the change of the surface material, the lower the accuracy of deformation.

Figure 11. The analysis of vertical deformation and height changes. (a) the UD deformation with topography and DEM change between April 2019 and August 2020. (b) and (c) the UAV 3-D real scene model acquired in April 2019 and August 2020, respectively. The location of profile X-X’ is also superimposed on (b) and (c).

5. Discussion

5.1. Underground mining-induced subsidence revealed by DEM difference

To explore the surface changes in the vertical direction, which are crucial to reveal the failure mechanism of mining-induced landslides, we calculated the changed slope volume of the Jianshanying landslide based on the DEM difference shown in Figure 6c. We first calculated the elevation changes of each pixel and then multiplied them by the pixel area 1 m² to obtain the change volume of each pixel. Thereafter, we analyzed the distribution characteristics of the changed volume (Figure 12). Otherwise, we also counted the average elevation of the pixels within a certain interval with the aim of finding the main changed area of the slope. In Figure 12, there was a significant volume difference between increased and reduced area, the total volume was −14.0 × 10⁵ m³ and 9.2 × 10⁵ m³, respectively. The difference in change reveals that the amount of material alteration on the surface within the landslide region is not conserved.

Figure 12. The comparison of increased and reduced slope volume. This graph is based on the statistics of DEM difference in Figure 6c. The X-axis shows the DEM difference value, the Y-axis represents the total area that fall within an interval of ±2 m, with the X-axis scale as the centre point. The size of the circle displays the size of the volume, with larger circle indicating larger volumes. The colour indicates the average elevation of pixels that fall within a certain interval. In the left of the grey line, the circles represent reduced slope volume, in the right, the circles show increased slope volume.

This phenomenon quantitatively revealed that, in addition to surface sliding, there was significant ground subsidence within the landslide area, with a total reduced volume of 4.8 × 10⁵ m³. Otherwise, as for the volume reduced area in the left of Figure 12, it was found that reduced volume caused by subsidence was mainly concentrated in high-elevation area, which was the rear area of the back edge of the Jianshanying landslide. This finding has important implications for revealing the failure mechanism of this type of mining-induced landslide.

5.2. Failure mechanism analysis

The 3-D surface deformation allows us to realize the evolutionary trend of different parts of the landslide to clarify the failure mechanism. We extracted horizontal, vertical deformation and DEM difference along seven profile lines, as shown in Figure 13(a–d) which spanned the major high deformation zone of the landslide. Since we have analyzed that the vertical deformation inverted by combining SAR and optical offset-tracking has low accuracy in the area where surface material changed greatly, we replaced the vertical deformation in Figure 13c with DEM difference in Figure 13d to analyze mechanism. Figure 13e demonstrates that little subsidence happened in area I, but enormous differential subsidence occurred in area II. For the horizontal direction in Figure 13f, the deformation was larger in area I, while it was relatively smaller in region II. Based on the characteristics of these differential deformations, the failure mechanism can be summarized. During the underground coal mining process, the rock layers directly connected to the coal seam lose their support and collapse after the formation of the goaf. With the expansion of coal seam, the lower rock layers continuously bend and sink, and the upper rock layers also begin to bend and sink under its own weight, resulting in the rapid transmission up to the surface and forming ground subsidence, as shown in area II of Figure 13e. During this process, the stress change in this area is characterized by vertical stretching and horizontal compression, which is fully verified by 3-D deformation in Figure 13(e,f). In contrast, as the front slope of the landslide, the stress change in area I is characterized by vertical compression and horizontal stretching. Due to the subsidence of the rear rock mass and the presence of a cliff in the front, this stress change reduces the strength of the slope, forms sliding surfaces in weak strata, and causes slope sliding with the effect of gravity. In addition, the difference in stress change between area I and area II causes tension cracks at the back edge of the landslide, through which precipitation enters the interior of the landslide, thereby accelerating the evolution of the landslide. This finding is basically consistent with the numerical analysis work presented by Li et al. (Li et al. 2022).

Figure 13. Failure mechanism of Jianshanying landslide revealed by 3-D surface deformation. (a) the locations of profile line in UAV image, the colours of each line are consistent with the colour in (b), (c) and (d). (b) and (c) are horizontal and vertical deformation of Figure 9 along profile, respectively, the horizontal deformation is calculated by synthesizing the projection of NS and EW deformation. (d) is DEM difference profile of Figure 6c. Gray arrows are display the points where deformation change significantly. (e) and (f) are failure mode of Jianshanying landslide, which superimpose DEM difference and horizontal deformation of line 3–3’ in (a), respectively, area I is landslide body, area II is subsidence area in landslide back scarp area.

6. Conclusions

In this study, multi-sensor RS approaches were applied to investigate and monitor the typical underground mining-induced Jianshanying landslide to trace the pre-event surface motions and reveal its failure mechanism. It was found that small-scale slope movement occurred during 2013 and 2018, and the landslide entered a high deformation stage after 2019, during which a large number of cracks developed and formed a landslide morphology after the overall slide in September 2020. In this study, InSAR was not able to effectively monitor landslide deformation because the wide gap between SAR acquisitions. In this case, only 19 archived ALOS/PALSAR-2 images over 3 years strongly caused the surface deformation to exceed the maximum limitation of InSAR. Therefore, the 3-D surface deformation was calculated by fusing optical and SAR offset-tracking observations, but the accuracy of vertical deformation is influenced by the change of surface material. Otherwise, in addition to surface sliding, there was a significant ground subsidence in the back scarp area with a total volume reduction of 4.8 × 10⁵ m³. Based on this phenomenon, we revealed the failure mechanism that the subsidence of the rear of the landslide body caused by underground mining resulted in differential stress changes in the horizontal and vertical directions between the rear slope and the front landslide body, which eventually triggered the landslide at the front slope along the free surface. This study provides a perspective for understanding the mechanism of mining-induced landslides from the perspective of deformation observations, and the proposed failure mechanism can help to prevent such disasters. In addition, this study proposes a novel method to achieve 3-D high deformation monitoring, and reveals the relationship between DEM changes and deformation changes. It provides a new tool for the monitoring of hazards in the world with high deformation characteristics, such as volcanoes, underground minerals, and oil and gas extraction, and further contributes to the research on the mechanism and early warning of geohazards.

Acknowledgments

The authors would like to thank JAXA, Japan for providing the ALOS/PALSAR-2 data and the European Space Agency for providing the Sentinel-2 optical images free of charge. In addition, thanks to MicMac software, which can be download from website https://github.com/micmacIGN/micmac.

Disclosure statement

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

Data availability statement

The Sentinel-2 images used in this paper are available from the website https://scihub.copernicus.eu/dhus/#/home. The ALOS/PALSAR-2 data can be acquired from https://gportal.jaxa.jp/gpr/search?tab=1 after applying to JAXA, Japan. Other data are available from the corresponding author [Chaoying Zhao] upon reasonable request.

Additional information

Funding

This work was funded by the National Key R&D Program of China (No.2022YFC3004302), the National Natural Science Foundation of China (Grant No. 41929001), and the Fundamental Research Foundation of the Central Universities (No. 300203211262).