Planet Image-Based Inventorying and Machine Learning-Based Susceptibility Mapping for the Landslides Triggered by the 2018 Mw6.6 Tomakomai, Japan Earthquake

Abstract

The 5 September 2018 (UTC time) Mw6.6 earthquake of Tomakomai, Japan has triggered about 10,000 landslides with high density, causing widespread concern. We attempted to establish a detailed inventory of this slope failure and use proper methods to assess landslide susceptibility in the entire affected area. To this end we applied the logistic regression (LR) and the support vector machine (SVM) for this study. Based on high-resolution (3 m) optical satellite images (planet image) before and after the earthquake, we delineated 9295 individual landslides triggered by the earthquake, occupying an area of 30.96 km². Ten controlling factors were selected for susceptibility analysis, including elevation, slope angle, aspect, curvature, distances to faults, distances to the epicenter, Peak ground acceleration (PGA), distance to rivers, distances to roads and lithology. Using the LR and SVM, two landslide susceptibility maps were produced for the study area. The results show that in the LR model, the success rate is 84.7% between the landslide susceptibility map and the training dataset, and the prediction rate is 83.9% shown by comparing the test dataset and the landslide susceptibility map. In the SVM model, a success rate of 90.9% exists between the susceptibility map and the test samples, and a prediction rate of 87.1% from comparison of the test dataset and the landslides susceptibility map. In comparison, the performance of the SVM is slightly better than the LR model.

1. Introduction

The 5 September 2018 Mw6.6 earthquake of Tomakomai, Japan has triggered about 10,000 landslides with high density, causing serious economic losses and casualties. In some local areas, the damage caused by this event even exceeded that by greater earthquakes such as the 2008 Wenchuan, China Mw7.9 [1] and 2015 Nepal Mw7.8 [2]. According to reports, this earthquake killed at least 41 persons, including 36 dead by landslides [3]. Hundreds of landslides have broken the ridges and changed the landform there. However, so far no comprehensive landslide inventory and landslide susceptibility assessment related to this event have been reported. To fill this gap, this work attempted to establish a detailed landslide inventory of the Tomakomai earthquake and use effective methods to assess landslide susceptibility in the entire affected area.

Currently, the methods for landslide susceptibility mapping can be classified into two types: Qualitative and quantitative approaches. In the qualitative approach, expert opinions are very vital to estimate landslide potential from intrinsic factors, assigning weights of variables and rating susceptibility assessment are subjective [4]. Thus, the resulting susceptibility assessment is not ideal. In contrast, the quantitative approach can be considered as more objective due to their data-dependent characteristics.

The quantitative approaches include statistical and advanced machine learning approaches. The statistical approaches use a landslide inventory to establish a relationship between landslides and their controlling factors for a sampled area, which is then extended to the total research area to attain the susceptibility assessment [5]. Many scholars have employed these methods to carry out landslide susceptibility assessment with the weight of evidence [6], discriminant analysis [7], binary statistical analysis [8,9,10], and the multivariate statistical model [11,12]. Among them, the LR method has been proved to be more applicable and have achieved good assessment results in several earthquakes cases [13,14,15,16,17]. Meanwhile, with the development of machine learning, the landslide susceptibility assessment based on machine learning is gradually favored by scholars [18], such as neural network [13,19] and SVM [20,21,22]. Among them, the SVM is more powerful to handle high-dimensional and nonlinear problems, and its prediction accuracy is superior to the statistical models in most cases [23,24,25].

At present, most statistical models and advanced machine learning models have been used in non-seismic landslide areas with a small number of landslide samples. While there are few applications in the specific earthquake events with a large number of seismic landslides. And these studies mostly rely on incomplete landslide database to conduct the landslides susceptibility assessment in local quake-affected area, which cannot fully reflect the overall distribution of landslides triggered by a single earthquake [25,26,27].

The aim of this study is to establish a detailed landslide inventory for the 2018 Mw6.6 Tomakomai, Japan earthquake based on high-resolution Planet images (3 m) before and after the earthquake and test the applicability of the LR and SVM to landslide susceptibility mapping in the affected area. Then the application effects of the two models were compared. This study provides successful experience and scientific reference for the selection of remote sensing images, landslide interpretation, selection of landslide influencing factors, and selection of evaluation model in similar studies.

2. Study Area

The Japanese archipelago lies at the interaction between the Pacific, North American, Eurasia, and Philippine plates, where there exists the deepest trench on the Earth’s surface. The North American and Eurasia plates here can be further subdivided into the Okhotsk Sea and the Amur plates. On 5 September 2018 (UTC Time), an Mw 6.6 earthquake occurred in the Oshima Belt region, east of Tomakomai on the island of Hokkaido, Japan (Figure 1). The strata of the whole tectonic area include from Holocene to Jurassic. The detailed description of lithology of different ages can be seen in Table S1. The study area extends on a surface of about 600 km², spanning from 42.6°N to 42.9°N of latitude and from 141.8° to 142.2°. The elevation of the area ranges between 0 and 624 m, with an average value of 166 m. The study area is made up of Miocene sedimentary rocks as its basement, and covered mostly up with air-fall lapilli-sized pumice layers. Surface soil layers covering low to middle mountain ranges are inter-bedded with the pumice and ash. Total thickness of the surface layer is about 4–5 m [3].

At the location of this earthquake, the Pacific plate moves northwestward at a velocity of about 87 mm/yr relative to the North America plate and collides with the North American Eurasia and Philippine plates along the West Pacific subduction zone, causing frequent earthquakes. The focal mechanism solutions of the Tomakomai earthquake (

Table 1. Parameters of the 2018 Mw6.6 Tomakomai, Japan earthquake.

Latitude (°)	Longitude (°)	Dip (°)	Depth (km)	Rake (°)	Mo (Nm)	Var. Red.	Strike (°)
42.6908	142.0067	30; 65	35	59; 107	1e + 19	89.47	134; 349

) indicate that faulting occurred on either a moderately dipping reverse fault striking Northwest, or on a shallow-to-moderately dipping fault striking southeast. In view of these focal mechanism solutions and focal depth (35 km), the earthquake likely represents rupture of a fault within the upper North America plate or Okhotsk microplate, instead of the subduction zone plate boundary interface between the Pacific and North America plates (which is 100 km deep at the location of this epicenter) (https://earthquake.usgs.gov).

Figure 2 and Figure 3 are photos and satellite images of the landslides triggered by the 2018 Tomakomai earthquake. As shown in Figure 2, the density of landslides triggered by this earthquake is high, most of them are shallow landslides, and distributed continuously. Substantial landslides have destroyed the ridges and caused the landform changes (Figure 2a,c). Figure 3a,b show the landslides damaged many houses and farmland, and blocked roads and valleys. Figure 3c shows that the debris area of the landslide was a mixture of the uprooted trees and the pumice and ash precipitate, and spreads over croplands and villages.

3. Data and Method

To map the landslide susceptibility, the following preparations were made: (1) Constructing database: The landslide inventory of 2018 Tomakomai earthquake and influencing factors of landslide susceptibility; (2) model construction: Samples for training and testing LR and SVM modeling; and (3) model validation: Preparing landslide susceptibility maps of LR and SVM, and compare the assessment results. Figure 4 shows the flow chart of these efforts.

3.1. Data Source

3.1.1. Landslide Inventory

A detailed and accurate earthquake-triggered landslide inventory is an important basis for landslide susceptibility assessment [31,32]. With the advance of sensor and space technology, remote sensing is able to acquire detailed temporal and spatial information on landslides on the Earth’s surface. This study used the satellite images of post-earthquake with high-resolution (3 m) (Planet satellite), which were acquired within 5 days after the earthquake [33]. These images cover the entire earthquake affected area, and have low cloud coverage (Figure 5b). The series of pre-earthquake ortho-images also come from the high-resolution (3 m) optical satellite images (Planet satellite) to ensure that the landslides existing before the earthquake are not recognized as co-seismic landslides (Figure 5a). The processed images are overlaid on the Google Earth Platform. Thus, the images can provide a stereo-like perspective view, and are very convenient to compare them before and after the seismic event [34].

We followed the principles of image interpretation: (1) About the landslide scale, we interpret it as long as it is a landslide that can be distinguished from the image. (2) About the landslides before the earthquake, we think that they were caused by this earthquake if their shapes are not consistent before and after the earthquake. (3) In general, the areas with flow-textures, light tones, and non- or destroyed-vegetable coverage are identified as co-seismic landslides. Due to the thick vegetation coverage in the quake-affected area, the tone of non-seismic landslides is often dimmer than co-seismic landslides. Therefore, we mostly judge the co-seismic landslides based on whether the sliding surface is fresh or not. (Figure 6). (4) In addition, some of the detailed features can also be used, including the direction of the landslide, damaged plants on a landslide, the destroyed house and roads. These features allow us to build a co-seismic landslide inventory based on the satellite images. Figure 6b shows the distribution of co-seismic landslides near the epicenter (acquired on 11 September 2018) and images before the earthquake on 3 August 2018 (Figure 6a) in the same area. From the comparison of pre- and post-earthquake images, we can observe that a large number of co-seismic landslides were triggered by the 2018 Tomakomai earthquake (Figure 6b).

In order to keep the objectivity of the landslide number, congregate and patched landslides were divided into individual landslides expressed as multi-polygons, rather than an individual landslide as a single polygon. Finally, we established a landslide inventory of the Tomakomai earthquake, which contains at least 9295 landslides, of which 6682 pieces are larger than 1000 m² (Figure 7). These landslides are mostly distributed in mountainous areas in Atsuma, Abira-chō, and Mukawa of Hokkaido. These landslides are mostly contiguous and shallow, several meters deep-seated [3]. The study area was defined as an irregular elliptical area according to the distribution of these landslides (Figure 7), which covers 588.06 km². The total area of 9295 landslides is 30.96 km², the proportion of landslide area in the whole study area is 5.26%, and the density of landslide points is 15.80/km².

3.1.2. Influencing Factors of Landslide Susceptibility

In landslide susceptibility assessment, there is no widely accepted standard for selecting landslide controlling factors [35]. Based on our experiences and literature [17], combining the topographical, geological and seismic aspects, here we selected elevation, slope angle, aspect, curvature, distances to fault, distances to epicenter, PGA, distance to rivers, distances to roads, and lithology to mapping landslide susceptibility in the study area.

In areas of similar elevation, seismic landslides are usually equally distributed. Therefore, the elevation was selected as an influencing factor of landslide susceptibility. From a 30 m resolution DEM of the study area, the elevation was divided into seven classes: <100 m, 100~200 m, 200~300 m, 300~400 m, 400~500 m, 500~600 m, and >600 m (Figure 8a).

Another factor affecting earthquake-induced landslides is the slope features, including its angle, aspect (or facing direction), and curvature. These data were derived from DEM. Based on previous studies [36] and local terrain, the slope angle was divided into five categories: <10°, 10°~20°, 20°~30°, 30°~40°, and >40° (Figure 8b), the aspect was divided into nine categories: Flat, East, South east, North, Northeast, West, Northwest, South, and Southeast (Figure 8c), and the curvature was divided into six categories of <−2, −2~−1, −1~0, 0~1, 1~2, and >2 (Figure 8d).

Because of under-cutting action of rivers or roads may affect the stability of natural slopes, distance to rivers and roads were also selected as influence factors. Distances to rivers and roads were classified by 100m buffer interval. Then, the maps were converted into raster format, respectively (Figure 8h,i).

Lithology is a major influence factor on landslide occurrences. The lithology data was from a 1:200,000 geological map of the study area [37], which was divided into 10 categories according to the stratigraphic age groups, including Late Pleistocene lower terrace (Q3tl), Middle Pleistocene higher terrace (Q2th), Middle Pleistocene marine and non-marine sediments (Q2sr), Late Eocene to Early Oligocene marine and non-marine sediments (PG3sr), Middle Eocene marine and non-marine sediments (PG2sr), Late Miocene to Pliocene marine and non-marine sediments (N3sn), Middle to Late Miocene marine and non-marine sediments (N2sn), Early Miocene to Middle Miocene marine and non-marine sediments (N1sr), Late Cretaceous marine sedimentary rocks (K2sm), and Late Pleistocene to Holocene marine and non-marine sediments (Hsr).

Strong ground motion is an immediate trigger of seismic landslides. Because the seismic intensity interval is larger compared to PGA, or the same intensity, the upper and lower limits of PGA can be several times different due to the different magnitudes. Therefore, we selected PGA as an index to represent ground shaking. The PGA data used in this work was from the USGS (https://earthquake.usgs.gov) The range of PGA for the study area is from 0.16 to 0.82g with a 0.1g interval (Figure 8g). The area ranging from 0.5g to 0.7g accounted for 63.8% of the total area. Most landslides are concentrated in this area with PGA from 0.5g to 0.7g. Besides, the epicenter was made up of buffers with an interval of 1 km, the outer belts were truncated where they intersect with the boundary of the study area and the study area was divided into 28 categories (Figure 8f).

The seismogenic fault is strongly correlated with seismic landslides. Many studies have proved that seismogenic faults strongly controlled distributions of seismic landslides in many earthquake events [38,39] and the landslides were mostly distributed along seismogenic faults [39]. Therefore, we selected the distances to fault as an influencing factor. As the exact ground projection of the seismogenic fault is unknown yet, we referred the focal mechanism solution of this earthquake, in which the causative fault of the earthquake trends in 349° (http://www.fnet.bosai.go.jp/). Using a parallel line with a dip angle of 349° at 1 km intervals to represent the distance to the seismogenic structure, the study area was divided into 19 categories with respect to the seismogenic fault (Figure 8e).

All of the ten maps were converted into raster format with a grid cell size of 30 × 30 m. These controlling factor maps and their classes are shown in Figure 8.

3.1.3. Sampling Method for Landslide Data

For model construction, both training and testing samples are needed. Training samples are used for calculating the coefficients of landslide influence factors, and testing samples are applied to test the predictive capability of the model.

Considering the accuracy of data, we only considered the landslides with an area greater than 1000 m², thus 6682 landslides (area > 1000 m²) were selected, with a total area of 29.48 km². The median value of each landslide elevation was used to roughly distinguish the landslide source area and debris deposits. The area larger than the median elevation is the landslide source area; otherwise it is the landslide debris deposits (Figure 9). Finally, the source area of these landslides is 15.79 km².

In order to equally train the model toward positive and negative observations, we balanced the samples by randomly selecting a number of non-sliding points that are equal to the number of the real landslide observations (6682), resulting in 50% landslide points and 50% non-sliding points. Non-landslide points chosen from the outside the buffer zone of the landslide source area (buffer radius = 100 m). 6682 random points were employed as the landslide sample points which were randomly selected from the entire source area. Finally, we got 13,364 samples including 6682 sliding points and 6682 non-sliding points. Then we randomly selected 9364 samples for training modeling (including 4682 landslide points and 4682 non-sliding points), and the remaining 4000 samples were used for tests (including 2000 landslide points and 2000 non-landslide points).

3.2. Methodology

In this work, we selected LR and SVM to assess landslide susceptibility in the entire affected area of the 2018 Tomakomai earthquake. Below we will introduce these two methods. In addition, how to use the above data in these methods to construct models will be explained, representatively.

3.2.1. Logistic Regression (LR)

The Logistic Regression model (LR) is a regression analysis method, in which the dependent variable is a binary categorical one, belonging to the nonlinear multivariate statistical model [15,40]. A LR model describes the relationship between a dichotomous dependent variable P, coded to take the values 1 (presence of landslide) or 0 (absence of landslide), and j independent variables $x_1$, $x_2$, ………, $x_j$. It predicts a dependent variable on the basis of continuous or categorical explanatory variables. The selection of the logistic regression model in this study mainly is for the following three reasons: (1) The independent variables do not necessarily have normal distributions. (2) The LR model carries out a comprehensive evaluation between various control factors and landslide samples, and can better solve the problem of interdependence between factors. (3) It accepts both binary and scalar values as the independent variables, which allows for the use of independent variables that are not continuous or qualitatively derived.

The algorithm of logistic regression applies maximum likelihood estimation after transforming the dependent variable into a logit variable. It means that the LR coefficients that make the observed results most ‘‘likely’’ are selected. The relationship between the explanatory variables and the dependent variable P is nonlinear in the LR model. The relationship between the probability of landslide occurrence and the independent variables can be written as follows:

$$\mathrm{Z}=a+{\sum }_j{b}_j{x}_j \mathrm{j}=1,2,\dots \mathrm{m}$$

(1)

$$P=\frac{1}{1+{\mathrm{e}}^{-\left(\mathrm{z}\right)}}$$

(2)

where $P$ is the probability of a landslide occurrence, varying from 0 to 1 on a S-shaped curve; $Z$ is the linear combination, varying from −∞ to +∞; $a$ is the intercept of the model; $\mathrm{j}$ represents the number of independent variables; $b_j$ ($j$ = 1, 2, 3, …,m) is the slope coefficient of the model; and $x_j$ ($j$ = 1, 2, 3, …, m) represents the independent variable.

3.2.2. Support Vector Machine (SVM)

The support vector machine is a new generation of machine learning algorithms based on the nonlinear theory of covariate transformation to higher dimensional feature space and the principle of structural risk minimization, which was presented firstly by Cortes and Vapnik in 1995 [41]. It can solve nonlinear and high-dimensional pattern recognition problems with fewer samples. Generally, the characteristics of the SVM can be briefly stated as: (1) Significantly precise and robust; (2) able to model complex nonlinear decision boundaries [20]; (3) less prone to over fitting in comparison with other models such as artificial neural network (ANN) [25]; and (4) potential of implementation in pattern recognition, regression, and classification [20].

In recent years, SVM modeling approaches have been implemented extensively in landslide automatic identification and landslide susceptibility assessment [25,42]. The key to the application of the support vector machine is (1) to find an optimal classification hyperplane that is out of the data model; and (2) to use a reasonable kernel function to map linearly inseparable data into high-dimensional feature space, making it linearly separable [20].

A set of linear separable training vectors $x_i$ (i = 1, 2, …, n) was considered. The training vectors consist of two classes, which are defined as $y_i$ = ± 1. The purpose of support vector machine is to search an n-dimensional hyperplane differentiating the two classes by their maximum gap. Mathematically, it can be expressed as:

$$\frac{1}{2}{\parallel w\parallel }^2$$

(3)

The constraint is:

$$y_i\left(w•x_i+b\right)\ge 1$$

(4)

where $w$ represent the norm of the normal of the hyperplane, $b$ is a scalar base, and • is the scalar product operation, x_i is ith sample value of the input vector, and y_i is the corresponding output.

With the Lagrangian multiplier, the cost function is:

$$L=\frac{1}{2} w^2-{\displaystyle \sum }_{i=1}^n{\lambda }_i\left(y_i\left(\left(w·x_i\right)+b\right)-1\right)$$

(5)

where ${\lambda }_i$ is the Lagrangian multiplier, $w$ and $b$ can be achieved by the dual minimizing Equation (5).

For a non-separable case, one can modify the constraints Equations (6) and (7):

$$y_i\left(w·x_i+b\right)\ge 1-{\xi }_i$$

(6)

$$L=\frac{1}{2}{w}^2-\frac{1}{vn}{\displaystyle \sum }_{i=1}^n{\xi }_i$$

(7)

where ${\xi }_i$ is slack variables, and v (0, 1] represents the misclassification. In the present study, +1 and −1 indicate the failed and stable cases, respectively.

There are four commonly used kernel function types including Linear, polynomial, radial basis, and sigmoid. Research shows that the radial basis function (RBF) is one of the most powerful kernel functions [43], which has been applied successfully to remote sensing classification and landslide susceptibility mapping. Therefore, this work chose the RBF in the assessment of earthquake landslide susceptibility. The mathematical expression of the RBF $K_{RBF}\left(x_i-x_j\right)$ is

$$K_{RBF}\left(x_i-x_j\right)=exp(-\gamma {\parallel x_i-x_j\parallel }^2)$$

(8)

where γ is the parameter of the kernel functions.

Before constructing LR and SVM models, we need to prepare the training and testing sample datasets which include both the presence and absence of landslide occurrences and associated controlling factors. These values were derived and rasterized by GIS. We extracted cell values at locations specified in a point feature class (training and testing samples we selected in Section 3.1.3) from 30 × 30 raster derived from the influence factor maps, and recorded the values to the attribute table of the point feature class.

4. Results and Analyses

In this section, we will show how to use the sample dataset to perform regional susceptibility evaluation using LR and SVM, and conduct the quantitative analysis of the assessment results.

4.1. LSM of LR

For the LR model, the sample dataset includes both the presence and absence of landslide occurrences and associated controlling factors. Then the sample dataset was inserted into the SPSS software to model the relationship of the probability of landslide occurrences with the associated controlling factors [44]. In this case, the last category was used as the reference category for each categorical variable [10]. Using the LR model, the spatial relationship between landslide occurrence and controlling factors of landslides was assessed.

Table 2. Classifications and LR coefficients of ten controlling factors of landslides.

Classification	Regression Coefficients	Classification	Regression Coefficient	Classification	Regression Coefficients
<DEM>		11	4.321	7: >0.7g	0
1: <100 m	5.51	12	3.238	<Distances to roads>
2:100–200 m	6.524	13	2.892	1:<100 m	−1.005
3:200–300 m	6.364	14	2.68	2:100–200 m	−0.644
4:300–400 m	4.532	15	2.027	3:200–300 m	−0.496
5:400–500 m	−11.78	16	1.607	4:300–400 m	−0.361
6:500–600 m	0	17	1.429	5:400–500 m	−0.444
7:>600 m	0	18	0.755	6:500–600 m	−0.366
<Slope angle>		19	0	7:600–700 m	−0.231
1:<10°	−25.717	<Distances to epicenter>		8:700–800 m	−0.56
2:10–20°	−24.894	1	−27.93	9:800–900 m	−0.587
3:20–30°	−24.426	2	−8.91	10:900–1000 m	−0.173
4:30–40°	−25.102	3	−6.511	11: >1000 m	0
5:>40°	0	4	−5.07	<Distance to rivers>
<Aspect>		5	−5.185	1: <100 m	−0.028
1: Flat	−18.45	6	−4.388	2:100–200 m	−0.287
2: North	0.416	7	−4.125	3:200–300 m	0.116
3: Northeast	0.92	8	−2.974	4:300–400 m	−0.552
4: East	1.454	9	−2.353	5:400–500 m	−0.511
5: Southeast	1.423	10	−2.379	6:500–600 m	−0.614
6: South	1.021	11	−2.185	7:600–700 m	−0.386
7: Southwest	0.652	12	−2.474	8:700–800 m	−0.222
8: West	0.187	13	−2.956	9:800–900 m	−0.209
9: Northwest	0	14	−3.254	10:900–1000 m	−0.6
<Curvature>		15	−3.59	11:>1000 m	0
1: <−2	0.14	16	−3.692	<Lithology>
2: −2~−1	0.732	17	−3.897	1: Hsr	−2.517
3: −1~0	0.926	18	−4.128	2: K2sm	1.048
4:0~1	1.061	19	−5.292	3: N1sr	−3.386
5:1~2	0.834	20	−6.396	4: N2sn	0.486
6: >2	0	21	−7.523	5: N3sn	−2.831
<Distances to fault>		22	−9.391	6: PG2sr	−0.852
1	−17.689	23	−26.766	7: PG3sr	−15.327
2	0.766	24	−23.596	8: Q2sr	−0.889
3	2.891	25	0	9: Q2th	−5.309
4	4.263	<PGA>		10: Q3tl	0
5	4.229	1: <0.2g	1.369	Constant	17.594
6	3.889	2:0.2–0.3g	1.604
7	5.634	3:0.3–0.4g	1.22
8	5.249	4:0.4–0.5g	0.8
9	6.142	5:0.5–0.6g	0.798
10	5.794	6:0.6–0.7g	2.089

and Figure 10 show the extracted coefficients of all categorical variables by LR and the classified area of each categorical variable. In this paper, the last class was used as the reference class for each class variable (i.e., the coefficient of the last class is 0) [10].

Generally, the greater the regression coefficient is, the higher the significance of landslide occurrence is. Based on the regression coefficients, we can explain the statistical relationship between each control factor and the occurrence of landslides. Figure 10a shows that the highest LR coefficients are in the range of elevations 100~200 m and 200~300 m (Figure 10a). Figure 10b shows that the LR coefficients of slope angel change little within 0~40°, then reach the highest value in >40° interval. Figure 10c shows that the east and southeast aspects are easier to produce landslides. The coefficient curve of the parallel seismogenic fault shows that the LR coefficients decrease with the increase of fault distances as a whole (Figure 10e). Figure 10g shows that the areas with 0.6~0.7g are more prone to landslides. The lithology of Late Cretaceous marine sedimentary rocks (K2sm), Middle to Late Miocene marine and non-marine sediments (N2sn) have the highest LR coefficients (Figure 10j).

The corresponding weight values were given to each controlling factor according to the LR coefficients (Table 2). Based on Formula (1), the Z values and the weighted linear combination of the independent variables were calculated based on GIS platform and the superposition operations of each controlling factor layer (Figure 11a). By virtue of Formula (2) and Z and P values, the probability of landslide occurrence, was calculated subsequently. Based on the commonly used five-class interval method, the landslide susceptibility index was ranked into five classes: (1) Very low (0 to 0.2), (2) Low (0.2 to 0.4), (3) Moderate (0.4 to 0.6), (4) High (0.6 to 0.8), and (5) Very high (0.8 to 1). Figure 11b shows the susceptibility class map.

Figure 12 and

Table 3. Statistics in different landslide susceptibility classes (LR).

	Area/km2	Area of Classification/(%)	Number of Landslides	LND/km2	LAD/(%)
Very low	303.96	51.68	209	0.68	0.18
Low	66.06	11.23	283	4.28	1.19
Moderate	55.52	9.44	554	9.97	2.33
High	63.91	10.86	1307	20.44	4.43
Very high	98.60	16.76	4329	43.90	10.37

Note: LND is landslide number density; LAD is landslide area density.

are the statistical results of the susceptibility zoning area, the number of landslides, landslide area density (LAD) and the landslide number density (LND) in different susceptibility classes of two models.

As shown in Figure 12 for the LR, the area of very high susceptibility class is 98.60 km², which accounts for 16% of the entire study area; the number of landslides in this class is 4329, accounting for 64.78% of the total landslides, and the LAD and LND are 10.37% and 43.9/km², respectively. 51% of the study area is designated as the very low susceptibility classes, where the number of landslides is 209, accounting for 3.12% of the total number of landslides, and the LAD and LND are 0.18% and 0.68/km², respectively. On the whole, most landslides are concentrated in the very low and very high susceptibility classes.

4.2. LSM of SVM

Software LibSVM3.22, which is a widely used SVM software library, was employed in this study [45]. Currently, there is a lack of information in most related studies concerning the SVM learning parameter estimation of C and Gaussian kernel width γ [46]. In this study, the optimization of parameters was conducted by the cross validation algorithm of software LibSVM3.22 using matlab2012b [45].

The maps of ten influence factors with a grid cell size of 30 × 30 m were adopted to produce the landslide susceptibility maps. The training dataset included 9364 grid cells and each grid consisted of ten layers representing the landslide-influencing factors. As the prediction accuracy of the SVM depends on sample size, a model test was carried out by varying the sample size. Using the RBF kernel parameters, the best C and γ were determined as 2 and 0.5, respectively. The SVM model was constructed. Then the landslide susceptibility map (Figure 13a) for the study area was created, the closer the value is to 0, the smaller the possibility of landslide occurrence is; and the closer the value is to 1, the greater the possibility of landslide occurrence is. The landslide susceptibility index was ranked into five classes: (1) Very low (0 to 0.2), (2) Low (0.2 to 0.4), (3) Moderate (0.4 to 0.6), (4) High (0.6 to 0.8), and (5) Very high (0.8 to 1) (Figure 13b).

Figure 14 and

Table 4. Statistics of different landslide susceptibility classes (SVM).

	Area/km2	Area of Classification/(%)	Number of Landslides	LND/km2	LAD/(%)
Very low	5.92	1.00	3	0.50	0.19
Low	332.96	56.61	103	0.30	0.12
Moderate	133.98	22.78	587	4.38	1.83
High	108.31	18.41	5484	50.62	10.72
Very high	6.88	1.17	505	73.29	17.41

Note: LND is landslide number density; LAD is landslide area density.

are the statistical results of SVM, the number of landslides, landslide area density (LAD) and the landslide number density (LND). In the SVM model, the area of middle susceptibility class is 133.98 km², which accounts for 22.78% of the entire study area. The actual number of landslides is 587, accounting for 8.78% of the entire landslides, and the LAD and LND are 0.12% and 0.30 km², respectively. The area of low susceptibility class is 332.96 km², which accounts for 56.61% of the entire study area. The number of landslides is 103, accounting for 1.54% of the entire number of landslides, and the LAD and LND are 1.19% and 4.28/km², respectively.

4.3. Model Validation and Quantitative Analysis

In this study, validations of the mapping results were performed by the area under the curves of the two models. The landslide susceptibility index values were sorted in descending order and divided into 30 classes of accumulated area ratio percentages. The success rate curve and the predictive rate curve were created with the landslide training dataset and the testing dataset, respectively.

To compare the result quantitatively, the AUC (Area Under Curve) were re-calculated, in which the total area is 100 means perfect success rate and prediction rate. The results of the training dataset and the validation dataset were calculated, and the success rate curve and the predictive rate are shown in Figure 15. The results show that the two models have good prediction accuracy. For training set, the AUCs of LR and SVM are 0.847 and 0.909, respectively. For validating set, the AUC is 0.839 for LR, 0.871 for SVM. Overall, the prediction accuracy of the SVM is higher than LR.

In order to determine the relationship between the landslide susceptibility index and landslide area percentage (LAP), the range of susceptibility index was divided into 32 equal intervals, and then the classification area within each level and its LAP were counted. Overall, the LAP increases rapidly with the increase of the susceptibility level for two models. As shown in Figure 16a, the classified area based on the LR decreases rapidly with the increasing of susceptibility level, most of the study area is in the interval of <0.1 and >0.9. The results of SVM show that most of the research area is in the intervals of 0.3–0.45 and 0.65–0.8 (Figure 16b). This is why the distribution area of very high (0.8–1) and very low (0–0.2) based on LR is larger (Figure 11b), but the assessment result of SVM concentrates in the classification area of the high (0.6–0.8) and low (0.2–0.4) (Figure 13b).

5. Discussion

This work attempted to establish a detailed landslide inventory and used the LR and SVM to carry out the susceptibility mapping for the 2018 Tomakomai, Japan Mw6.6 earthquake. Using the LR and SVM, two landslide susceptibility maps were produced for the study area. The results show that the SVM slightly outperformed LR method. There are two reasons. First, the SVM method is not oriented to linear features or even features that interact linearly. It can solve nonlinear and high-dimensional pattern recognition problems better than LR [47]. Secondly, SVM is less prone to over fitting in comparison with LR. While LR is a regression analysis method, which is based on logarithmic transformation of a generalized linear model. Thus, LR does not perform well when there is more complicated relationship between the landslides and the influencing factors [47].

In the application of the LR model to landslide susceptibility mapping (LSM), many scientists have tried to exploit the data using continuous variables [14,48] or dummy binary variables. When using dummy binary variables, the complexity of the data structure will increase. It means that many parameters are included, the regression equation will be very long, and it may even introduce numerical problems [48]. Thus, some studies extended the application of logistic regression by using the continuous data as they can attain good assessment results [40,48]. Although this can decreases the complexity of the data structure and increase flexibility of the LR model, for some independent variables such as elevation, distances to epicenter and distances to faults, they may not have the direct positive or negative correlation with the prediction values (susceptibility index). If the independent variables are not converted into categorized variables, it may influence the prediction accuracy. Thus, based on previous studies [17], we used dummy binary variables for LR [17]. By classifying the independent variables we calculated the LR coefficients for each class.

Based on the resulting susceptibility index map of LR (Figure 11a) and SVM (Figure 13a), the distribution of susceptibility index are basically the same, and the high-susceptibility areas are mainly distributed along the river valleys, which is basically consistent with the actual landslide distribution. Meanwhile, the susceptibility index was divided into 32 equal intervals, and then the classification area and its LAP within each susceptibility level were counted. From the statistical results (Figure 16), the LAP curves of two models are basically consistent. The LAP increases rapidly with the increase of the susceptibility level. However, the classified area distribution of two models is dissimilar. For LR, the peak values of classified area were mainly concentrated in <0.1 and >0.9 intervals (Figure 16a). So the distribution area concentrates on the interval of very high (0.8–1) and very low (0–0.2) (Figure 11b). For SVM, by contrast, the peak values of classified area were distributed in the interval of 0.3–0.45 and 0.65–0.8 (Figure 16b). Thus, the classification area of the high (0.6–0.8) and low (0.2–0.4) is larger (Figure 13b). In this paper, the assessment results of the two models are classified by the commonly used five-class interval method [27]. As the absolute values of the susceptibility index are different, causing different class results, which can also be seen in previous studies [11,49,50]. Therefore, this reminds us that when we conduct the landslide susceptibility assessment, different classification criteria will result in different susceptibility mapping results. Thus, how to establish reasonable and objective classification standards is very important in this aspect.

In terms of optical images, the planet images have three advantages. (1) Planet Satellite can provide high frequency, medium-high-resolution (3–4 m) remote sensing images, which ensure to acquire detailed temporal and spatial information of landslides. Especially for the landslides with high density, we can follow the principles of image interpretation to prepare the landslide inventory based on Planet images. (2) High data coverage efficiency: The Planet small satellite constellation has more than 170 satellites, which can achieve global daily coverage. In the past, it was impossible to obtain images covering the entire earthquake zone in one month or even few months after the earthquake (such as the Wenchuan, Nepal, and Lushan earthquakes), which brought great difficulties to the construction of the seismic-landslide database [31,32,51]. The Planet satellite changed this status, and the images can be acquired within a few days after the earthquake, which greatly increases the acquisition speed of landslides. (3) Planet images are merit for the high precision and slight noise. This ensures the accuracy and completeness of landslides interpretation.

Otherwise, this study further proves the feasibility and reliability of the two types of machine learning methods for large-scale, high-density seismic landslide susceptibility mapping. Previous studies [11,52] have shown that machine-learning methods are more accurate than Newmark and human empirical models. However, the machine-learning method has high dependence on the comprehensive and objective landslide dataset, which makes the method have low efficient, and the assessment results are mostly years after the earthquake, which led to the fact that scientific research results could not be applied to practice [27]. But, the planet images can be acquired within a few days after the earthquake, a comprehensive and objective earthquake landslide database can be obtained quickly, then the machine learning method was used to conduct accurate landslide susceptibility mapping, which can be applied to landslide disaster prevention and mitigation in earthquake-stricken areas in the short-term after the earthquake. In this study, we completed the landslide inventory within 1 month after the earthquake, so that the assessment results can be applied to the medium-term resettlement and reconstruction, and are helpful to earthquake landslide prevention and mitigation.

6. Conclusions

A detailed and complete landslide inventory is the basis for landslide susceptibility assessment, and landslide susceptibility mapping is of great significance for disaster prevention and mitigation of earthquake-triggered landslide hazard. This work attempted to establish a detailed landslide inventory of the 2018 Tomakomai, Japan Mw6.6 earthquake and used the LR and SVM to carry out the susceptibility mapping for the affected area. Based on the planet images acquired within a few days after the earthquake, a comprehensive and objective earthquake landslide database was obtained within 1 month, so that the results can be applied to the medium-term resettlement and reconstruction.

To compare the result quantitatively, the modeling results of landslide susceptibility were evaluated by comparing the maps with known landslide locations. For LR, more than 68.44% of landslides concentrate in areas of very low and very high susceptibility classes. For SVM, more than 89.6% of landslides concentrate in areas of high and very high susceptibility classes. The AUCs show that the SVM works slightly better than the LR method, the success rates of the LR and SVM are 84.7% and 90.9%, respectively. The prediction rates of the LR and SVM are 83.9% and 87.1%, respectively. Overall, the LND and LAD of the two models increase rapidly with the increase of the susceptibility level. This study explores the applicability of the logistic regression (LR) model and support vector-based model (SVM) in the study area, and the results provide a useful reference for earthquake disaster prevention and mitigation in the quake-affected area. Otherwise, the results can provide a scientific reference for the quick establishment of co-seismic landslide inventories and landslide susceptibility mapping based on landslide databases.