Scale and spatial distribution assessment of rainfall-induced landslides in a catchment with mountain roads

This study focused on landslides in a catchment with mountain roads that were caused by Nanmadol (2011) and Kong-rey (2013) typhoons. Image interpretation techniques were employed to for satellite images captured before and after the typhoons to derive the surface changes. A multivariate hazard evaluation method was adopted to establish a landslide susceptibility assessment model. The evaluation of landslide locations and relationship between landslide and predisposing factors is preparatory for assessing and mapping landslide susceptibility. The results can serve as a reference for preventing and mitigating slope disasters on mountain roads.

1 Introduction

Taiwan is an island with three quarters of its land area consisting of slope land that is 100 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$, or less but has an average gradient of 5 % above (SWCB, 2017). Much of this sloped land has a steep gradient and fragile geological formations. Taiwan is hit by an average of 3.4 typhoons every year during the years 1911 to 2016 (Central Weather Bureau, 2017). In addition, average annual rainfall reaches 2502 mm in the years 1949 to 2009 (Water Resources Agency, 2017). Typhoons usually occur between July and October, and 70–90 % of the annual rainfall is composed of heavy rain directly related to typhoons (SWCB, 2017). Concentrated rainfall causes heavy landslides and debris flows every year (Dadson et al., 2004). The threat of disaster currently influences industrial and economic development and the road networks in endangered areas, thus establishing disaster evaluation mechanisms is imperative.

Landslide susceptibility can be evaluated by analysing the relationships between landslides and various factors that are responsible for the occurrence of landslides (Brabb, 1984; Guzzetti et al., 1999, 2005). In general, the factors that affect landslides include predisposing factors (e.g. geology, topography, and hydrology) and triggering factors (e.g. rainfall, earthquakes, and anthropogenic factors) (Chen et al., 2013a, b; Chue et al., 2015). Geological factors include lithological factors, structural conditions, and soil thickness; topographical factors include slope, aspect, and elevation; and anthropogenic factors include deforestation, road construction, land development, mining, and alterations of surface vegetation (Chen et al., 2013a, b; Chue et al., 2015). The method used to assess landslide susceptibility can be divided into qualitative and quantitative. Qualitative methods are based completely on field observations and an expert's prior knowledge of the study area (Stevenson, 1977; Anbalagan, 1992; Gupta and Anbalagan, 1997). Some qualitative approaches incorporate ranking and weighting, and become semi-quantitative (Ayalew and Yamagishi, 2005). For example the analytic hierarchy process (AHP) (Saaty, 1980; Barredo et al., 2000; Yoshimatsu and Abe, 2006; Kamp et al., 2008; Yalcin, 2008; Kayastha et al., 2013; Zhang et al., 2016) and the weighted linear combination (WLC) (Jiang and Eastman, 2000; Ayalew et al., 2005; Akgün et al., 2008). Quantitative methods apply mathematical models to assess the probability of landslide occurrence, and thus define hazard zones on a continuous scale (Guzzetti et al., 1999). Quantitative methods developed to detect the areas prone to landslide can be divided mainly into two categories: deterministic approach and statistical approach. The deterministic approach is based on the physical laws driving landslides (Okimura and Kawatani, 1987; Hammond et al., 1992; Montgomery and Dietrich, 1994; Wu and Sidle, 1995; Gökceoglu and Aksoy, 1996; Pack et al., 1999; Iverson, 2000; Guimarães et al., 2003; Xie et al., 2004) and are generally more applicable when the underground conditions are relatively homogeneous and the landslides are mainly slope dominated. The statistical approach is based on the relationships between the affecting factors and past and present landslide distribution (Van Westen et al., 2008). Statistical methods analyse the relation between predisposing factors affecting the landslide which include bivariate statistical models (Van Westen et al., 2003; Süzen and Doyuran, 2004; Thiery et al., 2007; Bai et al., 2009; Constantin et al., 2011; Yilmaz et al., 2012), multivariate statistical approaches as discriminant analysis (Baeza and Corominas, 2001; Carrara et al., 2003, 2008; Pellicani et al., 2014), and linear and logistic regression (Dai and Lee, 2002; Ohlmacher and Davis, 2003; Ayalew and Yamagishi, 2005; Yesilnacar and Topal, 2005; Greco et al., 2007; Carrara et al., 2008; Lee et al., 2008; Pellicani et al., 2014), as well as non-linear methods such as artificial neural networks (ANN) (Lee et al., 2004; Yesilnacar and Topal, 2005; Kanungo et al., 2006; Wang and Sassa, 2006; Li et al., 2012) and multivariate hazard evaluation method (MHEM) (Su et al., 1998; Lin et al., 2009). The MHEM is a non-linear mathematical model that presents an instability index to indicate landslide susceptibility (Lin et al., 2009). In addition, in some studies, landslide susceptibility analyses have focused on man-made facilities such as roads and railroads and have examined the landslide susceptibility of the surrounding environments (Das et al., 2010, 2012; Pantelidis, 2011; Devkota et al., 2013; Martinović et al., 2016; Pellicani et al., 2016, 2017). The aforementioned studies on the landslide susceptibility of areas surrounding man-made facilities have not investigated characteristics such as the location and scale (area) of landslides occurring in upper or lower slopes, and these thus constitute one of the objectives of the present study.

Technological progress has provided various advanced tools and techniques for land use monitoring. In recent years, aerial photos or satellite images have been commonly used in post-disaster interpretations and assessments of landslide damage on large-area slopes (Erbek et al., 2004; Lillesand et al., 2004; Nikolakopoulos et al., 2005; Lin et al., 2005; Chen et al., 2009; Otukei and Blaschke, 2010; Chen et al., 2013a). Satellite images offer the advantages of short data acquisition cycles, swift understanding of surface changes, large data ranges, and being low cost, particularly for mountainous and inaccessible areas. With the assistance of computer analysis and geographic information system (GIS) platforms, researchers can quickly determine land cover conditions. Thus, satellite images are suitable for investigating large areas and monitoring temporal changes in land use (Liu et al., 2001). Satellites can capture images of the same area multiple times within a short period. Studies have indicated that land surface change detection is the process of exploring the differences between images captured at different times (Liu et al., 2001; Chadwick et al., 2005; Chen et al., 2009; Chue et al., 2015). With multispectral satellite images, land surface interpretations involve comparisons of multitemporal images that are completely geometrically aligned.

We selected part of the catchment area of Laonung River which include Provincial Highway 20 in southern Taiwan as our study area. Regarding time, we focused on periods before and after landslides that occurred in the study area as a result of Typhoon Nanmadol (2011) and Typhoon Kong-rey (2013). We applied the maximum likelihood method to interpret and categorize high-resolution satellite images, thereby determining the land surface changes and landslides in the study area before and after the rainfall events. By using a GIS platform, we constructed a database of the rainfall and natural environment factors. Subsequently, we developed a landslide susceptibility assessment model by using the MHEM. The model performance was then verified by historical landslides. In addition, we extracted the locations of landslide areas to explore the relationship between the natural environment and the spatial distribution of the scale of these areas.

2 Methodology

2.1 Maximum likelihood

The maximum likelihood classifier is a supervised classification method (SCM). SCMs include three processing stages: training data sampling, classification, and output. The underlying principle of supervised classification is the use of spectral pattern recognition and actual ground surface data to determine the types of data required and subsequently select a training site, which has a unique set of spectral patterns. To accurately estimate the various spectral conditions, the spectral patterns of the same type of feature are combined into a coincident spectral plot before the class of the training site is selected. Once training has been completed, the entire image is classified based on the spectral distribution characteristics of the training site by using statistical theory for automatic interpretation (Lillesand et al., 2004).

To facilitate the calculation of probability in the classification of unknown pixels, the maximum likelihood method assumes a normal distribution in the various classes of data. Under this assumption, the data distribution can be expressed using covariance matrices and mean vectors, both of which are used to calculate the probability of a pixel being assigned to a land cover class. In other words, the probability of X appearing in class i is calculated using Eq. (1), and the highest probability is used to determine the feature of each pixel (Lillesand et al., 2004).

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }p\left(X|C_i\right)=\\ \text{(1)}& {\displaystyle }& {\displaystyle }{\left(\mathrm{2}\mathit{\pi }\right)}^{d/\mathrm{2}}{\left|\sum _i\right|}^{-\mathrm{1}/\mathrm{2}}\exp \left[-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\left(X-{\mathit{\mu }}_i\right)}^T\sum _i^{-\mathrm{1}}\left(X-{\mathit{\mu }}_i\right)\right]{\displaystyle }& {\displaystyle }\mathrm{X}=\left[\begin{array}{c}{\mathrm{X}}_{\mathrm{1}}\\ {\mathrm{X}}_{\mathrm{2}}\\ \mathrm{⋮}\\ {\mathrm{X}}_d\end{array}\right]{\mathit{\mu }}_i=\left[\begin{array}{c}{\mathit{\mu }}_{\mathrm{1}i}\\ {\mathit{\mu }}_{\mathrm{2}i}\\ \mathrm{⋮}\\ {\mathit{\mu }}_{\text{di}}\end{array}\right]\sum _i=\left[\begin{array}{cccc}{S}_{\mathrm{11}}& S_{\mathrm{12}}& \mathrm{\cdots }& S_{\mathrm{1}d}\\ S_{\mathrm{21}}& S_{\mathrm{22}}& \mathrm{\cdots }& S_{\mathrm{2}d}\\ \mathrm{⋮}& \mathrm{⋮}& \mathrm{\ddots }& \mathrm{⋮}\\ S_{d\mathrm{1}}& S_{d\mathrm{2}}& \mathrm{\cdots }& S_{\text{dd}}\end{array}\right]\end{array}$$

In this equation, d denotes the number of features, X denotes a sample expressed using features and has d dimensions, p(X|C_i) denotes the probability that X originates from class i, Σ_i denotes the covariance matrix of class i, ${\mathrm{\Sigma }}_i^{-\mathrm{1}}$ denotes the inverse matrix of Σ_i, $\left|{\mathrm{\Sigma }}_i\right|$ denotes the determinant of Σ_i, μ_i denotes the mean vector of classification i, (X−μ_i)^T denotes the transpose matrix of (X−μ_i), and S_ij denotes the covariance of classes i and j.

During classification, the maximum value of the probability density functions of sample X in each class is used to determine which class the sample belongs to. The maximum likelihood classification decision is shown in Eq. (2).

$${\displaystyle }{\displaystyle }X\in C_m,m\subset \left\{\mathrm{1},\mathrm{2},\mathrm{\cdots },k\right\}$$

if

$$\begin{array}{}\text{(2)}& p\left(\left.X\right|C_m\right)=max\left\{p\left(X|C_j\right),j=\mathrm{1},\mathrm{2},\mathrm{\cdots }k\right\},\end{array}$$

in which k denotes the number of classes. The question regarding classification is how to effectively separate the classes in the feature space, or in other words, how to divide the feature space. Maximum likelihood is a common approach that offers fairly good classification accuracy (Bruzzone and Prieto, 2001; Chen et al., 2004). Thus, we adopted maximum likelihood to interpret and classify the satellite images.

2.2 Accuracy assessment

This study employed the aforementioned maximum likelihood method to classify satellite images. To determine whether the accuracy of image classification was acceptable, we adopted an error matrix to test for accuracy. An error matrix is a square matrix that presents error conditions in the relationship between ground surface classification results and reference data (Verbyla, 1995). It contains an equal number of columns and rows, and the number is determined by the number of classes. For example, Table 1 contains four classes. The columns show the reference data, and the rows show the classification results. The various elements in the table indicate the quantity of data corresponding to each combination of classes.

Table 1 Relationship table of error matrix (Verbyla, 1995).

In the Table 1, X₁₂ represents the amount of data that were interpreted as Class A but actually belong to Class B, whereas X₂₁ indicates the amount of data that were interpreted as Class B but actually belong to Class A. X₁₁ and X₂₂ represent the amount of data accurately classified as Class A and Class B. An error matrix is generally used to check the quality of classification results in statistics (Congalton, 1991; Verbyla, 1995). In the present study, we evaluated the accuracy of the classification results based on the overall accuracy and kappa value (Cohen, 1960), which is the coefficient of agreement derived from the relationship between the classification results and training data. These two parameters are explained as follows.

2.2.1 Overall accuracy (OA)

OA is the simplest method of overall description. For all classes, OA represents the probability that any given point in the area will be classified correctly.

$$\begin{array}{}\text{(3)}& {\displaystyle }{\displaystyle }\text{OA}=\left[{\displaystyle \frac{\mathrm{1}}{N}}{\sum }_{i=\mathrm{1}}^n{X}_{\text{ii}}\right]\times \mathrm{100}\mathit{\%}\end{array}$$

In Eq. (3), N denotes the total number of classifications, n denotes the total number of rows in the matrix, and X_ii is the number of correctly classified checkpoints.

2.2.2 Kappa coefficient

The kappa ($\widehat{K}$) coefficient indicates the degree of agreement between the classification results and reference values and shows the percentage reduction in the errors of a classification process compared with the errors of a completely random classification process. Generally, the kappa coefficient ranges from 0 to 1, and a greater value indicates a higher degree of agreement between the two sets of results, as shown in Eq. (4):

$$\begin{array}{}\text{(4)}& {\displaystyle }{\displaystyle }\widehat{K}={\displaystyle \frac{N{\sum }_{i=\mathrm{1}}^n{X}_{\text{ii}}-{\sum }_{i=\mathrm{1}}^n\left(X_{i+}\times X_{+i}\right)}{N^{\mathrm{2}}-{\sum }_{i=\mathrm{1}}^n\left(X_{i+}\times X_{+i}\right)}}\times \mathrm{100}\mathit{\%},\end{array}$$

in which X_i+ is the total number of pixels for a given class on the actual ground surface and X_+i is the number of pixels in that class. As reported by Landis and Koch (1977), a kappa coefficient greater than 0.8 signifies a high degree of accuracy, whereas a coefficient between 0.4 and 0.8 or less than 0.4 indicates moderate or poor accuracy.

2.3 Rainfall analysis method

In previous studies regarding the influence of rainfall on landslides, rainfall intensity and accumulated rainfall have been most commonly used as predisposing factors of landslides (Giannecchini, 2006; Chang et al., 2007; Giannecchini et al., 2012; Ali et al., 2014). Therefore, we adopted effective accumulated rainfall and intensity of rolling rainfall as rainfall indices and predisposing factors of landslides in the present study. These two indices are explained as follows.

2.3.1 Effective accumulated rainfall (EAR)

Generally, rainfall is considered the trigger of slope collapse, whereas previous rainfall can be regarded as a potential factor of a landslide. Previous rainfall influences the water content of the soil, which in turn affects the amount of rainfall required to trigger a landslide (Seo and Funasaki, 1973).

Figure 1 shows an illustration of rainfall events defined based on EAR (Seo and Funasaki, 1973). The diagram shows a concentrated rainfall event with no rainfall in the preceding or subsequent 24 h; thus it can be considered a continuous rainfall event. A continuous rainfall event that occurs simultaneously with a landslide is the main rainfall event. The beginning of the main rainfall event is defined as the time point when the rainfall first reaches 4 mm. The calculation of accumulated rainfall ends at the time when the landslide occurs. However, because the exact time of a landslide cannot be precisely determined, we regarded the hour with the maximum rainfall during the main rainfall event as the time at which the landslide occurred in this study.

In accordance with previous studies, we defined EAR as the sum of direct and previous indirect rainfall. Previous indirect rainfall is the rainfall accumulated during the 7 days prior to the main rainfall event and can be expressed as follows (Seo and Funasaki, 1973; Crozier and Eyles, 1980):

$$\begin{array}{}\text{(5)}& {\displaystyle }{\displaystyle }{\sum }_{n=\mathrm{1}}^{\mathrm{7}}{k}^n{P}_n=P_b,\end{array}$$

where P_b denotes the previous indirect rainfall, P_n denotes the rainfall during the n days prior to the main rainfall event (mm), and k denotes a diminishing coefficient set as 0.9 in this study (Chen et al., 2005). Direct rainfall encompasses the continuous rainfall accumulated during the rainfall events, starting from the first rainfall to the time of landslide occurrence. Direct rainfall has a direct and effective impact on landslide occurrence and is thus not diminished. Therefore, EAR could be expressed as follows in this study:

$$\begin{array}{}\text{(6)}& {\displaystyle }{\displaystyle }\text{EAR}=P_r+P_b,\end{array}$$

where P_r (mm) represents the rainfall accumulated during the main rainfall event from the first rainfall to the time of landslide occurrence, and P_b (mm) represents the previous indirect rainfall.

Figure 1 Definition of rainfall events based on effective accumulated rainfall (modified from Seo and Funasaki, 1973).

2.3.2 Intensity of rolling rainfall (I_R)

Rainfall intensity refers to the amount of rainfall within a unit of time. It is considered a crucial index for evaluating disasters because greater intensity or longer durations have considerable impacts on slope stability. Furthermore, rainfall-induced landslides may be triggered by several hours of continuous rainfall. The raw rainfall data in this study were hourly precipitation; thus I_R (mm h⁻¹) can be expressed as follows:

$$\begin{array}{}\text{(7)}& {\displaystyle }{\displaystyle }{I}_{m\text{R}}={\sum }_{t-m+\mathrm{1}}^{m}I=I_{t-m+\mathrm{1}}+I_{t-m+\mathrm{2}}+\mathrm{\cdots }+I_t,\end{array}$$

where I denotes rainfall intensity, m denotes the number of rolling hours of rainfall (set as 3 h in this study), I_mR denotes the I_R during m hours, and I_t denotes the rainfall intensity during hour t.

2.4 MHEM

The MHEM is a diverse non-linear mathematical model. Based on relative relationships, the MHEM presents an instability index (D_t) to indicate susceptibility in different areas. The objective is to estimate the variance of predisposing factors and then to determine the weight of each factor according to the value of variance, finally to derive a suitable landslide susceptibility assessment model (Su et al., 1998; Lin et al., 2009; Chue et al., 2015).

The predisposing factors in the MHEM are rated based on the frequency of landslide occurrence, which is calculated as follows:

$$\begin{array}{}\text{(8)}& {\displaystyle }{\displaystyle }{R}_i={\displaystyle \frac{r_i}{r_{\mathrm{T}}}},\end{array}$$

where R_i represents the landslide pixel ratio of the various factors in class i, r_i represents the number of landslide pixels in class i, and r_T represents the total number of pixels. Thus, landslide percentage X_i is expressed as

$$\begin{array}{}\text{(9)}& {\displaystyle }{\displaystyle }{X}_i={\displaystyle \frac{R_i}{\sum R_i}},\end{array}$$

where X_i denotes the landslide percentage of class i and ΣR_i denotes the sum of the landslide pixel ratios.

Based on the landslide percentages of the various classes for each predisposing factor, the normalized score value of classes for each factor (d_n) can be calculated using Eq. (10), and presented in relative values ranging from 1 to 10.

$$\begin{array}{}\text{(10)}& {\displaystyle }{\displaystyle }{d}_n={\displaystyle \frac{\mathrm{9}\left(X_i-X_{min}\right)}{\left(X_{max}-X_{min}\right)}}+\mathrm{1}\end{array}$$

In Eq. (10), X_i represents the causal rate of the sample region, and X_max and X_min represent the maximum and minimum landslide percentages of the factor in the various sample regions.

To estimate the weight of influence of each predisposing factor, the coefficient of variation (V) of the landslide ratios derived from the class of the predisposing factors is used to represent the sensitivity of landslide ratios in different predisposing factor classes. A smaller coefficient of variation denotes higher similarity among the landslide probabilities in the various classes, which indicates that this factor grading method cannot determine which areas have higher or lower landslide probabilities. By contrast, a greater coefficient of variation denotes that this factor grading method can be used to describe the influence of factor classes on landslides. Thus, the coefficient of variation among the predisposing factors can indicate the factor weights. The coefficient of variation is calculated as shown in Eq. (11):

$$\begin{array}{}\text{(11)}& {\displaystyle }{\displaystyle }V={\displaystyle \frac{\mathit{\sigma }}{X}}\times \mathrm{100}\mathit{\%},\end{array}$$

where σ is the standard deviation and X is the mean landslide percentage of the various factor classes.

We divided the coefficient of variation of each individual factor by the total coefficient of variation of all factors to derive the factor weight, which represented the degree of influence of the factor on landslide occurrence. The factor weight can be calculated as shown in Eq. (12), where W is the factor weight and V is a coefficient of variation.

$$\begin{array}{}\text{(12)}& {\displaystyle }{\displaystyle }{W}_i={\displaystyle \frac{V_i}{V_{\mathrm{1}}+V_{\mathrm{2}}+\mathrm{\cdots }+V_n}}\end{array}$$

Finally, the weight (W_i) of each factor is determined by the rank of its variance (V), and each factor is assigned a different weight. Subsequently, a non-linear mathematical model can be derived as follows:

$$\begin{array}{}\text{(13)}& {\displaystyle }{\displaystyle }{D}_t=d_{\mathrm{1}}^{W_{\mathrm{1}}}\times d_{\mathrm{2}}^{W_{\mathrm{2}}}\times d_{\mathrm{3}}^{W_{\mathrm{3}}}\times d_{\mathrm{4}}^{W_{\mathrm{4}}}\times d_{\mathrm{5}}^{W_{\mathrm{5}}}\mathrm{\cdots }\mathrm{\cdots }\times d_n^{W_n},\end{array}$$

where D_t is the instability index of the samples, expressed using relative values ranging from 1 to 10. A cumulative value closer to 10 indicates greater landslide potential, whereas a cumulative value closer to 1 indicates lower landslide potential.

By using the concept of log-normal distribution in statistics, we converted the levels of instability index derived using the MHEM into probabilities of landslide occurrence. The calculation formula of the log-normal distribution is shown in Eq. (14):

$$\begin{array}{}\text{(14)}& {\displaystyle }{\displaystyle }P\left(F\right)={\displaystyle \frac{\mathrm{1}}{x\mathit{\sigma }\sqrt{\mathrm{2}\mathit{\pi }}}}{e}^{-\frac{\mathrm{1}}{\mathrm{2}}{\left[\left(\ln x-\mathit{\mu }\right)/\mathit{\sigma }\right]}^{\mathrm{2}}},\end{array}$$

where x denotes the level of the instability index and μ and σ denote the mean and standard deviation of the level of the instability index. After calculating the probabilities of landslide occurrence by using the log-normal distribution, we normalized the probabilities to range from 0 to 1 for convenience. The normalization formula is shown in Eq. (15).

$$\begin{array}{}\text{(15)}& {\displaystyle }{\displaystyle }P{\left(F\right)}^{\prime }={\displaystyle \frac{\left(X_i-X_{min}\right)}{\left(X_{max}-X_{min}\right)}}\end{array}$$

In Eq. (15), X_i represents the factor being normalized and X_max and X_min represent the maximum and minimum values of the factor, respectively.

3 Study area

We referred to the historical data on road disasters from the NCDR (National Science and Technology Center for Disaster Reduction, 2017) and considered road sections where rainfall-induced landslides occurred frequently in southern Taiwan. We focused on the periods before and after Typhoon Nanmadol (2011) and Typhoon Kong-rey (2013) hit southern Taiwan, and we selected part of the catchment area of Laonung River in southern Taiwan as our study area (Fig. 2), which includes areas from three districts in Kaohsiung city (Jiashian, Liouguei, and Taoyan). The Laonung River flows SW across the south of the study area and originates from the Jade Mountain. The study area is located in a tropical monsoon climate zone. According to the climate statistics (1983–2012) recorded by the Central Weather Bureau, the average annual rainfall is approximately 2758 mm. Provincial Highway 20 is in an east–west direction, the starting point of the highway is Tainan city in southern Taiwan, and the ending point is Degao Community in Guanshan town, Taitung County, with a total length of 203.982 km. Within the study area, Provincial Highway 20 starts from Liouguei District (76 K + 000) in the west and goes to Taoyan Village (87 K + 500) in the Taoyan District. According to the survey data from the Directorate General of Highways (2017), the road width of Provincial Highway 20 passing through the study area is about 8.8 m. The average traffic flow and the total number of vehicles carried per day for both directions are 2260 PCU (passenger car unit) and 1434. In the study area, most of the traffic vehicles are sedans, followed by trucks and buses.

Figure 2 Study area in the south of Taiwan. Purple line depicts the mountain road.

4 Image interpretation and classification

4.1 Preprocessing of satellite images

This study employed and interpreted satellite images taken by FORMOSAT-2 (FM2). FM2 images have been extensively used to identify natural disasters and land use (e.g. Lin et al., 2004, 2006, 2011; Liu et al., 2007; Chen et al., 2009, 2013a). The FORMOSAT-2 satellite has a circular and sun synchronous orbit. With its high torque reaction wheels for all axes, the FORMOSAT-2 is able to point to a ±45^∘ along track and ±45^∘ across track and is thus able to capture any scene each day in all of Taiwan if necessary (Liu et al., 2007). FORMOSAT-2 images are available in 2 m resolution in panchromatic (pan) and 8 m in multispectral (ms) from visible to near-infrared with a coverage of 24 km×24 km. In the present study, prior to interpretation, the satellite images underwent spectral fusion, coordinate positioning, cropping, and cloud removal. The images taken by FM2 are multispectral with blue, green, red, and near-infrared (NIR) wavelengths (Chen et al., 2013a; Chue et al., 2015). Image fusion and coordinate positioning were conducted using the import data and coordinate positioning tool of ERDAS IMAGINE (2013). Then, we used the image analysis tool of ArcGIS to remove clouds from the images.

Table 2 Error matrix of interpretation results of satellite images after Typhoon Kong-rey in 2013.

Figure 3 Interpretation and classification results of satellite images before (a, c) and after (b, d) Typhoon Nanmadol (a, b) and Typhoon Kong-rey (c, d).

4.2 Training site selection and mapping

To map the sample areas required for image interpretation, we overlapped the high-resolution, preprocessed satellite images of the study area before and after the typhoons and mapped the training sites by using a GIS platform. Based on field investigations and relevant studies (Chen et al., 2009, 2013a; Chue et al., 2015), we selected areas with water, roads, buildings, crops, vegetation, river channels, and bare land within the study area as the sample area factors for interpretation training.

4.3 Image interpretation and accuracy assessment

Image interpretation and classification were conducted using the maximum likelihood module in ERDAS IMAGINE. The interpretation and classification results of the satellite images before and after Typhoon Nanmadol in 2011 and Typhoon Kong-rey in 2013 are shown in Fig. 3. The different colours in the images represent different interpretation factors.

To verify the accuracy of the results, we randomly extracted 25 points from the satellite images for each training factor as checkpoints and tested the accuracy by using the aforementioned error matrix approach. With the satellite images before and after Typhoon Kong-rey in 2013 as an example, Table 2 shows the error matrix and accuracy assessment results of the satellite image interpretation and classification processes. Table 3 presents the kappa values and OA results of the satellite images captured before and after the two typhoons. As mentioned, kappa values ranging from 0.4 to 0.8 indicate moderate accuracy, and thus the interpretation results had moderate to high accuracy.

Taking the landslide inventory after Typhoon Kong-rey in 2013 as an example, there were 291 landslides, which totaled a landslide area around 135 ha. This was equivalent to an average concentration of 5.1 % throughout the whole study area of 2659 ha. Among the 291 landslides, 256 landslides were recognized with an area smaller than 1 ha, the remaining 35 landslides were recognized in the range of 1 to 8 ha. The biggest landslide, with an area around 7 ha.

5 Landslide susceptibility assessment

To evaluate the landslide susceptibility of slopes within the study area, we constructed 8 m×8 m grids by using the GIS platform along with the interpretation results of the two typhoons. We also constructed an 8 m×8 m digital elevation model (DEM) (SWCB, 2011) and input the classification results, thematic map of predisposing factors, and rainfall data into the pixel to aid subsequent landslide susceptibility assessments.

5.1 Predisposing factor selection and factor correlation test

5.1.1 Predisposing factor selection

Referring to Chen et al. (2009), we divided the predisposing factors of landslides into three categories: natural environment, land disturbance, and rainfall.

Natural environment factors

Elevation

The influence of elevation varies with the climate and thus affects the distribution of vegetation on the slope and type of weathering. In addition, elevation reflects the influence of geological structure, stress, and time. The highest and lowest elevations in the study area were 1481 and 365 m, respectively. Using the GIS platform, we extracted the elevation data from the DEM of the study area to estimate the mean elevation of each grid. We divided the elevation data into seven classes at intervals of 300 m.

Slope gradient

A slope's gradient generally exerts significant impact on slope stability. By using the DEM and gradient analysis of the GIS platform, we calculated the mean gradient of each pixel in the study area; subsequently, we divided the gradient values in the pixels within the study area into seven classes.

Aspect

Rainfall-induced landslides are subject to the influence of seasonal changes such as those related to rainfall and wind direction. Thus, the direction of the slope must be considered. As described, we used the DEM and aspect analysis function of the GIS platform to calculate the average aspect of the pixels in the study area. According to their direction, we divided them into six classes from windward to flat ground.

Table 3 Interpretation results of satellite images before and after typhoons Nanmadol and Kong-rey.

Geology

Referring to the digital file of the Geologic Map of Taiwan, scale 1 : 50 000, Chiahsien, which was compiled by the Central Geological Survey of the Ministry of Economic Affairs in 2000, we determined that the geology of the study area includes five types of rock: the upper part of Changshan Formation, the Tangenshan Formation, the Changchihkeng Formation from the Miocene period, and modern alluvium and terrace deposits from the Holocene period. We divided geological strength into six classes (Chen et al., 2009).

Terrain roughness

Terrain roughness refers to the degree of change in pixel height. Wilson and Gallant (2000) proposed the use of the standard deviation of height within a radius to measure the degree of change in height because of its indicative meaning in relation to changes in regional height. Using the Neighborhood (focal statistics) of the Spatial Analyst Toolbox in ArcGIS, we calculated the terrain roughness of the DEM. Statistical cluster analysis was used to automatically divide terrain roughness into six classes.

Slope roughness

Slope roughness refers to the fluctuations in slope gradient in the pixels. High slope roughness means that the slope gradient varies considerably (Wilson and Gallant, 2000). Slope roughness is calculated through the same method as terrain roughness, except with the original elevation values being replaced with the slope gradient values obtained using ArcGIS. Just as terrain roughness was graded, we first used the Spatial Analyst Toolbox in ArcGIS to estimate the slope roughness of each pixel, after which we used cluster analysis to automatically divide slope roughness into six classes.

Distance to water

Streams will cause soil erosion and riparian erosion, which directly or indirectly affect the stability of the slope. We calculated the distances to water using the Buffer tool in ArcGIS and divided the distances into seven classes. Distance to road

The construction of the roads will also have an influence on the stability of the slope. Therefore, we also calculated the distances to road using the Buffer tool in ArcGIS and divided the distances into seven classes.

Land disturbance factors

Land disturbance varies with space and time. Based on the tendency to promote landslides, the index of land disturbance was developed, and we made some revisions to the qualitative approach proposed by Chen et al. (2009, 2013b) to calculate land disturbance and selected roads, buildings, crops, bare land, and vegetation as the land disturbance factors of landslides in the study area. We extracted the disaster and ground surface data from previous satellite image interpretation and classification results and input the land disturbance factors into the pixels by using the GIS platform. Referring to Chen et al. (2009, 2013b), the scores of the index for the disturbance condition (I_DC) in the pixels are assigned from five to one, corresponding to bare land, roads, buildings, crops, and vegetation.

Table 4 Effective accumulated rainfall and intensity of rolling rainfall observed at weather stations during Typhoon Nanmadol and Typhoon Kong-rey.

Rainfall factors

We collected precipitation data from weather stations of the Central Weather Bureau, including Guanshan, Biaohu, Hsiao Guanshan, Gaojhong, Sinfa, Jiashian, and Xinan. We then calculated the EAR and 3 h I_R (I_3R) levels observed at each station. The results from Typhoon Nanmadol in 2011 and Typhoon Kong-rey in 2013 are compiled in Table 4. By using the inverse distance weighting (IDW) function of ArcGIS and the EAR and maximum I_3R values of the weather stations, we estimated the rainfall of each pixel throughout the study area and then used cluster analysis to divide the results into six classes.

Table 5 Correlation test results between the predisposing factors.

Table 6 Paired sample t test results between elevation and terrain roughness and slope gradient and slope roughness.

5.1.2 Factor correlation test

To establish a landslide susceptibility assessment model, we selected elevation, slope gradient, aspect, geology, terrain roughness, slope roughness, distance to water, distance to road, I_DC, and rainfall as landslide-predisposing factors. Rainfall included EAR and maximum I_3R.

We employed the Pearson correlation test tool in SPSS software (2005) to examine the correlation among these factors. The correlation coefficients ranged from −1 to +1, with +1, −1, and 0 indicating complete positive correlation, complete negative correlation, and no correlation between two variables, respectively. Factors with high correlation were then subjected to a paired sample t test conducted using SPSS to examine the significance of the correlation between them. Those with high correlation were eliminated.

Table 5 presents the test results regarding the correlation between the predisposing factors. As shown, the degree of correlation between most factors was moderate to low. A high degree of correlation was found only between elevation and terrain roughness and between slope gradient and slope roughness. Thus, we administered paired sample t tests to these two factor pairs to test the significance of the correlation. The results in Table 6 show that the significance was 0 (<0.05) for the correlation between both pairs, indicating no correlation; thus these factors were not eliminated.

5.2 Landslide susceptibility assessment and hazard map

To apply the MHEM in order to establish a landslide susceptibility assessment model, we input the natural environment, land disturbance, and rainfall factors into the pixels by using the GIS platform. By using the changes in bare land between the images before and after the typhoons and applying image subtraction aided by manual checking, we obtained the pixel data of the rainfall-induced landslide locations in the study area. With the study area after Typhoon Nanmadol in 2011 as an example, we examined EAR during the rainfall period and rated the classes by using the factor weights derived using the MHEM, as shown in Table A1 of Appendix A.

The calculation process is explained in this paper using elevation as an example. In accordance with factor selection, the elevation factor was divided into seven classes. Aided by the GIS platform, we calculated the total number of pixels, total number of landslides, and landslide percentage within each elevation level in the study area by using Eqs. (8) and (9). Based on the landslide percentages of the elevation factor and the minimum and maximum landslide percentages, we subsequently obtained the scores of the factors by using Eq. (10). We then calculated the standard deviation, coefficient of variation, and weight values by using Eqs. (11) and (12); the results are listed in Table A1 of Appendix A. The presented results show that the standard deviation (σ), coefficient of variation (V), and factor weight (W) of the landslide percentage were 0.021, 0.764, and 0.087, respectively. Finally, we calculated the instability indices by using the weight values and scores of the factors through Eq. (13). Furthermore, the results in Table A1 of Appendix A indicate that the degrees of land disturbance (I_DC), geology (G_s), slope gradient (S_s), and slope roughness had the greatest influence on landslides in the study area, followed by distance to water (D_s), EAR, and elevation (E_l).

Table 7 Intervals of instability index and landslide probability of rainfall factors.

For EAR and I_3R we used an instability index to determine the level of landslide susceptibility of slopes throughout the study area. The derived instability index intervals (Table 7) for EAR and I_3R ranged from 2.05 to 9.59 and 2.02 to 9.96. By using Eqs. (14) and (15), the landslide probability intervals calculated based on EAR and I_3R are presented in Table 7.

We employed the mean probability of landslide occurrence to differentiate between high and low landslide susceptibility. Landslides were considered more likely to occur in areas where the probability of landslide occurrence was greater than the mean. By contrast, landslides were considered less likely to occur in areas where the probability of landslide occurrence was lower than the mean. With rainfall factor EAR as an example, we determined the mean probability of landslide occurrence to be 0.46. We further divided landslide susceptibility into four levels: high (0.731–1), medium high (0.461–0.73), medium low (0.23–0.46), and low (0–0.23). The results showed that the mean probability of landslide occurrence varied little, regardless of whether it was calculated using EAR or I_3R.

Table 8 Accuracy of landslide susceptibility map for different rainfall factors and typhoons.

Table 9 Numbers of landslide pixels in the study correspond to different D_t levels under different rainfall factors after typhoons.

By using the GIS platform, we considered the landslide susceptibility calculated using EAR for Typhoon Nanmadol in 2011 as an example. As illustrated in Fig. 4, we included an overlay created by the NCDR showing the locations of historical disasters within the study area. The results revealed a total of 24 historical disasters, 17 of which were situated in areas of medium high or high landslide susceptibility. Therefore, the estimation accuracy in this study was approximately 71 %. Regarding Typhoon Kong-rey in 2013, 18 historical disasters occurred within areas of medium high or high landslide susceptibility, thereby yielding 75 % accuracy. Table 8 presents the accuracy levels associated with using different rainfall factors to calculate landslide susceptibility for different typhoons.

Figure 4 Landslide susceptibility in the study area, in which cross symbols represent the historical disasters collected from NCDR (2017).

Figure 5 Relationships among instability index, effective accumulated rainfall, and landslide occurrence in the study area after typhoons Nanmadol (2011) and Kong-rey (2013).

5.3 Investigation of rainfall factors and instability index

To understand the relationship between the rainfall factors and the degree of instability on the slopes in the study area after the typhoons, we first removed the cloud cover grids from post typhoon images and subsequently employed cluster analysis to divide the instability indices of the pixels into three levels: high, medium, and low. We then collected random samples based on the proportions of landslide and non-slide pixels at each level (50 landslide and 50 non-landslide pixel points) and plotted their relationship. Table 9 and Fig. 5a–d present the relationships between the rainfall factors (EAR and I_3R), instability index, and landslide occurrence in the pixels following Typhoon Nanmadol in 2011 and Typhoon Kong-rey in 2013. Figure 5a and b show EAR, whereas Fig. 5c and d show I_3R. The presented results indicate that the typhoon events increased the degree of slope instability (D_t) and landslide occurrence, regardless of whether EAR or I_3R was considered. Furthermore, significantly more landslide points were situated in areas of high instability than in areas of other levels of instability, and landslides rarely occurred in areas of low instability. Moreover, areas of high slope instability were prone to landslides even if their EAR or I_3R was low. By contrast, areas of low instability required more rainfall for the occurrence of landslides. The results (Table 9) further showed that the EAR and I_3R levels of Typhoon Kong-rey in 2013 were greater than those of Typhoon Nanmadol in 2011. Thus, at any D_t level, the proportion of landslides that occurred in the study area after Typhoon Kong-rey was higher than that after Typhoon Nanmadol. Figure 5e and f present the relationships between EAR×I_3R, the instability index, and landslide occurrence; EAR×I_3R is the index of rainfall-induced landslide (ILR), with a higher value indicating higher susceptibility to a landslide. The figures show that for a high instability index, even a small rainfall event could trigger a landslide (lower-right corners of the figures). By contrast, for a low instability index, a larger rainfall event could not easily trigger a landslide (upper-left corners of the figures).

Figure 6 Relationships between landslide predisposing factors and number of landslide pixels in the study area.

Figure 7 Diagrams of landslide area on slope, in which d_r represents the distance between the highest point of a landslide area and the nearest ridge, d_s the distance between the lowest point of the landslide area and the nearest stream, and d_t the distance between ridge and stream.

6 Landslide location analysis

We analysed the spatial characteristics of landslides by using landslide locations collected from before and after the two typhoons and the land surface interpretation results of the study area.

6.1 Investigation of landslide predisposing factors and landslide area

The influence of predisposing factors on landslides varies. In this study, we examined the relationships between landslide area and various predisposing factors. By using the area of landslides (i.e. the number of landslide pixels) induced by Typhoon Nanmadol in 2011 as an example, we investigated the influences of the predisposing factors (elevation, slope gradient, aspect, geology, slope roughness, terrain roughness, distance to water, distance to road, and degree of land disturbance) on landslides. The various factor classes and corresponding numbers of landslide pixels are shown in Fig. 6a–i.

Figure 6a presents the relationship between different classes of elevation and the number of landslide pixels (landslide area). As shown in the figure, the number of landslide pixels in the study area peaked at elevations between 450 and 750 m and then declined as the elevation increased. Figure 6b displays the relationship between different classes of slope gradient and the number of landslide pixels (landslide area). As shown in the figure, the number of landslide pixels in the study area increased with the slope gradient and peaked between 30 and 55^∘. Landslides rarely occurred on slopes steeper than 55^∘. Figure 6c illustrates the relationship between aspect and the number of landslide pixels, with aspect divided into eight categories: north, north-east, east, south-east, south, south-west, west, and north-west. As shown in the figure, the number of landslide pixels was highest on slopes facing south, followed by those on slopes facing east and south-east. We speculate that this is because rainfall during the typhoon season in Taiwan promotes poor cementation and high weathering on slopes along rivers, which consequently prompts these slopes to develop toward low-lying rivers (which run from the north-east to the south-west) after rainfall events.

Figure 6d shows the relationship between geology and the number of landslide pixels. As shown in the figure, the Sanhsia Group and its stratigraphic equivalence lead to landslides more easily than the Lushan Formation in the study area. The Sanhsia Group and its stratigraphic equivalence mainly comprise sandstone, shale, and interbedded sandstone and shale. Shale has weaker cementation, lower strength, and a greater tendency to weather and fracture. By contrast, the Lushan Formation consists of argillite, slate, and interbedded argillite and sandstone, and its strength is controlled by cleaving; some areas are prone to weathering and fracturing. Thus, both rock types are more likely to collapse, but on the whole, the Sanhsia Group and its stratigraphic equivalence collapse more easily than the Lushan Formation. Furthermore, this result indicates that the locations of landslide areas within the study area are associated with geology. Figure 6e presents the relationship between slope roughness and the number of landslide pixels. The number of landslide pixels within a level of slope roughness first increased with the slope roughness and then began to decline once a certain level of slope roughness (35–40) was reached. This result is similar to that of the influence of the slope gradient on the number of landslide pixels. Figure 6f displays the relationship between terrain roughness and the number of landslide pixels. As shown in this figure, the results are similar to those regarding the influence of elevation on the number of landslide pixels: the number of pixels declined when the terrain roughness was greater than 500 and was very small when the terrain roughness was greater than 1200.

Figure 6g illustrates the relationship between distance to water and the number of landslide pixels. The presented results show a significantly greater number of landslide pixels within 300 m of water. The width of the river channel within the study area was determined to range from 100 to 200 m, revealing that the development of landslide areas near water in the study area is caused by rainfall significantly raising the water level in the river, which scours the slope toe, affects slope stability, and triggers landslides. Figure 6h presents the relationship between distance to road and the number of landslide pixels. The presented results reveal that areas between 100 and 300 m from roads had the greatest number of landslide pixels. Further examination of the relationship between distance to road and the area and number of landslides revealed that most landslides between 0 and 100 m from roads were small collapses, whereas those between 100 and 300 m from roads were larger in area. The number of landslides 0–100 m from roads was greater than that 100–300 m from roads.

Figure 8 Spatial distribution of bare land in the study area before and after the typhoons Nanmadol (top) and Kong-rey (bottom). The sizes of the bubbles reflect the areas of bare land.

The degree of land disturbance can represent the changes in surface conditions including roads, buildings, crops, bare land, and vegetation. A greater degree of land disturbance likely indicates a greater degree of surface changes, which can yield a greater number of landslide pixels. Figure 6i shows the relationship between the degree of land disturbance and the number of landslide pixels. The presented results indicate that the number of landslide pixels increased with the degree of land disturbance.

6.2 Landslide scale and spatial distribution

We employed the terrain tool in ERDAS IMAGINE and the DEM to identify the ridges and valleys in the study area. Following the methods in previous studies (Meunier et al., 2008; Chue et al., 2015), we extracted the distances between the highest point of a landslide area and the nearest ridge (d_r), between the lowest point of the landslide area and the nearest stream (d_s), and between the ridge and the stream (d_t) (Fig. 7). Furthermore, in Taiwan, many slopes are visible on developed mountain roads built between ridges and streams. Therefore, we explored the spatial distribution of landslides above and below mountain roads. Similarly to Fig. 7a, to explore the spatial distribution of landslides, we extracted the distances between the highest point of a landslide area on a slope above a road and the nearest ridge (d_r), between the lowest point of the landslide area and the nearest mountain road (d_mu), and between the ridge and the mountain road (d_tu) (Fig. 7b). We also investigated this distribution by extracting the distances between the highest point of a landslide area on a slope below a road and the nearest mountain road (d_md) between the lowest point of the landslide area and the nearest stream (d_s) and between the mountain road and the stream (d_td) (Fig. 7c).

This study examined the spatial distribution of landslides in the region along Provincial Highway 20 before and after Typhoon Nanmadol in 2011 and Typhoon Kong-rey in 2013. Using the approach shown in Fig. 7a, we mapped the bare land in the study area, as shown in Fig. 8a–d. Of these figures, Fig. 8a and c show the conditions before the typhoons, whereas Fig. 8b and d present the conditions after the typhoons. The presence of bare locations near the Y axis ($d_{\text{r}}/d_{\text{t}}\approx \mathrm{0}$) denotes that the bare land originated near the ridge. By contrast, the presence of bare locations near the X axis ($d_{\text{s}}/d_{\text{t}}\approx \mathrm{0}$) denotes that the bare land progressed toward the stream. Thus, the presence of bare locations near the origin denotes that the bare land originated near the ridge and progressed toward the stream.

The results in Fig. 8a–d show more bare locations in the lower-right halves of the graphs, some of which are larger in area. The figures indicate fewer bare locations in the upper-left halves of the graphs, and the ones that are present are smaller in area. These spatial distribution characteristics are similar to those derived by Meunier et al. (2008). We speculate that this is because the frequency of rainfall-induced landslides increases significantly because of bank erosion, which is shown in the lower-right half of Fig. 8 ($d_{\text{r}}/d_{\text{t}}\ge \mathrm{0.5}$ and $d_{\text{s}}/d_{\text{t}}\le \mathrm{0.5}$). Furthermore, the bare locations before and after typhoons Nanmadol and Kong-rey show that the bare land does not increase in number but increases significantly in area, implying that old landslides may result in more collapses or expansions of the affected area. In addition, the number of old landslides is greater than that of new landslides.

Figure 9 Spatial distribution of bare land on slopes above mountain roads in the study area before and after Typhoon Kong-rey in 2013. The sizes of the bubbles reflect the areas of bare land.

We explored the spatial distribution of landslides on slopes above (Fig. 9) and below (Fig. 10) mountain roads in the study area before and after Typhoon Kong-rey in 2013. Figures 9a and 10a present the spatial distribution of bare land before the typhoon, whereas Figs. 9b and 10b present the spatial distribution of bare land after the typhoon.

Figure 10 Spatial distribution of bare land on slopes below mountain roads in the study area before and after Typhoon Kong-rey in 2013. The sizes of the bubbles reflect the areas of bare land.

As shown in Fig. 9, most landslides on the slopes above the mountain roads occurred close to the roads, most likely because road construction involves cutting the slope toe and increasing the gradient. After the typhoon, the bare locations on the slopes above the roads in the study area did not increase in number significantly; thus rainfall did not exert a substantial impact on the slopes above the roads. The results in Fig. 10 show bare locations on the slopes below the mountain roads developing from near the roads to the streams. The bare locations near the streams may also have been affected by rainfall-induced bank erosion. However, the bare land near the roads may have been a result of roads being constructed in the study area, which affects slope stability and increases the probability of landslides. Furthermore, the bare locations near the roads slightly increased in number after the typhoon, likely because the roads changed the routes of surface run-off. The area of bare land near the streams also increased, possibly because the water flow scours the slope toe and causes continual bank collapses. Thus, typhoons have a significant impact on the stability of slopes below mountain roads.

7 Conclusions

This study applied the maximum likelihood method to interpret and classify satellite images before and after two typhoons in 2011 and 2013. We extracted landslide and land use information from the areas surrounding roads and then compiled the rainfall and DEM data from the typhoon events. By using the MHEM, we established a landslide susceptibility assessment model and examined the relationships between predisposing factors and the area and number of landslides within the study area, as well as the relationships between roads and the spatial distribution of landslides. The results show that the kappa coefficients associated with the use of the maximum likelihood method to interpret and classify satellite images before and after Typhoon Nanmadol in 2011 and Typhoon Kong-rey in 2013 ranged from 0.53 to 0.66, whereas the OA ranged from 61 to 71 %, indicating moderately high accuracy. According to the results of the instability index-based landslide susceptibility assessment model, the degree of land disturbance, geology, slope gradient, and slope roughness had the greatest impacts on landslides. A comparison of historical landslides triggered by the typhoons and the results of the hazard map revealed 71 % accuracy for Typhoon Nanmadol in 2011 and 75 % accuracy for Typhoon Kong-rey in 2013. Regarding the influence of the predisposing factors, an elevation of 450–750 m, a slope gradient of 30–55^∘, and distances within 300 m of water or roads were associated with a larger scale of landslides. The scale of landslides also increased with the degree of land disturbance. The relationships between the ILR, instability index, and landslide occurrence indicate that for a high instability index, even a smaller rainfall event could trigger a landslide. By contrast, for a low instability index, a larger rainfall event could not easily trigger a landslide. Thus, the instability index can effectively reflect landslide susceptibility. Comparisons of the distribution of bare land before and after typhoon events showed that most landslides in the study area were caused by stream water scouring away the toes of bank slopes. Although bare locations did not significantly increase in number after the typhoon events, they increased significantly in area, implying that the number of old landslide areas holding more collapses or expansions was greater than that of new landslide areas developing. In addition, the results obtained from observing changes on slopes above and below mountain roads after the typhoon events indicate that the number of bare locations on the slopes above the roads in the study area did not increase significantly, whereas the bare locations near the roads on the slopes below the roads slightly increased in number after the typhoon events, likely because of the roads changing the routes of surface run-off. The amount of bare land near streams also increased, possibly because the water flow scours the slope toe.

Appendix A

Table A1 Weights and scores of predisposing factors after rainfall brought by Typhoon Nanmadol in 2011.