The Mid-Niigata earthquake occurred on 23 October 2004. The Japan Meteorological Agency (JMA) measured the magnitude of the mainshock as 6.8, the epicenter is located at 37°18′16.56″N, 138°50′10.32″E, and the focal depth is about 13.1 km (Chigira and Yagi, 2006; Kokusho et al., 2011). Within 3 days after the mainshock, more than 900 landslides were induced by the earthquake (Chigira and Yagi, 2006; Kokusho et al., 2014). After the earthquake sequences, a very small surface rupture was found along a previously unmapped northern extension fault zone. The length of the surface rupture was about 1 km (Maruyama et al., 2007). The surface slip of the Mid-Niigata earthquake event was also very small (<20 cm of vertical displacement). In addition, the surface rupture zone is also far away from the epicenter zone, where seismic landslides have concentrated distribution (Sato et al., 2005), i.e., the study area of seismic landslide susceptibility did not contain the surface rupture zone. So in this study, we consider that the surface rupture zone has little effect on the formation of seismic landslides and regard the earthquake as a buried-rupture earthquake.

In this study, the assessment of seismic landslide susceptibility is performed on two scales, the large area and the epicenter area. As shown in Figure 1, the large-scale area is 22 km wide (east to west) and 40 km long (north to south). The total area of the large-scale area is about 880 km². The epicenter area is 7 km long (north to south) and 9 km wide (west to east). The total area is about 56 km². The epicenter area is located in the bordering area between Nagaoka City and Ojiya City.

Many methods have been utilized to set up landslide inventory maps, including satellite image interpretation, aerial photography, field survey, and historical landslide records (Vařilová et al., 2015). In this research, the landslide inventory map was interpreted from satellite image data and then checked by field survey data (Kokusho 2008; Kokusho et al., 2009). As shown in Figure 1A, a total of 957 landslide locations were recorded in the large-scale area, most of which are distributed in the mountainous area around the epicenter area and spread to the northeastern mountainous area. There are also some landslides located in the eastern and southern mountain areas. The landslide inventory map of the epicenter area is also shown in Figure 1B.

The factors that affect the occurrence of earthquake-triggered landslides usually include geology, topography, hydrology, climate, human activities, and earthquake-related parameters. Based on the availability of data and impacted factors used in previous studies (Reichenbach et al., 2018), seven landslide influencing factors (lithology, elevation, slope, slope aspect, surface curvature, peak ground acceleration, and the distance from the road) were taken into consideration for landslide susceptibility analysis for the large-scale area. In the later analysis in the epicenter area, coseismic ground deformation was added as an influencing factor. The scales of the landslide influencing factors are shown in Table 1.

Lithology directly determines the physical and mechanical properties of the slope, which have a direct impact on slope stability. The lithology maps in the large area and the epicenter area are shown in Figure 2A and Figure 3A. The elevation also affects the occurrence of seismic landslides (Hasegawa et al., 2009). The elevation has been regarded as a key factor determining the gravitational potential energy of the terrain. Elevation also indicates the amplification effect of seismic-induced landslides. The elevation maps of the large-scale area and the epicenter area are shown in Figure 2B and Figure 3B, respectively. The slope angle has a direct impact on slope stability that determines the ratio of the anti-sliding force to sliding force. The slope angle maps of the large-scale area and epicenter area are shown in Figure 3C and Figure 3C, respectively. According to previous studies (Xu et al., 2012a; Pham et al., 2016; Hong et al., 2017), the slope aspect is divided into nine groups. The slope aspect maps of the large study area and the epicenter area are shown in Figure 2D and Figure 3D, respectively, and the P and FL mean the flat area. The surface curvature distributions in the large-scale and epicenter areas are shown in Figure 2E and Figure 3E, respectively.

The peak ground acceleration (PGA) of an earthqauke is the maximum absolute value of the acceleration of the surface soil in the earthquake (Tian et al., 2019; Xu and Xu, 2013; Xu et al., 2013; Li et al., 2013; Xu et al., 2012a,b). The PGA was obtained from inversion analysis based on the K-NET and KiK-net observations. The resolution of PGA data is 350 m. The distribution of peak accelerations in the large-scale area and the epicenter area is shown in Figure 2F and Figure 3F, respectively. Human activities have also greatly impacted the topography features. Road construction not only produced a new steep cutting slope but also caused great disturbance to the original slope. Therefore, the distance to the road is taken into account in the assessment of landslide susceptibility. In this study, the locations of high-grade roads like expressways were interpreted from the satellite image. The distances to the road map were divided into seven classes (0–50, 50–100, 100–200, 200–300, 300–400, 400–500, and >500 m). The distances to road maps of the large-scale area and the epicenter area are shown in Figure 2G and Figure 3G, respectively.

For earthquakes with surface ruptures, previous research studies show that there is a clear connection between landslide distribution and the distance to the rupture zone (Xu et., al 2012b; Xu, 2014), which means the distance to the surface rupture could be used as an influencing factor. However, for buried-rupture earthquakes, as very short or no surface rupture is exposed, it is difficult to establish the relationship between the distribution of landslides and surface rupture. Therefore, it is necessary to introduce new influencing factors to improve the accuracy of landslide susceptibility analysis for buried-rupture earthquakes.

Coseismic ground deformation characterizes the absolute permanent ground deformation before and after the earthquake, and it has been demonstrated that there is a good correlation between landslide distribution and the values of coseismic ground deformation (Chang et al., 2005; Zhao et al., 2014). Therefore, coseismic ground deformation could make up for the disadvantage of losing surface rupture in the assessment of seismic landslide susceptibility to a certain extent. Coseismic ground deformation can be obtained by decomposed high-resolution DEM before and after the earthquake (Zhang et al., 2010; Zhao et al., 2012). Figure 4 shows the description of landform changes in Lagrangian and Eulerian manners. Supposing that a small patch i of the ground surface with one particular node mapped on it is inclined in East-West (x) and North-South (y) directions, Δz _i ^e is expressed in terms of the Lagrangian vector {Δx _i ^l Δy _i ^l Δz _i ^l } of the movement of the patch as Δ Z i e = { t X , i   t y , i   1 } ⋅ { Δ X i l   Δ y i l   Δ Z i l } T , (1) where t _x,i and t _y,i are tangents of the patch plane in the x and y directions, respectively. Taking three adjacent patches, i1, i2, and i3 in a triangle and using the displacement of its center {Δx _i ^l Δy _i ^l Δz _i ^l } as the representative displacement vector of the triangle, the following simultaneous equations are to be satisfied: { Δ z i 1 e Δ z i 2 e Δ z i 2 e } = [ t x , i 1     t x , i 1     1 t x , i 1     t x , i 1     1 t x , i 1     t x , i 1     1   ] { Δ x i l Δ x i l Δ x i l } = T • { Δ x i l Δ x i l Δ x i l } . (2)

An assumption that the triangle undergoes a rigid body translation is used in the aforementioned formulation. The inclination of the moving plane (plane i1) is essential for calculating t _x,i and t _y,i . Suppose the equation of the moving plane is expressed as z = a x + b y + c , (3) where a= t _x,i1=tanθ _x,i1 and b= t _y,i1=tanθ _y,i1

Zhao et al. (2012) provided a more rigorous solution method, including the definition of a nominal plane, the improvement of DEM comparability, and matrix condition test. In this study, we used the method proposed by Zhao et al. (2012) to calculate coseismic ground deformation. It is to be noted that the decomposition algorithm requires high-resolution (2 m) DEM. Thus, coseismic ground deformation is added as an influencing factor for the epicenter area only. The DEMs before the earthquake were derived from aerial photos shot by the Geospatial Information Authority of Japan in 1975 and 1976. Triangulation points prepared for road construction in 1986 were then used to orthogonalize and digitize these photos by Aero Asahi (2004). The standard deviation of the digitized aero photos from the triangulation data is 0.589 and 0.517 m in the horizontal and vertical directions, respectively. The DEMs after the earthquake were prepared by airborne LiDAR scanning conducted by Nakanihon Air Service on 24 Oct 2008, the second day of the mainshock and three major aftershocks. After scanning and post-processing standard courses, calibration courses were performed. The average differences for the x coordinate, y coordinate, and elevation are 0.09, 0.2, and 0.06 m respectively, showing good repeatability of scanning. Both sets of DEMs have a resolution of 2 m×2 m. Since the two sets of DEMs were prepared in different ways, it is not appropriate to directly compare them in the calculation. The contour line derived from the 1975 DEMs is quite smooth while that from the 2004 DEMs is curly and tilted owing to the high spatial frequency components. Therefore, we used smoothed elevations instead of the original ones for the 2004 DEMs by substituting the x and y coordinates of each point into the equation of its nominal plane. Given that three adjacent points undergo a rigid body translation movement, their coseismic deformation can be decomposed by solving simultaneous equations. A cut-off window was selected according to regional geology to obtain the tectonic displacements. The properly defined nominal plane could also fully or partially eliminate terrain changes owing to human activities during the time gap between two sets of DEMs since the size of the smooth window is larger than all man-made changes. In addition, to verify the accuracy of the calculated coseismic ground deformation, the calculated displacements were also compared with those at points of triangulation. In total, 11 available triangulation points buried on roads were used. The compared results showed that the difference between the observed displacement and calculated displacement was small, which demonstrated that the calculated coseismic ground deformation is accurate. Details can also be found in the studies by Konagai et al. (2009) and Zhao et al. (2012). The distribution of coseismic ground deformation in the epicenter area is shown in Figure 5.

The orientation of the computed coseismic ground deformation could be divided into two directions: lateral components and vertical components. Zhao et al. (2012) compared the location of earthquake-triggered landslides with the displacement field of lateral components and vertical components, respectively. The results showed that landslide clusters were found within large lateral deformation regions, while landslides seem to be off where the vertical displacement is large. Therefore, in this study, only lateral deformation is used. Figure 5 also shows the distribution of the absolute value of the lateral ground deformation.

In this study, the numbers of landslide points and non-landslide points are sampled at a ratio of 1:1.2 for the large-scale area. A total of 1,117 non-landslide points’ data were randomly selected in the non-landslide area. Subsequently, 70% of the landslide points and non-landslide points were selected randomly from the landslide inventory map as the training dataset, with the rest as the testing dataset. To obtain optimum results, we randomly selected the sample points (landslides points and non-landslide points) 10 times, respectively. For different selections, the training and testing samples are different, but the numbers of sample points are the same. In the epicenter area, as the used method to calculate coseismic ground deformation needs high-resolution DEM, the whole epicenter area was converted into 2 m pixels. The total number of pixels is 555,324, and the number of seismic landslide pixels is 45,852. Similarly, 70% of the landslide pixels and non-landslide pixels were selected randomly as the training dataset, with the remaining 30% as the testing dataset.

Logistic regression is suitable for describing the relationship between the categorical outcome (landslide or non-landslide) and input variables (landslide-affecting factors). The principle of the LR is to analyze the spatial relationship between the landslide-affecting factors and the occurrence of a landslide. The results of the regression usually can be interpreted as the probability which is constrained in the interval between 0 and 1.

The LR is indicated by an equation of the form Y = f ( P ) = ln ( P 1 − P ) = β 0 + β 1 X 1 + β 2 X 2 + ⋅ ⋅ ⋅ + β n X n , (4) where Y represents outcome variables (landslide or non-landslide), X = X ₁, X ₂ … Xn represents input variables, n is the nth landslide-affecting factor, β ₀ is the intercept condition, and β ₁, β ₂ … β _n are the regression coefficients (Tu, 1996).

The SPSS 10.0 was used to conduct the LR analysis to predict the correlation between the occurrence of landslides and landslide-affecting factors. The regression coefficients were then obtained.

The probability of a landslide event (P) can be determined from the following equation: P = P ( Y / X ) = e β 0 + β 1 X 1 + β 2 X 2 + ⋅ ⋅ ⋅ + β n X 1 + e β 0 + β 1 X 1 + β 2 X 2 + ⋅ ⋅ ⋅ + β n X . (5)

The probability values change from 0 to 1, with 0 indicating a 0% probability of landslide occurrences and one indicating a 100% probability.

An ANN model has many advantages compared with other models (Yi̇lmaz, 2009a). An ANN could process imprecise and fuzzy data without any assumptions. The ANN model with the most frequently used back-propagation BP algorithm (Pradhan and Lee, 2010b) is used in this study.

The model mainly consists of one input layer, several hidden layers, and an output layer. There are usually two stages for using an ANN, the training stage and classifying stage. During the training stage, the hidden and the input layer neurons handle their inputs by a corresponding weight, sum the product, and then deal with the sum using a nonlinear transfer function to generate a result. During the classification period, the ANN predicts a target value by adjusting the weights in accordance with the errors between the actual output values and the target output ones and makes the difference minimum.

In this study, the number of hidden layer nodes is calculated by the following equation (Yi̇lmaz, 2009a): N h = 2 N i + 1 , (6) where N _i is the number of input nodes and N _h is the number of hidden nodes.

Then, a three-layer network with one input layer (7 neurons), one hidden layer (15 neurons), and one output layer was used in the large-scale area. In the epicenter area, a three-layer network consisting of one input layer (8 neurons), one hidden layer (17 neurons), and one output layer was utilized. It is important to decide the initial weight range influencing the convergence of the model. In this study, the initial weights were randomly selected from a small range of [-0.25 to 0.25] as proposed by Yi̇lmaz, 2009b. For the hidden layer, the activation function is the hyperbolic tangent, and the activation function is the softmax function for the output layer. We did not use an optimization algorithm for the ANN model.

The SVM model employs nonlinear transformations of the covariates into a higher dimensional feature space. The two main principles of the SVM are the optimal classification hyperplane and the use of a kernel function (Yao et al., 2008).

The details of a two-class SVM model are described as follows. Given a set of linear separable training vectors x _i (i=1,2 … n) that consist of two categorical outcomes (landslide or non-landslide denoted as y= ±1), the purpose of the SVM is to find an n-dimensional hyperplane differentiating the two categories by the maximum gap.

Mathematically, the gap 1 2 ‖ w ‖ 2 could be a minimized subject to the following constraints: y i ( ( w ⋅ x i ) + b ) ≥ 1 , (7) where ||w|| is the norm of the normal of the hyperplane, b is a scalar base, and (·) denotes the scalar product operation. Using the Lagrangian multiplier, the cost function can be defined as L = 1 2 ‖ w ‖ 2 − ∑ i = 1 n λ i ( y i ( ( w • x ) + b ) ) ≥ 1 , (8) where λ _i is the Lagrangian multiplier. The solution can be obtained by the dual minimization of Eq. (8) with respect to w and b

In this study, the two-class SVM method was used due to its good performance in landslide susceptibility analysis (Yao et al., 2008; Yi̇lmaz, 2010).

In this study, the performances of the three models (LR, ANN, and SVM) for the large-scale area were validated using a receiver operating characteristic (ROC) curve. The area under the curve (AUC) indicates how good the statistical model is. It means the model has perfect performance when the AUC value equals to 1. A higher AUC value indicates better performance of the statistical model.

Because each sample dataset is selected randomly, the landslide susceptibility calculated by the same model is not the same. To determine the best model, the models are utilized 10 times for analysis of randomly selected datasets, respectively. For different analyses, the training and testing samples are different. For the same analyses, the training samples and testing samples are the same for all three models. The area under the ROC curve (AUC) of each analysis was compared to explore the difference between the three methods. The results are shown in Table 2.

Table 2 shows that the ANN model performed the best among the three models with the highest AUC value and the accuracy of the SVM model was the worst. Based on the maximum AUC values of 10 simulations, the ANN simulation result was also the best (83.6%). The average value and variance of the ANN model were 82.5 and 0.44%, respectively, which was better than those of the LR and SVM models. It means the robustness of the ANN model is better than that of the LR and SVM models.

Yi̇lmaz (2009a) used three models including frequency ratio (FR), ANN, and LR to generate the landslide susceptibility maps of Kat County (Tokat–Turkey). The result showed that the ANN model performed better than the other models. Yi̇lmaz (2010) utilized four different models, namely, conditional probability (CP), LR, ANN, and SVM models to assess the landslide susceptibility of Koyulhisar (Sivas, Turkey). The results also showed that the performance of the ANN model was the best. Some other research studies also showed that the ANN model performed more accurately than the other models (Gómez and Kavzoglu, 2005; Yesilnacar and Topal 2005). We consider the ANN model performed better than the other models because it has a good global searching ability and can learn the near-optimum solution without the gradient information of error functions. As there are about a total of 2,000 samples in the calibration and validation sets, large numbers of samples in the calibration stage will lead to sufficient training of the model and establish an appropriate structure of the ANN model. So the ANN model performs well on the condition that those large numbers of samples were available. For any algorithm, the quantity and quality of samples have key impacts on the accuracy of the predicted results the algorithm makes.

In this study, all three models have been used to calculate the landslide susceptibility index (LSI) and then generate the landslide susceptibility maps. There are several mathematical methods including quantiles, natural breaks, standard deviation, equal intervals, and descending area percentage to be reclassified in the LSI (Ayalew et al., 2004). Among the aforementioned methods, the descending area percentage technique is the most widely used. In this study, the descending area percentage technique was used. The landslide susceptibility maps were constructed into four classes: low (40%), moderate (30%), high (20%), and very high (10%). The landslide density was used to assess the performance of landslide susceptibility maps. The landslide density (LD) is defined as the ratio of the number of landslides and the area of each susceptible class.

The calculated landslide densities by using the three different models are shown in Table 3. It can be observed that all maps present good spatial predictions of landslides as landslide density is ascending from very low to very high class (Yi̇lmaz, 2009b). The results using the ANN model show that the very high class contains 42.01% of the total landslides; however, it only covers 9.95% of the total study area, and the LD of the very high class was 4.59. In comparison, the low classes only contain 3.34% landslides; however, it covers 40.15% area and the LD of the low class was 0.09. This indicates that the ANN model performed well in susceptibility classification as it fits well with the landslide inventories.

The landslide susceptibility maps of different methods are shown as Figure 6. The analysis results of LR, SVM, and ANN models are very close. The epicenter area is a very high susceptible area, the northeast and the southwest mountain areas are high and very high susceptible areas, respectively, and the northern plains area is distributed with the low susceptible class. The susceptibility map of the ANN model shows that the high susceptible areas and low susceptible areas are more concentrated into blocks, and zonation produced by the SVM and LR are more dispersed. Overall, all three models could generate reasonable landslide susceptibility maps.

From the landslide susceptibility map of the large-scale area, it is known that the susceptibility level in the epicenter area is generally high. Since it is still too costly to remediate all slopes in the approximately 60 km² area, it is necessary to further evaluate the landslide susceptibility in the epicenter area. It can be seen in Section 4.1 that the ANN model is the most suitable model for landslide susceptibility assessment in this area. Therefore, we only use the ANN model to analyze and evaluate the landslide susceptibility in the epicenter area.

First, to evaluate the significance of landslide susceptibility analysis by considering different scales, the values of the AUC for the epicenter area are calculated in two different conditions. First, we calculate the values of the AUC based on the corresponding calculated LSI of the epicenter area from the large-scale (whole affected area) datasets. Then, the values of the AUC are calculated based on the calculated LSI from the epicenter area datasets. The values of the AUC of the two different conditions are shown in Figure 7. The results show that the AUC is 56.2% based on the calculated LSI from the large-scale datasets, in contrast, the AUC is 72.3% based on the calculated LSI from the epicenter area solely. The results show that it is necessary to assess landslide susceptibility under different scales.

Subsequently, the landslide density and landslide susceptibility map of the epicenter area were obtained as shown in Table 4 and Figure 8. The results show that the very high class contains 40.44% of the total landslides; however, it only covers 8.6% of the epicenter area and the LD of the very high class was 26.54. In comparison, the low classes contain only 5.12% landslides; however, it covers 40.15% area and the LD of the low class was only 0.73. The landslide density increases gradually between the low class and very high class. This indicates that the landslide susceptibility map fits well with the landslide inventories.

As shown in Figure 8, the very high-class area is mainly distributed along the long axis of the ellipse in the east of the study area, and a large amount of deep-seated landslides occurred in this area. The high-susceptibility area is also distributed in the northwestern area. The occurrence possibility of landslides in the central area and southwest plain area is relatively low. Compared with the epicenter area parts in the landslide susceptibility map of the large-scale area, the landslide susceptibility maps obtained by the epicenter area research have a better discrimination degree, which can meet the key prevention and control requirements in the small area.