River network center line and catchment area automatic extraction method based on river network density maintenance

Abstract

The invention discloses a river network center line and catchment area automatic extraction method based on river network density maintenance, and belongs to the field of digital hydrological modeling and geographic information, and the method comprises the steps: taking two-dimensional river network density distribution as spatial variable input, combining flow direction data, achieving the collaborative extraction of a river network center line and a catchment area, and obtaining a river network center line and catchment area extraction result; and keeping the density of the extracted river network in each local watershed unit consistent with the input density field. The method provides high-precision river network data support for water regimen simulation, water and soil loss evaluation and intelligent water conservancy construction.

Description

Technical Field

[0001] The present invention belongs to the fields of digital hydrological modeling and geographic information, and particularly relates to a method for automatically extracting river network centerlines and catchment areas while maintaining river network density. Background Art

[0002] River network extraction based on digital elevation models is a key link in hydrological simulation, and its accuracy directly affects the simulation effect of the basin hydrological process. Currently, the mainstream river network extraction methods generally follow the following technical processes: (1) calculating the flow direction based on DEM data; (2) calculating the cumulative runoff or the cumulative runoff area; (3) setting a fixed threshold to extract the river channel centerline; (4) dividing the catchment area. Among them, the implementation of steps (3) and (4) both rely on a unified cumulative runoff threshold set artificially.

[0003] The prior art has the following significant technical defects: Using a fixed threshold will result in the river network density (the total length of river channels per unit area) extracted being inconsistent with the actual situation. This is because there is significant spatial heterogeneity in the river network density within a natural basin, which is affected by multiple factors such as topography, geology, and climate, and the fixed threshold cannot reflect this spatial variation characteristic.

[0004] To solve the above problems, existing improved methods adopt a multi-level threshold strategy. However, the multi-level threshold method requires multiple iterative calculations, significantly increasing the consumption of computing resources. Especially when applied at the large basin scale, the computing cost increases exponentially and cannot be achieved. Summary of the Invention

[0005] To solve the above technical problems, the present invention provides a method for automatically extracting river network centerlines and catchment areas while maintaining river network density, including:

[0006] Calculating the cumulative runoff length according to the flow direction and river network density;

[0007] Calculating the cumulative runoff area according to the flow direction;

[0008] Extracting all river reaches and their centerlines according to the flow direction, cumulative runoff length, and cumulative runoff area;

[0009] Extracting the catchment area according to the flow direction and centerline.

[0010] Preferably, the process of calculating the cumulative runoff length according to the flow direction and river network density includes:

[0011] Initializing the number of inflow grids of each grid to zero;

[0012] Traversing the flow direction grids and adding one to the number of inflow grids of its downstream grids;

[0013] Initializing the cumulative runoff length of each grid to be equal to the grid area multiplied by the river network density;

[0014] Traverse all grids, find the grids with zero inflow grids, and use them as the river sources;

[0015] Select one of the river source grids in turn, set it as the current grid, search downstream according to the flow direction, and calculate the cumulative length of confluence.

[0016] Preferably, the process of searching downstream according to the flow direction and calculating the cumulative length of confluence includes:

[0017] Find the downstream grid of the current grid according to the flow direction. If there is no downstream for the current grid, the search ends;

[0018] If there is a downstream grid, update the cumulative length of confluence of the downstream grid, and its value is equal to the sum of the cumulative length of confluence of the downstream grid and the current downstream grid;

[0019] If the number of inflow grids of the downstream grid is greater than or equal to two, subtract one from the number of inflow grids of the downstream grid, and the search ends. Otherwise, set the downstream grid as the current grid, and re-find the downstream grid of the current grid according to the flow direction for judgment.

[0020] Preferably, the process of calculating the cumulative confluence area according to the flow direction includes:

[0021] Initialize the number of inflow grids of each grid to zero;

[0022] Traverse the flow direction grids and add one to the number of inflow grids of their downstream grids;

[0023] Initialize the cumulative confluence area of each grid to be equal to its grid area;

[0024] Traverse all grids, find the grids with zero inflow grids, and use them as the river sources;

[0025] Select one of the river source grids in turn, set it as the current grid, search downstream according to the flow direction, and calculate the cumulative confluence area.

[0026] Preferably, the process of searching downstream according to the flow direction and calculating the cumulative confluence area includes:

[0027] Find the downstream grid of the current grid according to the flow direction. If there is no downstream for the current grid or the downstream is outside the search area, the search ends;

[0028] If there is a downstream grid, update the cumulative confluence area of the downstream grid, and its value is equal to the sum of the cumulative confluence area of the downstream grid and the current downstream grid;

[0029] If the number of inflow grids of the downstream grid is greater than or equal to two, subtract one from the number of inflow grids of the downstream grid and terminate the search; otherwise, set the downstream grid as the current grid, and re-find the downstream grid of the current grid according to the flow direction for judgment.

[0030] Preferably, the process of extracting all river reaches and their centerlines according to the flow direction, cumulative length of confluence, and cumulative area of confluence includes:

[0031] Initialize an empty centerline temporary queue and an empty centerline completed queue. Both queues store the centerlines of river reaches, and each centerline of a river reach is represented by a grid sequence, and the front-back order of the grids represents the upstream-downstream relationship;

[0032] Traverse the flow direction grids. If a grid has no downstream grid or the downstream is outside the search area, then the grid is the outlet of the basin. Take each outlet grid of the basin as a centerline grid sequence with a length of one and save it to the centerline temporary queue;

[0033] Select the centerline of the river reach at the front end of the centerline temporary queue and delete it from the queue;

[0034] Calculate the selected centerline of the river reach to obtain the data of all river reaches and their centerlines.

[0035] Preferably, the process of calculating the selected centerline of the river reach to obtain the data of all river reaches and their centerlines includes:

[0036] Calculate the target length of the centerline of the river reach. If the centerline of the river reach has no upstream, that is, there is no centerline in the temporary queue and the completed queue flowing to this centerline, then the target length is equal to the cumulative length of confluence at the location of the most downstream grid of the centerline of the river reach; otherwise, the target length is equal to the cumulative length of confluence at the location of the most downstream grid of the centerline of the river reach minus the sum of the cumulative lengths of confluence at the locations of the most downstream grids of all its upstream centerlines. If the length of the centerline is greater than or equal to the target length, then this centerline is saved to the completed queue and the loop is terminated;

[0037] If the length of the centerline is less than the target length, then search for the set of upstream grids of this centerline according to the flow direction;

[0038] Compare the cumulative areas of confluence based on the set of upstream grids, and find the grid with the largest cumulative area of confluence as the grid to be determined;

[0039] If the grid to be determined flows to the most upstream grid in the centerline, then add the grid to be determined to the front of the river channel centerline sequence, re-calculate the target length of the centerline of the river reach and make a judgment;

[0040] If the grid to be determined does not flow towards the most upstream grid in the center line, find the downstream grid of the grid to be determined according to the flow direction. Taking the downstream grid of the grid to be determined as the boundary, divide the center line into two sections: the upstream section and the downstream section. The upstream section does not include the downstream grid of the grid to be determined, and the downstream section includes the downstream grid of the grid to be determined. Save both the upstream and downstream sections of the center line to the end of the center line temporary queue.

[0041] Take the grid to be determined as a river reach center line grid sequence with a length of one, and save it to the end of the center line temporary queue.

[0042] Take a river reach from the front end of the temporary queue again, delete it from the queue, calculate the target length of the center line of this river reach and make a judgment; repeat the judgment process until the center line temporary queue is empty. At this time, the center line complete queue saves the center lines of all river reaches in the river network.

[0043] Preferably, the process of extracting the catchment area according to the flow direction and the center line includes:

[0044] Establish an empty exclusion set. If there is an upstream center line of the center line, traverse the upstream center line of the center line, and add the most downstream grid of the upstream center line to the exclusion set.

[0045] Initialize an empty catchment area temporary queue and an empty catchment area complete set. Both the temporary queue and the complete set save the grids covered by the catchment area.

[0046] Push the grids in the center line sequence into the end of the temporary queue in the order from downstream to upstream.

[0047] If the temporary queue is not empty, take a grid from the front end of the temporary queue, save it to the complete set, and delete it from the temporary queue. According to the flow direction, find all the grids upstream of the taken grid and not included in the complete set or the exclusion set, and push these grids into the end of the temporary queue. Repeat until the temporary queue is empty to complete the extraction of the catchment area.

[0048] Preferably, the river network center line is represented by the grid where the river network center line is located, and the river network center line vector is obtained by vectorizing the grid where the river network center line is located.

[0049] Compared with the prior art, the present invention has the following advantages and technical effects:

[0050] In view of the problem of river network density distortion caused by the global unified threshold in the existing river network extraction methods, the present invention proposes an automated extraction algorithm based on dynamic river network density constraints. By taking the two-dimensional river network density distribution as a spatial variable input and combining the flow direction data, the collaborative extraction of the river network centerline and catchment area is realized, and the density of the extracted river network within each local watershed unit is kept consistent with the input density field. This method provides high-precision river network data support for flood situation simulation, soil erosion assessment, and the construction of intelligent water conservancy. BRIEF DESCRIPTION OF THE DRAWINGS

[0051] The drawings constituting a part of this application are used to provide a further understanding of this application. The schematic embodiments of this application and their descriptions are used to explain this application and do not constitute an improper limitation to this application. In the drawings:

[0052] Figure 1 It is a schematic flowchart of the method according to an embodiment of the present invention. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0053] It should be noted that, without conflict, the embodiments in this application and the features in the embodiments can be combined with each other. The following will refer to the drawings and combine the embodiments to detail this application.

[0054] It should be noted that the steps shown in the flowchart of the drawings can be executed in a computer system such as a set of computer-executable instructions, and although the logical order is shown in the flowchart, in some cases, the steps shown or described can be executed in a different order than here.

[0055] Embodiment 1

[0056] As Figure 1 shown, this embodiment provides an automated extraction method for the river network centerline and catchment area with river network density maintained, including:

[0057] 1. Calculate the cumulative length of flow concentration according to the flow direction and river network density.

[0058] The cumulative length of flow concentration is defined as the total length of all upstream main and tributary rivers of a certain grid cell (or river network node), that is, the sum of the lengths of the water flow paths flowing towards this grid.

[0059] The cumulative length of flow concentration can be calculated through the following steps:

[0060] (1) Calculate the number of inflowing grid cells for each grid. 1. Initialize the number of inflowing grid cells for each grid to zero.

[0061] 2. Traverse the flow direction grid and add one to the number of inflowing grid cells of its downstream grid.

[0062] (2) Calculate the accumulated length of flow concentration. 1. Initialize the accumulated length of flow concentration for each grid to be equal to the grid area multiplied by the river network density. 2. Traverse all grids to find the grids with zero number of inflowing grids, which are the river sources. 3. Select one of the river source grids in turn and set it as the current grid, and search downstream according to the flow direction: (a) Find the downstream grid of the current grid according to the flow direction. If there is no downstream for the current grid, the search terminates; (b) If there is a downstream grid, update the accumulated length of flow concentration of the downstream grid, and its value is equal to the sum of the accumulated length of flow concentration of the downstream grid and the current downstream grid; (c) If the number of inflowing grids of the downstream grid is greater than or equal to two, subtract one from the number of inflowing grids of the downstream grid and the search terminates; otherwise, set the downstream grid as the current grid and repeat (a) to (c).

[0063] 2. Calculate the accumulated area of flow concentration according to the flow direction.

[0064] The accumulated area of flow concentration is defined as the total area of all upstream flow concentration regions of a certain grid cell (or river network node), that is, the sum of the catchment areas covered by the water flow paths flowing to this grid.

[0065] The accumulated length of flow concentration can be calculated through the following steps:

[0066] (3) Calculate the number of inflowing grids for each grid. 1. Initialize the number of inflowing grids for each grid to zero.

[0067] 2. Traverse the flow direction grids and increment the number of inflowing grids of their downstream grids by one.

[0068] (4) Calculate the accumulated area of flow concentration. 1. Initialize the accumulated area of flow concentration for each grid to be equal to the grid area. 2. Traverse all grids to find the grids with zero number of inflowing grids, which are the river sources.

[0069] 3. Select one of the river source grids in turn and set it as the current grid, and search downstream according to the flow direction: (a) Find the downstream grid of the current grid according to the flow direction. If there is no downstream for the current grid or the downstream is outside the search area, the search terminates; (b) If there is a downstream grid, update the accumulated area of flow concentration of the downstream grid, and its value is equal to the sum of the accumulated area of flow concentration of the downstream grid and the current downstream grid; (c) If the number of inflowing grids of the downstream grid is greater than or equal to two, subtract one from the number of inflowing grids of the downstream grid and the search terminates; otherwise, set the downstream grid as the current grid and repeat (a) to (c).

[0070] 3. Extract the centerline of the river network according to the flow direction, the accumulated length of flow concentration, and the accumulated area of flow concentration.

[0071] In this embodiment, the centerline of the river network is represented by the grid where the centerline of the river network is located, and the vector of the centerline of the river network can be further obtained by vectorizing the grid where the centerline of the river network is located.

[0072] The river network centerline is extracted by gradually searching upstream from the basin outlet raster. The specific steps are as follows.

[0073] (1) Initialize an empty centerline temporary queue and an empty centerline completed queue. Both queues store the river reach centerlines, and each river reach centerline is represented by a raster sequence. The front-to-back order of the rasters represents the upstream and downstream relationships.

[0074] (2) Traverse the flow direction rasters. If a raster has no downstream raster or the downstream is outside the search area, then this raster is the basin outlet. Each basin outlet raster is used as a centerline raster sequence of length one and saved to the centerline temporary queue.

[0075] (3) Select the river reach centerline at the front end of the centerline temporary queue and delete it from the queue.

[0076] (4) For the selected river reach centerline, perform the following steps in a loop:

[0077] a) Calculate the target length of the river reach centerline. If the centerline has no upstream, that is, no centerline in the temporary queue and the completed queue flows to this centerline, then the target length is equal to the accumulated flow length at the location of the most downstream raster of the river reach centerline; otherwise, the target length is equal to the accumulated flow length at the location of the most downstream raster of the river reach centerline minus the sum of the accumulated flow lengths at the locations of the most downstream rasters of all its upstream centerlines. If the centerline length is greater than or equal to the target length, then this centerline is saved to the completed queue and the loop is terminated; otherwise, execute the following steps.

[0078] b) According to the flow direction, search for the set of upstream rasters of this centerline. The rasters in this set meet three conditions: First, they are adjacent to any raster in the centerline; second, they flow to the centerline; third, they are not included in any centerline in the temporary queue and the completed queue. If the set of upstream rasters is empty, then this river channel centerline is saved to the completed queue and the loop is terminated; otherwise, execute the following steps.

[0079] c) In the set of upstream rasters, compare the accumulated flow areas and find the raster with the largest accumulated flow area as the pending raster.

[0080] d) If the pending raster flows to the most upstream raster in the centerline, then add the pending raster to the front of the river channel centerline sequence, go back to step (a), and repeat the loop; otherwise, execute the following steps.

[0081] e) Find the downstream grid of the grid to be determined according to the flow direction. Taking the downstream grid of the grid to be determined as the boundary, divide the centerline into two sections: the upstream section and the downstream section. The upstream section does not include the downstream grid of the grid to be determined, and the downstream section includes the downstream grid of the grid. Save both the upstream and downstream sections of the centerline to the end of the centerline temporary queue.

[0082] f) Take the grid to be determined as a river reach centerline grid sequence with a length of one, save it to the end of the centerline temporary queue, and terminate the loop.

[0083] (5) Repeat steps (3) and (4) until the centerline temporary queue is empty. At this time, the centerline complete queue has saved the centerlines of all river reaches in the river network.

[0084] 4. Extract the catchment area according to the flow direction and the centerline.

[0085] In this embodiment, the catchment area is represented by the grids covered by the catchment area, and the catchment area boundary can be further obtained by vectorizing the catchment area grids.

[0086] In this embodiment, the catchment area corresponds to the river reach centerline one by one. For the centerline of each river reach in the river network, its corresponding catchment area can be extracted through the following steps:

[0087] (1) Establish an empty exclusion set. If there is an upstream centerline for the centerline, traverse the upstream centerline of the centerline and add the most downstream grid of the upstream centerline to the exclusion set;

[0088] (2) Initialize an empty catchment area temporary queue and an empty catchment area complete set. Both the temporary queue and the complete set save the grids covered by the catchment area.

[0089] (3) Push the grids in the centerline sequence into the end of the temporary queue in the order from downstream to upstream.

[0090] (4) If the temporary queue is not empty, take a grid from the front end of the temporary queue, save it to the complete set, and delete it from the temporary queue, and loop to execute the following steps.

[0091] (5) According to the flow direction, find all the grids upstream of the taken grid and included in the complete set or the exclusion set, and push these grids into the end of the temporary queue.

[0092] (6) Repeat steps (4) and (5) until the temporary queue is empty. At this time, the complete set saves the catchment area corresponding to the river reach centerline, and all the grids in the catchment area complete set flow to this centerline.

[0093] The above is only a preferred specific embodiment of the present application, but the protection scope of the present application is not limited thereto. Any changes or substitutions that can be easily conceived by those skilled in the art within the technical scope disclosed in the present application should be covered by the protection scope of the present application. Therefore, the protection scope of the present application should be subject to the protection scope of the claims.

Claims

1. An automatic extraction method for river network centerlines and catchment areas to maintain river network density, characterized in that Including: Calculating the cumulative length of flow concentration according to the river flow direction and river network density; Calculating the cumulative area of flow concentration according to the river flow direction; Extracting all river reaches and their centerlines according to the said flow direction, the said cumulative length of flow concentration and the said cumulative area of flow concentration; Extracting the catchment area according to the said flow direction and the said centerline; 2. The method according to claim 1, characterized in that The process of calculating the cumulative length of flow concentration according to the flow direction and river network density includes: Initializing the number of inflow grids of each grid to zero; Traversing the flow direction grid and adding one to the number of inflow grids of its downstream grid; Initializing the cumulative length of flow concentration of each grid to be equal to the grid area multiplied by the river network density; Traversing all grids to find the grids with zero inflow grid number as the river sources; Selecting one of the river source grids in turn, setting it as the current grid, and searching downstream according to the flow direction to calculate the cumulative length of flow concentration.

3. The method according to claim 2, characterized in that, The process of searching downstream according to the flow direction to calculate the cumulative length of flow concentration includes: Finding the downstream grid of the current grid according to the flow direction. If there is no downstream for the current grid, the search terminates; If there is a downstream grid, updating the cumulative length of flow concentration of the downstream grid, and its value is equal to the sum of the cumulative length of flow concentration of the downstream grid and the current downstream grid; If the number of inflow grids of the downstream grid is greater than or equal to two, subtracting one from the number of inflow grids of the downstream grid and terminating the search. Otherwise, setting the downstream grid as the current grid and re-finding the downstream grid of the current grid according to the flow direction for judgment.

4. The method according to claim 1, wherein The process of calculating the cumulative area of flow concentration according to the flow direction includes: Initializing the number of inflow grids of each grid to zero; Traversing the flow direction grid and adding one to the number of inflow grids of its downstream grid; Initializing the cumulative area of flow concentration of each grid to be equal to the grid area; Traversing all grids to find the grids with zero inflow grid number as the river sources; Selecting one of the river source grids in turn, setting it as the current grid, and searching downstream according to the flow direction to calculate the cumulative area of flow concentration.

5. The method according to claim 4, characterized in that, The process of searching downstream according to the flow direction to calculate the cumulative area of flow concentration includes: Finding the downstream grid of the current grid according to the flow direction. If there is no downstream for the current grid or the downstream is outside the search area, the search terminates; If there is a downstream grid, updating the cumulative area of flow concentration of the downstream grid, and its value is equal to the sum of the cumulative area of flow concentration of the downstream grid and the current downstream grid; If the number of inflow grids of the downstream grid is greater than or equal to two, subtracting one from the number of inflow grids of the downstream grid and terminating the search; otherwise, setting the downstream grid as the current grid and re-finding the downstream grid of the current grid according to the flow direction for judgment.

6. The method according to claim 1, wherein The process of extracting all river reaches and their centerlines according to the flow direction, the cumulative length of flow concentration and the cumulative area of flow concentration includes: Initializing an empty centerline temporary queue and an empty centerline completed queue. Both queues save the centerlines of river reaches, and each centerline of a river reach is represented by a grid sequence, and the front-back order of the grids represents the upstream-downstream relationship; Traversing the flow direction grid. If a grid has no downstream grid or the downstream is outside the search area, then this grid is the basin outlet. Taking each basin outlet grid as a centerline grid sequence with a length of one and saving it to the centerline temporary queue; Selecting the centerline of the river reach at the forefront of the centerline temporary queue and deleting it from the queue; Calculate the center line of the selected river section to obtain the data of all river sections and their center lines.

7. The method according to claim 6, wherein The process of calculating the center line of the selected river section to obtain the data of all river sections and their center lines includes: Calculate the target length of the center line of this river section. If the center line of the river section has no upstream, that is, there is no center line flowing to this center line in the temporary queue and the completed queue, the target length is equal to the cumulative length of flow concentration at the grid where the most downstream grid of the center line of the river section is located. Otherwise, the target length is equal to the cumulative length of flow concentration at the grid where the most downstream grid of the center line of the river section is located minus the sum of the cumulative lengths of flow concentration at the grids where the most downstream grids of all its upstream center lines are located. If the length of the center line is greater than or equal to the target length, save this center line to the completed queue and terminate the loop; If the length of the center line is less than the target length, search for the set of upstream grids of this center line according to the flow direction; Compare the cumulative areas of flow concentration based on the set of upstream grids to find the grid with the largest cumulative area of flow concentration as the grid to be determined; If the grid to be determined flows to the most upstream grid in the center line, then add the grid to be determined to the front of the river center line sequence, recalculate the target length of the center line of this river section and make a judgment; If the grid to be determined does not flow to the most upstream grid in the center line, find the downstream grid of the grid to be determined according to the flow direction. Taking the downstream grid of the grid to be determined as the boundary, divide the center line into two sections: the upstream section does not include the downstream grid of the grid to be determined, and the downstream section includes the downstream grid of the grid to be determined. Save the center lines of the upstream and downstream sections to the end of the center line temporary queue; Save the grid to be determined as a river center line grid sequence with a length of one to the end of the center line temporary queue; Take a river section from the front end of the temporary queue again, delete it from the queue, calculate the target length of the center line of this river section and make a judgment; repeat the judgment process until the center line temporary queue is empty. At this time, the center line completed queue saves the center lines of all river sections in the river network.

8. The method according to claim 1, wherein The process of extracting the catchment area according to the flow direction and the center line includes: Establish an empty exclusion set. If the center line has upstream center lines, traverse the upstream center lines of the center line and add the most downstream grid of the upstream center line to the exclusion set; Initialize an empty catchment area temporary queue and an empty catchment area completed set. Both the temporary queue and the completed set save the grids covered by the catchment area; Push the grids in the center line sequence into the end of the temporary queue in the order from downstream to upstream; If the temporary queue is not empty, take a grid from the front end of the temporary queue, store it in the completed set, and delete it from the temporary queue. According to the flow direction, find all the grids upstream of the taken grid and not included in the completed set or the exclusion set, and push these grids into the end of the temporary queue. Repeat until the temporary queue is empty to complete the extraction of the catchment area.

9. The method according to claim 1, wherein The river network center line is represented by the grid where the river network center line is located, and the river network center line vector is obtained by vectorizing the grid where the river network center line is located.

Images