Next Article in Journal
Sensitivity of Green-Up Date to Meteorological Indicators in Hulun Buir Grasslands of China
Next Article in Special Issue
Coupling Data- and Knowledge-Driven Methods for Landslide Susceptibility Mapping in Human-Modified Environments: A Case Study from Wanzhou County, Three Gorges Reservoir Area, China
Previous Article in Journal
Analysis of Canopy Gaps of Coastal Broadleaf Forest Plantations in Northeast Taiwan Using UAV Lidar and the Weibull Distribution
Previous Article in Special Issue
Transfer Learning for Improving Seismic Building Damage Assessment
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

A Simple Method of Mapping Landslides Runout Zones Considering Kinematic Uncertainties

1
State Key Laboratory of Resources and Environmental Information System, Institute of Geographic Sciences and Natural Resources Research, Chinese Academy of Sciences, Beijing 100101, China
2
University of Chinese Academy of Sciences, Beijing 100049, China
3
National Disaster Reduction Center of China, MEM, Beijing 100124, China
*
Author to whom correspondence should be addressed.
Remote Sens. 2022, 14(3), 668; https://doi.org/10.3390/rs14030668
Submission received: 15 December 2021 / Revised: 25 January 2022 / Accepted: 29 January 2022 / Published: 30 January 2022

Abstract

:
Landslides can be triggered by natural and human activities, threatening the safety of buildings and infrastructures. Mapping potential landslide runout zones are critical for regional risk evaluation. Although remote sensing technology has been widely used to discover unstable areas, an entire landslide runout zone is difficult to identify using these techniques alone. Some simplified methods based on empirical models are used to simulate full-scale movements, but these methods do not consider the kinematic uncertainties caused by random particle collisions in practice. In this paper, we develop a semi-empirical landslide dynamics method considering kinematic uncertainties to solve this problem. The uncertainties caused by the microtopography and anisotropy of the material are expressed by the diffusion angle. Monte Carlo (MC) simulations are adopted to calculate the probability of each cell. Compared with the existing Flow-R model, this method can more accurately and effectively estimate runout zones of the Yigong landslide where random particle collisions are intense. Combining the D-InSAR technique, we evaluate the runout zones in the Jinsha River from June 2019 to December 2020. This result shows that the method is of great significance in early warning and risk mitigation, especially in remote areas. The source area of the landslide and DEM resolution together affect the number of MC simulations required. A landslide with a larger volume requires a larger diffusion angle and more MC simulations.

1. Introduction

Changes in the shear stress and shear strength of slope, which can be triggered by several factors acting alone or together, can cause the slope to become unstable and lead to a landslide [1]. Natural causes include earthquakes [2,3], rainfall [4,5], and repeated freezing and thawing [6]. Human activities, such as deforestation and mining, can also trigger landslides [7]. As a consequence of their unpredictability [8,9], landslides have always caused a large number of casualties and property losses and threatened basic infrastructure [10,11,12,13]. Mapping the runout zone of potentially unstable slopes is critical for regional hazard evaluations and the maintenance of major engineering facilities in alpine countries.
The most significant features in landslide monitoring research are slope deformation and displacement. The identification of deformation areas previously relied on time-consuming and laborious fieldwork [14]. With the development of remote sensing satellite sensors and related technology, the landslide mapping method [15,16,17] has been rapidly developed, which allows us to monitor large and inaccessible mountainous regions. Landslide interpretation is based on shape, size, pattern, color, and texture with characteristic contrast on optical images with different resolutions [18,19]. This method has advantages in identifying regional landslides, especially those with large deformations. However, the quality and resolution of the image limit the scope of its application. The landslide mapping result strongly depends on the experience and skills of the investigator. Satellite-based synthetic aperture radar interferometry (InSAR) provides an effective tool for landslide deformation monitoring from regional to local scales [20,21]. InSAR has been increasingly explored and successfully applied to identify and monitor different types of slope deformation [22]. In contrast with optical remote sensing technology, InSAR can identify small deformations on the ground (cm or less) [23]. However, InSAR can only measure the ground deformation, which is insufficient to reveal the runout zone of a landslide.
Numerical modeling, including the finite element method [24], discrete element method [25], and depth-integrated continuum mechanics method [26], provides an alternative way to identify runout landslide areas and can simulate full-scale movements of an entire landslide body [27,28]. Pure numerical modeling methods have certain advantages in the simulation of single landslides [29,30], but their application in regional landslide simulation is relatively rare due to the difficulty of obtaining parameters and the time-consuming computation [31,32].
To solve the above problems, some simplified methods based on empirical [33,34] or semi-empirical [35,36] approaches have been proposed. The Flow-R, a distributed empirical model for the regional assessment of gravitational hazards [37], has been widely used in landslide susceptibility mapping by scholars due to its low data requirements and openness [38,39]. The propagation area is simulated by frictional laws and flow direction algorithms. The authors developed a new spreading algorithm that avoids over-channelization and thus produces more realistic extents. To facilitate simulations, the conditions or parameters in the actual propagation process must be simplified and abstracted, and sources diffused along a single direction or multiple directions must be determined with a set ratio, which ignores the irregular diffusion caused by particle propagation after random collisions.
The kinematic uncertainties mentioned above are considered probabilistic problems in statistics. In the interest of developing an efficient tool that considers the uncertainties caused by random collisions, we provide an accurate and efficient semi-empirical method for mapping the almost full-scale landslide runout zones. In this method, source areas may be derived from previous studies or by other means. The propagation direction is given by aspect or D-infinity, and the kinematic uncertainties are assessed through Monte Carlo (MC) simulation and the diffusion angle. A complex landslide process is regarded as a probabilistic problem and the particle’s maximum travel distance is calculated by frictional laws in this method. The susceptibility of runout zones is obtained according to the arrival times of the MC simulation. Then, the proposed methodology is applied to two case studies. In the first case, we map the runout zone of the Yigong landslide and compare the result with those from the same type of Flow-R software. In the second case, we interpret source areas based on D-InSAR technology and apply our method for mapping regional landslide runout zones in a section of the Jinsha River.
Specific methods are introduced in the second chapter. The application of this method to individual landslides and regional landslides is presented in the third chapter. Finally, several conclusions are given.

2. Methods

2.1. Computational Approach

In this study, we simulate the propagation process of landslide sources based on landslide dynamics to obtain landslide runout areas. The source area regions and digital elevation model (DEM) data with a certain resolution are important input data. The model deals with the propagation of the landslide downslope, which is a very relevant landslide hazard component, especially for materials displaced by fluidized landslides, but does not correspond to a complete landslide susceptibility assessment. Thus, the assessment method of the source areas is not presented in detail, and it can be customized or derived from existing research and other methods such as D-InSAR technology. Vectorized landslide source regions (shapefile format) are used as inputs. We transfer the vectors into the grids for computation by rasterization. Grids within the source area are assigned a value of 255, and the resolution is consistent with the input DEM. This transfer process is realized by the GDAL package [40]. The simulation of the propagation process is designed in Python language. The main implementation process is shown as follows (Figure 1):
  • The cells within the vector boundary in the landslide source area are marked as source cells, and each cell is taken into a propagation loop.
  • The first source cell is marked as the first active cell. Then, the kinetic energy   E k i n ( i , j ) in the first active cell is calculated by energy conservation theory. If E k i n ( i , j ) > 0, the cell in the propagation direction calculated by MC simulation is marked as the next active cell. The loop continues until E k i n ( i , j ) ≤ 0.
  • Each cell in the source cell list executes step 2 many times, which is determined by preset MC simulation iteration. When the active cell propagates to cell ( i , j ) , 1 is added to T ( i , j ) in the arrival times matrix.
  • Every cell in the source cell list executes steps 2 and 3, and an arrival times matrix is obtained. Finally, the susceptibility map is calculated by the method mentioned in the later section.

2.2. Theoretical Model

In this section, we introduce the specific theoretical models or algorithms of the main steps. First, the whole simulation process is carried out based on the DEM, and two coordinate systems based on grid data are established. Subsequently, we describe the simulation process of landslide runout from three main aspects: kinematic processes, selection of propagation paths, and hazardous zone analysis. As shown in Figure 2, there are four specific steps.
1.
Step 1: Establishment of coordinate systems in grid data
The DEM is a powerful tool to describe regional surface physical characteristics and is a key foundation for simulating dynamic processes. In this paper, DEM data are used to describe the topography and related properties, such as slope and direction, which are important parameters for the calculation of energy changes during movement. Before the simulation, the depressions and flat regions in raw DEMs that prevent the DEM from making unreasonable mistakes that lead to inaccurate slopes and aspects need to be removed. This process is mainly implemented through the Python package Pysheds [41].
We first establish two main coordinate systems (Figure 3), a local coordinate system based on a single cell, and a basic (global) coordinate system. In the local coordinate system, the energy conversion during the landslide propagation process is mainly calculated. Each cell in the DEM is a separate local coordinate system. The basic coordinate system is a two-dimensional coordinate system, mainly used to perform the propagation direction and diffusion algorithm.
The movement trajectory of the simulated material in the source area mainly considers two aspects. The first aspect is the propagation direction of the material in a basic (global) coordinate system, and the second aspect is the movement process of the landslide along the movement direction in a local coordinate system based on a single cell.
2.
Step 2: Minimum energy trajectory calculation
The change in external conditions breaks the stability of the source area of the landslide. The rock and soil mass of the landslide moves along the sliding surface in the direction of low potential energy. In the process of landslide movement, the material potential energy of the landslide body changes into kinetic energy. Friction occurs between the bottom of the landslide body and the slide bed, resulting in internal friction energy, which is released in the form of heat energy, sound energy, and so on. Finally, the velocity potential energy variable is consumed as 0. The whole process is the mutual transformation of gravitational potential energy, kinetic energy, and the internal energy of friction:
E k i n ( i , j ) = E k i n ( i 1 , j 1 ) + Δ E p o t ( i , j ) E f ( i , j )
where E k i n ( i , j ) is the kinetic energy of cell ( i , j ) , E k i n ( i 1 , j 1 ) is the kinetic energy of the previous cell ( i 1 , j 1 ) , Δ E p o t ( i , j ) is the change in potential energy to cell ( i , j ) , and E f ( i , j )   is the energy lost in friction to cell ( i , j ) .
The change in potential energy to cell ( i , j ) is set as the height difference between the current position and the previous position, as shown in the equation:
Δ E p o t ( i , j ) = m g ( H ( i , j ) H ( i 1 , j 1 ) )
where H ( i , j )   is the elevation value of cell ( i , j ) , H ( i 1 , j 1 )   is the elevation value of the previous cell ( i 1 , j 1 ) , m is the mass of a single cell in the source area, and g is the acceleration due to gravity.
The loss of internal friction is not included. The friction energy can be obtained by a simplified friction-limited model (SFLM) [37], which is characterized by the minimum travel angle (or reach angle) [42]:
E f ( i , j ) = m g Δ x   t a n φ
where Δ x   is the increment of horizontal displacement, and t a n φ is the gradient of the energy line [43]. We focused on mapping the landslide runout areas in a simple way, so the specific energy in the propagation process was not obtained. The mass m mentioned in Equations (2) and (3) which is a parameter assumed for the convenience of explaining the algorithm, was not involved in the actual budgeting process.
This approach may result in improbable runout distances in steep catchments due to unrealistic kinetic energy amounts reached during propagation. To keep the energy within reasonable values, a maximum limit can be introduced to ensure that it does not exceed realistic kinetic energy:
E k i n ( i , j ) = m v ( i , j ) 2
E k i n ( i , j ) = m i n {   E k i n ( i , j ) ,   E k i n m a x }
where v ( i , j )   is the velocity value of cell ( i , j ) . The maximum kinetic energy   E k i n m a x   can be estimated by the propagation maximum velocity as the mass m is eliminated from these formulas. The maximum velocity of different types of landslides varies. The maximum speed limit for debris flow is generally 11–35 m/s [44,45]. Other types of landslides may have even greater velocities, especially long-runout landslides. Thus, the speed limit must be set according to the regional situation. In this paper, the maximum speed is generally set at 36 m/s.
3.
Step 3: Selection of propagation path
The calculation of propagation direction is essential to determine the propagation path of the landslide. In addition, the uncertainties caused by the collisions and frictions between particles are also not ignored. In this step, the direction calculation and diffusion calculation are used to implement the selection of the propagation path.
4.
Direction calculation
Under the action of no external force, the debris or rock in the source area will propagate in the direction of the fastest descent surface. There are many algorithms applied to the direction of landslide propagation, such as D8 and D-infinity. In our paper, the D-infinity and aspect propagation direction determination methods are applied to landslide simulation.
5.
D-infinity:
In a 3 × 3 grid window (Figure 4), eight spatial triangles can be formed by the central grid and its 8 adjacent grids, and the slope of each triangle can be calculated. The triangle with the largest slope is selected and taken as the slope of the central grid, and the aspect of the triangle is taken as the propagation direction of the central grid. The specific algorithm is as follows:
For any triangular surface, its slope can be expressed by a vector ( s 1 , s 2 ). The calculation formula for s 1 and s 2 is:
s 1 = H 0 H 1 d 1
s 2 = H 0 H 2 d 2
H i (i = 0,1,2) represents the elevation of the grid, d1 and d2 represent the distance between adjacent grids, and the direction r (ranging from 0° to 360°) of the triangular surface is:
r = a r c t a n ( s 2 s 1 )
6.
Aspect:
Each pixel in the input grid is accessed by a moving 3 × 3 window, and the slope value of each pixel in the centre of the window is calculated using an algorithm that includes the values of eight adjacent pixels. These pixels are identified using the letters e1 through e8, where e0 represents the pixel whose slope direction is currently being calculated. The calculation method of slope aspect is adopted in ArcGIS, and the specific algorithm is as follows:
The rate of change in e0 in the x-direction is calculated by the following algorithm:
  Slope   w e = ( e 8 + 2 e 1 + e 5 ) ( e 7 + 2 e 3 + e 6 ) 8 ×   Cellsize  
The rate of change in e0 in the y-direction is calculated by the following algorithm:
  Slope   s n = ( e 7 + 2 e 4 + e 8 ) ( e 6 + 2 e 2 + e 5 ) 8 ×   Cellsize  
The slope and aspect are calculated by the following algorithm:
Slope   = t a n   Slope   we 2 +   Slope   s n 2
Aspect   =   a r c t a n   ( Slope   s n /   Slope   w e )
In the above formula,   Slope   w e represents the slope in the X-direction, and   Slope   s n represents the slope in the Y-direction. According to the value of the propagation direction, it is divided into 8 different ranges representing different propagation directions, which point to the next propagation grid (Figure 5).
7.
Diffusion calculation
If the propagation process only follows the above propagation direction method, the simulation results will be mainly composed of multiple curves. However, in the actual process of propagation, debris or rock can spread due to collisions and other causes [46]. Each chip unit is an anisotropic discrete block, the kinematics of which are affected by its position, size, and shape and are not easily predicted by deterministic models. According to the above propagation algorithm, these uncertainties cannot be described. Wu et al. [28] established a statistical model of particles by describing the diffusion term and the offset term through differential equations.
In our model, we achieve this by setting the diffusion angle; that is, the landslide does not propagate along a single direction but along a certain range, as shown in Figure 6. MC simulation is used to obtain the random propagation direction within the range of the diffusion angle (Figure 6). Then, the next cell is evaluated according to the value range of cells in each direction shown in Figure 5. After several MC simulations, each cell records the number of times it has been activated.
8.
Step 4: Hazardous zone analysis
After several MC simulations, we obtain a matrix of propagation times, representing the possibility of the landslide impact range, that is, the sensitivity of the landslide. With additional simulation iterations, the arrival times of the highly sensitive grid cells increase, and the difference between them and the low-sensitivity areas becomes much larger. To better display the results, we utilize the method of detecting outliers with boxplots and mark outliers greater than the upper limit of normal as the maximum normal value (see Equation (13)).
T ( i , j ) = m i n { T ( i , j ) , Q U + 1.5 I Q R }
where T ( i , j ) is the arrival times of cell ( i , j ) , Q U is the top one-fourth quartile, and   I Q R is the interquartile range. Then, the arrival times matrix is normalized to express the susceptibility of the landslide (see Equation (14)).
P ( i , j ) = T ( i , j ) T m i n T m a x T m i n
where P ( i , j ) is the susceptibility value of cell ( i , j ) , and the range of the susceptibility value is from 0 to 1. T m a x is the maximum value in the arrival times matrix, and   T m i n is the minimum value in the arrival times matrix.

3. Results

3.1. Case Study 1

In this case, we simulate the runout zone of the Yigong landslide and compare the result with those from the same type of model Flow-R to express the advantages of our method in dealing with the uncertainty of the landslide movement process and to verify the accuracy of the model. The Yigong landslide, which is regarded as a typical long-runout landslide by scholars that occurred in April 2000, in the southern part of the Nyainqentanglha Range in southeastern Tibet (Figure 7) is selected as a case study. The landslide mass traveled from its source area at 5000 m to the sediment area at 2190 m. The horizontal displacement is about 6.7–7 km. It took only 3 min to travel down a horizontal distance of about 8 km at an average speed of about 37–39 m/s. The volume of material in the source area of the Yigong landslide is estimated to be 30–300 Mm2 [47,48]. Based on the interpretation of the Landsat-5 image taken after the landslide, the landslide covered an area of about 12.75 km2. Zhamu Creek, which is approximately 2.7 km long, is the drainage channel of the Yigong landslide. The horizontal distance from the source area to the Yigong riverbed is 10.7 km with a mean longitudinal slope of 31° [48].
The setting range of the landslide source area is based on previous studies [49]. Our method and Flow-R are both used to perform this case on a 30 m global DEM of the Yigong landslide area as acquired from Massflow software case 5 [50]. The D infinity method is chosen to determine the propagation direction, and the MC simulation iteration number is set to 300 in our method. According to previous reports, we set the friction angle to 13° and the maximum speed to 44 m/s [51,52].
The large source material volume and obvious elevation variation make the mass move at high speed, and the collision and entrainment effect is strong in the upper and middle parts of the creek [53]. In our method, we aim to expand the diffusion angle to describe the uncertainties of the block motion.
The Flow-R model contains various algorithms that are accurate for many kinds of gravitational landslides. The model mainly includes three parts: a direction algorithm, a persistence algorithm, and a friction algorithm. The direction algorithm determines the propagation direction of the landslide unit, the persistence algorithm is mainly used to express the repeated inertia behavior in the propagation process, and the friction algorithm determines the maximum propagation distance. We choose three types of direction algorithms including the Holmgren algorithm, the D-infinity algorithm, and the Quinn algorithm [54]. The different weights models are selected as persistence algorithms. The travel angle model and the Perla algorithm [55] are selected as friction algorithms. The detailed parameter settings for the different algorithms are shown in Figure 8. For a description of the Flow-R parameters, see the article by Horton [37]. We obtain a result (the red box in Figure 8) with the default parameter settings to compare the results with different settings. The results obtained with the different direction determination algorithms show small changes. Satisfactory results are not obtained by changing the persistence algorithm. Finally, we choose the Perla algorithm as the friction model. When the friction coefficient mu decreases or the mass-to-drag ratio m/d increases (the recommended value ranges from 103 to 105 [55]), the propagation length increases, but still stops at Zhamu Creek.
The simulation results based on the method in this paper with a 45° diffusion angle are not satisfactory, as with Flow-R. Due to the lack of consideration of uncertainty or smaller diffusion angle in the landslide transmission process, the propagation cells in the simulation process gather along the line of Zhamu Creek. Deposition or a backslope along the river line that causes the runout process to stop.
When the diffusion angle expands to 60° and 70°, the cells travel down the creek to the toes, resulting in a better simulation result especially at 70°. The runout cells start to spread along the two sides of the creek, and the effects of depositions and backslope are weakened or eliminated. However, the results based on our method with a diffusion angle of 70° show that the runout cells cross the actual runout path even at high slopes or erosion/transportation zones. With the increase of diffusion angle, the range of propagation direction angle is extended. The landslide mass may spread over the erosion zone on either side of the creek. Meanwhile, due to the energy loss of internal friction which is not considered in our method, the kinetic energy is overestimated resulting in an amplified runout zone. If the difference between the propagation direction of the next cell and the current propagation direction is greater than 90°, it is marked as a backslope. Accurately judging the propagation direction when the landslide passes through a backslope is difficult. In this method, the runout cells after a backslope are stopped to reduce the possibility of incorrect propagation. Therefore, due to the effect of the backslope on the other side of the Yigong River, the movement stops in the middle of the river, and as a result, the accumulation zone and the surge zone are not completely identified. The unstable region is abstracted as a single grid point by these two methods, which are based on DEM data. The volume change during the moving process, which is an important factor to map the accumulation zone and surge zone, is difficult to simulate.
Fluid mechanics methods [57,58] can be used to simulate the movement process of a landslide, and readily identify its accumulation zone and surge zone. However, these models often require a long calculation process, especially for regional landslide simulations. Existing empirical models rarely consider the uncertainty caused by the random collision of blocks. Therefore, it is necessary to identify the landslide hazard area by simulating the process based on some simplified dynamic theories and empirical parameter models. Our method combines a semi-empirical dynamics model and a statistical model to simulate the runout zone of a landslide and create a susceptibility distribution map. We focus on mapping landslides runout zones and attempt to design a quick and efficient model that requires simple parameters and considers kinematic uncertainties. Therefore, the state of motion can be ignored to speed up the computation.

3.2. Case Study 2

In this case, we apply our method to the identification of regional landslide runout zones in a section of the Jinsha River. D-InSAR stacking technology [59] with Sentinel-1 SAR images is used to identify the location and scope of the unstable slopes. The Sentinel-1 dataset is composed of 48 ascending-pass scenes (Path 99, Frame 1280), VV polarized, acquired between June 2019 and December 2020. The acquisition mode is interferometric wide (IW). The 12.5 m ALOS DEM is used to perform the landslide dynamics simulation process in this study area.
The coverage of Sentinel-1 data selected in this study is nested between Ganzi City, Sichuan Province and Qamdo City, Xizang Province approximately between 97.7° E–100.8° E and 29.6° N–31.3° N (as shown in Figure 9). The study area covers approximately 57,940 km2, where the G318 road, the Jinsha River, and Yalong River and its branches pass. The elevation varies from 5886 m down to 2431 m, bringing about extreme topographic gradients. The region has a high-altitude monsoon in which the climate influences microclimate phenomena, with sufficient rainfall lasting from June to September [60]. This region is a transition zone of tectonic activity from the Qinghai-Tibet Plateau to the Sichuan Basin, which includes multiple fault zones [21]. Long-term tectonic activities have broken rock masses and reduced the stability of slopes [61]. The Baige landslide located in the study area dammed the Jinsha River twice and has received extensive attention in China and worldwide. The identification of landslide runout zones in this area is conducive to the prevention of landslide disasters and corresponding, reductions in casualties and property losses.

3.2.1. Assessment of Landslide Deformation Area by D-InSAR

The D-InSAR analysis results can help identify the main deformation area of the landslide, which will be used as one of the main input data in the following study. The specific process of assessing the deformation area using the D-InSAR technique is as follows: (1) Sentinel-1 SAR data with short temporal and spatial orbital separations are selected to minimize decorrelation. After coregistration, every two images that are adjacent in time are used to generate interferograms, and then high-quality interferograms are selected. (2) The topographic phase is removed using the SRTM DEM (30 m), and Goldstein interferogram filtering [62] for the turbulent atmosphere is performed. (3) The average processing of multiphase interference results is carried out by the stacking tool in GAMMA, which can eliminate the random phase error and enhance the deformation phase. The obtained smooth phase gradients (in radians) allow the establishment of clear boundaries between potentially stable and unstable areas, which are helpful for manual interpretation and evaluation of deformation properties. Examples of the mapped potential slope instabilities are shown in Figure 10. (4) According to a comparison with Google images, the deformation areas are mapped manually. Finally, based on interpretation experience, the deformations unrelated to landslides are removed.
In the process of interpreting the landslide likelihood, considering the landslide type, the following kinds of deformation should be excluded: (1) Glacial surge and glacial front moraine flow deformation and freeze-thaw surface deformation in high and gentle areas; (2) Seasonal deformation caused by the change in moisture content in the loose bottom layer, phase change caused by the multipath refraction of the radar wave by crops, and the change in lake ice; (3) Phase characteristics are caused by surface elevation changes by human factors, such as site levelling and cultivated land reclamation. In addition, the interpretation process of landslides needs to consider the features of landforms, optical remote sensing images, and interference results to identify whether the deformation is an active landslide likelihood. As shown in Figure 10A, there is prominent overall deformation with a clear boundary, and compared with the optical image, a rockfall or slide may occur in this area. In Figure 10B, the deformation pattern is fuzzy and discontinuous, but an active landslide is found by combining the geomorphic features and optical image features, which may be caused by repeated freezing and thawing.

3.2.2. D-InSAR Assessment Results

Figure 11 presents the spatial pattern of D-InSAR stacking deformation, which represents the average displacement value during the period of data acquisition time. Overall, the D-InSAR stacking results are complete and smooth, and the ground displacement region is clear in partial Figure 11i–iii. However, there are partial gaps in low-altitude areas near the Jinsha River, which may be caused mainly by slope-induced shadow and layover problems. High-elevation areas also found gaps caused mainly by vegetation-induced decorrelation effects. The positive and negative deformation in Figure 11 indicates movement towards the east and west, respectively. We find a large deformation zone in Figure 11ii, which occurred in the Xiongba ancient landslide and has been reported by scholars [63]. After human interpretation and comparison with Google images, a total of 39 potential slope deformations are identified. In the process of identification, deformation is caused by glacial movement and artificial reclamation is removed. These deformations occurred between June 2019 and December 2020. Since we focused on analysing the potential runout area of the landslide, there is no specific time of occurrence.
The deformation is mainly concentrated on both sides of the Jinsha River and Yalong River. The total area of deformation is 20.68 km2, of which the largest landslide area is 3.59 km2 and the smallest area is 0.04 km2. The average slope of the deformation area is 25.72°. The deformation occurs in the range of 2700 m–4800 m, and the distribution of each interval is uniform (Figure 12A). The main slope direction of the deformation area is concentrated in the east direction (0–180°) (Figure 12B).

3.2.3. Landslide Runout Modeling

Based on our experience in the simulation process, when the number of MC simulation iterations is more than 50, the change in the simulation results is small or even negligible. Therefore, the number of MC simulations in the whole study area is also set at 50. The aspect-based method is selected for the simulation, the friction coefficient is set at 0.2, the diffusion angle is set at 30°, and the maximum velocity is set at 36 m/s. The landslide runout zones of the entire study area are obtained, and the simulation results (the same areas as shown in Figure 11) are shown in Figure 13. We identify the runout zones of 39 landslides in this study, with a total area of 39.73 km2 and an average danger area of 1.02 km2. The largest runout zone is 4.13 km2, and the smallest is 0.06 km2. In our model, the susceptibility index is obtained from the matrix of runout cell arrival times. The region with a high susceptibility index is concentrated in the gullies because the source cells often spread along the gullies. The susceptibility index on both sides of the creek, especially in the steep slopes and erosion zones, is low because of the lower number of arrival times.
Through a comparison of Google images in different periods and considering topographic factors and image features, the potential runout zones of landslides are drawn using the expert visual interpretation method. To verify the accuracy of the simulation results, we compared the simulation results with the visual interpretation results. We evaluated the model by calculating the precision, recall, and IoU (Intersection over Union) in the three groups of unstable zones in Figure 13. The Precision is the ratio of the True Positive (TP) area to the simulated area from our model. The TP area is the overlap region between the simulation results and the results from the expert visual interpretation. The Recall is the ratio of the TP area to the results from expert visual interpretation. The IoU is the ratio of the intersection area to the union area of the simulation results and the visual interpretation results. As shown in Table 1, the average Precision in the three zones is 0.813, the average Recall is 0.847, and the average IoU is 0.706. The average Precision value indicates that the effective value in our simulation results is at a high level, and the Recall results show that the simulated results cover most of the visual interpretation results. IoU results greater than 0.5 indicate that the gap between the simulation results and the visual interpretation results is acceptable. In summary, our results can approximately reflect the potential Runout region of the unstable region.
The lower right-hand corner of Figure 13i shows a simulation result that travels along two channels, which differs from the visual interpretation results. Landslides travel along two channels. The reason for this phenomenon may be related to the parameter setting in the simulation process; thus, we focus on the influence of parameter setting on the result in the discussion.
Identifying the regional landslide runout zone in the study area reveals that the method of combining the D-InSAR technique and landslide dynamics, which considers energy conservation and uncertain diffusion, can effectively represent the uncertainty of the landslide propagation process and substantially improve the reliability of regional landslide runout zone results. During the past few decades, scholars have carried out many analyses on the failure mechanism, process, and sensitivity of unstable landslide areas based on InSAR technology [64,65]. Through field investigation, optical image discrimination, and geomorphic feature reference, scholars can obtain an accurate unstable landslide area. However, In-SAR technology can only identify the area where ground deformation occurs and cannot accurately identify landslide danger zones, which should include transportation zones, deposition zones, and surge zones. Scholars have identified unstable areas based on InSAR technology, described landslide movement by numerical simulation methods, and obtained landslide hazardous zones. Landslide simulations based on Newton’s second law of rigid bodies are widely used, such as Flow-R and Rockfall analysis [66]. Due to the superimposed influence of multiple factors and irregular collisions during the landslide propagation process, there is some uncertainty in the landslide propagation process [67]. Therefore, based on the traditional landslide propagation mechanism and considering the uncertainty in the process of landslide propagation, it is meaningful to accurately identify the danger area after the landslide.

4. Discussion

The diffusion angle, friction angle, and DEM spatial resolution are the main parameters in our method. To illustrate the effect of some choices in the parameters, a sensitivity analysis of the input parameters is carried out. The friction angle, a fundamental physical constant, is always obtained by apparent friction angle or physical experiments, so we focus on the influence of the input DEM spatial resolution and diffusion angle to obtain the results in this paper. In addition, MC simulation iterations related to the completeness of the runout zone are also discussed. Finally, we provide a suggested range of parameter values. The Mangkam deformation zone on the bank of the Jinsha River in the city of Mangkam (Figure 14), which may threaten community safety and the capacity of the Jinsha River, is a typical case used to discuss the effect of parameter variation on the results. The deformation area is mapped on the D-InSAR stacking result in case study 2.

4.1. Influence of the Elevation Model Resolution

We select DEMs (descriptions of which are shown in Table 2) with resolutions of 12.5 m, 30 m, and 90 m to illustrate the impact of data resolution and the role of the modification of MC simulation iterations. Other parameters are set as shown in Figure 15. The method of determining the propagation direction is based on aspect direction.
As shown in Figure 15, the results of 10 MC simulations are relatively incomplete, which is most obvious on the 90 m DEM. In these simulations, the propagation at every resolution is satisfied with more than 50 simulation iterations. Although all the simulation results in the 12.5 m and 30 m resolutions DEM are propagated along these two channels that are located on the two sides of the Mangkam likelihood, the former simulation results are more detailed. Due to the roughness of the 90 m resolution DEM, the simulated results show three propagation channels. Finer resolutions can simulate an appropriate runout zone, while rougher resolutions artificially enlarge the extent. This suggests that higher resolution grids are needed for researchers who focus on risk analysis in a single unstable region. According to previous studies [68,69], a 10 m resolution DEM is a reasonable compromise between improved resolution and data volume in the simulation of geomorphological and hydrological processes, and a larger resolution result does not bring more important information.
In addition, we find that the simulation results of 12.5 m and 30 m grids do not change significantly after 50 cycles, except for some changes in the boundary. In general, DEMs with higher resolutions yield more runout paths, and more simulation iterations are needed to obtain all the propagation results. However, only 10 simulation iterations in the 12.5 m DEM are needed to obtain the basic propagation range, and 50 simulation iterations are needed in the 90 m DEM. The number of grid elements in the landslide likelihood increases the times of their propagation to the danger area, which is equivalent to increasing the number of simulation iterations. Hence, when the likelihood of a landslide is larger or the grid is more refined (which means much more grids in the likelihood), final results can be obtained with fewer simulation iterations.

4.2. Influence of the Diffusion Angle

The migration process of landslides is mainly affected by convective terms and diffusion terms in the depth-integrated continuum model. Based on the random walk mobility, we simplify the diffusion term into the diffusion angle to obtain an accurate runout zone. We perform MC simulations at different times with three different diffusion angles of 15°, 30°, and 45° based on the DEM grid of 12.5 m to illustrate the influence of diffusion angles.
Figure 16 shows the results with different diffusion angles and different simulation iteration numbers. Only one channel appears in the simulation results of a diffusion angle of 15°, while two channels appear in the simulation results of 30° and 45°. As the diffusion angle increases, the widths of the two channels also increase. With additional MC simulation iterations, the simulation results of the 15° diffusion angle do not change significantly, the results of the 30° diffusion angle do not change significantly after 50 cycles, and the simulation results of the 45° diffusion angle do change. This phenomenon occurs because the larger diffusion angle increases the number of transmission paths of lattice points in the unstable region. Thus, for a larger diffusion angle, more simulation iterations are needed.
As the landslide volume increases and the block size decreases, the opportunities for impact among the blocks increases as the diffusion effect during the propagation process is enlarged [28,70]. We find that the 15° diffusion angle result is consistent with the visual interpretation in case study 2. D-InSAR results show that the Mangkam deformation is at a low level, which indicates that the volume is small, so we need to set a smaller diffusion angle. In case 1, due to the large volume of the Yigong landslide, it is necessary to set a larger diffusion angle. The selection of diffusion angle should be determined according to different types of landslides. In regional landslide hazard simulations, we need to set a compromise value to ensure a more accurate range.

4.3. The Suggested Parameter Setting Range

We propose a simple landslide runout zone mapping method with fewer input parameters. However, parameterization is still an issue that needs to be discussed. When modeling a phenomenon for which information is limited, the subjectivity of choosing parameters should be reduced. In this section, we provide a method for forecasting values (as shown in Table 3), which consequently allows us to hypothesize reliable runout zones without having to resort to back analysis processes.
The diffusion angle should be selected according to different types of landslides. The diffusion angle can affect the shape of the landslide danger zone. In previous studies, the most common landslide shape measure is the relationship between the length and width of a landslide, which can be used to describe the diffusion angle [71,72]. Taylor et al. [73] presented an ellipse as a reasonable model for >80% of landslide shapes in different geomorphic settings and proved the feasibility of using ellipticity and the length-to-width ratio to describe the shape of landslides. The elliptic length-to-width ratio ranged from 1.2 to 15.1. The diffusion angle can be calculated by determining the arctangent of the elliptic length-to-width ratio which ranges from 4 to 39.55°. However, due to uncertainties, further research is required to estimate the diffusion angle using the landslide shape. In landslide simulations with a high degree of kinematic uncertainties (e.g., long-runout landslides), the diffusion angle should be set to 45°–90°. In other cases, the angle is usually set to 15°–45°, which is a compromise range. A smaller diffusion angle may lead to incomplete simulation results due to the strong effect of channelization [37], while a larger diffusion angle may lead to simulation results beyond the real range.
Choosing parameters for SFLM models is not easy, but some simple methods can be used to define friction angles. Jaboyedoff [74] proposed a method to obtain the friction angle using easily available parameters. The minimum friction angle of debris flows has been set to 11° in previous studies. We generally recommend a range of friction angles from 10° to 30°.
Following the upper limit of 70 m/s for “extremely rapid” landslides published in 1995 by the international geological union landslide working group [75], the recommended range for the maximum speed is 11–70 m/s. In the early years, scholars studied the prediction and assessment of landslide speeds [35,76]. Intuitively, greater landslide driving forces are derived from steeper slopes. Researchers have shown that the speed of a landslide seems to be related to the slope angle according to a power function [77]. Although this relationship is subject to further discussion, it can serve as an important reference factor.

5. Conclusions

Identifying hazardous landslide runout zones in certain areas, especially near cities, rivers, and roads, is necessary to predict landslide disasters and reduce economic losses and casualties. This paper proposes a semi-empirical landslide runout zones mapping method considering uncertainty in propagation for assessing landslide runout zones. Instead of considering propagation in a certain direction or multiple directions, we used the diffusion angle to describe the landslide propagation process. MC simulation was applied to express the probability of propagation to each position. The simulation method was used in two case studies. The following major conclusions can be drawn from this work:
  • Compared with other models of the same type, our method obtains a better Yigong landslide simulation result. The diffusion angle concept provides the possibility to express the uncertainty of landslide simulation in grid data and solves the false stops problem due to obstacles in the process of propagation.
  • Although there are partial gaps in certain areas, the D-InSAR stacking method can effectively identify unstable zones in the case study 2 area. Thirty-nine deformation zones were successfully detected. Compared with our results and expert interpretation results, we found that the combined method is simple and effective for regional landslide hazardous area identification.
  • We discussed the influence of DEM resolution and diffusion angle on the results and their relationship with the number of MC simulations required. Among the three kinds of DEMs tested, the simulation effect of 12.5 m resolution is the optimal choice, and a higher resolution DEM is needed if high-quality simulation results are needed. A larger landslide volume requires larger diffusion angles because of the increased opportunities for collisions between blocks. Under the same conditions, due to the influence of the size of the source area, a DEM with higher resolution does not necessarily need additional MC simulation iterations. The larger diffusion angle requires more MC simulation iterations.
Although the proposed method brought some improvements, some issues must be remarked. The first is that the parameters of the landslide propagation model are difficult to determine accurately because of its movement process type, terrain influence, provenance type, and triggering factors. Thus, empirical model parameters were applied to calculate the friction angle, diffusion angle, and maximum propagation velocity. Although these equations have been commonly used in previous studies, accuracy tests need to be further improved to develop theoretical models. In addition, we suggest that the compromise value should be selected in the process of regional landslide hazard area simulation. The second is that the surge zone and accumulation area were not involved in this landslide simulation model. These two aspects will be further discussed in the future.

Author Contributions

Conceptualization, Y.W.; methodology, J.L., X.G. and Y.W.; software, J.L.; validation, J.L.; resources, X.G.; writing—original draft preparation, J.L.; writing—review and editing, J.L. and X.Z.; visualization, X.G.; supervision, X.G. and X.Z.; project administration, Y.W. All authors have read and agreed to the published version of the manuscript.

Funding

This research was supported by the Strategic Priority Research Program (Class A) of the Chinese Academy of Sciences (No. XDA23090503); by the National Key Research and Development Plan of China (No. YS2018YFGH000001).

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Data sharing is not applicable to this article.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Cui, P.; Guo, C.X.; Zhou, J.W.; Hao, M.H.; Xu, F.G. The mechanisms behind shallow failures in slopes comprised of landslide deposits. Eng. Geol. 2014, 180, 34–44. [Google Scholar] [CrossRef]
  2. Cui, P.; Chen, X.Q.; Zhu, Y.Y.; Su, F.H.; Wei, F.Q.; Han, Y.S.; Liu, H.J.; Zhuang, J.Q. The Wenchuan Earthquake (May 12, 2008), Sichuan Province, China, and resulting geohazards. Nat. Hazard. 2011, 56, 19–36. [Google Scholar] [CrossRef]
  3. Yin, Y.P.; Wang, F.W.; Sun, P. Landslide hazards triggered by the 2008 Wenchuan earthquake, Sichuan, China. Landslides 2009, 6, 139–152. [Google Scholar] [CrossRef]
  4. Hu, X.; Burgmann, R.; Schulz, W.H.; Fielding, E.J. Four-dimensional surface motions of the Slumgullion landslide and quantification of hydrometeorological forcing. Nat. Commun. 2020, 11, 1–9. [Google Scholar] [CrossRef] [PubMed]
  5. Kirschbaum, D.; Adler, R.; Adler, D.; Peters-Lidard, C.; Huffman, G. Global Distribution of Extreme Precipitation and High-Impact Landslides in 2010 Relative to Previous Years. J. Hydrometeorol. 2012, 13, 1536–1551. [Google Scholar] [CrossRef] [Green Version]
  6. Zhou, J.W.; Cui, P.; Hao, M.H. Comprehensive analyses of the initiation and entrainment processes of the 2000 Yigong catastrophic landslide in Tibet, China. Landslides 2016, 13, 39–54. [Google Scholar] [CrossRef]
  7. Watkinson, I.M.; Hall, R. Impact of communal irrigation on the 2018 Palu earthquake-triggered landslides. Nat. Geosci. 2019, 12, 940–945. [Google Scholar] [CrossRef]
  8. Li, W.; Chen, Y.; Liu, F.; Yang, H.F.; Liu, J.L.; Fu, B.H. Chain-Style Landslide Hazardous Process: Constraints From Seismic Signals Analysis of the 2017 Xinmo Landslide, SW China. J. Geophys. Res.: Solid Earth 2019, 124, 2025–2037. [Google Scholar] [CrossRef]
  9. Pudasaini, S.P.; Miller, S.A. The hypermobility of huge landslides and avalanches. Eng. Geol. 2013, 157, 124–132. [Google Scholar] [CrossRef]
  10. Petley, D. Global patterns of loss of life from landslides. Geology 2012, 40, 927–930. [Google Scholar] [CrossRef]
  11. Haque, U.; Blum, P.; da Silva, P.F.; Andersen, P.; Pilz, J.; Chalov, S.R.; Malet, J.P.; Auflic, M.J.; Andres, N.; Poyiadji, E.; et al. Fatal landslides in Europe. Landslides 2016, 13, 1545–1554. [Google Scholar] [CrossRef]
  12. Huang, R.Q. Some catastrophic landslides since the twentieth century in the southwest of China. Landslides 2009, 6, 69–81. [Google Scholar] [CrossRef]
  13. Haque, U.; da Silva, P.F.; Devoli, G.; Pilz, J.; Zhao, B.X.; Khaloua, A.; Wilopo, W.; Andersen, P.; Lu, P.; Lee, J.; et al. The human cost of global warming: Deadly landslides and their triggers (1995–2014). Sci. Total Environ. 2019, 682, 673–684. [Google Scholar] [CrossRef] [PubMed]
  14. Guzzetti, F.; Mondini, A.C.; Cardinali, M.; Fiorucci, F.; Santangelo, M.; Chang, K.T. Landslide inventory maps: New tools for an old problem. Earth Sci. Rev. 2012, 112, 42–66. [Google Scholar] [CrossRef] [Green Version]
  15. Dini, B.; Manconi, A.; Loew, S. Investigation of slope instabilities in NW Bhutan as derived from systematic DInSAR analyses. Eng. Geol. 2019, 259, 105111. [Google Scholar] [CrossRef]
  16. Moragues, S.; Lenzano, M.G.; Moreiras, S.; Lo Vecchio, A.; Lannutti, E.; Lenzano, L. Slope instability analysis in South Patagonia applying multivariate and bivariate techniques on Landsat images during 2001–2015 period. Catena 2019, 174, 339–352. [Google Scholar] [CrossRef]
  17. Yang, W.T.; Wang, Y.J.; Wang, Y.Q.; Ma, C.; Ma, Y.H. Retrospective deformation of the Baige landslide using optical remote sensing images. Landslides 2020, 17, 659–668. [Google Scholar] [CrossRef]
  18. Soeters, R.; Van Westen, C.J. Slope instability recognition, analysis and zonation. Landslides Investig. Mitig. 1996, 247, 129–177. [Google Scholar]
  19. Coe, J.A.; Bessette-Kirton, E.K.; Geertsema, M. Increasing rock-avalanche size and mobility in Glacier Bay National Park and Preserve, Alaska detected from 1984 to 2016 Landsat imagery. Landslides 2018, 15, 393–407. [Google Scholar] [CrossRef] [Green Version]
  20. Lauknes, T.R.; Shanker, A.P.; Dehls, J.F.; Zebker, H.A.; Henderson, I.H.C.; Larsen, Y. Detailed rockslide mapping in northern Norway with small baseline and persistent scatterer interferometric SAR time series methods. Remote Sens. Environ. 2010, 114, 2097–2109. [Google Scholar] [CrossRef]
  21. Xu, Q.; Zheng, G.; Li, W.; He, C.; Dong, X.; Guo, C.; Feng, W. Study on successive landslide damming events of jinsha river in baige village on Octorber 11 and November 3,2018. J. Eng. Geol. 2018, 26, 1534–1551. [Google Scholar]
  22. Barboux, C.; Delaloye, R.; Lambiel, C. Inventorying slope movements in an Alpine environment using DInSAR. Earth Surf. Processes Landf. 2014, 39, 2087–2099. [Google Scholar] [CrossRef] [Green Version]
  23. Bejar-Pizarro, M.; Notti, D.; Mateos, R.M.; Ezquerro, P.; Centolanza, G.; Herrera, G.; Bru, G.; Sanabria, M.; Solari, L.; Duro, J.; et al. Mapping Vulnerable Urban Areas Affected by Slow-Moving Landslides Using Sentinel-1 InSAR Data. Remote Sens. 2017, 9, 876. [Google Scholar] [CrossRef] [Green Version]
  24. Zienkiewicz, O.C.; Taylor, R.L.; Taylor, R.L.; Taylor, R.L. The Finite Element Method: Solid Mechanics; Butterworth-Heinemann: Oxford, UK, 2000; Volume 2. [Google Scholar]
  25. Lo, C.-M.; Lee, C.-F.; Chou, H.-T.; Lin, M.-L. Landslide at Su-Hua highway 115.9 k triggered by typhoon Megi in Taiwan. Landslides 2014, 11, 293–304. [Google Scholar] [CrossRef]
  26. Wang, F.; Sassa, K. Landslide simulation by a geotechnical model combined with a model for apparent friction change. Phys. Chemistry Earth Parts A/B/C 2010, 35, 149–161. [Google Scholar] [CrossRef]
  27. Ma, P.F.; Cui, Y.F.; Wang, W.X.; Lin, H.; Zhang, Y.Z. Coupling InSAR and numerical modeling for characterizing landslide movements under complex loads in urbanized hillslopes. Landslides 2021, 18, 1611–1623. [Google Scholar] [CrossRef]
  28. Wu, Y.M.; Lan, H.X. Debris Flow Analyst (DA): A debris flow model considering kinematic uncertainties and using a GIS platform. Eng. Geol. 2020, 279, 105877. [Google Scholar] [CrossRef]
  29. Ouyang, C.J.; He, S.M.; Xu, Q.; Luo, Y.; Zhang, W.C. A MacCormack-TVD finite difference method to simulate the mass flow in mountainous terrain with variable computational domain. Comput. Geosci. 2013, 52, 1–10. [Google Scholar] [CrossRef]
  30. Mergili, M.; Fischer, J.T.; Krenn, J.; Pudasaini, S.P. r.avaflow v1, an advanced open-source computational framework for the propagation and interaction of two-phase mass flows. Geosci. Model Dev. 2017, 10, 553–569. [Google Scholar] [CrossRef] [Green Version]
  31. Ouyang, C.J.; Zhou, K.Q.; Xu, Q.; Yin, J.H.; Peng, D.L.; Wang, D.P.; Li, W.L. Dynamic analysis and numerical modeling of the 2015 catastrophic landslide of the construction waste landfill at Guangming, Shenzhen, China. Landslides 2017, 14, 705–718. [Google Scholar] [CrossRef]
  32. Carrara, A.; Crosta, G.; Frattini, P. Comparing models of debris-flow susceptibility in the alpine environment. Geomorphology 2008, 94, 353–378. [Google Scholar] [CrossRef]
  33. Iverson, R.M.; Schilling, S.P.; Vallance, J.W. Objective delineation of lahar-inundation hazard zones. Geol. Soc. Am. Bull. 1998, 110, 972–984. [Google Scholar] [CrossRef]
  34. Hungr, O. A model for the runout analysis of rapid flow slides, debris flows, and avalanches. Can. Geotech. J. 1995, 32, 610–623. [Google Scholar] [CrossRef]
  35. Scheidegger, A.E. On the prediction of the reach and velocity of catastrophic landslides. Rock Mech. 1973, 5, 231–236. [Google Scholar] [CrossRef]
  36. Hunter, G.; Fell, R. Travel distance angle for "rapid" landslides in constructed and natural soil slopes. Can. Geotech. J. 2003, 40, 1123–1141. [Google Scholar] [CrossRef]
  37. Horton, P.; Jaboyedoff, M.; Rudaz, B.; Zimmermann, M. Flow-R, a model for susceptibility mapping of debris flows and other gravitational hazards at a regional scale. Nat. Hazards Earth Syst. Sci. 2013, 13, 869–885. [Google Scholar] [CrossRef] [Green Version]
  38. Zou, Q.; Jiang, H.; Cui, P.; Zhou, B.; Jiang, Y.; Qin, M.Y.; Liu, Y.G.; Li, C. A new approach to assess landslide susceptibility based on slope failure mechanisms. Catena 2021, 204, 105388. [Google Scholar] [CrossRef]
  39. Nguyen, B.Q.V.; Kim, Y.T. Regional-scale landslide risk assessment on Mt. Umyeon using risk index estimation. Landslides 2021, 18, 2547–2564. [Google Scholar] [CrossRef]
  40. GDAL. Available online: https://pypi.org/project/GDAL/2.4.0/ (accessed on 15 December 2021).
  41. Bartos, M. Pysheds: Simple and Fast Watershed Delineation in Python. Available online: https://github.com/mdbartos/pysheds (accessed on 11 November 2021).
  42. Corominas, J. The angle of reach as a mobility index for small and large landslides. Can. Geotech. J. 1996, 33, 260–271. [Google Scholar] [CrossRef]
  43. Li, L.P.; Lan, H.X.; Strom, A. Automatic generation of landslide profile for complementing landslide inventory. Geomat. Nat Haz. Risk. 2020, 11, 1000–1030. [Google Scholar] [CrossRef]
  44. Kang, H.; Yun-Tae, K. Parameter Analysis of Flow-R Model for Physical Vulnerability Assessment of Debris Flow Disaster in Regional Scale. J. Korean Soc. Hazard Mitig. 2017, 17, 233–245. [Google Scholar] [CrossRef]
  45. Rickenmann, D.; Zimmermann, M. The 1987 Debris flows in switzerland-documentation and analysis. Geomorphology 1993, 8, 175–189. [Google Scholar] [CrossRef]
  46. Gnyawali, K.R.; Xing, A.G.; Zhuang, Y. Dynamic analysis of the multi-staged ice-rock debris avalanche in the Langtang valley triggered by the 2015 Gorkha earthquake, Nepal. Eng. Geol. 2020, 265, 105440. [Google Scholar] [CrossRef]
  47. Xu, Q.; Shang, Y.J.; van Asch, T.; Wang, S.T.; Zhang, Z.Y.; Dong, X.J. Observations from the large, rapid Yigong rock slide-debris avalanche, southeast Tibet. Can. Geotech. J. 2012, 49, 589–606. [Google Scholar] [CrossRef]
  48. Delaney, K.B.; Evans, S.G. The 2000 Yigong landslide (Tibetan Plateau), rockslide-dammed lake and outburst flood: Review, remote sensing analysis, and process modelling. Geomorphology 2015, 246, 377–393. [Google Scholar] [CrossRef]
  49. Zhou, G.G.D.; Roque, P.J.C.; Xie, Y.X.; Song, D.R.; Zou, Q.; Chen, H.Y. Numerical study on the evolution process of a geohazards chain resulting from the Yigong landslide. Landslides 2020, 17, 2563–2576. [Google Scholar] [CrossRef]
  50. Massflow. Available online: http://www.massflow-software.com/ (accessed on 5 November 2021).
  51. Shen, Y.; Shen, W.; Li, T.; Guo, J.; Lei, Y. Numerical investigation on impact of bed entrainment to the mobility of rapid flow-like landslides. J. Eng. Geol. 2019, 27, 1405–1414. [Google Scholar]
  52. Huang, R.Q.; Xu, Q. Catastrophic landslides in China; Science Press: Beijing, China, 2008. [Google Scholar]
  53. Shang, Y.J.; Yang, Z.F.; Li, L.H.; Liu, D.; Liao, Q.L.; Wang, Y.C. A super-large landslide in Tibet in 2000: Background, occurrence, disaster, and origin. Geomorphology 2003, 54, 225–243. [Google Scholar] [CrossRef]
  54. Quinn, P.; Beven, K.; Chevallier, P.; Planchon, O. The prediction of hillslope flow paths for distributed hydrological modeling using digital terrain models. Hydrol. Processes 1991, 5, 59–79. [Google Scholar] [CrossRef]
  55. Perla, R.; Cheng, T.; McClung, D.M. A two–parameter model of snow–avalanche motion. J. Glaciol. 1980, 26, 197–207. [Google Scholar] [CrossRef] [Green Version]
  56. Flow-R. Available online: https://www.terranum.ch/en/products/flow-r/ (accessed on 19 July 2021).
  57. Wu, Y.M.; Lan, H.X. Landslide Analysta landslide propagation model considering block size heterogeneity. Landslides 2019, 16, 1107–1120. [Google Scholar] [CrossRef]
  58. Pastor, M.; Blanc, T.; Haddad, B.; Petrone, S.; Morles, M.S.; Drempetic, V.; Issler, D.; Crosta, G.B.; Cascini, L.; Sorbino, G.; et al. Application of a SPH depth-integrated model to landslide run-out analysis. Landslides 2014, 11, 793–812. [Google Scholar] [CrossRef] [Green Version]
  59. Yao, X.; Deng, J.; Liu, X.; Zhou, Z.; Yao, J.; Dai, F.; Ren, K.; Li, L. Primary Recognition of Active Landslides and Development Rule Analysis for Pan Three-river-parallel Territory of Tibet Plateau. Adv. Eng. Sci. 2020, 52, 16–37. [Google Scholar]
  60. Fan, X.M.; Xu, Q.; Alonso-Rodriguez, A.; Subramanian, S.S.; Li, W.L.; Zheng, G.; Dong, X.J.; Huang, R.Q. Successive landsliding and damming of the Jinsha River in eastern Tibet, China: Prime investigation, early warning, and emergency response. Landslides 2019, 16, 1003–1020. [Google Scholar] [CrossRef]
  61. Ding, C.; Feng, G.C.; Liao, M.S.; Tao, P.J.; Zhang, L.; Xu, Q. Displacement history and potential triggering factors of Baige landslides, China revealed by optical imagery time series. Remote Sens. Environ. 2021, 254. [Google Scholar] [CrossRef]
  62. Goldstein, R.M.; Werner, C.L. Radar interferogram filtering for geophysical applications. Geophys. Res. Lett. 1998, 25, 4035–4038. [Google Scholar] [CrossRef] [Green Version]
  63. Guo, C.; Yan, Y.; Zhang, Y.; Zhang, X.; Zheng, Y.; Li, X.; Yang, Z.; Wu, R. Study on the Creep-Sliding Mechanism of the Giant Xiongba Ancient Landslide Based on the SBAS-InSAR Method, Tibetan Plateau, China. Remote Sens. 2021, 13, 3365. [Google Scholar] [CrossRef]
  64. Li, L.J.; Yao, X.; Yao, J.M.; Zhou, Z.K.; Feng, X.; Liu, X.H. Analysis of deformation characteristics for a reservoir landslide before and after impoundment by multiple D-InSAR observations at Jinshajiang River, China. Nat. Hazard. 2019, 98, 719–733. [Google Scholar] [CrossRef]
  65. Schlogel, R.; Doubre, C.; Malet, J.P.; Masson, F. Landslide deformation monitoring with ALOS/PALSAR imagery: A D-InSAR geomorphological interpretation method. Geomorphology 2015, 231, 314–330. [Google Scholar] [CrossRef]
  66. Lan, H.X.; Martin, C.D.; Lim, C.H. RockFall analyst: A GIS extension for three-dimensional and spatially distributed rockfall hazard modeling. Comput. Geosci. 2007, 33, 262–279. [Google Scholar] [CrossRef]
  67. Eidsvig, U.M.K.; Papathoma-Kohle, M.; Du, J.; Glade, T.; Vangelsten, B.V. Quantification of model uncertainty in debris flow vulnerability assessment. Eng. Geol. 2014, 181, 15–26. [Google Scholar] [CrossRef]
  68. Zhang, W.H.; Montgomery, D.R. DIgital elevation model grid size, landscape representation, and hydrologic simulations. Water Resour. Res. 1994, 30, 1019–1028. [Google Scholar] [CrossRef]
  69. Quinn, P.; Beven, K.; Lamb, R. The in (a/tan/β) index: How to calculate it and how to use it within the topmodel framework. Hydrol. Processes 1995, 9, 161–182. [Google Scholar] [CrossRef]
  70. Okura, Y.; Kitahara, H.; Sammori, T.; Kawanami, A. The effects of rockfall volume on runout distance. Eng. Geol. 2000, 58, 109–124. [Google Scholar] [CrossRef]
  71. Hovius, N.; Stark, C.P.; Allen, P.A. Sediment flux from a mountain belt derived by landslide mapping. Geology 1997, 25, 231–234. [Google Scholar] [CrossRef] [Green Version]
  72. Milledge, D.G.; Bellugi, D.; McKean, J.A.; Densmore, A.L.; Dietrich, W.E. A multidimensional stability model for predicting shallow landslide size and shape across landscapes. J. Geophys. Res.-Earth Surf. 2014, 119, 2481–2504. [Google Scholar] [CrossRef] [Green Version]
  73. Taylor, F.E.; Malamud, B.D.; Witt, A.; Guzzetti, F. Landslide shape, ellipticity and length-to-width ratios. Earth Surf. Processes Landforms 2018, 43, 3164–3189. [Google Scholar] [CrossRef]
  74. Jaboyedoff, M.; Rudaz, B.; Horton, P. Concepts and parameterisation of Perla and FLM model using Flow-R for debris flow. In Proceedings of the 5th Canadian Conference on Geotechnique and Natural Hazards, Kelowna, BC, Canada, 15–17 May 2011. [Google Scholar]
  75. IUGS–International Working Group. A suggested method for describing the rate of movement of a landslide. Bull. Int. Assoc. Eng. Geol. 1995, 52, 75–78. [Google Scholar] [CrossRef]
  76. Hungr, O.; Corominas, J.; Eberhardt, E. Estimating landslide motion mechanism, travel distance and velocity. In Landslide Risk Management; CRC Press: Boca Raton, FL, USA, 2005; pp. 109–138. [Google Scholar]
  77. Wen, B.P.; Wang, S.J.; Wang, E.Z.; Zhang, J.M. Characteristics of rapid giant landslides in China. Landslides 2004, 1, 247–261. [Google Scholar] [CrossRef]
Figure 1. Workflow of the computational approach.
Figure 1. Workflow of the computational approach.
Remotesensing 14 00668 g001
Figure 2. Schematic diagram of the main steps of this method.
Figure 2. Schematic diagram of the main steps of this method.
Remotesensing 14 00668 g002
Figure 3. Definition of the coordinate system and cell plane.
Figure 3. Definition of the coordinate system and cell plane.
Remotesensing 14 00668 g003
Figure 4. Definition of the coordinate system and cell plane.
Figure 4. Definition of the coordinate system and cell plane.
Remotesensing 14 00668 g004
Figure 5. A 3 × 3 window and definition of the different propagation directions.
Figure 5. A 3 × 3 window and definition of the different propagation directions.
Remotesensing 14 00668 g005
Figure 6. Schematic diagram of diffusion angle and MC simulation. Grey cells are the cells where landslides propagate.
Figure 6. Schematic diagram of diffusion angle and MC simulation. Grey cells are the cells where landslides propagate.
Remotesensing 14 00668 g006
Figure 7. Yigong landslide map. The background image was taken by Landsat5-TM on 12 May 2000.
Figure 7. Yigong landslide map. The background image was taken by Landsat5-TM on 12 May 2000.
Remotesensing 14 00668 g007
Figure 8. Yigong landslide danger zone simulation results for different methods. The value of the susceptibility index represents the possibility of landslide occurrence in every grid. In the results of Flow-R, the red box shows the result with the default parameter setting, and the black boxes show the results with the modified parameters, which are marked with black text. The level of landslide susceptibility index represents the possibility that a landslide will occur in the region. Flow-R software version is 2.0.7 [56].
Figure 8. Yigong landslide danger zone simulation results for different methods. The value of the susceptibility index represents the possibility of landslide occurrence in every grid. In the results of Flow-R, the red box shows the result with the default parameter setting, and the black boxes show the results with the modified parameters, which are marked with black text. The level of landslide susceptibility index represents the possibility that a landslide will occur in the region. Flow-R software version is 2.0.7 [56].
Remotesensing 14 00668 g008
Figure 9. Study site. The grey polygon represents the coverage of the Sentinel-1 data.
Figure 9. Study site. The grey polygon represents the coverage of the Sentinel-1 data.
Remotesensing 14 00668 g009
Figure 10. Two examples of objects mapped on single wrapped interferograms. (A) Observation of a rock slope deformation. (B) Displacement on a rock glacier. Within the white dashed line, the area with displacements.
Figure 10. Two examples of objects mapped on single wrapped interferograms. (A) Observation of a rock slope deformation. (B) Displacement on a rock glacier. Within the white dashed line, the area with displacements.
Remotesensing 14 00668 g010
Figure 11. D-InSAR Stacking results. The red points represent the deformation area. (iiii) show the results in detail. Red polygons represent unstable regions.
Figure 11. D-InSAR Stacking results. The red points represent the deformation area. (iiii) show the results in detail. Red polygons represent unstable regions.
Remotesensing 14 00668 g011
Figure 12. Variation in (A) elevations and (B) aspects of deformations in the study area. Aspects are shown on a 360° rose diagram. Numbers on circles are counts within that range.
Figure 12. Variation in (A) elevations and (B) aspects of deformations in the study area. Aspects are shown on a 360° rose diagram. Numbers on circles are counts within that range.
Remotesensing 14 00668 g012
Figure 13. Landslide danger zone simulation results. The black polygon regions represent the runout areas of visual interpretation.
Figure 13. Landslide danger zone simulation results. The black polygon regions represent the runout areas of visual interpretation.
Remotesensing 14 00668 g013
Figure 14. D-InSAR result and optical image of Mangkam deformation zone. The blue polygon represents the Mangkam likelihood. The red polygon represents the propagation area by our method (diffusion angle = 30°, tangent of friction angle = 0.3).
Figure 14. D-InSAR result and optical image of Mangkam deformation zone. The blue polygon represents the Mangkam likelihood. The red polygon represents the propagation area by our method (diffusion angle = 30°, tangent of friction angle = 0.3).
Remotesensing 14 00668 g014
Figure 15. Effect of variation of DEM resolution and MC simulation iterations on the Mangkam deformation runout zone (diffusion angle = 30°, tangent of friction angle = 0.3).
Figure 15. Effect of variation of DEM resolution and MC simulation iterations on the Mangkam deformation runout zone (diffusion angle = 30°, tangent of friction angle = 0.3).
Remotesensing 14 00668 g015
Figure 16. Effect of variation of the diffusion angle on the Mangkam deformation runout zone (tangent of friction angle = 0.3, DEM resolution = 12.5 m).
Figure 16. Effect of variation of the diffusion angle on the Mangkam deformation runout zone (tangent of friction angle = 0.3, DEM resolution = 12.5 m).
Remotesensing 14 00668 g016
Table 1. The evaluation indicator results in different zones of Figure 13.
Table 1. The evaluation indicator results in different zones of Figure 13.
Area_IDLA_IDIoUPrecisionRecall
i10.5450.5550.969
20.5850.7730.706
30.6830.8530.774
ii40.7750.8780.868
50.7970.9280.849
60.7170.7670.917
70.8980.9730.921
80.6140.7010.833
90.7310.9060.791
iii100.7600.8350.894
110.7460.8590.851
120.6160.7300.798
Table 2. The DEMs used in this study.
Table 2. The DEMs used in this study.
DEMResolutionAcquisition Time
ALOS12.5 mJune 2008
ASTER GDEM30 m2000–2009
SRTM90 m2000
Table 3. The setting range of input parameters.
Table 3. The setting range of input parameters.
Parameter NameSetting Range
Diffusion Angle15°–90°
Friction Angle10°–30°
Maximum Speed11–70 m/s
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Liu, J.; Wu, Y.; Gao, X.; Zhang, X. A Simple Method of Mapping Landslides Runout Zones Considering Kinematic Uncertainties. Remote Sens. 2022, 14, 668. https://doi.org/10.3390/rs14030668

AMA Style

Liu J, Wu Y, Gao X, Zhang X. A Simple Method of Mapping Landslides Runout Zones Considering Kinematic Uncertainties. Remote Sensing. 2022; 14(3):668. https://doi.org/10.3390/rs14030668

Chicago/Turabian Style

Liu, Jia, Yuming Wu, Xing Gao, and Xuehua Zhang. 2022. "A Simple Method of Mapping Landslides Runout Zones Considering Kinematic Uncertainties" Remote Sensing 14, no. 3: 668. https://doi.org/10.3390/rs14030668

APA Style

Liu, J., Wu, Y., Gao, X., & Zhang, X. (2022). A Simple Method of Mapping Landslides Runout Zones Considering Kinematic Uncertainties. Remote Sensing, 14(3), 668. https://doi.org/10.3390/rs14030668

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop