Identification of Dominant Controlling Factors and Susceptibility Assessment of Coseismic Landslides Triggered by the 2022 Luding Earthquake

1. Introduction

Landslides are widespread natural disasters around the world, and earthquakes are considered one of the main triggers of landslides. The damage caused by coseismic landslides is often more severe than the direct damage caused by the earthquake itself [1]. China is situated at the junction of two major global seismic belts—the Pacific Ring of Fire and the Eurasian seismic belt. Due to the compression from the Pacific, Indian, and Philippine Sea plates, the seismic fault zones in the region are highly developed, making China a country prone to frequent earthquakes [2]. For example, the globally recognized Wenchuan earthquake (Ms 8.0) in 2008, the Yushu earthquake (Ms 7.1) in 2010, the Ya'an earthquake (Ms 7.0) in 2013, the Ludian earthquake (Ms 6.5) in 2014, and the Jiuzhaigou earthquake (Ms 7.0) in 2017 have all resulted in widespread coseismic landslides. China is geographically characterized by numerous mountains and hills, with mountainous areas covering about 70% of the country's land. This vast mountainous terrain provides abundant geographical conditions for coseismic landslides. In particular, the mountainous canyon regions of western China feature highly complex geological conditions, with densely distributed active faults and frequent strong earthquakes, leading to severe coseismic landslides. On September 5, 2022, at 12:52 PM, a magnitude 6.8 earthquake struck Luding County, Ganzi Prefecture, Sichuan Province, with the epicenter located at 102.08°E, 29.59°N and a focal depth of 16 km, reaching a maximum intensity of IX. The earthquake triggered numerous landslides, resulting in significant casualties and property damage. Therefore, conducting regional coseismic landslide susceptibility assessments is of great significance for earthquake disaster prevention, mitigation, and future urban planning in seismic zones.

Coseismic landslide susceptibility assessment refers to the evaluation of the likelihood of landslides occurring in a specific area after an earthquake. An objective, detailed, and accurate coseismic landslide inventory map forms the basis of coseismic landslide susceptibility assessments [3,4]. Several research teams have used multi-source remote sensing satellite images or drone aerial images from before and after the earthquake to interpret the coseismic landslides induced by the Luding earthquake and have constructed corresponding coseismic landslide inventories. However, due to cloud cover, these coseismic landslide inventories have significantly missed the high-altitude mountainous areas above 2500 m on the western side of the seismogenic fault [5,6,7,8]. Shao et al. used cloud-free post-earthquake remote sensing images to conduct a comprehensive interpretation of landslides in the entire Luding earthquake-affected area [9]. They provided a more complete supplement to previous landslide inventories, resulting in the most comprehensive and reliable coseismic landslide inventory for Luding to date. Therefore, We use this updated landslide inventory as the fundamental data for the coseismic landslide susceptibility assessment.

The methods for assessing coseismic landslide susceptibility can be broadly divided into three categories: engineering geological analysis, statistical regression models, and mechanics-based analysis methods. The engineering geological analysis method is a comprehensive evaluation approach based on an understanding of slope stability and engineering experience [10]. On the one hand, the historical analysis method can be used to study the developmental history of the slope and predict its future evolution. On the other hand, based on extensive field geological surveys, the structural analogy method evaluates the stability of slopes under similar geological conditions. Statistical regression models summarize and predict patterns based on the analysis of factors influencing the development of coseismic landslides. By studying the distribution patterns of past landslides, the relationships between landslides and factors like seismic activity, topography, lithology, and geological structures are explored. Susceptibility trends are extrapolated based on existing statistical rules or through semi-quantitative assessments using expert knowledge. For example, methods like the comprehensive index method [11], or quantitative assessments using statistical techniques such as cluster analysis [12], frequency ratio (FR) [13], information method [14], and weights of evidence model [15], as well as machine learning algorithms like logistic regression [16], artificial neural networks [17], and support vector machines [18], can also be employed. The third method is mechanics-based analysis, which includes pseudo-static methods [19], limited sliding displacement method [20,21,22], numerical simulations [23], and physical model testing [24,25,26,27]. These methods are based on the physical mechanisms of coseismic landslides and use mathematical or physical models to quantitatively assess slope stability under seismic conditions.

The development of coseismic landslides is influenced by various factors such as seismic activity, topography, and geological conditions, making it a multi-criteria decision-making problem. The Analytic Hierarchy Process (AHP), first proposed by Saaty, is a classic and powerful method used for multi-criteria analysis [28,29,30,31]. The AHP method has been widely applied in assessing typhoon damage [32], landslides [33,34], spontaneous coal combustion hazards [35], avalanches [36], flood hazards [37,38,39], and other natural disasters. The AHP method is also one of the effective approaches for assessing coseismic landslide susceptibility [40]. However, the traditional method has three main limitations: first, it relies on expert judgment for pairwise comparisons, introducing subjectivity and failing to quantify the weight of each factor. Second, in the absence of a disaster database, the results of the AHP method are expressed as a range of scores rather than the actual probability of disaster occurrence. Third, in seismic landslide susceptibility assessment, including more influencing factors does not always enhance the analysis. On the contrary, low-correlation factors may compromise the accuracy of the evaluation.

In response to the aforementioned issues, this paper focuses on the landslides induced by the Luding Ms 6.8 earthquake. Based on remote sensing interpretation and field validation, the landslide inventory for the 2022 Luding earthquake was refined. On this basis, the highest point of the landslide scarp (the highest elevation of the landslide's back scarp) was extracted as the primary data, and the dataset was divided into a training set (70%) and a test set (30%). Using the FR-AHP-Pearson coupling algorithm, the sensitivity of each influencing factor to coseismic landslides was analyzed. Factors with strong autocorrelation and low sensitivity were excluded, and a coseismic landslide susceptibility assessment model for the Luding earthquake was constructed, leading to the establishment of landslide susceptibility zones. Finally, 30% of non-landslide samples and 30% of landslide samples were randomly selected to form a test set, and the model's accuracy was validated using the area under the receiver operating characteristic (ROC) curve (AUC). The findings of this paper provide valuable scientific guidance for post-earthquake reconstruction and the strategic planning of future disaster prevention and mitigation in seismic zones.

2. Study Area

At 12:52 PM on September 5, 2022, a Ms 6.8 earthquake occurred in Luding County, Ganzi Prefecture, Sichuan Province. According to data from the China Earthquake Networks Center, the epicenter of the Luding earthquake was located in the Hailuogou Glacier Forest Park near Moxi Town (102.08°E, 29.59°N), with a focal depth of approximately 16 km and a maximum seismic intensity of IX. This earthquake triggered numerous landslides, resulting in 93 fatalities, 25 missing persons, over 270 injuries, and widespread damage to houses and roads [5,41]. The seismogenic fault of the earthquake was the Moxi segment of the Xianshuihe Fault, located on the southeastern margin of the Tibetan Plateau, with a left-lateral strike-slip fault mechanism [42]. In this study, the area on both sides of the Moxi segment of the Xianshuihe Fault, where coseismic landslides are densely distributed, was selected as the study area to assess landslide susceptibility. The study area covers approximately 963.77 km².

The study area, located in the Hengduan Mountains on the southeastern margin of the Tibetan Plateau, features typical alpine gorge terrain. The topography generally slopes from west to east and from north to south, with significant elevation changes. The highest elevation in the study area is 4668 m at Wanglangbao, in the northwest, while the lowest elevation is 874 m, located in the lower reaches of the Dadu River. The average elevation is 2312 m, with a relative height difference of 3736 m. The study area lies in a typical subtropical monsoon climate zone, influenced by the southeast monsoon, southwest monsoon, and cold air from the Tibetan Plateau. The annual average temperature is 15.5°C, and the average annual precipitation is 664.4 mm [8]. The largest river flowing through the study area is the Dadu River, which runs from north to south, with an overall flow direction of 171°. The river channel within the area is approximately 50 km long, with a vertical drop of about 266 m and an average slope of 5.32‰. The river is rugged and wide. The Dadu River is a primary tributary of the Minjiang River, joining it at the southern part of Leshan City. There are three major tributaries of the Dadu River within the study area: the Moxi River, Wandong River, and Tianwan River, all of which are primary tributaries on the right bank of the Dadu River.

The exposed strata in the study area, from oldest to youngest, include the Sinian, Silurian, Devonian, Carboniferous, Permian, Triassic, and Quaternary systems. According to the 2022 geological cloud 1:200,000 regional vector geological map, granite has the highest proportion in the area and is mainly distributed on the eastern side of the seismogenic fault zone. Quartzite and marble of the Permian system form the second-largest lithological component in the area, predominantly distributed on the western side of the study area. Additionally, limestone and marble of the Permian system are primarily distributed near the earthquake epicenter, while limestone and slate of the Devonian system are mainly found in the southern and central parts of the study area, south of the Xianshuihe Fault. Quaternary sediments are mainly distributed along both banks of the Dadu River and its tributaries, including Holocene alluvial layers, Pleistocene alluvial layers, gravel, and clay. Other strata in the area are less widely distributed. Due to prolonged intense tectonic activity and weathering, the rock masses are fractured, and structural surfaces are well-developed.

Since the Neogene, due to the uplift of the Tibetan Plateau, neotectonic movements have been intense. The uplift of the Tibetan Plateau is generally divided into four main stages [43]. The first stage is Himalayan stage (45 to 8 Ma BP), during which the Tibetan Plateau primarily underwent progressive northward compression and uplift, causing the overall uplift of the Sichuan-Yunnan block and resulting in eastward compression [44]. The second stage is Qiangtang stage (3.6 to 3.4 Ma BP), the Hengduan Mountains experienced significant uplift, and the plateau began to move eastward, leading to strong tectonic activity in the Sichuan-Yunnan region [45]. The third stage occurred Between 2.6 and 2.5 Ma BP, the Jinsha River-Xianshuihe Fault zone gradually transitioned from compressive-shear to left-lateral strike-slip [46,47]; and The fourth stage is Kunhuang stage (1.2 to 0.5 Ma BP), the plateau experienced accelerated uplift, shaping the current tectonic and geomorphological features [48]. Over the past 300 years, 17 earthquakes with magnitudes greater than 6.5 have occurred along the Xianshuihe Fault zone. The spatially closest earthquake to this event was the 1786 Kangding earthquake [49], while the most recent earthquake in time was the 1973 Luhuo earthquake [50]. Based on the data recorded by seismic stations within 150 km of the epicenter, a peak ground acceleration (PGA) zoning map was generated for the study area, with values ranging from 0.23 to 0.65g.

The visual interpretation of landslides was performed using a human-machine interactive approach, based on field surveys and optical satellite imagery. The optical satellite images used in this study were obtained from Planet satellite data captured in July 2022 and on September 29, 2022. Contrasting high-resolution remote sensing images before and after the earthquake, a total of 13,717 landslides were identified within the study area, with a total landslide area of 39.27 km². The largest landslide covered an area of 120,747 m², while the smallest was 16 m², with an average landslide area of approximately 3,451 m² [51]. The distribution of the interpreted coseismic landslides is illustrated in Figure 1.

3. Materials and Methodology

3.1. Sources of Data

The datasets required for the landslide susceptibility assessment of the 2022 Luding earthquake are listed in Table 1. They include: ① a remote sensing-interpreted seismic landslide inventory database, ② a 30m resolution digital elevation model (DEM), ③ a 1:200,000 scale vector geological map, ④ a satellite imagery map at a 1:200,000 scale, and ⑤ a vector PGA distribution map of the 2022 Luding earthquake.

3.2. Coseismic Landslide Susceptibility Assessment Methods

3.2.1 Analytical hierarchy process

AHP was proposed by American operations researcher Saaty in 1973. It is a decision-making method that combines qualitative and quantitative analysis to solve multi-criteria problems [31]. AHP is widely used in landslide susceptibility assessments. The basic steps of the AHP method include:

(1)

Establish a hierarchical structure model, clarifying the relationships between influencing factors. The hierarchical structure model is established by the recursive relationships among the influencing factors in the criterion layer and sub-criterion layer.

(2)

Construct the judgment matrix. The judgment matrix is the core of the AHP method. It involves pairwise comparisons of the relative importance of influencing factors in the sub-criterion layer under the criterion layer, forming a judgment matrix. In traditional AHP, the relative importance of two factors is qualitatively expressed as equally important, slightly important, moderately important, strongly important, and extremely important, and is quantified using a scale of 1, 3, 5, 7, and 9. In this paper, the modified Landslide Relative Frequency (mLRF) value is used to assign a quantitative index to each factor to construct the judgment matrix.

(3)

Conduct a consistency check. To ensure that the weight distribution obtained from the judgment matrix is reasonable, consistency checks are performed using equations (1) and (2).

$$CI={\displaystyle \frac{{\lambda }_{max}-n}{n-1}}$$

(1)

$$CR=CI/RI$$

(2)

In this equation, CI represents the consistency index, λ_max is the maximum eigenvalue of the judgment matrix, and n is the order of the matrix. If CI =0, it indicates perfect consistency in the matrix. Conversely, the larger the CI value, the worse the consistency of the matrix. Typically, a CI value below 0.1 is generally deemed acceptable in practice. When n is greater than 3, to eliminate the influence of n on CI, the average random consistency index (RI) is introduced, with values taken from Table 2. CR is the consistency ratio, and a CR value of less than 0.1 is generally required. If the CR value is below 0.1, the matrix passes the consistency test; otherwise, the matrix needs to be revised.

(4)

Calculate the overall weight of each influencing factor. After passing the consistency test, normalize the weight of each influencing factor at all levels to compute the comprehensive weight value (${\omega }_i$), and create a ranking table of the comprehensive weights of the influencing factors.

3.2.2 Frequency ratio

The FR statistical approach is a research approach that calculates the probability of landslides occurring within different classification intervals of each landslide influencing factor, thereby determining the correlation between the influencing factors and the location of landslides [52,53]. In this study, the FR method was employed to analyze the correlations between coseismic landslides and 15 influencing factors. The FR is defined by Equation (6). An FR value less than 1 indicates that the factor interval is less prone to coseismic landslides, FR > 1 indicates that the factor interval is more prone to coseismic landslides, and FR = 1 suggests that the result is indeterminate [54]. However, the single FR model overlooks the weight of each influencing factor. To address this limitation, mLRF was adopted to quantify the correlations, as defined in Equation (4).

$${FR}_{ij}={\displaystyle \frac{N_{ij}/N}{S_{ij}/S}}$$

(3)

$$mLRF=\sum {\displaystyle \frac{{(N_{ij}/N)}^2}{S_{ij}/S}}$$

(4)

In the equation: $i$ represents the influencing factor; $j$ represents the classification under the influencing factor; $F{R}_{ij}$ is the frequency ratio of the $j$th division of the $i$th influencing factor; $mLFR$ is the modified landslide relative frequency; $N_{ij}$ is the number of landslides in the $j$th division of the $i$th influencing factor; $N$ is the total number of landslides in the study area; $S_{ij}$ is the number of grid cells occupied by the $j$th division of the $i$th influencing factor; and $S$is the total number of grid cells in the study area.

3.2.3 Pearson correlation coefficient

The Pearson correlation coefficient is a statistical measure used to assess the strength of the linear relationship between two continuous variables. Given two variables, X and Y, where variable X contains n sample observations$(x_1,x_2,x_3,\cdot \cdot \cdot {,\mathrm{x}}_n)$ and variable Y contains n observations$({\mathrm{y}}_1,y_2,y_3,\cdot \cdot \cdot {,\mathrm{y}}_n)$, the Pearson correlation coefficient is defined as follows [55]:

$$\mathrm{r}={\displaystyle \frac{{\sum }_{i=1}^n\left(X_i-\overline{x}\right)\left(y_i-\overline{y}\right)}{\sqrt{{\sum }_{i=1}^n{\left(X_i-\overline{x}\right)}^2}\sqrt{{\sum }_{i=1}^n{\left(y_i-\overline{y}\right)}^2}}}$$

(5)

The value of r ranges between [-1, 1]. The larger the $r$ value, the stronger the linear relationship between variables $X$and $Y$. When $r$ = 1, $X$ and $Y$ are perfectly positively correlated; when $r$ = -1, $X$ and $Y$ are perfectly negatively correlated; and when $r$= 0, $X$ and $Y$ are uncorrelated. Generally, if the correlation coefficient $r$ > 0.5, variables $X$ and $Y$ are considered to have a high correlation.

3.3. Methodology Flow

(1) Collect and organize data such as the DEM, remote sensing images, roads, rivers, and landslide information of the study area. Use the elevation raster layer as the base map, extract the highest elevation point from the landslide surface vector data obtained by remote sensing interpretation as the landslide initiation point. This will serve as the fundamental data for subsequent analysis. The data is then split into training and testing sets in a 7:3 ratio.

(2) Extract 15 factors including elevation, slope gradient, slope aspect, plan curvature, profile curvature, surface cutting degree, topographic relief, elevation coefficient variation, lithology, distance to faults, epicentral distance, PGA, distance to rivers, FVC and distance to roads. Except for lithology, which is a discrete variable, all other factors are continuous data and are used for subsequent correlation analysis between influencing factors.

(3) Use the natural breaks method in the GIS platform to classify the 15 influencing factors, apply the FR method to obtain the FR values for each classification, and analyze the spatial distribution characteristics of coseismic landslides.

(4) Calculate the mLRF values for the 15 influencing factors, use these mLRF values as quantitative indicators for expert scoring in the AHP method to determine the weights of the 15 influencing factors, and perform a correlation test between the factors using Pearson's coefficient method, excluding factors with strong correlation but low weight values.

(5) After filtering the factors, the remaining factor FR layers and their filtrated weights are substituted into Equation (6) for weighted overlay. This yields the coseismic landslide susceptibility index for the study area, which is then classified using the natural breaks method to create the landslide susceptibility zonation map for the 2022 Luding earthquake.

$$OLSI=\sum {\omega }_i·F{R}_i$$

(6)

Note: OLSI is the susceptibility index; ${\omega }_i$ is the normalized weight value of the filtrated influencing factors; and $F{R}_i$ is the filtrated frequency ratio layer of the influencing factors.

(6) Slope units were delineated using hydrological analysis. A number of slope units equal to the size of the testing set were randomly selected. The highest elevation points of the slopes (initiation points) are extracted, and the ROC curve is used to test the rationality of the coseismic landslide susceptibility zonation results, validate the accuracy of the coseismic landslide susceptibility assessment, and analyze the main controlling factors of coseismic landslides in the study area.

The specific process of coseismic landslide susceptibility assessment is illustrated in Figure 2.

4. Results and Analysis

4.1. Selection of Coseismic Landslide Influencing Factors

When assessing susceptibility to adverse geological conditions, the selection process for influencing factors must be rigorous [56]. Earthquake-triggered landslides, as a highly destructive geological disaster, are influenced by a combination of factors, including topographical conditions, geological conditions, seismic motion, rivers, and human engineering activities. In terms of topographical conditions, elevation, slope gradient, slope aspect, plan curvature, profile curvature, surface cutting depth, topographic roughness, and elevation variation coefficient reflect surface morphology and topographic relief, which are crucial for assessing susceptibility to coseismic landslides. Regarding geological conditions, the lithology of strata significantly affects the formation and distribution of coseismic landslides. Distance to faults, epicenter distance, and PGA illustrate the correlation between coseismic landslides and their inducing factors. Additionally, distance to rivers and FVC are also important factors affecting the susceptibility to coseismic landslides. Finally, the terrain in the Luding area is steep, and highways are often built in valleys; excavation of roads may disturb the surrounding slopes. Landslides triggered by the Luding earthquake were mostly distributed on both sides of the roads, while slopes without road construction rarely experienced landslides. This indicates that human engineering activities, particularly road excavation, may play a critical role in exacerbating landslide occurrences.

Based on a comprehensive consideration of the above factors, we selected the relevant influencing factors based on existing data, including elevation, slope gradient, slope aspect, plan curvature, profile curvature, surface cutting degree, topographic relief, elevation coefficient variation, lithology, distance to faults, epicentral distance, PGA, distance to rivers, FVC and distance to roads [3,57,58]. The selection of these factors aims to comprehensively reflect the various influencing factors of coseismic landslide occurrence, providing a scientific basis for subsequent susceptibility assessments of coseismic landslides.

4.2. Analysis of the Correlation Between Influencing Factors and Coseismic Landslide Distribution

After determining the influencing factors of coseismic landslide susceptibility, to enhance the reliability of subsequent correlation testing between the factors, the factors should be converted into continuous data wherever possible, as shown in Figure 3. The lithology, a discrete variable, was classified according to its actual state, while the remaining continuous data were classified using the reclassification tool. Secondly, a total of 13,717 coseismic landslides were identified in the study area. The landslide samples were split into a training set (9,602 landslides) and a testing set (4,115 landslides) in a 7:3 ratio. Finally, using Equation (1), the frequency ratio$F{R}_{ij}$for each classification of influencing factors was calculated. This metric quantifies the occurrence frequency of coseismic landslides per unit area and visually reflects the susceptibility of landslides in different regions. Further, a statistical table of FR values for each classification of influencing factors was created (Table 3). Spatial distribution layers of coseismic landslide susceptibility factors (Figure 3) and corresponding relationship curves (Figure 4) were plotted to analyze the spatial distribution patterns of coseismic landslides. The results are as follows:

(1) Topographic Conditions: As shown in Figure 5A, within the elevation ranges of 874~1994m and 3279~4668m, the FR value decreases significantly with increasing elevation. In the elevation range of 1994~3279m, the FR value shows a slight decrease with increasing elevation, but the change is not pronounced. The susceptibility of coseismic landslides is strongly positively correlated with slope, surface dissection depth, terrain roughness, and elevation variation coefficient, with the highest FR values and slopes in the last two grades of these factors (Figures 5B, 5F~5H). Plan curvature is calculated based on slope direction, describing the terrain characteristics in the horizontal direction [59]. Within the plan curvature range of -23.60~0.17, the FR value slightly decreases as the curvature increases, while in the range of 0.17~23.57, the FR value increases with curvature. Profile curvature describes the complexity of the terrain [60,61]. Within the profile curvature range of -24.10~0.07, the FR value decreases as the curvature increases, while in the range of 0.07~23.88, the FR value increases with curvature.

(2) Geological Conditions: Lithology is a crucial influencing factor, as different lithologies exhibit varying strengths, playing a decisive role in slope stability [62]. As shown in Figure 5E, granite occupies the largest proportion of the study area, accounting for 52.77% of the total area, followed by quartzite, which covers 19.60%. Approximately 69% of coseismic landslides occur evenly distributed across granite and quartzite formations. Pyroxene peridotite and marble regions show the highest FR values, with the greatest density of coseismic landslide development.

(3) Seismic motion conditions: Within 5 km of the fault, the FR value slightly decreases with increasing distance from the fault, with FR values remaining above 1; beyond 5 km from the fault, the FR value shows a clear negative correlation with the fault distance (Figure 5J, Table 3). The Xianshuihe Fault zone, the source fault zone of the Luding earthquake, shows a concentration of large landslides as depicted in Figure 5. Figure 6K demonstrates that the FR value increases with the epicentral distance initially and then shows reduced volatility; the FR value peaks in the 4-6 km range, with a maximum value of 1.92, indicating the highest density of coseismic landslide occurrences within 4-6 km of the fault. Analysis of PGA shows that approximately 85% of coseismic landslides occur within the 0.23-0.42 g PGA range, with a significant positive correlation between FR and PGA. The FR value is highest at 0.37-0.42 g PGA, reaching 1.58, indicating the greatest density of coseismic landslides. Within the 0.42-0.65 g PGA range, the FR value decreases initially and then increases with stronger ground motion, with the minimum FR value corresponding to the 0.51-0.55 g PGA range.

(4) Rivers and FVC: Table 3 shows that the number of coseismic landslides and the FR values both exhibit a significant negative correlation with distance to rivers. Within 0-800 m of rivers, the number of coseismic landslides is the highest, accounting for 25.48% of the total landslides; in the 400-800 m range, the maximum FR value is 1.68. This indicates that river erosion at the base of slopes plays a controlling role in the development of coseismic landslides. FVC reflects the degree of vegetation cover, ranging from 0 to 1, with higher FVC values indicating denser vegetation cover. Table 3 and Figure 5N show that, overall, the number of coseismic landslides is positively correlated with FVC; within the 0-40.23% FVC range, there is a clear positive correlation between FR and FVC. When FVC exceeds 40.23%, the FR value increases slightly with FVC, but the change is minimal; within the 62.50-71.88% FVC range, the FR value is highest, indicating the greatest density of coseismic landslides.

(5) Human Engineering Activities: As shown in Figure 5O, within a distance of 0-2800 m from roads, the FR value decreases slightly with increasing distance from roads, but the change is minimal. When the distance from roads exceeds 2800 m, the FR value gradually decreases with increasing distance, and approaches zero when the distance exceeds 5200 m. The density of coseismic landslides is highest within the 0-400 m range from roads.

4.3. Determination of initial weights for influence factors

After the relationship between each influence factor and coseismic landslide susceptibility was established, the mLRF values for all 15 factors were calculated based on their corresponding FR raster layers, as described in Equation (4). The results are presented in Table 4.

Constructing the judgment matrix is a crucial step for determining the weights of influence factors. The core step in constructing the judgment matrix is to compare the influence factors pairwise to determine their relative importance. Traditional AHP methods rely on expert judgment to determine the relative importance of influence factors, which is highly subjective. This study uses the mLRF values of the influence factors to quantitatively represent their relative importance, thereby effectively reducing the subjectivity of the AHP method. The pairwise comparison of the 15 influence factors resulted in the judgment matrix shown in Figure 5.

Figure 5. Judgment matrix of influence factors (a. Elevation；b. Slope gradient；c. Slope aspect；d. Plan curvature；e. Profile curvature；f. Surface cutting degree；g. Topographic relief；h. Elevation coefficient variation； i. Lithology；j. Distance to faluts；k. Epicentral distance；l. PGA；m. Distance to rivers；n. FVC；o. Distance to roads)

Figure 5. Judgment matrix of influence factors (a. Elevation；b. Slope gradient；c. Slope aspect；d. Plan curvature；e. Profile curvature；f. Surface cutting degree；g. Topographic relief；h. Elevation coefficient variation； i. Lithology；j. Distance to faluts；k. Epicentral distance；l. PGA；m. Distance to rivers；n. FVC；o. Distance to roads)

The eigenvector of the maximum eigenvalue of the judgment matrix, after normalization, can be used as the weight vector. The maximum eigenvalue of the judgment matrix was calculated to be 15, with the corresponding eigenvector being {0.264, 0.314, 0.315, 0.240, 0.244, 0.274, 0.298, 0.332, 0.228, 0.239, 0.142, 0.278, 0.210, 0.248, 0.174}. After normalization, the weight ranking of the influence factors for coseismic landslides was obtained (Figure 6). The CR value of the judgment matrix, calculated using Equations (1) and (2), was found to be, which is well below 0.1, indicating that the weight allocation of the influence factors is reasonable.

Figure 6. Initial weight ranking of influence factors (ECV. Elevation coefficient variation；SA. Slope aspect；SG. Slope gradient；TR. Topographic relief；PGA. Peak ground acceleration；SCD. Surface cutting degree；EL. Elevation；FVC. Fractional vegetation cover；PRC. Profile curvature；PLC. Plan curvature；DTFA. Distance to faluts；LI. Lithology；DTRI. Distance to rivers；DTRO. Distance to roads；ED. Epicentral distance)

Figure 6. Initial weight ranking of influence factors (ECV. Elevation coefficient variation；SA. Slope aspect；SG. Slope gradient；TR. Topographic relief；PGA. Peak ground acceleration；SCD. Surface cutting degree；EL. Elevation；FVC. Fractional vegetation cover；PRC. Profile curvature；PLC. Plan curvature；DTFA. Distance to faluts；LI. Lithology；DTRI. Distance to rivers；DTRO. Distance to roads；ED. Epicentral distance)

4.4. Correlation analysis of influence factors and normalization of filtrated weights

In assessing coseismic landslide susceptibility, high autocorrelation among influence factors can lead to redundancy and reduce modeling accuracy [63]. Using the multi-value extraction tool on the GIS platform, raster values from the 15 layers shown in Figure 4 were extracted to the attribute table of the landslide initiation points in the testing set. This attribute table was then imported into SPSS software, where Pearson correlation analysis was performed to examine the relationships between the 15 influence factors, as shown in Figure 7. If the correlation coefficient r satisfies$\left|r\right|\le 0.5$, the correlation between the two factors is considered low [64]. The results indicate that elevation has significant correlations with elevation variation coefficient, distance to road, PGA, and distance to fault. Distance to epicenter shows strong correlation with PGA, plan curvature with profile curvature, and surface cutting depth with terrain roughness. Based on the obtained weights in Section 3.3, factors with smaller weights were removed. Consequently, elevation, distance to epicenter, planar curvature, and surface cutting depth were excluded. Finally, the weights of the remaining factors were renormalized, resulting in the final influence factor weight ranking for the Luding study area (Figure 8). Figure 8 shows that the elevation variation coefficient is the most important influence factor with a weight of 0.1152, followed by slope aspect and slope degree, both with weights around 0.11. Compared to other influence factors, distance to road is the least important factor, with a weight of 0.0605.

4.5. Mapping of coseismic landslide susceptibility zoning

After the filtrated weights of the influence factors were obtained, the FR raster layers derived in Section 3.2, together with the weights presented in Figure 8, were applied in Equation (6). A weighted overlay of the filtered FR raster layers and their corresponding weights was then performed to compute the coseismic landslide susceptibility index for the study area. Based on the values of the coseismic landslide susceptibility index, the study area is classified into five categories using the natural breaks method: non-susceptible, low susceptibility, moderate susceptibility, high susceptibility, and very high susceptibility (Figure 9).

As shown in Figure 9, the areas of very high and high landslide susceptibility are mainly distributed along both sides of the Dadu River and near the Xianshuihe Fault zone. The very high susceptibility area accounts for 12.97% of the total study area, while the high susceptibility area covers 23.84% of the study area (Figure 10). According to Figure 10, as landslide susceptibility increases, both the number of landslides and the FR value also increase. The landslide susceptibility assessment results are highly consistent with the actual distribution of coseismic landslides.

4.6. Validation of coseismic landslide susceptibility zoning results

The accuracy of the coseismic landslide susceptibility zoning results was evaluated using the Area Under the ROC Curve. The ROC curve plots the true positive rate against the false positive rate under various threshold values. AUC values ranges from 0.5 to 1.0, where values greater than 0.7 indicate good model performance and values exceeding 0.9 suggest excellent predictive accuracy [65].

In this study, 30% of the landslide initiation points (4,115 in total) were randomly selected as the testing dataset. Slope units were extracted using the hydrological analysis method in the ArcGIS (Figure 2), resulting in 44,018 units across the study area. Excluding the slopes containing all landslide initiation point from the training and testing sets, the study area's elevation raster layer was used as the base map. After excluding units containing all training or testing landslides, raster center points with the highest elevation in each slope were extracted as landslide initiation points, yielding 30,513 samples. For slope units lacking raster center point coverage, a total of 8,846 slope center points were extracted. During data preprocessing, 180 slope units near the boundary of the study area were identified as lacking valid raster center points and were thus excluded. The remaining 43,838 slope units were retained, from which 4,115 landslide initiation points were randomly selected to form the testing dataset for model validation. Using the multi-value extraction tool, susceptibility index values corresponding to the testing points were obtained. Subsequently, the ROC curve was plotted and the AUC value was calculated. As shown in Figure 11, the AUC value of 0.8445 demonstrates a high level of accuracy for the coseismic landslide susceptibility zoning model.

5. Discussion

5.1. Controlling factors of the coseismic landslides

The high and very high susceptibility zones of landslides triggered by the Luding earthquake show a strong consistency with the actual distribution of coseismic landslides interpreted through remote sensing (Figure 10 and Figure 11), validating the accuracy of the coseismic landslide susceptibility zones presented in this study. The results indicate that the elevation variation coefficient is the primary controlling factor for coseismic landslides. The elevation variation coefficient, which is the ratio of the standard deviation of elevation to the mean elevation, reflects the degree of surface erosion and incision within the region. A higher elevation variation coefficient indicates more significant surface erosion and incision [66]. There is a clear positive correlation between the susceptibility of coseismic landslides and the elevation variation coefficient (Figure 4H), meaning that regions with more pronounced surface erosion and incision typically reflect higher tectonic activity, intense river downcutting, and the development of steep slopes. Additionally, it suggests that these slopes have experienced long-term unloading deformation, with highly developed rock mass fractures, resulting in poor overall rock quality and lower strength, conditions that are conducive to the occurrence of coseismic landslides. Furthermore, Yang et al. proposed using two quantitative geomorphic parameters, slope gradient and elevation variation coefficient, to identify active faults [67]. They discovered that in the high mountain gorge region of the Yarlung Tsangpo River, east of Namcha Barwa, the elevation variation coefficient map revealed discontinuous, linear low-value bands along active fractured fault zones. A similar phenomenon was observed in this study at the location of the Xianshuihe Fault zone, reinforcing the applicability of using the elevation variation coefficient as a quantitative geomorphic parameter to identify active faults. Additionally, the high elevation variation coefficient of the slopes on both sides of the Xianshuihe Fault at this location indicates severe slope erosion, providing material sources for the development of large landslides near the fault zone. This also reflects the basic fact that the rocks along the fault's fractured zone are less resistant to erosion than the surrounding bedrock.

Slope aspect emerged as the second most important factor, likely due to three reasons. First, research on the 2008 Wenchuan earthquake found that slopes aligned with the direction of fault movement were more susceptible to landslides [68]. A similar pattern was observed in landslides triggered by the Luding earthquake, where the initial slope directions of most landslides in the western mountain gorge region were S, SE, and E, consistent with the SE strike-slip direction of the fault (Figure 3C), further supporting this observation. This reflects that the SE extension of the seismogenic fault was a primary cause of landslides, as seismic surface waves encountered free faces, where wave reflection and other effects caused tension cracks and ejection failures in the slope’s surface rock. Second, remote sensing images (Figure 1) show that in the western mountain gorge region of the study area, most coseismic landslides developed on sunward-facing slopes, with very few on shaded slopes. Snow cover was observed on shaded mountain tops, indicative of glacial geomorphology. Glaciers can play a role in moderating seismic activity and controlling landslide scale. Where glacier thickness approaches local slope height, glaciers can reduce the topographic amplification effect of seismic shaking [69]. Sunward-facing slopes, exposed to solar radiation, had melted snow cover, and glacier retreat exposed steeper, higher terrain, enhancing the topographic amplification of seismic shaking [66]. This amplification effect is most pronounced near steep slopes, mountain tops, and ridges [70,71]. Third, as shown in Figures 4N and 4I, in the western part of the study area, vegetation cover was found to be lower on shaded slopes and higher on sunward-facing slopes. Coseismic landslides primarily developed on hard rocks, such as quartzite and granite. Roots of plants can secrete organic acids that absorb minerals from rocks, altering their composition. Roots can infiltrate rock fissures, gradually expanding them. This process, particularly in already fractured rocks, can cause rock disintegration, providing ample material for coseismic landslide development. In summary, the role of slope aspect as the second most important factor is driven by a complex mechanism that integrates multiple factors influencing the development of coseismic landslides.

Slope gradient is the third most important factor, with a clear positive correlation between the susceptibility to coseismic landslides and slope gradient (Figure 4C). This finding aligns with the distribution patterns of landslides triggered by the 2008 Wenchuan earthquake, the 2010 Yushu earthquake, and the 2017 Jiuzhaigou earthquake, where slope gradient also held significant weight in the susceptibility studies of coseismic landslides [28,70,71]. Steeper slopes are more unstable, and thus, have higher susceptibility to coseismic landslides. Distance to rivers is the factor with the lowest weight. As shown in Figure 3M, in the western part of the study area, the Hailiu River and Sala Pond River, two tributaries of the Dadu River, are far from the seismogenic fault. Both the elevation variation coefficient and slope gradient are small, indicating gentle slopes and low topographic relief, which result in a sparse distribution of coseismic landslides, reducing the influence of rivers on landslides. The elevation variation coefficient, derived from elevation factors, reflects geomorphological differences. Roads and faults are located near valleys, and from the valley to the mountain peak, elevation increases with distance from roads and faults. Furthermore, seismic waves have an elevation amplification effect: the higher the elevation, the greater the PGA. Therefore, elevation, the elevation variation coefficient, distance to roads, PGA, and distance to faults are highly correlated. Seismic energy is typically concentrated near the epicenter. In the central and northern parts of the study area, PGA generally decreases with increasing distance from the epicenter. Plan curvature and profile curvature describe the terrain features in the horizontal and vertical directions. Coseismic landslides tend to develop on convex or concave surfaces in the horizontal direction and at curvature transition points in the profile (Figures 3D and 3E). Surface dissection depth and terrain relief describe the degree of surface erosion and topographic variation, both of which are factors reflecting geomorphological differences. The greater the surface dissection, the more steep and tall slopes develop. These slopes tend to have a long history of unloading deformation, with highly fractured rock masses, which promote the occurrence of coseismic landslides. Therefore, epicentral distance is strongly correlated with PGA, plan curvature, profile curvature, surface cutting degree, and topographic relief.

5.2. Influence of the seismogenic fault

The spatial distribution of coseismic landslides in Figure 1 reveals that the orientation of the Dadu River’s tributary valleys on the western side of the study area is nearly perpendicular to the Xianshuihe Fault. Furthermore, the "locked section" of the Xianshuihe Fault, where large landslides are densely concentrated, is clearly observable. This is highly consistent with the "locked section effect" proposed by Xu and Li based on the distribution characteristics of landslides from the 2008 Wenchuan earthquake [72]. Earthquakes typically release accumulated energy through strong shaking or surface deformation, with the energy most concentrated near the seismogenic fault [1]. The "locked section effect" occurs when the locked section is displaced during the earthquake and further sheared and ruptured, releasing additional energy. This results in more intense local shaking, forming secondary seismic sources and leading to the concentrated development of large landslides [72]. Therefore, the phenomenon of densely developed large landslides induced by the Luding earthquake in the locked section of the Xianshuihe Fault further corroborates that the Xianshuihe Fault is the seismogenic fault for this event. The 1999 Chi-Chi Mw 7.6 earthquake in Taiwan, the 2005 Kashmir Mw 7.6 earthquake in Pakistan, and the 2008 Wenchuan earthquake all exhibited the hanging wall effect, where more large landslides developed on the hanging wall of the thrust fault than on the footwall. The number of landslides on the western side of the Xianshuihe Fault is greater than on the eastern side [9]. This observation aligns with the fact that most aftershocks of the Luding earthquake occurred on the western side of the Xianshuihe Fault [73]. This is considered an indication of the hanging wall effect, leading to the hypothesis that the Xianshuihe Fault exhibits both strike-slip and thrust motion. To verify this hypothesis, we conducted field investigations and discovered traces of thrust-fold contact between the hanging wall of the Xianshuihe Fault and the fault zone (Figure 12). The measured strike and dip of the fault zone were 256°∠42°. After the Luding earthquake, some research teams processed InSAR data and obtained coseismic surface displacements. The results showed a maximum uplift of approximately 15 cm on the west side of the Xianshuihe Fault and a maximum subsidence of 14 cm on the east side [73,74]. Furthermore, the results of coseismic fault modeling indicated that the deformation field in the epicentral area and to the north was dominated by sinistral strike-slip motion, with a minor thrust component [75,76]. In light of the above findings, we propose that the Moxi segment of the Xianshuihe Fault may include a thrust component in its kinematics. This hypothesis requires further verification by future research teams.

5.3. Research Prospects and Disaster Mitigation Recommendations

Landslide inventories are crucial data for landslide susceptibility assessment, especially the quality of the landslide inventory [77]. The quality of a landslide inventory depends on various factors, such as scale, data collection, and the quality of remote sensing imagery [78]. However, most existing landslide inventories are incomplete, as they do not distinguish between landslide initiation zones and deposition zones [79]. In this study, we updated the landslide inventory for the Luding coseismic landslides based on remote sensing imagery and field investigations. We extracted the highest elevation points of the landslide polygons as the initiation zones of the coseismic landslides [69], using these as input for the model to improve the quality of the Luding coseismic landslide inventory.

From an administrative unit perspective, the zones with the highest susceptibility to landslides induced by the Luding earthquake are primarily concentrated in the townships of Tianwan, Caoke, Detuo, and Moxi. These townships show high landslide susceptibility, and efforts should be made to enhance InSAR-based hazard identification and dynamic monitoring of deformed bodies to prevent further exacerbation of potential landslide hazards during rainfall. In the Wandong Township area, the fractured slope rock mass, influenced by the Xianshuihe Fault Zone, has led to a high density of landslides. A large number of landslide deposits have formed a small barrier dam in the upper reaches of the Wandong River. Although the dam has breached, a substantial amount of loose material remains in the riverbed, which could easily trigger debris flows during rainfall events. This poses a significant challenge for post-earthquake reconstruction, highlighting the need for special attention to the potential geological disaster chain risks in this area. In contrast to the monotonic use of a single method for evaluating coseismic landslide susceptibility by other scholars [28,71,80], this study applied the AHP-FR-Pearson coupled algorithm to assess the landslides triggered by the Luding earthquake. The FR method reduced the subjectivity inherent in AHP, while the Pearson coefficient method minimized the correlation between influencing factors. The results indicate that the susceptibility map obtained using the AHP-FR-Pearson coupled algorithm achieved a prediction accuracy of 84.45%. The database used in this study is relatively large. In future research, we plan to apply more complex algorithms, such as random forests and convolutional neural networks.

6. Conclusions

This study applied the AHP-FR-Pearson coupled algorithm to conduct a susceptibility assessment of coseismic landslides triggered by the 2022 Luding Ms 6.8 earthquake, based on landslide initiation points data. In addition, the controlling factors affecting the spatial distribution of coseismic landslides were systematically analyzed. The main conclusions are as follows:

(1) An updated landslide inventory associated with the Luding earthquake was developed, identifying 13,717 landslides across the study area, with a cumulative area of 39.27 km².

(2) The study area was classified into five susceptibility levels, with very high and high susceptibility zones concentrated along the Dadu River and the Moxi segment of the Xianshuihe Fault. Particularly, towns such as Tianwan, Caoke, Detuo, and Moxi fall within the very high susceptibility zone, warranting focused landslide hazard investigations.

(3) The elevation variation coefficient, slope aspect, and slope gradient are the main controlling factors of coseismic landslide distribution. Coseismic landslide susceptibility is highest when the elevation variation coefficient is between 0.1 and 0.134, the slope aspect is southeast, and the slope gradient ranges from 70° to 76.63°.

(4) The coseismic landslide susceptibility model established using the FR-AHP-Pearson models achieved a prediction accuracy of 0.8445, indicating high accuracy, and can be widely applied in coseismic landslide susceptibility assessment.