Three-dimensional slope stability analysis using laser scanning and numerical simulation

ABSTRACT Detailed surface topography is important when analyzing the stability of slopes. Recent advances in new technologies such as interferometric synthetic aperture radar and light detection and ranging have allowed us to obtain high-precision profiles of distant landscapes and objects including detailed slope information for three-dimensional (3D) slope stability analysis. However, techniques of reconstructing 3D numerical models from scanned data of slope geometry have not been well investigated or tested. This paper proposes a comprehensive approach that integrates laser scanning and finite element method for slope stability analysis, particularly failure prediction under precipitation scenarios. The methodology is applied to a slope in the Wenchuan earthquake-stricken mountainous region that failed in 2013, triggered by severe rainfall. The modelling results show that the surface sampling resolution can affect the prediction accuracy of the potential failure zones. When constructing the slope model, the selected surface grid should be fine enough to capture the important topographic features of the slope while minimizing computation demand. This paper confirms that the proposed method can be successfully used to identify the potential failure zones of the investigated slope under severe rainfall conditions.


Introduction
Effective numerical models for landslide stability analysis require a good set of input parameters. Accurate topology data are essential for a reliable landslide stability model; however, this requires considerable effort involving various technologies and methods such as field surveys, global positioning system (GPS) mapping, and digital elevation model (DEM) data. Field surveys of landslide sites or slopes can be challenging if they are located in areas with poor accessibility arising from steep topography or slippery and loose land cover. Moreover, areas that experienced recent earthquakes (such as the mountainous areas struck by the 2008 Wenchuan earthquake) can be very dangerous, because many of the slopes that need to be assessed for hazard level and possible mitigation measures are at high risk of sliding.
In the last decade, new technologies such as interferometric synthetic aperture radar (InSAR) and light detection and ranging (LiDAR) have been intensively studied and used in topographic surveys; the results confirmed these methods as effective tools for acquiring detailed topographic information. Using the digital topography data retrieved by InSAR and LiDAR with relatively high resolution (millimetre to metre), a wide range of geoscience and geoengineering applications were developed, from real-time monitoring to numerical modelling . Although the techniques and scheme designs of InSAR and LiDAR are very different, they share the same principle of physics that the travel time of a signal reflects the distance it travelled. Both InSAR and LiDAR can be used to measure changes in distant objects with high resolution (of mm scale). The different applications (terrestrial, airborne, and mobile) make InSAR and LiDAR suitable for location-specific measurements as well as large-scale morphology scanning. For example, the slow movement of slopes can be monitored by using ground-based InSAR while spotlight interferometry data can be used to estimate the movement rate of the landslides (Leva et al. 2003;Milillo et al. 2014). Similarly, a terrestrial laser scanner (TLS) was used to generate displacement maps (Teza et al. 2007), and the time-series of the TLS point clouds enabled researchers to estimate the movement progression of the investigated landslides (Kasperski et al. 2010). The use of three-dimensional (3D) laser scanning provides accurate measurements of landslides, and the results can be used to estimate earthwork volume (Du & Teng 2007). Based on LiDAR-extracted DEMs of the study area from different time periods, landslides can be identified, and the landslides and non-landslide information can be used to identify and characterize the important factors that affected the landslides .
Integrated approaches combining multiple sources and data acquisition technologies to generate models that characterize both geomorphological and subsurface structures were also explored (Merritt et al. 2014). For example, coal-mining subsidence can be monitored using a combination of 3D TLS and GPS technology (Zhou et al. 2014). High-resolution terrestrial, airborne, and satellite photogrammetry, as a complementary technology, was used to estimate landslide volumes (Martha et al. 2010) and to detect changes in landslides and slopes (Tsutsui et al. 2007;Gonzalez-Diez et al. 2014). 3D models were developed by Ozbay et al. (2015) using two DSLR digital cameras and a laser range finder. Brideau et al. (2012) combined terrestrial photogrammetry and airborne LiDAR data for stability analysis of landslides.
Although the new technologies (InSAR, LiDAR, high-resolution photogrammetry, etc.) provide effective and diversified support in landslide and slope reconstruction, researchers have not yet taken full advantage of the data quality available to improve landslide modelling and simulation and understand the mechanisms that govern landslides, rockfall, and debris flow . Various methods have been used in stability analysis of slopes. These include the finite element method (FEM), discrete element method, limit equilibrium method, finite difference method (FDM), and various combinations of these methods (Wang et al. 2007;Kanungo et al. 2013;Zhang et al. 2014;Ozbay et al. 2015). These numerical methods have recently been incorporated with new technologies to generate suitable 2D/3D models for analysis. For example, LiDAR technology was used to generate a 3D model, upon which 2D FEMs were built to estimate the safety factor at specific sites (Hu et al. 2012). TLS technology was applied to reconstruct the geometry of slopes in Italy, and the generated 3D model was used to simulate the landslide and rockfall processes (Marsella et al. 2015). However, because of the complexity of the problem, especially when rainfall infiltration is of interest, big uncertainties still exist in the slope stability analysis, which makes the failure estimation difficult (Zhang et al. 2011).
Analysis based on 3D modelling is advantageous over 2D modelling because it provides much more detailed topology; moreover, Vassallo et al. (2015) found that the results generated by the 3D model may differ significantly from those of the 2D model. The 3D geometry of a slope can now be obtained from InSAR and LiDAR data; however, processing the large data-set (i.e. the laser cloud points) is challenging. First, the detailed geometry data-set usually contains some redundant and false information, for example, trees and temporary ground cover, which may lead to distorted models. Models with distorted geometry, even in a small, very localized area, can often generate errors in the analysis and misleading results. Second, the creation of 3D models that include detailed geometry may generate a very large number of elements; this will create large errors in most numerical methods (i.e. FEM, FDM) which have limits on the number of elements they can handle, especially when rainfall infiltration is considered. Therefore, careful processing of 3D LiDAR or InSAR data-sets is necessary to eliminate noise data that may cause distortion, and to determine an appropriate geometry model that balances precision with computational capacity.
In this study, we propose a method that integrates 3D laser scanning and 3D FEM analysis for slope stability analysis, particularly failure prediction under scenarios of extreme precipitation. The main steps of the integrated approach are (1) acquisition of a 3D point cloud, (2) reconstruction of the geometry model, (3) creation of a 3D FEM model, and (4) safety analysis and failure prediction. The model is applied in a case study of a failed slope in 2013, triggered by severe rainfall in the Wenchuan earthquake-stricken mountainous region. The investigated slope experienced high-position landslides and is at risk of undergoing deeper large-scale sliding during severe rainfall. The proposed model was tested using different geometry precisions to determine the effect of the precision on the results. In this paper we demonstrate the techniques that integrate 3D cloud points with TSL and 3D FEM analysis and discuss the practical issues that arise during data post-processing and model conversion.
The rest of this paper is structured as follows. Section 2 describes the study site and introduces the overall methodology of the integrated approachobtaining the 3D point cloud data and the 3D FEM modelling and analysis. In Section 3 we apply the proposed methodology to the investigated slope, and the results of the numerical analysis are presented in Section 4, followed by the Conclusions.

The Hongxi landslide
The 2008 Wenchuan earthquake and the subsequent severe rainfall triggered a large number of landslides in the disaster-affected areas. Because of the unstable mountainous environment after the earthquake, numerous slopes located along rivers and highways are currently at high risk of emerging and expanding sliding, which poses a big challenge for post-earthquake reconstruction and recovery. Therefore, prompt assessment of the stability of these slopes is essential before engineering and ecological measures can be taken in these regions. However, obtaining the profiles of the slopes can be difficult because of the complex mountainous topography and changing terrains after the earthquake. The large landslides induced by the earthquake generated a large amount of loose material (debris) that makes geological surveying extremely difficult and dangerous. Many landslides are still very active even without significant man-made (e.g. farming and mining) or natural (e.g. severe rainfall and seismic motion) disturbances. Moreover, traditional surveying and measurement techniques for obtaining the profiles of a large number of slopes can be time-consuming and involves high labour costs. An automated and easy-to-implement technique that integrates slope information acquisition and slope modelling/analysis is urgently required to identify high-risk slopes.
On the night of 11 July 2013, Pingwu County experienced extreme rainfall which lasted 12 h. This severe rainfall triggered a large number of landslides and debris flow in Pingwu, disrupting transport across the county. The village of Hongxi was a post-earthquake resettlement site located on the north bank of the Hongxi river. More than 40 villagers reconstructed their houses in the village, which was completely destroyed by the rainfall-induced landslides, debris flow, and flooding. During the July 11 severe rainfall, the steep mountain slope near the south bank of the river became unstable and a deep landslide formed at an elevation of about 100 m. The landslides generated over 30,000 m 3 of loose soil and rocks that plunged into the river at high speed, and the loose material and water continued flowing onto the north bank, burying and destroying 80% of the houses there. After the disaster, the village of Hongxi was abandoned and the surviving villagers were resettled again in other villages. Bare soil and rock with very limited vegetation were exposed on the mountain near the south bank after the landslide; this site poses a high risk of further landslides, which may disrupt road transport along the south bank and the water flow of the Hongxi river. Figure 1 shows the geographic location of the investigated slope. Figure 2 presents the overview of the slope (Google Earth) before and after the landslide and a photograph taken during a recent field survey.

Methodology
Based on a geological survey conducted in the study area, the Hongxi landslide has a high potential of further sliding. Every summer this area experiences seasonal precipitation, particularly in July and August, raising concern for the stability of this slope in future heavy rain scenarios. A detailed 3D FEM model to determine the conditions of slope failure under rainfall is therefore needed. The methodology of this research is outlined in Figure 3 and described below.

Acquisition of slope configuration
A 3D laser scanner was used to obtain point cloud data for the 3D representations of the investigated slope. The 3D laser scanner records a large number of points on the surface of the measured object and outputs the data as point cloud data. A point cloud is a set of vertices in a 3D coordinate system. The high accuracy of the point cloud data reflects the accurate profile of the slope; however, it also generates some unwanted points, for example, the trees and temporary ground cover, which clearly are not part of the slope profile. Therefore, careful post-processing of the point cloud data is needed to eliminate these points. Another reason for processing the points cloud is that the farther objects contain fewer points than the closer ones; hence, we use an interpolation algorithm in the low-resolution areas to generate an evenly gridded 3D surface covering the whole investigated slope. Furthermore, very high-resolution slope topography data will require extensive numerical computation and will have difficulty converging; therefore, a resampling of the high-resolution point cloud data is needed to create a data-set of lower resolution which still contains the necessary slope profile information. The resampled points are then converted to a triangular-meshed network surface to create the numerical model, and the 3D smooth surface is input into the FEM software (PLAXIS 3D) to construct the 3D slope geometry model. The detailed steps in the post-processing of the data are summarized as follows:  Step 1: Remove manually the obvious abnormal points (i.e. electric wire, tall trees, etc.) from the obtained point cloud data.
Step 2: Resample the high-resolution point cloud data by using spatial sampling method with a given spacing interval (i.e. 80 cm).
Step 3: Extract the coordinates of the resampled point cloud data and import them into MATLAB; Step 4: Create the boundaries of grids in x-y (horizontal) plane and define the grid size (i.e. 1.5 m, 2.0 m, 3.0 m, etc.).
Step 5: Apply interpolation (i.e. linear or triangle interpolation) to find every z (vertical) value corresponding to the grids in x-y plane. This creates a new set of spatial points representing the original point cloud.
Step 6: Import the new set of spatial points into AutoCAD, and then connect the surface points that are closest together to create triangulation lines. The triangulation lines form the triangles that make up the surface triangulation-meshed network.
After the six steps, a smooth surface that approximates the original point cloud data is generated according to the defined grid boundaries and size. The generated smooth surface can be read and further analyzed in PLAXIS.

Testing of soil properties
The soil parameters are crucial for the slope stability analysis and can be derived by tests such as X-ray diffraction tests, compression tests, triaxial tests, and hydraulic conductivity tests. The X-ray diffraction test determines the mineral composition of the soil samples collected at different locations. If the mineral composition is significantly different, more soil samples should be collected, and the soil difference should also be modelled. The compression tests, triaxial tests, and hydraulic conductivity tests are performed to determine the physical properties of the soil, including the Young's modulus, cohesion, friction angle, and hydraulic conductivity. The soil properties are then assigned to the FEM model.

Generation of precipitation scenarios
Historical precipitation records from the study site can be analyzed to generate the precipitation scenarios. Data recorded by rain gauges or a meteorological station located close to the study site are preferred. Brief intensive rainfall and long-lasting rainfall can threaten slope stability. The historical maximum daily rainfall and maximum prolonged rainfall data were extracted from 60-year precipitation records and used as the simulated extreme rainfall events in the FEM model for slope failure prediction.
The numerical analysis was conducted based on the slope configuration data, the soil properties, and the simulated precipitation scenarios; the change of factor of safety (FS) and the potential failure zone of the slope was examined.
3. Slope geometry and 3D model 3.1. 3D scanning and data acquisition The Trimble GX system was used to scan the surface of the Hongxi slope in the study area. The scanner is based on time-of-flight principle and has a laser source that emits light pulses with a wavelength of 532 nm. The distance is determined using the time it takes for the laser pulse to travel between the point of emission and the surface being scanned. The scanner can measure distances up to 350 m, with a scanning speed of 5000 points per second. The field of view is 360 horizontally and 60 vertically. The single point accuracy is 12 mm at 100 m range and the vertical and horizontal angular widths are 0.0014 and 0.0012 , respectively.
The scanning was performed with 30-cm intervals from the top to the bottom of the slope in a clockwise direction. The results of the scanning contain 70,000 cloud points with detailed geometric information of the slope (Figure 4). The scanned horizontal area is 299 m £ 144 m with a maximum height of 110 m. The angle of the sloping surface varies from 30 to 60 .

Reconstruction of slope geometry
The scanned cloud data were used to create the numerical model. The surface generated directly from the point cloud data has irregularities which complicate the numerical analysis (Figure 4(c)). Therefore, we post-processed the original point cloud to create a smooth surface. An ASCII file was exported containing information about the geometric coordinates of each point X, Y, and Z. The extracted data were then imported using the MATLAB software and the interpolated function was used to fit a surface formulated as Z 0 = F(X 0 ,Y 0 ) to the scattered data X and Y, where X 0 and Y 0 are defined between the maximum and minimum values of X and Y. To determine the appropriate interpolation interval, different intervals (1.5, 2, 3, 4, and 5 m) were selected. Too large a point interval would cause information loss while too small a point interval would create a surface with many small curved surfaces and requires long calculation times in the numerical analysis. Figure 5 shows the reconstructed surfaces with an interpolation interval of 1.5 m.

3D FEM model
The numerical model of the slope was then created in PLAXIS which is a geotechnical finite element software. Since the reconstructed surface generated by MATLAB cannot be directly read by PLAXIS, the interpolated point data were then imported to AutoCAD to form continuous triangles using a triangulated irregular network mesh as shown in Figure 6. The continuous triangles were then  exported into PLAXIS to create the 3D FEM models. As shown in Figure 7, the 3D FEM models generated using various surface mesh sizes (intervals of 1.5-5 m) present different geometry details. In the high-resolution model (with 1.5-m point interval), details of the gullies running from the top to the toe of the slope can be observed; these gullies were formed as debris flow gullies during a severe rainfall in July 2013. The model with 2-m point intervals is similar to the one with 1.5-m intervals; however, as the interpolation intervals increase (3, 4, and 5 m), the details of the gully geometry diminish. It should be noted that PLAXIS may face difficulty generating mesh with too much geometry details although the details are important for accurate slope reconstruction. The model with 1.0-m point intervals was used in PLAXIS but it failed to converge when the 3D slope data from our study site were applied. Therefore, when constructing a 3D FEM model from the acquired geometry point cloud, we need to find a balance between more detailed geometry and computational capacity.
The slope model consists of a 20-m-thick soil layer and the slope bedrock. The soil depth of the slope is referred to geological survey (Ding et al. 2009;Hou et al. 2011) and consultation with local geological agents. The interaction between soil and rock is considered as fully coupled in this study. The numerical model is meshed with wedge elements. At the bottom of the finite-element mesh z = 0 plane, the displacements are set to zero in the three directions x, y, and z. Figure 8 shows the hydraulic boundary conditions of the model. The sides and bottom of the slope were set as non-flow boundaries. Rainfall infiltration was simulated by applying a unit flux equal to the rainfall intensities on the sloping surface. The bedrock was assumed to be an impermeable layer. The water table was set as a straight line 2 m from the ground surface at the interface between the bedrock and soil and 0 at the toe of the slope.

Soil properties
The soil was modelled as a homogeneous material using the Mohr-Coulomb failure criterion. The shear strength of the soil t in the Mohr-Coulomb model is given by where s n is the effective normal stress, ' is the angle of internal friction of the soil, and C is the cohesion force of the soil. To obtain the basic parameters for the Mohr-Coulomb model and other soil properties, comprehensive tests (X-ray diffraction tests, compression tests, triaxial tests, and hydraulic conductivity tests) were performed. Nine soil samples (three at the bottom, three in the middle, and three at the top of the slope) were collected from a depth of around 1 m. These samples were analyzed with X-ray diffraction tests. The test results from the different soil samples were similar, indicating that the soil was homogeneous. The soil at the landslide zone analyzed from nine samples is 47%-62% quartz, 15%-24% illite, 8%-12% chlorite, 6%-8% plagioclase, 4%-6% potassium feldspar, 3%-8% dolomite, and 2%-6% montmorillonite.
In the triaxial tests (Figure 9(a)), the soil samples were subjected to vertical stresses while horizontal stresses were applied to the specimens through the surrounding fluid pressures. The soil specimens were 300 mm in height and 150 mm in diameter. The tests were conducted according to the test methods of soils for highway engineering, China (JTG E40-2007, 2007). Once the angle of internal friction and the cohesion of the soil were determined, the shear strength and the Mohr circles (Figure 9(b)) could be calculated. The physical properties used in the numerical model are summarized in Table 1.

Rainfall data
The rainfall pattern parameters are crucial for the study of rainfall-induced slope failure. The rainfall pattern (rainfall intensity, variation, duration, etc.), is highly variable, depending on the geographic  location, climatic condition, and season. In this study, precipitation time-series from historical records were used to track the precipitation variations in Pingwu County. The rainfall in Pingwu is characterized by short and intense or less intense and prolonged rainfall throughout the year, while only prolonged rainfall occurs during the summer, from May to September. The largest rainfall volume in a 6-day period found in the historical records from day 1 (259.5 mm), day 2 (72.3 mm), day 3 (0.4 mm), day 4 (9.1 mm), day 5 (5.8 mm), and day 6 (20.6 mm) was used for the slope stability analysis.

FEM analysis results
The slope stability can be assessed by calculating the FS. In PLAXIS, the FS of the slope was calculated using the strength reduction method. In the strength reduction method, the friction angle and cohesion are reduced incrementally by the same factor to determine the safety factor as shown in Equation (2): where C and ' are the cohesion and the friction angle of the soil, respectively, and C r and ' r are the reduced cohesion and friction angle, respectively, which are controlled by the by the total multiplier P Msf . This parameter is increased in a step-by-step procedure until failure occurs. The FS is then determined as the value of P Msf at failure. The initial (without rainfall) FS of the slope was calculated. The total displacement increment of the slope obtained from the numerical analysis is shown in Figure 10 for five models with different sampling intervals (1.5-5.0 m). The displacements do not represent a physical meaning but give an indication of the potential failure mechanism. Different results were obtained from models with different slope geometry details. In the model with the finest slope geometry (interpolation interval of 1.5 m), the large displacements are in the lower left side of the slope and the potential failure zone is on the left side of the debris flow gully. The displacement pattern of the model with a 2-m interpolation interval is similar to the 1.5-m interval model. However, the potential failure zone changes as the resolution becomes coarser from 3 m to 5 m. The model with 3-m intervals has a large and flat potential failure zone that crosses the gully; hence, the full geometry has already been averaged out in the 3-m model. The potential failure zone moves to the bottom left corner of the slope in the 4-m and 5-m models, implying that in the models with smoother surface features the gully geometry has disappeared. During the field survey, the debris flow gully was clearly observed; the left part of the gully slope appeared unstable with loose soil and rocks. This important topographic feature of the slope surface was lost in the coarser-resolution models. The initial FS also varies with model geometry, ranging from 1.722 to 1.808. The 1.5-m model gives the minimum FS and the 5-m model gives the maximum value.
The 1.5-m model was adopted for slope stability analysis under various rainfall scenarios. Figure 11 shows a 6-day rainfall scenario that represents the largest total precipitation over a 6-day period. The FS of the investigated slope drops as the rain continues, from an initial value of 1.722 to 1.425 on the sixth day. This is due to the increase of pore water pressure during precipitation which leads to a decrease of effective stress of the drained soil and further decreases the shear strength of the soil according to the effective stress principle. The failure surfaces estimated for this scenario by our model for day 1 to day 6 are shown in Figure 12, indicating the highest failure potential is on the left side of the debris flow gully. The potential failure surface extends from the top to the toe of Figure 10. Total displacement increment of the slope for the five models with sampling intervals of 1.5-5.0 m (The total displacement increment does not represent realistic displacement and only shows the potential failure mechanism).
Note: The stress state development in a structure is regarded as an incremental process and a global iterative procedure is required to achieve both equilibrium conditions and constitutive relation. The incremental displacement is the summation of sub-increments from all iterations in the current increment step. the slope, representing a landslide that cuts through the whole vertical side of the slope, as in the landslide that occurred in 2013. This result confirms that the proposed approach integrating both 3D laser scanning and 3D FEM simulation can successfully estimate the safety of the investigated slope and identify the potential failure surface. In the simulation, only an extreme 6-day rainfall scenario is input into the model for the stability analysis, while in reality the slope may have experienced additional precipitation before the modelled scenario, in which case the soil of the slope may already be partially saturated. Therefore, an extreme rainfall scenario in a partially saturated soil may easily lead to slope failure (FS = 1). Although only a case study of one slope is analyzed, the results indicate that many slopes that have already experienced sliding during the 2008 Wenchuan earthquake or in subsequent severe rainfall events may still be unstable. Engineering measures should be taken to monitor the stability of those slopes and implement mitigation measures to reduce disaster risks in the earthquake-stricken regions.

Conclusions
In this study, we present a comprehensive framework for a 3D slope stability analysis in a severe precipitation scenario by integrating laser scanning and FEM. First, a 3D laser scanner was used to acquire point cloud data containing detailed geometric information of the distant slope. Postprocessing of the original point cloud dataset was performed to create a smooth surface with different grid sizes from 1.5 m to 5 m. The surface was then converted into a triangular-mesh surface in AutoCAD and used to construct the finite element model in PLAXIS. Second, a series of laboratory tests (X-ray diffraction tests, compression tests, triaxial tests, and hydraulic conductivity tests) were performed to obtain the physical properties of the slope soil field samples. The mechanical soil properties obtained in the tests were used in the FE model for more realistic simulation and analysis. Third, historical precipitation data were analyzed to generate representative extreme rainfall scenarios; these scenarios were applied to estimate the slope failure potential and the critical condition. The paper takes a large slope as a case study to test the proposed approach. The slope is located in the Wenchuan earthquake-stricken area where numerous landslides were triggered by seismic forces or subsequent severe rainfall. The investigated slope experienced a high-elevation landslide that destroyed the adjacent Hongxi village during an extreme rainfall event in the summer of 2013; however, it still holds a high risk of further sliding and requires detailed analysis. The safety factor and the potential failure surface of the slope under an extreme rainfall scenario were analyzed.
The major conclusions obtained from this study are summarized below.
(1) We confirmed that the proposed approach can be successfully used to obtain detailed 3D geometry information of distant slopes and to generate a computable FE model for analysis. This method can be used to estimate the safety of slopes and identify the potential failure surfaces, providing insights into the physical process of slope failure and tools to examine the safety of specific sites.
(2) The precision of the modelled slope geometry influences the slope stability analysis results. A geometry capturing the main characteristics of the slope profile is required for an accurate slope stability analysis; the balance between detailed geometry and computation cost should also be considered. (3) The results of the case study analysis indicate that under prolonged rainfall, the slope will approach failure state. Engineering measures should be taken to monitor and provide means of mitigating the stability limitations of such slopes to reduce disaster risks in earthquakestricken regions.