Using an Unmanned Aerial Vehicle-Based Digital Imaging System to Derive a 3D Point Cloud for Landslide Scarp Recognition

Abstract

Landslides often cause economic losses, property damage, and loss of lives. Monitoring landslides using high spatial and temporal resolution imagery and the ability to quickly identify landslide regions are the basis for emergency disaster management. This study presents a comprehensive system that uses unmanned aerial vehicles (UAVs) and Semi-Global dense Matching (SGM) techniques to identify and extract landslide scarp data. The selected study area is located along a major highway in a mountainous region in Jordan, and contains creeping landslides induced by heavy rainfall. Field observations across the slope body and a deformation analysis along the highway and existing gabions indicate that the slope is active and that scarp features across the slope will continue to open and develop new tension crack features, leading to the downward movement of rocks. The identification of landslide scarps in this study was performed via a dense 3D point cloud of topographic information generated from high-resolution images captured using a low-cost UAV and a target-based camera calibration procedure for a low-cost large-field-of-view camera. An automated approach was used to accurately detect and extract the landslide head scarps based on geomorphological factors: the ratio of normalized Eigenvalues (i.e., λ1/λ2 ≥ λ3) derived using principal component analysis, topographic surface roughness index values, and local-neighborhood slope measurements from the 3D image-based point cloud. Validation of the results was performed using root mean square error analysis and a confusion (error) matrix between manually digitized landslide scarps and the automated approaches. The experimental results using the fully automated 3D point-based analysis algorithms show that these approaches can effectively distinguish landslide scarps. The proposed algorithms can accurately identify and extract landslide scarps with centimeter-scale accuracy. In addition, the combination of UAV-based imagery, 3D scene reconstruction, and landslide scarp recognition/extraction algorithms can provide flexible and effective tool for monitoring landslide scarps and is acceptable for landslide mapping purposes.

1. Introduction

High urbanization rates and unplanned settlements expose the residents of urban areas to natural hazards-related risks. Earthquakes, floods, and landslides are among the major threats to urban areas [1,2]. During the 1990s, nearly nine percent of the world’s natural disasters were landslides, which are difficult to predict [3,4,5]. Landslides are a natural geological phenomenon that is widely recognized as an important process in the transport of sediment. They occur on a wide variety of spatial and temporal scales in many mountainous areas [6,7,8,9,10,11,12], and identifying the smallest detectable area is important because small landslides are very likely to broaden under heavy rainfall conditions. The accurate detection and quick identification of small landslides are crucial for adopting appropriate mitigation measures and efficient decision-making strategies [13,14,15]. In recent years, an increasing number of studies have been conducted worldwide regarding geohazard susceptibility mapping and risk assessment. Traditional landslide mapping methods are typically based on interpretations of aerial photography and field survey inspections, which are used to better understand the characteristics of landslides [11,16]. The photographic, or image interpretation, approach can be adapted and implemented manually, automatically, or semi-automatically [17,18]. Depending on the scale and spatial resolution of the image, the details of the extracted geomorphologic features can be significantly affected. Singhroy [19] asserted that an aerial image scale of 1:25,000 is considered the largest scale at which one can properly interpret slope instability phenomena from aerial photographs. Mantovani et al. [20] suggested that the best image scale is 1:15,000 because the disrupted topography of a landslide scarp can be clearly identified at this scale. To achieve a landslide map with high accuracy and reliability, manual interpretation requires a well-trained geologist to delineate landslides in a stereoscopic environment. This process is time- and labor-intensive [17,18]. Alternatively, ground surveys with GPS and a total station are also time-consuming and have sparse spatial coverage, which results in the omission of fine-scale terrain structures in the resulting digital surface models (DSMs) [21]. Aside from being inaccessible and time-consuming, ground surveys and fieldwork are also subjective and prone to human and instrumental errors. Traditional strategies are unlikely to provide a satisfactory solution, especially for small-area investigations where mass movement rates are slower [16,22]. For these reasons, researchers have adopted remote sensing techniques in their research methodology [8,9,20,23,24,25,26].

A number of studies have explored the feasibility of an automated approach using satellite or airborne imagery and object-based image analysis [27,28,29,30]. The performance of this approach is limited due to the method’s inability to distinguish the Earth’s surface changes. In addition, improper processing strategies impact the adjacent pixels, resulting in low efficiency, and require machine learning theories, such as image classification and segmentation. This process requires large and remarkably different learning samples or a priori knowledge. It is also difficult to detect the causes of changes in information, which can lead to incorrect identification. Li et al. [31] used multi-resolution segmentation and object features from a digital terrain model and high-resolution satellite images. They encountered problematic results due to the presence of landslide and non-landslide pixels present along the landslide borders, which can directly affect classification accuracy.

Others have attempted to use a combination of satellite imagery and digital elevation models (DEMs). Barlow et al. [32] combined Landsat enhanced thematic mapper plus (ETM+) imagery and the use of DEM-derived geomorphometric data to detect and classify fresh translational landslide scarps within an area of the Cascade Mountains in British Columbia, Canada. They attempted to overcome the problem of unreliable detection of most types of landslides by segmenting the images and using the geomorphometric data. This method had a 75% overall accuracy in the detection of landslides that were over 1 ha². The abovementioned methods are less effective than Light Detection and Ranging (LiDAR) at detecting large landslides that have experienced significant historical activity [33]. The use of LiDAR technology for landslide identification has become increasingly popular in digital terrain modelling [12,34,35,36,37,38,39,40]. LiDAR is a valuable tool in geology [41], geomorphology [12], and hazard reduction efforts. Only a few studies have attempted to develop computer-assisted methods for extracting landslides from single or multiple pulse LiDAR data via pixel-based analysis [4,12,42]. Schulz [33,43] used a LiDAR DEM to identify the topography of the Puget Sound region of Seattle, USA. They found LiDAR data to be useful in identifying complex large-scale landslides and for locating potential landslides in the area. Schulz [33] used a LiDAR DEM with a 2 m resolution to produce slope and aspect-shading maps with azimuth angles of 45°, 135°, and 315° to identify landslides. Miner et al. [44] used a LiDAR DEM to interpret old landslides and identify potential landslides in the Victoria Otway Ranges near the Johanna area of Australia and successfully extracted vegetation from real ground elevation and shaded maps with different aspect directions to verify the landslide detection. Airborne LiDAR data and imagery are safe, accurate, and able to achieve a valuable top-view. However, the extremely high cost associated with the use of aircraft and its time-consuming nature makes this strategy an impractical solution, especially for investigations of small areas. Leshchinsky et al. [45] used LiDAR DEMs and head scarps to perform a semi-automated landslide inventory comparing three different study areas. They noted an increase in computational cost associated with post-processing data, which became especially prominent when one or more of the following parameters increased: study area size, input parameters, and resolution of the datasets.

To overcome the previous drawbacks of airborne LiDAR, the development of other active remote sensing techniques, such as terrestrial laser scanning (TLS), are making it possible to collect dense 3D point cloud for sites of interest. TLS has attracted interest for use in landslide studies, including (1) landslide detection and characterization [46,47]; (2) hazard assessment and susceptibility mapping [41]; and (3) modelling and monitoring [40,48,49,50,51,52]. A comprehensive review of laser scanning technology and its applications in landslide investigations can be found in the work of Jaboyedoff et al. [53].

This technology is not without limitations. These limitations include orientation biases that occur when the TLS line of sight is nearly parallel to the orientation of the discontinuity, occlusions that occur when parts of a rock face cannot be sampled because it is obscured by protruding surfaces, and truncation that occurs when the exposure of the discontinuity is less than the available resolution of the TLS point cloud. The TLS results can have occluded areas that leave gaps in the data due to shadows, truncated data, and orientation biases from sensor positioning. TLS does not eliminate the need for visiting the landslide site because these devices have a limited range, and a larger ranges (i.e., the greater distances between the scanner and the object) lead to lower accuracies.

The precision of the measurements, the amount of time, and money required to conduct these measurements are important considerations when creating a landslide inventory [54]. It is not necessary to reach sub-centimeter precision when monitoring shallow landslides because the effort is directed towards evaluating the entire landslide body [55]. The need for high-accuracy measurements cannot be satisfied by traditional mapping methods due to the limitations in financial and technical resources. High-resolution topographic data are necessary for the morphological analysis of landslides [4]. Identifying the smallest detectable area possible is important because a small landslide has a high probability of broadening as a result of heavy rainfall. Recent advances in low-cost digital imaging, navigation systems, and software development have made it possible to accurately reconstruct 3D surfaces without using costly mapping-grade data acquisition systems. Furthermore, advancements in mobile mapping systems, such as unmanned aerial vehicles (UAVs), have made accurate 3D surface reconstruction more feasible whenever and wherever it is required. When small target areas need to be examined, a UAV is a better and more cost-effective choice than other platforms because UAVs are highly portable and are dynamic acquisition platforms. UAVs combine the advantages of vertical aerial photogrammetry and the high resolution of ground-based images. UAVs have the ability to fly at low altitudes and slow speeds and to reach areas that are not accessible from the ground or with manned aircrafts [56]. Moreover, these systems can be stored and deployed at a minimal expense [56,57,58]. These characteristics make UAVs an optimal platform for affordable rapid responses in mapping applications; consequently there is great interest in utilizing them in natural hazard applications.

Recently, 3D object reconstruction using UAV systems have become a popular area of research. This technology is currently used in a wide variety of applications that span across many fields such as, soil erosion [59,60], landslides or rocky surfaces inventory [56,61], landslide detection [62], landslide dynamics and monitoring [57,63,64,65,66,67], and natural disaster monitoring and evaluation [68,69]. A comprehensive review of remote sensing applications based on UAVs equipped with specific sensor-based technologies has been accomplished to varying degrees across a range of applications. A detailed summary, specifically focused on the progress made in environmental science applications, is provided in Pajares [56] and Colomina et al. [58].

The purpose of this research is to present a comprehensive system that uses a low-cost UAV, Semi-Global dense Matching (SGM) techniques, and automated approaches to detect and extract the geomorphological features of landslide scarps. This system reduces the many limitations that decrease the accuracy, completeness, and reliability of the previously described methods, thereby increasing the effectiveness of the landslide mapping.

2. Study Area

2.1. Location of the Study Area

This study was conducted in a 266 m × 185 m area situated along the main highway from Amman to Irbid in the Salhoub/Al-Juaidieh area of north-central Jordan. The study area is in a mountainous region that receives an average annual precipitation of 200 mm and contains creep landslides induced by heavy rainfall, [70]. The map in Figure 1 illustrates the location and extent of the study area, which is 21 km northwest of Amman and located at the following coordinates: 32.11430°N–35.85689°E. The terrain elevation varies from 650 m to 720 m above sea level and has a slope inclination of between 10° and 70°. The study area is dominated by two significant landslide features, which are known as the Al-Juaidieh slides. Abderahman [71] determined that the Al-Juaidieh landslide was a reactivation of an old landslide. Various field observations across the slope body and deformation along the highway and existing gabions indicate that the slope is still active and scarp features across the slope will continue to open and develop new tension crack features; leading to further drift of rocks downward. Substantial first failure activity in the area dates back to August 1992 during the construction of the Irbid-Amman highway. This landslide affected approximately 250 m of the highway, destroyed two houses, caused the evacuation of 10 others, and delayed the opening of the highway for several months [71].

Figure 1. Location of study area and distribution of the landslides.

Figure 1. Location of study area and distribution of the landslides.

2.2. Geologic Setting of the Study Area

The land-cover types include landslide areas, roads, bare soil, and a few manmade buildings. The vegetation cover on the landslide is predominantly grass with occasional isolated bushes. The geology of the study area consists mainly of the Ajlun Group, which is subdivided into the five formations [72] shown in Figure 2. From bottom to top, these formations are the Na’ur, Fuheis, Hummar, Shu’ayb, and Wadi Es-Sir. The Na’ur limestone formation forms the majority of the outcrops in the study area, which starts with a few meters of a very hard limestone overlain by a medium-hard to somewhat medium-marly limestone, soft marl, very hard dolomitic limestone, soft marl, and intermixed gravel and boulders of limestone and marl at the surface [71]. The Hummar formation forms outcrops of a few meters of marly limestone followed by successive layers of medium-thick bedded grey to whitish crystalline limestone, dolomitic limestone, and dolomite that changes near the top to marly limestone [73]. The thickness of this formation is approximately 65 m within the Wadi Shu’ayb, Suweileh, and Hummar areas [74]. The Fuheis formation consists of soft marl and marly limestone interbedded with thin to medium-bedded shally to nodular limestone. As can be seen in Figure 3, the geological cross section taken along (A-A′ (in the study includes a marlstone layer the bottom and a dolomite layer, limestone layer, and caliche layer at the top [75]. Saket [76] concluded that most of Jordan’s landslides take place within the Lower Cretaceous Kurnub Sandstone, and the Upper Cretaceous Fuheis and Na’ur formations. These formations are characterized by the presence of marls and clays that act as weak interlayers, rendering the slope vulnerable to instability and sliding.

Figure 2. Geological map of the study area (modified after Sawariah and Barjous [73]).

Figure 2. Geological map of the study area (modified after Sawariah and Barjous [73]).

Figure 3. Schematic geologic cross section of the landslide ([71,77]).

Figure 3. Schematic geologic cross section of the landslide ([71,77]).

From a structural point of view, the outcropping rocks in the study area show well-developed joints and cracks. Their widths of these cracks range from a few centimeters to multiple meters, indicating that discontinuity sliding could take place. Furthermore, all the landslides within the area are above the highway ridge and the mass movement along these landslides is rotational in a westerly to northwesterly direction.

3. Methodology

Figure 4 illustrates the proposed framework of the implemented methodology in this research work. This research aims to analyze and examine the automated approaches for landslide scarp detection and extraction. This is achieved first through mission planning and data collection using a UAV equipped with an off-the-shelf large-field-of-view (LFOV) camera. The camera calibration and stability analysis performed for this step can be expressed through the analysis approach proposed by Habib et al. [78,79]. After calibration 3D surface reconstruction was completed using the SGM approach to generate dense image-based point cloud [80,81]. The next step in achieving the research goal is to develop automated approaches for detecting and extracting the geomorphological features of landslide scarps using a ratio of Eigenvalues based on principal component analysis (PCA), examining topographic surface roughness index, and measuring the variability in slope in the local neighborhood of dense 3D image-based point cloud. As proposed by Varnes [82] (see Figure 5), the detection of landslides is based on the identification of the following geomorphological features: (1) a crown with tension cracks; (2) a main scarp, which tends to be the easiest feature to recognize because it is semi-circular (long and narrow) with a steep slope and convex planform and has a main direction perpendicular to the flow direction; (3) a minor scarp; (4) related slide blocks; and (5) a toe bulge or accumulation area, which is higher than the surrounding area due to the accumulation of debris.

Root mean square error (RMSE) analysis and a confusion (error) matrix is performed in order to validate and accuracy examination of the assessment in order to quantify errors as well as to assess the quality difference between the two datasets obtained for the manually digitized landslide scarps (the reference data) and the automated extraction scrap segments using the proposed approaches.

Figure 4. Methodology flow chart.

Figure 4. Methodology flow chart.

Figure 5. Geomorphologic characteristics of large-scale landslide (modified after Varnes [82]).

Figure 5. Geomorphologic characteristics of large-scale landslide (modified after Varnes [82]).

3.1. Mission Planning and Data Acquisition

To apply a cost-effectively collect data and detect landslide detection in hazardous, unstable areas a low-cost remote sensing approach can be applied using a UAV and LFOV digital camera. In this study, a GoPro camera was chosen because the UAV’s gimbaled camera mount is specifically compatible with GoPro cameras. The Hero 3+ Black Edition was the highest quality camera available at the time of purchase. This camera had several innate advantages for this study, such as a fixed focal length and a robust design. GoPro cameras are built for use in rugged environments, which makes them an ideal choice for a UAV platform. The accompanying UAV is the DJI Phantom 2 (see Figure 6). The DJI Phantom 2 is designed for hobbyist use and has user-friendly controls. It also offers several other functions, such as an autopilot, no-fly zones, and auto-return home. This UAV offers approximately 25 min of flight duration on a single battery charge, can carry less than 1 kg of payload, and fly up to 1 km from the controller. This design is sufficient to cover a relatively small to medium area. According to Van Den Eeckhaut et al. [83] small landslide area are less than 0.02 km² in size and a medium to large landslide area size is greater than 0.02 km². However, it lacks several features that are necessary to obtain a photogrammetric product, such as time synchronization between the camera and the navigation unit, control of image exposure, and the overlap and sidelap of imagery. Furthermore, the camera has a short focal length that causes considerable barrel lens distortion.

Figure 6. The DJI Phantom 2 unmanned aerial vehicles (UAV) drone equipped with the GoPro camera used in this study.

Figure 6. The DJI Phantom 2 unmanned aerial vehicles (UAV) drone equipped with the GoPro camera used in this study.

The mission planning (flight and data acquisition carried out in September 2014) for the study area involved determining a detailed flight plan, including the flight direction, number of flight lines, flying height, and knowledge of the interior orientation parameters (IOPs) of the mounted digital camera. This configuration depends on the desired overlap, sidelap, and ground sampling distance (GSD) or footprint. The DJI Phantom 2 is equipped with an autopilot package that allows for pre-programmed flight paths, enabling the user to choose horizontal waypoints and the desired height of flight above the ground. The desired overlap (in flying direction) and sidelap (between adjacent flight lines) were set considering the array width to be parallel to the flying direction. The mission planning was designed to achieve the optimum configuration for the area of interest, as shown in Figure 7. Figure 7a illustrates the mission flight plan for image capture that covers the area of interest using north-south flight paths with 80% overlap and 60% sidelap. Figure 7b shows the east-west flight paths with 50% overlap and 30% sidelap which were used to obtain minimum data redundancy and to fill in data gaps caused by shadowed areas, occlusions, and blurry images. This step guarantees success in the automated point matching process and application of the SGM approach. Table 1 summarizes the designed north-south flight configuration for the area affected by the creeping landslide in the study area.

Figure 7. Graphical interface of the image capturing mission flight plan designed to cover the area of interest: (a) north-south flight path and (b) east-west flight path.

Figure 7. Graphical interface of the image capturing mission flight plan designed to cover the area of interest: (a) north-south flight path and (b) east-west flight path.

Table 1. North-south flight path configuration for area of interest.

Flight Configuration
Average flying height (AGL)/speed	50 m/5 m/s
Autopilot	Available
Camera Specs	GoPro Hero3 + Black
Image format (pixels)	3000 × 2250
Pixel size	1.55 µm
Focal length (nominal)	3 mm
Time lapse	~2 s
Image Block Specs
GSD (nominal)	~2.0 cm
Overlap/sidelap %	80/60
Image footprint	83 m × 62 m
Distance between images	14 m
Distance between lines	39.8 m
Number of strips	6
Number of Images	370
Total Area Covered
Study area	48,981 m2 (12.73 min)

Table 1. North-south flight path configuration for area of interest.

Flight Configuration
Average flying height (AGL)/speed	50 m/5 m/s
Autopilot	Available
Camera Specs	GoPro Hero3 + Black
Image format (pixels)	3000 × 2250
Pixel size	1.55 µm
Focal length (nominal)	3 mm
Time lapse	~2 s
Image Block Specs
GSD (nominal)	~2.0 cm
Overlap/sidelap %	80/60
Image footprint	83 m × 62 m
Distance between images	14 m
Distance between lines	39.8 m
Number of strips	6
Number of Images	370
Total Area Covered
Study area	48,981 m2 (12.73 min)

To overcome the difficulties of time synchronization, limited image exposure control, and barrel lens distortion, the following steps are proposed. (a) An accurate analytical camera calibration of a low-cost LFOV camera followed by stability analysis using the ROT (Rotation) procedure was performed to evaluate the similarities and differences between the derived camera calibration parameters [78]; (b) The camera was set to automatic exposure using the camera’s built-in timer. The selected time lapse is a function of the image footprint and flying speed, as shown in Equation (3). The GSD can be estimated based on the assumption that the terrain is relatively flat (see Equation (2)); (c) The exterior orientation parameters (EOPs) of the images can be indirectly recovered as part of the image-based point cloud generation process, which will be discussed later.

$$Pixel Size=\frac{Sensor width \left(mm\right)}{Image width \left(px\right)}$$

(1)

$$GSD=\frac{Pixel Size . Elevation above ground \left(H\right)}{Focal length \left(f\right)}$$

(2)

$$Time lapse=\frac{Sensor width .\frac{H}{f}. \left(1-overlap\right) }{Flying speed}$$

(3)

where time lapse is the time between two successive exposures, sensor width is the number of pixels, f is the camera focal length, H is the flying height, overlap is a user-defined value that represents the overlap percentage between two successive images, and flying speed is a user-defined value.

3.2. Automated Surface Reconstruction

To generate a dense 3D image-based point cloud both Structure from Motion (SfM) and SGM techniques are utilized. The proposed SfM approach is first performed for the automated recovery of image EOPs. Then, SGM [80,81] is adapted for the dense image matching using the derived image orientation from the SfM approach from the scanned scene using the captured images. An overview of the proposed procedure is shown in Figure 8.

Figure 8. Proposed dense 3D reconstruction procedure.

Figure 8. Proposed dense 3D reconstruction procedure.

This procedure was initiated by determining the camera’s EOPs at each imaging station via the SfM approach developed by He and Habib [84]. This procedure is based on a two-step linear solution for the initial recovery of the image EOPs First, point feature correspondences are identified, followed by calculation of the relative orientation parameters relating stereo-images from derived conjugate points. A local reference coordinate frame is established in the second stage. Then, the EOPs of the remaining images are sequentially recovered through an incremental augmentation process. In this approach, tie points are identified using the scale-invariant feature transform (SIFT) feature detection algorithm in individual images as described by Lowe [85]. The acquired images are then normalized based on epipolar geometry while considering the camera IOPs and estimated EOPs, to minimize the search space for corresponding features in overlapping imagery. The algorithm used to determine the epipolar geometry and to define the normalized image plane is described in detail by Cho et al. [86]. This algorithm effectively removes the y-parallax in each stereo-pair of images. The images are projected onto a normalized image plane in which the rows of pixels in one image lie on the same row as in the normalized stereo-pair image.

An SGM procedure is then utilized to find the pixel-wise correspondence between the stereo-paired images. The algorithm used here starts by minimizing a cost function along eight different directions (two horizontally, two vertically, and four diagonally) from a pixel to determine an initial disparity value. The cost function is the sum of all the costs of the disparity along a given path. The path (direction) that has the least cost for a given pixel is used to determine the final disparity value, which is used to relate the pixels in one image to the conjugate pixels in the stereo-pair image.

The next step is a tracking procedure that identifies the matched pixels in all images. In this step, individual pixels are tracked across multiple images using the disparity values from the previous step. Starting with a stereo-pair of images, which are adjacent, and their relative orientation, which is known (determined with the EOPs), a corresponding disparity image can be computed. This image can be used to determine matching pixels in other images. Once a pixel is tracked through the disparity images until a valid disparity value cannot be found, the tracking algorithm for that pixel stops and the coordinates are transformed back into the ground coordinate system for use in the next step.

Finally, a dense 3D point cloud is generated based on the spatial intersection of light rays from all matched pixels. Using the EOPs, IOPs, and pixel coordinates of the corresponding pixels, the collinearity equations can be used to compute the object coordinates of each point. This procedure is conducted as a least squares estimation process, which is updated until the coordinates of the point do not change between iterations. There is a possibility that a pixel may be incorrectly tracked among images (known as a blunder) or that the intersection contains very poor geometry (low precision), and these situations should be avoided. To avoid low precision, the ground coordinates that are not sufficiently precise are eliminated. Blunders can be eliminated by projecting the ground coordinates back onto an image plane and computing the residuals, which should be within one pixel of the original coordinates. In addition, the ground coordinates are not computed if fewer than three pixels contribute to the intersection because this will not create a reliable solution. The resulting point cloud uses an arbitrary ground coordinate system with an arbitrary scale. The next step is to transform the point cloud from the arbitrary coordinate system into a meaningful one using surveyed ground control points.

3.3. Automated Landslide Scarp Features Detection and Extraction

The rapid development of high-resolution DEM data has provided more detailed topographic information. Higher resolution topographic data result in the more accurate detection of landslide scarps, thereby assisting in the monitoring and assessment of landslides as well as the development of mitigation plans. In addition, an in-depth understanding of topographic variations within hazardous landslide terrains is vital for companies developing construction plans for new or existing projects. Identifying landslide-specific spatial features from single surface models is important because not all the changes can be detected through a temporal analysis of landslide-susceptible areas. Previous landslides are the key predictors of the distribution of future landslides [87]. The following methodological subsections focus on detecting and extracting landslide scarps from a single surface using an image-based point cloud that is generated using the SGM approach. High-resolution topographic data was used to extract data on landslide scarps by examining local topographic variability through an analysis of the Eigenvalue ratio, slope, and surface roughness.

3.3.1. Eigenvalue Ratios

Eigenvalue ratios represent the degree of three-dimensional roughness or the crease edge of land surfaces [88,89,90]. Point clouds that result from the SGM approach are comprised of massive amounts of 3D coordinates. The well-known KD-tree data structure [91] was used in this study to handle the point cloud. The KD-tree is a tool that organizes the point cloud and allows for different query processes in the 3D space during the Eigenvalue estimating procedure. The KD-tree is used to search for neighbors within a specified search radius from the query point. A simple method for establishing this local neighborhood is to select the closest points to the query point according to a fixed Euclidian distance [92]. In this study, a 0.5 m radius was selected because of the size of the relevant geomorphic features. Selecting a smaller radius would have resulted in an increase in recognition of non-scarp landscape features, such as shrub vegetation, tree stumps or boulders. The Eigenvalue ratio methodology begins by utilizing PCA to determine the geometric properties of the local neighborhood of image-based points. To check whether a certain point (query point) belongs to a rough surface or a crease edge, the following steps are taken. First, a local neighborhood $(P_n)$is defined to enclose the ($n$) neighbors nearest to the query point. Then, a covariance matrix (𝐶𝑜𝑣) is formed based on the dispersion of the points (P_n) from their centroid$({\overline{P}}_c$), as given by Equation (4). An Eigenvalue analysis is then performed to decompose the covariance matrix into two matrices (Equation (5)). The first matrix (W) is comprised of thee eigenvectors $({\overrightarrow{\mathrm{e}}}_1,{\overrightarrow{\mathrm{e}}}_2,{\overrightarrow{\mathrm{e}}}_3),$and the other matrix (Ʌ) provides their corresponding Eigenvalues (${\mathsf{\lambda }}_1,{\mathsf{\lambda }}_2,{\mathsf{\lambda }}_3$).

$${\mathrm{Cov}}_{3\times 3}=\frac{1}{\mathrm{n}}\mathrm{ }{\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{n}}\mathrm{ }\left(\left[\begin{array}{c}{\mathrm{P}}_{{\mathrm{i}}_{\mathrm{X}}}\\ {\mathrm{P}}_{{\mathrm{i}}_{\mathrm{Y}}}\\ {\mathrm{P}}_{{\mathrm{i}}_{\mathrm{Z}}}\end{array}\right]-\left[\begin{array}{c}{\overline{\mathrm{P}}}_{{\mathrm{c}}_{\mathrm{X}}}\\ {\overline{\mathrm{P}}}_{{\mathrm{c}}_{\mathrm{Y}}}\\ {\overline{\mathrm{P}}}_{{\mathrm{c}}_{\mathrm{Z}}}\end{array}\right]\right){\left(\left[\begin{array}{c}{\mathrm{P}}_{{\mathrm{i}}_{\mathrm{X}}}\\ {\mathrm{P}}_{{\mathrm{i}}_{\mathrm{Y}}}\\ {\mathrm{P}}_{{\mathrm{i}}_{\mathrm{Z}}}\end{array}\right]-\left[\begin{array}{c}{\overline{\mathrm{P}}}_{{\mathrm{c}}_{\mathrm{X}}}\\ {\overline{\mathrm{P}}}_{{\mathrm{c}}_{\mathrm{Y}}}\\ {\overline{\mathrm{P}}}_{{\mathrm{c}}_{\mathrm{Z}}}\end{array}\right]\right)}^{\mathrm{T}}$$

(4)

where

$${\overline{\mathrm{P}}}_{\mathrm{c}}=\frac{1}{\mathrm{n}}{\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{n}}\left[\begin{array}{c}{\mathrm{P}}_{{\mathrm{i}}_{\mathrm{X}}}\\ {\mathrm{P}}_{{\mathrm{i}}_{\mathrm{Y}}}\\ {\mathrm{P}}_{{\mathrm{i}}_{\mathrm{Z}}}\end{array}\right]$$

$${\text{Cov}}_{3\times 3}=\mathrm{W} Ʌ {\mathrm{W}}^{\mathrm{T}}=\left[\begin{array}{ccc}{\overrightarrow{\mathrm{e}}}_1& {\overrightarrow{\mathrm{e}}}_2& {\overrightarrow{\mathrm{e}}}_3\end{array}\right]\left[\begin{array}{ccc}{\mathsf{\lambda }}_1& 0& 0\\ 0& {\mathsf{\lambda }}_2& 0\\ 0& 0& {\mathsf{\lambda }}_3\end{array}\right]\left[\begin{array}{c}{{\overrightarrow{\mathrm{e}}}_1}^{\mathrm{T}}\\ {{\overrightarrow{\mathrm{e}}}_2}^{\mathrm{T}}\\ {{\overrightarrow{\mathrm{e}}}_3}^{\mathrm{T}}\end{array}\right]$$

(5)

The Eigenvectors/Eigenvalues are quite advantageous in determining the geometric nature of the established neighborhood. The Eigenvectors represent the orientation of the neighborhood in 3D space, whereas the Eigenvalues define the extent of the neighborhood along the directions of their corresponding eigenvectors [93]. The relative sizes of the Eigenvalues and the Eigenvectors’ directions can indicate the type of primitive feature. For a rough surface/crease edge point, two of the estimated Eigenvalues will be much smaller than to the third Eigenvalue, for which the conventional equations (Equations (6a) and (6b)) were used in this study. The three normalized Eigenvalues denoted by “λ” were sorted from largest to smallest values as λ₃, λ₂, and λ₁,

$${\lambda }_1{ }_{\approx } {\lambda }_2$$

(6a)

$$\frac{{\lambda }_1}{{\lambda }_2}\ge {\lambda }_3$$

(6b)

3.3.2. Slope

The slope angle is the most important parameter in a slope stability analysis [94] because the slope angle is directly related to landslide probability. Other parameters contributing to landslide probability include geology, soil type, and hydrology [95], but this study specifically utilizes sudden change detection via examination of changes in slope angle. This method can also be useful in the identification and extraction of cliff information for other site-specific needs. The variability in slope is used to detect and extract landslide scarps [96,97,98]. The slope is often initially employed to identify landslides [99], based on the assumption that the slope changes abruptly between two successive scarps and that scarps become more distinct from their surroundings as they evolve. The slope angle for each point in an image-based point cloud can be estimated using Equation (7) because the Eigenvector for each point is calculated based on the PCA analysis in the previous step.

$$\left|\mathsf{\theta }\right|=ta{n}^{-1}\frac{ \sqrt{N{x}^2+N{y}^2}}{Nz}* \frac{180}{\mathsf{\pi }}$$

(7)

where $\mathsf{\theta }$ is the slope angle, and Nx, Ny, and Nz are the normal vector components of the plane.

3.3.3. Surface Roughness Index

Surface roughness can be defined as the irregularity of a topographic surface [96]. The surface roughness index is based on the calculation of deviations between the elevation model’s surfaces fitted in the local range of a moving 0.5 m search radius window. The surface of most landslides is rougher, at the local scale of a few meters, than adjacent stable slopes, which are relatively smoother. Tarolli [100] states that as long as erosion associated with precipitation and time is not severe, old landslides can be identified because there will still be a distinct difference in surface roughness between the landslide area and the local terrain. On the other hand, low average rainfall and semiarid climatic conditions dominating the study area result in poor susceptibility to weathering effects, thus maintaining associated surface roughness features.

This characteristic can be exploited to automatically detect and map landslides captured in high-resolution 3D point clouds. As illustrated in Figure 9, the surface roughness of landslide terrain (bottom image) features higher topographic variability than stable terrain (top image). McKean and Roering [12] and Glenn et al. [4] examined surface roughness and confirmed that landslide surfaces are rougher than the neighboring stable terrain due to the landslide mechanics, surface deformation, and subsidence of material. The algorithm developed in this study utilizes surface roughness information to detect and extract landslide scarps based on the measurement of the variability in local topographic surfaces.

Figure 9. Surface normal representation of topographic surface roughness in a digital elevation model (DEM). The difference between smooth terrain (top) and rough terrain (bottom), illustrates greater variability in the latter [12,101].

Figure 9. Surface normal representation of topographic surface roughness in a digital elevation model (DEM). The difference between smooth terrain (top) and rough terrain (bottom), illustrates greater variability in the latter [12,101].

The surface roughness index is the standard deviation of the object height (h) within 0.5 m local sampling window was used to calculate the surface roughness in this study in Equation (8).

$$\gamma =\sqrt{\frac{\sum \text{}{(h- h\prime )}^2}{n-1}}$$

(8)

where $\mathsf{\gamma }$ is the surface roughness index, $h\prime$ is the mean height of all points within the local window, and n is the number of points within the local search window.

4. Results and Discussion

In this study, three sets of results were produced. The first set of results consists of the camera calibration parameters that were obtained by calibrating the LFOV digital camera. The second set of results consists of the image-based point cloud and orthophotos produced by implementing the image processing using the SGM algorithm for the study area. Finally, the third set of results consists of the landslide scarps detected and extracted using (a) the ratio of the Eigenvalues (λ1/λ2 ≥ λ3) based on PCA analysis; (b) the topographic surface roughness index; and (c) the approach based on measuring the variability in slope in the local neighborhood of the 3D image-based point cloud.

4.1. Dataset Description

The study area (Salhoub/Al-Juaidieh, Jordan), yielded a set of 370 nearly vertical images with overlap of 80% and sidelap of 60% along the north-south flight path direction. These images were captured with a calibrated GoPro Hero 3+ Black Edition camera mounted to a DJI Phantom 2 UAV over a total flight time of 12.73 min. To obtain the lowest data redundancy while filling in data gaps caused by shadowed areas, occlusions, and blurred images, a set of 160 images was collected along east-west flight paths with an overlap of 50% and a sidelap of 30%. The utilized GoPro camera was equipped with several imaging modes; however, in this study, only the medium field-of-view mode was utilized. Each image covered an area of approximately 83 m × 62 m on the ground, resulting in imagery with a resolution of approximately 2 cm (GSD). Figure 10 shows an example of a strip of aerial photographs in the study area.

Figure 10. Strip of aerial photographs of the study area captured from the areas indicated in red along of the east-west flight path shown in the upper image.

Figure 10. Strip of aerial photographs of the study area captured from the areas indicated in red along of the east-west flight path shown in the upper image.

4.2. Camera Calibration and Stability Analysis

The camera was calibrated and checked for stability over a two-month period from January 2014 to March 2014. It should be noted that, over the two-month the camera is mounted on a UAV that is flown in different environments, difficult terrain, and have been subjected to rough landings, significant vibrations from the UAV, and even a few crashes. For each image dataset, a total of 12 convergent images of a 2D calibration test field with 116 targets were captured with the same image network configuration. The ultimate goal was to decide whether the two sets of IOPs were equivalent. If the sets of IOPs are similar, then a stability analysis of the camera is ensured and the IOP parameters have not changed over time. Thus, the camera can be deemed stable, which is essential for achieving an accurate 3D model later on. The experimental results were examined using the United States Geological Survey (USGS) Simultaneous Multi-frame Analytical Calibration (SMAC) model with R₀ values of 0 and 3 mm with K₁, K₂, and K₃ as the distortion parameters. The derived square roots of posteriori variance (σ₀) values are all smaller than one pixel in size (1.55 µm) according to Habib et al. [102], and the image coordinate measurement had an accuracy range of one-half pixel, which indicates that no blunders were present, that the utilized SMAC model was appropriate and that the derived IOPs were acceptable.

The stability of the camera was analyzed using the ROT method [78]. In the case of the GoPro Hero 3+ Black Edition camera used in this study, the pixel size was 1.55 µm. Therefore, if the square root of the variance component was less than 1.55 µm, which is approximately around one pixel, then the IOP sets were considered similar and stable according to the ROT method [102,103,104] and was the correct choice for an accurate estimation of the 3D coordinates.

4.3. Automated Point Cloud Generation

To generate a point cloud using the SGM methodology, a total of 530 images were captured. The image orientation (EOP) is a fundamental prerequisite parameter in any image-based reconstruction. Orienting the images was accomplished by applying the SfM approach developed by He and Habib [83]. The estimated image position and orientation and the reconstructed sparse point cloud from the UAV image dataset are shown in Figure 11.

Once a set of images was oriented, the surface was digitally reconstructed by implementing the SGM algorithm starting from the known exterior orientation and camera calibration parameters. A dense point cloud with more than 13.65 million points was constructed for the study area using all of the captured images with an average point spacing of 1.5 cm and an average point density of approximately 4431.52 points/m². The dense point cloud obtained from SGM is illustrated in Figure 12.

Figure 11. Perspective view of the constructed sparse point cloud. The red (North to South direction) and blue (East to West direction) dots over the sparse points show the camera positions and orientations during image acquisition by the UAV.

Figure 11. Perspective view of the constructed sparse point cloud. The red (North to South direction) and blue (East to West direction) dots over the sparse points show the camera positions and orientations during image acquisition by the UAV.

Figure 12. The dense 3D point cloud generated from the UAV image dataset collected in September 2014. The blue inset square represents the point density within 1 m².

Figure 12. The dense 3D point cloud generated from the UAV image dataset collected in September 2014. The blue inset square represents the point density within 1 m².

4.4. Detection and Extraction of Landslide Scarp Features

The landslide geomorphological analysis in this study was based on a 1.5 cm point spacing in a dense 3D point cloud generated for the study area. The dense image-based point cloud proved to be a very useful tool for rapidly creating profiles along the slope direction. Using these profiles the vertical walls of the landslide scarps were measured and observed to have a range of 3 m to 5 m. The cross section along the profile A-A′ is shown in Figure 13.

Figure 13. Cross-section (left image) along the profile A-A′ (right image is a subset of the geological map) of the landslide in the central part of the study area.

Figure 13. Cross-section (left image) along the profile A-A′ (right image is a subset of the geological map) of the landslide in the central part of the study area.

4.4.1. Topographic Eigenvalue Ratios

Figure 14 illustrates the variation in the topographic parameters using the normalized Eigenvalue ratio of λ1/λ2 computed in a 0.5 m moving sampling window on a dense 3D image-based point cloud. The map in Figure 14 clearly shows that the boundary of the landslide scarps (yellow-red color in the map) were recognized in detail, which means the normalized Eigenvalue ratio values inside and surrounding the landslide scarps were much less clustered than those in the nearby terrain outside the slide.

Figure 14. Eigenvalue ratio (λ1/λ2) computed with a 0.5 m moving sampling window in the dense 3D image-based point cloud. Higher Eigenvalue ratios indicate landslide scarp features because the values of these areas are less clustered than the smooth terrain surface.

Figure 14. Eigenvalue ratio (λ1/λ2) computed with a 0.5 m moving sampling window in the dense 3D image-based point cloud. Higher Eigenvalue ratios indicate landslide scarp features because the values of these areas are less clustered than the smooth terrain surface.

The extracted features were filtered, and only the areas with λ1/λ2 Eigenvalue ratios greater than λ3 were selected as a candidate scarp points. The landslide scarps that were automatically extracted using the proposed threshold values (λ1/λ2 ≥ λ3) are shown in Figure 15. Some open areas associated with the spacing/opening between two successive landslide scarps were observed to have been filled with recently fractured rocks and sediments. Therefore, after the initial landslide scarp extraction, these areas were also classified as scarps because they presented similar characteristics (i.e., top and bottom scarps). Further refinement to filter out noisy extracted points (i.e., points enclosed between higher candidate scar points circled in red in Figure 16) based on the height profile (cross section) along the slide direction was required to obtain more accurate results. In Figure 16, the height profile along the slide direction is in blue, the preliminary scarp points are in green, and the accepted scarp points are circled in red.

Figure 15. Extraction of landslide scarp features for the of Salhoub/Al-Juaidieh landslide based on variation in the local topography’s Eigenvalue ratio.

Figure 15. Extraction of landslide scarp features for the of Salhoub/Al-Juaidieh landslide based on variation in the local topography’s Eigenvalue ratio.

Figure 16. Height profile (cross section) along the slide direction in blue, preliminary scarp points in green, and accepted scarp points circled in red.

Figure 16. Height profile (cross section) along the slide direction in blue, preliminary scarp points in green, and accepted scarp points circled in red.

The extracted features based on the proposed threshold were overlaid on the orthophoto image (Figure 17). Visual inspection indicates that the features related to landslide scarps (blue segments) were accurately extracted and accurately overlaid on the corresponding landslide scarps.

Figure 17. Final landslide scarp detection results (blue) based on Eigenvalues ratios of the topographic surface (left image) overlaid on the generated orthophoto using UAV images of the study area (right image).

Figure 17. Final landslide scarp detection results (blue) based on Eigenvalues ratios of the topographic surface (left image) overlaid on the generated orthophoto using UAV images of the study area (right image).

4.4.2. Topographic Slope Surface

The 3D topographic image-based point cloud data were used to generate maps of the local topographic slopes. The slope angle for each point of the image-based point cloud was derived using a 0.5 m moving sampling window and was determined by calculating the slope of the normal vector for each point. As shown in Figure 18, the resulting slopes of the study area were distributed between 0° and 90°. The slope calculation results show that the active landslide scarps had slopes up to 50° associated with the abrupt changes in the Earth’s surface compared with their local neighborhood and, to a lesser degree, they formed a ridge of the depletion area (toe area) that is comprised of marlstone sediments.

Figure 18. Slope map of the creeping Salhoub/Al-Juaidieh landslide calculated using a 0.5 m moving sampling window on a dense 3D dense image-based point cloud.

Figure 18. Slope map of the creeping Salhoub/Al-Juaidieh landslide calculated using a 0.5 m moving sampling window on a dense 3D dense image-based point cloud.

For landslide scarp detection using the slope method, the point cloud classified the ground objects into either landslide or non-landslide classes based on the assumption that the slope changes abruptly between two successive scarps and that the scarps become more distinguishable from their surroundings as they evolve. In the analysis of the slope distribution in this study, the classification of landslide scarps usually occurred above a mean slope value of 22°, which corresponds to the internal frictional angle at which most slope instabilities in the area occur [28]. The statistics for the slope distribution, derived from the point cloud, are shown in Figure 19.

Figure 19. Statistics of the slope distribution.

Figure 19. Statistics of the slope distribution.

An estimate of the threshold was based on the fact that the histogram of the slope value was set to a mean slope value of 22°. The obtained threshold provided a filtering process to extract the most likely landslide scarps in the study area. The extracted potential landslide scarps from the derived slope variation based on the defined threshold are shown in Figure 20.

Figure 20. Distribution of the identified Salhoub/Al-Juaidieh landslide scarp features based on local slope variations.

Figure 20. Distribution of the identified Salhoub/Al-Juaidieh landslide scarp features based on local slope variations.

The extraction smooth landslide head scarps from the previous results were based on the assumption that the detected candidate points should only decrease in height. In this study, the height profile (cross section) in the slide direction was analyzed to ensure that the smoothed head scarps fit the established criteria (i.e., Figure 16). As seen in Figure 21, the final landslide scarp detection results were overlaid on the orthophoto image. Visual inspection revealed that the detection results were accurately overlaid onto the corresponding landslide scarps.

Figure 21. Final landslide scarp detection results (blue) overlaid on the orthophoto generated used UAV images of the study area.

Figure 21. Final landslide scarp detection results (blue) overlaid on the orthophoto generated used UAV images of the study area.

4.4.3. Topographic Surface Roughness Index

The surface of a landslide is usually rougher on a local scale than that neighboring stable terrain, which means that the local vector orientations of the rougher surfaces more highly variable than those of the smooth topographic surfaces, which have similar orientations. The roughness of a landslide surface is calculated using point clouds without interpolation. To provide a more reliable roughness indicator, the standard deviation of the height point (h) within a 0.5 m local moving sampling window is used to calculate surface roughness of each point in the dense 3D point cloud, using Equation (8) in Section 3.3.3. As illustrated in Figure 22, the surface roughness is higher in the landslide terrain in this study than in the stable terrain. The surface roughness varies from approximately 5 cm to 14 cm, with greater roughness near active landslide scarps and along the higher steeper slopes due to landslide mechanics, surface deformation, and subsidence of material. The smoother topographic surface roughness, which covered the flat surfaces and stable areas, was within a maximum range of 5 cm.

Figure 22. Topographic surface roughness index map of the creeping Salhoub/Al-Juaidieh landslide calculated using a 0.5 m moving sampling window in the dense 3D point cloud.

Figure 22. Topographic surface roughness index map of the creeping Salhoub/Al-Juaidieh landslide calculated using a 0.5 m moving sampling window in the dense 3D point cloud.

The topographic surface roughness map (Figure 22) indicates that the dataset exhibits a good correlation with the landslide scarp boundaries using the topographic surface roughness index approach. To extract the candidate landslide scarps, the histogram of the topographic surface roughness index distribution (Figure 23) was used to estimate the thresholds.

A defined threshold, which was equal to 2σ (0.046 m), was used as a filtering process. The filtering process selected areas with a surface roughness value greater than 0.046 m. Figure 24 shows the extracted features based on the 2σ threshold value.

Figure 23. Statistics of the topographic surface roughness index distribution.

Figure 23. Statistics of the topographic surface roughness index distribution.

Figure 24. Extraction of landslide scarp features of the Salhoub/Al-Juaidieh landslide based on local topographic surface roughness values.

Figure 24. Extraction of landslide scarp features of the Salhoub/Al-Juaidieh landslide based on local topographic surface roughness values.

The final results of the landslide scarp extraction process using the surface roughness index approach are shown in Figure 25a. The extracted landslide scarps are overlaid on the orthophoto image for verification in Figure 25b, and is the result shows that the extraction results are accurately overlaid on the corresponding landslide scarps.

Figure 25. (a) Final landslide scarp detection results (blue) based on the topographic surface roughness index and (b) the scarps overlaid on the orthophoto using UAV images of the study area.

Figure 25. (a) Final landslide scarp detection results (blue) based on the topographic surface roughness index and (b) the scarps overlaid on the orthophoto using UAV images of the study area.

4.5. Accuracy Assessment

The accuracy assessment of the identified results consisted of a comparison between the automatically extracted scarp segments obtained from the different topographic surface analyses and manually digitized scarp segments. The manual measurements were accomplished using GIS software, whereby the active landslide scarps were manually digitized (screen digitizing) by visual recognition of the landslide scarp edges in the orthophoto image of the study area generated using the dense 3D image-based point cloud, and this manual image was treated as ground truth data (reference data) (Figure 26).

Figure 26. Manually digitized ground truth (reference) data for the landslide scarps overlaid on the orthophoto image of the study area.

Figure 26. Manually digitized ground truth (reference) data for the landslide scarps overlaid on the orthophoto image of the study area.

A per-pixel-based approach is used to validate and examine the accuracy of the extracted landslide scarps. The results were compared to the reference data by calculating the RMSE between the automatically extracted landslide scarp segments and the manually digitized landslide reference segments (ground truth). If the average RMSE was above a threshold, then the entire landslide scarp segment was regarded as an unmatched segment. The maximum RMSE tolerance was set to 6 pixels, which is equal to 30 cm (number of pixels multiplied by the pixel size equal to 5 cm) due to the horizontal accuracy in the orthophoto when visually and manually digitizing. This maximum tolerance buffer of 6 pixels was chosen to account for the manually digitized reference accuracy and the minimum spacing between neighboring scarps (~70 cm) observed in this dataset. Using lower tolerance values might miss true matches because of manual reference inaccuracies, and using a higher tolerance value might cause detected scarps to be double matched to neighboring manually referenced scarps.

The RMSE analysis was achieved using the proposed approaches shown in Figure 27a–c. No significant differences in accuracy were found between the different automatic extraction methods in this paper and the manually digitized landslide scarps. The RMSE values, which represent an accuracy assessment, of the Eigenvalue ratio, topographic surface slope, and topographic surface roughness index methods are 11.98 cm, 9.05 cm, and 10.45 cm, respectively. Thus, the positions of the extracted scarp segments may have shifted by approximately two pixels on average from their positions in the reference data. In general, the automatically extracted results were comparable to the results obtained through manual digitization of the orthophoto image of the study area for landslide scarp identification. Therefore, the scarps in this study were accurately extracted.

Figure 27. The average distance (RMSE) between the ground truth (red segments) and the extracted landslide scarp data (blue segments); Scarp segments extracted using (a) Eigenvalues ratios; (b) topographic surface slopes; and (c) topographic surface roughness index values.

Figure 27. The average distance (RMSE) between the ground truth (red segments) and the extracted landslide scarp data (blue segments); Scarp segments extracted using (a) Eigenvalues ratios; (b) topographic surface slopes; and (c) topographic surface roughness index values.

Further validation and examination of the accuracy of the assessment was accomplished by performing a confusion (error) matrix to quantify the errors as well as to assess the quality difference between the two datasets obtained by visual interpretation (the reference data) and the map data (extracted scarp segments) using a per-pixel approach. Each variable of the two datasets was organized in a matrix format populated by the number of agreements and disagreements according to the following four classes: (1) true positive (TP), or matched, which indicates that an extracted landslide scarp correspond to a reference landslide scarp; (2) true negative, which indicates that the extracted non-landslide pixels correspond to reference non-landslide pixels; (3) false positive (FP), or unmatched, which indicates an extracted landslide scarp correspond to reference non-landslide pixels; and (4) false negative, which indicates that an extracted non-landslide pixel corresponds to a reference landslide scarp Matching was performed based on the constant in the pre-defined buffer method. The maximum tolerance of the buffer around the reference landslide scarp data was set at 6 pixels. The parts of the extracted data within the buffer were considered to match the manual reference. The accuracies of the extracted landslide scarp segments based on the different approaches and the reference datasets are shown in Table 2, Table 3 and Table 4. The error matrix was obtained for the per-pixel assessments of quality based on a comparison of the extracted scarp segments as presented in Figure 17a, Figure 21a, and Figure 25a. The reference map is shown in Figure 26.

The overall accuracy of the confusion matrix was calculated by dividing the total number of agreements (i.e., the sum of the diagonal cells of the matrix) by the total number of samples. The user’s accuracies, producer’s accuracies, and Cohen’s kappa coefficient statistics were calculated from the confusion matrix to obtain a more in-depth perspective of the uncertainty analysis. The user’s accuracy (row values) was based on the agreement of a particular class to the summation of all classes in each row. Similarly, the producer’s accuracy (column values) was computed considering the agreement of a particular class to the summation of that column. In many cases, according to Zhan et al. [105] and Tuermer et al. [106], the user’s accuracy represents a measure of correctness (Equation (9)), and the producer’s accuracy represents as a measure of completeness (Equation (10)).

$$\text{Correctness}=\left(\frac{\text{TP}}{\text{TP}+\text{FP}}\right)\times 100$$

(9)

$$\text{Completness}=\left(\frac{\text{TP}}{\text{TP}+\text{FN}}\right)\times 100$$

(10)

Cohen’s kappa coefficient was calculated from the confusion matrix and was estimated for the performance evaluation of landslide extraction. This coefficient is a measure of the agreement between the extracted and reference data. In other words, the kappa statistics are a measure of true agreement, which is represented by the following relationship [107]:

$$\text{Kappa coefficient}\left(\mathrm{k}\right)=\frac{\mathrm{N}{\sum }_{\mathrm{i}=1}^{\mathrm{r}}{\mathrm{X}}_{\text{ii}}-{\sum }_{\mathrm{i}=1}^{\mathrm{r}}({\mathrm{X}}_{\mathrm{i}+}{* \mathrm{X}}_{+\mathrm{i}})}{{\mathrm{N}}^2-{\sum }_{\mathrm{i}=1}^{\mathrm{r}}({\mathrm{X}}_{\mathrm{i}+}{* \mathrm{X}}_{+\mathrm{i}})}$$

(11)

where (k) is the Kappa value, r is the number of rows in the confusion (error) matrix, ${\mathrm{X}}_{\text{ii}}$ is the number of observations in row i and column i on the major diagonal of the matrix, ${\mathrm{X}}_{\mathrm{i}+}$ is the total observations in row i, ${\mathrm{X}}_{+\mathrm{i}}$ is the total observations in the column, and N is the total number of observations included in the matrix.

Table 2. Confusion (error) matrix for assessing the quality of extracted landslide scarps based on the topographic Eigenvalue ratios approach.

		Ground Truth (Reference)
		Positive	Negative	Total	User’s Accuracy/Correctness (%)	Error of Commission (%)
Extracted Features	Positive	3917	846	4763	82.24	17.76
	Negative	1306	26698	28004	95.37	4.66
	Total	5223	27544	32767
	Producer’s Accuracy/Completeness (%)	75.0	96.93	Overall accuracy 93.43%; kappa 74.58%.
	Error of Omission (%)	25.00	3.07

Table 2. Confusion (error) matrix for assessing the quality of extracted landslide scarps based on the topographic Eigenvalue ratios approach.

		Ground Truth (Reference)
		Positive	Negative	Total	User’s Accuracy/Correctness (%)	Error of Commission (%)
Extracted Features	Positive	3917	846	4763	82.24	17.76
	Negative	1306	26698	28004	95.37	4.66
	Total	5223	27544	32767
	Producer’s Accuracy/Completeness (%)	75.0	96.93	Overall accuracy 93.43%; kappa 74.58%.
	Error of Omission (%)	25.00	3.07

Table 3. Confusion (error) matrix for assessing the quality of extracted landslide scarps based on the topographic surface slope approach.

		Ground Truth (Reference)
		Positive	Negative	Total	User’s Accuracy/Correctness (%)	Error of Commission (%)
Extracted Features	Positive	3325	781	4106	80.98	19.02
	Negative	1468	27193	28661	94.88	5.12
	Total	4793	27974	32767
	Producer’s Accuracy/Completeness (%)	69.37	97.20	Overall accuracy 93.14%; kappa 70.78%.
	Error of Omission (%)	30.63	2.79

Table 3. Confusion (error) matrix for assessing the quality of extracted landslide scarps based on the topographic surface slope approach.

		Ground Truth (Reference)
		Positive	Negative	Total	User’s Accuracy/Correctness (%)	Error of Commission (%)
Extracted Features	Positive	3325	781	4106	80.98	19.02
	Negative	1468	27193	28661	94.88	5.12
	Total	4793	27974	32767
	Producer’s Accuracy/Completeness (%)	69.37	97.20	Overall accuracy 93.14%; kappa 70.78%.
	Error of Omission (%)	30.63	2.79

Table 4. Confusion (error) matrix for assessing the quality of extracted landslide scarps based on the topographic surface roughness approach.

		Ground Truth (Reference)
		Positive	Negative	Total	User’s Accuracy/Correctness (%)	Error of Commission (%)
Extracted Features	Positive	3168	450	3618	87.56	12.44
	Negative	1669	27480	29149	94.27	5.73
	Total	4837	27930	32767
	Producer’s Accuracy/Completeness (%)	65.50	98.39	Overall accuracy 93.53%; kappa 71.31%.
	Error of Omission (%)	34.50	1.61

Table 4. Confusion (error) matrix for assessing the quality of extracted landslide scarps based on the topographic surface roughness approach.

		Ground Truth (Reference)
		Positive	Negative	Total	User’s Accuracy/Correctness (%)	Error of Commission (%)
Extracted Features	Positive	3168	450	3618	87.56	12.44
	Negative	1669	27480	29149	94.27	5.73
	Total	4837	27930	32767
	Producer’s Accuracy/Completeness (%)	65.50	98.39	Overall accuracy 93.53%; kappa 71.31%.
	Error of Omission (%)	34.50	1.61

The overall qualities of the landslide scarp segments extracted from the Salhoub/Al-Juaidieh landslide test site using different approaches were calculated to be 93.43% for the Eigenvalue ratio approach, 93.14% for the slope approach, and 93.53% for the surface roughness index approach. The producer’s accuracies (completeness) and user’s accuracies (correctness) were found to be 75.0% and 82.24%, respectively, for the Eigenvalue ratio approach; 69.37% and 80.98%, respectively, for the slope approach; and 65.50% and 87.56%, respectively, for the surface roughness index approach. The main causes of error, in terms of the extracted landslide scarp segments, were due to a few manmade buildings and isolated bushes within the study area. These segments were not extracted because they were directly connected to adjacent bushes and do not show the desired characteristics of landslide scarps. Additional incomplete extracted scarp segments were associated with errors in the visual and manual digitizing due to horizontal inaccuracy present in the orthophoto.

Kappa statistics arrange in from value of 1 or (100%), indicating a strong agreement, to 0 (or 0%), indicating that any agreement is entirely due to chance (i.e., incorrectly extracted). Kappa values greater than 75% indicate very good to strong agreement [108,109,110]. The values of the kappa coefficient associated with the confusion matrices in Table 2, Table 3 and Table 4 were 74.58%, 70.78%, and 71.31%, respectively, thereby denoting good agreement between the extracted results and the reference results.

5. Conclusions and Recommendations for Future Work

This paper presents a practical approach for the detection of landslide scarps using point clouds that have been derived from captured imagery by low-cost unmanned aerial vehicles (UAVs). The main advantages of the proposed methodology is allowing for the derivation of accurate information for landslide characterization while alleviating the inherent risk in surveying hazardous landslide-prone areas and reducing the incurred cost. Due to the inherent excessive lens distortions in the utilized imaging system, GoPro Hero 3+ Black Edition onboard a DJI Phantom 2, a camera calibration and stability analysis procedure is essential. Derived interior orientation parameters (IOPs) for the GoPro during the time period from January to March, 2014 revealed that the camera maintains the stability of its internal characteristics and is suitable for landslide mapping. The paper also introduces procedures for automated recovery of the exterior orientation parameters (EOPs) of the images that have been captured over the investigated landslide location as well as generation of dense point clouds representing the surface within the mapped area. More specifically, a point cloud comprised of more than 13.65 million points—whose average inter-point spacing and local point density values are 1.5 cm and 4431.5 points/m², respectively—have been derived from a set of 530 overlapping images.

The paper presented three automated approaches for the detection and extraction of scarps using the morphometric characteristics of the derived point cloud. The first approach uses ratios among the principal component analysis PCA-based Eigenvalues at local neighborhoods to extract salient geo-morphometric features that represent scarps. The second approach is based on slope variability within local neighborhoods, which are defined by a 0.5 m moving window, while assuming that abrupt slope changes will take place between two successive scarps. This slope change will become more distinct as the scarps evolve. Finally, the last approach is based on evaluating the surface roughness index for the derived point cloud at local neighborhoods represented by a moving window with 0.5 m search radius.

The experimental results show that the proposed approaches accurately identify and extract landslide scarps with a recognition accuracy of approximately 72%. The results also indicate that the scarp detection is 70% complete and 84% correct with an overall quality of 93.4%. These measures have been evaluated by calculating the Root Mean Square Error (RMSE) between the automatically extracted landslide scarp segments and manually-digitized ones. No significant differences in accuracy are observed among the different approaches. The RMSE analysis also revealed that the accuracy of the topographic surface slope, topographic surface roughness index, and Eigenvalue-based approaches are 9.05 cm, 10.45 cm, and 11.98 cm, respectively. Thus, the proposed approaches can accurately identify and extract landslide scarps at the decimeter-level accuracy. Such ability to detect landslide scarps will lead to better understanding of the landslide mechanisms for a given area. This in turn will lead to better identification of the most likely failure site within a landslide prone area and estimation of the volume of potential sliding rock mass. In summary, the developed approaches are fast, economical, labor-saving, and safe tools for detecting and recognizing landslide scarps. They can be also used to monitor and assess the rate of horizontal displacement between the extracted landslide scarps at different times.

The experimental results are based on datasets covering bare earth with little to no vegetation, where the proposed approaches will lead to accurate scarp detection and extraction. Future work will be extended to deal with areas covered by bushes, boulders, and tree stumps. More specifically, filtering techniques will be developed to remove such features before the application of the proposed approaches. In addition, rather than having three independent approaches, future research will integrate the three proposed approaches to explore the possibility of improving the scarp detection and extraction accuracy, especially in areas with higher levels of vegetation cover than the current study site. Finally, the proposed methodologies will be incorporated for monitoring and change detection evaluation while using temporal UAV-based datasets over a given site.