Risk assessment of rain-induced debris flow in the lower reaches of Yajiang River based on GIS and CF coupling models

2.1 Research scope

The study area is located in the east of Nyingchi City, Tibet Autonomous Region (Figure 3a). The overall elevation trend of Nyingchi City is higher in the north and lower in the west, resulting in higher terrain in the upper half of Yajiang and lower in the lower half of Yajiang (Figure 3b). Nyingchi area is a high-frequency area of debris flow due to its unique hydroclimatic conditions, plateau environment, complex stratigraphic lithology conditions, and vegetation coverage. The Yajiang River flows from Lang County to Nyingchi City and from Metuo County to India. The overall distribution is high in the east and low in the west, with a large elevation span of 5,110–86 m (Figure 3c). It passes through the Brahmaputra Canyon which is the largest and deepest canyon in the world. Steep cliffs on both sides of the canyon, narrow riverbeds, and rapid flow of water, for the occurrence of debris flow, nurtured good conditions.

Figure 3

Geography and disaster site information of the study area.

With the development and utilization of the water resources of the Yajiang River, it is found that natural geological disasters are easy to occur near the basin of the Yajiang River. Once the disasters occur, they will cause serious threats and losses to engineering and construction. In this paper, the width of the two banks of the Yajiang River in southeast Tibet is 2 km in the study area. The study area is a total area of about 4309.5 km², mud-rock flow point 42, roughly the information shown in Table 1; 42 non-debris flow points were randomly selected in the study area, and 84 debris flow research sample banks were obtained to carry out assessment work.

2.2 Disaster-inducing factors

2.3 Data collection

The debris flow points were provided by geological disaster survey data of counties and cities in the Tibet Autonomous Region. In this article, all debris flow events that occurred in the study area from 2001 to 2009 were selected as research objects. Through the random point creation function in Arcgis, 42 non-debris flow points were selected (there was no debris flows near these points for more than a decade).

Geological and geomorphic information such as elevation, slope aspect, and relief will not change greatly in decades, so there is little difference in the results caused by choosing different times. In this article, the data from 2020 are used to predict, and the data from 2000 that did not occur debris flow are used to study.

2.4 Correlation of factors

Based on the data of 42 debris flow points, the correlation analysis of 11 important index data was carried out using SPSS software (Table 2).

There is a weak positive correlation between elevation and relief, elevation and NDVI, and relief and SPI. The correlation coefficients between the other factors are small, which can be considered as having no correlation. The factors used met the factor selection conditions.

3.1 CF model

The CF model can well calculate the CF values of different categories of the same factor, and the calculation formula is as follows:

where PPa is the occurrence probability of debris flow points in different categories of single factor area, and PPs is the occurrence probability of debris flow points in the total study area. CF values range from −1 to 1, with the closer to 1 the more certain, and the closer to −1 the more uncertain.

The 10 disaster-causing factors of 42 debris flow sample points in the study area were classified. The CF value is calculated by the CF model, as shown in Table 3.

It can be seen from Table 3 that the fault distance is between 12,000 and 14,000 m; The elevation is between 2,900 and 3,100 m; NDVI ranges 0.4–0.6; annual rainfall in the 900–1,000 mm range; the profile curvature is between 0 and 5; the fluctuation is in the range of 0–10 m; silty soil content in the range of 30–40; TWI ranges from 10 to 12; SPI is in the range of −4 to 0; when the slope direction is northwest, debris flow is most likely to occur.

3.2 CF–AHP

AHP is an evaluation method combining qualitative and quantitative methods, which decomposes the elements related to decision-making into levels of objectives, criteria, schemes, etc. [24]. The quantization table of AHP is shown in Table 4. According to the quantization values of 1–9, all disaster-causing factors are compared in pairs and assigned, and the judgment matrix is listed. Scholar Wu Chenshuang, who has a certain in-depth study on the whole geological disaster of Tibet, was invited to score this time [6], and the opinions of Academician Cui Peng, Professor Tang Chuan, and other relevant personnel with great research achievements in this field were also referred to the study by Liu et al and Guo et al. [25,26]. The score is relatively reliable.

By multiplying the weight of each factor calculated by AHP by the CF value of the corresponding 11 factors of a grid, the debris flow occurrence weight w of the grid can be obtained. W is a positive number that represents the mudslide, w is a negative representative mudslide that does not occur, and the size of the |w| represents the influence degree. Finally, 84 debris flow points in the sample database were used to verify the accuracy of the model.

In order to verify the rationality of the scoring of quantified values, all evaluation matrices should be subjected to a consistency test, which can prove that the expert scoring is reasonable. Equation (5) is obtained by converting Aw = λw, which is used to calculate λ _max. According to formula (6), CI can be obtained. The smaller the value of CI, the more reliable the weight. According to formula (7), CR can be calculated. When CR <0.1, the matrix passes the consistency test.

where A is the constructed judgment matrix, w is the weight, n is the number of indicators, and RI can be obtained by referring to Table 4.

The disaster factors are divided into two layers. In the first criterion layer, A1 is the terrain, A2 is the solid source, and A3 is the water source. The matrix is shown in Table 5, CR = 0.058 < 0.1, passing the consistency test. The second criterion layer is divided into three matrices, as shown in Table 6, which are B1 elevation, B2 slope aspect, B3 profile curvature, and B4 relief; C1 fault distance, C2 silt content, C3 NDVI, C4 land utilization; D1 rainfall, D2 TWI, D3 SPI. CR was 0.033, 0.0582, and 0.004, respectively.

3.3 CF–AHP

The binary LR model is an algorithm that expresses or predicts trends with a linear relationship based on statistics. 70% of the sample points, including 30 groups of debris flow points and 30 groups of non-debris flow points, were randomly selected to construct the model. The remaining 30% of debris flow was used to validate the model. The calculated CF value was used to replace the data of 11 disaster-causing factors of debris flow and was used as the independent variables of the model. Whether debris flow occurs (“1” means debris flow occurs, “0” means no debris flow occurs) is used as the dependent variable of the model [27,28], and the results are shown in Table 7.

Therefore, the occurrence probability of debris flow in the cell can be expressed as follows:

3.4 CF–RF

RF is a non-parametric statistical technique based on regression or classification of a decision tree set (forest) [29]. The workflow flow chart is shown in Figure 4. The importance of features is obtained by using 70% debris flow sample points. The remaining 30% verifies model accuracy.

Figure 4

Schematic diagram of random forest.

4.1 Analysis of calculation results

In the AHP, the weight of the first criterion layer and the second criterion layer are multiplied to obtain the weight value of the disaster factor. The variable coefficients and constants of disaster factors were obtained by binary LR. The characteristic importance of disaster factors is obtained by random forest after data training (Table 8).

It can be seen from Figure 5 that elevation, silty soil, profile curvature, and SPI in the CF–LR model have a large weight. In the CF–AHP model, the weight of rainfall, land use, and TWI are larger. In the CF–RF model, elevation, slope aspect, and relief have greater weight.

Figure 5

Schematic diagram of weight.

Although the maximum weight factors of the three models are different, they all show that topography has an important impact on the development of debris flow. The difference in weight factors may be due to the lack of local attention to the debris flow database in the lower reaches of the Yajiang River, resulting in insufficient debris flow sample points and insufficient research data. The lack of sufficient sample data will affect the accuracy of the LR model and the random play of RF. At the same time, AHP is highly subjective, which exaggerates the influence of rainfall on debris flow when assigning factors.

4.2 Comparison of model accuracy

Confusion matrix is a standard way to evaluate the accuracy of a model. Its function is to judge the probability of correct judgment and the probability of wrong judgment. Its main indicators are precision, Recall, specificity, and F1 Score (the closer to 1, the better the effect). It can be seen from Table 9 that the specificity of the RF model is slightly higher than that of the LR model, and other indicators of the LR model are better than other models. Therefore, the LR model is applicable as the research method in the research area at the present stage.

Receiver operating characteristics (ROC) graphs are useful for organizing classifiers and visualizing their performance. The ROC curve of the three methods was verified by the established debris flow database (Figure 6). According to the area under the curve (AUC) (Figure 7), the CF–LR coupling model has the highest accuracy, followed by the CF–RF coupling model, and finally CF–AHP coupling model.

Figure 6

ROC curve validation of the three models.

Figure 7

AUC values for different models.

4.3 Susceptibility map of study area

GIS can well collect, analyze, and process complex spatiotemporal data and provide good help in processing geological disaster information. The calculated CF value is used to assign values to the raster graph, as shown in Figure 8. The CF–LR coupling model is used to make the debris flow disaster susceptibility distribution map. The evaluation chart of debris flow susceptibility in the lower reaches of Yajiang River is shown in Figure 9.

Figure 8

CF value raster layer. (a) Fault distance, (b) elevation, (c) NDVI, (d) average rainfall, (e) aspect, (f) profile curvature, (g) relief, (h) silt content, (i) TWI, (j) SPI, and (k) land utilization.

Figure 9

Susceptibility map of rainfall-type debris flow in the study area.

As can be seen from Figure 10, the probability of most grid points is close to 0 or 1, with a roughly uniform distribution in the middle. Therefore, the equidistant classification method was used to divide the susceptibility into four categories: low susceptibility area, mild susceptibility area, moderate susceptibility area, and high susceptibility area. It can be seen from Figure 7 that most of the study area belongs to the low-risk area of debris flow, accounting for 54.57% of the total study area. mild, moderate, and high-risk areas accounted for 6.1, 5.61, and 33.76%, respectively. With the Grand Canyon as the boundary point, the areas prone to debris flow are basically distributed on both sides of the upstream river bank of the Grand Canyon and a small amount of them are distributed downstream.

Figure 10

Grid occurrence probability distribution map.

Through satellite images (https://www.earthol.com/) and field investigation, several places with potential debris flow hazards in high-prone areas were found, mainly in the “n” shaped watershed. In addition to the steep terrain, these areas also have fast-flowing water and complicated flow direction, and the corrosion and scouring of the water greatly increased the potential debris flow hazards (Figures 11–13).

Figure 11

“n” Shape watershed topography (93.1°E, 29.1°N).

Figure 12

The first multi-channel bend terrain (93.5°E, 29.16°N).

Figure 13

The second multi-channel bend terrain (94.5°E, 29.45°N).

4.4 Study of different rainfall

By introducing CF values of different annual rainfall into the model, the influence of annual rainfall on the study area is studied. The CF value of grade A rainfall (600–700 mm) is −0.7581, that of grade B rainfall (700–800 mm) is 0.2136, that of grade C rainfall (800–900 mm) is 0.2582, and that of grade D rainfall (900–1,000 mm) is 0.4092. The susceptibility of debris flow at different gears in the study area was obtained through grid calculation and processing, as shown in Figure 8 (Figure 14).

Figure 14

Comparison of vulnerable areas under four grades of annual rainfall.

The results show that when the annual rainfall is between 600 and 1,000 mm, the proportion of low-prone areas increases with the increase in annual rainfall, while the proportion of high-prone areas decreases with the increase in annual rainfall. The number of debris flows is not much when the rainfall is 600–700 mm, but suddenly increases when the rainfall of B-grade is 700–800 mm. This may be due to the influence of the geographical environment and geological factors in the study area. When the annual rainfall reaches 600–700 mm, it is at the critical point of debris flow outbreak, and when the annual rainfall exceeds the critical point, debris flow will be triggered. Debris flow takes away effective solid material sources after the eruption, making it difficult to repeat debris flow. However, due to the concave terrain in some areas, drainage is difficult, so more water sources have accumulated. With the increase in rainfall, the accumulated water source can drive the movement of heavier solid materials, which will cause a small amount of debris flow disaster.