A New Approach Using AHP to Generate Landslide Susceptibility Maps in the Chen-Yu-Lan Watershed, Taiwan

This paper proposes a new approach of using the analytic hierarchy process (AHP), in which the AHP was combined with bivariate analysis and correlation statistics to evaluate the importance of the pairwise comparison. Instead of summarizing expert experience statistics to establish a scale, we then analyze the correlation between the properties of the related factors with the actual landslide data in the study area. In addition, correlation and dependence statistics are also used to analyze correlation coefficients of preparatory factors. The product of this research is a landslide susceptibility map (LSM) generated by five factors (slope, aspect, drainage density, lithology, and land-use) and pre-event landslides (Typhoon Kalmaegi events), and then validated by post-event landslides and new landslides occurring in during the events (Typhoon Kalmaegi and Typhoon Morakot). Validating the results by the binary classification method showed that the model has reasonable accuracy, such as 81.22% accurate interpretation for post-event landslides (Typhoon Kalmaegi), and 70.71% exact predictions for new landslides occurring during Typhoon Kalmaegi.


Introduction
Landslides are disasters which often occur in hilly or mountainous places in the rainy season. They are classified as a category of disasters triggered by hazards related to geological processes (earthquakes, volcanic eruptions, floods in mountainous areas) and meteorology (heavy rains, storms, or typhoons) [1,2]. According to a statistical report on disasters in the last decade (2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015) produced by the International Federation of Red Cross (IFRC), landslides account for 4.7% (177 events) of the total number of reported natural hazards [1]. The highest density of landslide occurrences is in Asia, accounting for 118 of the 177 reported landslides occurring worldwide [1]. A report by the The World Bank, "Natural Disaster Hotpots: A Global Risk Analysis" (2005), noted that Taiwan is a small country with an area of approximately 36,000 km 2 and a population of about 23 million people, but is classified as a high-risk country for such disasters [3]. Taiwan ranked first of the top 15 countries (based on land area) with regard to being the most exposed to multiple hazards (affected by three or more disasters) [3]. In other statistics examining countries which are at relatively high mortality risk from multiple hazards, Taiwan also ranked first of the top 35, with 95.1% of its population living in areas at risk [3]. Besides earthquakes and typhoons, landslides are the most common natural disasters in Taiwan (from 1986 to 2007), and found that four major groups of parameters were used to analyze the landslide susceptibility index, including geological parameters (geology/lithology, fault/lineaments, strata-slope interactions), topographical parameters (drainage, surface roughness, elevation, slope angle, aspect, curvature, slope length) geotechnical parameters (soil thickness, soil texture), and environmental parameters (land-use/land-cover, anthropogenic parameters, rainfall, positions within catchment) [30]. Of these, the five most commonly used parameters were slope angle (97.9%), geology/lithology (92.31%), land-use/land-cover (75.52%), drainage (72.73%), and aspect (59.44%) [30]. This suggests that the above five categories are important factors for landslides. In addition, N. Casagli et al. also suggested that the three key factors causing landslides were topography (slope failure), material composition (geology/lithology), and land-use [2]. Based on the above analysis, five parameters, including lithology, slope angle, aspect, drainage density and land-use, were used to analyze the landslide susceptibility index in this study.

Study Area
The Chenyulan river watershed is located in the central part of Taiwan, in the island's major fault zone (slate belt) ( Figure 1). This river is a tributary of the Zhoushui river. There are important fault systems in this area, such as Chenyulan, Dili, Kuling Chiao, Shenmu, and so on. The average elevation is 1680 m above sea level, and the average slope is 34.5 • . The erosion activities here are complex, and occur at high densities. The area of landslides is 16,598,445 m 2 , or 4.2% of the total of 392,444,312 m 2 in Taiwan. This area is also affected by the tropical monsoon climate of Taiwan, which is warm and humid all year. The rainy season is from May to October, and typhoons often occur in the summer and autumn, from June to October, more frequently in August and September, with on average four direct hits per year [31]. The annual precipitation is more than 2500 mm. The rainfall in Taiwan often concentrates in mountainous areas, and therefore the Chen-Yu-Lan river watershed is also one of the areas with the highest rainfall in Taiwan. Summer rainfall accounts for around 80% of the annual precipitation.
The study area's lithology can be divided into the following: metamorphic rocks are distributed in the east, sedimentary rocks are located in the west, terrace deposits (gravel, sand, clay) are scattered from the central to the south, and alluvium is concentrated along the main river in the north. Geological structures here are very complicated, with many faults and fold systems which are distributed from the north to the south of the area. The fault systems are in two main directions, including the north-south (Chenyulanchi fault, Shanshilchia fault, Dili fault, and Kulingchiao fault) and the east-west (Shenmu fault, Shihpachekeng fault, and Erhyu fault). The fold systems are arranged in the directions of the north-south (Pinglin syncline) and the north northeast-south southwest (Hoshe anticline and Tungfushan syncline). The Chen-Yu-Lan River is located on a large fault (Chenyulanchi), where a recent significant earthquake (Chi-Chi) occurred in 1999. For all these reasons, this area is one of the most vulnerable to landslides in Taiwan.
Based on a statistical report of the IFRC, the highest number of landslides occurred in the period from 2008 to 2010 (accounting for 72 of the 177 events occurring worldwide) [1]. This period also corresponds to some very strong typhoons in Taiwan (Ex: Typhoon Kalmaegi 2008, Typhoon Morakot in 2009). As typhoons that landed on the main island are often accompanied by enormous amounts of rainfall, they often trigger shallow landslides, such as debris flows or slides, and these are two common categories of landslides in the Chen-Yu-Lan river watershed. Typhoon Kalmaegi, which formed on 13 July, brought heavy rainfall of approximately 1200 mm over a three-day period (16 July to 18 July 2009), caused a loss of lives (25 fatalities and one missing) and destruction of about 332.3 million USD of property [31]. Typhoon Morakot, which formed on 2 August and dissipated on 12 August 2009, hitting Taiwan from the 8 August to 9 August, brought much heavier rainfall of 3000 mm caused a greater loss of lives (677 deaths, four severely injured, 22  Kalmaegi: 120 km/h; Typhoon Morakot: 140 km/h), the Regional Specialized Meteorological Center Tokyo (RSMC Tokyo) ranked Typhoon Kalmaegi and Typhoon Morakot in the typhoon group of the Tropical Cyclone Intensity Scale. However, the impacts (damages and fatalities) of these two typhoons in Taiwan are very different (as mentioned above) because of the huge difference in their rainfall (Typhoon Kalmaegi: 1200 mm; Typhoon Morakot: approximate 3000 mm). Due to this difference, in this study, we considered Typhoon Kalmaegi as a normal event and Typhoon Morakot is an extreme event which triggered landslides. Landslide data related to these events are thus used to generate a landslide susceptibility map and evaluate the accuracy of the predictions of our combined model. Morakot in the typhoon group of the Tropical Cyclone Intensity Scale. However, the impacts (damages and fatalities) of these two typhoons in Taiwan are very different (as mentioned above) because of the huge difference in their rainfall (Typhoon Kalmaegi: 1200 mm; Typhoon Morakot: approximate 3000 mm). Due to this difference, in this study, we considered Typhoon Kalmaegi as a normal event and Typhoon Morakot is an extreme event which triggered landslides. Landslide data related to these events are thus used to generate a landslide susceptibility map and evaluate the accuracy of the predictions of our combined model.  Figure 2 shows the flowchart of this research, the types of data and the steps to conduct the study are as follows:

Producing Landslide Inventory Maps Related to Events
Unclouded Formosat-2 satellite images taken from 2008/06/10 to 2008/06/22 (the period before Typhoon Kalmaegi occurred), from 2008/08/24 to 2009/02/03 (the period after Typhoon Kalmaegi), and from 2009/10/18 to 2009/10/21 (the period after Typhoon Morakot) were used to produce landslide inventory maps (LMs) related to the event (Typhoon Kalmaegi), which are called, in turn, the pre-event LM and the post-event LM , using the Formosat-2 Automatic Image Processing System (F-2 AIPS) [32][33][34]. The product of this processing is called a spectral summation intensity modulation (SSIM) image. The SSIM images thus obtained were combined with the digital elevation model (DEM) data to produce 2D ( Figure 3a) and 3D (Figure 3b) composite images [32]. The DEM is presented as a raster graphics image with a resolution of 20 meters, and is able to show the height information of the surface of the study area, which was used to assist in the preparation of the landslide inventory [34] (Figure 3a [32][33][34]. The product of this processing is called a spectral summation intensity modulation (SSIM) image. The SSIM images thus obtained were combined with the digital elevation model (DEM) data to produce 2D ( Figure 3a) and 3D (Figure 3b) composite images [32]. The DEM is presented as a raster graphics image with a resolution of 20 m, and is able to show the height information of the surface of the study area, which was used to assist in the preparation of the landslide inventory [34] (Figure 3a

Producing Landslide Inventory Maps Related to Events
Unclouded Formosat-2 satellite images taken from 2008/06/10 to 2008/06/22 (the period before Typhoon Kalmaegi occurred), from 2008/08/24 to 2009/02/03 (the period after Typhoon Kalmaegi), and from 2009/10/18 to 2009/10/21 (the period after Typhoon Morakot) were used to produce landslide inventory maps (LMs) related to the event (Typhoon Kalmaegi), which are called, in turn, the pre-event LM and the post-event LM , using the Formosat-2 Automatic Image Processing System (F-2 AIPS) [32][33][34]. The product of this processing is called a spectral summation intensity modulation (SSIM) image. The SSIM images thus obtained were combined with the digital elevation model (DEM) data to produce 2D ( Figure 3a) and 3D (Figure 3b) composite images [32]. The DEM is presented as a raster graphics image with a resolution of 20 meters, and is able to show the height information of the surface of the study area, which was used to assist in the preparation of the landslide inventory [34] (Figure 3a    The Expert Landslide and Shaded Area Delineation System (ELSADS) was then used, along with the differentiation thresholds of the Normalized Difference Vegetation Index (NDVI) and the normalized green red difference indices (NGRDI) to automatically detect and localize vegetation coverage and non-coverage areas, as shown by the green lines in Figure 3a,b [34,35]. The resulting flexible, composite 3D images ( Figure 3b) can enable users to easily recognize and re-detect actual landslides, as well as erase buildings, roads, rivers, farms (vegetation non-coverage areas) and so on to detect the real landslides, as shown by the yellow lines in Figure 3a,b. This work was used to produce event-based landslide inventory maps in the whole study area (Figure 4). Most of the landslides in the Chen-Yu-Lan watershed are shallow ones that are known as debris flows and debris slides. Overlapping two kinds of LMs (the pre-event LM and the post-event LM) generated the during-event LM (this LM showed new landslide positions which occurred during the event) ( Figure 4). The Expert Landslide and Shaded Area Delineation System (ELSADS) was then used, along with the differentiation thresholds of the Normalized Difference Vegetation Index (NDVI) and the normalized green red difference indices (NGRDI) to automatically detect and localize vegetation coverage and non-coverage areas, as shown by the green lines in Figure 3a,b [34,35]. The resulting flexible, composite 3D images ( Figure 3b) can enable users to easily recognize and re-detect actual landslides, as well as erase buildings, roads, rivers, farms (vegetation non-coverage areas) and so on to detect the real landslides, as shown by the yellow lines in Figure 3a,b. This work was used to produce event-based landslide inventory maps in the whole study area (Figure 4). Most of the landslides in the Chen-Yu-Lan watershed are shallow ones that are known as debris flows and debris slides. Overlapping two kinds of LMs (the pre-event LM and the post-event LM) generated the during-event LM (this LM showed new landslide positions which occurred during the event) ( Figure 4).

Building a Hierarchy Tree of Landslide Factors
To build a hierarchical tree, factors (or variables) were classified into components (classes) based on experts' experience. The DEM in 20 m resolution ( Figure 5a) is supported by the Forestry Bureau, Council of Agriculture in Taiwan was used to extract slopes (SL), drainage densities (DD),

Building a Hierarchy Tree of Landslide Factors
To build a hierarchical tree, factors (or variables) were classified into components (classes) based on experts' experience. The DEM in 20 m resolution ( Figure 5a) is supported by the Forestry Bureau, Council of Agriculture in Taiwan was used to extract slopes (SL), drainage densities (DD), and aspects (AS) (Figure 5d-f). Besides, a geological map (GM) (scale 1:250,000 was used to measure the strength of the material. Based on the information about the rocks' resistances and distributions, the lithology in this area was divided into six main groups, namely deposits, alluvium, sandstone-shale, argillite-slate, quartzite slate, and phillite-slate ( Figure 5b). This GM data was then converted to the raster data in 20 m resolution. A land-use map (LU), which was published by Council of Agriculture and Central Geological Survey in 2006, was used to determine the coverage of the surfaces. This data was processed by SPOT 5 satellite images which were taken in 2004. The LU is classified into six major groups, including agriculture, forest, road and building, river body and wetland, grassland, and bare land (Figure 5c). The LU data was also converted to the raster data in 20-m resolution. and aspects (AS) (Figure 5d-f). Besides, a geological map (GM) (scale 1:250,000 was used to measure the strength of the material. Based on the information about the rocks' resistances and distributions, the lithology in this area was divided into six main groups, namely deposits, alluvium, sandstone-shale, argillite-slate, quartzite slate, and phillite-slate ( Figure 5b). This GM data was then converted to the raster data in 20 m resolution. A land-use map (LU), which was published by Council of Agriculture and Central Geological Survey in 2006, was used to determine the coverage of the surfaces. This data was processed by SPOT 5 satellite images which were taken in 2004. The LU is classified into six major groups, including agriculture, forest, road and building, river body and wetland, grassland, and bare land (Figure 5c). The LU data was also converted to the raster data in 20-m resolution.

Computing the Landslide Density of Factors
Instead of based on expert opinion analyses to assess the important level of class pairs in each factor, we used the landslide density of each class to compare and calculate this important level. The landslide density of class i (Di) is calculated in Equation (1) with the pre-event LM data, is ratio of landslide area in class i (Li, pixel unit) and distribution area in class i (Ai, pixel unit).

Computing the Landslide Density of Factors
Instead of based on expert opinion analyses to assess the important level of class pairs in each factor, we used the landslide density of each class to compare and calculate this important level. The landslide density of class i (D i ) is calculated in Equation (1) with the pre-event LM data, is ratio of landslide area in class i (L i , pixel unit) and distribution area in class i (A i , pixel unit).

Calculating the Weight of Factors by Combining of AHP and Bivariate Analysis
Then, this comparison is made between pairs of indicators (classes) and is combined into a matrix of n lines and n columns (n is the number of indicators or classes in each factor) ( Figure 6). The element D ij represents the important level of the index i (line) versus the index j (column). The relative importance of indicator i versus j is calculated in D i /D j ratio, and vice versa of indicator j versus i is D j /D i , and D ii = 1 ( Figure 6). Each matrix has been set up for each factor that can set the priorities of the elements (classes) in the hierarchy tree. The priority is a nonnegative number that fluctuates from 0 to 1 representing the weight association in each element at each level. By definition, the priority of the target is 1, and the total priority of a level is also 1. In this study, the priority is also the weight (W) of each class, and computed in Equation (2), where n is number of classes in each factor, D i is the landslide density of class i (Equation (1)), and W i is the weight of class i.

Calculating the Weight of Factors by Combining of AHP and Bivariate Analysis
Then, this comparison is made between pairs of indicators (classes) and is combined into a matrix of n lines and n columns (n is the number of indicators or classes in each factor) ( Figure 6). The element Dij represents the important level of the index i (line) versus the index j (column). The relative importance of indicator i versus j is calculated in Di/Dj ratio, and vice versa of indicator j versus i is Dj/Di, and Dii = 1 (Figure 6). Each matrix has been set up for each factor that can set the priorities of the elements (classes) in the hierarchy tree. The priority is a nonnegative number that fluctuates from 0 to 1 representing the weight association in each element at each level. By definition, the priority of the target is 1, and the total priority of a level is also 1. In this study, the priority is also the weight (W) of each class, and computed in Equation (2), where n is number of classes in each factor, Di is the landslide density of class i (Equation (1)), and Wi is the weight of class i.
After executing steps of calculating the weights of classes, one question was set forth whether any method evaluates the validity of the significance values of the indicators (classes). According to Saaty [36], the consistency ratio (CR) can be used and calculated in Equations (3)-(5), where CI is the consistent index, RI is the random number, and n is number of classes in each factor. This ratio compares the consistency with the random objectivity of the data. The consistency ratio (CR) of less than 0.1 is acceptable. If it is greater than 0.1, the decision maker needs to reduce the heterogeneity by changing the significance level between the pair of indicators (classes).

 
For each of the n-level comparison matrices, Saaty [36] tested the creation of random matrices and calculated the RI (random index) corresponding to matrix levels as Table 1.
After executing steps of calculating the weights of classes, one question was set forth whether any method evaluates the validity of the significance values of the indicators (classes). According to Saaty [36], the consistency ratio (CR) can be used and calculated in Equations (3)-(5), where CI is the consistent index, RI is the random number, and n is number of classes in each factor. This ratio compares the consistency with the random objectivity of the data. The consistency ratio (CR) of less than 0.1 is acceptable. If it is greater than 0.1, the decision maker needs to reduce the heterogeneity by changing the significance level between the pair of indicators (classes).
For each of the n-level comparison matrices, Saaty [36] tested the creation of random matrices and calculated the RI (random index) corresponding to matrix levels as Table 1.

Finding Out Coefficients of Landslide Factors by a Combination of AHP and Correlation Model
In order to measure and compare the importance or influence level of pairs of factors on landslides, the correlation and dependency model is used. The weighting correlation of landslide factors is analyzed on each pixel, then aggregated and output, and is called the correlation (r). According to Pearce, correlation (r) is a statistical indicator that measures correlation between two variables (x and y) [37]. The correlation (r) fluctuates from (−1) to 1. If the r equals or approximately 0, two variables have no relation, and vice versa if the r is −1 or 1, two variables have an absolute relation. The process of calculating correlation coefficient in turn is performed in the order of steps as explained in formulas 6 to 10.
Firstly, the average weight value of factor k (µ k ) must be calculated throughout the study area, as in Equation (6), where W ik is the weight value of pixel i of factor k, and n is number of pixels in the study area. Then the variance (Var k ) and standard deviation (σ k ) of factor k are calculated, as in Equations (7) and (8). The covariance (Covar kj ) between the weight values of factors k and j is calculated as in Formula (9). The variance and covariance of factors are represented by a variance-covariance matrix, and shown in a symmetric matrix ( Table 2). The correlation (r kj ) between the weighting values of the pair of factors (k & j) is calculated by the Formula (11), where σ k is the standard deviation of factor k and σ j is the standard deviation of factor j. These correlations are given the absolute value also combined into a symmetric matrix (Table 3) which is applied for calculating coefficients (C) by the AHP model (Section 2.3).

The Weight of Landslide Factors Calculated by the Combination of AHP and Bivariate Analysis Model
Low-resistance petrographic groups due to being compressed and crumbled by tectonic activities such as phyllite-slate-interbedded sandstone or quartzite-slate-coaly shale correspond to high and the highest weights, whereas high-strength rocks formed in less severe environments such as sandstone, shale, and argillite are low in weight. Although deposits and alluvium are of low-strength materials, their weights are not high. This inadequacy can be explained by the fact that these types of materials are distributed in areas with flat terrain and low gradient.
From the results of calculating the weights of the aspect factors, it is found that the south and west landslide risk is higher than the other directions. The distribution of the aspects in the study area is fairly uniform (in the 10-17% range), so the difference in weight of these classes is not as great as in the case of slopes, land use. The case of DD is similar to that of aspects that the difference in weight values of classes is not large. It is also unreasonable for smaller DD values (0.002-0.003) to have greater weights in a comparison to larger one (0.004-0.005) ( Table 2). This can be explained by two reasons: the influence of the distribution area on the weight calculation of DD classes and almost DDs 0.002-0.003 are distributed in locations with steep slopes.
The result of the computation of the consistency ratio shows that all CRs less than 0.1 (Table 5), which mean that the values of significance between pairs of indicators are acceptable.  Table 6 depicts the symmetric matrix of landslide factors' coefficients that calculated by the AHP. This shows that GM and LU are two important factors in computing LSI because their coefficients are higher than the others'. The two factors with the best correlation are LU and GM, whereas DD has a poor correlation with the other factors.

Generating the LSM by the Combination Model of the Geometric Multivariate and Factor Coefficients
A landslide susceptibility map (LSM) was produced by as in the Equation (11), where the significance of landslide factors is also the factor's coefficient (C) (Figure 7a). The LSM showed that most of landslide areas (black color) which occurred in events (Typhoon Kalmaegi and Typhoon Morakot) distributed in the highest indices (red color) (Figure 7b-e). The safe regions are distributed in low-index positions corresponding to the river channel and its neighborhood with fairly flat terrain and low gradient slopes. Figure 8 shows how to determine the thresholds of the LSM indices. Two cumulative percentage curves of landslide and non-landslide are set and viewed in the same graph, the intersection of which is the threshold (the dotted red line) between the safe zone (accounting for more than 80% of study area) and the dangerous zone (approximate 20% of study area) (Figure 8a). Areas in the danger zone are highly coincident with real landslides. After comparing and reviewing the landslide data of different years, we found that the landslides that occurred afterwards (new landslides) often occurred in adjacent locations along the old landslides. Precise warnings about areas of landslide are very significant, therefore, in order to increase the accuracy in warning, we create additional buffer zones for warning areas and then two sub-thresholds are set up (Figure 8b). Figure 8b shows the distribution percentage of LSI corresponding to actual landslides and non-landslides. This figure shows that the landslide curve is variable (at the threshold LSI (0.1458), the landslide percentage increases suddenly from less than 10% to greater than 60%), while non-landslide does not fluctuate (Figure 8b). Therefore, based on the variation of landslide curve (Figure 8b), two sub-thresholds are also defined (the dotted brown and green lines). As a result, the LSM indices were divided into four levels, including low (index: 0.0148-0.0927), low medium (0.0928-0.1161), high (0.1162-0.1457), and very high (LSI: 0.1458-0.4123), with distribution areas of 30.25%, 30.52%, 19.92%, and 19.31%, respectively (Figures 8 and 9).  Compared to the other factors in the LSI modeling calculations, the GM here does not have a high resolution (scale: 1:200,000), so the rock classification was not detailed, which are grouped together for example the argillite-slate-phyllite group (three rock categories were grouped together). In addition, because GM's landslide coefficient (or significant level for landslides) is higher than other factors' ones, high LSI indices are roughly aligned with the geological boundaries of the GM classes.
As result, the very high LSI (0. 39-0.41), which occupies about 46.11 ha, is mainly distributed in the southern part of the study area (Yu-Shan mountain peak) (Figure 7a). Although the LSI value is high, there are some positions that do not have landslide occurrences (non-landslides). This can be explained by the fact that the LSI values of these positions are dominated by three factors including slope, lithology and land-use. All of these factors have the highest weighting in above positions, thus making the LSI values very high. Although these positions are located in steep cliffs and are in the bare-land group, the exposed surface is bedrock rather than loose materials, so landslides do not occur here.
in examining the correlation between the during-event precipitation and the LSI at severe landslides. The largest landslide scar in Typhoon Morakot (111.07 hectares) was located in the southwestern, which corresponds to the highest rainfall in Morakot but LSI indices range from low (0.06) to very high (0.30). In contrast, the largest landslide scar in Typhoon Kalmaegi (49.09 hectares) was distributed in the south near Mount Yushan, corresponding to a high LSI and low-moderate rainfall. This again shows that the precipitation factor in extreme events is extremely important, because it can cause landslides in all indices from low to very high.

Validating the LSM with Types of Related-Typhoon Kalmaegi LMs
The binary classification and Kappa index (K) methods were used to validate the LSM. The "binary classification" method is a method of assessing the accuracy of interpretative models, which is widely used in research topics related to statistics in the fields of economics, medicine, or GIS [7]. The binary classification is a model, in which objects were classified into two groups based on its attributes. In this research, the attributes were classified in particular are landslide occurrence. Therefore, the LSI is divided into two groups consisting of dangerous and safety groups. A threshold of division between danger and safety is established.
Here, the correct interpretation is called accuracy, which is the sum of TP and TN. TP is the area of real landslides corresponding to LSI area that is divided into the dangerous group. TN is non-landslide areas and corresponds to LSIs that are partitioned in safe areas. In contrast, the error is the sum of FN and FP, where FN is the area to be interpreted as having a safe LSI (safe group) but with landslide occurrences, and FP is the area to be interpreted as a high risk of landslides (dangerous group) but no landslide occurs. The study area is a total of TP, TN, FP, and FN (Equation (13)).
The Kappa index method developed by Cohen Kappa (1960) is used to evaluate the similarity between the statements. Specifically, in this study, there was a claim of landslide or no landslide. According to Cohen Kappa, the main interpretations of K can be divided into the following: (K < 0.  Table 7 showed that the prediction accuracy of the LSM is high if overlapping to the normal events (Typhoon Kalmaegi: the post-event 77.35%, the during-event 70.26%; Typhoon Morakot: the post-event 65.64%, the during-event 58.49%). The calculating K index results are also the consensus of predictions of Typhoon Kalmaegi (post-event: K = 0.559; during-event: K = 0.412) was higher than Typhoon Morakot (post-event: K = 0.313, during-event: K = 0.170). These K values are at moderate consensus. This can be explained by the fact that the rainfall of these two events is very different. Kalmaegi is classified as a normal event (maximum total of precipitation: 774 mm), while Morakot (maximum total of precipitation: 1974 mm) is classified as an extreme event. To increase interpreting agreement, a rainfall factor should be considered to add into the computation model. Most of the previous studies [21,38], researchers have often considered rainfall as one of the landslide factors and they add it to the other LSI computational model. However, in our project, finding thresholds between the rainfall and the LSI indices is very significant in the warning of landslides before the storm occurs. Therefore, the rainfall is considered as one of the trigger factors, so this parameter is not included in this computational model (Phase 1), it will be considered in phase 2. In the upcoming work, the relationship between the rainfall and LSI indices will be analyzed from which to find their corresponding thresholds (Phase 2). This work is thus done in phase 2. Besides, success rate curves have been used in many studies to assess the validity of predictions produced under different scenarios, and they have also been applied in some works that obtained LSI maps using different factors [39][40][41][42]. The "success rate curve" is a method which has been used to measure the accurate proportion of prediction by plotting cumulative percentages of high LSI risk and real landslide during typhoon. Therein, the data of overlapping LSI values and actual landslides were analyzed, and based on the resulting statistics, the cumulative percentages of LSI and landslide were calculated. These results were then plotted on the same graph, known as a "success rate curve". By observing the distributions of the curves in the graph, it is possible to assess the accuracy of the related predictions. The LSM's success rate curves which are evaluated by events (Typhoon Kalmaegi and Typhoon Morakot) showed a high accuracy of interpretation rates with approximately 86% of the real landslides (post Typhoon Kalmaegi) occurring in areas with high and very high LSI (dangerous distribution area: 39%), and over 75% of new landslides occurring during this event. Also, 61.73% of actual landslide areas occurred post Typhoon Morakot and 49.37% of new landslide positions also occurred during the typhoon ( Figure 9, Table 8). These curves once again show the importance of rainfall parameters in the LSI calculation model for extreme events in the comparison of the results between Typhoon Kalmaegi and Typhoon Morakot. The above conclusions are further strengthened in examining the correlation between the during-event precipitation and the LSI at severe landslides. The largest landslide scar in Typhoon Morakot (111.07 hectares) was located in the southwestern, which corresponds to the highest rainfall in Morakot but LSI indices range from low (0.06) to very high (0.30). In contrast, the largest landslide scar in Typhoon Kalmaegi (49.09 hectares) was distributed in the south near Mount Yushan, corresponding to a high LSI and low-moderate rainfall. This again shows that the precipitation factor in extreme events is extremely important, because it can cause landslides in all indices from low to very high.

Conclusions
Producing a disaster prediction map (landslide susceptibility map) is one of the important tasks that disaster researchers are interested in and need to carry out. They always try to develop statistical analysis models so that they can have the most accurate predictions of landslides in the future. Landslide researchers have constructed many models for landslide prediction by various methods, such as multivariate statistics, binary statistics, the AHP, neutral networks, etc. Especially, in using the AHP method for landslide interpretation, researchers often use the statistics of experts' opinion and experience to establish weights or compare the importance for pairwise indicators (classes of each landslide factor).
The purpose of this paper is to present a new approach and via AHP modeling to generate landslide prediction maps (LSMs). In this study, a combination of the AHP, bivariate, and correlation statistics model was developed. Instead of using experts' opinions to build a weighting matrix for the classes of factors, the bivariate analysis method (considering the landslide density of each class) is used to analyze this weight. In addition, correlation statistics were also used to evaluate the correlation coefficients of the landslide factors (C) by the AHP model. Finally, a multivariate analysis model incorporating the correlation coefficients of the landslide factors was applied to produce the LSM map.
To assess the accuracy of the above integration model, the LM data (the post-event landslide inventory maps, the during-event landslide inventory maps) associated with the events was used to evaluate the binary classification method, Kappa index, and success rate curves. The results showed that for normal events, such as Kalmaegi, LSM's success rate is high (over 77% for post-event, and over 70% for during-event), whereas for extreme events, such as Morakot, it is less accurate (over 65% for post-event, and approximately 60% for during-event). This assessment proposed that further consideration of the rainfall data in landslide predictive models is needed in future work. Therefore, in the forthcoming work, we will continue to use this approach in developing a model for calculating landslide hazard indices, where the correlation coefficient between LSI (built on preparatory factors) and rainfall (a triggering factor) will be considered in the calculation model.