3.1 Fuzzy C-means clustering (FCM)

The fuzzy C-means method belongs to soft clustering, which is widely used at present. Its core idea is to map data points of a multi-dimensional space to different clustering sets in the form of membership degree so as to determine C cluster centres in such a manner that the intercluster associations are minimised and the intracluster associations are maximised (Bezdek, 1981). For every group, each point is assigned a membership degree between 0 and 1. The membership values indicate the probability of each point belonging to the different groups (Eke et al., 2019). The steps of the FCM algorithm are as follows (Fig. 6).

1. The membership matrix μ_ij is initialised with random numbers between 0 and 1, which are used to represent the membership degree of x_i to the cluster j; it satisfies the constraint conditions

   $$\begin{array}{}\text{(1)}& \sum _{i=\mathrm{1}}^C{u}_{ij}=\mathrm{1},j=\mathrm{1},\mathrm{2},\mathrm{\dots },n,\end{array}$$

   where C represents the number of clusters.

2. For calculating clustering centres C_i, the formula is as follows (Hammah and Curran, 1998):

   $$\begin{array}{}\text{(2)}& C_i=\sum _{j=\mathrm{1}}^n{u}_{ij}^m{x}_j/\sum _{j-\mathrm{1}}^n{u}_{ij}^m,\end{array}$$

   where m controls the degree of fuzziness and m=2 is deemed to be the best for most applications (Bezdek, 1981). X_j represents the jth sample.

3. For determining the number of clustering centres, the clustering number C of the FCM algorithm is not clearly given, which is one of the key factors affecting the clustering effect. So this paper combines the non-distance-based FCM clustering effectiveness index proposed by Chen and Pi (2013) to determine the value of C. The exponent (V_cs) consists of the compactness index and separation index. The definition of compactness is as follows:

   $$\begin{array}{}\text{(3)}& C_{ij}=\left\{\begin{array}{ll}{u}_{ij}^{\mathrm{2}},& u_{ij}\ge \frac{\mathrm{1}}{c},\\ \mathrm{0},& u_{ij}<\frac{\mathrm{1}}{c},\end{array}\right.\end{array}$$

   where C_ij is the compactness of the jth sample with the ith. When u_ij is greater than or equal to 1∕c, it avoids being meaningless for being too small. When $u_{ij}<\mathrm{1}/c$, this indicates that the J sample is unlikely to belong to the ith class. When all samples clearly belong to a certain class, the compactness degree is the maximum – that is, the clustering result is compact. We define the whole compactness of sample data as follows:

   $$\begin{array}{}\text{(4)}& C=\sum _{i=\mathrm{1}}^c\sum _{j=\mathrm{1}}^n{C}_{ij}.\end{array}$$

   The definition of the separation index is as follows:

   $$\begin{array}{}\text{(5)}& S_{ij}=\min \left(u_{ik},u_{jk}\right),k=\mathrm{1},\mathrm{2},\mathrm{\dots },n,\end{array}$$

   namely, the minimum value of the membership degree of samples belonging to these two categories. When the division of the two categories is relatively clear, this indicates that the membership degree of samples belonging to a certain category must be greater than other values. Therefore, the better the clustering result is, the smaller S_ij should be, and the total separation is defined as

   $$\begin{array}{}\text{(6)}& S=\underset{i=\mathrm{1},j=\mathrm{1},i\ne j}{\overset{c}{max}}{S}_{ij}.\end{array}$$

   The smaller the dispersion is, the greater the difference between the two classes is and the better the clustering result is.

   Based on this, the clustering effectiveness V_cs index is defined as follows:

   $$\begin{array}{}\text{(7)}& V_{\mathrm{cs}}={\displaystyle \frac{C}{S}}.\end{array}$$

   In conclusion, when C is larger and the S value is smaller, V_cs is larger and the clustering effect is better.

4. Calculating the value function J,

   $$\begin{array}{}\text{(8)}& J=\sum _{j=\mathrm{1}}^N\sum _{i=\mathrm{1}}^C{u}_{ij}^m{d}^{\mathrm{2}}\left(X_j,V_i\right),\end{array}$$

   where N is the total number of observations, and j is the fuzzy objective function; d² is the Euclidean distance between the ith clustering centre and the jth data point (Wang et al., 2008).

   The operation is stopped when J is less than a certain threshold.

5. Calculating the new matrix U_ij and returning to step 2,

   $$\begin{array}{}\text{(9)}& u_{ij}={\displaystyle \frac{\mathrm{1}}{\sum _{k=\mathrm{1}}^C{\left(\frac{d_{ij}}{d_{kj}}\right)}^{\mathrm{2}/(m-\mathrm{1})}}}.\end{array}$$

Figure 6A flowchart of FCM.

Download

3.2 Factor analysis

FA is a multivariate statistical analysis method which studies the internal dependence of variables and reduces some variables with intricate relations to a few comprehensive factors (Li et al., 2016). FA is the inferred decomposition of observed data into two matrices. One matrix represents a set of underlying unobserved characteristics of the subject which give rise to the observed characteristics and the other explains the relationship between the unobserved and observed characteristics (Tolkoff et al., 2018). The mathematical formula can be expressed as follows:

where X(x₁, x₂, …, x_p) is the original factor; F(F₁, F₂ …, F_m) is the common factor; A=(akj), p×m, is the factor-loading matrix; akj represents the load of the k original factor on the J common factor; and ε=(ε₁, ε₂, …, ε_p) is a special factor.

The main calculation steps of the factor analysis method can be divided into six steps (Fig. 7).

1. Test the feasibility of FA of original evaluation index variables. In this paper, SPSS was used to provide a Bartlett sphericity test to determine whether variables are suitable for FA.

2. Standardised calculation of original data. In order to eliminate the numerical differences of different variables in the order of magnitude and dimension, the original data should be standardised. This paper adopted the Z standardisation method in SPSS software.

3. Construct a common factor F. In the study, the first m factors for which the cumulative variance contribution rate is no less than 85 % were selected as common factors to represent the original data.

4. Factor rotation. In this paper, varimax orthogonal rotation was used to realise factor rotation.

5. Calculating factor scores. The most common method for calculating factor scores is the Thomson regression method (Tolkoff et al., 2018), and the formula is as follows:

   $$\begin{array}{}\text{(11)}& F=A^{\prime }{R}^{-\mathrm{1}}X,\end{array}$$

   where $A^{\prime }{R}^{-\mathrm{1}}$ is factor-scoring coefficient matrix and A is the factor-loading matrix after rotation.

6. Calculating weight. The product of the factor-scoring coefficient and variance contribution rate is the contribution of each factor in the sample, and the sum of the contribution of each factor divided by the contribution of all indices is the weight of each factor. It is expressed by the formula

   $$\begin{array}{}\text{(12)}& {\mathit{\omega }}_i={\displaystyle \frac{\sum _{j=\mathrm{1}}^m{\mathit{\beta }}_{ji}{e}_j}{\sum _{i=\mathrm{1}}^p\sum _{j=\mathrm{1}}^m{\mathit{\beta }}_{ji}{e}_j}},\end{array}$$

   where β_ji is the coefficient score of each index in principal component F_j;  i=1, 2, …, p; j=1, 2, …, m; and e is the contribution rate of factor variance.

Figure 7A flowchart of FA.

Download

3.3 Combination weighting method

Considering the defects of the current method for determining the weight of factors, the combination of a analytic hierarchy process and factor analysis method is used to determine the weight of each influencing factor of debris flow.

Table 1The random average consistency index.

Download Print Version | Download XLSX

Table 2Definition of comparative importance.

Download Print Version | Download XLSX

3.3.1 Analytic hierarchy process (AHP)

The AHP was first proposed by Saaty (1978), a famous American mathematician. It decomposes the factors related to decision-making into multiple layers, such as the target layer, criterion layer and scheme layer. The AHP is a subjective weighting method and has obvious advantages in determining the weight of each factor. The specific steps are as follows.

1. Establishing a hierarchical structure model. The hierarchical structure is mainly divided into three layers: the target layer, criterion layer and scheme layer.

2. Establishing the judgement matrix. For the same level, the judgement matrix is established by pairwise comparison. The formula is as follows:

   $$\begin{array}{}\text{(13)}& A={\left(a_{ij}\right)}_{n\times n},a_{ij}>\mathrm{0},a_{ij}={\displaystyle \frac{\mathrm{1}}{a_{ji}}},(i,j=\mathrm{1},\mathrm{2},\mathrm{\dots },n),\end{array}$$

   where a_ij is the ratio of relative importance between element B_i and B_j, which is usually expressed by the scoring method from 1 to 9 (Saaty, 1978), as shown in Table 2.

3. Consistency testing. The consistency test is divided into three steps. Calculating the consistency index (CI) (Saaty, 1977a, b), the expression is

   $$\begin{array}{}\text{(14)}& \mathrm{CI}={\displaystyle \frac{{\mathit{\lambda }}_{\max }-n}{n-\mathrm{1}}},\end{array}$$

   where λ_max is the largest eigenvalue of the judgement matrix A.

4. Average random consistency RI. RI is associated with the order of the judgement matrix, and their relationship is shown in Table 1.

5. Obtaining the test coefficient CR. This can be calculated by the following equation:

   $$\begin{array}{}\text{(15)}& \mathrm{CR}={\displaystyle \frac{\mathrm{CI}}{\mathrm{RI}}}.\end{array}$$

   If CR < 0.1, the judgement matrix has a good consistency with reasonable judgement. Otherwise, the judgement matrix needs to be revised until the consistency test is satisfied.

Table 3Factors frequently used in susceptibility analysis of debris flow.

Download Print Version | Download XLSX

3.4 Efficiency coefficient method

Based on the principle of multi-objective programming, the efficiency coefficient method transforms each factor into a measurable evaluation score through the efficiency function and combines the weight of factors to make a comprehensive evaluation. The specific steps are as follows.

1. Select the evaluation factors.

2. Determine the satisfactory value and the unallowable value: the satisfactory value is a value based on years of experience, while the unallowable value is the lowest or highest acceptable value of the evaluation index.

3. Calculate the single efficacy coefficient. The single efficacy coefficient was calculated by the corresponding efficacy function based on the sensitivity of each factor. It is mainly divided into three variables: the extremely large variable (the higher the factor, the higher the efficiency coefficient), the infinitesimal variable (the smaller the index value, the larger the efficiency coefficient value) and the interval variable (the value reaches its highest in a certain interval). The specific formula is as follows:

   $$\begin{array}{}\text{(19)}& g_{\mathrm{1}i}=\left(\begin{array}{ll}\frac{x_i-x_{ni}}{x_{yi}-x_{ni}}\times \mathrm{40}+\mathrm{60},& x_i<x_{yi},\\ \mathrm{100},& x_i\ge x_{yi},\end{array}\right.\end{array}$$

   where g_1i is the single efficacy coefficient value of the ith extremely large factor, X_i is the actual value of the ith factor, X_yi is the satisfactory value of the ith factor and X_ni is the unallowable value of the ith factor.

   The infinitesimal variable is calculated as follows:

   $$\begin{array}{}\text{(20)}& g_{\mathrm{2}i}=\left(\begin{array}{ll}\frac{x_i-x_{ni}}{x_{yi}-x_{ni}}\times \mathrm{40}+\mathrm{60},& x_i>x_{yi},\\ \mathrm{100},& x_i\ge x_{yi}.\end{array}\right.\end{array}$$

   The interval variable is calculated as follows:

   $$\begin{array}{}\text{(21)}& \begin{array}{rl}& g_{\mathrm{3}i}=\\ & \left\{\begin{array}{ll}\left(\mathrm{1}-\frac{x_{\min }-x_i}{x_{\min }-x_{n\min }}\right)\times \mathrm{40}+\mathrm{60},& x_i<x_{\min },\\ \mathrm{100},& x_{\min }<x_i<x_{\max },\\ \left(\mathrm{1}-\frac{x_i-x_{\max }}{x_{n\max }-x_{\max }}\right)\times \mathrm{40}+\mathrm{60},& x_i>x_{\min }.\end{array}\right.\end{array}\end{array}$$

4. Calculating the total efficiency coefficient,

   $$\begin{array}{}\text{(22)}& G=\sum _i^m\left(g_i{\mathit{\omega }}_i\right),\end{array}$$

   where G is the total efficacy coefficient, g_i is the single efficacy coefficient and ω_i is the weight of the ith factor.

   The flowchart for the method used for our classification and susceptibility analysis is shown in Fig. 8.

Figure 8Flowchart used for classification and susceptibility assessment.

Download

Table 4The values for the 13 factors of the 21 debris flow catchments.

Download Print Version | Download XLSX

3.5 Influencing factors

The topographical, geological and climatic factors play a critical role in the distribution and activities of debris flows (Di et al., 2008). Table 3 shows the influencing factors selected by research in debris flow susceptibility assessment in recent years. Rainfall is one of the most pivotal external factors inducing debris flow disasters, but the meteorological data in our area are all from the same station, which cannot reflect the differences between each catchment. Therefore, rainfall was not included in this study. In addition, the frequency of debris flow and the size of soil particles are difficult to obtain accurately. The loose-material volume reflects the lithological characteristics and fault length to some extent, so lithology and fault length were not taken into account. The basin area, main channel length, drainage density, average slope angle, average gradient of the main channel, vegetation coverage, maximum elevation difference and curvature of the main channel, which were cited and available, were selected in this paper. As source conditions, the loose-material volume and the loose-material supply length ratio were also considered. As the study area is located in a tourist area with a relatively dense population, population density is selected as the factor of human activities. A total of 13 influencing factors were selected based on the previous research findings to reflect the characteristics of the watershed. All these factors were acquired in our field survey or calculated in ArcGIS, as described below.

• Basin area (F1) (km²). The basin area reflects the scale of debris flow. Generally, the larger the basin area is, the greater the risk of debris flow will be. It was obtained by geometric operations in ArcGIS and corrected by the remote-sensing images in Google Earth.

• Main channel length (F2) (km). Main channel length reflects the potential for increasing loose sources along the route. This value was measured from ArcGIS by combining RS technology and topographic maps.

• Drainage density (F3) (km km⁻²). Drainage density is the ratio of the total drainage length to the watershed area, and it is an important index for describing the development of gullies in the watershed.

• Average gradient of the main channel (F4). This is the ratio of the maximum elevation difference of the main channel to its linear length. The larger the value, the better the hydrodynamic condition is. This value is obtained from the digital elevation model (DEM).

• Average slope angle (F5) (^∘). As F5 increases, the erosion capacity and intensity of precipitation increase. The value was obtained by the ArcGIS slope analysis tool.

• Maximum elevation difference (F6) (m). The difference between the maximum and minimum elevation values in the basin provides the kinetic energy condition of disaster. This value is also obtained from the DEM.

• Curvature of the main channel (F7). F7 is the ratio of the main channel length to its linear length, which reflects the degree of channel blockage.

• The loose-material volume (F8) (×10⁴ m³). The loose material is one of fundamental factors triggering debris flows. This factor is obtained through field investigation with tape and a laser rangefinder. The thickness was obtained by field estimation and a trench test.

• The loose-material supply length ratio (F9). F9 is the ratio of loose-material length along a channel to total channel length, which reflects the successive supplied sediments. It was obtained through field surveys and RS technology.

• Vegetation coverage (F10). The lower the vegetation coverage, the more serious the soil erosion. It was estimated from field surveys and SPOT5 imaging.

• Population density (F11) (number of people per km²). With the development of social economy, human activities have gradually become an important factor affecting debris flow. Population density reflects the intensity of human activities, which is estimated according to the number of buildings through field survey and RS technology.

• Roundness (F12). Roundness is the morphological statistical element of a gully, and the plane shape of a gully differs from its developmental stage. F12 is the ratio of the length of the main channel of debris flow to its area.

• The highest volume of one flow (F13) (×10⁴ m³). Liu et al. (1993) selected F13 as the main factor in the risk assessment of debris flow, which is one of the most important factors in evaluating the degree of debris flow hazard.

4.1 Fuzzy C-means clustering analysis

The curve of the clustering effectiveness index V_cs with the number of clustering centres is shown in Fig. 9, and the optimal number of clustering centres is four. Based on the basic data of 21 debris flow catchments, the FCM was carried out and set the fuzzy weighted index to m=2. The results are shown in Table 5.

Figure 9Clustering validity function V_cs.

Download

Thus 21 debris flow catchments in the study area are divided into four categories. The data of each catchment belonging to the same category have a certain internal similarity and vary greatly among different categories. In other words, data of different influencing factors have different effects on different types of debris flows, which provide a favourable basis for us to analyse the main influencing factors of debris flows, and also point out the direction for monitoring and prevention of debris flows.

4.2 Factor analysis

Based on the clustering results of 21 debris flow catchments, FA was used to analyse each type of debris flow. Tables 6–9 are the results of the first, second, third and fourth categories, respectively.

Table 5Clustering results of 21 debris flow catchments.

Download Print Version | Download XLSX

As shown in Table 2, in the first category, the accumulative contribution rate of the first three factors (C1–C3) reaches 86.40 %, which retains most information of the 13 original variables. For the first group, the load values of the main factors 1–3 are relatively large in the basin area, the highest volume of one flow, the maximum elevation difference, the main channel length and curvature of the main channel, and population density and drainage density. Similarly, in the second type, the load values of the main factors 1–3 are relatively large in the basin area, the main channel length and population density, loose-material volume and drainage density, and maximum elevation difference. In the third category, the load values of the main factors 1–3 are relatively large in the basin area, main channel length, the highest volume of one flow, loose-material volume, and the loose-material supply length ratio and vegetation coverage. In the fourth category, the load values of the main factors 1–3 are relatively large in main channel length, drainage density, loose-material volume, the highest volume of one flow, and the loose-material supply length ratio and population density.

Table 6The factor-loading matrix after rotation and contribution ratios for the first category.

Download Print Version | Download XLSX

Among the 13 factors, the basin area and the highest volume of one flow reflect the scale of debris flow eruption. The main channel length, drainage density, average gradient of the main channel, the average slope, maximum elevation difference, curvature of the main channel and roundness reflect the topographical condition. The loose-material volume and the loose-material supply length ratio are the material sources for debris flow. Vegetation coverage reflects the geomorphologic condition. Population density reflects the impact of human activities on nature to some extent. Therefore, four types of debris flows can be named according to the results of the FCM and FA.

Table 7The factor-loading matrix after rotation and contribution ratios for the second category.

Download Print Version | Download XLSX

The first category can be defined as debris flow closely related to scale–topography–human activities. Considering the situation, monitoring and control of basic material sources is recommended. Similarly, the second, third and fourth categories can be defined as topography–human activities–provenance, scale–provenance–topography and topography–scale–provenance–human activities, respectively. In the same way, corresponding prevention measures can be proposed according to the characteristics of each type of debris flow.

Table 8The factor-loading matrix after rotation and contribution ratios for the third category.

Download Print Version | Download XLSX

4.3 Weights of major factors

Based on the FA of each category of debris flow in the previous section, the main influencing factors were obtained. However, the repeatability of evaluation information should be reduced. Average slope angle and average gradient of the main channel are both indicators of potential energy, so the average gradient of the main channel is omitted. Similarly, curvature of the main channel, the loose-material supply length ratio and roundness were omitted. Thus nine factors, including basin area (F1), main channel length (F2), drainage density (F3), average slope angle (F5), maximum elevation difference (F6), the loose-material volume (F8), vegetation coverage (F10), population density (F11) and the highest volume of one flow (F13) were selected. On the other hand, a reduction in the number of indicators facilitates the allocation of weight values.

Table 9The factor-loading matrix after rotation and contribution ratios for the fourth category.

Download Print Version | Download XLSX

Figure 10Hierarchical structure for debris flow susceptibility analysis.

Download

4.3.1 Subjective weights

The AHP was applied to calculate the subjective weight in this paper. The hierarchical structure (Fig. 10) was constructed, and the 1–9 scale method was used to grade each factor. The judgement matrices $\mathbf{A}-{\mathbf{A}}^{\prime }$ (Table 10) and $\mathbf{B}-{\mathbf{B}}^{\prime }$ (Table 11) were constructed, and the consistency test was conducted. The weight values of each factor are shown in Table 12.

Table 10Comparison matrix elements for geology condition.

CR = 0.0045<0.1 met the conformance inspection requirements.

Download Print Version | Download XLSX

Table 11Comparison matrix elements of the criterion level factors.

CR = 0.02<0.1 met the conformance inspection requirements.

Download Print Version | Download XLSX

Table 12The weighted values of the factors obtained by AHP.

Download Print Version | Download XLSX

4.3.2 Objective weights

FA was applied to calculate the objective weight in this paper. The weight values of each factor are shown in Table 13.

Figure 11Tree diagram obtained by hierarchical clustering.

Download

Table 13The weighted values of the factors obtained by factor analysis.

Download Print Version | Download XLSX

Table 14The combined weighted values of the factors.

Download Print Version | Download XLSX

4.4 The efficacy coefficient of factors

Among the nine factors, basin area, main channel length, drainage density, maximum elevation difference, the loose-material volume, the highest volume of one flow and population density are all extremely large variables. Vegetation coverage is the infinitesimal variable. The average slope angle is an interval variable. Table 15 shows the efficacy coefficient scores of 21 debris flow catchments after being combined with the weight calculation.

Table 15The efficacy coefficient scores of 21 debris flow catchments.

Download Print Version | Download XLSX

4.5 Susceptibility assessment of debris flow

Taking the total efficiency coefficient of each catchment as the evaluation standard (the larger the value, the higher the possibility of debris flow), the FCM was conducted for 21 debris flow catchments in the study area. The result showed that the susceptibility of debris flow was divided into three grades: high (H), moderate (m) and low (L). Combined with the classification of each debris flow mentioned above, the final results are shown in Table 16.

Table 16The qualitative description and susceptibility class for each debris flow catchment.

Download Print Version | Download XLSX

As shown in Table 16, susceptibility for the 10th and 13th catchments was high, and both of them belong to the debris flow with a close relationship between topography, human activities and provenance. Susceptibility for six catchments, including the 1st, 4th, 6th, 17th, 20th and 21st, was medium. The other 13 had low susceptibility.

Normative scoring, the k-means clustering algorithm and hierarchical clustering were determined to validate susceptibility analysis methods used in this paper.

Based on the field investigation, the 10th catchment is located in Huangsongyu National Mining Park, where a large amount of slag accumulated. With low vegetation coverage and steep terrain, the gully was in its prime, which directly threatened the safety of villagers and tourists. What is more, there are several warning boards for natural disasters and corresponding monitoring equipment in the scenic spot (as shown in Fig. 5). The 13th catchment is located in the Lishugou village scenic spot. Part of the pedestrian passageway was built, but a lot of stones were piled up in the trench and the road was broken and steep (as shown in Fig. 6). However, there is no obvious accumulation of loose materials in the catchments with low susceptibility. The gully was in its old stage, with high vegetation coverage and little human interference. The quantitative comprehensive evaluation results of debris flow susceptibility are shown in Table 17, which are divided into two levels: low (L) and moderate (M). Among them, the susceptibility levels of the 10th and the 13th catchments were moderate and the others were low.

Table 17Comparison of susceptibility analyses based on different algorithms.

Download Print Version | Download XLSX

The k-means clustering algorithm (k) (Hartigan and Wong, 1978) and hierarchical clustering (H) (Kimes et al., 2017) were used for the classification of our data to measure the classification performance in this paper. The results are shown in Table 17. The susceptibility results obtained by k and the FCM are exactly the same. The susceptibility assessment of 17th and 21st were high when based on H and moderate from the FCM and k. However, such minor differences are acceptable. On the other hand, the susceptibility results obtained by the FCM and normative scoring are different. This is mainly because the number of categories is different, and the level was generally higher when obtained by the FCM. In addition, it can be seen from the tree graph (Fig. 11) obtained by hierarchical clustering that the clustering results are more reasonable when divided into three categories, which is consistent with V_cs. Therefore, the susceptibility model established in this paper is suitable and reasonable.