Landslide susceptibility evaluation method based on terrain unit
1. The landslide susceptibility assessment method based on the terrain units is characterized by comprising the following steps:
s1, determining the position of the landslide point and related disaster-causing factor data;
s2, discretizing the disaster causing factor data to obtain discrete intervals of the disaster causing factor data, wherein each discrete interval corresponds to one classification type, and the distribution area occupation ratio AR of each disaster causing factor in each discrete interval is obtained;
s3, obtaining a weight ratio exponential curve according to the area occupation ratio AR of each disaster causing factor in each discrete interval, taking the intersection point and inflection point of the area ratio exponential curve and the weight ratio exponential curve as a judgment standard, and taking the inflection point or intersection point of the weight ratio exponential curve and the area ratio exponential curve as the endpoint value of the interval to obtain an optimized interval;
the weight ratio index QW is:
wherein H is the number of landslides in each classification category, H0The total number of the landslide points;
s4, counting the area of each discrete interval and the corresponding landslide point number after optimization, and obtaining a landslide inclination index as follows:
wherein TTsRepresenting the ratio of the total number of landslide points to the total area occupied by the disaster-causing factor data, TTaThe ratio of the total number of the landslide points existing in the optimized interval corresponding to the corresponding disaster-causing factor to the area of the interval a is as follows:
if OS >0, the slip inclination index OS is a positive value in the optimized interval range, and the possibility of occurrence of landslide is high;
if the OS is less than 0, the slip inclination index OS is a negative value in the optimized interval range, and the possibility of occurrence of landslide is low;
if OS is approximately equal to 0, the slip inclination index OS is close to zero in the optimized interval range, and whether the landslide occurs or not cannot be determined;
s4, constructing a river network by utilizing a terrain surface model DEM;
s5, calculating connection points in a river network to obtain water outlet points of a certain catchment basin, extracting catchment basins under positive and negative terrains respectively by combining water flow direction data, determining positions of all grid units of the catchment basins, and combining the catchment basins under the positive and negative terrains to obtain a terrain unit;
s6, assigning the slip inclination index of each optimized interval of each disaster causing factor to each terrain unit;
and S7, constructing a random forest model, inputting data of each terrain unit, and outputting a landslide probability value of each terrain unit.
2. The method for evaluating the landslide susceptibility according to claim 1, wherein in S2, the continuous disaster factor data is discretized by using an equidistant method, and the area integral analogy AR is:
wherein, YiArea corresponding to the ith classification category of the disaster-causing factor data, YsThe overall area occupied by the disaster-causing factor data is shown.
3. The method for evaluating the vulnerability of landslide based on terrain unit of claim 2, wherein in the S4, constructing a river network by using a terrain surface model DEM comprises:
s41, setting different simulated depression depth values in the river region according to the terrain characteristics according to the terrain surface model DEM, and filling depressions in the terrain surface model DEM by using the simulated depression depth values;
s42, simulating the water flow direction in the terrain surface model DEM:
the flow direction average index model is used as a real water flow direction judgment basis, and the flow direction average index formula is as follows:
wherein G isiFor the straight-line distance, G, between the centers of two grid elements in the terrain surface model DEMzFor the elevation difference between two grid cells in the terrain surface model DEM:
if (Max FI) >0, it indicates that the flow direction between two grid cells is uncertain;
if (Max FI) ≦ 0 and the maximum is only one, then the direction corresponding to this maximum is taken as the flow direction between the two grid elements;
if max (fi) is 0 and there is more than one 0 value, the direction values corresponding to these 0 values are added together to form the flow direction between two grid cells;
if Max (FI) <0 and more than one maximum value exists, selecting the direction corresponding to one maximum value as the flow direction between two grid units;
s43, in the terrain surface model DEM, simulating the confluence:
giving water levels of different depths to each grid node of the terrain surface model DEM, and calculating the water volume of each grid unit in all directions according to the regional terrain and the water flow direction according to the natural law that water flows from high to low to obtain confluence multidirectional cumulant;
and S44, setting different catchment area values according to the confluence multidirectional cumulant of each grid unit obtained in S43 to simulate river networks under different conditions, comparing the catchment area values with real river networks to find out the catchment area value closest to the real river networks, and generating the river networks by using the catchment area value.
4. The terrain unit-based landslide susceptibility assessment method according to claim 3, wherein the confluent multi-directional accumulation amount is:
DU＝D1+D2+D3+...+DN(ii) a Wherein D isnN is 1,2,3, and N denotes the number of water flow directions of a central grid unit.
5. The terrain-unit-based landslide susceptibility assessment method of claim 4, wherein Dn＝DnAnd M, wherein M is 1 or 0, if the flow direction points to the grid unit which needs to calculate the confluence multidirectional accumulation amount at present, M is 1, and otherwise M is 0.
6. The terrain cell-based landslide susceptibility assessment method according to claim 1, wherein in S1, the disaster factor data comprises continuous data and discrete data, the continuous data is elevation, slope direction, distance from river, distance from road, section curvature, NDVI, and the discrete data is rock-soil body type, land use type, and slope structure.
Background
The landslide susceptibility evaluation is the basis for realizing landslide geological disaster prevention and control, and the landslide disaster-causing factor and the landslide evaluation unit division are two most important components in the landslide susceptibility evaluation. The landslide disaster-causing factor analysis is the primary work of risk assessment and prevention and control of regional landslide disasters, and meanwhile, the landslide susceptibility evaluation method is poor in reliability due to the fact that some problems still exist in the existing method for dividing landslide susceptibility evaluation units and no uniform division standard exists.
The landslide occurrence is a relatively complex process, and the establishment of a scientific incidence evaluation model is very important for the landslide incidence evaluation result.
Disclosure of Invention
The invention provides a landslide susceptibility evaluation method based on a terrain unit, aiming at the problem that the existing landslide susceptibility evaluation method is poor in reliability.
The invention discloses a landslide susceptibility evaluation method based on a terrain unit, which comprises the following steps:
s1, determining the position of the landslide point and related disaster-causing factor data;
s2, discretizing the disaster causing factor data to obtain discrete intervals of the disaster causing factor data, wherein each discrete interval corresponds to one classification type, and the distribution area occupation ratio AR of each disaster causing factor in each discrete interval is obtained;
s3, obtaining a weight ratio exponential curve according to the area occupation ratio AR of each disaster causing factor in each discrete interval, taking the intersection point and inflection point of the area ratio exponential curve and the weight ratio exponential curve as a judgment standard, and taking the inflection point or intersection point of the weight ratio exponential curve and the area ratio exponential curve as the endpoint value of the interval to obtain an optimized interval; the weight ratio index QW is:
wherein H is the number of landslides in each classification category, H0The total number of the landslide points;
s4, counting the area of each discrete interval and the corresponding landslide point number after optimization, and obtaining a landslide inclination index as follows:
wherein TTsRepresenting the ratio of the total number of landslide points to the total area occupied by the disaster-causing factor data, TTaThe ratio of the total number of the landslide points existing in the optimized interval corresponding to the corresponding disaster-causing factor to the area of the interval a is as follows:
if OS >0, the slip inclination index OS is a positive value in the optimized interval range, and the possibility of occurrence of landslide is high;
if the OS is less than 0, the slip inclination index OS is a negative value in the optimized interval range, and the possibility of occurrence of landslide is low;
if OS is approximately equal to 0, the slip inclination index OS is close to zero in the optimized interval range, and whether the landslide occurs or not cannot be determined;
s4, constructing a river network by utilizing a terrain surface model DEM;
s5, calculating connection points in a river network to obtain water outlet points of a certain catchment basin, extracting catchment basins under positive and negative terrains respectively by combining water flow direction data, determining positions of all grid units of the catchment basins, and combining the catchment basins under the positive and negative terrains to obtain a terrain unit;
s6, assigning the slip inclination index of each optimized interval of each disaster causing factor to each terrain unit;
and S7, inputting the data in each terrain unit and outputting the landslide probability value of each terrain unit to the constructed random forest model.
Preferably, in S2, the continuous disaster factor data is discretized by an equidistant method, and the area fraction analogy AR is:
wherein, YiArea corresponding to the ith classification category of the disaster-causing factor data, YsThe overall area occupied by the disaster-causing factor data is shown.
Preferably, in S4, the constructing a river network by using the terrain surface model DEM includes:
s41, setting different simulated depression depth values in the river region according to the terrain characteristics according to the terrain surface model DEM, and filling depressions in the terrain surface model DEM by using the simulated depression depth values;
s42, simulating the water flow direction in the terrain surface model DEM:
the flow direction average index model is used as a real water flow direction judgment basis, and the flow direction average index formula is as follows:
wherein G isiFor the straight-line distance, G, between the centers of two grid elements in the terrain surface model DEMzFor the elevation difference between two grid cells in the terrain surface model DEM:
if (Max FI) >0, it indicates that the flow direction between two grid cells is uncertain;
if (Max FI) ≦ 0 and the maximum is only one, then the direction corresponding to this maximum is taken as the flow direction between the two grid elements;
if max (fi) is 0 and there is more than one 0 value, the direction values corresponding to these 0 values are added together to form the flow direction between two grid cells;
if Max (FI) <0 and more than one maximum value exists, selecting the direction corresponding to one maximum value as the flow direction between two grid units;
s43, in the terrain surface model DEM, simulating the confluence: giving water levels of different depths to each grid node of the terrain surface model DEM, and calculating the water volume of each grid unit in all directions according to the regional terrain and the water flow direction according to the natural law that water flows from high to low to obtain confluence multidirectional cumulant;
and S43, setting different catchment area values according to the confluence multidirectional cumulant of each grid unit obtained in S43 to simulate river networks under different conditions, comparing the catchment area values with real river networks to find out the catchment area value closest to the real river networks, and generating the river networks by using the catchment area value.
Preferably, the confluence multidirectional accumulation amount is:
DU＝D1+D2+D3+...+DN(ii) a Wherein D isnN is 1,2,3, and N denotes the number of water flow directions of a central grid unit.
Preferably, D isn＝DnAnd M, wherein M is 1 or 0, if the flow direction points to the grid unit which needs to calculate the confluence multidirectional accumulation amount at present, M is 1, and otherwise M is 0.
Preferably, in S1, the disaster factor data includes continuous data and discrete data, the continuous data is elevation, slope direction, river distance, road distance, section curvature, NDVI, and the discrete data is rock-soil body type, land use type, and slope structure.
The invention has the beneficial effects that: the method constructs a landslide susceptibility evaluation model based on the terrain units, and has the following characteristics: firstly, the method of the invention fully considers the factors of various disaster-causing influences and ensures the comprehensive enrichment of the invention content to a certain extent. Secondly, the method can carry out quantitative analysis on the landslide disaster-causing influence factors, thereby providing a foundation for landslide prevention and control work. In addition, the method selects the evaluation unit based on the terrain on the basis of fully considering the terrain and the landform, and ensures the relatively reliable research result.
Drawings
FIG. 1 is a flow chart of the present invention;
FIG. 2 is a plot of the area of the example zone positive topography simulated catchment area values versus the corresponding valley density;
FIG. 3 is a plot of an embodiment regional random forest model ROC based on terrain units;
FIG. 4 is a plot of an embodiment grade of landslide susceptibility based on terrain cells.
Detailed Description
The technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the drawings in the embodiments of the present invention, and it is obvious that the described embodiments are only a part of the embodiments of the present invention, and not all of the embodiments. All other embodiments, which can be derived by a person skilled in the art from the embodiments given herein without making any creative effort, shall fall within the protection scope of the present invention.
It should be noted that the embodiments and features of the embodiments may be combined with each other without conflict.
The invention is further described with reference to the following drawings and specific examples, which are not intended to be limiting.
The landslide susceptibility evaluation method based on the terrain unit of the embodiment includes:
step one, determining a landslide point position and related disaster-causing factor data;
and step two, dividing the selected disaster causing factors into continuous data and discrete data, and uniformly substituting the discrete data into the operation. Discretizing the disaster causing factor data by using an equidistant method to obtain discrete intervals of the disaster causing factor data, wherein each discrete interval corresponds to one classification type to obtain the distribution area occupation ratio AR of each disaster causing factor in each discrete interval;
the area integral analogy AR is:
wherein, YiArea corresponding to the ith classification category of the disaster-causing factor data, YsThe overall area occupied by the disaster-causing factor data is shown.
Step three, acquiring a weight ratio exponential curve according to the area occupation ratio AR of each disaster causing factor in each discrete interval, taking the intersection point and inflection point of the area ratio exponential curve and the weight ratio exponential curve as a judgment standard, and taking the inflection point or intersection point of the weight ratio exponential curve and the area ratio exponential curve as an end point value of an interval to obtain an optimized interval;
the weight ratio index QW is:
wherein H is the number of landslides in each classification category, H0The total number of the landslide points;
s4, counting the area of each discrete interval and the corresponding landslide point number after optimization, and obtaining a landslide inclination index as follows:
wherein TTsRepresenting the ratio of the total number of landslide points to the total area occupied by the disaster-causing factor data, TTaThe ratio of the total number of the landslide points existing in the optimized interval corresponding to the corresponding disaster-causing factor to the area of the interval a is as follows:
if OS >0, the slip inclination index OS is a positive value in the optimized interval range, and the possibility of occurrence of landslide is high;
if the OS is less than 0, the slip inclination index OS is a negative value in the optimized interval range, and the possibility of occurrence of landslide is low;
if OS is approximately equal to 0, the slip inclination index OS is close to zero in the optimized interval range, and whether the landslide occurs or not cannot be determined;
fourthly, constructing a river network by utilizing a terrain surface model DEM;
calculating connection points in the river network to obtain water outlet points of a certain catchment basin, extracting catchment basins under positive and negative terrains respectively by combining water flow direction data, determining positions of all grid units of the catchment basins, and combining the catchment basins under the positive and negative terrains to obtain a terrain unit;
in this step, the water outlet point of a certain catchment basin, i.e. the lowest point of the catchment basin, can be obtained by calculating the river connection point, and all the grids flowing through the water outlet port at the upstream of the water outlet point are analyzed and searched by combining the water flow direction data until all the grids of the catchment basin determine the position, which represents the position of the boundary and watershed of the catchment basin. And inputting water flow direction data and river connection point data to output a catchment basin grid, and obtaining a vector file of the catchment basin by using a grid vector conversion function in ArcGIS. And (4) reversing the DEM to obtain a negative terrain, extracting a catchment basin by the same steps, and merging catchment basins obtained by two times of analysis to obtain a terrain unit.
Assigning the slip inclination index of each optimized interval of each disaster causing factor to each terrain unit;
and seventhly, inputting the data of each grid unit in the terrain units and outputting the landslide proneness probability value of each grid unit to the constructed random forest model. And evaluating the classification effect of the model by using the AUC value of the area under the ROC (receiver Operating Characteristic curve), wherein the larger the AUC value is, the better the classification effect is.
In the fourth step of the embodiment, the method for constructing the river network by using the terrain surface model DEM comprises the following steps:
fourthly, according to the terrain surface model DEM, different simulated depression depth values are set in the river region due to terrain characteristics, and depression filling is carried out on depressions in the terrain surface model DEM by utilizing the simulated depression depth values;
the DEM is a smoother terrain surface model, due to DEM errors and real terrains, the surface of the DEM has some sunken areas, unreasonable and even wrong water flow directions are often caused when water flow direction calculation is carried out, and therefore depression filling is carried out on original DEM data before the water flow direction is calculated, and the DEM without depressions is obtained. Carrying out depression filling on the original DEM data by adopting a simulated depression depth value, wherein the judgment formula of the simulated depression depth value is
Where B is a simulated depression profile obtained by setting different depression depths in consideration of topographic features of the study area, and B is a depression profile of real topography of the study area. When the value of N approaches 1, the simulated depression depth is closest to the actual depression distribution, thereby obtaining a simulated depression depth value capable of filling the depression;
step four, simulating the water flow direction in the terrain surface model DEM:
the flow direction average index model is used as a real water flow direction judgment basis, and the flow direction average index formula is as follows:
wherein G isiFor the straight-line distance, G, between the centers of two grid elements in the terrain surface model DEMzFor the elevation difference between two grid cells in the terrain surface model DEM:
if (Max FI) >0, it indicates that the flow direction between two grid cells is uncertain;
if (Max FI) ≦ 0 and the maximum is only one, then the direction corresponding to this maximum is taken as the flow direction between the two grid elements;
if max (fi) is 0 and there is more than one 0 value, the direction values corresponding to these 0 values are added together to form the flow direction between two grid cells;
if Max (FI) <0 and more than one maximum value exists, selecting the direction corresponding to one maximum value as the flow direction between two grid units;
step four, in the terrain surface model DEM, simulating convergence:
giving water levels with different depths at each grid node of the terrain surface model DEM, calculating the water volume of each grid unit in all directions according to the regional terrain and the water flow direction according to the natural law that water flows from high to low, and obtaining confluence multidirectional cumulant as follows:
DU＝D1+D2+D3+...+DN(ii) a Wherein D isnN is 1,2,3, and N denotes the number of water flow directions of a central grid unit.
Dn＝DnAnd M, wherein M is 1 or 0, if the flow direction points to the grid unit which needs to calculate the confluence multidirectional accumulation amount at present, M is 1, and otherwise M is 0.
And fourthly, setting different catchment area values according to the obtained confluence multidirectional cumulant of each grid unit to simulate the river network under different conditions, comparing the catchment area values with the real river network, finding out the catchment area value closest to the real river network, and generating the river network by using the catchment area value.
Wherein C is analog digital river network data obtained by setting different catchment area values, and C is real digital river network data. When the value of Q is close to 1, the corresponding digital river network under the catchment area value is close to the real digital river network, and the catchment area value at the moment can be used as the input data of the generated river network.
In the first step of the embodiment, the disaster-causing factor data includes continuous data and discrete data, the continuous data includes elevation, slope direction, distance from river, distance from road, section curvature, and NDVI, and the discrete data includes rock-soil body type, land utilization type, and slope structure.
The specific embodiment is as follows: in order to achieve the above-mentioned prediction and planning objectives, a landslide susceptibility assessment method based on a terrain unit is provided. By simultaneously considering the evaluation unit and the disaster-causing factor, the random forest model is used as the landslide susceptibility evaluation model, and reasonable prediction of future landslide susceptibility sections is provided. The calculation scheme adopted by the present invention is further explained in several steps with reference to the attached figures and specific examples:
1. data acquisition:
the data used in this embodiment includes 1:1 ten thousand county level second tone data in a certain city and a certain province, 1: 650000 geological disaster distribution diagram, GDEMDEM 30m resolution digital elevation data of geographic space data cloud, GF-1 remote sensing image and rainfall site data of Hubei provincial weather bureau.
2. The selected disaster-causing factors are divided into continuous data and discrete data, and the discrete data are used in a unified mode for substitution operation. In order to be convenient and then uniformly substituted into a model, firstly, discretizing the continuous disaster causing factors in the selected disaster causing factors by using an equal-spacing method, substituting discretization data into a surface integral analogy formula to count the distribution area occupation ratio of each disaster causing factor in each discrete interval, and using a surface integral analogy AR;
3. the intersection point and inflection point of the surface integral analog curve and the weight ratio exponential curve are used as a judgment standard to further optimize the discretization interval, and the landslide data and the surface integral analog value data calculated by each disaster causing factor are substituted into a surface integral analog formula;
4. and (4) counting the area of each discrete interval and the corresponding number of landslides after the area classification ratio and the weight ratio index are optimized, and calculating a slip inclination index to substitute the slip inclination index into a random forest evaluation model, so that the disaster tendency of each disaster-causing factor is calculated.
5. First, according to the embodiment, the relative height difference of the area is between 500 m and 1300 m. And comparing the embodiment area depression depth distribution map calculated in the ArcGIS platform with the depression depth distribution map of the real terrain, so as to calculate the depression depth of the positive terrain. And subtracting the elevation value of the positive terrain from 0 to obtain the negative terrain, and calculating the filling threshold of the negative terrain by referring to a positive terrain hollow filling threshold calculation method. Finally, the hole filling threshold for positive topography of the example zone is set to 500, and the hole filling threshold for negative topography is set to 8000.
6. The GDEMDEM 30m resolution digital elevation data of the embodiment area is substituted into the flow direction average index model for calculation
Therefore, the water flow direction of each grid unit is determined, and the real water flow direction is simulated.
Wherein G isiIs the linear distance between the centers of two grid cells, GzIs the difference in elevation between two grid cells:
if (Max FI) is greater than 0, the grid unit flow direction is undetermined;
if (Max FI) ≦ 0 and the maximum is only one, then the direction value corresponding to this maximum is taken as the flow direction value at the center grid cell;
if max (fi) is 0 and there is more than one 0 value, the direction values corresponding to these 0 values are added to obtain the flow value.
7. In order to simulate the real confluence flow, adopting confluence multidirectional cumulant and counting the water quantity of each grid unit in all directions by using the obtained water flow direction, namely inputting the water flow corresponding to the grid adjacent to each central unit into the confluence multidirectional cumulant;
based on the obtained confluence multidirectional cumulant, different catchment area values are set to simulate river networks under different conditions, and compared with real river networks, when the value of a catchment area value judgment formula Q approaches to 1, the corresponding digital river network under the catchment area value approaches to the real digital river network, and the catchment area value at the moment can be used as input data for generating the river networks.
In this embodiment, firstly, based on an ArcGIS platform and in combination with the above operations, a filled non-hollow DEM is obtained, and 9 values such as a normal terrain catchment area threshold value of 500, 1000, 1500 … 12000 are respectively taken; and taking 9 values of inverse terrain catchment area threshold values of 500, 1000, 1500 … 12000 and the like, calculating corresponding valley density and comparing the valley density with the density of the real terrain river network, wherein the specific calculation result is shown in fig. 2. Finally, the present example sets the catchment area threshold for the positive terrain to 12000. Similarly, the catchment area threshold of the negative topography is selected to be 8000. Inputting the positive and negative terrain simulated catchment area value obtained by the calculation into corresponding modules in ArcGIS, respectively extracting catchment watersheds of the positive and negative terrains, and combining the catchment watersheds and the catchment watersheds to obtain a terrain unit corresponding to the embodiment area.
And taking the obtained terrain units as basic units for landslide susceptibility evaluation, assigning the sensitivity value of each optimized interval of each disaster-causing factor to each terrain unit, bringing the sensitivity value into a constructed random forest model, and evaluating the classification effect of the model by utilizing the AUC value of the area under the ROC curve. The AUC value is between 0.5 and 1.0, and the classification effect is better when the AUC value is larger. The AUC value of the model in the embodiment is 0.82 (figure 2), and the value meets the requirement of spatial distribution verification accuracy.
In this embodiment, based on 9 landslide disaster-causing factors in a county, a random forest model is adopted to perform landslide easiness evaluation on an embodiment area through repeated adjustment of key model parameters, so that a probability value of easiness of occurrence of each grid of the embodiment area is determined, and accordingly, a landslide easiness distribution map of the embodiment area is drawn by using a natural discontinuity method, and the result is shown in table 1 and fig. 4. The results show that the example zone susceptibility rating is 77.84% for the area below medium susceptibility and that the landslide occurs 55.92% for the very high and higher susceptibility zones. The higher incidence areas and the high incidence areas are mostly distributed in the areas of the Yangtze river and the two sides of the tributary of the Yangtze river, such as the river town of the Guojia, the dam town of the Xiangxi, the village town of the deserts and the like. It can be known from on-site investigation and comparison of historical landslide data that the result of landslide development evaluation based on topographic unit is reasonable and reliable, and has higher degree of coincidence with ziguo county historical landslide distribution characteristics.
TABLE 1 number of landslide hazard points and grading units corresponding to different susceptibility grades in example area
