Pre ‐ Event Deformation and Failure Mechanism Analysis of the Pusa Landslide, China with Multi ‐ Sensor SAR Imagery

: The Pusa landslide, in Guizhou, China, occurred on 28 August 2017, caused 26 deaths with 9 missing. However, few studies about the pre ‐ event surface deformation are provided because of the complex landslide formation and failure mechanism. To retrieve the precursory signal of this landslide, we recovered pre ‐ event deformation with multi ‐ sensor synthetic aperture radar (SAR) imagery. First, we delineated the boundary and source area of the Pusa landslide based on the coherence and SAR intensity maps. Second, we detected the line ‐ of ‐ sight (LOS) deformation rate and time series before the Pusa landslide with ALOS/PALSAR ‐ 2 and Sentinel ‐ 1A/B SAR imagery data, where we found that the onset of the deformation is four months before landslide event. Finally, we conceptualized the failure mechanism of the Pusa landslide as the joint effects of rainfall and mining activity. This research provides new insights into the failure mechanism and early warning of rock avalanches.


Introduction
On 28 August 2017, the long-runout collapse initiated by the ridge-top rockslide in Pusa Village, Zhangjiawan Town, Guizhou Province, China, buried residential areas and caused 26 deaths with 9 missing [1]. This catastrophic disaster is a typical rock avalanche caused by combined effects of natural and anthropogenic factors in the Yunnan-Guizhou Plateau and its surrounding areas, China, which is the largest karst mountain area in China, and even in the world, occupying around 6.2 × 10 5 km 2 [2]. Due to the vulnerable geological environment and frequent anthropogenic activities in this area, more than 2800 people have been killed by landslide disasters in recent decades. The landslide showed complex geological structures and uncertain triggering factors, which makes it difficult to invert the failure mechanism. The surface deformation before the landslide and its correlation with external factors can hopefully reveal the failure mechanism and triggering factors.
Synthetic aperture radar (SAR) intensity and coherence maps have been widely used to study the landslide disasters owing to the area being free from clouds. Zhao et al. [3] detected the landslide affected area with SAR intensity maps. Mondini et al. [4] detected worldwide rapid landslides by Sentinel-1 SAR intensity maps. Raspini et al. [5] revealed the pre-event movement of the slope in Montescaglioso village of South Italy by SAR intensity information. Yun et al. [6] acquired damage proxy maps through changes in coherence maps, which were generated by pre-and post-event SAR images. Dai et al. [7] combined interferograms and their corresponding coherence maps to identify the source region and boundary of the landslide.
In addition, interferometric SAR (InSAR) can monitor the surface deformation with centimeter to millimeter precision [8,9]. It has been extensively used to monitor high-risk slopes and restore their deformation history [10][11][12][13]. Further, landslide types and failure modes can be deduced based on the pre-event deformation, which will provide great guidelines for landslide prevention and hazard mitigation [10,[14][15][16]. To this end, Kang et al. [10] recovered the failure mode of the Guanling landslide in China and found that the shear failure in the locking segment of a key body triggered by heavy precipitation based on the pre-event deformation features was the cause. Dong et al. [14] and Intrieri et al. [17] captured a sudden acceleration of deformation dozens of days before the landslide of Xinmo landslide in Maoxian County, Sichuan Province through InSAR technique, which was triggered by torrential rainfall. Besides, Zhao et al. [3] addressed the mining-induced Jiweishan landslide failure mechanism through L-band ALOS/PALSAR-1 datasets. As for the loess landslides' failure modes, Shi et al. [18] and Zhao et al. [19] revealed the irrigation triggering factor by using fusing InSAR and groundwater measurements. In addition, four loess landslides types were summarized by Liu et al. [16] based on two-dimensional surface deformation generated from ascending and descending TerraSAR datasets. Tantianuparp et al. [20] uncovered landslide failure caused by the temporal variation of water level in the Three Gorges, China, using multi-sensor SAR imagery. Bru et al. [21] applied InSAR techniques to the operational geotechnical monitoring of the Leintz Gatzaga landslide in Basque Country, Spain. Hu et al. [22] detected the deformation of seasonal landslides in Washington by time-series SAR imagery. Xu et al. [23] revealed the movement of a rainfall-driven landslide in southwestern Oregon based on soil moisture and InSAR time series.
However, as for the triggering factors and failure mechanism of rock avalanches of southwestern China, it was only analyzed based on the geotechnical engineering technique [24]. How can the landslides be identified and delineated based on the SAR images? Is there any precursory signal that can be captured? And what is the failure mechanism of the rock avalanches? To answer these questions, we took the Pusa landslide as an example, where intensity, coherence, and phase-based information were all involved. And multi-sensor SAR imagery was employed. In addition, both ascending and descending SAR data and in-situ data were collaborated to analyze the triggering factors and the failure mechanism of the Pusa landslide. Figure 1 shows the location of Pusa landslide and the coverage of SAR images. The average temperature in a year is approximate 13.6 ℃ ranging from 4 ℃ to 32 ℃. And the annual rainfall is abundant from May to September, which accounts for 74% of the yearly cumulative rainfall ranging from 1200 mm to 1300 mm. The bird's eye view of the Pusa landslide is shown in Figure 2a, where A, B, and C indicate the source area, track-abrasion area, and deposition area, respectively [25]. The depth of the sliding surface is 85 m. The on-site photos of this landslide are shown in Figure 2b-d, where the surface of the back scarp area is fragmented, and the tension crack is visible at the back edge of the remaining blocks. The scarp after the landslide is about 2100 m a.s.l.. The elevation of the deposition area is 1833 m to 1930 m. The volumes of rock material in the source area and deposition area are about 5.0 × 10 5 m 3 and 8.23 × 10 5 m 3 , respectively [1].

Study area
The geological settings of the Pusa landslide is shown in Figure 3, among which the exposed lithologic formations from the top to the bottom are the Yelang Formation of the Lower Triassic (T1y), about 100 m thick, the Changxing-Dalong Formation of the Upper Permian (P2c+d) and the Longtan Formation of the Upper Permian (P2l). The thickness of P2c+d is about 50 m. The T1y consists of limestone, silty mudstone, and argillaceous siltstone. The P2c+d consists of silty mudstone, limestone, argillaceous siltstone, and coal seams. The P2l consists of argillaceous siltstone, carbonaceous mudstone, and coal seams [1]. Besides, there are six coal seams under the foot of the landslide, where the mining behaviors were mainly conducted in CS16 coal seams, resulting in several goafs underground. The formation of these goafs changed the stress state of the slope, which had the potential to cause the deformation of the slope.

Data
Two SM3 ALOS/PALSAR-2 images acquired on 4 September 2016 and 22 January 2017, and three SM1 ALOS/PALSAR-2 images from 14 May 2017 to 6 August 2017 were used to detect preevent deformation area of the Pusa landslide ( Figure 4a). Two post-event SM1 ALOS/PALSAR-2 images acquired on 12 November 2017 and 24 December 2017 were used to identify the boundary and source area of the landslide. Twenty ascending and eighteen descending Sentinel-1A/B images were used to calculate the surface deformation of the Pusa landslide before it occurred (Figure 4b).
The parameters of SAR images are listed in Table 1. The shuttle radar topography mission (SRTM) digital elevation model (DEM) with a resolution of 30 m was employed to remove the topographic phase in InSAR processing.

Methods
The flow chart for Pusa landslide analysis with multi-sensor SAR imagery is shown in Figure 5, which includes three main functions. Coherence map and intensity maps were applied to identify the landslide boundary, source area, and buried villages, while deformation maps from both L-band ALOS/PALSAR-2 and Sentinel-1 were used to retrieve pre-event deformation.

Coherence Estimation
Coherence can be used as a measure for the accuracy of the interferometric phase [26]. The coherence magnitude is the most common measure of coherence with the range of [0,1]. Under good SAR imaging conditions, it is close to 1 for bare surfaces, but close to 0 in dense or thick vegetation, such as forests or crops. Based on these characteristics, we investigated the region where there was a vegetation-covered area before the landslide occurred but was covered by deposition after the landslide. Accordingly, we could delineate the boundary of the landslide.
can be calculated commonly using Equation (1) based on SAR image amplitude [26,27]: where is the number of pixels in a window, and are the complex value of master and slave images at pixel , respectively, * denotes complex multiplication, |•| denotes the absolute value of the complex.

Surface Change Detection with Intensity Maps
The SAR intensity is affected by factors such as dielectric, surface roughness, and terrain slope [28]. The ratio of SAR intensity maps reflects the surface change during the period of the repeat-pass SAR images [3,29,30]. Therefore, we delineated the extent of the source area of the landslide by using two intensity maps acquired before and after the event.
To detect the small scale landslide, we adopted single-look intensity maps with a spatial resolution of 3 m. However, to decrease the serious speckle noises in the intensity maps, we filtered the intensity maps with the modified Lee filter algorithm based on statistically homogeneous pixels (SHPs). First, we selected SHPs by the modified pixel clustering algorithm based on the hypothesis test and confidence interval method (HTCI) [31]. Then the filtering weight is calculated by the local variance, mean, and coefficient of variation within the detected SHPs, which was based on a span image. The span image was defined as the sum of all intensity maps. Finally, the images were filtered by using the modified Lee filter [32] as Equation (2).
where I , , … , is a vector of intensity maps, I ̅ ̅ , ̅ , … , ̅ is the mean value of intensity maps, which is calculated by averaging the intensity of the detected SHPs, is the number of SAR intensity maps, represents the filtering weight, I is the intensity vector after filtering. After removing the effects of speckle noises, the ratio between two intensity maps can be calculated by using Equation (3) where and are backscatter coefficient values for pixel in two intensity maps, and represents the window size.

Stacking Interferograms
The average deformation rate map is often used to identify the potential landslide with subtle deformation. The stacking interferograms method can be used to acquire the average deformation rate with high accuracy. As the atmospheric delay often causes temporally random errors in interferograms, it can be reduced by averaging the unwrapped interferograms as follows [10,33,34]: where ℎ is the average deformation rate in phase, and represent the unwrapped phase and the temporal interval of -th interferogram, respectively, and is the number of interferograms. In this study, the unwrapped interferograms, which have high coherence and no unwrapping errors at the landslide region, are selected to calculate average deformation.

SBAS-InSAR
Small baseline subsets (SBAS) InSAR combines SAR images to form interferograms within the limits of spatial and temporal baselines to limit the spatial and temporal decorrelation phenomenon [35]. In this paper, the spatial baseline and temporal baseline thresholds were set to 150 m and 100 days, respectively, for Sentinel-1A/B data. The multi-look ratio was 4:1 in range and azimuth directions, and the coherence threshold was 0.3 to get more monitoring targets. The minimum cost flow (MCF) method was applied to unwrap the interferograms [36]. And the ramp phases were removed by a quadratic polynomial, shown as Equation (5). ℎ , where ℎ , is the unwrapped phase at coordinates , , , , , , , and are the unknown coefficients. We subtracted the unwrapped atmospheric phase by estimating the linear model between the residual phase and elevation. Finally, we selected the high quality unwrapped interferograms to calculate the accumulative deformation and deformation time series.

The Boundary and Source Area of the Landslide Identification
We used coherence and intensity maps to identify the landslide boundary, source area, and buried villages. Figure 6 shows the coherence and intensity maps from ALOS/PALSAR-2 images of the Pusa landslide region. As Pusa village is usually covered by dense vegetation, InSAR coherence often gets lost in summer for both C-band and L-band SAR images. When Figure 6a,b are compared, the pre-event coherence in the Pusa landslide area was low, while the post-event interferograms showed high coherence. The reason is that this area was covered by heavy vegetation originally but changed to bare soil and rock mass after the occurrence of the landslide on 28 August 2017, which can be verified in Figure 6c. Therefore, through the changes in surface coherence before and after the landslide occurred, the landslide boundary was delineated, shown by white and yellow lines in Figure 6, which is well consistent with the boundary (landslide boundary in Figure 7) depicted by post-event in-site image [1]. Furthermore, the source area of the landslide was detected by the two closest intensity maps crossing the event acquired on 6 August 2017 and 12 November 2017, respectively, as shown in Figure 6d and Figure 6e. It can be seen that the areas marked by black lines in Figure 6d show a brighter tone than those in Figure 6e. It can be inferred that these areas were covered by objects which showed strong backscattering characteristics before the landslide occurred. By comparing the Google Earth image, we confirmed these areas were the houses of the local village, which were buried by the deposition. Furthermore, the delineated range of the source area is shown as green lines in Figure 6d-f. It can be seen that the slope in the green line (Figure 6e) had slid relative to the pre-event intensity map (Figure 6d). Meanwhile, the landslide back wall can be clearly seen in Figure 6e. We then calculated the ratio between pre-and post-event intensity maps with the average ratio in the green line ( Figure 6e) up to 0.90, which manifested that the source area of landslide could be well delineated by intensity maps.

Pre-event Deformation from ALOS/PALSAR-2 Datasets
As only two SM3 modes and three SM1 modes ALOS/PALSAR-2 images are available before the landslide, it is unfavorable to calculate the deformation rates and time series, so we executed conventional differential interferometric processing to detect pre-event deformation area. As shown in Figure 7, the local fringes were visible over and around the landslide source region, which indicates that the large deformation in the sliding slope arose before the landslide occurred. The areas of three regions shown as black dash lines in Figure 7b-d were 0.123 km 2 , 0.242 km 2 , 0.310 km 2 , respectively. Furthermore, the extent and the number of local fringes over landslide source area increased significantly with time, especially from 14 May 2017 to 11 June 2017, which indicates that the deformation accelerated after 14 May 2017. This deformation characteristic may verify that the mining activity beneath the landslide changed the stress state within the slope body [1], which would lead to instability in the source area of the landslide. To further reveal the temporal evolution of surface deformation before the landslide occurred, we applied both ascending and descending Sentinel-1 datasets with the SBAS-InSAR technique in the next section.

Pre-event Deformation from Sentinel-1 Datasets
Annual LOS deformation rate maps calculated with ascending and descending Sentinel-1A/B images are shown in Figure 8a,b, respectively. It is worth noting that the negative values (red points) represent the surface moves away from the satellite, while the positive values (blue points) indicate the surface moves towards the satellite sensor. As the Pusa landslide is usually covered by dense vegetation, the interferograms within short temporal baselines were used to calculate the deformation rate by the stacking interferograms method. As shown in Figure 8, obvious deformation before the landslide was visibly seen at the trailing edge and eastern side of the sliding body (source area A). In source area A, many coherent points suffered large deformation, where the maximum deformation rate reached −106 mm/year in the line-of-sight (LOS) direction at the back scarp area of the landslide. It can be evidently seen that the eastern side of source area A showed a large deformation rate with the maximum rate of +83 mm/year. Besides, there was a large land subsidence area at the bottom of the whole slope. Furthermore, the average descending LOS deformation rates, shown in Figure 8b,d, were consistently larger than those from ascending datasets, but the signs were opposite in the eastern and bottom of the slope. Unfortunately, few points could be captured in source area A due to the layover for descending SAR geometry.

Decomposition of the Pre-event Deformation
As we can see, different regions show different deformation values in ascending and descending results. It can be further divided into two cases from Figure 8: The deformations were both negative relative to two flight directions, and the deformations showed opposite signs relative to the different flight directions. Considering the different imaging geometries of the ascending and descending satellites, we deduced the relationships between the LOS deformation and actual slope direction deformation, shown in Figure 9. In this case, the average azimuth direction of the descending and ascending satellites was 190°39 18 " and 350°06 54 " , respectively. As for the first case, a point with deformation as negative in both ascending and descending deformation rate maps, the range of the actual slope deformation direction as for descending track result is shown as in Figure 9a, while the ascending track result is shown as .
is the intersection of and , which is shown as a red shadow, this range is the actual deformation direction for this case. After calculation, the moving direction of this area was approximately between 9°53 06 "~1 0°39 18 " , which approximates north. Similarly, the actual deformation direction for the second case could also be inferred, as approximately between 190°39 18 "~3 50°06 54 " , shown as in Figure 9b. Finally, we put deformation in both cases in Figure 8a to decompose the pre-event deformation of the Pusa landslide and its neighbors. As for the points a1, a2, a3, the arc segments represent the range of actual deformation direction. The arrows represent the probable direction in these areas by combining it with the detected ALOS/PALSAR-2 pre-event deformation. The actual deformation direction of area a1 and a2 moved toward the center of the goaf, while the deformation of area a3 was well consistent with the factual sliding direction of the Pusa landslide. Therefore, the deformation inferred by InSAR measurements consists of two parts: coal mining-induced deformation and landslide deformation. and were the range of moving azimuth angle for the descending track. and were the range of moving azimuth angle for the ascending track. and were the intersections for the descending track and ascending track, respectively.

Pre-event Deformation Time Series
Some previous research has shown that the sliding body suffered constant deformation before the landslide occurred [1,25,37,38]. To study the sliding mechanism of the Pusa landslide, we obtained the cumulative time series deformation of the landslide with the SBAS-InSAR technique. It can be seen in Figure 10 that the ascending InSAR measurements can be applied to capture the precursory deformation, as one SAR image was acquired just five days before the landslide occurred. To analyze the temporal evolution of the landslide deformation, we extracted the deformation time series at three points, which are marked as P1, P2, P3 in Figure 8a,b. P1 and P2 are located on the sliding body and back scarp area of the landslide, while the P3 is located on the top of the landslide. To better uncover the deformation trend, we conducted a curve fitting by using the Gaussian function for the deformation time series (Figure 10). For P1 in Figure 10a, the cumulative deformation can better reflect the deformation of the slope before the sliding. The sliding body showed slight LOS deformation since 3 October 2016, but the cumulative deformation of the slope body turned to a positive value after 13 April 2017. Moreover, it increased continuously and reached +116 mm by 23 August 2017, just 5 days before the landslide. For P2, it can also be seen that the deformation time series increased continuously for both ascending track (Figure 10b) and descending track (Figure 10d).
The maximum deformation reached −106 mm and −120 mm, respectively. For P3 in Figure 10c, we can see the sliding body moved continuously before the landslide occurred with the maximum deformation of −86 mm.

Discussion
As there are no in-situ measurements one year before the failure of Pusa rock avalanche, the preevent InSAR deformation time series can hardly compare quantitatively with external measurements. Even so, the InSAR deformation field can be verified by comparing different modes of ALOS/PALSAR-2 interferograms (Figure 7) and the stacking Sentinel-1 interferograms (Figure 8), where the deformed regions were consistent. Moreover, the deformation rate can be visually compared with ascending and descending Sentinel-1A images, shown in Figure 8, where green points indicate a stable region in ascending and descending results. And deformed regions can be verified by the Google earth image shown in Figure 2(a) before the event. Further, the deformation time series can be validated by ascending and descending results at P2 shown in Figure 10 (b) and (d), where the main deformation trends are consistent, and the difference has resulted from different SAR geometry.

The Triggering Factors
From Figure 3, six coal seams are located beneath the source region of the Pusa landslide, whose cumulative thickness is about 8 m [37]. To analyze the relationship between slope deformation, underground mining, and weather variations, we first compared the ascending deformation time series and deformation rate at point P1 with cumulated rainfall amounts in five days. The daily rainfall was logged by the Zhangjiawan Town rainfall station [1]. To synchronize the temporal resolution of InSAR measurements and rainfall, the deformation time series from 7 January 2017 to 23 August 2017 was selected to conduct the comparison (Figure 11). After April 2017, the Zhangjiawan Town area entered a period of sunny and rainy. The rainy season occurred in late June and the middle of July annually. It can be seen from Figure 11, the deformation of landslide body was affected by climate change, and the deformation between two adjacent heavy rainfall times could be roughly divided into three stages: I, II, III. The deformation of the source area accelerated after a large rainfall, namely stage I. Then, it entered the decreased rate, i.e., stage II, due to the sunny weather. Finally, the deformation rate kept constant under the continuous sunny weather, that is, stage III, which lasted until the next heavy rainfall. During late June and early mid-July 2017, persistent heavy rainfall took place in this area, and the cumulative rainfall increased quickly from 147.9 mm on 21 June to 522.2 mm on 21 July. Followed by persistent sunny days until the heavy rainfall occurred again on 24 August 2017. Accordingly, the slope deformation accelerated significantly after 15 July 2017. Then it began to decrease. Unlike the case in the former two periods, heavy rainfall occurred again on 24 August 2017, which would accelerate the slope deformation again according to the deformation law mentioned above. On the same day as the Pusa landslide occurred, it rained again and increased the weight of the slope body [15,39]. When stress exceeded the critical locking stress, the landslide eventually occurred at 10:30 local time, 28 August 2017. Furthermore, mining may lead to surface deformation because the stress state was changed within the slope. According to the relevant research [25], the mining tunnels shown in Figure 8 were located 200-300 m beneath the landslide. The mining resumed production in early April 2017. Thereafter, mining was suspended four times before the landslide occurred, and the total closed time was 40 days. Finally, it was shut down on 17 August 2017. It can be seen from Figure 11 that the source area was stable before April 2017. Then it entered into the deformation stage, which indicates that underground coal mining was another triggering factor of the Pusa landslide.

The Failure Mechanism
Based on the aforementioned analysis, the mechanism of the Pusa landslide can be conceptualized in Figure 12 as follows. The original profile of the Pusa landslide is shown in Figure  12a. Figure 12b shows that the underground mining of the coal seams leads to the formation of goafs. The roof of the goaf loses the support from the rock mass below, so it will deform and subside by its gravity, which results in the surface settlement and tension cracks within a certain range. Once it enters the rainy days, the rainwater infiltrated into rock mass through these cracks in Figure 12c, which will increase the hydrostatic pressures, and decrease the strength of the soft rocks in the slope and undermine slope stability [25,40,41]. Besides, the study area is located on the east side on the Yungui Plateau, where the change in the temperature and long-term solarization will aggravate the fragmentation of rock mass, favorable to the rainfall infiltration. The factors of rainfall and temperature change lead to wetting-drying alternation and heating-cooling alternation, which accelerates the rock fragmentation [42]. Finally, the landslide occurs at the source region, based on the local geomorphologic conditions shown in Figure 12d.
Some previous studies have focused on rainfall-induced landslides [14,22], the similarity between these landslides and the Pusa landslide is that the rainfall infiltration will increase the hydrostatic pressures and decrease the strength of the soft rocks in the slope. But the particularity of Pusa landslide was the joint effect of underground mining and rainfall, which accelerated the process of deformation, especially under the karst environment.

Conclusions
The Pusa landslide was fully analyzed with SAR imagery in terms of the boundary and source area delineation, pre-event deformation monitoring, and triggering factors' explanation. First, we manifested that coherence maps and SAR intensity maps generated by ALOS/PALSAR-2 images can be applied to identify the boundary and source area of the landslide. Second, the surface deformation can be uncovered by differential interferograms and the SBAS-InSAR technique by using both L-band ALOS/PALSAR-2 and C-band Sentinel-1A/B datasets. Moreover, the pre-event deformation time series was calculated by both ascending and descending Sentinel-1A/B datasets, which revealed that both underground mining-induced subsidence and unstable slope deformation jointly occurred. Third, two triggering factors, namely rainfall and underground mining, are uncovered by correlation analysis among the pre-event deformation time series, deformation rate variation, cumulative rainfall, and underground mining information. Moreover, the precursory deformation was clearly captured four months before the landslide event. Lastly, the failure mechanism of the Pusa landslide is summarized as four stages, which can be referred to as the rock avalanches in the southwestern of China.
Our research provides a systematic strategy for landslide early-warning in similar areas, such as the unstable slope identification by stacking interferogram method, deformation time-series monitoring by both ascending and descending SAR images. Moreover, the source area and the deposition area can be well measured with InSAR coherence maps and intensity maps for loss assessment. Otherwise, as for the combined effects of natural and anthropogenic factors to the rock avalanches, we suggest prevention measures should be developed in terms of rainfall data collection, real-time surface deformation monitoring, and geotechnical survey, especially caused by underground mining.
Author Contributions: L.C. and C.Z. performed the experiments and produced the results. L.C. drafted the manuscript. C.Z., C.Y., and Y.L. finalized the manuscript. Y.K., H.C., B.L., and A.X. contributed to the discussion of the results. All authors conceived the study and reviewed and approved the manuscript. All authors have read and agreed to the published version of the manuscript.