Received 15 April 2016; accepted 23 May 2016; published 27 May 2016
Landslides are naturally occurring complex geological phenomena that cause significant damage in mountainous regions. Landslide mitigation and risk reduction require mapping of susceptible areas and estimating the likelihood of landslide occurrences  . Landslide susceptibility (LS) deals with the likelihood of landslide occurrence in an area on the basis of local terrain conditions  . Most LS studies follow a simple principle: the past and the present are the keys to the future. The conditions leading to the past and present failures will help in estimating the style, frequency, extent, and consequences of failures in the future  . Over the years, the number of studies concerning methods and the progress in landslide susceptibility have grown rapidly. They involve either qualitative or quantitative modeling  . Many pioneer works in this field concern qualitative studies where the judgment established by experts, based on the data investigated, was used to produce susceptibility maps  . The subjectivity of these methods was addressed by the adoption of quantitative assessment methods such as bivariate or multivariate statistical analysis, logistic regression, likelihood ratio, and weight-of-evidence  -  .
Recently, various machine-learning (ML) techniques have been used in landslide research, more often because of their robustness in handling large complex data. These approaches include Artificial Neutral Networks (ANN), Support Vector Machine (SVM), and Decision Trees  -  . Random Forests (RF) is a relatively new ML technique in the field of landslide research, and is utilized in this study. RF, a bagged trees ensemble, is considered as a superior classification algorithm and is widely used in various fields of data mining such as gene classification  -  . However, the use of RF in landslide research is still limited to a few examples  -  .
Different intrinsic and extrinsic parameters are used to analyze LS, and many of them such as geology, soil depth, soil type, and land use usually have limitations of availability and scale  . Therefore, LS evaluation solely based on a digital elevation model (DEM) has been conducted assuming that topography reflects other factors such as geology and land use  -  . Increased availability of high-resolution global DEMs (e.g., SRTM 1 Arc-Second Global and ASTER GDEM) and recent advances in DEM acquisition techniques encourage this approach.
However, the selection of an appropriate DEM scale is necessary to achieve high precision in LS research  . At coarser scales, terrain presentation may be too smoothed. Therefore, Keijsers et al.  suggest that slope morphology and hydrological patterns are better represented with fine resolution DEMs. However, Tarolli and Tarboton  found that LS prediction performance decreases at finer resolutions because too localized topography does not represent the processes governing landslide initiation. Catani et al.  found that the importance of landside predicting parameters changed with spatial scale, and concluded that for some parameters, scale representing not local values but their trends should be evaluated. However, they did not conduct a concrete study to incorporate the variability of parameter importance at different scales for LS mapping (LSM).
This DEM-based study proposes a novel approach to identify the optimal resolution of each topographic parameter and use those parameters at multiple-optima for an LS study using an RF model. The method was applied to two areas with different geo-environmental settings and landslide density to examine how regional differences affect the scale of each parameter.
2. Study Area
The analysis was carried out at two equal-sized areas in Japan that differ in geological and environmental settings and landslide density (Figure 1). The study area in Niigata Prefecture, Honshu, was selected as a representative of an area with frequent landslides. It has an area of 625 km2 (25 × 25 km), and the location has the highest landslide density in the region. An equal sized area in Ehime Prefecture, Shikoku, was selected as a representative of an area with fewer landslides.
The landslides in Niigata reflect specific geotectonic and climatic settings   . The area lies in a large graben called North Fossa Magna  with active neotectonics  -  . According to a 10 m DEM (see Section 3 for details), elevations in the area range from 5 to 1284 m, with a mean of 369 m. The mean slope is 17.4˚. The mountains in the area mainly consist of Tertiary to Quaternary sedimentary rocks, including so-called “Green Tuff”, and Quaternary volcanic rocks and their deeply weathered materials (Figure 2)  -  . The study area and its surroundings are characterized by heavy snowfall, whose meltwater contributes to rich groundwater and frequent landslides  . The closest meteorological station (Tsunan, Japan Meteorological Agency) receives an average annual precipitation of about 1900 mm (1981-2010) with an average annual snowfall of about 1.35 m (1989-2010). The study area in Ehime is located in central Shikoku. Its elevation ranges from 4 to 1895 m with a mean of 825 m. Mean slope is 31.7˚, much larger than in Niigata. Most of the area is underlain by crystalline schist of the Jurassic complex from the Sanbagawa belt. Low-grade metamorphic greenstones from the Chichibu belt, a Jurassic accretionary complex zone, also dominate the southern section of
Figure 1. Map showing density of landslides in (a) Japan, (b) Niigata, and (c) Ehime, prepared using the land- slide data from the NIED (National Research Institute for Earth Science and Disaster Prevention, Japan).
Figure 2. Geological map of the study areas in (a) Niigata and (b) Ehime.
the study area (Figure 2)   . A meteorological station in the southern part of the study area (Hongawa, Japan Meteorological Agency) and one in the northern part (Niihama, Japan Meteorological Agency) receive average annual precipitation of 3077 and 1305 mm (1981-2010), respectively. The study area is steep and affected by major tectonic lines  , favoring landslides, but their density is lower than that in Niigata.
A 10-m DEM obtained from the GSI (Geospatial Information Authority of Japan) was used for topographic analysis. The minimum elevation for both study areas is about 5 m, whereas the mean and maximum elevations are, respectively, 369 and 1284 m for Niigata and 825 and 1895 m for Ehime. The landslide inventory used was provided by the NIED (National Research Institute for Earth Science and Disaster Prevention, Japan; available online at http://lsweb1.ess.bosai.go.jp/gis-data/index.html). The inventory covers the whole of Japan and was prepared by interpreting 1:40,000 aerial photographs. While the inventory includes many historical landslides identified from topographical discontinuities, small slope disturbances were not included. The use of historical (geomorphological) landslide inventories, which summarize past multiple landslide events  , may enable a robust LSM because it reflects various environmental conditions, and the number of available data tends to be large. However, because of the uncertainties associated with the identification of large failures as unique (single) events  , and lack of differentiation among landslide types (possible inclusion of slow moving earth flows), exceptionally large landslides greater than the 95th percentile of the land-slide area (Figure 3) were not used in this study. Table 1 shows the properties of the investigated landslides.
Topographic parameters were derived from the 10 m DEM. Random Forest was used to construct LS models using the topographic parameters and the landslide inventory. It was implemented in JMP Pro 11.2 (SAS Institute Inc., Cary, NC), and GIS based calculations were performed using ArcGIS 10.3 (ESRI, Redlands, CA) and Python 2.7.
4.1. Topographic Parameters
This study employs 16 DEM-derived topographic parameters (Figure 4, Table 2) previously used in landslide research   -  . Elevation (El) is a measure of the height of a surface above mean sea level. Slope (Sl) indicates the degree of inclination of the surface and shows the rate of elevation change. Slope aspect (Asp) represents the direction of the maximum slope. Drop (Dr), equivalent to the hydrologic slope   , shows the ratio of the maximum change in elevation along the direction of flow between cell-centers. It provides an exact measure of surface inclination in relation to flow. Profile curvature (Pfc) is surface curvature in the direction of slope and affects the acceleration and deceleration of surface flow. Plan curvature (Plc) is surface curvature perpendicular to slope direction and affects the convergence and divergence of surface flow. Total curvature (Cr) reflects both the plan and profile curvatures and represents the overall surface curvature.
Ridges and channels are fundamental features of terrain morphology, and therefore are used in various terrain analyses  . In this study, drainage density (Dd), the total length per unit area, was computed using a circular
Figure 3. Percentile distribution of landslide area.
Table 1. Statistical properties of the studied landslides in Niigata and Ehime.
Figure 4. Interrelationship of the topographic parameters.
moving window of radius equaling the length of 10 cells. This unit area changing with DEM resolution allowed the mapping of drainage texture ranging from purely local to the regional average drainage density  . The shortest distance to a drainage line (Dtd) and shortest distance to a ridge (Dtr) were also computed. In areas of high seismicity, Dtr is a significant LS parameter reflecting the amplified motion observed at mountain tops  . Cells with flow accumulation higher than a threshold value were identified as drainage networks, while ridges were defined as lines of cells with zero flow accumulation.
Relative relief or internal relief (Ir) is the maximum elevation difference in a unit area  . The elevation- relief ratio (Er) describes the area distribution at different elevations and is defined as:
Table 2. Topographic parameters used for susceptibility modeling.
features with the 10 m DEM and coarser-scale form of hillslope with the coarser DEMs)  .
The sediment transport capacity index (STCI) is equivalent to the length-slope factor of the Revised Universal Soil Loss Equation  . Therefore, it accounts for the effects of topography on both sediment transport and erosion  . STCI is calculated as:
where A is the upslope contributing area (m2), β is the local slope gradient (in degrees), and m and n are constants. Because the sensitivity of erosion predictions is not strongly affected by the values of the constants  , the values m = 0.4 and n = 1.3 suggested by Chen and Yu  for Taiwan, an area of landslide activity comparable to our study areas, were used in this study.
The stream power index (SPI) also describes the potential of channel erosion and sediment transport processes  . It is defined as:
where As is the specific catchment area (upslope contributing area per unit contour length).
TCI is related to the spatial variability of soil depth and sediment transportation capacity   and is defined as:
TWI has been used to describe soil moisture distribution and has been found useful for landslide susceptibility studies; higher TWI values are often found in landslide bodies   . It is defined as:
In addition to the 16 topographic parameters, random integers rand was also used to assess the performance of other parameters according to the parameter ranking provided by the RF model  .
Although the DEM-derived parameters represent distinct terrain properties and processes, their interrelationship may lead to multicollinearity. However, for LSM, Harrell  suggests that multicollinearity does not influence the predictions from training and testing datasets if both have the same type of collinearities. This applies to our study because all parameters used with the training and testing datasets are mathematical derivatives of the same 10 m DEM.
4.2. Random Forest
RF is an ensemble learning method of classification using regression trees which combines the idea of bagging with random feature selection  -  . RF utilizes bootstrap and random techniques to select the subsample of data and predictor parameters while growing an ensemble of trees (hence called “forest”). In addition to constructing each tree using a different bootstrap sample of the data, RF changes how the classification or regression trees are constructed. In classical decision trees, each node is split using the best split among all parameters. In RF, by contrast, each node is split using the best among a subset of predictors randomly chosen at that node. This strategy leads to higher performance than many other classifiers such as discriminant analysis, SVM, and ANN; it also makes RF robust against overfitting   . RF has several other advantages: 1) it does not require assumptions on the distribution of explanatory parameters; 2) it allows for the mixed use of categorical and numerical parameters without using dummy parameters; 3) it can account for interactions and nonlinearities between parameters  .
RF produces multiple outputs to aid the interpretation of results, including out-of-bag (OOB) accuracy estimates and parameter importance measures  . OOB errors from RF classifications provide an alternative to cross-validation. For each tree in the forest, a random third of all observations are held out from the training set, and are referred to as OOB. The OOB error is, thus, the proportion of misclassified observations. The other crucial output is the measure of parameter importance, i.e., the statistical weight of each predictor variable. This study employs this measure to analyze the influence of scale on landslide causative parameters. OOB accuracy estimates provide the predictive efficacy of RF models. The change in generalized R-square (R2), a measure of variance in the dependent variable explained by the independent variables, was analyzed to identify the required number of trees.
Classification data used in an RF model for LSM should contain information about both landslides and no- landslide areas. A single-pixel sampling strategy using the centroid of a landslide polygon was employed. No- landslide points, whose number is equal to the number of landslides, were randomly created in no-landslide areas. Values of parameters for the landslide and no-landslide points were extracted. The data (50% landslide and 50% no-landslide) for Niigata and Ehime consist of 21,324 and 5086 entries, respectively. The data were randomly divided into training (50%) and testing (50%) datasets.
4.3. Multi-Resolution LSM
The seven DEM-scales (10 to 300 m) were applied to construct RF models classifying landslide presence and absence to identify the optimal resolution for each parameter. Figure 5(a) outlines the process. First, the values of the 16 parameters are computed for all the scales, and the values are used in an RF model to classify the landslide data. The scale of each parameter with the highest importance in the classification, determined as an average of 10 iterations, is considered optimal. This process is repeated for all the parameters, and finally, a combination of all parameters at their optimal scales is used to create a multi-resolution LS model and an LS map. The finest grid size among the parameters at their optimal scales is selected as the mapping unit. Figure 5(b) outlines the determination of parameter importance in a multi-resolution LS model. In the hypothetical example shown in Figure 5, Sl at 30 m contributes most to the classification, and therefore is the most important parameter for the LS study.
To identify the minimum number of trees required for a stable RF, the number of trees was gradually increased to 5000. Figure 6 shows that the RF model stabilizes with fewer than 100 trees. This study uses 500 trees to accommodate unseen inconsistencies, following the strategy of Díaz-Uriarte and Alvarez de Andrés  .
Figure 5. Outline of multi-resolution technique for: (a) selection of optimal parameter scale and (b) determination of parameter importance.
Figure 6. Generalized R-squared and number of in an RF model.
Table 3 and Table 4 show the contribution of each parameter in the RF model at each scale for Niigata and Ehime, respectively. The results show that the optimal scale with the highest contribution differs among parameters.
For Niigata, the finest resolution (10 m) is optimal for parameters Dtd, Er, Ir, and Sl, while most of the other parameters exhibit optima at higher resolutions, preferably at 30, 60, and 300 m. Two parameters, Asp and Dd,
Table 3. Contribution of each parameter and each scale to modeling landslide distribution in Niigata. The highest contributor is selected as optimal (averaged over 10 iterations, RF T# 500).
Table 4. Contribution of each parameter and each scale to modeling landslide distribution in Ehime. The highest contributor is selected as optimal (averaged over 10 iterations, RF T# 500).
contribute the most at 300 m, the coarsest resolution. The results from Ehime are similar. For Ehime, parameters Asp, Dtd, Dtr, El, Er, Ir, and Sl are optimal at the finest scale, whereas the other parameters exhibit optima at coarser scales. Similar to Niigata, Dd in Ehime was found to be optimum at the coarsest scale. Figure 7(a) compares the optimal parameter scale between the two study areas. The optimal scales for most parameters show similarity in both areas except Asp, for which the coarsest resolution in Niigata and the finest resolution in Ehime are optimal.
The relative importance of each optimal parameter is shown in Figure 7(b). Parameters Dr, Sl, El, Ir, STCI, and TWI for Niigata and Dr, TCI, Plc, Cr, Sl, and TWI for Ehime constitute the top six parameters explaining the distribution of landslides in each area. Parameters with higher relative importance (≥8; Figure 7(b)) in both areas are Cr, Dtd, Dr, Ir, Sl, TCI, and TWI. Some parameters were found to be of greater importance in one area but not in the other. The parameters with large differences in relative importance between the two areas (≥5; Figure 7(b)) are Asp, El, Ir Plc, SPI, STCI, and TCI; Parameters El, Ir, SPI, and STCI have higher importance in Niigata while Asp, Plc, and TCI have higher importance in Ehime. The random variable “rand” was correctly identified as the least important parameter in the analysis.
Figure 8 presents the testing and training accuracies of the RF models (averaged over 10 iterations) constructed with parameters at various scales. In both areas, the training accuracies (85.08% for Niigata and 95.44% for Ehime) and testing accuracies (79.70% and 78.62%) were found to be highest for the models with parameters at the optimal scales. Among the model iterations, the models closest to the mean testing accuracy were used to produce LS maps. Figure 9 shows 25 km2 parts of the LS maps for the combined optimal scale and the 10, 30, 90, and 150 m scales.
The effectiveness of the model with parameters at the optimal scales was evaluated using receiver operating characteristics (ROC). The area under the ROC curve (AUC) characterizes the quality of a prediction model  . AUC values range from 0 to 1, and values closer to 1 represent a better classification model. Figure 10 illustrates AUC values on test data for RF models at different scales. The model with the parameters at optimal
Figure 7. The optimal scales (upper) and relative importance (lower) of parameters for Niigata and Ehime.
Figure 8. Training and testing accuracies of LS models for different parameter scales for Niigata and Ehime (averaged over 10 iterations, RF T# 500) (Nii-train = training samples at Niigata; Nii-test = testing samples at Niigata; Ehi-train = training samples at Ehime; and Ehi-test = training samples at Ehime).
scales shows an AUC value of 0.877 for Niigata and 0.870 for Ehime. The highest AUC values show that multi-resolution LS modeling outperforms the conventional single scale modeling.
6.1. Parameter Scale
Different terrain parameters vary in different ways when the DEM resolution changes  . This study has shown that the optimum scales for LS modeling also differ according to parameter (Table 3 and Table 4), and they tend to be common for the two study areas (Figure 7(a)). Although it is expected that the finest DEM can describe detailed topography and is hence suitable for LSM, several parameters are more significant at relatively coarse scales (≥30 m). This suggests that the smallest-scale variabilities of these parameters do not well re- present the physical processes of landslide triggering, as suggested by some previous studies   . For example, all curvature-related variables (Cr, Pfc, and Plc) and the composite topographic indices (SPI, STCI, TCI, and TWI) show meaningful influences on landslide susceptibility at scales equal to or larger than 30 m (Figure 7(a)). This suggests that a 3 × 3 moving window for parameter computation properly encompasses a meaningful topographic unit, including both the detachment and deposition areas of a landslide, only at relatively coarse resolutions  . The landslide size distribution of the two study areas (Figure 3, Table 1) also suggests that most landslides are larger than the terrain represented at 10 m resolution. Scaling of Dr with the optimal 60 or 150 m resolution also seems to correspond to the landslide size distribution. In contrast, for Ir and Er, the finest resolution (10 m) is the optimal scale for both areas, and it is ascribable to a larger 10 × 10 moving window used for their computation, which can represent a relatively large area even at the finest resolution. However, although computed using a 3 × 3 moving window, Sl is optimal at the finest scale. This seems to reflect the high sensitivity of slope calculation to DEM resolution. It is widely known that lower resolution DEMs result in lower Sl values for the same terrain  . Because Sl is directly related to gravitational force triggering landslides, its accurate computation using fine DEMs is important. A similar explanation can be given for Dtd, which is also optimal at the finest 10 m scale. Fluvial activity such as channel erosion tends to induce landslides along the river course. The accurate location of rivers is better represented if the finest DEM is used  .
For parameter Dtr, a coarser scale (60 m) for Niigata and the finest scale (10 m) for Ehime were found to be optimal. Ridges extracted at coarser scales usually correspond to major ridge lines, while at finer scales they include local topographic highs  . Dtr for Niigata at a coarser scale could therefore include the amplified motion observed along the major outstanding ridges during seismic events  . Indeed, some landslides in Niigata
Figure 9. Landslide distribution map (grey) and part of LS maps (25 km2 and probability ≥ 0.5) from RF models with different parameter scales for Niigata (upper) and Ehime (lower).
Figure 10. AUC values of ROC curves for test data for the models at various scales (left). ROC curves for the LS model with the parameters at the optimal scales (right).
were due to high seismicity, as was the case of the 2004 Chuestu earthquake  .
The optimal scales of Dd and Asp differ significantly from those of the other parameters; the coarsest scale (300 m) is optimum for Dd, while Asp shows the largest deviation between the two areas, 10 and 300 m. Dd in this study is estimated over a unit area dependent on the scale of analysis, thus at finer scales itmight be related only to the local presence/absence of drainage lines. However, at coarser scales, it can reflect the known relationship between general relief characteristics and landslide occurrences    . For Asp, the finest resolution (10 m) is optimal for Ehime, while the coarsest resolution (300 m) is optimal for Niigata. Local meteorological conditions and their relationship with LS may explain this variation. Both the study areas receive large amounts of precipitation; however, Niigata receives a significant portion of its precipitation as snowfall (1.35 m per year). The increased overburden due to accumulated snow and the increased soil moisture from snowmelt are responsible for landslides there  . Asp at a coarser scale indicates the overall direction of a hillslope and suggests that the difference in the deposition thickness of snow on the windward and leeward sides is crucial for LS in Niigata. In contrast, Asp at a finer scale could depict local variations in micro-climate, such as insolation and related groundwater conditions, which is related to rock weathering  . The close relationship between weathering and distribution of landslides has been reported in Shikoku Island including Ehime  , indicating the effect of finer scale Asp on local climate, weathering, and landslides.
6.2. Parameter Importance
Among the values of relative parameter importance (Figure 7(b)), higher values for Cr, Dtd, Dr, Ir, Sl, TCI, and TWI in both study areas suggest that these parameters are instrumental in landslide occurrences. Landslide probability generally increases with terrain slope because of increased shear stress, and slope is considered very important in LS studies   . Therefore, the higher importance of Dr and Sl is reasonable. The higher importance of Dr compared to Sl confirms that Dr is a more direct representation of local maximum slope. Claessens et al.  provided a similar observation on the different effects of these slope parameters. The higher importance of Ir in both study areas also suggests the importance of topographic steepness.
The importance of Dtd is explained by the bidirectional relationship between fluvial processes and slope failures. While landslides contribute to channel initiation, stream incision also contributes to landslides   -  . The higher relative importance of TCI and TWI may reflect the significance of hydrological variations related to rock weathering and soil properties    . There is a general consensus regarding Cr that landslides are more likely to occur on concave slopes because of groundwater concentration  . In contrast, earthquake-induced landslides may be more likely on convex slopes with higher ground acceleration. The importance of Cr in both study areas therefore hints to such mechanisms controlling LS.
Parameters Asp, Plc, and TCI have markedly higher importance in Ehime than in Niigata. A combination of geological and environmental variables may explain this observation. As noted, the importance of Asp in Ehime may be due to local micro-climatic differences that lead to differential weathering. By contrast, the higher importance of Plc in Ehime indicates the positive influence of horizontal flow movement in LS, as suggested by Nefeslioglu et al.  (Table 2). The area in Ehime receives a larger amount of rainfall than in Niigata, hence increased water concentration may contribute more to landslides. The higher importance of TCI in Ehime can be explained similarly because it is a parameter strongly related to terrain curvature.
For Niigata, the relative importance of El, Ir, SPI, and STCI are evidently higher than in Ehime. The effects of El and Ir may be related to the effects of earthquakes because higher or high-relief areas tend to undergo accelerated seismic movement. SPI and STCI provide measures of erosion by water flowing from the entire upstream area (Table 2). They are hydrological parameters but are different from curvatures that reflect much more localized water concentration and divergence. The importance of curvature in Ehime and that of SPI and STCI in Niigata suggest that hilly and gentler terrain in Niigata requires not local but wider-scale water concentration for active landslides.
6.3. Assessment of LS Models
The performance of LS models analyzed at eight different scales were compared based on their accuracy estimates (Figure 8) and AUC values (Figure 10). The scale dependency of input parameters was also observed in the accuracy estimates of the models (Figure 8). Except for the training accuracies of LS models for Ehime and an LS model for Niigata at 300 m resolution, the accuracy estimates decreases with an increasing analytical scale beyond 30 m. The slight increase in the model accuracy for Niigata at 300 m is the effect of the two parameters that are optimal at that scale (Figure 7(a)). The discrepancy observed with the training accuracies for Ehime might be due to the smaller number of training samples. Higher testing accuracies were obtained for models at coarser scales (≥30 m) than at the finest scale (10 m). The proposed multi-resolution LS technique resulted in a boost of testing accuracies according to the obtained AUC values (Figure 10).
Figure 9 shows examples of local LS maps at different scales. The figure indicates that the usability of an LS map depends on the mapping scale as well as the model used. Selection of mapping scale of each parameter for the proposed multi-resolution approach enables us to achieve higher accuracy LSM.
LSM provides the relative likelihood of future landslides, conditional on local geomorphic and topographic characteristics  . Results of our LS study suggest that a single parameter-scale analysis falls short in accommodating the heterogeneity of geomorphological characteristics of the landslides and their surrounding area. This study proposes a multi-resolution LSM technique to incorporate such variabilities. The method requires an identification of optimum scales for all parameters to best represent the conditions of slope failure. The parameters at different optimum scales are then brought together for the final LSM.
The study has also demonstrated the usefulness of a DEM-based LS analysis in areas without other sets of high-quality thematic data. The analysis of scale and importance of the DEM-derived parameters reveal that while some parameters show similar importance and scale dependency for different regions, environmental differences result in variability between regions. The performance of LS models also suggests that the finest scale of analysis is not always the best. The proposed multi-resolution LS analysis permits higher accuracy LSM than any single-scale analysis. Further study of different areas is necessary to confirm the usefulness of the multi- resolution approach.