Adaptive slope cell extraction method based on watershed mechanism

Abstract

The adaptive slope unit extraction method based on the watershed mechanism comprises the following steps: 1) preprocessing a digital elevation map to generate a feature map; 2) improving an initial coarse slope unit generated by a watershed algorithm; 3) extracting corresponding features of the slope unit through the feature map, constructing an adjacency relationship and a feature table; 4) traversing the slope unit to find a line segment table that needs to be dissolved; 5) updating the line segment table and marking the line segments that can be dissolved as a uniform value; 6) marking the line segments on the slope unit, dissolving the slope units with the same mark, returning to step three for iteration until the iteration termination condition is met. The present application introduces a flow convergence accumulation to control the dissolution of the slope unit, so that the dissolution result is closer to the geographical features; the multiple dissolution is used to improve the efficiency of the slope unit dissolution and reduce the number of iterations; the number of parameters in the dissolution process is reduced, so that the slope unit is no longer affected by multiple parameters.

Description

Technical Field

[0001] This invention belongs to the field of geological hazard zoning, specifically involving an adaptive slope unit extraction method based on the watershed mechanism. Background Technology

[0002] Commonly used unit types in geological hazard zoning include grid units, regional units, homogeneous condition units, sub-basin units, and slope units. Among these, grid units have relatively regular shapes, facilitating rapid subdivision, and the resulting matrix data is beneficial for further calculations; however, they cannot fully reflect topographic relief and lack a close connection to geological conditions. Homogeneous condition units do not consider the differences in geological environmental conditions across different regions. Sub-basin units are suitable for debris flow hazard zoning but not for landslides and collapses. Slope units are the basic units for the development of geological hazards such as landslides and collapses. Furthermore, among various controlling or influencing factors, the developmental stages of rivers and valleys have a significant controlling effect on the formation of landslides and collapses. Therefore, using slope units as evaluation units allows for close correlation with geological environmental conditions, comprehensively reflecting the effects of various controlling or influencing factors, and making the evaluation results closer to reality. Therefore, provided that DEM accuracy requirements are met, slope unit subdivision is more suitable for geological hazard zoning.

[0003] The following are some literatures in this field that study the division of this slope unit: [1] Turel M, Frost JD. Delineation of Slope Profiles from Digital Elevation Models for LandslideHazard

[0004] Analysis[J]. American Society of Civil Engineers, 2011:829-836; [2] Yan Ge, Liang Shouyun, Zhao Hongliang. Improvement and implementation of slope unit division method based on GIS[J]. Geographical Science, 2017, 37(11):1764-1770; [3] Alvioli M, Marchesini I, Reichenbach P, et al. Automatic delineation of geomorphological slope units with r.slopeunits v1. 0 and their optimization for landslide susceptibility modeling[J]. Geoscientific Model Development, 2016, 9(11): 3975; [4] Wang K, Zhang SJ, Ricardo DelgadoTéllez, et al. A new slope unit extraction method for regional landslide analysis based on morphological image analysis[J]. Bulletin of Engineering Geology and the Environment, 2018:1-13.

[0005] Among them, Reference [1] extracts slope units by generating a depression-free DEM, extracting flow direction and runoff accumulation, generating a river network, extracting positive and negative catchment basins, and merging catchment basins. Reference [2] proposes a curvature-based method, first using the DEM to obtain positive and negative curvature, then using the positive and negative curvature to perform watershed segmentation and watershed calculation, and extracting convex and concave landform elements. The convex and concave landform elements are vectorized using the "raster to surface" tool provided by the GIS toolbox, and then merged to form slope units. Reference [3] first uses the r.watershed hydrological module embedded in GrassGIS to perform runoff accumulation analysis, dividing the DEM raster data into a few larger sub-basins. Then, by reducing the runoff accumulation threshold, each sub-basin is divided into two parts, which Alvioli calls "half basin" (HB). The process of repeatedly reducing the threshold of the cumulative flow will subdivide HB into many HBchild small regions. Each HBchild small region is composed of several grid cells. The boundary line of the small region is the watershed or the watershed. During the subdivision process, if the area of ​​a certain HBchild small region is less than the minimum area threshold a set by the user, then the small region is regarded as a slope cell. If the area of ​​the HBchild small region is less than the maximum area threshold m set by the user, and the circular squareness of the grid cells inside the small region meets the set threshold c, the subdivision process stops, and the HBchild small region is regarded as a slope cell. Reference [4] first discretizes DEM into many grid center points, and uses the fourth surface model method to calculate the average curvature of each grid center point. According to the physical meaning of the average curvature, the convex and concave shape of the grid cell is judged, and the terrain is divided into two categories: ridge and valley regions. The extraction results of the ridge and valley regions are binarized, and the morphological skeleton lines of the ridge and valley are extracted using morphological imaging algorithms. The morphological skeleton lines can reflect the undulation characteristics of the actual micro-terrain and can be used as the boundary lines of the small regions. Then, an iterative search method is used to connect the morphological skeleton lines of ridges and valleys into a closed morphological skeleton network. The actual terrain covered by each small region in the network has a uniform slope and aspect. Principal component analysis is used to extract the fitting plane of the closed small regions, and vector similarity theory is used to merge adjacent small regions whose normal vector similarity meets a set threshold to form slope units.

[0006] However, the conventional slope element division has the following problems:

[0007] (1) It is not difficult to apply to areas lacking high-precision DEMs. The higher the DEM accuracy, the finer the extracted slope units, and the better the slope unit boundary lines can match the obvious turning points of the actual landform. For areas lacking high-precision DEMs, the number of slope units extracted is small, and the extraction results will be relatively coarse, affecting subsequent landslide stability analysis;

[0008] (2) The extraction results are affected by the threshold set by the user. In all existing methods, users need to set some empirical thresholds as the conditions for stopping the iteration. The increase of user-defined parameters will increase the subjectivity of the extraction results and reduce the reproducibility of the extraction results.

[0009] (3) There is a lack of simple and easy-to-implement optimization methods for the extracted results. Existing methods lack simple and easy-to-implement optimization steps. For example, the hydrological process analysis method based on forward and reverse DEMs produces a large number of long, unreasonable strip-shaped units, requiring extensive and tedious manual correction later. When the study area is large and the number of such unreasonable units is high, the workload of manual correction will increase dramatically. The r.slopeunits method uses a relatively complex method to optimize the extracted results. Each extraction requires running a landslide judgment model and plotting the ROC curve.

[0010] Calculating the area of ​​AUCROC takes a considerable amount of time.

[0011] Therefore, proposing a method with a reasonable threshold and superior extraction results has profound significance in the field of slope unit partitioning. Summary of the Invention

[0012] To address the problems existing in the prior art, this invention proposes an adaptive slope element extraction method based on the watershed mechanism.

[0013] The adaptive slope cell extraction method based on the watershed mechanism proposed in this invention is applicable to slope cell division and regional early warning system application scenarios. The specific steps are described as follows:

[0014] Step 1: First, preprocess the digital elevation map to generate a reverse elevation map. Then, fill in depressions in both elevation maps. Using the Grass GIS Python library functions, calculate the aspect raster, slope raster, and flow direction raster. Calculate the cosine and sinine values ​​of the aspect (i.e., the horizontal component sinno and the vertical component cosno) using a raster calculator. These values ​​are used for subsequent aspect consistency calculations. The confluence accumulation algorithm is used to calculate the confluence accumulation values ​​su_acc and re_acc for both the forward and reverse elevation maps. The confluence accumulation algorithm is as follows:

[0015] Input flows to the raster map dir.

[0016] (1.1) Initialize the inflow matrix NIDP as a zero matrix of the same size as dir, and set the nodata part of dir to nodata in NIDP as well;

[0017] (1.2) Using the NIDP-base algorithm, traverse the dir to calculate NIDP. For each cell that flows into the traversed cell, the value of that cell is incremented by 1.

[0018] (1.3) Initialize the flow accumulator matrix FlowAccu, a matrix of all ones with the same size as dir;

[0019] (1.4) Traverse the NIDP matrix, find the cell with an NIDP value of 0 as the starting point, and according to the flow direction of the dir matrix, flow to the position in NIDP that is not 1, and decrement the NIDP value at that position by one. Add the value of the outflow position to the value of the FlowAccu matrix.

[0020] (1.5) Export FlowAccu tif

[0021] Step 2: Using the Watershed algorithm in Grass GIS, a minimum area threshold is set (an optimal value is found through bisection and a predetermined slope unit density). Initial coarse slope units (SLU) are extracted. In most practical applications, flat slope areas are not considered, so flat slope areas with a slope of less than 5 degrees are extracted as masks. The fragmented parts of the masks are removed, and the masked parts of the SLU are cropped. If there is an administrative boundary vector map, it is overlaid on the SLU. Each slope unit is numbered according to its administrative division. If there is no administrative boundary vector map, the administrative division number of each slope unit is set to 0. Fragmented slope units caused by calculation errors are dissolved. Finally, the slope unit boundary lines are smoothed to obtain the initial version of the slope units. The slope unit density of this version is slightly larger than the target density to facilitate subsequent refinement and dissolution.

[0022] Step 3: Extract slope unit features and line segment features. First, extract the surface features corresponding to the slope unit based on the feature raster map in Step 1, such as area, slope aspect, administrative division, etc. Then, construct the adjacency relationship table of the slope unit to extract the intersection line of the adjacent unit. Extract the confluence accumulation of the intersection line through su_acc and re_acc. Quantize the values ​​in su_acc and re_acc and take the larger value of the two to obtain the confluence accumulation feature of the intersection line.

[0023] Step 4: Traverse the slope cells. Based on the area features obtained in Step 3, determine whether the slope cells need to be dissolved. Calculate the slope aspect consistency of each edge of the slope cells that need to be dissolved, and select the edges that meet the slope aspect consistency threshold. The slope aspect consistency threshold is calculated as follows:

[0024]

[0025] Here, is the threshold of a certain slope unit, k is the threshold control coefficient, the smaller the value, the easier it is to be dissolved, and is the result of the area standardization of a certain slope unit. Then, from these edges, edges that are consistent with the administrative division of adjacent slope units are selected, and the edge with the largest cumulative flow is selected as the edge to be dissolved and stored in flush_line.

[0026] Step 5: Input flush_line to find the slope cells that can be continuously dissolved, reducing the number of iterations. left represents the slope cell to the left of the line segment, and right represents the right.

[0027] (5.1) Initialize is_fush to -1, indicating that it has not been modified;

[0028] (5.2) Traverse the fush_line, skipping the line segments where is_fush is 1;

[0029] (5.3) Add the unmodified line segments from the left neighbor of the line segment to the queue;

[0030] (5.4) Dequeue the line segment in the queue, set the fusion flag to the id of left, set is_fush to 1, and add the neighboring line segments of the line segment to the queue until the queue is empty;

[0031] (5.5) After traversal, update flush_line;

[0032] Updating `fush_line` allows some consecutive slope cells to dissolve simultaneously, reducing the number of dissolution iterations and increasing the dissolution rate.

[0033] Step 6: Set the dissolution markers of the left and right ramp cells of the line segment in `fush_line` to the dissolution markers of that line segment. Then dissolve the ramp cells with the same dissolution markers and calculate the density of the dissolved ramp cells. If the density meets the target density or reaches the maximum number of iterations, output the ramp cells; otherwise, return to Step 3, extract the features of the dissolved ramp cells, and dissolve them again. Note that to improve the calculation speed, the intersection features that have not changed do not need to be extracted again.

[0034] This invention proposes a slope unit extraction method based on existing watershed methods. This invention is applicable to regional risk prediction systems based on slope units, particularly for landslide early warning in mountainous areas. It extracts a large number of coarse slope units using the watershed algorithm, and iteratively dissolves the slope unit data to obtain slope unit results with consistent slope aspect that do not disrupt the geographical structure. This method reduces terrain damage caused by the threshold setting of the watershed algorithm and improves the uniformity among slope units. Compared to existing methods, this invention makes the following improvements:

[0035] 1. An improvement based on Alvioli's r.slopeunits method is made. The drawback of r.slopeunits is that the initial threshold is set too high when using the watershed algorithm, causing slope units to cross valleys and ridgelines. This method finds a smaller threshold that meets the requirements through binary search. Although more slope units are extracted at this time, the number of slope units can be reduced through subsequent dissolution. When dissolving, the flow accumulation at the slope unit boundary is considered, which can prevent slope units from crossing valleys or ridgelines.

[0036] 2. To improve the rate of slope unit dissolution, in each iteration, the line segments to be dissolved are traversed, and breadth-first search is used to dissolve some connected slope units. The dissolved slope units are then marked, so that when extracting features at the beginning of each iteration, only the features of these marked slope units need to be extracted, thereby avoiding repeated feature extraction and reducing the dissolution time.

[0037] 3. In addition, in order to reduce the influence of artificial thresholds on slope units, this method reduces the use of parameters and sets the area and aspect consistency as an aspect consistency threshold to control the dissolution of slope units. The threshold is controlled by a coefficient k, the smaller the value, the easier it is to be dissolved. Furthermore, the value of the confluence accumulation is not set as a threshold, but is selected by selecting the maximum confluence accumulation. This method can also prevent the ridges and valleys from being dissolved.

[0038] The advantages of this invention are: (1) It introduces the accumulation of flow to control the dissolution of slope units, making the dissolution results closer to the geographical features. (2) It uses multiple dissolution to improve the efficiency of slope unit dissolution, which can dissolve a large area of ​​slope units that meet the conditions, reducing the number of iterations. (3) It reduces the number of parameters in the dissolution process, so that the slope units are no longer affected by multiple parameters. Attached Figure Description

[0039] Figure 1 This is a flowchart of the overall algorithm of the present invention.

[0040] Figure 2 This is a flowchart of the confluence accumulation algorithm.

[0041] Figure 3 This is a flowchart of the dissolution process for the slope unit.

[0042] Figure 4 shows an example of the ramp element obtained by the present invention, wherein... Figure 4a It is a slope aspect feature map superimposed on slope units. Figure 4b It is a map of slope units superimposed with contour lines. Figure 4c This is a diagram of a slope element. Figure 4dThis is a diagram of the key areas where slope units are superimposed. Detailed Implementation

[0043] The preferred embodiments of the present invention are given below with reference to the accompanying drawings to illustrate the technical solution of the present invention in detail. The embodiments described herein are for illustrative purposes only and do not limit the scope of the invention.

[0044] like Figure 1 As shown, the adaptive slope unit extraction method based on the watershed mechanism proposed in this invention is applicable to slope unit division and regional early warning system application scenarios. The specific steps are as follows:

[0045] Step 1: First, preprocess the digital elevation map to generate a reverse elevation map. Then, fill in depressions in both elevation maps. Using the Grass GIS Python library functions, calculate the aspect raster, slope raster, and flow direction raster. Calculate the cosine and sinine values ​​of the aspect (i.e., the horizontal component sinno and the vertical component cosno) using a raster calculator. These values ​​are used for subsequent aspect consistency calculations. The confluence accumulation algorithm is used to calculate the confluence accumulation values ​​su_acc and re_acc for both the forward and reverse elevation maps. The confluence accumulation algorithm is as follows:

[0046] Input flow to raster map dir

[0047] (1.1) Initialize the inflow matrix NIDP as a zero matrix of the same size as dir, and set the nodata part of dir to nodata in NIDP as well;

[0048] (1.2) Using the NIDP-base algorithm, traverse the dir to calculate NIDP. For each cell that flows into the traversed cell, the value of that cell is incremented by 1.

[0049] (1.3) Initialize the flow accumulator matrix FlowAccu, a matrix of all ones with the same size as dir;

[0050] (1.4) Traverse the NIDP matrix, find the cell with an NIDP value of 0 as the starting point, and according to the flow direction of the dir matrix, flow to the position in NIDP that is not 1, and decrement the NIDP value at that position by one. Add the value of the outflow position to the value of the FlowAccu matrix.

[0051] (1.5) Export FlowAccu tif

[0052] Step 2: Using the Watershed algorithm in Grass GIS, a minimum area threshold is set (an optimal value is found through bisection and a predetermined slope unit density). Initial coarse slope units (SLU) are extracted. In most practical applications, flat slope areas are not considered, so flat slope areas with a slope of less than 5 degrees are extracted as masks. The fragmented parts of the masks are removed, and the masked parts of the SLU are cropped. If there is an administrative boundary vector map, it is overlaid on the SLU. Each slope unit is numbered according to its administrative division. If there is no administrative boundary vector map, the administrative division number of each slope unit is set to 0. Fragmented slope units caused by calculation errors are dissolved. Finally, the slope unit boundary lines are smoothed to obtain the initial version of the slope units. The slope unit density of this version is slightly larger than the target density to facilitate subsequent refinement and dissolution.

[0053] Step 3: Extract slope unit features and line segment features. First, extract the surface features corresponding to the slope unit based on the feature raster map in Step 1, such as area, slope aspect, administrative division, etc. Then, construct the adjacency relationship table of the slope unit to extract the intersection line of the adjacent unit. Extract the confluence accumulation of the intersection line through su_acc and re_acc. Quantize the values ​​in su_acc and re_acc and take the larger value of the two to obtain the confluence accumulation feature of the intersection line.

[0054] Step 4: Traverse the slope cells. Based on the area features obtained in Step 3, determine whether the slope cells need to be dissolved. Calculate the slope aspect consistency of each edge of the slope cells that need to be dissolved, and select the edges that meet the slope aspect consistency threshold. The slope aspect consistency threshold is calculated as follows:

[0055]

[0056] Here, is the threshold of a certain slope unit, k is the threshold control coefficient, the smaller the value, the easier it is to be dissolved, and is the result of the area standardization of a certain slope unit. Then, from these edges, edges that are consistent with the administrative division of adjacent slope units are selected, and the edge with the largest cumulative flow is selected as the edge to be dissolved and stored in flush_line.

[0057] Step 5: Input flush_line to find the slope cells that can be continuously dissolved, reducing the number of iterations. left represents the slope cell to the left of the line segment, and right represents the right.

[0058] (5.1) Initialize is_fush to -1, indicating that it has not been modified;

[0059] (5.2) Traverse the fush_line, skipping the line segments where is_fush is 1;

[0060] (5.3) Add the unmodified line segments from the left neighbor of the line segment to the queue;

[0061] (5.4) Dequeue the line segment in the queue, set the fusion flag to the id of left, set is_fush to 1, and add the neighboring line segments of the line segment to the queue until the queue is empty;

[0062] (5.5) After traversal, update flush_line;

[0063] Updating `fush_line` allows some consecutive slope cells to dissolve simultaneously, reducing the number of dissolution iterations and increasing the dissolution rate.

[0064] Step 6: Set the dissolution markers of the left and right ramp cells of the line segment in `fush_line` to the dissolution markers of that line segment. Then dissolve the ramp cells with the same dissolution markers and calculate the density of the dissolved ramp cells. If the density meets the target density or reaches the maximum number of iterations, output the ramp cells; otherwise, return to Step 3, extract the features of the dissolved ramp cells, and dissolve them again. Note that to improve the calculation speed, the intersection features that have not changed do not need to be extracted again.

[0065] Figure 4 shows an example of the slope element obtained by this invention. It is the result of slope element extraction using the DEM of a town in Huangshan after 15 iterations. The top left image shows the result of overlaying the slope element with the slope aspect of the area, the top right image shows the result with contour lines overlaid, the bottom left image shows the result without overlay, and the bottom right image shows the result with key areas such as roads and houses overlaid. Since there were no actual key areas, these key areas were set by the experiment.

[0066] This invention designs the following slope unit extraction criteria: using the watershed algorithm with a relatively small initial threshold, initial coarse slope units are obtained; these coarse slope units are then dissolved according to rules to obtain better slope units. Compared to existing methods, this invention employs an extraction-then-dissolution approach, which can maximize the number of slope units with reasonable geographical features. Furthermore, this invention incorporates confluence accumulation into the dissolution process, ensuring that the slope units do not disrupt ridgelines or valley lines.

[0067] The above are preferred embodiments of the present invention, but the scope of protection of the present invention is not limited thereto. Any modifications or substitutions conceived by those skilled in the art within the scope of the technology disclosed in the present invention without inventive effort are included within the scope of protection of the present invention. The scope of protection of the present invention is determined by the scope defined in the claims.

Claims

1. An adaptive slope element extraction method based on the watershed mechanism, comprising the following steps: Step 1: First, preprocess the digital elevation map to generate a reverse elevation map. Then, fill in depressions in both elevation maps. Using the Grass GIS Python library functions, calculate the aspect raster, slope raster, and flow direction raster. Calculate the cosine and sinine values ​​of the aspect (i.e., the horizontal component sinno and the vertical component cosno) using a raster calculator. These values ​​are used for subsequent aspect consistency calculations. The confluence accumulation algorithm is used to calculate the confluence accumulation values ​​su_acc and re_acc for both the forward and reverse elevation maps. The confluence accumulation algorithm is as follows: Input flow to raster map dir (1.1) Initialize the inflow matrix NIDP as a zero matrix of the same size as dir, and set the nodata part of dir to nodata in NIDP as well; (1.2) Using the NIDP-base algorithm, traverse the dir to calculate NIDP. For each cell that flows into the traversed cell, the value of that cell is incremented by 1. (1.3) Initialize the flow accumulator matrix FlowAccu, a matrix of all ones with the same size as dir; (1.4) Traverse the NIDP matrix, find the cell with an NIDP value of 0 as the starting point, and according to the flow direction of the dir matrix, flow to the position in NIDP that is not 1, and decrement the NIDP value at that position by one. Add the value of the outflow position to the value of the FlowAccu matrix. (1.5) Export FlowAccu tif Step 2: Using the watershed watershed algorithm in Grass GIS, a minimum area threshold is set. This minimum area threshold is optimized using a bisection method and a set slope unit density. Initial coarse slope units (SLU) are extracted. Flat slope areas with a slope less than 5 degrees (mask) are extracted, and the fragmented parts of the mask are removed. The mask parts in the SLU are also trimmed. If there is an administrative boundary vector map, it is overlaid on the SLU. Each slope unit is numbered according to its administrative division. If there is no administrative boundary vector map, the administrative division number of each slope unit is set to 0. The fragmented slope units caused by calculation errors are dissolved. Finally, the slope unit boundary lines are smoothed to obtain the initial version of the slope unit. The slope unit density of this version is slightly larger than the target density to facilitate subsequent refinement and dissolution. Step 3: Extract slope unit features and line segment features. First, extract the surface features corresponding to the slope unit based on the feature raster map in Step 1. The surface features are area, slope aspect or administrative division. Then, construct the adjacency relationship table of the slope unit to extract the intersection line of the adjacent unit. Extract the confluence accumulation of the intersection line through su_acc and re_acc. Quantize the values ​​in su_acc and re_acc and take the larger value of the two to obtain the confluence accumulation feature of the intersection line. Step 4: Traverse the slope cells. Based on the area features obtained in Step 3, determine whether the slope cells need to be dissolved. Calculate the slope aspect consistency of each edge of the slope cells that need to be dissolved, and select the edges that meet the slope aspect consistency threshold. The slope aspect consistency threshold is calculated as follows: in It represents the threshold value for a specific slope unit, where k is a threshold control coefficient; the smaller the value, the easier it is to dissolve. The result is the area of ​​a certain slope unit standardized; then, from these edges, the edges that are consistent with the administrative division of the adjacent slope units are selected, and the edge with the largest cumulative flow is selected as the edge to be dissolved and stored in flush_line; Step 5: Input fush_line to find the slope cell that is continuously dissolved, reduce the number of iterations. left represents the slope cell to the left of the line segment, and right represents the right. (5.1) Initialize is_fush to -1, indicating that it has not been modified; (5.2) Traverse the fush_line, skipping the line segments where is_fush is 1; (5.3) Add the unmodified line segments from the left neighbor of the line segment to the queue; (5.4) Dequeue the line segment in the queue, set the fusion flag to the id of left, set is_fush to 1, and add the neighboring line segments of the line segment to the queue until the queue is empty; (5.5) After traversal, update flush_line; Updating fush_line allows some consecutive ramp cells to dissolve simultaneously, reducing the number of dissolution iterations and increasing the dissolution rate. Step 6: Set the dissolution markers of the left and right ramp cells of the line segment in `fush_line` to the dissolution markers of that line segment. Then dissolve the ramp cells with the same dissolution markers and calculate the density of the dissolved ramp cells. If the density meets the target density or reaches the maximum number of iterations, output the ramp cells; otherwise, return to Step 3, extract the features of the dissolved ramp cells, and dissolve them again. Note that to improve the calculation speed, the intersection features that have not changed do not need to be extracted again.

Images