A method and device for constructing a water system spatial map and a medium

Abstract

The application discloses a method and device for constructing a water system space atlas and a medium. The method comprises the following steps: preprocessing high-precision DEM data of the water system space atlas by using a depression filling analysis method to obtain final corrected DEM; calculating the total flow distribution proportion of the water system based on the final corrected DEM; determining a catchment area threshold based on the total flow distribution proportion; performing water system grading based on the catchment area threshold to obtain a water system grading result; and performing basin extraction and grading based on the final corrected DEM and the total flow distribution proportion to generate basin vector boundary data and obtain a complete water system space atlas.

Description

Technical Field

[0001] This invention relates to the field of water conservancy information technology, and more specifically, to a method, apparatus and medium for constructing a water system spatial map. Background Technology

[0002] Existing technologies for constructing spatial water system maps mainly rely on steps such as DEM depression filling, flow direction calculation, catchment area threshold segmentation, water system classification, and watershed extraction. However, each step has significant shortcomings: depression filling cannot distinguish between pseudo-depressions caused by data errors and real hydrological depressions such as natural lakes and wetlands; applying a globally fixed small uplift to all depressions results in poor adaptability to mountainous and plain landforms; and relying solely on elevation data without constraints from water conservancy projects such as river embankments and dams easily leads to distortion of water system morphology and misalignment between artificial and natural water systems; the multi-direction flow (MFD) algorithm's flow allocation is based solely on downslope gradient and ignores the inertia of water flow, failing to simulate the physical law that water bodies in nature cannot easily turn sharply due to the conservation of momentum, leading to… The existing methods suffer from several shortcomings. First, the current methods often deviate from the actual water catchment processes on slopes. Second, the determination of catchment area thresholds is often based on single statistical indicators such as river network density or relies on manual trial and error, lacking utilization of the two-dimensional spatial distribution information of the river network. This results in strong subjectivity and difficulty in ensuring optimal matching between the extracted river network and the actual water system. Third, in water system classification, the Shreve method treats all source rivers equally, ignoring objective differences in their cumulative runoff, gradient, catchment area, and river length. The Gravelius method requires pre-defined relationships between main streams and tributaries but lacks quantitative criteria. Fourth, in watershed extraction, there is a misalignment between the measured coordinates of hydrological stations and other outlets and the DEM grid positions. Simple adsorption positioning lacks river network constraints and elevation co-optimization, easily leading to distorted watershed boundaries. These combined deficiencies make it difficult for existing methods to generate water system spatial maps that meet the refined requirements of smart water conservancy and digital watershed analysis in terms of fidelity, accuracy, objectivity, and engineering applicability. Summary of the Invention

[0003] To address the shortcomings of existing technologies, this invention provides a method, apparatus, and medium for constructing a spatial map of a water system.

[0004] According to one aspect of the present invention, a method for constructing a spatial map of a river system is provided, comprising: The depression-filling analysis method was used to preprocess the high-precision DEM data of the water system spatial map to obtain the final corrected DEM. Based on the final corrected DEM, the total flow distribution ratio of the water system is calculated. The catchment area threshold is determined based on the total flow distribution ratio; The water system is classified based on the catchment area threshold, and the water system classification results are obtained. Based on the final corrected DEM and the total flow distribution ratio, watershed extraction and classification are performed to generate watershed vector boundary data, resulting in a complete spatial map of the water system.

[0005] According to another aspect of the present invention, an apparatus for constructing a spatial map of a water system is provided, comprising: The preprocessing module is used to preprocess high-precision DEM data of water system spatial map using depression filling analysis method to obtain the final corrected DEM; The calculation module is used to calculate the total flow distribution ratio of the water system based on the final corrected DEM; The determination module is used to determine the catchment area threshold based on the total flow allocation ratio; The classification module is used to classify water systems based on a catchment area threshold and obtain the water system classification results. The generation module is used to extract and classify watersheds based on the final corrected DEM and the total flow distribution ratio, generate watershed vector boundary data, and obtain a complete water system spatial map.

[0006] According to another aspect of the present invention, a computer-readable storage medium is provided, the storage medium storing a computer program for performing the methods described in any of the above aspects of the present invention.

[0007] According to another aspect of the present invention, an electronic device is provided, the electronic device comprising: a processor; a memory for storing executable instructions of the processor; the processor being configured to read the executable instructions from the memory and execute the instructions to implement the method described in any of the preceding aspects of the present invention.

[0008] Therefore, this invention addresses the shortcomings of existing DEM (Digital Elevation Model) depression-filling techniques, such as the inability to distinguish between pseudo-depressions and genuine hydrological depressions, poor adaptability of fixed uplift amounts to landforms, and distortion of water system morphology due to the lack of hydraulic engineering constraints. It also addresses the deficiencies of multi-flow direction algorithms, which ignore flow inertia leading to flow distribution deviations, the reliance on a single indicator for determining catchment area thresholds resulting in strong subjectivity, the lack of objective criteria for determining tributaries in water system classification, and the lack of collaborative optimization of river network elevation for hydrological station coordinate offsets in watershed extraction. The invention proposes a method for constructing a water system spatial map based on a digital elevation model. The technical solution is as follows: First, a depth-area dual-threshold discrimination mechanism is used to distinguish between pseudo-depressions and natural genuine depressions; the uplift amount is dynamically calculated using an adaptive slope formula; and high-fidelity DEM preprocessing is achieved by combining the NDWI water body index, river channel and embankment vectors, and hydraulic structure data as triple external constraints. Then, based on traditional slope weighting... The flow inertia weight is increased—the momentum of the upslope cells is used to determine the current flow direction, and the cosine of the angle between the downslope direction and the flow direction is used as the inertia weight (negative values ​​are truncated to 0)—so that the flow distribution is controlled by both slope and inertia. Then, the river network is transformed into a two-dimensional density distribution field through the core density. The fitting effect is comprehensively evaluated by the geometric weighted average of three indicators: RMSE, KGE, and MSSIM, and the catchment area threshold is automatically determined by optimization. Next, based on the cumulative runoff, the source river depth is corrected by introducing three parameters: gradient, catchment area, and river length. The relationship between the main stream and tributaries is objectively determined based on the depth at the confluence, and then Gravelius classification is performed. Finally, the outlet location is optimized under the constraint of the river network with the objective function of minimizing the sum of squares of the elevation deviation of the hydrological stations. Based on the optimized outlet, the catchment area is divided and the vector boundary is generated. This invention, through a five-stage integrated architecture, achieves precise removal of pseudo-depressions and complete preservation of natural water bodies, more natural flow direction calculation, objective and automatic threshold selection, quantitative and reproducible determination of main and tributary streams, and outlet positioning constrained by both river network and elevation. It significantly improves the fidelity, accuracy, automation, and engineering practicality of the entire process of constructing a water system spatial map. Attached Figure Description

[0009] Exemplary embodiments of the present invention can be more fully understood by referring to the following figures: Figure 1 This is a flowchart illustrating a method for constructing a spatial map of a water system provided in an exemplary embodiment of the present invention; Figure 2 This is another flowchart illustrating a method for constructing a water system spatial map provided in an exemplary embodiment of the present invention; Figure 3 This is a flowchart illustrating a high-precision DEM data preprocessing method for water system spatial maps based on depression filling analysis, provided in an exemplary embodiment of the present invention. Figure 4This is a schematic flowchart of a method for allocating water flow direction and flow rate in a spatial map of a water system provided by an exemplary embodiment of the present invention. Figure 5 This is a flowchart illustrating a method for determining the water catchment area threshold in the construction of a water system spatial map, provided by an exemplary embodiment of the present invention. Figure 6 This is a schematic diagram of the structure of a device for constructing a water system spatial map provided in an exemplary embodiment of the present invention; Figure 7 This is the structure of an electronic device provided in an exemplary embodiment of the present invention. Detailed Implementation

[0010] Hereinafter, exemplary embodiments according to the present invention will be described in detail with reference to the accompanying drawings. Obviously, the described embodiments are merely some embodiments of the present invention, and not all embodiments of the present invention. It should be understood that the present invention is not limited to the exemplary embodiments described herein.

[0011] It should be noted that, unless otherwise specifically stated, the relative arrangement, numerical expressions, and values ​​of the components and steps described in these embodiments do not limit the scope of the invention.

[0012] Those skilled in the art will understand that the terms "first," "second," etc., in the embodiments of the present invention are only used to distinguish different steps, devices, or modules, and do not represent any specific technical meaning, nor do they indicate a necessary logical order between them.

[0013] It should also be understood that in the embodiments of the present invention, "multiple" can refer to two or more, and "at least one" can refer to one, two or more.

[0014] It should also be understood that any component, data or structure mentioned in the embodiments of the present invention can generally be understood as one or more unless explicitly defined or given contrary instructions in the context.

[0015] Furthermore, the term "and / or" in this invention is merely a description of the relationship between related objects, indicating that three relationships can exist. For example, A and / or B can represent: A existing alone, A and B existing simultaneously, or B existing alone. Additionally, the character " / " in this invention generally indicates that the preceding and following related objects have an "or" relationship.

[0016] It should also be understood that the description of the various embodiments in this invention emphasizes the differences between the various embodiments, and the similarities or similarities can be referred to each other. For the sake of brevity, they will not be described in detail.

[0017] At the same time, it should be understood that, for ease of description, the dimensions of the various parts shown in the accompanying drawings are not drawn according to actual scale.

[0018] The following description of at least one exemplary embodiment is merely illustrative and is in no way intended to limit the invention or its application or use.

[0019] Techniques, methods, and equipment known to those skilled in the art may not be discussed in detail, but where appropriate, they should be considered part of the specification.

[0020] It should be noted that similar labels and letters in the following figures indicate similar items; therefore, once an item is defined in one figure, it does not need to be discussed further in subsequent figures.

[0021] The embodiments of this invention can be applied to electronic devices such as terminal devices, computer systems, and servers, and can operate together with a wide range of other general-purpose or special-purpose computing system environments or configurations. Well-known examples of terminal devices, computing systems, environments, and / or configurations suitable for use with electronic devices such as terminal devices, computer systems, and servers include, but are not limited to: personal computer systems, server computer systems, thin clients, thick clients, handheld or laptop devices, microprocessor-based systems, set-top boxes, programmable consumer electronics, network PCs, minicomputer systems, mainframe computer systems, and distributed cloud computing environments including any of the above systems, etc.

[0022] Electronic devices such as terminal devices, computer systems, and servers can be described in the general context of computer system executable instructions (such as program modules) executed by a computer system. Typically, program modules can include routines, programs, object programs, components, logic, data structures, etc., which perform specific tasks or implement specific abstract data types. Computer systems / servers can be implemented in distributed cloud computing environments, where tasks are executed by remote processing devices linked through communication networks. In distributed cloud computing environments, program modules can reside on local or remote computing system storage media, including storage devices.

[0023] Exemplary methods Figure 1 This is a schematic flowchart illustrating a method for constructing a spatial map of a water system according to an exemplary embodiment of the present invention. This embodiment can be applied to electronic devices, such as... Figure 1 As shown, the method 100 for constructing a spatial map of a water system includes the following steps: Step 101: The high-precision DEM data of the water system spatial map is preprocessed using the depression filling analysis method to obtain the final corrected DEM; Step 102: Calculate the total flow distribution ratio of the water system based on the final corrected DEM; Step 103: Determine the catchment area threshold based on the total flow distribution ratio; Step 104: Classify the water system based on the catchment area threshold to obtain the water system classification results; Step 105: Based on the final corrected DEM and the total flow distribution ratio, watershed extraction and classification are performed to generate watershed vector boundary data and obtain a complete water system spatial map.

[0024] Specifically, refer to Figure 2 As shown, the specific implementation process is as follows: This invention addresses the shortcomings of traditional DEM depression-filling formulas, such as the inability to distinguish between true and false depressions, poor adaptability to fixed small uplift amounts, and the tendency to distort water system morphology due to a lack of external data constraints. It designs an integrated optimization technology solution through threshold constraint discrimination, adaptive parameter adjustment, and multi-source data coupling control. This solution aims to eliminate DEM false depressions, ensure continuous water flow throughout the entire area, and maximize the restoration of the true topography, preserving the ecological functions of natural water systems and the original water flow patterns of artificial water conservancy projects. The overall technology is implemented step-by-step according to the process of depression identification and threshold determination → adaptive parameter dynamic calculation → multi-source data-constrained depression-filling execution. The specific technical roadmap is shown below. Figure 3 As shown.

[0025] The specific implementation steps are as follows: (1) Depression-area dual threshold discrimination mechanism Traditional depression-filling algorithms uniformly perform elevation raising operations on all depression grids, failing to distinguish between pseudo-depressions caused by data errors and genuine hydrological depressions such as naturally formed lakes, wetlands, ponds, and river bends, easily damaging the original structure of the water system. This invention adds a depression identification module and dual-threshold judgment logic to the original depression-filling formula, employing a selective depression-filling strategy to accurately distinguish between the two types of depressions and match them with different processing rules.

[0026] 1) Preliminary identification of depressions across the entire region First, the entire DEM grid dataset is scanned grid by grid, and all local depression areas are identified using 3×3 neighborhoods as the basic unit. Then, a region connectivity analysis algorithm is used to divide the spatially contiguous depression grids into independent depression units. For each independent depression unit, the following characteristic indicators are calculated: First, the maximum depth H of the depression is calculated, which is the difference between the lowest grid elevation inside the depression and the elevation of the depression's overflow outlet. Second, the contiguous area A of the depression is calculated, which is the actual surface area corresponding to all the depression grids contained in the depression.

[0027] 2) Setting the discrimination threshold and optimizing the formula system Based on the geomorphological features and hydrological control requirements of the study area, two sets of fixed discrimination thresholds were preset: H 0 To determine the minimum depth threshold of a natural depression,A 0 To determine the minimum contiguous area threshold for natural depressions, a segmented depression-filling calculation formula is constructed based on the threshold determination results: in, The grid elevation after the depression filling treatment; The original grid elevation; This is the set of elevations of the current grid's surrounding neighborhood. This represents the lowest elevation value within the current grid neighborhood. For a fixed small lift (values ​​range from 0.001 to 0.01 m); H This represents the actual maximum depth of a single depression unit; A This represents the actual contiguous area of ​​a single depression unit; Set a preset threshold for judging depression depth; The preset threshold for judging the area of ​​depressions is set.

[0028] 3) Execute gap-filling logic according to different scenarios When a depression meets either condition that its depth is greater than H0 or its contiguous area is greater than A0, the area is determined to be a real depression with hydrological and ecological functions, such as a natural lake, wetland, or pond. The original elevation of the grid is directly preserved, and no elevation lifting operation is performed on Zorig, thus preserving the natural water accumulation space and water system ecological elements in their entirety. When a depression simultaneously meets the conditions of a depth less than or equal to H0 or a contiguous area less than or equal to A0, it is determined to be a micro-pseudo-depression formed by DEM acquisition and interpolation errors. The traditional depression filling formula is then used to complete the elevation correction, eliminating the water flow obstruction problem caused by local depressions.

[0029] This mechanism only performs depression-filling operations on small-scale, shallow-depth error-type pseudo-depressions, which not only removes local depressions and ensures water flow connectivity, but also preserves the real hydraulic elements in the water system.

[0030] (2) Construct an adaptive micro-lift based on terrain slope Calculation method Traditional methods use a globally fixed value to set a small uplift ε, which cannot adapt to the diverse terrains such as mountains, hills, and plains. In plains, the terrain is gentle, and a fixed ε can easily lead to an excessively small hydraulic gradient, scattered water flow, and the emergence of parallel pseudo-channels. In mountainous areas, the terrain is undulating, and a fixed ε is insufficient to eliminate continuous deep depressions, still posing a risk of water flow interruption. This invention abandons the globally fixed value model and introduces the grid terrain slope S as a dynamic factor to construct an adaptive ε solution formula, allowing the uplift to automatically adjust according to terrain features.

[0031] 1) Definition of basic parameters and construction of formulas First, the terrain slope is calculated grid-by-grid based on the original DEM data. SThen, the minute uplift is dynamically solved using a linear function. ε Adaptable to different terrains: in, ε min The minimum lift is fixed at 0.001m, providing the basic hydraulic gradient for the entire grid. k The slope correction coefficient is pre-calibrated based on the overall landform type of the region (Table 1) to control the extent of the slope's influence on the uplift. S This represents the terrain slope value of the current DEM grid.

[0032] Table 1. Recommended values ​​for uplift of different landforms 2) Adaptive control logic based on different landforms In steep slope areas (mountainous and hilly areas), the grid slope S value is large, and the terrain itself already has sufficient natural hydraulic slope. The ε calculated by the formula only increases slightly, mainly due to the foundation uplift. While filling the depressions and opening up the water flow path, the original terrain undulations are preserved to the maximum extent to avoid excessive modification of the surface morphology. In flat areas (plains, floodplains, irrigation areas), the grid slope S approaches 0, and the natural hydraulic gradient is insufficient. The formula will simultaneously increase the value of ε to artificially construct a reasonable micro hydraulic gradient, constrain the direction of water flow, and solve problems such as disordered flow direction, disappearance of small tributaries, and straightening of river morphology in flat areas.

[0033] This method can solve the problem of insufficient adaptability of fixed parameters under different landforms, and ensure that the direction of water flow is reasonable and the connectivity is stable throughout the entire area.

[0034] (3) Multi-source data coupling correction of artificial water systems and key landforms Traditional depression-filling algorithms rely solely on DEM elevation data, neglecting constraints from water conservancy projects such as rivers, embankments, dams, and artificial canals. This can easily alter existing cross-sections and deviate from the main flow line, leading to misalignment between artificial and natural water systems. This invention integrates three types of multi-source spatial data: water body remote sensing data, water system vector data, and thematic data on hydraulic structures. It establishes multi-level constraint rules, embedding external control boundaries into the entire depression-filling process, achieving precise matching between DEM correction and measured water system and water conservancy project morphology.

[0035] 1) Delineate the entire area as a no-fill zone based on water quality index. The Normalized Difference Water Index (NDWI) was extracted from remote sensing imagery to automatically identify and vectorize the main channels of perennially water-bearing rivers, the main bodies of lakes, and large wetlands in the study area, generating a permanent water body extent layer. Before the depression-filling calculation, all DEM grids within this layer were marked as no-fill grids. During the algorithm's traversal calculation, no no-filling areas were skipped, and no elevation lifting operations were performed, strictly protecting the perennial water space and preventing the problems of water area reduction and water body morphology changes caused by depression filling.

[0036] 2) Establish elevation constraint lines based on river channel and embankment vectors. Import measured vector elements such as the river centerline, estuary line, and levee edge line as topographic fault constraint lines. For the river flood passage section and levee protection area within the constraint line range, add elevation change limits. Strictly control the grid elevation rise during the filling process, prohibiting significant increases in the elevation of the river bottom and levee toe, ensuring that the designed river cross-section and levee cross-section dimensions remain unchanged, and guaranteeing that the original functions of flood control, water conveyance, and drainage, as well as the water flow path, do not deviate.

[0037] 3) Add boundary control for hydraulic structures such as sluices, dams, and canals. Thematic vector data for sluice gates, dams, irrigation and drainage canals, and structures crossing dikes are extracted to delineate hydraulic structures and their impact ranges. Rigid elevation boundary constraints are then applied within these areas. On the one hand, arbitrary filling of depressions is prohibited on the structures themselves and their upstream and downstream connections to avoid altering the sluice gate / dam opening and closing water levels and the longitudinal slope of the canal bottom. On the other hand, for linear artificial water systems such as main canals, branch canals, and field ditches, filling and correction of depressions are strictly carried out according to the original slope and cross-sectional shape of the canal system to ensure that the connectivity and water conveyance capacity of the artificial canal network are consistent with the design standards.

[0038] The multi-constraint rules primarily target the entire depression-filling process of the initial DEM correction for precise correction and boundary control, effectively improving the shortcomings of traditional depression-filling methods that are unconstrained and result in uniform elevation across the entire area. Three types of differentiated constraints achieve refined terrain correction: by shielding the dynamic elevation calculation of water body grids, the problem of conventional algorithms mistakenly raising the elevation of natural water bodies such as rivers and lakes is avoided, ensuring the original topography of natural water bodies and preventing distortion of water system topography; river and levee fracture constraints, based on topographic fracture constraint lines, delineate the boundaries of key flood control areas, strictly control the upper limit of grid elevation within the area, solve the problem of excessive depression filling, effectively avoid defects such as the disappearance of river slope and weakening of levee topography, and stabilize the original trend of the water system and the cross-sectional shape of flood control; rigid boundary constraints for hydraulic structures, by delineating rigid elevation boundaries, strictly limit the correction range and magnitude of the DEM around water conservancy facilities such as dams, canals, and structures penetrating levees, preserving the original topographic structure of artificial water conservancy facilities. The overall process of correcting pseudo-depressions is precisely controlled by zoning and grading, eliminating topographic distortions caused by unconstrained depression filling, and fully preserving the true topographic features of natural depressions, water systems and artificial water conservancy facilities, thus achieving high-precision DEM depression filling correction that conforms to the actual water system topography.

[0039] Therefore, this invention addresses three major shortcomings of traditional DEM depression-filling algorithms: inability to distinguish depression types, weak adaptability to fixed-parameter terrain, and lack of multi-source data constraints. It employs an integrated optimization scheme that combines depth-area dual-threshold discrimination, slope adaptive parameter calculation, and multi-source spatial data coupling constraints. This completes a technological upgrade across the entire process, from depression identification and parameter adjustment to boundary control. Significant technical effects are achieved in areas such as depression removal, water flow connectivity, terrain fidelity preservation, and water system function preservation, comprehensively improving the accuracy of DEM preprocessing and the reliability and practicality of water system spatial maps. Detailed technical effects are described below: (1) Implement classified treatment of depressions, taking into account both the removal of pseudo-depressions and the complete preservation of natural water systems and ecological spaces. Traditional depression-filling algorithms rely solely on elevation relationships to determine depression areas, applying indiscriminate elevation increases to both pseudo-depressions and natural depressions. While eliminating data-related depressions, they also fill in hydrological and ecological units within the watershed, such as natural lakes, ecological wetlands, small ponds and dams, and river bends that serve to regulate floods, conserve water resources, and maintain biodiversity. This directly leads to damage to the ecological structure of the water system spatial map, loss of key water system elements, and a significant reduction in the integrity and authenticity of the map data.

[0040] Furthermore, this invention introduces a dual-threshold discrimination mechanism for depression depth and contiguous area. First, it performs depression traversal identification and characteristic index calculation on the entire DEM to achieve automated differentiation of depression types. For shallow, small-scale pseudo-depressions caused by surveying errors, data interpolation, and terrain noise, normal depression-filling operations are performed to eliminate local "stagnant water" areas, addressing the problems of water flow obstruction and confluence path breakage caused by depressions at their source. For natural waterlogged depressions with depths and areas exceeding the set thresholds, the original grid elevation is retained without further elevation correction. This method not only achieves the basic objectives of DEM preprocessing but also fully preserves the original water bodies, wetlands, and water-retaining spaces, effectively protecting the regional natural water system pattern and ecological structure. This ensures that the final generated water system spatial map objectively reflects the true surface hydrological features, fundamentally improving the integrity and hydrological rationality of the map information.

[0041] (2) Parameters are dynamically and adaptively adjusted to adapt to complex terrain across the entire area, ensuring continuous water flow and a natural and reasonable river network morphology. Traditional algorithms use a globally fixed, small lift amount ε without differentiating it based on terrain slope and geomorphological features. This reveals significant shortcomings in different terrain regions: In plains, where the terrain is flat, a fixed ε can easily result in an excessively small hydraulic gradient, leading to scattered water flow direction and path deviation, which in turn can cause problems such as parallel pseudo-channels, straightening of main channels, and disappearance of small tributaries. In mountainous and hilly areas, where the terrain is undulating and there are many contiguous deep depressions, a fixed lift amount cannot meet the requirements for filling depressions. Incomplete depression removal still poses the risk of water flow interruption and confluence failure, seriously affecting the accuracy of hydrological analysis and river network extraction results.

[0042] This invention constructs an ε-adaptive calculation model based on grid topographic slope, enabling minute uplift amounts to dynamically change with the landform, achieving precise matching between parameters and geomorphological features. In steep slopes and hilly areas, where the terrain itself possesses sufficient natural hydraulic slope, the algorithm automatically reduces the ε value to correct local depressions with minimal impact, preserving the original topographic undulations to the maximum extent and avoiding excessive modification of the landform. In gentle areas such as plains, floodplains, and irrigated areas, where the natural hydraulic slope is insufficient, the algorithm synchronously and appropriately increases the ε value, artificially constructing a stable and orderly hydraulic slope to constrain the direction of water flow, effectively preventing phenomena such as turbulent flow, pseudo-channel formation, and tributary absence. After optimization, a continuous, unidirectional, and logically sound confluence system can be formed throughout the entire watershed, regardless of whether it is mountainous, hilly, or plain. Water flow connectivity fully meets the standards, and the automatically extracted river network direction and branch structure highly match the actual terrain, making the river network morphology closer to natural reality.

[0043] (3) Introduce rigid constraints from multi-source data to protect artificial water conservancy projects and achieve seamless connection between natural and artificial water systems. Traditional depression-filling algorithms rely solely on single DEM elevation data for calculations, without combining measured water conservancy thematic vector data such as river channels, dikes, sluice gates, and irrigation and drainage systems for boundary control. The elevation correction process is not limited by engineering boundaries, which can easily alter the cross-section of dikes, the cross-section of river channels, the longitudinal slope of canals, and the original flow patterns upstream and downstream of sluice gates. This disrupts the design and operation pattern of water conservancy projects, ultimately leading to misalignment of the spatial location and flow direction of artificial and natural water systems in the water system spatial map. This fails to meet the actual business needs of river and lake management, canal network scheduling, and flood control planning.

[0044] This invention integrates multi-source spatial data, including measured water system vectors, thematic data on hydraulic structures, and remote sensing water indices, to establish multi-level rigid constraints throughout the entire process of filling depressions. On the one hand, water indices are used to delineate areas where filling is prohibited due to year-round water accumulation, strictly protecting the scope and morphology of natural water bodies. On the other hand, based on river, embankment, and canal system vectors, as well as thematic data on sluice gates, dams, and river-blocking structures, elevation control lines and management boundaries are delineated to limit the range and magnitude of elevation rise within the water conservancy project and its impact area, strictly prohibiting changes to the original water passage cross-section, canal bottom slope, and designed water flow path. This constraint mechanism embeds water conservancy project management rules and actual water system morphology into the DEM preprocessing stage, effectively avoiding the distortion of artificial water system morphology caused by filling depressions. It ensures that the original functions and water flow patterns of water conservancy facilities such as dikes, dams, irrigation and drainage canals remain unchanged, and allows the processed DEM data and water system spatial map to match the distribution of natural water systems and the layout of artificial water conservancy projects. This achieves the coordination and unity of the two types of water system spatial morphology and water flow logic, and the data results can be directly applied to water conservancy operations such as river management, canal system operation and maintenance, and flood control scheduling.

[0045] (4) Overall application effect and overall value In summary, this invention effectively addresses many pain points of traditional iterative depression-filling algorithms in water system spatial map construction and watershed hydrological analysis through the combined application of three optimization techniques. The entire technical solution retains the advantages of traditional algorithms, such as simple logic, stable operation, and strong versatility, while also overcoming their shortcomings, such as insufficient depression identification ability, poor terrain adaptability, and lack of engineering constraints.

[0046] At the data level, it achieves triple effects of terrain fidelity, water body shape preservation, and engineering state preservation, significantly improving the accuracy of DEM data and the realism of water system maps. At the functional level, it simultaneously meets multiple requirements such as depression removal, water flow connectivity, ecological protection, and engineering management, and is adaptable to various water conservancy scenarios such as natural watersheds, plain river networks, artificial irrigation areas, and flood storage and detention areas. At the application level, the processed results can stably support a series of water conservancy GIS tasks such as river network extraction, confluence simulation, urban flooding risk assessment, river and lake shoreline planning, and soil and water conservation analysis. The technology's applicability and versatility are greatly enhanced, enabling it to better serve the industry development needs of smart water conservancy, water system digitization, and refined management of river and lake spaces.

[0047] This invention improves upon the existing classic multi-flow direction algorithm (MFD) by considering the influence of water flow inertia on flow distribution. In addition to the slope factor of the downhill cells, it adds the angle factor between the direction of the downhill cell and the current cell's water flow direction. By comprehensively considering the effects of slope and water flow inertia, the flow distribution more closely aligns with the physical water flow process in nature. (Reference) Figure 4 As shown, the specific implementation steps are as follows: 1. Calculation of the first flow distribution ratio under the influence of slope factor In the 3x3 DEM grid formed by the current cell and its adjacent cells, the eight adjacent cells are numbered 1-8, starting with the cell facing east (numbered 1) and increasing clockwise. The calculation of the first flow allocation ratio under the influence of slope follows the classic multi-direction flow (MFD) algorithm. The flow allocation ratio for each downhill cell is related to the downhill slope of that cell. The specific calculation method is as follows: Where d gr (j) represents the proportion of the flow distribution from the downhill cell to the j-th adjacent cell due to the slope factor; tanβ is the slope; p is the flow distribution weight, taken as 1.1; L j L is the weighting factor for the contour length of the j-th adjacent cell. j =0.5, when j=1, 3, 5, 7 and the elevation is lower than the current cell; L j = When j=2, 4, 6, 8 and the elevation is lower than the current cell; L j =0, when the elevation is not lower than the current cell; 2. Calculation of the second flow distribution ratio under the influence of water flow inertia. 2.1 Calculation of water flow direction in the current cell Establish a coordinate system with the center point of the current cell as the origin, with east as the positive x-axis and north as the positive y-axis. The azimuth of the vector direction starts at north and increases clockwise.

[0048] The current cell's water flow direction is the vector sum of the water flow from adjacent uphill cells flowing into the current cell. In this invention, the flow rate of all adjacent uphill cells is uniformly set to 1, and the flow velocity is proportional to the square root of the slope of the uphill cell according to the Manning formula.

[0049] The projections of the total water flow in the current cell onto the x-axis and y-axis are as follows: Projection of water flow along the x-axis Projection of water flow along the y-axis Where tanβ j This represents the slope of the adjacent cell.

[0050] K j Used to determine whether the adjacent cell is an uphill cell or a downhill cell. For uphill cells, K j =1; for downhill cells, K j =0.

[0051] D u (j) represents the azimuth angle from the adjacent cell to the current cell, D u (j) = 1.5π, when j = 1; D u (j) = 1.75π, when j = 2; D u (j)=0, when j=3; D u (j) = 0.25π, when j = 4; D u (j) = 0.5π, when j = 5; D u (j) = 0.75π, when j = 6; D u (j) = π, when j = 7; D u (j) = 1.25π, when j = 8.

[0052] Based on the projection of water flow rate in the x and y directions, pw x ,pw y Calculate the azimuth angle α of the water flow in the current cell, with a value in the range [0, 2π).

[0053] when pw y >0 and pw x ≥0, α=arctan(| pw x | / | pw y |); when pw y >0 and pw x <0, α=2π-arctan(| pw x | / | pw y |); when pw y <0 and pw x ≥0, α=π-arctan(| pw x | / | pwy |); when pw y <0 and pw x <0, α=π+arctan(| pw x | / | pw y |); when pw y =0 and pw x >0, α=π / 2; when pw y =0 and pw x <0, α=3π / 2.

[0054] 2.2 Calculation of Downhill Cell Direction The direction of the downhill cell is the direction in which the current cell flows into the adjacent cell, and its azimuth angle D is... d (j) Determined based on the numbers of adjacent cells, D d (j) = 0.5π, when j = 1; D d (j) = 0.75π, when j = 2; D d (j) = π, when j = 3; D d (j) = 1.25π, when j = 4; D d (j) = 1.5π, when j = 5; D d (j) = 1.75π, when j = 6; D d (j)=0, when j=7; D d (j) = 0.25π, when j = 8.

[0055] 2.3 Calculation of included angle The angle γ(j) between the downhill cell direction and the current cell's water flow direction is obtained by subtracting their azimuth angles: γ(j)=D d (j)-α The value range of γ(j) is (-2π, 1.75π).

[0056] 2.4 Calculation of Water Flow Inertia Weight The inertial weight of the water flow in each downhill cell is the cosine of the angle between the downhill cell direction and the current cell's water flow direction. Since the angle ranges from -2π to 1.75π, the cosine value can be negative. However, negative values ​​have no meaning in the actual physical water flow process, so the negative values ​​are replaced with 0. That is, downhill cells located in the direction perpendicular to or opposite to the current cell's water flow direction do not participate in flow distribution.

[0057] The calculation method for the water flow inertia weight W(j) of each downhill cell is as follows: W(j)=max(cos(γ(j))*P j ,0) Where P j To determine whether an adjacent cell is an uphill or downhill cell, for an uphill cell, P j =0, for downhill cells, P j =1.

[0058] 2.5 Calculation of the Second Flow Allocation Ratio The proportion of flow allocated to each downhill cell under the influence of water flow inertia is related to the angle between the direction of that downhill cell and the direction of the current cell's flow. The flow allocation proportion d for each downhill cell... in (j) The calculation method is as follows: 3. Total flow allocation ratio Under the combined influence of slope and flow inertia, the total flow distribution ratio allocated to each downhill cell is the first flow distribution ratio d under the influence of the downhill cell slope factor. gr (j) The second flow distribution ratio d under the influence of water flow inertia factors in The product of (j).

[0059] The total flow distribution ratio d of each downhill cell t (j) The calculation method is as follows: Therefore, in existing multi-flow direction algorithms (MFD), the MFD-md algorithm (also known as the MFD-fg algorithm) considers the influence of uphill cells on flow distribution. However, this algorithm comprehensively uses data from all adjacent cells, including the uphill cells, to calculate the local maximum downhill slope in the 3*3 structure DEM grid, which is used to improve the flow distribution weight parameter p, thereby affecting the flow distribution in each downhill cell. This algorithm does not directly consider the factor of flow inertia. This invention, in terms of downhill cell slope factors, adopts the classic multi-flow direction algorithm (MFD). The advantage of this invention is that it calculates the flow inertia factor separately and then combines it with the downhill cell slope factor, resulting in a flow distribution calculation method that simultaneously reflects the influence of both downhill cell slope and flow inertia.

[0060] Furthermore, many existing methods for comparing extracted river network data with actual water system data (blue lines) are fundamentally based on river network density analysis. For example, the river network density method treats river network density as a constant, while the watershed width distribution method uses river network density as a one-dimensional mathematical description; however, river networks are distributed in two dimensions in actual watersheds, and current methods have not fully utilized the distribution pattern of river network density in a two-dimensional plane. This invention's technique for determining the catchment area threshold in water system spatial mapping constructs a system by numerically representing the core density parameter and providing a two-dimensional description of the spatial distribution of river network density in the corresponding watershed. A core density fitting index, obtained by geometrically weighting the root mean square error, Kling-Gupta efficiency coefficient, and average structural similarity index, is used to verify the fitting effect of the extracted river network core density distribution on the actual watershed core density distribution. The catchment area threshold with the best fitting effect is then defined as the optimal catchment area threshold. Figure 5 As shown, the specific steps are as follows: 1. Water system extraction with different catchment area thresholds Based on the actual natural geographical conditions of the watershed, a lower limit Amin, an upper limit Amax, and a step size ΔA are set for the catchment area threshold, resulting in a list of catchment area thresholds A1, A2, A3...An. Different catchment area thresholds from A1 to An are compared with the cumulative runoff volume of each cell in the DEM. If the cumulative runoff volume is greater than or equal to the catchment area threshold, the cell is classified as a water system; if the cumulative runoff volume is less than the catchment area threshold, the cell is classified as a slope.

[0061] Generate river network links for all cells identified as water systems, convert them into vector line features, and extract the water systems under different catchment area thresholds.

[0062] 2. Matching actual water system data The actual water system data is converted into vector line features using geographic information system tools and georeferenced with the DEM data used to extract the water system. This ensures that the actual water system data is correctly matched with the DEM data and the extracted water system, facilitating the calculation of the water system core density D in the next step.

[0063] 3-line kernel density calculation This section calculates the river network core density (D) for each DEM cell in the watershed under different catchment area thresholds and actual river network conditions. Core density is an indicator that characterizes the density of a certain type of line feature in a local DEM cell. River network core density can describe the spatial distribution of river network density in the watershed in a two-dimensional way. The calculation method for river network core density (D) is as follows: D: Core density of the river system; R: Search radius, taken as 1.5 times the average channel spacing; l i : Water system features covered within the search radius of the cell center; Weight i The weight field for the line is set to 1 for all values. d i ( s ): Distance from the center of the cell to the smallest arc segment in the water system line element; K(u): Double square kernel function.

[0064] In the formula, u As an intermediate variable in the composite function, u= .

[0065] 4. Calculation of Fitting Performance Index The fitting effect of the extracted core density of the water system on the actual core density of the water system was examined using the core density fitting index G. The core density fitting index G was calculated from the root mean square error (RMSE), the Kling-Gupta efficiency coefficient (KGE), and the mean structural similarity index (MSSIM).

[0066] The core density of the water system extracted under a certain catchment area threshold is: D m ( x , y The actual core density of the water system is... D z ( x , y The calculation methods for each fitting effect index are as follows: (1) The root mean square error (RMSE) is calculated as follows. The smaller the value, the better the fit: In the formula, M Total number of cells in the grid, which is the number of rows multiplied by the number of columns; x i ,y i : No. i The spatial coordinates of each cell;D z ( x i , y i ): The core density of the actual water system at this cell location; D m ( x i , y i ): This cell contains the core density of the water system.

[0067] (2) The Kling-Gupta efficiency coefficient (KGE) is calculated as follows, and its value ranges from -∞ to 1. The closer it is to 1, the better the fitting effect is: r Correlation coefficient; α : Variability ratio; β : Deviation ratio; where the correlation coefficient r is calculated as follows: The method for calculating the variability ratio α is as follows: The deviation ratio β is calculated as follows: The method for calculating the mean is as follows: In the formula, M Total number of cells in the grid, which is the number of rows multiplied by the number of columns; Dz : The core density of the actual water system at each cell; : The average core density of the actual water system; Dm Extract the core density of the water system at each cell; Extract the average linear core density of the water system; σ Dz Standard deviation of actual hydrographic core density; σ Dm Extract the standard deviation of the core density of the water system line; (3) The mean structural similarity index (MSSIM) is calculated as follows, with a value range of -1 to 1. The closer the value is to 1, the better the fit. In the formula, M Total number of cells in the raster, calculated as number of rows × number of columns; SSIM (window j ): The structural similarity index of a given cell.

[0068] Let z be the core density of the actual water system in a cell, and m be the core density of the extracted water system in a cell. The Structural Similarity Index (SSIM) of this cell is calculated as follows: In the formula, N : The cells within an 11×11 area surrounding the center cell; μ z The average core density of the actual water system within a 11×11 range; μ m Extract the average core density of the water system from cells within an 11×11 range; σ z : The variance of the core density of the actual water system within a 11×11 range; σ m Extract the core density variance of the water system from cells within an 11×11 range; σ zm The covariance between the actual core density of the water system and the extracted core density of the water system within a 11×11 range; z i The core density of the actual water system in each cell within an 11×11 range; m i Extract the core density of the water system from each cell within an 11×11 range; L C1: The difference between the maximum and minimum values ​​of the linear core density of the water system; K1: a constant; K2: a constant; K2: a constant.

[0069] (4) The fitting effect of the extracted core density of the water system on the actual core density of the water system is tested by the core density fitting index G. The core density fitting index G is the geometrically weighted average of the root mean square error (RMSE), the Kling-Gupta efficiency coefficient (KGE), and the average structural similarity index (MSSIM). The calculation method is as follows, and its value ranges from 0 to 1. The closer it is to 1, the better the fitting effect is: Among them, RMSE norm The normalized root mean square error: In the formula, RMSE: original root mean square error; σ Dz: Standard deviation of the actual core density of the water system.

[0070] Where KGE' is the truncated Kling-Gupta efficiency coefficient: KGE'=max(0,KGE) Where MSSIM' is the truncated average structural similarity index: MSSIM'=max(0,MSSIM) 5. Determination of the optimal catchment area threshold The line kernel density fitting index G was used to analyze the fitting effect of the water system extracted by different catchment area thresholds. The closer the line kernel density fitting index G is to 1, the better the fitting effect. The catchment area threshold with the best fitting effect was defined as the optimal catchment area threshold.

[0071] Therefore, this invention describes the spatial distribution of river network density in the corresponding watershed by numericalizing the core density parameter and using two dimensions. The core density fitting index, obtained by geometrically weighting the root mean square error, Kling-Gupta efficiency coefficient, and average structural similarity index, is used to test the fitting effect of the extracted core density distribution of the river network on the actual watershed core density distribution. The catchment area threshold with the best fitting effect is defined as the optimal catchment area threshold.

[0072] Furthermore, the water system spatial map water system classification technology of this invention improves the existing Gravelius classification method by coupling the Shreve classification method: based on parameters such as cumulative runoff, gradient, drainage area, and river length, the source river level in the Shreve classification method is improved to a depth value. The quantitative judgment of the relationship between tributaries and main streams is realized based on the cumulative river depth of each river segment, which serves as the basis for Gravelius classification, making the Gravelius classification results more objective and consistent.

[0073] This invention relates to a water system spatial map water system classification technology, which improves upon the existing Gravelius classification method by coupling it with the Shreve classification method. Based on parameters such as runoff accumulation, gradient, drainage area, and river length, the source river level in the Shreve classification method is modified to a depth value. The quantitative judgment of the relationship between tributaries and main streams is achieved based on the accumulated river depth of each river segment, serving as the basis for Gravelius classification and making the Gravelius classification results more objective and consistent.

[0074] 1. Calculation of the depth of the source river The original Shreve method defines the smallest river without any inflowing rivers as a Class 1 river. When calculating the class of a new river formed by the confluence of two rivers, the classes of the inflowing rivers are added together. However, this method ignores the differences among Class 1 rivers in terms of cumulative runoff, gradient, drainage area, and river length.

[0075] In the process of obtaining a river system spatial map through topographic analysis based on regular grid DEM data, the runoff accumulation FA, gradient J, drainage area A, and river length L of each first-order river are easily obtained river parameters. This invention defines depth D as a correction to the runoff accumulation based on gradient, drainage area, and river length. Depth is positively correlated with these three parameters; that is, the larger the gradient, drainage area, and river length of a first-order river, the greater its depth value.

[0076] Assuming there are m Class 1 rivers in the water network, the depth D of each Class 1 river is calculated as follows: 2. Depth calculation for each river section The depth of each downstream river section is calculated using the Shreve method, which is the sum of the depths of the rivers flowing into them. This process is repeated until the depth of all river sections has been calculated.

[0077] 3. Determining the relationship between main streams and tributaries For each river confluence, the relationship between the main stream and tributaries is determined based on the depth of the rivers flowing into them; the deeper river is designated as the main stream, and the shallower river as the tributary. This method establishes the main stream and tributary relationships of the complete river system.

[0078] 4. Gravelius method for river classification Once the relationship between the main stream and tributaries of a river system is clarified, the Gravelius method is used to classify the rivers. The largest main stream in the basin is defined as level 0, the secondary streams flowing into the main stream are defined as level 1, the tributaries flowing into level 1 main streams are defined as level 2, and so on, until all river systems are defined.

[0079] Furthermore, this invention provides a method for watershed extraction and classification based on hydrological station constraints in spatial maps of river systems, specifically including: Step S1: Data Acquisition. Acquire digital elevation model (DEM) data, flow direction raster data, runoff accumulation raster data, and location coordinate data of known hydrological stations for the study area.

[0080] The water flow direction grid data and the confluence accumulation grid data can be obtained in advance through the following methods: filling depressions in the DEM data to obtain a depression-free DEM; calculating the water flow direction of each grid cell using the D8 single-flow direction algorithm based on the depression-free DEM to generate water flow direction grid data; and counting the number of upstream grid cells flowing into each grid cell based on the water flow direction data to generate confluence accumulation grid data. The aforementioned depression filling, water flow direction calculation, and confluence accumulation calculation are all well-known techniques in the art, and this invention does not limit them.

[0081] Step S2: River network extraction. Based on the runoff accumulation raster data and a threshold, extract the river network raster data. Mark rasters with runoff accumulation greater than or equal to the threshold as river rasters, and rasters with runoff accumulation less than the threshold as non-river rasters, generating binarized river network raster data.

[0082] Step S3: Determine the location of the outlet based on the location constraints of the hydrological station. This is the core step of the invention. For each known hydrological station location coordinate... Determine the search window centered on it. In the search window Internally, construct the objective function for optimizing the outlet positioning: in, For search window The set of all candidate raster cells located on the river network raster data; Candidate raster cells Elevation value; Location of the hydrological station The elevation value at the location; This represents the number of candidate raster cells within the search window.

[0083] The physical meaning of the above objective function is: to search the river network for the grid cell whose elevation is closest to the measured location of the hydrological station, and determine it as the optimal outlet location. This objective function simultaneously satisfies two constraints: (1) The outlet must be located on the river network; (2) The sum of squares of the deviations between the elevation of the outlet and the measured elevation of the hydrological station is the smallest.

[0084] When there are multiple hydrological stations in the study area, the above-mentioned optimized positioning process is performed on each hydrological station to obtain multiple sets of outlet locations. .

[0085] Step S4: Watershed delineation based on the determined outlet location. Starting from the outlet location determined in Step S3, perform reverse tracing of the raster data along the flow direction: recursively trace all upstream raster cells that ultimately converge at the outlet. Mark all traced raster cells as belonging to the same watershed, thus obtaining the watershed delineation raster data. : When there are multiple outlets, the following division method is used: For each grid cell, calculate its distance from the water flow path to each outlet; The grid cell is assigned to the watershed corresponding to the outlet whose flow path is closest to it; When a grid cell belongs to the confluence range of multiple outlets, the attribution is determined according to the principle of prioritizing the shortest flow path distance. When the distance of the water flow path is the same, the decision is made according to the principle of lowest elevation priority.

[0086] Step S5: Watershed Classification and Identification. All watershed units obtained in Step S4 are classified, and a unique identifier (ID) is generated for each unit. Upstream and downstream hierarchical relationships between watershed units are also encoded. Each watershed unit contains the following attribute information: unique ID, watershed level, list of upstream watershed IDs, list of downstream watershed IDs, outlet coordinates, and watershed area.

[0087] Step S6: Watershed Vector Boundary Generation. Extract the outer boundary of each watershed unit from Step S4, convert the raster-formatted watershed boundary into a vector polygon, and generate watershed vector boundary data. : For shared boundaries between adjacent watershed units, a topological consistency processing method is adopted to ensure that the boundaries of adjacent watersheds completely coincide, without gaps or overlaps. The final generated hydrological spatial map contains complete elements such as watershed vector boundaries, watershed hierarchical information, and watershed topological relationships.

[0088] Therefore, this invention addresses the shortcomings of existing DEM (Digital Elevation Model) depression-filling techniques, such as the inability to distinguish between pseudo-depressions and genuine hydrological depressions, poor adaptability of fixed uplift amounts to landforms, and distortion of water system morphology due to the lack of hydraulic engineering constraints. It also addresses the deficiencies of multi-flow direction algorithms, which ignore flow inertia leading to flow distribution deviations, the reliance on a single indicator for determining catchment area thresholds resulting in strong subjectivity, the lack of objective criteria for determining tributaries in water system classification, and the lack of collaborative optimization of river network elevation for hydrological station coordinate offsets in watershed extraction. The invention proposes a method for constructing a water system spatial map based on a digital elevation model. The technical solution is as follows: First, a depth-area dual-threshold discrimination mechanism is used to distinguish between pseudo-depressions and natural genuine depressions; the uplift amount is dynamically calculated using an adaptive slope formula; and high-fidelity DEM preprocessing is achieved by combining the NDWI water body index, river channel and embankment vectors, and hydraulic structure data as triple external constraints. Then, based on traditional slope weighting... The flow inertia weight is increased—the momentum of the upslope cells is used to determine the current flow direction, and the cosine of the angle between the downslope direction and the flow direction is used as the inertia weight (negative values ​​are truncated to 0)—so that the flow distribution is controlled by both slope and inertia. Then, the river network is transformed into a two-dimensional density distribution field through the core density. The fitting effect is comprehensively evaluated by the geometric weighted average of three indicators: RMSE, KGE, and MSSIM, and the catchment area threshold is automatically determined by optimization. Next, based on the cumulative runoff, the source river depth is corrected by introducing three parameters: gradient, catchment area, and river length. The relationship between the main stream and tributaries is objectively determined based on the depth at the confluence, and then Gravelius classification is performed. Finally, the outlet location is optimized under the constraint of the river network with the objective function of minimizing the sum of squares of the elevation deviation of the hydrological stations. Based on the optimized outlet, the catchment area is divided and the vector boundary is generated. This invention, through a five-stage integrated architecture, achieves precise removal of pseudo-depressions and complete preservation of natural water bodies, more natural flow direction calculation, objective and automatic threshold selection, quantitative and reproducible determination of main and tributary streams, and outlet positioning constrained by both river network and elevation. It significantly improves the fidelity, accuracy, automation, and engineering practicality of the entire process of constructing a water system spatial map.

[0089] Exemplary device Figure 6 This is a schematic diagram of the structure of a device provided in an exemplary embodiment of the present invention. For example... Figure 6 As shown, the device 600 includes: Preprocessing module 610 is used to preprocess high-precision DEM data of water system spatial map using depression filling analysis method to obtain the final corrected DEM; Calculation module 620 is used to calculate the total flow distribution ratio of the water system based on the final corrected DEM; Module 630 is used to determine the catchment area threshold based on the total flow distribution ratio; The classification module 640 is used to classify water systems based on a catchment area threshold and obtain water system classification results. The generation module 650 is used to extract and classify watersheds based on the final corrected DEM and the total flow distribution ratio, generate watershed vector boundary data, and obtain a complete water system spatial map.

[0090] Exemplary electronic devices Figure 7 This is the structure of an electronic device provided in an exemplary embodiment of the present invention. For example... Figure 7 As shown, the electronic device 70 includes one or more processors 71 and a memory 72.

[0091] The processor 71 may be a central processing unit (CPU) or other form of processing unit with data processing capabilities and / or instruction execution capabilities, and may control other components in the electronic device to perform desired functions.

[0092] The memory 72 may include one or more computer program products, which may include various forms of computer-readable storage media, such as volatile memory and / or non-volatile memory. The volatile memory may include, for example, random access memory (RAM) and / or cache memory. The non-volatile memory may include, for example, read-only memory (ROM), hard disk, flash memory, etc. One or more computer program instructions may be stored on the computer-readable storage medium, and the processor 71 may execute the program instructions to implement the methods of the software programs of the various embodiments of the present invention described above, and / or other desired functions. In one example, the electronic device may also include an input device 73 and an output device 74, these components being interconnected via a bus system and / or other forms of connection mechanisms (not shown).

[0093] In addition, the input device 73 may also include, for example, a keyboard, a mouse, etc.

[0094] The output device 74 can output various information to the outside. The output device 74 may include, for example, a display, a speaker, a printer, and a communication network and its connected remote output devices, etc.

[0095] Of course, for the sake of simplicity, Figure 7 Only some of the components of this electronic device relevant to the present invention are shown, omitting components such as buses, input / output interfaces, etc. In addition, the electronic device may include any other suitable components depending on the specific application.

[0096] Exemplary computer program products and computer-readable storage media In addition to the methods and apparatus described above, embodiments of the present invention may also be computer program products, which include computer program instructions that, when executed by a processor, cause the processor to perform the steps in the methods according to various embodiments of the present invention described in the "Exemplary Methods" section above.

[0097] The computer program product can be written in any combination of one or more programming languages ​​to perform the operations of the embodiments of the present invention. The programming languages ​​include object-oriented programming languages ​​such as Java and C++, as well as conventional procedural programming languages ​​such as C or similar languages. The program code can be executed entirely on the user's computing device, partially on the user's computing device, as a standalone software package, partially on the user's computing device and partially on a remote computing device, or entirely on a remote computing device or server.

[0098] Furthermore, embodiments of the present invention may also be computer-readable storage media storing computer program instructions thereon, which, when executed by a processor, cause the processor to perform the steps of the methods according to various embodiments of the present invention described in the "Exemplary Methods" section above.

[0099] The computer-readable storage medium may be any combination of one or more readable media. A readable medium may be a readable signal medium or a readable storage medium. A readable storage medium may be, for example, an electrical, magnetic, optical, electromagnetic, infrared, or semiconductor system, system, or device, or any combination thereof. More specific examples of readable storage media (a non-exhaustive list) include: an electrical connection having one or more wires, a portable disk, a hard disk, random access memory (RAM), read-only memory (ROM), erasable programmable read-only memory (EPROM or flash memory), optical fiber, portable compact disk read-only memory (CD-ROM), optical storage devices, magnetic storage devices, or any suitable combination thereof.

[0100] The basic principles of the present invention have been described above with reference to specific embodiments. However, it should be noted that the advantages, benefits, and effects mentioned in the present invention are merely examples and not limitations, and should not be considered as essential features of each embodiment of the present invention. Furthermore, the specific details disclosed above are for illustrative and facilitative purposes only, and are not limitations. These details do not limit the present invention to the necessity of employing the aforementioned specific details.

[0101] The various embodiments in this specification are described in a progressive manner, with each embodiment focusing on its differences from other embodiments. Similar or identical parts between embodiments can be referred to interchangeably. For system embodiments, since they largely correspond to method embodiments, the description is relatively simple; relevant parts can be referred to the descriptions in the method embodiments.

[0102] The block diagrams of devices, systems, devices, and systems involved in this invention are merely illustrative examples and are not intended to require or imply that they must be connected, arranged, or configured in the manner shown in the block diagrams. As those skilled in the art will recognize, these devices, systems, devices, and systems can be connected, arranged, and configured in any manner. Words such as “comprising,” “including,” “having,” etc., are open-ended terms meaning “including but not limited to,” and are used interchangeably with them. The terms “or” and “and” as used herein refer to the terms “and / or,” and are used interchangeably with them unless the context clearly indicates otherwise. The term “such as” as used herein refers to the phrase “such as but not limited to,” and is used interchangeably with it.

[0103] The methods and systems of the present invention may be implemented in many ways. For example, they may be implemented by software, hardware, firmware, or any combination of software, hardware, and firmware. The above-described order of steps for the methods is for illustrative purposes only, and the steps of the methods of the present invention are not limited to the order specifically described above unless otherwise specifically stated. Furthermore, in some embodiments, the present invention may also be implemented as a program recorded on a recording medium, the program comprising machine-readable instructions for implementing the methods according to the present invention. Thus, the present invention also covers recording media storing programs for performing the methods according to the present invention.

[0104] It should also be noted that in the systems, apparatus, and methods of the present invention, the components or steps can be disassembled and / or recombined. These disassemblies and / or recombinations should be considered equivalents of the present invention. The above description of the disclosed aspects is provided to enable any person skilled in the art to make or use the invention. Various modifications to these aspects will be readily apparent to those skilled in the art, and the general principles defined herein can be applied to other aspects without departing from the scope of the invention. Therefore, the invention is not intended to be limited to the aspects shown herein, but rather to be carried out within the widest scope consistent with the principles and novel features disclosed herein.

[0105] The above description has been given for purposes of illustration and description. Furthermore, this description is not intended to limit the embodiments of the invention to the forms disclosed herein. Although numerous exemplary aspects and embodiments have been discussed above, those skilled in the art will recognize certain variations, modifications, alterations, additions, and sub-combinations thereof.

Claims

1. A method for constructing a spatial map of a water system, characterized in that, include: The depression-filling analysis method was used to preprocess the high-precision DEM data of the water system spatial map to obtain the final corrected DEM. Based on the final corrected DEM, the total flow distribution ratio of the water system is calculated. The catchment area threshold is determined based on the total flow distribution ratio. Based on the aforementioned catchment area threshold, the water system is classified to obtain the water system classification result; Based on the final corrected DEM and the total flow allocation ratio, watershed extraction and classification are performed to generate watershed vector boundary data, resulting in a complete spatial map of the water system.

2. The method according to claim 1, characterized in that, The high-precision DEM data of the water system spatial map was preprocessed using the depression-filling analysis method to obtain the final corrected DEM, including: The DEM grid dataset was scanned grid by grid to identify all local depressions, and the spatially connected depression grid was divided into independent depression units using a region connectivity analysis algorithm. Based on the preset depression depth and depression area thresholds, a segmented depression filling calculation formula is constructed to calculate the maximum depression depth and contiguous depression area for each independent depression unit, thereby generating depression type discrimination results. Based on the depression type identification results, combined with the terrain slope calculated grid by grid, the small uplift amount of each grid is dynamically solved to generate the dynamic uplift amount; Based on the dynamic uplift, the elevation correction is performed on the areas identified as pseudo depressions in the depression type discrimination results, while the original elevation is retained for the areas identified as natural depressions, thus obtaining a preliminary corrected DEM; By integrating remote sensing data of water bodies, vector data of water systems, and thematic data of hydraulic structures, multi-level constraint rules are constructed, constraint marker layers are generated, and the preliminary modified DEM is modified based on the constraint marker layers to generate the final modified DEM. The maximum depth of the depression is the difference between the lowest grid elevation inside the depression and the elevation of the depression's overflow outlet; the contiguous area of ​​the depression is the actual surface area corresponding to all the depression grids contained in the depression; the depth discrimination threshold is 0.5 to 2.0 m, and the area discrimination threshold is 50 to 200 m².

3. The method according to claim 2, characterized in that, The maximum depth of the depression is the difference between the lowest grid elevation inside the depression and the elevation of the depression's overflow outlet; the contiguous area of ​​the depression is the actual surface area corresponding to all the depression grids contained in the depression; the depth discrimination threshold is 0.5 to 2.0 m, and the area discrimination threshold is 50 to 200 m². The calculation formula for segmented depression filling is as follows: in, The grid elevation after the depression filling treatment; The original grid elevation; This is the set of elevations of the current grid's surrounding neighborhood. This represents the lowest elevation value within the current grid neighborhood. 'To fix a small lifting amount;' H This represents the actual maximum depth of a single depression unit; A This represents the actual contiguous area of ​​a single depression unit; Set a preset threshold for judging depression depth; Set a threshold for judging the area of ​​depressions; The dynamic lifting amount ε The calculation expression is: in, ε min Minimum lift based on the base amount; k This is the slope correction factor; S This represents the terrain slope value of the current DEM grid.

4. The method according to claim 1, characterized in that, Based on the final corrected DEM, the total flow distribution ratio of the water system is calculated, including: Based on the elevation relationship between the current cell and its adjacent cells in the final corrected DEM grid, the adjacent cells on the downhill slope are determined, and the first flow distribution ratio under the influence of the slope factor of each adjacent cell on the downhill slope is calculated. Based on the water flow vectors of adjacent cells on the uphill slope that flow into the current cell, calculate the vector synthesis direction of the water flow in the current cell, and use it as the water flow direction of the current cell. Based on the current cell's water flow direction and the inflow direction of each adjacent downhill cell, calculate the angle between the direction of each adjacent downhill cell and the current cell's water flow direction. Calculate the flow inertia weight of each adjacent cell on the downhill slope based on the directional angle, and calculate the second flow distribution ratio under the influence of the flow inertia factor of each adjacent cell on the downhill slope based on the flow inertia weight. Calculate the total flow distribution ratio of each adjacent cell on the downhill slope based on the first flow distribution ratio and the second flow distribution ratio of each adjacent cell on the downhill slope.

5. The method according to claim 4, characterized in that, The expression for calculating the first traffic allocation ratio is: In the formula, d gr ( j The slope factor of the downhill cell affects the flow of water to the first... j The flow distribution ratio between adjacent cells; tan β Slope; p Assign a weight of 1.1 to the water flow; L j As a weighting factor for the contour length of the j-th adjacent cell, when j=1, 3, 5, 7 and the elevation is lower than the current cell, L j =0.5; when j=2, 4, 6, 8 and the elevation is lower than the current cell, L j = When the elevation is not lower than the current cell, L j =0, j is the number of 8 adjacent cells, starting with the cell number 1 in the due east direction and increasing clockwise; In the calculation of the water flow vector, a coordinate system is established with the center point of the current cell as the origin, with due east as the positive x-axis direction and due north as the positive y-axis direction. The azimuth angle of the vector direction starts at due north and increases clockwise. The projection of the water flow rate along the x-axis is: Projection of water flow along the y-axis: In the formula, tan β j Adjacent cells j The slope; K j Used to determine whether the adjacent cell is an uphill cell or a downhill cell. For uphill cells, K j =1; for downhill cells, K j =0; where D u (j) represents the azimuth angle from the adjacent cell to the current cell. When j=1, D u (j)=1.5π; when j=2, D u (j) = 1.75π; when j = 3, D u (j)=0; when j=4, D u (j) = 0.25π; when j = 5, D u (j)=0.5π; when j=6, D u (j) = 0.75π; when j = 7, D u (j)=π; when j=8, D u (j) = 1.25π; The current water flow direction α in the cell is: when pw y >0 and pw x ≥0, α=arctan(| pw x | / | pw y |); when pw y >0 and pw x <0, α=2π-arctan(| pw x | / | pw y |); when pw y <0 and pw x ≥0, α=π-arctan(| pw x | / | pw y |); when pw y <0 and pw x <0, α=π+arctan(| pw x | / | pw y |); When pw y =0 and pw x >0, α=π / 2; when pw y =0 and pw x <0, α=3π / 2; The inflow direction D of adjacent cells on each downhill slope d (j) is: This j=1,D d (j) = 0.5π; When j=2, D d (j) = 0.75π; When j=3, D d (j) = π; When j=4, D d (j) = 1.25π; This j=5,D d (j) = 1.5π; When j=6, D d (j) = 1.75π; When j=7, D d (j)=0; When j=8, D d (j) = 0.25π; The expression for calculating the angle γ(j) between the direction of each adjacent cell on the downhill slope and the direction of the water flow in the current cell is as follows: γ(j)=D d (j)-a In the formula, γ(j) takes values ​​in the range of (-2π, 1.75π], and α is the direction of water flow in the current cell; Water flow inertia weights for each downhill cell W ( j The calculation expression is: W ( j )=max(cos(γ(j))*P j ,0) Among them, P j To determine whether an adjacent cell is an uphill or downhill cell, for an uphill cell, P j =0, for downhill cells, P j =1; Second traffic allocation ratio d in ( j The calculation method is as follows: Total flow distribution ratio of each downhill cell d t ( j The calculation method is as follows: 。 6. The method according to claim 1, characterized in that, Determining the catchment area threshold based on the total flow allocation ratio includes: Based on the preset catchment area threshold list and the total flow distribution ratio of each cell in the DEM, all water system cells are extracted and converted into vector line features to obtain extracted water system vector data. Acquire actual water system data and convert it into vector line features. Then perform georegistration between the actual water system vector line features and the DEM data to obtain actual water system vector data. The core density of the extracted water system vector data and the actual water system vector data at each DEM cell in the watershed is calculated respectively, so as to obtain the core density distribution of the extracted water system and the core density distribution of the actual water system under different catchment area thresholds. The core density distribution of each extracted water system and three indicators of the actual water system are calculated respectively, and the geometrically weighted average of the three indicators is used as the fitting effect index. The three indicators include: root mean square error between core density distributions, Kling-Gupta efficiency coefficient and average structural similarity index. Based on the fitting effect index, the water catchment area threshold with the best fitting effect is determined as the optimal water catchment area threshold.

7. The method according to claim 1, characterized in that, Based on the aforementioned catchment area threshold, a water system classification is performed to obtain the water system classification results, including: Obtain the spatial map data of the water system in the study area, and extract the parameters of cumulative runoff, gradient, drainage area and river length of each river segment based on the water catchment area threshold; All source rivers within the study area were identified, and the river depth of each source river was calculated based on parameters such as cumulative runoff, gradient, drainage area, and river length of each river segment, using a pre-defined improved river depth calculation formula. Based on the river depth of each source river, the Shreve method is used to add up the river depths of each confluence point in the direction from upstream to downstream to obtain the depth value of the new river segment downstream of each confluence point. This process is repeated until the depth values ​​of all river segments have been calculated. For each confluence of rivers in the water system, compare the depth values ​​of the river segments flowing into that confluence point, identify the river segment with the larger depth value among the two confluence river segments as the main stream, and the river segment with the smaller depth value as the tributary, thus determining the main stream and tributary relationship of the complete water system. The main rivers connecting the entire water system within the basin are defined as Class 0 rivers, tributaries that flow directly into Class 0 rivers are defined as Class 1 rivers, tributaries that flow directly into Class 1 rivers are defined as Class 2 rivers, and so on, until all river sections are classified, thus obtaining the water system classification result.

8. The method according to claim 1, characterized in that, Based on the final corrected DEM and the total flow allocation ratio, watershed extraction and classification are performed to generate watershed vector boundary data, including: Based on the final corrected DEM of the study area, the flow direction raster data of the total flow distribution ratio, the runoff accumulation raster data, and the location coordinate data of known hydrological stations; Based on the accumulated runoff raster data and a preset threshold, river network raster data is extracted, and based on the river network raster data and the location coordinate data of the hydrological stations, the outlet location of each hydrological station is determined. Starting from the outlet location, the grid data is traced in reverse along the water flow direction, and all upstream grid cells that eventually converge to the outlet location are marked as the same watershed, thus obtaining watershed division grid data. The watershed units in the watershed division raster data are hierarchically processed, a unique identifier is generated for each watershed unit, and the upstream and downstream hierarchical relationship encoding between watershed units is established to generate watershed hierarchical information. Extract the outer boundary of each watershed unit, convert the raster-format watershed boundary into a vector polygon, and generate watershed vector boundary data.

9. A device for constructing a spatial map of a water system, characterized in that, include: The preprocessing module is used to preprocess high-precision DEM data of water system spatial map using depression filling analysis method to obtain the final corrected DEM; The calculation module is used to calculate the total flow distribution ratio of the water system based on the final corrected DEM; The determination module is used to determine the catchment area threshold based on the total flow distribution ratio; The classification module is used to classify the water system based on the water catchment area threshold and obtain the water system classification result. The generation module is used to extract and classify the watershed based on the final corrected DEM and the total flow allocation ratio, generate watershed vector boundary data, and obtain a complete water system spatial map.

10. A computer-readable storage medium, characterized in that, The storage medium stores a computer program for performing the method described in any one of claims 1-8.

Images