Giant landslide displacement analysis using a point cloud set conflict technique: a case in Xishancun landslide, Sichuan, China

ABSTRACT Landslides, threatening millions of human lives, are geological phenomena on earth, occurred frequently. An increasing number of techniques are being used to monitor landslide deformation. Among these, Light Detection and Ranging (LiDAR) stands out for its high efficiency and accuracy in displacement detection, particularly for giant landslides. In this work, we collected two temporal datasets of terrain laser scanning and proposed a flowchart for giant landslide displacement analysis using the point cloud set conflict（PCSC） technique. First, the terrestrial points were obtained by performing registration and off-terrain point filtering. Second, the landslide displacement field was acquired using the proposed method based on its surface roughness. The displacement results from our established methodological system are comparable with the ones of Interferometric Synthetic Aperture Radar (InSAR)-derived deformations. The differences estimated from two systems are at the centimetre level. Cross-analysis on the trigger factor with landslide occurred mechanism could be achieved based on the results as well. Therefore, this work provides a novel system to analyse the displacement of a giant landslide in the future study.


Introduction
Landslides, geological hazardous phenomena occurring frequently worldwide, have been threatening millions of human lives and various types of infrastructures (Guzzetti et al. 2012). Monitoring towards it is particularly important for understanding the characteristics and evolutionary trends (Huang et al. 2015). Light Detection and Ranging (LiDAR) is a rapidly developing technique for landslide monitoring. Compared with traditional in-situ surveying techniques, such as Global Positioning System (GPS) (Travelletti et al. 2012;Malet, Maquaire, and Calais 2002), or remote sensing techniques, such as Interferometry Synthetic Aperture Radar (InSAR) (Herrera et al. 2009;Greif and Vlcko 2012), multi-temporal LiDAR stands out for its low cost and high effectiveness as millions of points could be obtained with only several surveys (Tapete et al. 2015;Szabó et al. 2016). Because of the development of full waveform laser scanning, multi-temporal LiDAR can be used to extract terrain points even in some vegetation-dense areas, particularly mountainous regions (Ferraz, Mallet, and Chehata 2016). Therefore, surface change detection of landslides based on LiDAR is becoming prevalent in landslide movement studies (Jaboyedoff et al. 2012;Ghuffar et al. 2013).
Methods for change detection of LiDAR are still under development. Comparisons of pre-and post-event digital terrain models are initially carried out on the LiDAR datasets (Dewitte et al. 2008). Models with detailed information could be acquired due to abundant point data available from LiDAR; thus, minor changes could be detected after comparisons (Petschko, Bell, and Glade 2016;Pradhan and Abdulwahid 2017). However, object identification in this method still requires extensive labour to construct the entire model, which leads to a high price and consumes large amounts of time particularly for giant landslides (Hsiao et al. 2004). On the other hand, direct three-dimensional comparisons between two point clouds provide a non-time-consuming topographic change detecting (Barnhart and Crosby 2013). Theoretically, for each point of the referenced point cloud, its distance to the nearest neighbour being compared is computed. Based on this theory, such methods as Cloud-to-Cloud (Girardeau-Montaut et al. 2005) and Cloud-to-Mesh (Cignoni, Rocchini and Scopigno 1998) are proposed to measure the deformation. But they are either limited in situations with highly dense points or only work well on flat surfaces (Monserrat and Crosetto 2008;Young et al. 2010), which is not suitable for rough natural objects (Travelletti, Malet, and Delacourt 2014). Therefore, in this study, a method using the Point Cloud Set Conflict (PCSC) is used to detect the surface changes of a giant landslide. Surface roughness is considered in PCSC in order to estimate different objects. Objects in a different temporal point cloud are considered as sets and surface changes are regarded as distances of the sets conflicting.

Methodology
As abovementioned, terrain laser scanning techniques have become an important approach for analysing deformation of a landslide. However, there still exist some problems in using two temporal point clouds. For example, rough natural surfaces would bring uncertainties to uniform distribution of studied points, which leads to a relatively lower accuracy of results. In this paper, we proposed a new system with PCSC technique to acquire and analyse landslide deformation.

PCSC technique
The fundamental theory of PCSC technique is described in Figure 1. Point Cloud Set a (PCS a ) and Point Cloud Set b (PCS b ) are two point cloud sets, representing two sets of a local object (depicted by a dark triangle) yet observed at different time. Here, PCS a is set as a referenced point cloud and PCS b is set as a compared point cloud. The object's location has changed from a(a x , a y ) to b(b x , b y ) during this time interval. Hence, the changing distance should be calculated precisely by a function: dist a; b ð Þ where the direction of it is denoted as in Figure 1 and magnitudes of it are denoted according to Equation (1).
However, as can be seen from Figure 1, there are always observation errors in point clouds, namely σ a and σ b in two results, respectively, hence the precise position of the local object could not be achieved. In practice, we would introduce a certain calculation rule to project the point cloud sets to virtual points, as Equation (2) denoted: In Equation (2), H (•) is the calculating rule to project the points to a virtual point and it can be extended in Figure (2): PCS a and PCS b are two observed point cloud results and their actual distance is calculated based on virtual points. To find the visual points and calculate the distance, four steps should be conducted: • A semi-diameter M/2 is set at every point in the referenced point cloud (PCS a ) and points (grey squares) within its reach are selected to fit a local plane (a black circle). • The normal (the black arrow) of the plane is constructed toward the compared point cloud (PCS b ), and a circular truncated cone is made with another semidiameter of m/2. • Points from each point cloud set have projections along this normal vector. Their average positions are recorded as two virtual points: a i and b i (dark solid dots). Thus, the distance between two virtual points is calculated as z i ¼ dist a i ; b i ð Þ according to Equation (1). • Iterate procedures above on another point and record z i values. The change distance between two point cloud sets is calculated as the average distance between virtual points through Equation (3).

Parameters determination in PCSC
Determining diameters M and m is a key step in PCSC procedure. As studied before, the surface roughness may lead to errors in the point cloud change detection (Heritage and Milan 2009;Hodge 2010;Schürch et al. 2011): if M or m is of similar scales as the point cloud roughness, the normal vector would fluctuate, which results in incorrect estimation on point cloud set conflicting. The roughness can be defined as in Equation (4).
D i is the distance between a point and a locally fitted plane. N is the number of the points included. r is the roughness indicator of the point cloud. Roughness is calculated as r a; i and r b;i , so diameters M and m can be empirically determined as in Equation (5). 3. Case study and LiDAR datasets 3.1. Geographic location of Xishancun landslide Xishancun landslide, located in the south-west of the Chinese Mainland, is selected as the study case in Figure 3(a). It is located on the northern bank of Zagunao river in Lixian County in Figure 3(b) and is only 22 km away from Wenchuan county which is well known for Wenchuan Earthquake in 2008. This regions are considered to be hazardous according to previous studies (Liu et al. 2013;Liu et al. 2016). As shown in Figure 3(c), Xishancun landslide is a giant soil accumulation body forming a 'U' shape. Landslide boundaries are defined by the topography and geomorphology, i.e. cracks, scarps, surface deformation, and bedrock (Xu et al. 2016). Facing south-west, it lies with a leading-edge elevation H 1 being 1,500 m and a trailing edge elevation H 2 being 2,900 m, the relative relief is 1,400 m. Its maximum length L 1 is 4,000 m, and its maximum width L 2 is 1,700 m. The general thickness is 55 m, and the earth volume is approximately 170 million cubic metres.

Field investigation and data collection
As can be seen in Figure 4(a), a GPS receiver was used to locate while investigating. Four field investigations, located at the point of f 1f 4 in Figure 4(c), have been further discussed in section 4. Among them, a debris flow has been found at f 1 in Figure 4(a). Original vegetation flowed away in the debris accident. The slope surface was covered with mixture of bare soil and breaking stones rolling down from a depletion area. And in Figure 4(b), a differential road sedimentation was found at f 2 . The demarcation could be identified clearly, and the uplifting amount was measured around six centimetres as can be seen from an enlarged view. The field investigation results indicate the landslide movement and studying it is of great value.
Data collecting process is presented in Figure 4(c) and (d). The boundary of Xishancun landslide is labelled with dark dash lines and the study area is circled by a dark solid line. As abovementioned in subsection 3.1, Xishancun landslide has a length of about 4 km and active regions are mainly located in the middle and toe of the whole landslide body as studied before (Qu et al. 2016;Shi et al. 2016). A long-distance terrain laser scanning system (Riegl VZ4000) was used to collect point cloud data. It has a maximum scan distance of 4 km, thus could acquire redundant points for study. Equipped with a full waveform module, terrain points could be detected even at vegetation dense regions.
Two phases of datasets were acquired on 13 April 2014 and 11 April 2015. Basic information of each scan is listed in Table 1. Scanning angles were 57.410°horizontally and 92.404°vertically in 2014. The corresponding parameters were 49.985°and 55.008°i n 2015. Horizontal resolutions were 0.010°and 0.013°and vertical resolutions were 0.013°and 0.006°, respectively. Scan frequency was set 30 Hz for both surveys. Sixhundred and ten million and 530 million points were recorded in two surveys. In

Data processing and results
The flowchart of applying the PCSC technique to obtain landslide displacement is shown in Figure 5. Original LiDAR datasets were pre-processed and two terrain point cloud sets (PCS a and PCS b ) were acquired then. After applying the PCSC technique, the displacement field Z PCSC is calculated with magnitude and direction. Finally, the performance of PCSC technique was estimated by comparing its results with InSAR's (Z LOS,InSAR ) . Differences (ΔZ) between them were evaluated at last.

4.1.
Data pre-processing 4.1.1. Registration of two point cloud sets First, aimed at aligning two point clouds, a manual registration was performed using the first-survey result as the referenced point cloud set. Buildings were clearly identified from dense points at the bottom of the mountain. They were considered stable with respect to the entire slope. Therefore, such features as corners and edges belonging to them were selected as to correspond two point cloud sets and points of these features were corresponding points. Standard Deviation Error (SDE) of the manual registration is 9.70 cm in the end. To minimize the alignment error, Multi-Station Adjustment (MSA) was performed then. This tool modifies the orientation and position of each scan position in several iterations to calculate the best overall fit. MSA uses tie points, tie objects, and poly data objects to detect the closest point in the other point cloud set so as to further align scan positions through iterating. Finally, the calculated SDE degraded to 2.18 cm.

Off-terrain objects filtering
Off-terrain objects filtering were performed to obtain terrain points and detect the surface displacements. In this paper, full waveform LiDAR surveying helped to filter off-terrain objects. Generally, the last back-scattered pulse represents terrain points, which is the most relevant for deformation analysis (Jaboyedoff et al. 2012). Therefore, the points belonging to first and other return are generally assigned to the off-terrain class and will be removed during filtering. After filtering processing, the average distribution of point density of two surveys is shown in Figure 6.
As shown in Figure 6, surface point density generally decreases while the scanning range increases. After off-terrain filtering, the surface point density is not evenly distributed. It reaches the highest value of 10 points per square metre in the lower part of the study area, which has a relatively short range with scanner's positions referred to Figure 4(c), and decreases while the range grows. However, there are four representative areas found with extremely low point density (lower than one) marked as Area 1-4 in Figure 6. In Area 1 and Area 2, there are some buildings, trees, and farmland, so most of the acquired points are targeted as off-terrain points. Low density in Area 3 and Area 4 where the vegetation is not widely distributed may result from landforms which degrade the performance of filtering.

Results of PCSC
4.2.1. Calculating the surface roughness As discussed in section 2, determining parameters M and m is a key step of PCSC. These parameters are related to surface roughness. According to Equation (4), roughness values are calculated first. Their average distribution can be seen from Figure 7. These values are classified to five intervals. Lower values are mainly distributed in the northern part of the study area, and higher ones are in the southern parts.
From Table 2 we can see that roughness is related to landforms. Values between 5 and 10 account for 71.3% of the entire study area and the representative object of this interval is bare slope with relatively flat surfaces. Values between 10 and 15 are related to areas distributing along roads. These surfaces are a little rougher because of some human activities. Areas with values between 0 and 5 are distributed with a flat cliff, which mostly lies in the north-east of the study area. Areas with values beyond 15 have the least occupation (0.6%). According to the representative object's picture, the surficial structure is rough and complicated. There seem to be many stones on these areas. Based on the investigation at f 1 in Figure 4(a), the slope surface is roughly covered by breaking stones and bare soil after a debris flow. Thus, high roughness values may be related to post-debris-flow events. Therefore, the parameter of roughness could be used to describe the terrain surface structures.

Displacement field of Xishancun landslide
Displacement field of the landslide is acquired according to the proposed flowchart in Figure 8. First, the point cloud set of the first survey in 2014 is taken as the reference cloud (referred to as PCS a in Figure 1), and the one in 2015 is taken as the conflicted cloud (referred to as PCS b in Figure 1). Parameters in related to roughness, k 1 and k 2 in Equation (4), are set as 10 empirically. Then, the distance is calculated at every core point (virtual point). To obtain the main deforming values, an up-scaling operation on the PCSC result is conducted: a raster with grids of 10 m × 10 m is created and displacement values within a grid are used to calculate a mean value. It is set as the grid value. In addition, a convolutional median filtering with a window of 10 is performed to mitigate the noise effects and smooth the corresponding result. Finally, a displacement field including deforming magnitude and directions is obtained.
As shown in Figure 8, there are certain deforming principles in the displacement results. First, displacements in the northern and eastern areas are considerably larger than the ones in the southern. Areas with the largest displacement values are located in the north-east of the study area, and the annual displacement is up to 246 mm. Second, a compound pattern can be seen in the deforming direction. In the northern and north-eastern areas, the sliding patterns are principally from north-east to southwest. And in the north-western area, they are from north-west to south-east. Then, the directions converge and turn mainly to the south in the middle and toe of the landslide.

Comparison of the results with InSAR
InSAR techniques could be used to detect the landslide displacement with a millimetre accuracy (Herrera et al. 2009;Greif and Vlcko 2012). Qu (Qu et al. 2016) has used nine descending TerraSAR images to calculate the annual changes through a small baseline subset method. As shown in Figure 9, the InSAR results cover most of the study area. The significant deformed areas are mostly located in the northern and eastern parts. However, the InSAR technique does not provide deformation results in some areas. This kind of result may come from two reasons: the influence of vegetation and the displacement magnitude this technique can reach. For example in Figure 9, dense vegetation results in low coherence of the interferometric pairs in Area 5, and the deformation magnitude is beyond the detecting ability of the InSAR technique in Area 6.
The displacement results through the proposed PCSC system are compared with the ones through InSAR system. Difference comparisons are performed in six parts with typically active deformation from the InSAR results, denote as region (a)-(f) in Figure 9. However, the displacements from two systems are calculated in different directions. The InSAR deformation was calculated from a Line Of Sight (LOS) direction, so the PCSC displacement should be projected to the LOS direction through a simple equation according to Qu (Qu et al. 2016) as follows:  where Z LOS;PCSC is magnitude of PCSC displacement along the LOS direction, Z PCSC is magnitudes of measured PCSC displacements, and θ is set 33°according to the look angle of TerraSAR satellites. After Equation (6) correcting PCSC displacements to the LOS direction, comparison of PCSC with InSAR results can be conducted through Equation (7): where Z LOS,InSAR (already in LOS direction) is displacement form InSAR. Z LOS,PCSC (projected into the LOS direction) is displacement from PCSC. ΔZ is their relative deviation. As shown in Figure 10, the PCSC displacement results differ notably from InSAR results. Figure 10 (a)(-(f) describes the deviation (ΔZ) histogram of correspondingly six selected areas in Figures 9 and 10(g) exhibits overall deviation histogram of six areas together. The average ΔZ has been calculated through Equation (8) (d) and (e) they are beyond 3.00 cm and no values less than zero, which means Z LOS;PCSC is all larger than Z LOS;InSAR . Additionally, only in area (c) is ΔZ less than zero (−1.15 cm) and more than a half of ΔZ values are negative. Totally speaking, even though histograms of ΔZ differ in distinct regions, there is a particularly concentrated distribution overall. In Figure 10(g), more than 85% of abstract ΔZ are less than 3.00 cm, less than 1% of them are larger than 5.00 cm, and ΔZ is 1.18 cm. The differences between PCSC and InSAR results are supposed to come from two aspects. 1) Data involved to calculate displacement comes from different time at different time. 2) Projecting the PCSC displacements to the LOS direction of SAR may bring in errors because the projecting equation (Equation 6) is too simple here. However, the comparison above still implicates that deviations between PCSC and InSAR results are at centimetre level. Totally average deviation is 1.18 cm. Aimed at studying the movement features of an oversized landslide, a confidence of centimetre could be acceptable. Therefore, an analysis combined with deforming modes could be further extended based on the displacement field results in the following section.

Trigger factor of results in landslide mechanism
Trigger factors of the displacement result could be originated from the potential landslide movement mechanism. The Xishancun landslide is shown to be a quite complex landslide exhibiting different deforming modes as reported in former studies (Qu et al. 2016;Shi et al. 2016). In this subsection, two types of deformation patterns in this landslide, known as compound deforming and multi-platform deforming, respectively, are further discussed based on the abundant displacement information through the proposed technique.
An area with compound deforming modes is found in the most active part, located at the north-east of the study area as discussed in subsection 4.2. This region has a displacement range of 130-246 mm, and it could be divided into a small part and a large part by cyan dash lines based on displacement characteristics in Figure 11(a). First, in the small part, displacement values are below 190 mm and deforming directions has a tendency of turning south-west from the southeast. On the contrary in the large part, displacement values are all beyond 290 mm and have consistent deforming directions from north-east to south-west. Additionally, there is a minor scarp and a larger scarp forming the landslide platforms marked in blue lines. Details of the minor scarp, in the neighbour of investigation point f 3 , are presented in Figure 11(b). An obvious trailing edge is exhibited, and the scarp bares out with an estimated relative fall of about 10 m in an enlarged view. And the main scarp, with a length of more than 200 m from west to east, is in the lower part of this region. Displacements below the scarp have the largest magnitude, and their directions are relatively west-south-west in Figure 11(a). The main scarp can be identified in an orthophoto image in Figure 11(c). There is dense vegetation distributed above and below the scarp. As a contrast, the scarp looks white because of little vegetation covering.
As described above, there are different deforming parts in this region (i.e. one small part and one large part with two scarps). Their deforming modes could be further analysed through acquired displacement results combined with the topographic interpretation in the view of the landslide mechanism. Therefore, this region is supposed to be compound deformed.
There are multi-sliding platforms in neighbour of investigation point f 4 . Two typical scarps, marked in blue lines, could be recognised through mapping slope degrees in Figure 12(a). They both have the largest slope angle grade of 44-86°. The main scarp is located lower than the secondary scarp. Deforming just below two scarps is consistent in this region in Figure 12(b). Most of it has magnitudes between 150 and 185 mm (coloured by yellow) and has directions from north-east to south-west. The main scarp has a length of more than 100 m. An orthophoto image and an enlarged field photo of the main scarp are  shown in Figure 12(c). The platform, covered by vegetation, could be clearly identified according to the enlarged picture according to it. Therefore, a multi-platform deforming mode can be recognised in this region based on the displacement results.
In this subsection, the trigger factor of results originated from the landslide movement mechanism, has been discussed through further analysing the sliding modes. On the one hand, the compound deforming mode and the multi-platform sliding mode have been successfully identified through abundant displacement information including magnitude and directions. On the other hand, the analysed sliding modes could be explained with the landslide movement mechanism and corresponding field pictures as well. Therefore, cross analysis on deforming modes combined with the trigger factor in landslide movement mechanism could make the displacement results more convincible.

Conclusion and discussion
The LiDAR-based technique stands out for its high accuracy and effectiveness in surveying giant landslides. This study applies a method called PCSC to acquire the surface changes from two temporal point clouds. In this method, registration and off-terrain objects filtering are adopted in the pre-processing. Then, surface displacement results, including magnitude and directions, are calculated through PCSC technique based on terrain roughness. These results are comparable with the ones of InSAR-derived deformations. The differences estimated from two systems are at centimetre level (1.18 cm) in relatively active areas. In addition, based on the displacement results obtained, two landslide deformation patterns are analysed that may improve the understanding of the landslide movement mechanism.
This study indicates that the PCSC technique could provide an informative displacement field for the study of giant landslides. The results have an acceptable difference at centimetre level compared with that of InSAR and they are consistent with the landslide movement mechanism as well. As result, the proposed technique may provide a useful guide for future research.
However, uncertainties originated from certain operations would still not be neglected in evaluating the results, although the corresponding results have been compared with the ones from the InSAR technique. Uncertainties may come from following three factors: • Registration uncertainties. To detect the change between two surveys, strict registration must be performed on these two-point cloud sets. In our study, efforts have been achieved to ensure it in two operations. First, the manual registration method is adopted by selecting certain discernible stable features (i.e. corners and edges of building) as the corresponding points. Then, the MSA method is performed to decrease the errors through iterating. Finally, the SDE decreases to 2.18 cm from 9.70 cm. However, these errors still have negative influences on the PCSC results. They may be the principal reason that the PCSC results have larger differences in comparison with InSAR results. • Displacement calculation uncertainties. PCSC models the calculated distances interactively: two fitting planes are calculated at every core point (virtual point) from the reference cloud and sub-core point from the compared cloud. Parameters as M and m are used to model the fitting planes. Since the fitting planes are used to describe the landslide surface, their parameters should vary distinctly with surface roughness changes. In our study, a simple linear relationship 10 times larger than the local roughness is used to determine M and m empirically. However, this relationship may be more complex (Sagy, Brodsky, and Axen 2007). And identifying each local geomorphology and providing a detailed description is time-consuming. Therefore, the simplification was used to lower the labour cost but it causes uncertainties in the results as well. • Filtering uncertainties. Some filters have been performed to calculate and map the surface roughness and displacement field of the landslide, for example, upscaling the point cloud sets to 10 m × 10 m grids and smoothing the displacement results with a median filter. These operations help to improve results in dealing with issues like noisy, low point density, etc.. However, certain detailed displacement information may be lost within these operations as well. Therefore, they may bring filtering uncertainties in results.

Disclosure statement
No potential conflict of interest was reported by the authors.

Funding
This work was supported by the National Basic Research Program of China [2013CB733204].