A method for automatic extraction of gully networks based on surface concavity index

Abstract

The present invention discloses a method for automatically extracting a gully network based on a surface concavity index, and relates to the field of terrain analysis. The key points of the technical solution include the following steps: S1, calculating potential gully network data based on DEM data; S2, calculating the surface concavity index based on DEM data; S3, obtaining the gully starting position information by integrating the potential gully network data and the surface concavity index data; S4, starting from the gully starting position, by finding the grid position of the cumulative amount of convergence in the neighborhood window, the generation of the gully network in the study area is realized. The surface concavity index proposed by the present invention can accurately describe the concave morphological characteristics of the terrain, has a clear definition and is easy to calculate based on grid DEM data, can be used to assist in the accurate extraction of the gully network, and can also provide technical references and method references for scientific research fields such as digital terrain analysis.

Description

Technical Field

[0001] The invention relates to the field of terrain analysis, and in particular to an automatic extraction method for a gully network based on a surface concavity index. Background Art

[0002] A gully is a narrow depression formed by water erosion on the earth's surface. In areas with severe water erosion, gullies of different levels are inlaid and interconnected to form a gully network in the area. In areas with severe water erosion, such as the Loess Plateau in my country, accurate extraction of gully networks is of great significance for agricultural production planning, soil and water conservation, sand control, and other scientific research, such as landform development and evolution.

[0003] The main method of extracting valley networks at present is to simulate the flow of surface runoff on the surface based on DEM data using the surface runoff overland flow model to generate a water system, that is, a valley network. The basic process is to first fill the DEM data to eliminate local depressions in the DEM data, and then determine the direction of water flow on each grid unit based on the principle of maximum elevation difference, thereby simulating the confluence route of water flow on the surface of a region; then, based on this confluence route, calculate the number of water catchment units passed by each grid unit, and the total number obtained by statistics is called the confluence accumulation of the grid unit; when the confluence accumulation of a grid reaches a certain numerical range, it can be considered that the grid becomes part of the water system or valley, that is, the location of the valley is determined by judging the confluence accumulation value of the grid. This method discriminates the flow path and simulates surface runoff based on the concept of hydrological confluence. It has a certain model basis, and the algorithm is simple, and it can directly generate a continuous flow path. It is currently a commonly used valley network extraction method.

[0004] A key step in the above method is to determine a runoff accumulation threshold, which determines the starting position of a valley water system and the number and distribution of valley networks in a region. Therefore, the main research on valley water system extraction focuses on how to determine this runoff accumulation threshold. Comprehensive analysis of existing research shows that there are two methods that are currently used more frequently: the first is the experimental method, which uses different thresholds to extract the valley network, and then overlays and compares with existing water system maps, remote sensing images and other data to determine the most suitable runoff accumulation threshold; the second is the method of finding mutation points. The main idea is to calculate the valley network under different runoff accumulation thresholds, then count its total length, draw a scatter plot between different threshold parameters and the total length of the valley network, fit the scatter plot into a curve, and find the mutation point of the curve as the optimal threshold. The typical method for finding mutation points is the mean change point analysis method.

[0005] The above two methods have obvious defects: first, the workload is large. Both methods need to set different threshold parameters to extract the valley network. At the same time, the subsequent result comparison and mutation point calculation are also time-consuming. Second, all methods currently try to determine an optimal runoff accumulation threshold when extracting the valley network, but ignore the variability of the threshold. Even in the same area, the degree of erosion of different surfaces is different, and the runoff accumulation threshold parameter corresponding to the starting position of the valley may be different. Therefore, if a threshold parameter is fixed, even if an appropriate parameter value is determined, the extracted valley network will have a large error.

[0006] In order to solve the above technical problems, the present invention starts from the definition of valleys, utilizes the "narrow depression" characteristics presented on the surface, and constructs a surface concavity index calculated based on DEM data. The index first determines the starting position of the valley on the regional surface, and then combines the surface runoff overflow model to follow the confluence path from the starting position of the valley to the watershed outlet. This method no longer fixes a single convergence accumulation threshold parameter in the process of extracting the valley network, but first finds and determines the starting position of the valley from the perspective of surface morphology, and then combines the surface runoff overflow model to obtain a continuous and complete valley network. The method of the present invention is easy to operate, and the extracted valley network is more consistent with the valley network on the actual surface. It has important application value in practical projects such as agricultural planning and soil and water conservation planning and in geomorphology-related scientific research fields. Summary of the invention

[0007] In view of the above problems, the present invention designs a surface concavity index that can effectively reflect the concave characteristics of the surface, and realizes the calculation of the index based on regular grid DEM data, and then uses the index to realize the rapid and accurate extraction of the valley network. The method includes two key technical links: the calculation of the surface concavity index based on DEM data and the extraction of the valley network assisted by the surface concavity index, which can effectively solve the technical problem that the threshold parameter of the accumulated flow is difficult to effectively determine in the current valley network extraction.

[0008] The technical solution adopted by the present invention is as follows:

[0009] The automatic extraction method of valley network based on surface concavity index includes the following steps:

[0010] S1, input DEM data, calculate the runoff accumulation value of each grid through depression filling, flow direction calculation, runoff accumulation calculation and other processes, and reclassify the runoff accumulation data according to whether the runoff accumulation value is greater than the specified parameter value. If it is greater, the value is 1, otherwise it is 0, and the first binary grid data of the study area is obtained;

[0011] S2, based on the DEM data, calculate the terrain concavity index of the study area, and reclassify the concavity index according to whether the terrain concavity index is greater than the specified parameter value. If it is greater, the value is 1, otherwise it is 0, and the second binary raster data of the study area is obtained;

[0012] S3, performing a superposition multiplication operation on the first binary data of the study area obtained in step S1 and the second binary grid data obtained in step S2 to obtain the third binary data of the study area;

[0013] S4, converting the binary data obtained in step S3, selecting the grid with a value of 1, into a vector line format, and then extracting the two endpoints of each vector line;

[0014] S5, using the endpoints of the vector lines obtained in step S4, obtain the cumulative flow value on the endpoints, and through window analysis, search for grids with increased cumulative flow values ​​until the watershed outlet, and finally obtain the valley network of the study area.

[0015] Furthermore, the specific process of step S2 is as follows:

[0016] S21, for the DEM data of the study area, read the basic information of the DEM data such as grid resolution h, number of rows m, number of columns n, etc.;

[0017] S22, starting from the first grid, determine the position of the current grid, if the current grid is located on the boundary, that is, its row number is 1 or m, or the column number is 1 or n, then assign 0 to the surface concavity index of the grid position, otherwise, execute step S23;

[0018] S23, if the current grid is not on the boundary, a 3×3 grid window is formed with the current grid as the center and nine grids including the current grid on its left, right, top, bottom, upper left, lower right, upper right, and lower left;

[0019] S24, according to the spatial position relationship of different grids in the window and the grid resolution h, respectively calculate the cosine value of the spatial angle between the central grid and its left and right grids, the central grid and its upper and lower grids, the central grid and its upper left and lower right grids, and the central grid and its upper right and lower left grids;

[0020] S25, comparing the four cosine values ​​calculated in step S24, and adding 1 to the maximum value as the surface concavity index of the central grid;

[0021] S26, repeating steps S22-S25, traversing all grids, and completing the calculation of the surface concavity index of each grid;

[0022] S27, setting a threshold parameter, judging the surface concavity index calculated in steps S21-S26, if the surface concavity index of a grid is greater than the threshold parameter, assigning the grid a value of 1, otherwise assigning the grid a value of 0, and obtaining binary data.

[0023] Furthermore, the specific process of step S5 is as follows:

[0024] S51, define a grid data B that is consistent with the basic information of the research data grid resolution, number of rows and columns, and assign an initial value of 0;

[0025] S52, read the point data obtained in step S4, obtain the cumulative value of the confluence at the spatial position corresponding to each point, and then execute steps S53-S55 in sequence for each point;

[0026] S53, reading a point, determining the grid cells at the same position as the point in the grid data A of the accumulated flow and the grid data B defined in step S51, and marking the corresponding grid cells in the grid data B as 1;

[0027] S54, taking the grid cell in the grid data A determined in step S53 as the current processing grid, searching for grid cells with a convergence accumulation value greater than the convergence accumulation value of the current processing grid cell in the surrounding grids (i.e., the left, right, top, bottom, upper left, upper right, lower left and lower right) thereof, if there are multiple grid cells, selecting the largest grid cell as the current processing grid cell, and marking the corresponding grid cell in the grid data B as 1;

[0028] S55, repeating step S54, during which the grid cells in the grid data B can be repeatedly marked until no grid cell with a flow accumulation value greater than the flow accumulation value of the currently processed grid cell is found in the grid data A, indicating that the watershed outlet has been reached, then other endpoints that have not been processed are selected, and steps S53-S55 are repeated;

[0029] S56, after all points are processed, the grids marked as 1 in the grid data B constitute the valley network of the study area, which is converted into real data output to obtain the valley network data of the study area.

[0030] Technical features and beneficial effects of the present invention:

[0031] (1) The surface concavity index proposed in the present invention can accurately describe the concave morphological characteristics of the terrain. The index is clearly defined and easy to calculate based on grid DEM data. It can provide technical reference and methodological reference for scientific research fields such as surface morphological analysis and digital terrain analysis.

[0032] (2) The valley network extraction process assisted by the surface concavity index proposed in the present invention combines the basic analysis process in the classical technical process on the one hand, and determines the starting position of the valley by calculating the concave shape of the surface on the other hand, thus solving the problem that the threshold parameter of the runoff accumulation amount is difficult to determine in the classical method.

[0033] (3) The valley network extraction process proposed in the present invention, which is assisted by the surface concavity index, no longer uses a fixed runoff accumulation threshold parameter as the only criterion for determining valleys. Instead, it incorporates the morphological characteristics of the valleys, significantly improving the accuracy of valley network extraction in areas with complex terrain changes. BRIEF DESCRIPTION OF THE DRAWINGS

[0034] Figure 1 It is a schematic diagram of the process of the present invention.

[0035] Figure 2 Schematic diagram of the local rectangular coordinate system constructed between grid cells in the method of the present invention.

[0036] Figure 3 Schematic diagram of vectors composed of center points of grid cells in the method of the present invention.

[0037] Figure 4 Schematic diagram of the process of finding the valley grid in the method of the present invention.

[0038] Figure 5 Schematic diagram of valley grid unit marking in the method of the present invention.

[0039] Figure 6 Schematic diagram of the valley generated in the method of the present invention.

[0040] Figure 7 Schematic diagram of the surface concavity index calculated in the experimental area in an embodiment of the present invention.

[0041] Figure 8 This is a schematic diagram of the candidate starting points of the valley calculated in the experimental area in an embodiment of the present invention.

[0042] Fig. 9 This is a schematic diagram of the valley network extracted from the experimental area in an embodiment of the present invention.

[0043] Fig.10 This is a schematic diagram of the valley starting points extracted from the valley network in the experimental area in an embodiment of the present invention.

[0044] Fig.11 This is a schematic diagram of the valley network extracted by the experimental area comparison method in an embodiment of the present invention. DETAILED DESCRIPTION

[0045] The principles and features of the present invention are described below in conjunction with the accompanying drawings. The examples given are only used to explain the present invention and are not used to limit the scope of the present invention.

[0046] like Figure 1 As shown, the specific steps of the method for automatically extracting valley networks based on the surface concavity index of the present invention are as follows:

[0047] Step 1: Read the DEM data of the study area, and implement the processes of depression filling, flow direction calculation, and runoff accumulation calculation in sequence to obtain the runoff accumulation value of each grid;

[0048] Step 2: Take 10 as the parameter, and judge the grid data of the runoff accumulation obtained in the step. If the runoff accumulation value of the grid unit is greater than 10, the grid is assigned a value of 1, otherwise it is assigned a value of 0, and the first binary grid data is obtained;

[0049] Step 3: Read the DEM data of the study area and calculate the surface concavity index of each grid cell in turn. If the grid to be calculated is located on the boundary, that is, it is in the first row, or the last row, or the first column, or the last column, the surface concavity index of the grid is assigned to 0, otherwise it is calculated according to step 4;

[0050] Step 4: For grid cells at non-boundary positions, the current processing grid is taken as the center, and the surrounding 8 grids form a 3×3 analysis window. The center point of the central grid is taken as the origin, and a local rectangular coordinate system is established in the horizontal and vertical directions passing through the origin (see Figure 2 );

[0051] Step 5: Based on the local rectangular coordinate system constructed in step 4, with the origin of the coordinate system as the center point, generate vectors of the center point and the eight points around it (see Figure 3 ), the figure is for displaying vectors, and the two coordinate axes are omitted. Then the cosine values ​​of the angles between vectors P1O and P8O, P2O and P7O, P3O and P6O, and P4O and P5O are calculated respectively, and the calculation formula is:

[0052]

[0053] Where θ represents the angle between two vectors, (x1, y1, z1) and (x2, y2, z2) are the spatial coordinates of the two vectors.

[0054] After the calculation is completed, take the maximum value of the cosine values ​​of the four angles and add 1 as the surface concavity index value of the current processing grid cell, that is:

[0055] Con=max(cosθ1+cosθ2+cosθ3+cosθ4)+1 (Formula 2)

[0056] Step 6: Repeat steps 4-5 to complete the calculation of the surface concavity index of the study area;

[0057] Step 7: For the calculated surface concavity index in the study area, 0.8 is used as the critical parameter, and the surface concavity index greater than 0.8 is assigned a value of 1, and the surface concavity index less than or equal to 0.8 is assigned a value of 0, to obtain the second binary raster data;

[0058] Step 8: Perform multiplication and superposition operation on two binary raster data. If the values ​​of the two binary raster data are both 1, the result of the superposition operation is 1, otherwise it is 0, and the third binary raster data is obtained;

[0059] Step 9: Select the grid with a value of 1 in the third binary grid data, and convert the grids at both ends of the continuously distributed grids into vector point data;

[0060] Step 10: For the two end point data obtained from each continuous grid in step 9, respectively obtain the cumulative flow value at the point position, and delete the point with the larger cumulative value;

[0061] Step 11: Generate result grid data and assign the initial value to 0;

[0062] Step 12: for the point data obtained in step 10, process each point in turn, calculate the position of the point in the grid data, determine the corresponding grid cells in the confluence accumulation grid data and the result grid data, and assign the result grid cell a value of 1;

[0063] Step 13: With the current processing grid determined in the runoff accumulation grid data as the center, search for grids around it whose runoff accumulation values ​​are greater than the runoff accumulation values ​​of the current grid, and use them as the new current processing grid, and mark the grid data corresponding to the position in the result grid data as 1 (see Figure 4 , Figure 5 );

[0064] Step 14: Repeat step 13 until no grid with a flow accumulation value greater than the current grid flow accumulation value is found, indicating that the watershed outlet has been reached and a flow path search is completed;

[0065] Step 15: Repeat steps 12-14. After all point data are processed, the points marked as 1 in the resulting grid data are the locations of the valleys (see Figure 6 );

[0066] Step 16: Convert the continuous grid with a median value of 1 in the result data into vector lines to obtain the valley network of the study area and output it.

[0067] Example

[0068] In order to verify the feasibility of the method of the present invention in extracting gully networks, the embodiment selects a small watershed in northern Shaanxi where the loess is higher than a certain area as the experimental area. The experimental area is a typical water erosion terrain with a very well-developed gully network system, which can be used as an ideal experimental area for differential analysis of gully network extraction results using different methods.

[0069] Step 1: For the sample area, after the processes of depression filling, flow direction calculation, and runoff accumulation calculation, the calculated runoff accumulation takes 10 as the critical parameter, and a value of 1 is assigned if it is greater than 10, otherwise it is assigned to 0, and the first binary raster data of the sample area is obtained;

[0070] Step 2: Calculate the surface concavity index for the sample area (see the results in Figure 7 ), and then take 0.8 as the critical parameter, assign a value of 1 if it is greater than 0.8, otherwise assign a value of 0, and obtain the second binary raster data of the sample area;

[0071] Step 3: Multiply the binary raster data obtained in step 1 and the binary raster data obtained in step 2 to obtain the third binary raster data of the sample area;

[0072] Step 4: Convert the binary raster data obtained in step 3 into vector lines, then take the two end points of the line segment, obtain the cumulative flow value on the endpoints, and delete the endpoint with the larger cumulative flow value from the same line segment to obtain the available point data ( Figure 8 );

[0073] Step 5: Based on the point data obtained in step 4, taking the point as the starting position, the water flow path line is obtained by searching for adjacent grids with increased water flow accumulation values;

[0074] Step 6: Convert the grid on the searched catchment path line into vector data to obtain the valley network result of the sample area ( Fig. 9 ), and then extract the valley starting point of the sample area based on the valley network ( Fig.10 );

[0075] Step 7: For the sample area, the valley network results obtained by the contrast method (the method of determining the threshold by the mean change point analysis method) are obtained (see Fig.11 ).

[0076] Through the above comparative analysis, there are two obvious problems in the valleys extracted by the runoff accumulation threshold calculated by the mean change point analysis method: the first is that a large number of pseudo valleys appear, mainly in the upstream area of ​​the small watershed, where the gullies are not well developed, there are many flat terrains, and pseudo valleys are easily formed; the second problem is that many short valleys are not extracted, because the catchment area of ​​these valleys is small, far from reaching the determined runoff accumulation threshold, so they cannot be effectively extracted. The first problem of the valleys extracted by the method of the present invention does not exist at all, because pseudo valleys are mainly distributed in flat terrain areas, and the surface concavity index of these areas is approximately 0, so the valleys extracted by the method of the present invention can effectively avoid the appearance of pseudo valleys; the second problem also exists partially for the method of the present invention, mainly because some valleys are not fully developed, or are in the early stage of valley development, and the concave shape of the surface is not obvious, so the method of the present invention cannot be effectively extracted, but the number of valleys that cannot be identified is significantly less than the runoff accumulation threshold method. Therefore, based on the above analysis, the valley network extraction result of the method of the present invention is more accurate.

Claims

1. A method for automatically extracting gully networks based on surface concavity index, characterized in that: The steps include: S1, input DEM data, calculate the runoff accumulation value of each grid through depression filling, flow direction calculation, and runoff accumulation calculation process, and reclassify the runoff accumulation data according to whether the runoff accumulation value is greater than the specified parameter value. If it is greater, the value is 1, otherwise it is 0, and the first binary grid data of the study area is obtained; S2, based on the DEM data, calculate the terrain concavity index of the study area, and reclassify the concavity index according to whether the terrain concavity index is greater than the specified parameter value. If it is greater, the value is 1, otherwise it is 0, and the second binary raster data of the study area is obtained; S3, performing a superposition multiplication operation on the first binary data of the study area obtained in step S1 and the second binary grid data obtained in step S2 to obtain the third binary data of the study area; S4, converting the binary data obtained in step S3, selecting the grid with a value of 1, into a vector line format, and then extracting the two endpoints of each vector line; S5, using the endpoints of the vector line obtained in step S4, obtaining the cumulative flow value on the endpoints, searching for grids with increased cumulative flow value through window analysis until the watershed outlet, and finally obtaining the gully network of the study area; The specific process of step S2 is: S21, for the DEM data of the study area, read the basic information of the DEM data grid resolution h, number of rows m, number of columns n; S22, starting from the first grid, determine the position of the current grid, if the current grid is located on the boundary, that is, its row number is 1 or m, or the column number is 1 or n, then assign 0 to the surface concavity index of the grid position, otherwise, execute step S23; S23, if the current grid is not on the boundary, a 3×3 grid window is formed with the current grid as the center and nine grids including the current grid on its left, right, top, bottom, upper left, lower right, upper right, and lower left; S24, according to the spatial position relationship of different grids in the window and the grid resolution h, respectively calculate the cosine value of the spatial angle between the central grid and its left and right grids, the central grid and its upper and lower grids, the central grid and its upper left and lower right grids, and the central grid and its upper right and lower left grids; S25, comparing the four cosine values ​​calculated in step S24, and adding 1 to the maximum value as the surface concavity index of the central grid; S26, repeating steps S22-S25, traversing all grids, and completing the calculation of the surface concavity index of each grid; S27, setting a threshold parameter, judging the surface concavity index calculated in steps S21-S26, if the surface concavity index of a grid is greater than the threshold parameter, assigning the grid a value of 1, otherwise assigning the grid a value of 0, and obtaining binary data.

2. The method for automatically extracting gully networks based on surface concavity index according to claim 1, characterized in that: The specific process of step S5 is as follows: S51, define a grid data B that is consistent with the basic information of the research data grid resolution and number of rows and columns, and assign an initial value of 0; S52, read the point data obtained in step S4, obtain the cumulative value of the confluence at the spatial position corresponding to each point, and then execute steps S53-S55 in sequence for each point; S53, reading a point, determining the grid cells at the same position as the point in the grid data A of the accumulated flow and the grid data B defined in step S51, and marking the corresponding grid cells in the grid data B as 1; S54, taking the grid cell in the grid data A determined in step S53 as the current processing grid, searching the surrounding grids for grid cells whose confluence accumulation value is greater than the confluence accumulation value of the current processing grid cell, if there are multiple grid cells, selecting the largest grid cell as the current processing grid cell, and marking the corresponding grid cell in the grid data B as 1; S55, repeating step S54, during which the grid cells in the grid data B can be repeatedly marked until no grid cell with a flow accumulation value greater than the flow accumulation value of the currently processed grid cell is found in the grid data A, indicating that the watershed outlet has been reached, then other endpoints that have not been processed are selected, and steps S53-S55 are repeated; S56, after all points are processed, the grids marked as 1 in the grid data B constitute the valley network of the study area, which are converted into vector data and output to obtain the valley network data of the study area.

Images