Digital twin basin river network high-precision extraction method

Abstract

The invention relates to the technical field of hydrology, and discloses a digital twin basin river network high-precision extraction method, which comprises the steps of obtaining DEM data, performing splicing and cutting processing, and judging whether to perform depression filling processing so as to generate DEM data after basin depression filling; calculating the water flow direction of each unit grid through unit grid division, generating flow direction data and calculating accumulated flow data of each unit grid; setting a water collection threshold value, judging whether accumulated flow data of each unit grid exceeds the water collection threshold value or not, generating a first water body unit, performing river linking, river network grading and vectorization processing, generating vectorized river network data for calculating an improved river network density index, and calculating an optimal water collection threshold value; generating a digital river network corresponding to the optimal water collection threshold value and calculating the river network density so as to establish a joint evaluation index based on the registration difference and the river network density for evaluating the precision of the digital river network; according to the method, the digital river network extraction precision and the extraction result evaluation accuracy are improved.

Description

Technical Field

[0001] The present invention relates to the technical field of hydrology, and in particular to a method for high-precision extraction of river networks in a digital twin watershed. Background Art

[0002] River networks play a vital role in hydrology and water resources management, and provide key data for multiple interdisciplinary fields, such as natural disaster monitoring and early warning (such as flood, drought and landslide monitoring), urban planning and management, etc. At present, the mainstream methods for extracting river networks are: river network extraction methods based on GIS and DEM and river network extraction methods based on remote sensing images. The river network extraction method based on remote sensing images is easily affected by factors such as mountain shadows and cloud cover, resulting in inaccurate extracted river information and difficulty in extracting small rivers. Therefore, most of the current river network extraction studies rely on GIS and DEM data, but the accuracy of river networks extracted based on GIS and DEM data is affected by many factors, among which the setting of water collection threshold is one of the key parameters affecting the accuracy results. The water collection threshold determines the accumulated amount of water flow, thereby affecting the identification and classification of rivers. If the water collection threshold is too large, pseudo rivers are likely to be generated, and if the water collection threshold is too small, small rivers may be lost. Therefore, determining the optimal water collection threshold is the key to extracting river networks.

[0003] Among them, the river network density method is the mainstream method for determining the optimal water collection threshold. Although it is widely used in various basin experiments, there is still a certain gap between the river network extracted by the river network density method and the actual river network. This gap is caused by the inaccurate determination of the optimal water collection threshold. In addition, the current important reference indicator for evaluating the accuracy of river network extraction is the fit difference. Although many researchers use the size of the fit difference to judge the quality of river network extraction results after extracting the digital river network, the fit difference is due to its calculation principle, that is, the area of ​​the fragmented area formed by the extracted river network and the actual river network water system lines occupies the entire basin. The ratio of the area cannot accurately evaluate the accuracy of river network extraction. The reason is that when the fitting error is the only indicator to evaluate the accuracy of river network, as the river length decreases, the river network density decreases, the intersection points of the extracted river network and the real river network decrease, the area of ​​the obtained fragmented polygons decreases, and the fitting error decreases. It is believed that the accuracy is improved, and it can be concluded that the smaller the river network density, the higher the accuracy of river network extraction. However, when the river network density is too small, the river network becomes sparse, and many small rivers are not extracted, and the river network accuracy is too low. This is inconsistent with the conclusion that the smaller the river network density, the higher the accuracy of river network extraction. Summary of the invention

[0004] In view of the above-mentioned deficiencies in the prior art, the present invention provides a high-precision river network extraction method for a digital twin watershed, which is used to solve the problems of insufficient accuracy of the water collection threshold and inaccurate extracted river network when extracting the water collection threshold using the existing river network density method. At the same time, an evaluation index of river network density is considered on the basis of the existing fitting difference index to solve the problem of inaccurate evaluation of river network extraction accuracy.

[0005] In order to achieve the above-mentioned object of the invention, the technical solution adopted by the present invention is:

[0006] A high-precision extraction method for a digital twin river network includes the following steps:

[0007] S1. Obtain the DEM data of the study area and perform splicing and cropping to generate the DEM data of the watershed. At the same time, determine whether there is a depression in the DEM data of the watershed. If so, fill the depression to obtain the DEM data of the watershed after filling. Otherwise, do not process it and directly use the DEM data of the watershed as the basic terrain data for river network extraction.

[0008] S2, divide the DEM data after the basin is filled or the basic terrain data extracted from the river network into unit grids, calculate the water flow direction of each unit grid, generate flow direction data, and calculate the cumulative flow data of each unit grid based on the flow direction data;

[0009] S3, by setting a water collection threshold, determining whether the accumulated flow data of each unit grid exceeds the water collection threshold, if so, identifying each unit grid exceeding the water collection threshold as a first water body unit, otherwise, not identifying it as a first water body unit;

[0010] S4, performing river linking, river network classification, and vectorization processing on all first water body units to generate vectorized river network data;

[0011] S5. Based on the vectorized river network data, the length of the final river and the length of the total river are obtained, and the improved river network density index is calculated;

[0012] S6. Establish a fitting function of the improved river network density index and the water collection threshold and perform a second-order derivative, take the point where the fitting function approaches zero as the optimal water collection threshold point, take the water collection threshold corresponding to the optimal water collection threshold point as the optimal water collection threshold, and again judge whether the cumulative flow data of each unit grid exceeds the optimal water collection threshold. If so, each unit grid exceeding the optimal water collection threshold is identified as a second water body unit, otherwise, it is not identified as a second water body unit;

[0013] S7. Link rivers, classify river networks, and vectorize all second water body units to generate a digital river network corresponding to the optimal water collection threshold and obtain the river network density of the digital river network. Establish a joint evaluation index based on the fit difference and river network density to evaluate the accuracy of the digital river network.

[0014] The present invention has the following beneficial effects:

[0015] A high-precision river network extraction method for a digital twin watershed proposed in the present invention takes the final river density as the main reference object. Compared with the traditional river network density method which takes the entire river network density as the reference object, it can determine a more accurate optimal water collection threshold and thus extract more accurate river network results. Compared with traditional river network extraction result accuracy evaluation indicators, the joint evaluation indicator proposed in the present invention simultaneously considers the influence of the fitting error and river network density on the river network extraction accuracy evaluation, thereby making the evaluation of the river network extraction effect more accurate. BRIEF DESCRIPTION OF THE DRAWINGS

[0016] Figure 1 A schematic diagram of the process of a high-precision extraction method of river network in a digital twin watershed proposed in the present invention;

[0017] Figure 2 Schematic diagram of the process of filling depression in the embodiment;

[0018] Figure 3 Schematic diagram of water flow direction process using D8 algorithm in the embodiment;

[0019] Figure 4 This is a schematic diagram of cumulative flow calculation in the embodiment;

[0020] Figure 5 This is a schematic diagram of river connection in the embodiment;

[0021] Figure 6 Schematic diagram of river network classification using the Strahler classification method and the Shreve classification method in the embodiment;

[0022] Figure 7 This is a schematic diagram of the 30m resolution river network density curve fitting in the embodiment;

[0023] Figure 8 This is a schematic diagram of the 90m resolution river network density curve fitting in the embodiment;

[0024] Fig. 9 This is a schematic diagram of the 450m resolution river network density curve fitting in the embodiment;

[0025] Fig.10 This is a schematic diagram of the 30m resolution river network extraction effect in the embodiment;

[0026] Fig.11This is a schematic diagram of the river network extraction effect at a resolution of 90m in the embodiment;

[0027] Fig.12 This is a schematic diagram of the river network extraction effect at a resolution of 450m in the embodiment;

[0028] Fig.13 It is a schematic diagram of local comparison of the 30m resolution river network extraction results in the embodiment. DETAILED DESCRIPTION

[0029] The specific implementation modes of the present invention are described below so that those skilled in the art can understand the present invention. However, it should be clear that the present invention is not limited to the scope of the specific implementation modes. For those of ordinary skill in the art, as long as various changes are within the spirit and scope of the present invention as defined and determined by the attached claims, these changes are obvious, and all inventions and creations utilizing the concept of the present invention are protected.

[0030] like Figure 1 As shown, a high-precision extraction method for river networks in a digital twin watershed includes the following steps S1-S7:

[0031] S1. Obtain the DEM data of the study area and perform splicing and cropping to generate the DEM data of the watershed. At the same time, determine whether there are depressions in the DEM data of the watershed. If so, fill the depressions to obtain the DEM data of the watershed after filling. Otherwise, do not process it and directly use the DEM data of the watershed as the basic terrain data for river network extraction.

[0032] In this embodiment, the ArcGIS hydrological analysis module is used to perform depression filling, flow direction calculation, runoff accumulation calculation, river network extraction and other operations on the DEM data, so as to obtain hydrological information such as the length of the river in the basin, the river network structure, and the accumulated flow, so as to realize the digital river network extraction based on the DEM data; however, due to the uncertainty factors such as the accuracy of the DEM data, "false depressions" will appear; and the existence of "false depressions" may cause interruptions or unreasonable redirections in the subsequent river network extraction, forming "pseudo-river channels", and further affecting the extraction of the digital river network in the basin; therefore, before performing subsequent processing, the DEM data needs to be corrected first, that is, the depression filling process, the principle of which is to modify the elevation of the depression point to the elevation of the pouring point, as shown in the following figure. Figure 2 As shown in the figure, the condition for the existence of depression points and pouring points is that in a unit confluence area, the elevation of the cell in the center is less than the elevation of all the surrounding boundary cells. In this case, it is considered that there are depression points and pouring points in this unit confluence area, and the central cell is a depression point, and the cell with the lowest elevation around it is a pouring point. In addition, a unit confluence area is as follows: Figure 3As shown in the figure, it is a nine-square grid composed of a cell and its surrounding cells; the elevation difference between the depression point and the pouring point is set to z, so the filling result of the unit confluence area can be determined by setting the size of z, that is, the threshold of the elevation difference. If a unit confluence area has a depression and the elevation difference between the depression point and the pouring point is less than the set elevation difference threshold, the elevation of the depression point is changed to the pouring point elevation to achieve depression filling. If it is greater than the elevation difference threshold, it will not be processed.

[0033] Specifically, the specific process of the depression filling process in step S1 is as follows:

[0034] The DEM data of the watershed is divided into unit grids to generate several unit confluence areas; each unit confluence area is an area composed of a unit grid and its adjacent unit grids.

[0035] Determine whether there are depressions and pouring points in each unit confluence area, and the elevation difference between the depression and the pouring point is less than the set elevation difference threshold. If so, adjust the elevation of the depression to the elevation of the pouring point to obtain the DEM data after the basin is filled. Otherwise, no processing is performed.

[0036] Specifically, the condition for the existence of depressions and pouring points is that in a unit confluence area, if the elevation of the central cell is less than the elevations of all surrounding boundary cells, then it is considered that there are depressions and pouring points in the unit confluence area, and the central cell of the unit confluence area is the depression point, and the cell with the lowest elevation of all surrounding boundary cells is the pouring point.

[0037] S2. Divide the DEM data after basin filling or the basic terrain data extracted from the river network into unit grids, calculate the water flow direction of each unit grid, generate flow direction data, and calculate the cumulative flow data of each unit grid based on the flow direction data.

[0038] In this embodiment, the D8 algorithm is used to calculate the river flow direction in ArcGIS, that is, the water flow of each cell can flow in one of the eight directions around it, and the adjacent cell with the largest slope with the central cell is selected as the next flow direction of the water flow, such as Figure 3 As shown, the specific operation process is as follows: step S21; in addition, the flow direction file (flow direction data) can be input into the flow hydrology analysis module of ArcGIS to calculate the cumulative amount of runoff. The basic principle is to use grids as units, each point represents a unit of water volume, and according to the natural law that natural water flows from high to low, the cumulative amount of water flowing through each grid is calculated according to the flow direction data, so as to obtain the cumulative amount of runoff in each unit runoff area; and when the runoff volume accumulates to a certain extent, water flow is generated on the surface, and all grids with runoff volume exceeding the critical value are potential water body cells, and these water body cells are connected to form a river network, as shown in the figure. Figure 4As shown, if the watershed area is a 16×16 grid area, after each grid obtains the flow direction through the D8 algorithm, assuming that each grid has a unit of water volume, then the cumulative flow of each grid is the cumulative water volume flowing through this grid in the entire watershed, and the operation process is as shown in the following steps S22-S23.

[0039] Specifically, step S2 includes S21-S23:

[0040] S21. Divide the DEM data after basin filling or the basic terrain data extracted from the river network into unit grids, and use the D8 algorithm to select the adjacent unit grid with the largest slope to the central unit grid as the next flow direction of the water flow, and calculate the flow direction of each unit grid, that is:

[0041]

[0042] Among them, D ij represents the flow direction of the cell grid (i, j), argmin represents the independent variable used to find the minimum value, N(i, j) represents the eight adjacent cell grids of the cell grid (i, j), and Z ij Represents the elevation of the cell grid (i, j), Z kl Represents the elevation of the cell grid (k,l).

[0043] S22. Generate flow direction data based on the water flow direction of each unit grid.

[0044] S23. Calculate the cumulative flow data of each unit grid according to the flow direction data, that is:

[0045]

[0046] Among them, A ij represents the cumulative flow data of the unit grid (i, j), U(i, j) represents the set of all units flowing to the unit grid (i, j), and A kl Represents the cumulative amount of water flowing from cell grid (k,l) to cell grid (i,j).

[0047] S3. By setting a water collection threshold, determine whether the accumulated flow data of each unit grid exceeds the water collection threshold. If so, identify each unit grid exceeding the water collection threshold as a first water body unit; otherwise, do not identify it as a first water body unit.

[0048] In this embodiment, a designated water collection threshold is set to obtain water body unit data results. If the cumulative flow exceeds the set water collection threshold, it is determined that a river exists in the area.

[0049] S4. Perform river linking, river network classification, and vectorization processing on all first water body units to generate vectorized river network data.

[0050] In this embodiment, the single rivers identified in the above steps (wherein the rivers can be obtained by connecting the water body units) are linked through rivers to obtain a river network. Figure 5 As shown, Figure 5 The red dots in the middle are the intersections of each river, and the black lines are the river links. At the same time, the river network is classified and each river is assigned a different level. The classification methods commonly used are Strahler classification and Shreve classification. Among them, Strahler classification defines the source river as level 1. The river network classification will only increase when rivers of the same level meet. Therefore, the intersection of level 1 connection lines and level 2 connection lines will retain level 2 connection lines. The Shreve classification also defines the source as a first-order river, but the river level is the sum of the levels of all previous rivers. For example, the intersection of two first-order connecting lines will create a second-order connecting line, the intersection of a first-order connecting line and a second-order connecting line will create a third-order connecting line, and the intersection of a second-order connecting line and a third-order connecting line will create a fifth-order connecting line. Strahler classification and Shreve classification are as follows: Figure 6 As shown; in addition, the present invention uses the Strahler classification method as the river network classification method for subsequent operations. Finally, the obtained river network classification raster data is converted into a river network in vector format through vectorization processing, thereby completing the river network data conversion from raster to vector; and then the river length of the vector data can be calculated to obtain the length of the first-level river, that is, the final-level river, and the sum of the lengths of all rivers (the total length of the river).

[0051] Specifically, step S4 includes S41-S43:

[0052] S41. Link all first water body units to rivers to generate a river network.

[0053] S42. The river network is classified using the Strahler classification method. By assigning different levels to each river in the river network, the classified river network raster data is obtained, specifically:

[0054] The source rivers in the river network are defined as level 1. Only when rivers of the same level meet, the level of the rivers in the river network is increased. That is, when a level 1 river meets a level 1 river, a level 2 river is generated. When a level 1 river meets a level 2 river, the level 2 river attributes are retained. When a level 2 river meets a level 2 river, a level 3 river is generated, and so on, until the main stream of the basin is generated and the river network classification is completed, that is, the classified river network raster data is obtained.

[0055] S43, vectorizing the classified river network raster data to generate vectorized river network data.

[0056] S5. Based on the vectorized river network data, the length of the final river and the length of the total river are obtained, and the improved river network density index is calculated.

[0057] In this embodiment, the length of the first-level river, i.e., the final river, and the sum of the lengths of all rivers are calculated through vector data (vectorized river network data), and an improved river network density index is calculated. Among them, the improved river network density index is proposed based on the principle of the river network density method, specifically: as the water collection threshold increases, the total length of the river decreases, and the change in the total length of the river is mainly caused by the decrease in the length of the final river, and the length of the main trunk of the river, such as the main stream, does not change much; based on this, the total length of the river L is divided into the final river length and the length of other rivers, and the latter is given a coefficient weight to adjust the relationship between the two, reduce the proportion of other river lengths, so that the curve change can better represent the changing trend of the final river length, thereby finding a more accurate turning point as the optimal water collection threshold. Among them, the calculation formula of the improved river network density index used in the improved river network density method is: I = (L a +B*(LL a )) / S.

[0058] Specifically, step S5 includes S51-S52:

[0059] S51. Based on the vectorized river network data, the length of the last-level river and the length of the total river are obtained; wherein the length of the last-level river is the length of the source river of the river network.

[0060] In this embodiment, the total river length also refers to the length of rivers of all levels and the sum of the lengths of all rivers. In the above steps, a vector file of the river network can be obtained. The attribute of the vector file is a line file. The total river length can be obtained by calculating the lengths of all lines in the line file. The calculation principle is as follows: the size of the cell is determined by the initial DEM data resolution. A 30m resolution is a 30m×30m cell. The position of each cell identified as a water body is determined, and the center points of the water body cells are connected to form a river network. The sum of the lengths of all connecting lines is the total river length. Generally, the calculation can be performed using geographic information software - ArcGIS software. At the same time, the length of each individual river can be calculated, or the length of each level of the river can be calculated according to the river level.

[0061] S52. Calculate the improved river network density index based on the length of the final river and the length of the total river, namely:

[0062] I=(L a +B*(LL a )) / S

[0063] Among them, I represents the improved river network density index, L a It represents the length of the last-stage river, L represents the total length of the river, B represents the first empirical coefficient, and S represents the total area of ​​the basin.

[0064] Among them, the value range of the first empirical coefficient is 0 to 1, and it is found through experimental verification that the improvement effect is best when the value of B is 0.

[0065] S6. Establish a fitting function between the improved river network density index and the water collection threshold and perform a second-order derivative. Take the point where the fitting function approaches zero as the optimal water collection threshold point. Take the water collection threshold corresponding to the optimal water collection threshold point as the optimal water collection threshold. Determine again whether the cumulative flow data of each unit grid exceeds the optimal water collection threshold. If so, each unit grid that exceeds the optimal water collection threshold is identified as a second water body unit. Otherwise, it is not identified as a second water body unit.

[0066] In this embodiment, the calculated improved river network density index and the water collection threshold are fitted with a function. Generally, the two satisfy a power function relationship. The second-order derivative of the fitting function is then solved, and the water collection threshold corresponding to the point where the second-order derivative approaches 0 is the optimal water collection threshold.

[0067] Specifically, step S6 specifically includes S61-S62:

[0068] S61. Establish a fitting function between the improved river network density index and the water collection threshold, and perform a second-order derivative on the fitting function to obtain the point where the second-order derivative approaches zero as the optimal water collection threshold point, and use the water collection threshold corresponding to the optimal water collection threshold point as the optimal water collection threshold.

[0069] In this embodiment, the calculation formula for taking the second-order derivative of the fitting function is: Where f′(x) is the first-order derivative, f″(x) is the second-order derivative, f(x) is the function value at x, f(x+h) is the function value at (x+h), and h is a value close to 0.

[0070] S62. Based on the optimal water collection threshold, determine again whether the cumulative flow data of each unit grid exceeds the optimal water collection threshold. If so, identify each unit grid that exceeds the optimal water collection threshold as a second water body unit. Otherwise, do not identify it as a second water body unit.

[0071] S7. Link rivers, classify river networks, and vectorize all second water body units to generate a digital river network corresponding to the optimal water collection threshold and obtain the river network density of the digital river network. Establish a joint evaluation index based on the fit difference and river network density to evaluate the accuracy of the digital river network.

[0072] In this embodiment, a joint evaluation index based on the fitting difference and the river network density is established to evaluate the accuracy of the river network density extracted by the optimal water collection threshold and the density of the actual river network. Specifically, the river network density difference is used as an auxiliary verification indicator for accuracy evaluation, and the smaller the difference, the better the extraction effect. In order to present this indicator digitally, a certain weight is given to the river network density difference to participate in accuracy verification. The specific calculation formula is: i=c*|Dd| / d+T.

[0073] Specifically, step S7 specifically includes S71-S73:

[0074] S71. Perform river linking, river network classification, and vectorization processing on all second water body units to generate a digital river network corresponding to an optimal water collection threshold, and obtain the river network density of the digital river network.

[0075] S72. Based on the river network density of the digital river network, a joint evaluation index based on the overlay difference and the river network density is established, namely:

[0076] i=c*|Dd| / d+T

[0077] Among them, i represents the joint evaluation index based on the fit difference and river network density, c represents the second empirical coefficient, D represents the river network density of the digital river network extracted using the optimal water collection threshold, d represents the actual river network density, and T represents the fit difference; among them, since the fit sum difference value range is 0 to 0.05, and the river network density difference range is 0 to 0.5, which differs by one order of magnitude, c can be set to 0.1.

[0078] In this embodiment, the river network fitting difference calculation formula is: Among them, T is the river network fitting difference, A i The area of ​​the fine-grained area is km 2 ; S is the total area of ​​the basin km 2 ; In addition, the total basin area in this embodiment refers to the area of ​​the basin where the river network is to be extracted, that is, the DEM data needs to be spliced ​​and cropped in the initial step, and the cropping process needs to determine the boundary of the study basin, such as the Jinsha River Basin. Once the boundary is determined, the calculated area of ​​the basin is the total basin area.

[0079] S73. Use a joint evaluation index based on the fit error and river network density to evaluate the accuracy of the digital river network.

[0080] In this embodiment, in order to verify the effectiveness of the high-precision extraction method of the river network of the digital twin river basin proposed in the present invention, the Jinsha River Basin was selected as the research area for the experiment, as follows:

[0081] Among them, the DEM products are: (1) ASTER GDEMV3 (30m resolution); (2) SRTM DEM (90m resolution); (3) GEBCODEM (450m resolution); the verification data is: the real river network and the map comes from the National Basic Geographic Information Center of the Ministry of Natural Resources; the research plan is set as follows: in order to discuss the advantages of the improved river network density method, DEM data with different resolutions are selected for river network extraction, specifically: first, the three types of DEM data and the traditional river network density method are used to calculate the optimal threshold and extract the digital river network; secondly, the improved river network density method is used to extract the optimal threshold and digital river network of the three types of DEM data; finally, the river network data published by the Ministry of Natural Resources is used to evaluate and verify the accuracy of the river network extracted by the two methods based on the joint evaluation index. Based on the above data foundation and experimental process, the experimental results are as follows:

[0082] (1) River network extraction results

[0083] The three types of DEM data of three basins were processed and analyzed, and appropriate water collection thresholds were set to extract the river network and obtain the river network density. At the same time, a curve diagram of the relationship between river network density and water collection threshold was drawn. The results are shown in the figure. Figure 7-Figure 9 As shown in Figure 2, it can be found that there is a power function relationship between the water collection threshold and the river network density, and the fitting correlation coefficient R 2 The values ​​of the two fitting curves are all greater than 0.9, indicating a good fitting effect. At the same time, the second-order derivative of the fitting curves obtained from different DEM resolution data is used to obtain the inflection point, and the optimal water collection thresholds for the three DEM data are 2.4 km 2 (30m resolution), 2.5km 2 (90m resolution), 2.6km 2 (450m resolution); In addition, the DEM data is processed in the same way. When calculating the river network vector file, the river network density index is calculated according to the improved river network density method, and the optimal water collection threshold is calculated with the water collection threshold fitting relationship curve. The river network extraction effect is as follows: Figure 10-12 As shown, it can be found that the result is 1.8km 2 (30m resolution), 1.9km 2 (90m resolution), 2.1km 2 (450m resolution), which is completely different from the optimal water collection threshold obtained by the river network density method and is generally smaller.

[0084] (2) Verification of river network extraction results

[0085] The accuracy of the extracted river network results is verified, and the results are shown in Table 1:

[0086] Table 1 Accuracy evaluation table of river network extraction results

[0087]

[0088] As shown in Table 1, the improvement of river network extraction evaluation index of GDEMV3 data in Jinsha River Basin is 0.56% (increasing river length by 38,000 km), the improvement of river network evaluation index of SRTM data is 0.74% (increasing river length by 45,000 km), and the improvement of river network evaluation index of GEBCODEM data is 0.79% (increasing river length by 46,000 km). Different DEM data show different degrees of accuracy improvement. Fig.13 The comparison of river network extraction results of 30m resolution data is specifically shown. Since the results of 90m and 450m resolution data are similar, only the 30m resolution data results are shown here. Fig.13 The river network of the Jinsha River Basin in China is divided into the traditional river network density method, the improved river network density method and the actual river network distribution from top to bottom. Among them, A and B show the extraction effect of the two terminal rivers in the Jinsha River Basin, while A1, A2 and B1, B2 are the enlarged display of the key areas. Figure 4 It can be observed that the main river channels extracted by the two methods are basically consistent with the actual river network, and the main difference is reflected in the extraction of river channels at the tributary level. Specifically, in area A, although both methods can identify the tributary, the length of the river channel extracted by the improved method proposed in the present invention is closer to the actual situation. In area B, the traditional method failed to identify the tributary, while the improved method proposed in the present invention successfully achieved the extraction. It can be confirmed from the local enlarged images of A1, A2, B1, and B2 that there are indeed real small river channels in these areas, rather than false "pseudo-river channels". The results show that the improved river network density method has stronger adaptability and can effectively process data with different terrain features and DEM resolutions, and its extraction accuracy is significantly improved compared with the traditional method.

[0089] To sum up, the high-precision river network extraction method for a digital twin watershed proposed in the present invention, first of all, compared with the traditional river network density method which takes the entire river network density as a reference object, the improved river network density method proposed in the present invention takes the final river density as the main reference object, and can calculate a more accurate optimal water collection threshold and obtain a more accurate river network extraction result; at the same time, compared with the traditional accuracy evaluation index, the joint evaluation index proposed in the present invention takes into account the influence of the fitting difference and river network density on the river network extraction accuracy, and can more accurately evaluate the river network extraction effect; therefore, the method proposed in the present invention improves the accuracy of the river network extraction results.

[0090] The present invention uses specific embodiments to illustrate the principles and implementation methods of the present invention. The description of the above embodiments is only used to help understand the method of the present invention and its core idea. At the same time, for those skilled in the art, according to the idea of ​​the present invention, there will be changes in the specific implementation methods and application scope. In summary, the content of this specification should not be understood as a limitation on the present invention.

[0091] Those skilled in the art will appreciate that the embodiments described herein are intended to help readers understand the principles of the present invention, and should be understood that the protection scope of the present invention is not limited to such specific statements and embodiments. Those skilled in the art can make various other specific variations and combinations that do not deviate from the essence of the present invention based on the technical revelations disclosed by the present invention, and these variations and combinations are still within the protection scope of the present invention.

Claims

1. A high-precision extraction method for river networks in digital twin watersheds, characterized in that: The following steps are involved: S1. Obtain the DEM data of the study area and perform splicing and cropping to generate the DEM data of the watershed. At the same time, determine whether there is a depression in the DEM data of the watershed. If so, fill the depression to obtain the DEM data of the watershed after filling. Otherwise, do not process it and directly use the DEM data of the watershed as the basic terrain data for river network extraction. S2, divide the DEM data after the basin is filled or the basic terrain data extracted from the river network into unit grids, calculate the water flow direction of each unit grid, generate flow direction data, and calculate the cumulative flow data of each unit grid based on the flow direction data; S3, by setting a water collection threshold, determining whether the accumulated flow data of each unit grid exceeds the water collection threshold, if so, identifying each unit grid exceeding the water collection threshold as a first water body unit, otherwise, not identifying it as a first water body unit; S4, performing river linking, river network classification, and vectorization processing on all first water body units to generate vectorized river network data; S5. Based on the vectorized river network data, the length of the final river and the length of the total river are obtained, and the improved river network density index is calculated; S6. Establish a fitting function of the improved river network density index and the water collection threshold and perform a second-order derivative, take the point where the fitting function approaches zero as the optimal water collection threshold point, take the water collection threshold corresponding to the optimal water collection threshold point as the optimal water collection threshold, and again judge whether the cumulative flow data of each unit grid exceeds the optimal water collection threshold. If so, each unit grid exceeding the optimal water collection threshold is identified as a second water body unit, otherwise, it is not identified as a second water body unit; S7. Link rivers, classify river networks, and vectorize all second water body units to generate a digital river network corresponding to the optimal water collection threshold and obtain the river network density of the digital river network. Establish a joint evaluation index based on the fit difference and the river network density to evaluate the accuracy of the digital river network.

2. The high-precision extraction method of river network in digital twin watershed according to claim 1 is characterized in that: The specific process of the depression filling process in step S1 is as follows: The DEM data of the watershed is divided into unit grids to generate several unit confluence areas; each unit confluence area is an area composed of a unit grid and its adjacent unit grids; Determine whether there are depressions and pouring points in each unit confluence area, and the elevation difference between the depression and the pouring point is less than the set elevation difference threshold. If so, adjust the elevation of the depression to the elevation of the pouring point to obtain the DEM data after the basin is filled. Otherwise, no processing is performed.

3. The high-precision extraction method of river network in digital twin watershed according to claim 2 is characterized in that: The condition for the existence of depressions and pouring points is that in a unit confluence area, if the elevation of the central cell is lower than the elevations of all the surrounding boundary cells, then it is considered that there are depressions and pouring points in the unit confluence area, and the central cell of the unit confluence area is the depression point, and the cell with the lowest elevation of all the surrounding boundary cells is the pouring point.

4. The high-precision extraction method of river network in digital twin watershed according to claim 3 is characterized in that: Step S2 specifically includes: S21. Divide the DEM data after basin filling or the basic terrain data extracted from the river network into unit grids, and use the D8 algorithm to select the adjacent unit grid with the largest slope to the central unit grid as the next flow direction of the water flow, and calculate the flow direction of each unit grid, that is: Among them, D ij represents the flow direction of the cell grid (i, j), argmin represents the independent variable used to find the minimum value, N(i, j) represents the eight adjacent cell grids of the cell grid (i, j), and Z ij Represents the elevation of the cell grid (i, j), Z kl Represents the elevation of the cell grid (k,l); S22, generating flow direction data based on the water flow direction of each unit grid; S23. Calculate the cumulative flow data of each unit grid according to the flow direction data, that is: Among them, A ij represents the cumulative flow data of the unit grid (i, j), U(i, j) represents the set of all units flowing to the unit grid (i, j), and A kl Represents the cumulative amount of water flowing from cell grid (k,l) to cell grid (i,j).

5. The high-precision extraction method of river network in digital twin watershed according to claim 4 is characterized in that: Step S4 specifically includes: S41, linking all first water body units to rivers to generate a river network; S42. The river network is classified using the Strahler classification method. By assigning different levels to each river in the river network, the classified river network raster data is obtained, specifically: The source river in the river network is defined as level 1. Only when rivers of the same level meet, the level of the river in the river network is increased. That is, when a level 1 river meets a level 1 river, a level 2 river is generated. When a level 1 river meets a level 2 river, the level 2 river attribute is retained. When a level 2 river meets a level 2 river, a level 3 river is generated. This process is repeated until the main stream of the basin is generated and the river network classification is completed, that is, the classified river network raster data is obtained. S43, vectorizing the classified river network raster data to generate vectorized river network data.

6. The high-precision extraction method of river network in digital twin watershed according to claim 5 is characterized in that: Step S5 specifically includes: S51. Based on the vectorized river network data, the length of the last-level river and the length of the total river are obtained; wherein the length of the last-level river is the length of the source river of the river network; S52. Calculate an improved river network density index based on the length of the final river and the length of the total river.

7. The high-precision extraction method of river network in digital twin watershed according to claim 6 is characterized in that: The formula for calculating the improved river network density index in step S52 is: I=(L a +B*(L-L a )) / S Among them, I represents the improved river network density index, L a It represents the length of the last-stage river, L represents the total length of the river, B represents the first empirical coefficient, and S represents the total area of ​​the basin.

8. The high-precision extraction method of river network in digital twin watershed according to claim 7 is characterized in that: Step S6 specifically includes: S61, establishing a fitting function of the improved river network density index and the water collection threshold, and taking the second-order derivative of the fitting function, obtaining the point where the second-order derivative approaches zero as the optimal water collection threshold point, and taking the water collection threshold corresponding to the optimal water collection threshold point as the optimal water collection threshold; S62. Based on the optimal water collection threshold, determine again whether the cumulative flow data of each unit grid exceeds the optimal water collection threshold. If so, identify each unit grid that exceeds the optimal water collection threshold as a second water body unit. Otherwise, do not identify it as a second water body unit.

9. The high-precision extraction method of river network in digital twin watershed according to claim 8 is characterized in that: Step S7 specifically includes: S71, linking rivers, classifying river networks, and vectorizing all second water body units to generate a digital river network corresponding to an optimal water collection threshold, and obtaining a river network density of the digital river network; S72. Based on the river network density of the digital river network, a joint evaluation index based on the overlay difference and the river network density is established; S73. Use a joint evaluation index based on the fit error and river network density to evaluate the accuracy of the digital river network.

10. The high-precision extraction method of river network in digital twin watershed according to claim 9 is characterized in that: The joint evaluation index based on the fit difference and the river network density in step S72 is: i=c*|Dd| / d+T Among them, i represents the joint evaluation index based on the fitting difference and river network density, c represents the second empirical coefficient, D represents the river network density of the digital river network extracted using the optimal water collection threshold, d represents the actual river network density, and T represents the fitting difference.

Images