Assessment of Rainfall-Induced Landslide Distribution Based on Land Disturbance in Southern Taiwan

Abstract

This study explores the impact of rainfall on the followed-up landslides after a severe typhoon and the relationship between various rainfall events and the occurrence, scale, and regional characteristics of the landslides, including second landslides. Moreover, the influence of land disturbance was evaluated. The genetic adaptive neural network was used in combination with the texture analysis of the geographic information system for satellite image classification and interpretation to analyze land-use change and retrieve disaster records and surface information after five rainfall events from Typhoon Morakot (2009) to Typhoon Nanmadol (2011). The results revealed that except for extreme Morakot rains, the greater the degree of slope disturbance after rain, the larger the exposed slope. Extreme rainfall similar to Morakot strikes may have a greater impact on the bare land area than on slope disturbance. Moreover, the relationship between the bare land area and the index of land disturbance condition (I_LDC) is positive, and the ratio of the bare land area to the quantity of bare land after each rainfall increases with the I_LDC. With higher effective accumulative rainfall on the slope in the study area or greater slope disturbance, the landslide area at the second landslide point tended to increase.

1. Introduction

Taiwan is prone to typhoons and heavy rainfall. Because of steep terrain and unfavorable geological conditions in mountainous areas of Taiwan, concentrated rainfall during typhoons or rainstorm attacks leads to landslide and debris flow [1]. In the last few decades, large-scale landslides, river siltation, and bank dike collapse caused by extreme rainfall have led to drastic changes in the natural environmental conditions of the watershed. For example, the Morakot Typhoon in August 2009 hit Taiwan with heavy rain, causing severe landslides, debris flow, and flooding in the southern, central, and eastern regions [2,3]. Such heavy rainfall that accompanies typhoons has caused frequent disasters in various hillside areas and affected people in those regions. Such disasters greatly affect the safety of people’s lives, property, and the living environment in the region, also presenting a threat to the construction of major public projects. The Taiwan government has invested large amounts of reconstruction funds in disaster-stricken areas. However, even after the reconstruction, these regions face the risk of repeated sediment disasters on the slopes. As a result, the overall economic development and transportation are threatened under the influence of disasters. Many settlements have been repeatedly affected by sediment disasters. Therefore, identifying strategies for the prevention and treatment of landslide disasters and repeated sediment disasters is necessary.

Commonly used satellite image classification methods can be categorized as supervised and unsupervised [4]. Unsupervised classification is based on the number of classifications (i.e., the predetermined number of clusters) set by the image analyst by using cluster analysis to achieve the predetermined number of classifications. By contrast, in supervised classification, the user must first specify the training sample area and then use the spectral grayscale value distribution of the training sample area as the basis for classification. Satellite imagery presents advantages, such as short surface-change–acquisition time, wide coverage, and low cost and is thus feasible for use in mountainous or remote areas. In addition, the acquired data can be analyzed on a computer and through the geographic information system (GIS), which can help quickly analyze the ground cover situation. Therefore, satellite imagery is a suitable tool for monitoring large-scale and time-series land-use changes [5,6].

Satellite images have recently been used for interpreting large-scale landslide damage [7,8,9,10,11,12,13,14]. In addition, some scholars have published relevant research on the application of artificial intelligence (AI) in the classification and interpretation of satellite images [10,14,15,16,17,18]. Texture, an aspect of an image, refers to changes in gray levels of the adjacent pixels in the image, correlation of the color space, or visual performance of the gray and color changes of the image in terms of spatial positions, such as edges, shapes, stripes, and color blocks. Therefore, texture analysis can help distinguish images in different groups [19]. Some scholars have used texture information as auxiliary information for satellite image interpretation to improve interpretation accuracy [14,20,21,22].

In general, landslides in each area are influenced by several factors, including latent factors (such as geology, soil, topography, hydrology, and land use) and incentives (such as rain and earthquake) [10,14,23]. Many relevant studies have investigated and discussed the hazards of landslides [10,14,24,25,26,27]. Human activities have affected slope stability through civil construction and slope reclamation. In some specific areas, human activities are the main influencing factor in landslides [28]. The sensitivity of different types of land-use has been evaluated by simulation analysis of spatial and temporal changes and suggested that land-use should be carefully planned to reduce the likelihood of disasters [29]. In addition, some scholars have studied the spatial distribution of landslides caused by rainfall and the impact of various influencing factors on the location of landslides [14,30,31,32]. The spatial distribution of landslides can reflect the occurrence of landslides [33,34,35]. According to the spatial distribution of landslides in different periods, Samia et al. (2017) [36] quantified the effects of earlier landslides on subsequent landslides.

The aim of this research is to investigate the effect of rainfall on the landslide after a catastrophic typhoon, the relationship between different rainfall scales and spatial distribution of the landslide and second landslide, as well as the influence of land disturbance. The artificial neural network (ANN) was used in combination with image texture analysis in the satellite image-based classification and interpretation of sediment disasters and land-use changes. Genetic algorithms, which automatically evolve, and train optimized neural network architectures to obtain disaster records and surface data, were used. In the past few years, the artificial neural networks (ANN) of artificial intelligence (AI) have been proven the advantages in satellite image interpretation [37,38]. The ANN extracts more parametric information through learning and interaction capabilities and has the ability to deal with linear and nonlinear relationships simultaneously [39]. In addition, ANN learning has considerable tolerance for error yet requires relatively low computation and memory demands [40]. The ANN method performs better than some traditional classifiers, such as maximum-likelihood classification [41,42].

The study also explored the potential hazards of exposure of the studied watershed area to sediment hazards by conducting a correlation analysis of the landslide history, land disturbance, and regional environmental characteristics. The mechanism of the second landslide in the site of the original landslide, the scale of the landslide, and mutual characteristics of the location was explored to obtain a reference for the prevention of future landslides and establishment of countermeasures.

2. Study Areas

In August 2009, Typhoon Morakot struck Taiwan, resulting in 704 deaths and 22 disappearances. The mountains in Southern Taiwan experienced the largest accumulative rain of more than 2900 mm, which caused the highest damage in Southern Taiwan. Landslide disasters caused by typhoons or heavy rain have been increasing. Therefore, after referring to the historical data on road disasters from the National Science and Technology Center for Disaster Reduction [43], Baolai Village and Jianshan Village in Liouguei District and Taoyuan District, respectively, in Lao-nong River Basin, Southern Taiwan, were selected as the research areas (Figure 1). By using satellite image records in the typhoon or rainstorm-affected areas in the study region, information regarding the 2009 Morakot typhoon attack up to 2011, five typhoon events during the 3-year period (2009 Morakot, 2010 Meranti, 2010 Fanapi, 2011 Meari, and 2011 Nanmadol), and a 24-h rainfall event with an accumulative rainfall of more than 130 mm (27 July 2010) were retrieved.

3. Methodology

3.1. Genetic Adaptive Neural Network

An ANN, an AI technology, presents a high-speed computing ability, high memory, and high learning ability. The ANN presents the following advantages: no requirement of any assumption of the type of equation in advance, simple training method, ability to process large amounts of data, high-precision, and high-speed. These advantages make the ANN suitable for processing highly complex nonlinear function problems. In this study, the backpropagation network (BPN) with gradient descent was used for training the neural network. BPN is one of the most representative and commonly used models in the current neural network learning model. The details of this algorithm were reported by Hagan et al. [44]. Based on the studies of Chen et al. [18,45], this study adopted tansig as the activation function, whereby the normalized output value was set between −1 and 1.

The application of ANN often presents challenges such as yielding the local minima rather than global minima, insufficient training or overtraining, and inability to converge. To improve the BPN performance, based on Adeli and Hung’s study [46], this study used genetic algorithms (GAs) to obtain optimal values of the parameters used in the BPN framework. The basic principle is to set the weighting matrix and other network parameters as the chromosomes in GAs, which is also the target solution to be searched. This study employed the genetic adaptive neural network (GANN) proposed by Chen et al. [45] and Chen et al. [18] and adopted the difference between the network predicted value and the actual value as the adaptive function of the algorithm; that is, the mean squared error was used as the evaluation guideline for chromosome adaptability. The actual value represents the classification of the factor to be interpreted, and the network prediction value is the spectral value of each sample area input. After simulation with GANN, the obtained classification prediction value is considered the output. In this study, GAs were employed to search for the optimal solution through a set of network parameters that minimize the error between the predicted and actual values.

For the GANN, various types of coding, crossover, and mutation have been proposed for different application areas. In this study, binary coding was adopted as the coding method for GAs, and elitism selection and uniform mating were the algorithms used. Based on the study of D’Ambrosio et al. [47], this study also selected 200 as the number of generations (number of chromosome groups) in the calculation process. The commonly used selection rate is approximately 0.10, and the general mating probability is mostly set between 0.5 and 0.8. Heng et al. [48] also proposed that the mutation rate should be between 0.1 and 0.001. Based on the study of Chen et al. [10], this study set the mating rate and mutation rate as 0.6 and 0.1, respectively. This study also referred to the studies of Chen et al. [18,45] and accordingly set the upper limit of hidden layers to 2, the upper limit of the number of neurons to 32, the upper limit of the learning rate to 3.2, and the upper limit of the number of learning times to 15,000. The evolution process is described in detail in Chen et al. [10,18,45].

3.2. Texture Analysis

This study used image texture analysis to establish and quantify the gray-level co-occurrence matrix (GLCM) of satellite images obtained from the study area [49]. GLCM is a second-order statistical method and was used to determine the frequency of each pair of grayscale values appearing at specific relevant locations and for the calculation of texture feature values. The calculation method is shown in Equation (1):

$$C_{ij}\left(d,\theta \right)=\frac{P_{ij}\left(d,\theta \right)}{{{\displaystyle \sum }}_{i=0}^N{{\displaystyle \sum }}_{j=0}^N{P}_{ij}\left(d,\theta \right)}$$

(1)

where P_ij represents the joint probability that grayscale values i and j appear in relative positions (d, θ) in the image, d is the distance, and θ denotes the direction.

Haralick et al. [49] suggested the use of a variety of texture statistics to quantify GLCM. Based on the studies of Chen et al. [50] and Chue et al. [14], this study used the following texture quantization formulas and used the resultant value as the input for the GANN during image interpretation training and classification:

$$Homogeneity={\displaystyle \sum }_{i=0}^N{\displaystyle \sum }_{j=0}^N\frac{1}{1+{\left(i-j\right)}^2}\cdot C_{ij}\left(d,\theta \right)$$

(2)

$$Contrast={\displaystyle \sum }_{i=0}^N{\displaystyle \sum }_{j=0}^N{\left(i-j\right)}^2\cdot C_{ij}\left(d,\theta \right)$$

(3)

$$Dissimilarity={\displaystyle \sum }_{i=0}^N{\displaystyle \sum }_{j=0}^N\left|i-j\right|\cdot C_{ij}\left(d,\theta \right)$$

(4)

$$Entropy={\displaystyle \sum }_{i=0}^N{\displaystyle \sum }_{j=0}^N{C}_{ij}\left(d,\theta \right)\cdot log C_{ij}\left(d,\theta \right)$$

(5)

$$\mathrm{Angular} \mathrm{second} \mathrm{moment}={\displaystyle \sum }_{i=0}^N{\displaystyle \sum }_{j=0}^N{C}_{ij}^2\left(d, \theta \right)$$

(6)

3.3. Accuracy Assessment

After satellite image interpretation, the image classification is considered suitable for monitoring if the classification results fulfill the requirements. After referring to the method of Chen et al. [10], this study used the error matrix [51] to evaluate the accuracy of the image interpretation classification. Take the error matrix in

Table 1. Example of error matrix.

		Actual Surface		Total
		Category A	Category B
Classification results	Category A	E11	E12	E1+
	Category B	E21	E22	E2+
Total		E+1	E+2	E++

as an example. Table presents four categories. The column represents the reference data category, and the row represents the classification category. All results can be derived from the matrix. The error matrix compares the reference category of the sampling area and the satellite image in the same sampling space to explain the difference between the classification results.

E₁₂ in the table indicates the amount of data that actually belongs to category B but is misclassified under category A in the interpretation result; E₂₁ indicates the amount of data that actually belongs to category A, but is misclassified under category B in the interpretation result; E₁₁ represents the number of correctly classified data of type A, and E₂₂ represents the number of correctly classified data of type B.

From the values in the error matrix, several primary and commonly used classification accuracy indicators [51] can be calculated; of them, overall accuracy (OA), the simplest indicator of the accuracy of a method, is calculated as follows:

$$OA=\left(\frac{1}{N}{\displaystyle \sum }_{i=1}^r{E}_{ii}\right)\times 100\%$$

(7)

where N and r represent the total number of classifications and the number of rows in the matrix, respectively.

In addition, the classification accuracy for a single category can be expressed by user accuracy (UA) and producer accuracy (PA). UA indicates the probability that the ground cover is classified correctly and is calculated as follows:

$$UA=\frac{E_{ii}}{E_{i+}}\times 100\%$$

(8)

where E_ii is the number of pixels actually present in a given class and E_i+ is the number of pixels classified to that class.

PA refers to the accuracy with which the ground reference materials can be correctly classified using a certain classification method and is calculated as follows:

$$PA=\frac{E_{ii}}{E_{+i}}\times 100\%$$

(9)

In addition, C [52] proposed a kappa ($\hat{K}$) index (Equation (10)), which indicates the extent to which the classification results are more accurate than those obtained from random classification. The kappa statistical index considers the difference between two consistency types: the consistency between automatic classification and reference data and the probability of consistency between sampling and reference classification. In general, the consistency coefficient (kappa index) is between 0 and 1. The larger the kappa value, the higher the classification accuracy is. According to Landis and Koch [53], the classification accuracy is high if the kappa value is greater than 0.8, moderate if the kappa value is between 0.4 and 0.8, and low if the kappa value is less than 0.4.

$$\hat{K}=\left(\frac{N{{\displaystyle \sum }}_{i=1}^r{E}_{ii}-{{\displaystyle \sum }}_{i=1}^r\left(E_{i+}\times E_{+i}\right)}{N^2-{{\displaystyle \sum }}_{i=1}^r\left(E_{i+}\times E_{+i}\right)}\right)\times 100\%$$

(10)

3.4. Accumulative Rainfall Analysis

The rainfall was calculated by the effective accumulative rainfall (EAR), which is defined as the sum of the direct rainfall from continuous rainfall and the rainfall received in the previous period [54]. Seo and Funasaki [54] proposed that concentrated rainfall is considered continuous (main rain field) when there is no rainfall 24 h before and after it. Direct rainfall is the accumulative rainfall from the first rainfall for the main rain field from the time of the landslide disaster (a time point when the first rainfall reaches ≥4 mm) to the time of the landslide disaster (the point of maximum rainfall in the main rain field). Previous indirect rainfall (P_b) was calculated as the amount of rainfall for the main rain field for 7 days as

$$P_b={\displaystyle \sum }_{n=1}^7{k}^n{P}_n$$

(11)

where P_n is the rainfall (mm) n days before the main rain field; k is the decreasing coefficient. Based on the studies of Chen et al. [10] and Chue et al. [14], this study set k as 0.9. Because direct rainfall (P_r) directly affects the occurrence of landslides, no reduction is allowed. Therefore, the EAR can be calculated as follows:

$$EAR=P_r+P_b$$

(12)

3.5. GIS

In the last decade, the GIS emerged as an essential tool for land-use planning and management and many other applications. GIS is a set of integrated systems that combine geographical information and computer technology. A GIS is composed of two parts: database and functional systems. The data are divided into two types, attribute and spatial data, mainly by analyzing information and establishing decision-making. Original geographical data can be transformed into spatial decision-making information by using forecasting models. In this study, thematic maps relevant to land-use and landslide occurrence were used to construct a vector-type spatial database using the ESRI ArcGIS [55]. The Data Management Tools and Conversion Tools in ArcToolbox of ArcGIS were used for the establishment of the basic grid, satellite image preprocessing, and interpretation, respectively. Moreover, the 3D Analyst Tools of ArcGIS were used to extract the topographic location of the bare ground.

4. Results

4.1. Construction of Original Cartographic Information

Some of the original cartographic information in this study was based on data from Chue et al. [14], including the satellite image maps, digital elevation maps (DEMs), and geological maps of Formosat-2 (FM2) before and after the six storms (typhoons) during 2009–2011. The cartographic information was combined with GIS to establish relevant attributes and spatial databases. The satellite images were first screened based on the date of the typhoon (rainstorm) event recorded by the Central Weather Bureau. From the screened images, those most suitable for use (low cloud cover rate) were selected, and finally, 10 satellite images were selected for interpretation. The data are obtained from [14] and shown in

Table 2. Basic data of satellite imagery

Event Number	Before/after Event (Date)	Image Shooting Date	Image Resolution
I	Before Typhoon Morakot (2009-08-05)	2009-05-09	2 m
	After Typhoon Morakot (2009-08-05)	2009-08-24
II	Before rainfall (2010-07-27)	2010-05-25
	After rainfall (2010-07-27)	2010-08-10
III	Before Typhoon Meranti (2010-09-09)	2010-08-10
	After Typhoon Meranti (2010-09-09)	2010-09-11
IV	Before Typhoon Fanapi (2010-09-17)	2010-09-11	8 m
	After Typhoon Fanapi (2010-09-17)	2010-11-21
V	Before Typhoon Meari (2011-06-23)	2011-05-08
	After Typhoon Meari (2011-06-23)	2011-08-17
VI	Before Typhoon Nanmadol (2011-08-27)	2011-08-17
	After Typhoon Nanmadol (2011-08-27)	2011-10-24

.

4.2. Establishment of Basic Grid

The DEM used in this research has a resolution of 40 × 40 m². Therefore, regardless of the FM2 image or DEM, the GIS software program ArcGIS was used to uniformly create a basic 40 × 40 m² grid, and the ArcGIS Spatial Analyst function was used to calculate the average slope, average elevation of each grid, and geological type of the grid.

4.3. Interpretation and Classification of Satellite Images

4.3.1. Satellite Image Preprocessing

In this study, the original FM2 satellite image data for each spectral band were obtained. We first used the ERDAS Imagine [56] function to merge the various spectral band images and then entered the coordinate position, pixel resolution size, and measurement units of the satellite images to correct the satellite image coordinates. In addition, in the interpretation of satellite images, to avoid the influence of cloud coverage on the interpretation results, this study used ArcGIS [55] in conjunction with the manual selection drawing method and used ERDAS Imagine removing the cloud part.

In this study, the aforementioned GANN was used to interpret and classify satellite images of the study area. The original spectral values of the satellite image in the study area and their spectral texture information was used as the input for network training. In this study, ArcGIS was used to capture the original spectral values of red light R (red), green light G (green), blue light B (blue), and near-infrared (NIR) light in satellite images. ENVI [57] was used to calculate the homogeneity, contrast, dissimilarity, entropy, and angular second moment of the texture information.

Furthermore, before drawing the training plots, this study first standardized the input values of network training. The normalization formula is as follows:

$$S=\frac{x-mean\left(X\right)}{std\left(X\right)}$$

(13)

where S is the normalized value of x, mean(X) is the average, and std(X) is the standard deviation of factor X.

4.3.2. Satellite Image Interpretation and Classification

Based on the studies of Chen et al. [10,13], this study used nine factors, namely waters, fruit trees, buildings, forests, grasslands, bare land, farmland, roads, and rivers, as interpretive classification factors. Before satellite image interpretation, the training sample area for the features required in this study should be selected. For selecting the training sample area, the study used aerial photographs as the base map together with the current data to select the various interpretation classification factors used in this research. For each factor, a training plot with approximately 5% of image grids (approximately 1000 image grids) in the study area must be selected, and the selected plots were required to be evenly distributed within the study area.

In this study, the relevant GANN parameter settings were modified after referring to the recommendations of Chen et al. [10,45]. The upper limits of the number of hidden layers, number of neurons, learning rate, and number of learning times were set as 2, 32, 3.2, and 15,000, respectively. After GANN training of the sample area, the optimal structure of the neural network was obtained through testing. This optimal structure was used to classify the full-scale satellite images. The optimal framework parameters of the 10 satellite images obtained after training are shown in

Table 3. Summary of genetic adaptive neural network (GANN) framework parameters and training results of satellite images.

Date	Hidden Layers	Neurons for 1st Hidden Layer	Neurons for 2nd Hidden Layer	Learning Rate	Learning Times	Training Accuracy (%)
2009-05-09	2	30	30	2.1	15,000	98.9
2009-08-24	2	30	32	2.3	5000	98.1
2010-05-25	2	31	28	2.7	8000	90.7
2010-08-10	2	30	15	3.1	10,000	92.8
2010-09-11 (2 m)	2	31	28	2.1	9000	98.1
2010-09-11 (8 m)	2	30	30	2.6	14,000	87.8
2010-11-21	2	30	29	2.3	14,000	86.1
2011-05-08	2	28	17	1.8	15,000	99.8
2011-08-17	2	30	29	2.3	15,000	89.3
2011-10-24	2	30	32	2.9	15,000	89.6

. Table 3 indicates that the average precision of the interpretation training of each satellite image is approximately 93%.

To confirm the accuracy of the interpretation results, based on the studies of Chen et al. [10,13], the study randomly selected 25 points on the satellite image for the verification of each interpretation factor, supplemented by aerial photographs or on-site survey data and used the aforementioned accuracy evaluation method to evaluate the accuracy of image interpretation and classification. For satellite images of the study area after the Morakot invasion in 2009 (24 August 2009), the image interpretation accuracy is shown in

Table 4. Error matrix of satellite imagery interpretation (24 August 2009).

	Building	Bare Land	Watershed	Road	Forest	River Course	Grassland	Orchard	Paddy Field	Total	UA (%)
Building	23	0	0	2	0	0	0	0	2	27	85.1
Bare land	1	25	1	4	0	2	0	3	4	40	62.5
Watershed	0	0	24	0	0	0	0	0	0	24	100.0
Road	1	0	0	14	0	3	0	0	3	21	66.6
Forest	0	0	0	0	25	0	2	3	0	30	83.3
River course	0	0	0	1	0	20	0	0	0	21	95.2
Grassland	0	0	0	0	0	0	22	0	0	22	100.0
Orchard	0	0	0	0	0	0	1	19	0	20	95.0
Paddy field	0	0	0	4	0	0	0	0	16	20	80.0
Total	25	25	25	25	25	25	25	25	25	225
PA (%)	92.0	100.0	96.0	56.0	100.0	80.0	88.0	76.0	64.0	Kappa index = 0.82 Overall accuracy = 83.6%

, and the results of the satellite image interpretation before and after Typhoon Morakot are shown in Figure 2. Table 4 shows that the consistency coefficient (kappa index value) of the interpretation results of the satellite images was 0.82, and the OA was 83.6%. In addition, this study interpreted a total of 10 satellite images before and after six typhoons or rainfall events in the study area from 2009 to 2011. As shown in

Table 5. Summary of interpretation results of satellite images before and after rain.

Rainfall Event	Date (before/after)	Overall Accuracy (%)	Kappa Index
I	2009-05-09 (before typhoon)	81.0	0.8
	2009-08-24 (after typhoon)	83.6	0.82
II	2010-05-25 (before rainfall)	76.4	0.75
	2010-08-10 (after rainfall, before typhoon)	75.2	0.73
III	2010-09-11—2 m (after typhoon)	80.4	0.79
IV	2010-09-11—8 m (before typhoon)	81.0	0.8
	2010-11-21 (after typhoon)	82.0	0.81
V	2011-05-08 (before typhoon)	76.9	0.75
VI	2011-08-17 (after typhoon Meari, before typhoon Nanmadol)	82.8	0.82
	2011-10-24 (after typhoon)	81.8	0.81

, the average kappa index value was approximately 0.8, and the OA of interpretation was approximately 80.1%, indicating that the image interpretation and classification results have a high degree of accuracy.

5. Discussion

To explore the relationship between the location, scale, and spatial distribution characteristics of the disaster caused by the landslide, this study explored the relationship between landslide and the slope area characteristics by using the relevant quantitative data on landslide and location data of the natural environment. The relationship between landslide and slope factor characteristics is described in the subsequent sections.

5.1. Relationship between the Change in Bare Land and Topographic Location

This study adopted the methods of Meunier et al. [30], Chue et al. [14], and Tseng et al. [31], which were based on the interpretation results of satellite images before and after six typhoons or rainstorms in different years. We extracted the location of the bare ground and calculated the distance between the highest point of the exposure range and the nearest ridge top (dr), between the lowest part of the exposure range and the nearest stream distance (ds), and between the ridge top and the stream distance (dt). This study evaluated the relationship of (dr/dt) and (ds/dt) before and after six rainfalls with the change in the size of the landslide area [14], and the results are shown in Figure 3a–f. The size of the circles in the figure represents the size of the exposed area.

Figure 3 presents the distribution of landslide points caused by six rainfall events from 2009 to 2011. A comparison of the number and area of bare sites before and after each rainfall event (

Table 6. Number and area of landslides before and after different rainfall events in the study area.

Rainfall Event (Number)	Number of Landslides		Total Landslide Area (m2)
	Before Rain	After Rain	Before Rain	After Rain
2009 Typhoon Morakot (I)	114	195	406,890	2,053,415
2010-07-27 rainfall (II)	115	121	1,309,168	1,365,362
2010 Typhoon Meranti (III)	121	134	1,365,362	1,697,533
2010 Typhoon Fanapi (IV)	154	168	1,831,335	1,887,852
2011 Typhoon Meari (V)	91	143	1,278,568	1,770,650
2011 Typhoon Nanmadol (VI)	143	175	1,770,650	2,188,420

) revealed that after each rainfall event, irrespective of the number or area of the exposed area, the landslide points were significantly increased in those areas compared with before the rainfall event. The difference in the distribution of bare land before (Figure 3a—left) and after (Figure 3a—right) the rainfall event of Typhoon Morakot in 2009 was the largest, with a large area of bare land after the event. Chue et al. [14] proposed that before each rainfall event, the landslide points can be predicted using the differences in the spatial distribution of the bare land before and after the rainfall events in different years. During the interval between the end of the rainy season in a year (Figure 3a, right, Figure 3d, right) and the beginning of the next year’s rainy season (Figure 3b, left, Figure 3e, left), the exposed areas of the slopes in the study area were decreased, and obvious vegetation restoration was observed.

The changes in each exposed point are illustrated in Figure 3. As shown in Figure 3a–right, after the rainfall before Typhoon Morakot, the exposed area in the study area is distributed more toward the stream. After the rain, the overall distribution of the bare land is fairly uniform, with a larger area of bare land near the stream. Figure 3b show that the small area close to the top of the ridge was bare land; however, after the rain on 27 July 2010, many places in the adjacent stream were exposed, and thus, the area of bare land increased. Figure 3c, right also shows that after Typhoon Meranti, although the bare land area increased, the position of the bare land did not change considerably. As seen from the spatial distribution of the exposed areas before and after Typhoon Fanapi (Figure 3d), after rainfall, the exposed areas were distributed toward the ridge top. This phenomenon is similar to that observed after Typhoon Meari (Figure 3e) and that observed after Typhoon Nanmadol (Figure 3f). Increasing bare land areas developed near the stream after Typhoons Meari and Nanmadol.

Figure 4 and Figure 5 show the relationship between the average EAR of each field from 2009 to 2011 and the increase in the numbers and areas of landslide. After six rainfalls in the study area, the number of landslides and the landslide area increased with an increase in the average EAR. The study used polynomial trend lines to fit the data. The results revealed that the coefficient of determination between the average EAR and the increase in the number of landslides and that between the average EAR and the increase in the landslide area was 0.8329 and 0.9171, respectively.

5.2. Relationship between Bare Area, Quantity, and Rainfall

The average EAR, bare land area, and the number of bare lands after each rainfall are shown in

Table 7. Number of bare lands and bare land area after each rainfall in the study area.

EAR (mm)	Bare Land Area (km2)	Number of Bare Ground after Each Rainfall							Total Exposed Area (m2)
		I	II	III	IV	V	VI	Total
0–100	0–0.001	0	0	18	0	0	0	18	12,741
	0.001–0.01	0	0	75	0	0	0	75	292,755
	0.01–0.05	0	0	36	0	0	0	36	887,110
	0.05–0.1	0	0	4	0	0	0	4	281,727
	0.1 or more	0	0	1	0	0	0	1	223,200
101–200	0–0.001	0	17	0	12	0	0	29	16,216
	0.001–0.01	0	56	0	26	0	0	82	9630
	0.01–0.05	0	30	0	19	0	0	49	1,239,202
	0.05–0.1	0	2	0	6	0	0	8	580,477
	0.1 or more	0	1	0	1	0	0	2	381,757
201–300	0–0.001	0	3	0	21	15	26	65	35,261
	0.001–0.01	0	10	0	31	34	101	176	660,141
	0.01–0.05	0	2	0	10	13	39	64	1,390,260
	0.05–0.1	0	0	0	0	2	7	9	652,919
	0.1 or more	0	0	0	0	1	2	3	703,239
301–400	0–0.001	0	0	0	13	6	0	19	10,792
	0.001–0.01	0	0	0	20	51	0	71	273,516
	0.01–0.05	0	0	0	7	17	0	24	466,976
	0.05–0.1	0	0	0	1	3	0	4	271,256
	0.1 or more	0	0	0	0	1	0	1	217,214
1001–1100	0–0.001	19	0	0	0	0	0	19	8132
	0.001–0.01	69	0	0	0	0	0	69	278,045
	0.01–0.05	30	0	0	0	0	0	30	720,793
	0.05–0.1	4	0	0	0	0	0	4	270,807
	0.1 or more	2	0	0	0	0	0	2	402,888
1101–1200	0–0.001	21	0	0	0	0	0	21	10,242
	0.001–0.01	43	0	0	0	0	0	43	156,679
	0.01–0.05	5	0	0	0	0	0	5	67,786
	0.05–0.1	2	0	0	0	0	0	2	138,043
	0.1 or more	0	0	0	0	0	0	0	0

. In addition to the extreme rainfall during Typhoon Morakot in 2009, this study divided the average EAR of the study area after the other five rainfalls in the numerical range 0–100, 101–200, and 201–300 mm and plotted the values against the corresponding total exposed area, as shown in Figure 6. The trend line of the exponential relationship was used to fit the data. The results revealed that the area of the exposed area in the study area increased with an increase in the average EAR after various rains, with R² = 0.98:

$$\mathrm{Total} \mathrm{exposed} \mathrm{area} \left(m^2\right)=10+0.6\cdot EXP\left(0.3534 \mathrm{average} \mathrm{EAR} \mathrm{level}\right)$$

(14)

5.3. Relationship between the Number and Area of Bare Lands after Each Rainfall and the Degree of Slope Disturbance

To explore the relationship between the number and area of bare lands after each rainfall and the degree of slope disturbance, this study obtained surface information through the aforementioned satellite image interpretation and classification and referenced the results of Chen et al. [10,13] to establish an index for land disturbance condition (I_LDC). First, six factors, namely vegetation cover (including grasslands and forests) rate, farmland (paddy fields and dry fields), planting rate, fruit tree planting rate, bare land density, building density, and road density, were selected as the slope disturbance factors affecting landslide occurrence. In this study, ArcGIS Spatial Analyst was used to input the results of the interpretation and classification of satellite images in the study area into the basic grid and calculate the proportion of the area occupied by each slope disturbance factor in each grid to be used as quantitative indicators of each slope disturbance factor.

Based on the studies of Chen et al. [10,13], this study defines I_LDC as follows:

$$I_{LDC}=\sum W_{LDC}\times R_{LDC},$$

(15)

where W_LDC is the evaluation index value of the slope disturbance factor in each basic grid and R_LDC is the proportion of the area of each slope disturbance factor in the basic grid area. The value of W_LDC is based on the results obtained and modified from Chen et al. [10] and is provided in

Table 8. Evaluation index values of the slope disturbance factors that affect landslide occurrence

Slope Disturbance Factor	Bare Density	Road Density	Building Density	Fruit Tree Planting Rate	Farmland Planting Rate	Vegetation Cover Rate
Score	6	5	4	3	2	1

.

After the rainfall in each field, the relationship between I_LDC and the exposed area of the overall landslide point is, as shown in Figure 7. The results revealed a positive relationship between the bare land area and I_LDC, except after the extreme Typhoon Morakot rain. The higher the value of I_LDC, the higher the exposed area is. This result indicated that the impact of extreme rainfall similar to that during Typhoon Morakot on the exposed area might be greater than the impact of slope disturbance.

The aggregate statistics of the slope disturbance in the study area and the area and number of bare lands after each rainfall are shown in

Table 9. Slope disturbance in the study area and the area and number of bare lands after each rainfall.

ILDC	Bare Area (km2)	Number of Bare Ground after Each Rain							Total Exposed Area (m2)
		I	II	III	IV	V	VI	Total
0–1	0–0.001	0	0	0	14	0	0	14	6010
	0.001–0.01	0	0	2	4	0	0	6	19,389
	0.01–0.05	0	0	0	0	0	0	0	0
	0.05–0.1	0	0	0	0	0	0	0	0
	0.1 or more	0	0	0	0	0	0	0	0
1.01–2	0–0.001	20	0	10	12	5	16	63	34,918
	0.001–0.01	45	0	24	18	9	25	121	356,405
	0.01–0.05	10	0	3	1	0	4	18	344,875
	0.05–0.1	1	0	0	0	0	0	1	67,160
	0.1 or more	0	0	0	0	0	0	0	0
2.01–3	0–0.001	14	4	8	9	15	9	59	33,895
	0.001–0.01	49	12	43	27	59	57	247	980,149
	0.01–0.05	16	5	19	11	13	11	75	1,602,908
	0.05–0.1	4	0	1	0	1	0	6	399,804
	0.1 or more	2	0	0	0	0	0	2	402,888
3.01–4	0–0.001	6	6	0	4	1	0	17	8811
	0.001–0.01	15	31	6	14	15	16	97	464,618
	0.01–0.05	7	16	11	12	17	16	79	1,891,087
	0.05–0.1	1	1	2	2	4	2	12	826,064
	0.1 or more	0	0	0	0	1	1	2	335,312
4.01–5	0–0.001	0	8	0	2	0	1	11	6356
	0.001–0.01	3	21	0	5	2	3	34	116,412
	0.01–0.05	2	10	2	5	0	8	27	709,678
	0.05–0.1	0	1	1	2	0	4	8	607,174
	0.1 or more	0	1	1	1	1	1	5	1,190,097
5.01–6	0–0.001	0	2	0	1	0	0	3	981
	0.001–0.01	0	2	0	8	0	0	10	25,280
	0.01–0.05	0	1	1	6	0	0	8	192,468
	0.05–0.1	0	0	0	2	0	1	3	196,710
	0.1 or more	0	0	0	0	0	0	0	0
6.01–7	0–0.001	0	0	0	4	0	0	4	2414
	0.001–0.01	0	0	0	2	0	0	2	11,939
	0.01–0.05	0	0	0	1	0	0	1	31,111
	0.05–0.1	0	0	0	1	0	0	1	98,318
	0.1 or more	0	0	0	0	0	0	0	0

. Figure 8 presents a plot of the I_LDC and the ratio of the area and number of bare lands after the six rainfalls. The trend line of the exponential relationship was used to fit the data. The results revealed that the ratio of the area of bare lands to the number of bare lands in the study area after various rainfall events increases with an increase in I_LDC, with R² = 0.72:

$$\frac{\mathrm{Bare} \mathrm{area} \left({\mathrm{km}}^2\right)}{\mathrm{Number} \mathrm{of} \mathrm{bare} \mathrm{ground}}=1681.2\cdot EXP\left(0.4343\cdot I_{LDC}\right)$$

(16)

5.4. Relationship between the Scale of New or Second Landslide, Rainfall, and I_LDC

After the rainfall in each field, new landslide sites were observed, and the related factors, such as the location of the landslide, scale of the landslide (including new landslide or second landslide), rainfall, and I_LDC were evaluated.

5.4.1. Relationship between the Landslide Location and the Corresponding Landslide Area

By using the interpretation of the satellite images and the extraction of the exposed area, as described in Section 5.1, after the rainfall in the study area during 2009–2011, we calculated the distance between the highest point of the exposed area and the nearest ridge top (dr), the distance between the lowest point of the bare area and its nearest stream (ds), and the distance between the ridge top to the creek (dt). The correlation between the landslide location and scale after each typhoon rainfall was evaluated, and the landslides (including new landslide and second landslide) induced by the rainfall event in the study area were assessed, as shown in Figure 9a–e. With the landslide points induced by the extreme rains of Typhoon Morakot in 2009 as the benchmark, we plotted the variation of the new landslide point (left) against the second landslide point (right) of the other five rainfall events (Figure 9). In Figure 9, the size of the circle indicates the size of the landslide area, and the plots represent the 3-year period from 2009 to 2011.

Figure 9 shows that the number of second landslide points is considerably greater than the number of new landslide points. This result indicates that the historical landslide points in the study area will experience rainfall events in the future and thereby experience landslides; moreover, the scale of the second landslide is greater than that of the new landslide. In addition, Figure 9 shows that most of the new landslide points are biased toward the ridge top (dr/dt ≈0), and the second landslide points of larger landslide scale are biased toward the stream (ds/dt ≈0).

5.4.2. Relationship between Landslide Area and Location, Rainfall, and Slope Disturbance

This study was based on the aforementioned calculation of the EAR and I_LDC. The relationship between the EAR, I_LDC, and landslide area (including new landslide and second landslide) after each rainfall in the study area is plotted in Figure 10a–f. The size of the circle in the figure represents the varying size of the landslide area; the black and red circles represent the newly added points and the second landslide points before and after each rainfall, respectively. Irrespective of the rainfall, when the I_LDC was large, the scale of the second landslide was also large. In addition to the extreme rain of Typhoon Morakot, the number and area of the second landslide after the rainfalls were greater than those of the new landslide.

5.4.3. Relationship between Variation in Second Landslide Scale, Rainfall, and Slope Disturbance

The relationship between the EAR, I_LDC, and variation in the area (the increased landslide area) of the second landslide after each rainfall from 2009 to 2011 is plotted in Figure 11a–f. In the figure, the increase in the landslide area at the second landslide point is represented by the size of the circle. The results revealed that irrespective of whether the area of the second landslide at the second landslide point was related to EAR or I_LDC, a positive relationship was noted overall. In the landslide sites of the study area, at a large EAR or I_LDC, the area of the second landslide site also increased.

6. Conclusions

We targeted the Lao-nong River Watershed in Southern Taiwan as the study area, and the study period was from 2009 to 2011, starting with the Typhoon Morakot invasion in 2009. During the 3-year study period, the study area experienced six rains, including five typhoons and one heavy rain. The study used the GANN and texture analysis on the GIS platform for satellite image classification and interpretation to analyze land-use change, characteristics of the landslide area, and location of second landslides to provide a reference for devising preventive and response countermeasures. In addition, this study explored the impact of rainfall on the second landslide and evaluated the relationship between the various rainfall events and the occurrence, scale, and regional characteristics of the landslide. The study conclusions are summarized as follows:

(1)

For the interpretation and classification of high-resolution satellite images, this study used the GANN combined with texture analysis. The OA and consistency coefficient values of the interpretation results revealed that the satellite image interpretation before and after each rainfall in the research area achieved medium to high accuracy;

(2)

A comparison of the number and area of the exposed areas before and after the six rainfalls revealed that the number or area of the bare land in the research area in each field significantly increased after the rainfall than before the rainfall. The distribution of bare land before and after Typhoon Morakot was the largest. In addition, after each rainfall, the number of bare lands and bare land areas increased with an increase in the average EAR. When the data were fitted with a polynomial trend line, the coefficient of determination between the average EAR and the increases in the number of bare lands and that between the average EAR and the increase in the landslide area was approximately 0.83 and 0.92, respectively;

(3)

In addition to the extreme rainfall during Typhoon Morakot in 2009, this study divided the average EAR after each rain into three levels in sequence and used the trend line of the exponential relationship to fit the bare land data. The results revealed that after each rainfall in the study area, the bare land area increased with an increase in the average EAR value, and the coefficient of determination of trend line reached 0.98;

(4)

The relationship between I_LDC and the bare land area after each rainfall indicated that except for the extreme Morakot rains, the greater the degree of slope disturbance was after rain, the greater the area of the exposed slope was. This result also indicated that when extreme rainfall similar to Typhoon Morakot strikes, the impact of rainfall on the bare land area may be greater than the impact of slope disturbance. In addition, the results of the joint mapping study after the rainfall in each field revealed a positive relationship between the bare land area and I_LDC;

(5)

The relationship between the I_LDC in the study area and the ratio of the area of bare land to the amount of bare land after each rainfall indicated that the ratio of the area of bare land to the number of bare lands after each rainfall increased with I_LDC;

(6)

The results of the rainfall-induced new landslide and second landslide in each field revealed that except for the number of new landslide points induced by the extreme rainfall event during Typhoon Morakot, which was considerably higher than the number of second landslide points, for the remaining landslides induced by rainfall, the number of second landslide points was higher than the number of new landslide points, and the area of the second landslide point was also greater than that of the new landslide point. In addition, despite the rainfall, the larger the slope disturbance, the larger the scale of the second landslide was. Consequently, more new landslide points were biased toward the ridge crest, whereas the second landslide points with larger landslide scales tended to develop toward the stream;

(7)

After rainfall in each field, the relationship between the EAR at the point of the second landslide, I_LDC, and re-increased area of landslide indicated that overall, a positive relationship was noted between the increased area of the landslide at the second landslide point and the EAR or I_LDC. With an increase in the EAR on the slope in the study area or the slope disturbance, the area of the landslide at the second landslide point also tended to increase.