Evaluation of landslide susceptibility mapping techniques using lidar-derived conditioning factors (Oregon case study)

ABSTRACT Landslides are a significant geohazard, which frequently result in significant human, infrastructure, and economic losses. Landslide susceptibility mapping using GIS and remote sensing can help communities prepare for these damaging events. Current mapping efforts utilize a wide variety of techniques and consider multiple factors. Unfortunately, each study is relatively independent of others in the applied technique and factors considered, resulting in inconsistencies. Further, input data quality often varies in terms of source, data collection, and generation, leading to uncertainty. This paper investigates if lidar-derived data-sets (slope, slope roughness, terrain roughness, stream power index, and compound topographic index) can be used for predictive mapping without other landslide conditioning factors. This paper also assesses the differences in landslide susceptibility mapping using several, widely used statistical techniques. Landslide susceptibility maps were produced from the aforementioned lidar-derived data-sets for a small study area in Oregon using six representative statistical techniques. Most notably, results show that only a few factors were necessary to produce satisfactory maps with high predictive capability (area under the curve >0.7). The sole use of lidar digital elevation models and their derivatives can be used for landslide mapping using most statistical techniques without requiring additional detailed data-sets that are often difficult to obtain or of lower quality.


Introduction
Landslides pose a great threat to lives and also lead to socio-economic losses across the globe. Worldwide, significant funds are spent mitigating the effects of landslides (Burns et al. 2013). With an increase in population and urban land use, landslide-triggered losses are expected to increase. In order to minimize economic and human losses, potential landslide areas need to be identified and communicated. Landslide susceptibility mapping (LSM) is considered as a critical component in disaster management and mitigation (Paudel et al. 2003;Pearce 2003;Chen, Liu, et al. 2006;Coppola 2006;Mansourian et al. 2006;Sassa et al. 2006;Burns & Madin 2009;Smith 2013). LSM techniques have evolved greatly over the last few decades and numerous techniques have been proposed and applied by many researchers. Each technique has its own unique advantages and disadvantages.
Landslides are triggered by a variety of factors, which makes them complex to analyse and predict. For example, topography, geomorphology, geology, and hydrological conditions are some of the most important factors to consider for LSM. However, the data-sets representing specific conditioning factors (e.g., slope, geology, aspect, terrain roughness, vegetation index, distance to rivers, distance to roads, etc.) included in the analysis for LSM are often not consistent between each study. Data source and quality can also vary widely.
With the emergence of light detection and ranging (lidar) technology, high-resolution digital elevation models (DEMs) have become increasingly important in landslide studies (Gutierrez et al. 2001;Hsiao et al. 2004;Chen, Chang, et al. 2006;Schulz 2007;Kasai et al. 2009). Jaboyedoff, Choffet, et al. (2012) have documented the importance of lidar DEMs in landslide research and suggest that the trend of using lidar data for landslide studies will increase substantially over the next 5À10 years given wider availability of the data. Programs such as the 3D Elevation Plan (Sugarbaker et al. 2014) are in place to collect seamless lidar data across the entire United States given its tremendous value. The potential information available from lidar when combined with advanced statistical techniques can strengthen the accuracy of LSM by improving data precision, resolution, accuracy, and completeness, especially in forested areas (Gutierrez et al. 2001;Haneberg et al. 2009;Lato et al. 2009;Ventura et al. 2011).
While these previous studies have shown the effective use of lidar data in landslide delineation and inventory Jaboyedoff, Oppikofer, et al. 2012), these studies also indicate that much work is still needed to evaluate the use of lidar in automated, predictive mapping. Given the rapid state of lidar technology evolution, many features of the data have not yet been exploited such as the ability to quantify topographic attributes at a fine scale and link those with hydrological parameters. Further, while many landslide mapping techniques have been applied on individual lidar data-sets with success (Schulz 2007;Burns 2010;Derron & Jaboyedoff 2010;Duplantis 2011), few studies have examined the same study area using multiple techniques to evaluate the consistency of results between techniques. Further, it is important to exploit the complete benefits of lidar data by performing LSM analysis on solely lidar-derived conditioning factors. Independent of our study, Jebur et al. (2014Jebur et al. ( , 2015 recently compared LSM using solely lidar-derived data-sets to those produced considering additional factors from other data sources. In their study, they used three different techniques for LSM. Their inventory contained 31 landslide samples (21 training and 10 testing samples) within a small catchment area (31 km 2 ). However, these landslide samples in the inventory were not mapped or derived from lidar data; rather they were derived from aerial photographs, satellite imagery, and field surveys. Their study concluded that solely using lidar-derived factors 'may possibly be sufficient' to produce LSM with acceptable accuracy and other factors did not significantly improve the maps. Jebur et al. (2015) performed similar work using more advanced methods.
Our study builds on this prior work to test the performance of LSM utilizing solely lidar-derived data-sets with a much expanded scope. First, our study considers a significantly larger area (134 km 2 ) with a detailed landslide inventory (581 landslide deposits and 672 scarps) generated directly from lidar data. The improved inventory, which contains polygon delineations of both scarps and deposits, enables the generation of robust training and test samples, which are independent of each other. Next, our study also considers several additional techniques (including advanced machine learning) beyond those considered in the Jebur et al. (2014) study. Finally, we analyse these techniques and the suitability of lidar-derived condition factors in a location with a much higher landslide susceptibility than the area considered in the Jebur et al. (2014Jebur et al. ( , 2015 study.
To this end, the primary objectives of this paper are to: 1. Investigate the applicability and performance of terrain factors from lidar (slope, slope roughness, terrain roughness, stream power index (SPI), and compound topographic index (CTI)) in LSM in a landslide prone location. Are these factors sufficient alone to produce a satisfactory landslide susceptibly map without requiring additional data-sets representing more complex conditioning factors? While Jebur et al. (2014Jebur et al. ( , 2015 provide preliminary results to indicate their suitability, additional studies are needed to confirm these results in a wide variety of terrain types as well as using a variety of techniques. 2. Evaluate the differences in LSM using six widely used and representative statistical techniques in a regional-level analysis. Can more advanced machine learning techniques provide higher predictive capability? While additional, high-quality data-sets representing other features can surely improve the results of any LSM effort, they were not included in this study so that the effectiveness of lidar data alone could be evaluated. It is hypothesized that the high-resolution topographic information fingerprinted into the lidar data encodes critical geologic and hydrologic properties that are useful for determining susceptibility. For this reason, this study did not directly incorporate the detailed surficial geologic mapping that was available for the study area. Additionally, a proper treatment of geology is not straightforward in many mapping techniques; sometimes it is incorporated as a numerical parameter and in other cases as a category. This, in turn, renders it difficult to compare the techniques themselves when it is included. In an evaluation of various factors for western Oregon, Sharifi-Mood et al. (2013) show that many of the topographic-based conditioning factors varied significantly by lithological unit and recommended that analyses be separated by geologic or lithological units.
Despite the wide array of techniques available, it should be noted that few studies utilize more than one technique for a given study area; hence, it is difficult to compare the effectiveness of each individual technique. More detailed information on each of the specific statistical techniques selected for this study is presented in the methodology section.

Conditioning factors
Based on prior literature, landslide conditioning factors are often grouped into two categories, namely intrinsic and extrinsic factors. Intrinsic factors are slope, aspect, hydrological characteristics, and geologic conditions, while extrinsic factors include precipitation, earthquakes, and human activities (Guzzetti 2000;Lee et al. 2004;Guzzetti et al. 2005;Raia et al. 2011;Santangelo et al. 2012;Ardizzone et al. 2013). Numerous studies have recorded the importance of slope, elevation, aspect, flow rates, flow direction, normalized difference vegetative index (NDVI), SPI, topographic wetness index, distance from streams, and distance from roads (Sidle 1984;Van Westen et al. 1997;Hewitt 1998 Extensive work is involved in the preparation of data layers to represent these conditioning factors before using them for landslide mapping. Collection, pre-processing, post-processing, and field validations are common protocols adopted in generation of the data-sets used in landslide mapping. Field collection of data is also dependent on the accessibility to the area and prevailing weather conditions. Often, for many types of data, field observations are transferred to a digital format through digitization or manual data entry. The manual nature of this stage can lead to errors, which can adversely affect the quality of the resulting maps (Ardizzone et al. 1999;Ardizzone et al. 2002;Booth et al. 2009).
Because of the difficulty in acquiring reliable data for many factors, information that can be extracted from remote sensing data is of high value. While such information may not be able to yield as accurate individual point measurements, remote sensing data provides broader spatial coverage. In particular, the high-resolution capabilities provided by lidar help infer information about these broad ranges of conditioning factors within a single data-set. For example, hydrological properties such as CTI and SPI can be derived from a DEM and provide information about the presence and activity levels of water shaping the topographic surface.

Study area and inventory
Landslides are one of the most common and most devastating geohazards in Oregon (North 1964;Swanson and Dyrness 1975;Palmer 1977;Bottom et al. 1985;Miles & Swanson 1986;Sessions et al. 1987;Frissell & Nawa 1992;Komar & Shih 1993;Robinson et al. 1999;Hoffmeister et al. 2002;Wang et al. 2002;Burns & Madin 2009). In Oregon, landslides contribute over $10 million of economic losses every year, and Oregon is one of the seven pilot states which derive funding from the U.S. Geological Survey for loss estimation. With the increase of population and urban land use, landslide-triggered losses are expected to increase (Wang et al. 2002). The Department of Geology and Mineral Industries (DOGAMI) has mapped approximately 37,226 landslides within the state of Oregon . Recent technological advances such as lidar, enabling more detailed topographic information have aided DOGAMI to map more landslides in the last 5 years than those mapped over the last 60 years (Burns, Duplantis, et al. 2012). Figure 1 shows the study area, the Gales Creek Quadrangle, which is located in north-west Oregon (123 07 0 30 00 W to 123 15 0 00 00 W and 45 30 0 00 00 N to 45 37 0 30 00 N). The 7.5 0 USGS quad covers approximately 144 km 2 ; however, it was cropped to a rectangular area of 134 km 2 since a small portion of the quad did not contain lidar data. Approximately one-third of the study area contains high topographic relief. The mean annual precipitation (over the last 30 years) ranged between 110 and 170 cm/year (PRISM Climate Group 2011) from east to west across the study area. Based on the Oregon Geologic Data Compilation (OGDC V5, Ma 2009), the generalized surficial geology of the study area consists of 41% general sediment deposits (predominately Missoula Flood deposits), 37% volcanics (predominately Tillamook volcanics), 14% marine sediments, and 8% intrusive rocks (predominately Columbia River Basalt). Landslides occur frequently in most of the generalized geologic units with the exception of the Missoula flood sediment deposits, which exhibits minimal topographic relief and shallow slopes.
An overview of the lidar elevation data used in this study can be found in Watershed Sciences (2009). The DEM for this quadrangle was part of a much larger data collection effort spanning 6420 km 2 . Because statistics were not computed by quadrangle, the results from the overall collection are described herein. The native spatial resolution of the bare earth lidar DEM for the collection was 0.91 m. The point cloud data-set was estimated to have an absolute (open-terrain) vertical accuracy of 4 cm (RMS) based on a closest-point comparison with nearly 17,000 RTK GPS coordinates. The total point density was typically around 7 points per m 2 , with approximately 0.7 ground points per m 2 . Data were acquired from a flight height of approximately 900 m above ground with a scan angle of §14 from nadir and 50% flight-line side-lap (100% overlap).
The original lidar DEM was resampled to 10 m using a bilinear interpolation technique. The data were resampled so that derivative products from the DEM such as slope are more representative of overall trends and topographic features of interest rather than smaller, localized features. For example, at high resolution small drainage ditches can be falsely identified as landslide prone because of the steep slope, but are generally small enough to not be of concern for landslide studies. Based on the findings of Mahalingam and Olsen (2015), a 10-m spatial resolution was selected for further analysis in this study. While higher resolution provides additional detail, it often results in many false positives for LSM across a larger scale. However, for inventory mapping, the higher resolution data are critical.
DOGAMI has performed detailed landslide mapping for this study area using lidar data . The landslide inventory used in this study was obtained from the Statewide Landslide Information Database of Oregon (SLIDO) v3.0 (Burns et al. 2008. Debris flow deposits were removed from the data-set because of their differing failure mechanism. Small deposits measuring less than 0.0007 km 2 were filtered from the database since 70% of these small deposits were subsequent failures mapped on the top of larger landslide deposits. After the processing and filtering, the SLIDO database for this area contains 499 complex movement type deposits, 45 earth flow deposits, 3 rotational earth slide deposit, and 34 rotational rock slide deposits. Of these deposits, 532 were classified as deep-seated landslides and 47 were classified as shallow landslides. Table 1 documents the details of the landslide inventory used in the analysis. Landslide scarps and deposits cover approximately 20% of the study area.

Methodology
The key stages and steps in the workflow to evaluate the landslide hazard mapping methodology consists of four primary stages: (1) prepare, (2) analyse, (3) produce, and (4) validate. ArcGIS 10.1 software was used for data processing throughout all stages.

Landslide conditioning factors
For this study, we evaluate the ability to derive representative intrinsic and extrinsic factors from the resampled (10 m) lidar bare earth DEM in order to simplify the data needed for landslide assessment. Table 2 describes each of these factors, how they are calculated, their influence on landslide susceptibility, their value ranges in the study area, and references that have utilized this parameter. Figure 2 shows the DEM-derived data layers representing these landslide conditioning factors for the study area. Of all of these factors, slope plays the most significant role in landslide assessment.

Training and test samples
In order to produce landslide susceptibility maps using bivariate, multivariate, and soft computing techniques, landslide training and testing samples are needed (Pradhan 2010b;Yalcin et al. 2011;Zare et al. 2012). Training samples are used to develop the model and testing samples are used for validation. These samples are chosen based on the completeness of the inventory, size of the study area, and the type of approach used. While there are currently no standard protocols for the selection of training and testing samples, it is important that the training and testing samples be independent of each other Yilmaz et al. 2012;Yilmaz et al. 2013). Basheer and Hajmeer (2000) have documented different splitting ratios used for various statistical methods. According to their summary, Swingler (1996) and Gopal et al. (1999) used an 80:20 ratio for training and testing of their models. More recently, Guzzetti et al. (2012) used a 70:30 ratio for delineating landslides from high-resolution images and Bui et al. (2012) used the same ratio for predicting landslides using soft computing techniques. Random sampling techniques are commonly employed in all of the above-mentioned studies. In order to make the models computationally robust, a 70:30 ratio was chosen.
The landslide inventory used in this study was created rigorously from high-quality lidar data and is a comprehensive inventory of where landslides have occurred. Training and testing samples were gathered from landslide deposit polygons, and scarp flank polygons, respectively (figure 3(a) and 3(b), respectively). The use of points sampled from the landslide deposits for training is expected to result in a conservative landslide susceptibility prediction model. In general, the landslide deposit polygons represent material which has reached equilibrium, forming a post-landslide terrain in which the slope angle is relaxed in comparison to the prelandslide terrain. Unlike the landslide deposits, the scarp flanks represent over-steepened portions of the post-landslide terrain. Given the anticipated conservative nature of the resulting landslide susceptibility prediction models, we expect the scarp flank samples to serve as a suitable candidate for independent validation.
Initially, one sample point per pixel occupied by each landslide deposit and scarp was created. However, the large data-set proved too cumbersome for processing using many of the techniques. Hence, randomly sampled points were taken from the deposits and scarp flanks with an approach to balance the representation of all individual landslide features as well as reflect the internal variability within larger landslide deposits. Therefore, landslide deposit polygons smaller than the mean deposit area (0.044 km 2 ) were represented with two sample points and polygons larger than the mean area were represented with 18 points per polygon. In total, 1720 training points were Table 2. Parameters derived from lidar bare earth model used in this study and their significance in landslide mapping (modified from Gorsevski et al. (2000)).
Conditioning factor As can be readily calculated from ArcGIS hydrological tools. generated from landslide deposits for the training data-set. The testing samples were similarly generated from the scarp flanks. However, because the range of areas is relatively small compared to the landslide deposits, each scarp flank was represented with three point features, resulting in a total of 837 testing samples for validation. Once the testing and training protocol was developed for samples representing the presence of landslides, an equal number of testing and training points were randomly generated outside of the landslide scarps and deposits to represent locations with the absence of landslides. These point features were then combined into a sample database. Values for all conditioning factors including slope, elevation, slope roughness, terrain roughness, CTI, and SPI were then normalized between 0 and 1, so the coefficients could be used to evaluate the relative importance of each factor in the derived models. For each sample, values for each conditioning factor were extracted from the raster layers using a GIS extract multi-valued routine.

Analyse
The following techniques were applied individually to generate a series of susceptibility maps.

Frequency ratio (FR)
FR provides the probabilities of presence and absence of an event for individual conditioning factors by generating weights based on the ratio of area which experienced landslides in the past to the total study area (Lee & Pradhan 2007;Pradhan 2010a;Umar et al. 2014). This technique helps in understanding the spatial relationship between landslides and individual conditioning factors because the values in the classes of each conditioning factor indicate the strength of correlation with the event. In order to compute the FR weights, the area ratio for landslide occurrence points was identified for the class of each conditioning factor considered in the current study. The landslide inventory deposit polygons were rasterized and overlaid with the conditioning factors. The area ratio for the class of each factor to the total area was obtained. The FR weights can be obtained by dividing the landslide occurrence ratio in a particular class to the area ratio in that class. Table 3 presents an example of calculated FR weights for slope. A final susceptibility map can be produced by summing up each factor's weights using equation (1). FR weights greater than one indicate higher correlation of that class in triggering landslides: (1)

Weights of evidence method (WofE)
WofE employs a log-linear form of Bayesian probability to estimate the spatial association of landslide conditioning factors and landslide locations. This is a popular method among bivariate analyses because of its flexibility in using thematic factors; in addition, it is often used in expert-based mapping techniques (Neuh€ auser & Terhorst 2007;Vijith et al. 2013;Tehrany et al. 2014). WofE functions on the principle that landslide occurrence can be calculated by understanding the event's prior probability and the occurrences of the conditioning factors. In this current study, WofE was performed by overlaying a binary map with pixels indicating the presence or absence of landslides with each of the conditioning factor classes. Positive and negative weights evaluating the ratios of the presence and absence of landslides for classes for each factor are computed as follows (Akgun 2012): W C D ln PðB=DÞ PðB=DÞ ; (2) W ¡ D ln PðB=DÞ PðB=DÞ ; where P denotes the probability of occurrence of an event (e.g., landslide), B and Bare the presence and absence of the conditioning factors, and D and D are the presence and absence of a landslide, respectively. The difference between these two weights is called contrast (C), which represents the overall weight, showing the association of landslide presence within each class for each conditioning factor. The landslide susceptibility (LS) value is then calculated for each pixel (j) as follows: where C represents the contrast value for the class of each conditioning factor (i) for pixel j, and n defines the number of conditioning factors used. An example of the weights calculated for the slope conditioning factor is shown in table 4. The sign of the C value for each class indicates whether there is positive or negative correlation between the conditioning factor and landslide presence. The weights from all the conditioning factors were normalized and summed in the raster calculator to generate the final LSM.

Logistic regression (LR)
LR is a multivariate regression between landslides (dependent variable) and the conditioning factors (independent variables) within a logit function. In LR, a dependent variable is dichotomous and independent variables can be categorical or continuous (Lee 2004). The regression is based on maximizing likelihood. In this current study, the dependent variable is dichotomous, indicating either the presence or absence of a landslide. However, the independent variables are continuous in nature. LR can be expressed with the logit function as follows: where P is the probability of landslide occurrence which varies between 0 and 1. z is defined with the linear equation (6), which is obtained through an LR, using a maximum likelihood estimation: where b 0 represents the coefficient of the model, andb 1 , b 2 , …, b n represent the regression coefficients. X 1 , X 2 , …, X n represent the independent variables used in the study (e.g., slope, elevation, slope roughness, terrain roughness, CTI, and SPI). If the values are normalized, the regression parameter indicates the relative contribution of each independent variable towards the event happening (e.g., landslide occurrence). The LR was performed using Stat Tools Ò . The resulting coefficients were then multiplied by conditioning factor values for each pixel using the GIS raster calculator to produce the final susceptibility map as follows: z D ¡ 1:340 C 0:1373£S C 0:0006£E C 0:0082£SR C 0:0021 £SPI C 0:6155£CTI C 0:4340£TR; 1where S is slope in degrees, E is elevation in metres, SR is slope roughness in degrees, SPI is the stream power Index, CTI is the compound topographic index, and TR is terrain roughness in metres.

Discriminant analysis (DA)
DA is a statistical method which explores the combination of independent variables that contribute to the presence and absence of landslide events (Reger 1979;Carrara 1983;Gorsevski et al. 2000;Dong et al. 2009). Stepwise DA enables the addition or removal of the independent variables based on the outcome of analysis. In the application for landslide mapping, two groups are established such as stable and unstable terrains. It is assumed that two groups are distinct in nature and mapping units can be categorized into one of these groups. In the context of LSM, DA can help determine the group membership of the mapping unit. This is achieved by identifying the linear combination of landslide conditioning factors which maximizes the difference between stable and unstable zones. Each conditioning factor plays a role in segmenting the mapping unit into its respective zones. The relative contribution of these factors in the identification process can be studied from equation (8). Application of this function for the unidentified units allows for identification of areas susceptible to landslides (Van Den Eeckhaut et al. 2009). DA optimizes a linear combination of the independent variables (e.g., conditioning factors) through multilinear regression based on a leastsquares error minimization approach. DA follows a similar process to LR; however, in LR, a logit function is used and a maximum-likelihood approach is utilized for the regression. Dichotomous points were selected representing the presence and absence of landslide events and normalized values for each conditioning factor at those locations were extracted. The following shows the linear form of the resulting model (Baeza & Corominas 2001): where z represents the DA score, C 1 ÀC n represent coefficients\weights for each conditioning factor, and a is an intercept coefficient. The DA regression was performed using Stat Tools Ò . Several models were generated for various combinations of each conditioning factor. Using testing samples for validation, a success rate curve is plotted to identify each model's ability to correctly identify the given input samples. The model with the highest success rate was chosen and the coefficients from equation (9) were used to produce the final susceptibility map: z D ¡ 0:0701 C 1:504£S C 0:859£TR C 0:2714£CTI C 0:055 £SR C 0:014£SPI C 0:003£E; where S is slope, TR is terrain roughness, CTI is compound topographic index, SR is slope roughness, SPI is stream power Index, and E is elevation.

Artificial Neural Network (ANN)
ANN is a 'computational mechanism able to acquire, represent and compute a mapping from one multivariate space of information to another, given a set of data representing that mapping' (Garrett, 1994). A back propagation algorithm is widely used and has proven successful in LSM. The model generates weight by generalizing and predicting the outcomes from a set of inputs that were not previously seen (Choi et al. 2010). ANN assigns weights to the layers by understanding the errors between the actual and target output values. A multilayered network used in landslide mapping consists of one input layer, hidden layers, and an output layer. Back propagation algorithms train a network until minimum error is reached between the desired and the targeted output values. When the network is trained, a feed-forward algorithm is used to classify the entire data-set. The prior stage is training and the latter is a classification stage. ANN has numerous interconnected nodes; each node represents a processing element which responds to the weights it receives from other input nodes. As such, each node receives a weighted value (Basheer & Hajmeer 2000;Kiymik et al. 2004;Ermini et al. 2005). In this study, the ANN technique has three layers, namely the input layer which represents the presence and absence of landslides, hidden layers which represent conditioning factor values, and an output layer. The back propagation algorithm employed in this technique enabled training by determining the weights of each factor. The trained model was then used for classification of the test data, which was omitted during training of the model. Weights are primarily determined by modifying the number of hidden nodes and the learning rate between the input and hidden layers, and between the hidden and output layers. RapidMiner Studio Ò was used to run neural network modelling by importing the data from GIS for the landslide presence and absence samples. All of the previously discussed conditioning factors were used independently in the model for analysis. When the RMSE goal (0.01) was not met, the model was set to terminate at 6000 epochs. Initially, a learning rate of 0.01 and random weights between 0.1 and 0.9 were used. An iterative process was implemented to refine those values and ensure an appropriate level of goodness of fit of the model. The final weights were then computed from the trained network with the highest accuracy and normalized to generate the final susceptibility map.

Support vector machine (SVM)
SVM is a popular machine learning technique which employs a supervised binary classifier based on statistical learning theory. Marjanovi c et al. (2011) and San (2014) documented that SVM implicitly maps the original input areas into a high-dimensional feature space, using the training data-set. The optimal hyperplane in the feature space is determined by maximizing the boundaries (Abebe et al. 2010;Pradhan 2013). Support vectors are defined as the training points that are closest to the optimal hyperplane. Classification of the new data is achieved after obtaining the decision surface, which is based on the hyperplane.
In a landslide application, consider a set of training samples with instance label pairs x i ; y i with x i 2 R n ; y i 2 f1; ¡ 1g and i D 1; :::; m, where x is a vector in the input section that encompasses other independent variables. The presence and absence of landslide points are denoted by f1; ¡ 1g in the input space. The two core concepts of SVM are to distinguish the landslide and non-landslide points by developing an optimum separating hyper plane; and the use of kernel functions (table 5) to convert originally nonlinear data to a linearly separable data-set (Yao et al. 2013). Yao et al. (2008) documented that two-class SVM methods produce higher accuracy landslide susceptibility maps when compared to results using more than two classes. The two-class SVM method used in this study utilizes a hyperplane to segment linearly separable data. Mathematically, the orientation of the hyperplane in the feature space is represented by w, a coefficient vector. The offset of the hyperplane from the origin is denoted by b. Lagrangian multipliers are used to determine an optimal hyperplane using equation (10) and is subject to the constraints mentioned in Table 5. Variable functions used in the support vector machine algorithm (modified from Pradhan, 2000).

Kernel Equation Variables Linear
Kðx i ; x j Þ D x i T :x j g; r; and d denote parameters of kernel functions. Polynomial Kðx i ; x j Þ D ðg:x i T :x j C rÞ d ; g > 0 Radial basis function (RBF) Kðx i ; x j Þ D e ¡ gðxi:xjÞ 2 ; g > 0 Sigmoid Kðx i ; x j Þ D tanhðg:x i T :x j C rÞ equation (11): Positive slack of the variables are used to form non-separable cases and is denoted by z 1 (Cortes and Vapnik 1995) in the following equation: Hence, for a non-separable case, equation (10) can be modified as follows: where n(0, 1) accounts for misclassification (Cortes & Vapnik 1995;Sch€ olkopf 2000;LaConte et al. 2005). Xu et al. (2012) reported that selection of the kernel function is very important when implementing the SVM method. Four commonly used kernel functions are linear, polynomial, radial, and sigmoid (table 5). Detailed explanation of these kernel functions can be found in many published research works (Yilmaz 2010;Marjanovi c et al. 2011;Ballabio & Sterlacchini 2012;Xu et al. 2012;Li & Kong 2013;Pradhan 2013;Yao et al. 2013;San 2014;Tehrany et al. 2014).
A complexity constant, C, is used to set the tolerance of misclassification and is defined for all types of kernels. Higher C values allow softer boundaries (over-generalization) and lower values of C create harder boundaries (over-fitting) (Pourghasemi et al. 2013).
Rapid Miner Studio Ò was used to perform the SVM operations and a cross-validation function was used to find the optimal kernel parameters. To begin, ranges were defined for each kernel before running the model. Based on the aforementioned research, each kernel has different impacts on the accuracy of the landslide susceptibility map. Hence, the current study applied the four above kernels to evaluate which kernel was most suitable for the lidar-derived conditioning factors. In addition to those four methods, a novel ensemble method was created by combining normalized FR values from the bivariate analysis with a radial basis function (RBF) kernel. The RBF used in this method utilized the following parameters: g D 0.3, C D 10, convergence epsilon D 0.001, and max iterations D 100,000.
The results of each of the five scenarios were subsequently analysed to determine the optimal kernel resulting in the highest prediction rate. In this study, the combined FR and RBF method was chosen for the development of the landslide susceptibility map.

Produce
Each technique resulted in a raster data-set stored in a geo-database. The landslide susceptibility maps from the six techniques were converted to a common scale of 0À1 and reclassified to very low, low, medium, high, and very high zones based on a natural break classification method to normalize all maps to a similar scale. Figure 4 shows the landslide susceptibility maps produced from six statistical techniques. The test samples from the scarp flanks were overlaid on the output maps to visually understand the distribution of landslide samples within the hazard zones.

Validate
In order to evaluate and compare the analysis and mapping techniques, several validation techniques were implemented including (1) receiver operating characteristic (ROC), (2) area under the curve (AUC), and (3) a histogram analysis. 3.4.1. Receiver operating characteristic (ROC) An ROC describes the true positive and false positive proportions along with the quality of the deterministic and probabilistic methods (Swets, 1988). An ROC curve is shown with sensitivity and specificity values. Sensitivity is the proportion of landslide points correctly classified as susceptible areas. Specificity is the proportion of the absence of landslide points correctly classified as very low susceptibility. An ROC can be constructed using training or test samples. Training samples are used to identify the model's goodness of fit and test samples are used for determining prediction skill (Van Den Eeckhaut et al. (2009). In this current study, ROC was calculated for validation samples to determine the prediction rate. For a model to classify the presence and absence point samples correctly, the proportion of true positives should be higher than the proportion of false positives.

Area under the curve (AUC)
The AUC is a quantitative measure of the predictive performance of the model. AUC ranges between 0 and 1, with higher values indicating better performance. Swets (1988) considers classification models with AUC values greater than 0.90 as highly accurate and those with AUC values less than or equal to 0.5 as poor classifiers. AUC values above 0.7 are generally considered satisfactory in landslide research Pradhan 2013;Althuwaynee et al. 2014a).

Histogram analysis
A histogram analysis was done using the test samples generated from the scarp flanks. The equivalent numbers of absence samples were generated from undisturbed morphological study sections. Values from the output susceptibility maps were extracted on the corresponding spatial location of presence and absence samples. The values were then reclassified using the same bins as used for the output map classification.

Results and discussion
The reclassified landslide susceptibility maps from the six statistical classification techniques are presented in figure 4. All maps, with exception to the one created using ANN, appear very similar. Figure 5 shows the ROC curves and AUC values for each technique. All curves are well above the 1:1 line and AUC values are well above 0.5, indicating that all the six methods could be successfully used with these conditioning factors alone and achieve satisfactory performance. The SVM (AUC D 0.89) and ANN (AUC D 0.72) techniques showed the best and worst performances, respectively. All other techniques shared similar ROC curves and AUC values.
For additional validation, the histogram analysis outputs (figure 6) plot the number of landslide presence (red) and absence (green) samples that fall within each hazard bin (e.g., very high). Ideally, all landslide presence points should plot within the high and very high bins while absence sample points should plot within the low and very low bins. These plots help compare the individual techniques and evaluate their ability to distinguish between landslide prone and stable locales. FR, WofE, LR, DA, and SVM show a clear separation of landslide samples placed in high and very high bins. Except for ANN, all the models classified more than 75% of the test samples as high and very high zones. Further, the ANN output more or less equally distributes landslide presence samples throughout all bins, leading to ambiguous results. A similar trend is noticed for absence samples. Visual analysis of figure 4 shows the ANN map is washed out and does not provide sufficient contrast to distinguish between susceptible and non-susceptible areas. Visual analysis of ANN results also indicate unstable zones derived from scarp test samples classified as stable zones. Therefore, it is noted that the ANN technique needs more ancillary data-sets for improved training and resulting LSM.
In addition to binning and comparing the test sample points, the original polygons can be used for validation. The total area of the scarp polygons in the study area is 4.6 km 2 (table 1). FR, WofE, LR, DA, and SVM identified more than 3.0 km 2 (65%) of the scarps as high hazard zones where as the ANN technique identified only 1.7 km 2 (37%) of scarps as high hazard.
All methods showed that slope was the dominant conditioning factor in landslide susceptibility for the study area. For example, results of using the FR technique indicate a contribution of each slope class in triggering landslides. It shows that landslide occurrences are typically between 21 and 45 . Although less than 4% of the total area has a slope greater than 45 , there is a positive correlation for failure of those steep slopes. C values (table 4) from WofE show similar results with a strong positive correlation between landslides and slope classes ranging from 20 to 51 (table 4). In the LR and DA analyses, the coefficients for slope were significantly larger in the regression analysis using the normalized conditioning factors.
After slope, the individual FR values of each factor indicated that elevation and SPI had the next most significant contribution. Other factors showed correlation with landslide occurrence; however, they were not as strong as those observed with slope. From the FR analysis, elevation values from 350 to 450 m show higher contribution to landslide occurrence. However, values above 450 m begin to show a negative correlation with landslides. SPI demonstrates a positive correlation from 0 to 12, which starts to decrease gradually in value after that range. Slope roughness demonstrates a positive correlation from 0 to 8 but shows a negative trend for greater than 8 . The entire range of CTI and terrain roughness values for the study area results in a negative correlation with landslide occurrence.
Results were similar with other techniques; however, the relative contribution of factors other than slope varied in order of importance. For example, the positively signed coefficients for LR indicated that slope, slope roughness, SPI, and elevation are correlated with landslide occurrence. For DA, ANN, and SVM, weights indicated slope and slope roughness as the strongest indicators when defining landslide prone areas. LR and DA also indicated that all six techniques were found to be statistically significant using both techniques. Each technique has its own advantages in terms of its ability to incorporate input parameters and information provided in interim stages or the final output that can be used for analysis. Bivariate (FR and WofE) analyses help provide a clear understanding of the details of classes in each of the conditioning factors, which were well reflected in the output maps. The multivariate (LR and DA) output maps had the ability to highlight susceptible regions near the streams even in regionally flat regions. The two soft computing techniques, ANN and SVM, had polarized results. On the one hand, the ANN technique produced a smooth output with limited details to delineate the hazard zones. On the other hand, in the SVM map, the details of other factors such as SPI and CTI can be seen in the output map, indicating that it captures these morphological features well in the output maps. When compared to the other models, the slope pattern and elevation dominate the susceptibility maps and show minimal contribution from the other factors. The SVM strength is attributed to the combination of using a RBF with the FR output to generate the output map.

Conclusion
The primary aim of this paper is to analyse powerful statistical techniques using only lidar-derived topographic variables for predictive mapping of areas susceptible to landslides. Most notably, this research shows the advantage of selecting fewer critical conditioning factors for LSM when using a lidar DEM similar to the results of Jebur et al. (2014Jebur et al. ( , 2015. A single-lidar DEM and its derivative products (e.g., slope) proved sufficient to be used effectively and reliably in LSM using bivariate, multivariate, or soft computing techniques. Additional detailed data-sets (e.g., hydrology) that may be difficult to obtain or poor in quality were not required. The ability of a lidar DEM to preserve morphological features enables other conditioning factors to be accounted for in the analysis indirectly, as evidenced by the satisfactory AUC values and ROC curves. Few studies are available in geomorphological literature to evaluate the reliability of various statistical analysis techniques. In this study, with the exception of ANN, all techniques appeared to produce statistically similar and reliable results, despite some minor differences. Based on a north-west Oregon study area, the SVM technique performed relatively well for LSM compared to FR, WofE, LR, DA, and ANN. The overall performance of ANN in LSM was low. The poor performance of the ANN method is attributed to the need for ancillary data to train the data-sets and generate weights.
It should be noted that all statistical techniques are sensitive to quality of the landslide inventory data used for calibration. The current study used a lidar-derived inventory, which may not always be available. Hence, in situations where such an inventory and high-resolution lidar are not available, data layers representing other conditioning factors may be necessary. In particular, geologic mapping can help distinguish areas susceptible to landslides. Geologic information pertaining to material type, age, orientation, and tectonic setting of in situ soil and rock can provide insight into the potential for slope instability.
The information derived from this study can assist researchers, planners, scientists, engineers, and government organizations to use a minimum selection of conditioning factors for LSM with the applied statistical techniques. It is envisioned that these techniques can be evolved in the future for disaster response to rapidly map susceptibility for specific events such as earthquakes, storms, or heavy rainfall provided that the lidar data and landslide inventory are available. While the current study indicates that for the study area, the use of a few conditioning factors from lidar can enable fast and reliable generation of susceptibility maps, it is recommended to use sufficient validation techniques in other locations before relying solely on this data-set for LSM for disaster mitigation purposes. It is also recommended to apply several techniques when mapping an area to compare and contrast the overall results since each technique has its own strengths and weaknesses. Combination of the results from each technique could prove to provide a more realistic assessment and minimize the artefacts of each technique.
Matt O'Banion has a background in geoscience and geotechnical engineering. He is currently working on his PhD in geomatics engineering at Oregon State University. Current Research interests include geohazards, lidar analysis, and immersive visualization.