Method for river network grading and automatic recognition of river topological relationship based on graph theory

Abstract

The application discloses a river network grading and river topological relation automatic identification method based on graph theory, and relates to the technical field of hydrological geographic information. The method comprises the following steps: acquiring vector river network, key points and DEM data of a target basin and cleaning; constructing an initial river network graph model based on line element connection relation; fusing DEM terrain evidence and graph theory connection characteristics, constructing and solving a global potential field equation containing terrain driving and boundary constraint, calculating the flow direction potential attribute of all network nodes to determine the flow direction relation; identifying topological abnormal structures such as strong connected components in the network based on the flow direction relation, performing breaking treatment by using direction confidence, and generating a loop-free directed network structure; performing river grade division based on topological transfer rules, and associating the key points to the river network skeleton. The application solves the problems of difficult flow direction identification in plain micro-landform areas and the problem that complex loops cannot be graded by solving and optimizing the global potential field, and realizes automatic construction of the river network topology.

Description

Technical Field

[0001] This invention belongs to the field of hydrological geographic information technology, and in particular to a method for automatic identification of river network hierarchy and river topological relationships based on graph theory. Background Technology

[0002] In the fields of computational hydrology and digital watershed construction, the accurate extraction of river network topology and the logical association of water conservancy project locations are the underlying foundation for constructing hydrodynamic models and water resource allocation systems. This technology aims to transform discrete spatial vector data into directed graph structures with clear flow direction, hierarchy, and connectivity, providing necessary boundary conditions and computational logic for flood evolution simulation, pollutant diffusion analysis, and joint watershed scheduling. It is a crucial link in realizing the mapping from physical watersheds to digital twin watersheds.

[0003] Existing river network topology construction techniques primarily rely on raster analysis of digital elevation models (DEMs) or geometric rule determination based on vector data. DEM-based methods typically employ the D8 or D-infinity algorithms to calculate confluence accumulation and extract the river network; vector data-based methods mostly utilize geometric connectivity to construct an initial network and determine river flow direction and tributary relationships by comparing elevation values ​​at river endpoints or using Strahler hierarchical rules for local recursion. Furthermore, for key locations such as reservoirs and pumping stations, existing technologies often simply project them onto the nearest river segment based on spatial proximity.

[0004] However, the aforementioned methods lack global hydraulic consistency constraints, making it difficult to generate reliable topological structures in plain micro-topographic areas or complex ring-shaped river systems. Specifically, the local elevation determination method is extremely sensitive to data noise. In areas with gentle terrain, the elevation error of the DEM is often greater than the actual natural slope. Relying solely on local endpoint elevation comparisons can easily lead to flow direction conflicts or dead-end rivers, failing to form a unified potential energy gradient across the entire network. Existing recursive algorithms require the network to be a strictly tree-like structure, making it difficult to effectively analyze complex water networks containing closed loops (strongly connected components). Once a loop structure is encountered, the hierarchical algorithm often gets stuck in an infinite loop or is forced to terminate. The fusion of point data and river network data only remains at the geometric level, lacking a mechanism to use hydraulic engineering as a hydraulic constraint to reverse the river network flow direction, making it difficult to achieve deep logical fusion of heterogeneous data. Summary of the Invention

[0005] The purpose of this invention is to provide a graph theory-based method for automatic identification of river network hierarchy and river topological relationships, in order to solve one of the problems mentioned above in the existing technology.

[0006] Technical solution: A graph theory-based method for automatic identification of river network hierarchy and river topological relationships, comprising:

[0007] Acquire vector river network data, key point data, and DEM elevation data for the target watershed, and perform spatial coordinate unification and data cleaning;

[0008] Based on the connection relationship of the centerline elements in the vector river network data, an initial river network map model containing node sets and edge sets is constructed.

[0009] Based on topographic evidence mapped from DEM elevation data, combined with the topological connectivity features of the initial river network model, the flow potential energy attributes of each node or edge in the network are calculated, and the flow direction relationships of the rivers are determined accordingly.

[0010] Based on flow relationship and graph theory connectivity analysis, we identify and process topological anomalies in the network to generate acyclic directed network structures.

[0011] Based on the directed network structure, the river level is classified according to the topology transfer rules, and key points are associated with the corresponding river network skeleton to output structured topology data.

[0012] Beneficial effects: This invention solves the problems of difficult flow direction identification and inability to classify complex loops in plain micro-topographic regions by solving the overall situation energy field and optimizing the topology, thereby realizing the automated construction of river network topology. Attached Figure Description

[0013] Figure 1 This is a flowchart illustrating the steps of the basic implementation scheme based on local rules and recursive algorithms in the embodiments of this application.

[0014] Figure 2 This is a functional module architecture diagram of the graph theory-based automatic identification system for river network hierarchies and river topology relationships in this application embodiment.

[0015] Figure 3 This is a flowchart illustrating the steps of the graph theory-based automatic identification method for river network hierarchy and river topology in this application.

[0016] Figure 4 This embodiment of the application shows the river network divided into two independent, unconnected water systems.

[0017] Figure 5 This is a DEM image for automatically identifying river sources and estuaries in this application embodiment.

[0018] Figure 6 In this embodiment of the application, the DEM is automatically identified and the river source and estuary map is manually adjusted based on human experience.

[0019] Figure 7 This is a river classification diagram in the embodiments of this application.

[0020] Figure 8 This is an automatic identification map of the spatial location mapping of key points in the embodiments of this application. Detailed Implementation

[0021] Example 1 describes an overall scheme for automatic identification of river network hierarchy and river topology based on graph theory. Addressing the problems of low automation in existing technologies for river network topology identification, difficulty in fusion of heterogeneous data, and excessive reliance on terrain data accuracy, this example provides an automated processing flow based on graph theory, compatible with both basic rule-driven methods and high-precision potential field-driven methods.

[0022] Step 101: Obtain vector river network data, key point data and DEM elevation data of the target watershed, and perform spatial coordinate unification and data cleaning.

[0023] In this embodiment, the vector river network data is typically a line feature file in Shp format, containing the geometric trajectory information of the river. Key point data includes point features with hydraulic engineering attributes such as reservoirs, pumping stations, and hydrological stations. DEM elevation data provides a digital elevation model of the watershed, with a preferred resolution of at least 30 meters. Spatial coordinate unification refers to unifying the above multi-source data into the same projected coordinate system, such as the WGS84 coordinate system, through affine transformation to ensure the accuracy of spatial overlay analysis. The data cleaning process aims to remove topological errors in the data, such as dangling points and duplicate points.

[0024] During the data cleaning process, if NoData missing values ​​exist in the DEM elevation data, interpolation is performed using a weighted average of the surrounding valid pixels to fill in the missing values. If self-intersecting line features exist in the vector river network data, the line features are broken at the intersection points. If the global situational energy field equations fail to converge due to the approximately singular coefficient matrix, the basic flow direction determination method is automatically applied.

[0025] Step 102: Based on the connection relationship of the line elements in the vector river network data, construct an initial river network map model containing a set of nodes and a set of edges.

[0026] Specifically, the initial river network graph model is an undirected graph structure G=(P,E) in graph theory. Here, the node set P consists of the endpoints and intersections of line features, and the edge set E consists of line segments connecting two nodes. The construction process involves extracting the geometric endpoints of all line features, identifying spatially overlapping endpoints as nodes of the graph, and establishing adjacency relationships between nodes. This embodiment transforms geospatial geometric data into a computer-processable mathematical graph structure.

[0027] Step 103: Extract the elevation values ​​of each node based on the DEM elevation data as topographic evidence. Combine the topological connection characteristics of the initial river network map model to calculate the flow direction potential energy attributes of each node or edge in the network. Determine the flow direction relationship of the river based on the difference in flow direction potential energy attributes between adjacent nodes.

[0028] In this embodiment, the flow direction potential energy attribute is a generalized physical quantity used to indicate the movement trend of water flow. In the basic implementation, this attribute can be the node elevation weight extracted based on the DEM; in the preferred implementation, this attribute is the physical potential energy field value obtained by solving the global Laplace equation. Determining the flow direction relationship follows the physical law that water flows downhill, that is, from nodes with high potential energy attributes to nodes with low potential energy attributes. By comparing the difference in potential energy attributes between the nodes at both ends of the connecting edge, an undirected edge can be given a clear directionality.

[0029] Step 104: Based on flow relationship and graph theory connectivity analysis, identify and process topological anomalies in the network to generate acyclic directed network structures.

[0030] Topological anomalies refer to structures that prevent the formation of strict hierarchical relationships in a network. They mainly include two types: first, pseudo-nodes with a degree of 2, which increase network redundancy; and second, closed loop structures, i.e., strongly connected components, which violate the properties of directed acyclic graphs (DAGs), causing hierarchical algorithms to fail to terminate. The processing involves deleting redundant nodes, disconnecting or flipping edges forming loops, and generating a standard DAG to ensure the unidirectional transitivity of the flow logic.

[0031] Step 105: Based on the directed network structure, perform river level classification according to the topology transfer rules, and associate key points with the corresponding river network skeleton to output structured topology data.

[0032] Topology propagation rules refer to the logic of hierarchical information propagation in a network. They can employ recursive backtracking algorithms or forward iterative algorithms based on topological sorting. Using these rules, the system can automatically identify the main stream and its tributaries at each level. The output structured topology data includes not only attributes such as river hierarchy, flow direction, and length, but also the precise locations of key points within the river network and their topological affiliations with the rivers; for example, a reservoir located upstream of a certain level tributary.

[0033] Example 2 describes the data cleaning and subgraph partitioning process before constructing the graph model. To address the potential geometric loops and multi-water system mixing issues in the original vector data, this example provides a cleaning mechanism based on geometric checks and a clustering mechanism based on disjoint-set data structures.

[0034] Step 201: Perform a geometric check on individual line features in the vector river network data to determine whether their starting coordinates and ending coordinates coincide. If there are line features whose starting coordinates and ending coordinates coincide, determine that the line feature forms a closed loop structure, and perform deletion or breakage processing on the line feature to ensure the geometric acyclicity of the input data.

[0035] In some embodiments, for identified closed loop structures, deletion is not simply done randomly, but rather a deletion priority score is calculated based on the attributes of the edges. del Edges that may be data noise are prioritized for deletion. The calculation formula is: Score del (e ij )=α1*1 / (l ij +ε)+α2*1 / (|z i -z j |+ε)+α3*1 / min(deg(p i ),deg(p j )); where Score del (e ij ) is edge e ij Deletion priority score, e ij Let be the edge connecting node i and node j; α1 is the weight coefficient for the edge length term (recommended value is 0.4); l ij For edge e ij The length of ; ε is to prevent small quantities with a denominator of 0, such as 10. -6 α2 is the weighting coefficient for the elevation difference term (recommended value is 0.4); z i z j They are nodes p respectively i p j DEM elevation values, |z i -z j | represents the DEM elevation value of node i; α3 is the weight coefficient of the node degree term (recommended value is 0.2); deg(p i ) is node p i degree (number of connected edges); deg(p) j ) is node p j The degree is denoted by ; min(·) is the minimum value function. This formula tends to delete edges that are short in length, have small elevation differences (uncertain flow direction), and low connectivity.

[0036] For the set of ring edges E cycle Each edge e in ij =(p i ,p j ), calculate the edge deletion priority score. del (e ij The formula is:

[0037] Score del (e ij )=α1·1 / (L ij +ε)+α2·1 / (|z i -z j |+ε)+α3·1 / min(deg(p i),deg(p j )); where L ij For edge e ij Length (in meters); z i z j They are nodes p respectively i p j DEM elevation values ​​(in meters); deg(p i ), deg(p j ) are nodes p i p j The degree of α1, α2, and α3 are weighting coefficients, with recommended values ​​of α1=0.4, α2=0.4, and α3=0.2; ε is used to prevent small quantities with a denominator of 0, with a recommended value of 10. -6 .

[0038] Deletion priority rule explanation: The shorter the edge length, the higher the Score. del The higher the score, the more likely it is to be deleted (shorter sides are more likely to be data noise); the smaller the elevation difference between the two ends of an edge, the higher the score. del The higher the score, the more likely it is to be deleted (low slope edges have higher flow direction uncertainty); the smaller the minimum degree of the edge endpoints, the higher the score. del The higher the degree, the more likely it is to be deleted (edges connecting low-degree nodes have less impact on network connectivity).

[0039] Select Score del The highest edge is deleted. After deletion, the cycle structure is re-checked. If the cycle still exists, the above process is repeated.

[0040] Step 202: Process the vector river network data using the disjoint-set data structure algorithm to construct an element grouping mechanism that includes a parent node dictionary and a rank dictionary; traverse the line elements in the vector river network data, and if two line elements have a spatial intersection relationship, perform a merge operation to cluster line elements with the same parent node into independent river systems; extract the node coordinates and connection relationships for each independent river system to construct the node set and edge set of the initial river network map model.

[0041] Specifically, the disjoint-set data structure (DFS) algorithm manages the element set by maintaining a parent node dictionary (parent) and a rank dictionary (rank). For the search operation `find(x)`, path compression is used. If `parent(x)` is not equal to `x`, `parent(x)` is searched recursively, and the parent nodes of all nodes on the path are directly pointed to the root node. For the merge operation `union(a,b)`, the root nodes `rootA` and `rootB` of `a` and `b` are searched respectively. If `rootA` is not equal to `rootB`, their ranks (`rank`) are compared. If `rank(rootA)` is less than `rank(rootB)`, the parent node of `rootA` is set to `rootB`; otherwise, `rootB` is pointed to `rootA`. If their ranks are equal, either one is chosen as the parent node, and its rank is incremented by 1. By traversing all line feature pairs, merge operations are performed on all line features with overlapping spatial endpoints. After traversal, the line features are grouped according to the parent node dictionary. Each group constitutes a spatially independent river system, decomposing the large-scale whole-basin computation task into multiple independent subgraph computation tasks.

[0042] Example 3, such as Figure 1 As shown, a basic implementation scheme based on local rules and recursive algorithms is described. This scheme has low computational complexity and is suitable for mountainous river networks with significant terrain slopes and no complex loops, or scenarios with limited computing resources.

[0043] Step 301: Calculate the degree of each node in the initial river network model, mark the node with a degree of 1 as an endpoint, mark the node with a degree of 2 as a pseudo node, and mark the node with a degree greater than 2 as a true node.

[0044] Wherein, the degree of the node is deg(p) i ) is defined as being related to node p i The number of directly connected edges. Endpoints usually correspond to river sources or estuaries; true nodes correspond to confluence or bifurcation points of rivers; pseudo nodes are usually segmentation points or unstructured turning points during data collection, which only play a role in transmission in the topology and have no function of branching or confluence.

[0045] Step 302: Traverse all nodes in the initial river network graph model and identify nodes with a degree of 2 as pseudo nodes; for each pseudo node, obtain its two connected adjacent nodes, delete the pseudo node and its two associated edges, and establish a new connecting edge between the two adjacent nodes. Integrate the attribute information of the original two edges into the new connecting edge to simplify the network topology.

[0046] Specifically, assume pseudo-node p k Connected to node p i and p j , side e ik The length is l ik , side ekj The length is l kj During processing, remove p from the graph's adjacency list. k Added edge e ij Its length attribute is set to l ik +l kj This process is repeated iteratively until there are no nodes with a degree of 2 in the network, which reduces the size of the graph and improves the efficiency of subsequent algorithms.

[0047] Step 303: Extract the elevation values ​​of the endpoints from the DEM elevation data and use the elevation values ​​as the initial flow direction potential energy attributes of the endpoints; determine the flow direction based on the magnitude relationship of the initial flow direction potential energy attributes, wherein the endpoint with the highest flow direction potential energy attribute is determined as the source of the main stream, and the endpoint with the lowest flow direction potential energy attribute is determined as the outlet of the basin; based on the potential energy difference between the source of the main stream and the outlet of the basin, determine the flow direction of each side along the connecting path in the initial river network map model as from high potential energy node to low potential energy node.

[0048] In some embodiments, considering potential errors in the DEM data or the influence of special terrain features, this embodiment provides a correction for manually experienced weights. The system allows the user to set the weight adjustment range, and the adjusted weight W... adj (p i ) Must satisfy: |W adj (p i )-W init (p i )|≦β*(z max -z min ); where W adj (p i ) is node p i Adjusted artificial weights; W init (p i ) is node p i Initial weights based on DEM; β is an adjustment coefficient, for example, 0.3; z max and z min These represent the maximum and minimum elevations for the entire basin. If the user specifies a river source as an endpoint, its weight is adjusted to the maximum value across the entire network; if it is specified as a river mouth, it is adjusted to the minimum value. The human-computer interaction mechanism ensures the ability to correct course when automatic identification fails.

[0049] Step 304: A recursive river network branching algorithm is adopted. Starting from the watershed outlet, the connected paths are searched and decomposed in reverse in the initial river network model until a true node or endpoint is encountered, forming a river segment path. According to the flow potential energy attribute of the endpoints from high to low, based on the principle of main stream priority, the decomposed river segment paths are assigned positive level values, and the paths that flow into the existing level river segments are marked as the next level tributaries.

[0050] In the specific implementation, define the recursive function TraceUpstream(p curr This function checks p. curr In-degree: If a true node (confluence point) is encountered, a branch is generated. The system will prioritize recursively tracing the branch connected to the high-weight endpoint (i.e., the source of the main stream) and mark it as part of the main stream. For other branches, their level is marked as the current level plus 1 when the recursion returns. This process goes upstream from the determined watershed outlet, traversing the entire tree network structure to complete the decomposition and level transfer of the river segment.

[0051] Example 4 describes a flow direction calculation method based on physical potential energy fields. Addressing the problem that traditional methods rely solely on local elevation to determine flow direction, resulting in low accuracy in plain micro-topographical regions, this example proposes an optimized implementation scheme that transforms river network topology and topographic evidence into a global Laplace linear equation system. Consistent optimization of the flow direction across the entire network is achieved through mathematical physics solutions.

[0052] Step 501: In this embodiment, based on the length of each side in the initial river network map model and the edge weight characteristics calculated from the DEM elevation data, a weighted Laplace matrix L is constructed. w Based on the boundary conditions determined by key point data and DEM elevation data, a constraint matrix R and boundary driving vector b are constructed; a structure of the form (L) is constructed. w The system of linear equations for the overall state energy field is given by +λR)Φ=b, where λ is the equilibrium coefficient and Φ is the node potential energy vector to be solved. The system of linear equations for the overall state energy field is solved by the sparse linear equations iterative solution method to obtain the overall state energy value of each node in the initial river network model, and the overall state energy value is used as the flow potential energy attribute.

[0053] In this embodiment, the global state energy field equation (L w +λR)Φ=b. Where, L w is an n x n weighted Laplace matrix, describing the connection strength and diffusion capacity between nodes in the river network diagram; λ is a balance coefficient used to adjust the weights of the internal smoothing term and the external constraint term, preferably a real number between 1.0 and 10.0; R is an n x n diagonal matrix; Φ is the n-dimensional node potential energy vector to be solved; b is an n-dimensional column vector representing the driving force generated by the terrain slope and boundary conditions. Essentially, this equation seeks a scalar field Φ such that, while satisfying the external constraints, the gradient of the internal potential energy is as consistent as possible with the terrain slope. The preferred solution method is the preprocessed conjugate gradient method (PCG) or the Cholesky decomposition method. The iteration termination condition is that the residual norm is less than a preset threshold η, for example, 10. -6 .

[0054] Define the edge weight function w(e) ij For edge e ij =(pi ,p j The weight calculation formula for )∈E is: w(e ij )=1 / (L ij ) γ Among them, L ij For edge e ij The Euclidean length of the potential function is given by γ (in meters), and γ is the length decay exponent. Recommended values ​​for γ are 1 (linear decay) or 2 (square decay). Shorter sides have greater weights, indicating a stronger influence of the shorter side on the propagation of the potential function.

[0055] For each undirected edge e uv =(u,v)∈E, calculate the evidence weight w' uv The formula is: w' uv =1 / (1+exp(-k·(|d uv |-τ d )));where, |·| is the absolute value function; exp(·) is the exponential function; τ d The recommended value range for the effective slope threshold is [0.001, 0.01]; k is the kurtosis parameter of the weighting function, and the recommended value range is [50, 200], with a default value of 100.

[0056] The Sigmoid function converts the descent evidence into weighted values ​​within the range [0,1]. The greater the absolute value of the descent is than τ, the more weighted the descent becomes. d The closer the weight is to 1, the more reliable the evidence of the edge flow.

[0057] Define the potential function value φ(p) of the boundary node. i )=W(p i ); where W(p i ) represents the elevation weight of the boundary node.

[0058] Let the boundary potential function vector φ B It is an m×1 column vector, and its i-th component is the potential function value of the i-th boundary node.

[0059] The matrix form of vector b is calculated as: b = -L IB ·φ B Among them, L IB For the weighted Laplace matrix L w The bottom left block after the block rearrangement.

[0060] The k-th component of vector b (corresponding to internal node p) m+k The element-wise calculation formula for ) is:

[0061] b k =Σ j w(e (m+k),j )·W(p j); where the summation is performed only on the internal node p m+k Adjacent boundary node p j The physical meaning of this formula is: the potential function tension exerted by the boundary node on the internal node through the adjacent edge.

[0062] Step 502: For each edge in the initial river network model, obtain the elevation values ​​of the nodes at both ends of the edge using DEM elevation data, and calculate the ratio of the elevation difference of the edge to the Euclidean distance to obtain slope evidence; based on the slope evidence and the decay function of the Euclidean distance, calculate the connection weight of the edge, where the connection weight is negatively correlated with the Euclidean distance; fill the weighted Laplace matrix L according to the connection weight. w The diagonal elements are the sum of the connection weights of all connected edges of the node, and the off-diagonal elements are the negative values ​​of the connection weights of the corresponding connected edges.

[0063] Specifically, for the edge e connecting nodes u and v uv Calculate slope evidence d uv The formula is: d uv =clip((z u -z v ) / (l uv +ε),-d max ,d max ); where d uv For edge e uv Evidence of slope descent; z u z v The DEM elevation value of the node; l uv For edge e uv The Euclidean distance length; ε is a zero-preventing decimal; d max This is the maximum gradient threshold; the clip function is used to limit the gradient to [-d]. max ,d max Within a certain range, abnormally steep data can be used to prevent computational instability.

[0064] Calculate edge weight w uv The Sigmoid function is preferred to enhance the effect of slope: w uv =1 / (1+exp(-k*(|d uv |-τ d )))*(1 / l uv γ ); where w uv For edge e uv Connection weights in the potential energy field equation; exp(·) is the exponential function; k is the kurtosis parameter; |d uv | represents the absolute value of the slope evidence; τ d The effective slope threshold; uv For edge euv The length of the slope; γ is the distance attenuation exponent, usually taken as 1 or 2. The physical meaning of this formula is: when the absolute value of the slope |d uv |Exceeds the effective threshold τ d The more edges there are and the shorter the distance between them, the greater their weight and the stronger the potential energy transfer.

[0065] Matrix L w The filling rules are as follows: For non-diagonal elements L w (u,v), if u and v are adjacent, then L w (u,v)=-w uv Otherwise, it is 0. For diagonal elements L w (u,u) is the sum of the weights of all connected edges to that node, i.e., ∑w uv (Sum of all adjacent v).

[0066] To illustrate the matrix construction process more clearly, a specific numerical calculation example is provided below: Assume a simple 3-node linear river network: P1 connects to P2, and P2 connects to P3. Node elevations: P1 is 100 meters, P2 is 90 meters, and P3 is 80 meters. Side length: l 12 For 1000 meters, l 23 The distance is 2000 meters. For simplicity, assume the edge weight is inversely proportional to the distance, w = 1 / l. Then w 12 =0.001, w 23 =0.0005. Matrix L w The elements of (3x3) are: L w (1,1)=0.001;L w (1,2)=-0.001;L w (1,3)=0L w (2,1)=-0.001;L w (2,2)=0.001+0.0005=0.0015;L w (2,3)=-0.0005L w (3,1)=0;L w (3,2)=-0.0005;L w (3,3)=0.0005.

[0067] Step 503: Based on DEM elevation data, candidate exit points are identified from the endpoints of the initial river network map model, and potential energy anchoring constraints with a target potential energy value of 0 are set for the candidate exit points; based on the type attributes of key point data, constraint points with upstream and downstream indication functions are selected, and the constraint points are mapped to the nodes of the initial river network map model, and corresponding target potential energy values ​​are set according to the point type; constraint reliability weights are assigned to each node with a set potential energy anchoring constraint or target potential energy value; based on the constraint reliability weights, the diagonal elements of the constraint matrix R are constructed, and the boundary driving vector b is assembled by combining the target potential energy value with the slope characteristics of each side in the initial river network map model calculated based on DEM elevation data.

[0068] In this embodiment, the constraint matrix R and vector b carry the fusion information of multi-source data. For a given watershed outlet node p... outlet Set the target potential energy Φ target =0, and gives extremely high reliability r outlet For example, 10 5 For the reservoir outflow node p res Based on its characteristics of controlling the upstream and releasing the downstream, relative potential energy constraints can be set.

[0069] Assembling vector b involves complex boundary condition handling. To improve solution efficiency, a matrix partitioning technique is employed in the preferred embodiment. The nodes are rearranged into a boundary node set P. B and internal node set P i At this point, the system of linear equations can be decomposed, and the driving vector b of the internal nodes... i The formula for calculating the kth component is: b k =∑w kj *Φ target (j); where b k ∑ represents the driving vector component value corresponding to the internal node k; ∑ is the summation symbol, and the summation is only performed on the boundary nodes j adjacent to the internal node k; w kj Φ represents the edge weight between internal node k and boundary node j; target (j) represents the target potential energy anchorage value of boundary node j (e.g., 0 at the estuary). Physically, this means that the fixed potential energy of the boundary node exerts a pulling force on the internal nodes through the edge weights. Block processing avoids invalid calculations of all-zero vectors, thus improving the computational speed of large-scale river networks.

[0070] Example 5 describes how to use the calculated global situation energy field to handle uncertainties and complex loops in river networks.

[0071] Step 801: For each edge in the initial river network model, calculate the difference in total potential energy values ​​between the two endpoints of the edge to obtain the potential energy difference. If the absolute value of the potential energy difference is greater than a preset effective judgment threshold, then mark the edge as a directed edge pointing from a high potential energy node to a low potential energy node according to the direction of potential energy drop. The typical range of the effective judgment threshold τ is 0.01 to 0.1, which can be adjusted according to the DEM resolution and watershed characteristics. For example, for a plain water network area with a 30-meter resolution DEM, it is recommended to take τ=0.05; for a mountain river network, τ=0.02 can be taken. Calculate the absolute value of the potential energy difference and define it as the direction confidence of the edge. Store the direction confidence as an attribute of the directed edge for subsequent topological anomaly processing.

[0072] Specifically, for edge e uv Potential energy difference ΔΦ uv =Φ u -Φ v ; where ΔΦ uv Φ represents the potential energy difference between the two ends of edge uv; u Let Φ be the total state energy of node u; v Let be the total state energy value of node v. If ΔΦ uv If the flow is greater than the positive threshold τ, the flow direction is u to v; if it is less than the negative threshold -τ, the flow direction is v to u. Direction confidence (conf) uv =|ΔΦ uv This indicator intuitively reflects the reliability of flow direction determination: the greater the potential energy difference, the clearer the flow direction and the higher the confidence level; when the potential energy difference approaches zero (such as in a still water area on a plain), the flow direction is easily disturbed and the confidence level is extremely low.

[0073] Step 901: Construct a directed graph model based on the flow relationship, and use the strongly connected component algorithm to identify the closed loop structure in the directed graph model; for each closed loop structure, traverse all the directed edges contained in it and obtain the direction confidence of each directed edge; identify the directed edge with the lowest direction confidence as the target edge to break the loop, and perform deletion operation or direction flipping operation on the target edge to break the loop until the closed loop structure is eliminated, and obtain a loop-free directed network structure.

[0074] In one aspect of this embodiment, it is assumed that an edge e is identified. AB e BC e CA The triangular closed loop has direction confidence levels of conf. AB =0.8、conf BC =0.6、conf CA =0.2. Based on the minimum cost destructive strategy, select the edge e with the lowest confidence. CAAs the target edge for breaking the loop, after the deletion operation is performed, the loop is broken, forming a one-way path A→B→C.

[0075] In this embodiment, the Tarjan algorithm or Kosaraju algorithm is preferably used to calculate the strongly connected components (SCCs) of the directed graph to identify closed loop structures. If an SCC contains more than one node, it is determined to be a loop. The strategy for handling loops is to break them at the lowest cost, that is, to modify the weakest link while maintaining the most reliable overall network structure.

[0076] Suppose a strongly connected component contains n c Each node and m c If there are 1 edge, then the number of edges r to be removed satisfies the lower bound constraint:

[0077] r≥m c -n c +1; This constraint stems from graph theory properties: n c A directed acyclic graph with n nodes contains at most n nodes. c-1 The edges form a tree structure. In practice, r is taken as the minimum value that breaks all loops, and this value is equal to the size of the minimum set of feedback arcs.

[0078] Since the minimum feedback arc set problem is NP-hard, a greedy edge-removal approximation algorithm is used for large-scale graphs:

[0079] Step 1: Calculate the deletion cost cost(e) for each edge, where cost(e) = α·ΔΦ(e) + β·L(e), ΔΦ(e) is the potential difference of the edge (negative for edges with opposing potential), L(e) is the edge length, and α and β are weight coefficients (α=1, β=0.01 recommended).

[0080] Step 2: Sort all edges in ascending order of cost (with the reverse edge at the front).

[0081] Step 3: Check the sorted edges one by one. If deleting the edge can reduce the number of cycles in the graph, then add it to the feedback arc set F.

[0082] Step 4: Repeat step 3 until there are no cycles in the diagram.

[0083] Depth-First Search (DFS) is used to detect loops. DFS is performed starting from any node, maintaining three node states: WHITE, GRAY, and BLACK. If a node in the GRAY state is visited during the traversal, a loop is detected.

[0084] After removing the selected edge set F, the Tarjan strongly connected component algorithm is run again on the processed subgraph for verification. If the size of each strongly connected component is 1 (i.e., there are no non-trivial strongly connected components), the verification passes; otherwise, additional edges need to be selected and removed.

[0085] Step 1501: Construct a disruption cost function based on directional confidence, where the disruption cost function is negatively correlated with directional confidence, i.e., the lower the directional confidence, the smaller the disruption cost; calculate the disruption cost function value of each directed edge in the closed loop structure, and determine the directed edge with the smallest disruption cost function value as the disruption target edge.

[0086] Specifically, the cost function for destruction can be defined as: cost(e uv )=1 / (conf uv +δ); where cost(e uv ) is edge e uv Removal cost in the minimum feedback arc set algorithm; conf uv For edge e uv The direction confidence is δ; δ is a small positive number to prevent the denominator from being zero. The formula logic is negatively correlated, meaning the lower the confidence, the smaller the cost. Alternatively, when searching for the minimum feedback arc set, we want to remove the edge with the smallest weight (confidence). Therefore, if the standard minimum feedback arc set algorithm (removing the set of edges with the smallest weight) is used, then the conf is directly... uv As the weight of the edge, we find the edge with the smallest weight and remove it. If expressed using a cost function, it should be expressed as prioritizing the removal of the edge with the smallest cost, where the cost should be proportional to the confidence (conf). Revised description: In the preferred embodiment, we transform the problem into finding the minimum set of feedback arcs, i.e., removing a set of edges that makes the graph acyclic and minimizes the total weight of that set of edges. Here, we directly define the edge weight as the direction confidence (conf). uv The system traverses all edges in SCC and identifies the conf. uv The smallest edge is the target for breaking the cycle.

[0087] As an alternative implementation to the aforementioned confidence-based edge destruction scheme, this embodiment also provides a greedy edge deletion strategy. In this strategy, the destruction cost function comprehensively considers the potential energy difference and edge length, defined as cost = α * ΔΦ + β * l; where cost is the cost score of an edge under the greedy edge deletion strategy; α is the weight coefficient of the potential energy difference term; ΔΦ is the driving force of the potential energy difference flow of the edge; β is the weight coefficient of the edge length term; and l is the edge length. The system prioritizes deleting the edge with the minimum cost, suitable for scenarios sensitive to geometric length. As another alternative implementation, the Minimum Spanning Tree (MST) strategy can be adopted. The SCC is treated as an undirected subgraph, with confidence... uvThe reciprocal of the edge weight is used to calculate the minimum spanning tree, retaining only the edges on the tree and reorienting them according to their potential energy. This method can simplify complex network structures into tree structures to the greatest extent.

[0088] Example 6 describes an efficient hierarchical algorithm based on a directed acyclic graph (DAG). Compared to traditional recursive algorithms, this scheme can handle large-scale river networks and is less affected by local breakpoints.

[0089] Step 1001: Perform topological sorting on the acyclic directed network structure to generate a topological sorting sequence with upstream and downstream dependencies; traverse each node according to the order of the topological sorting sequence; when traversing to a merging node with multiple merging edges, calculate the backbone determination score of each merging edge based on the difference in the overall situational energy value of the upstream nodes of each merging edge and the cumulative contribution of the upstream; identify the merging edge with the highest backbone determination score as the backbone branch, and identify the remaining merging edges as tributary branches; and based on the logical relationship between the backbone branch and the tributary branch, propagate downstream along the topological sorting sequence and calculate the river level of each river segment.

[0090] Specifically, the topological sorting preferably uses the Kahn algorithm to generate a linear sequence, ensuring that for any node v in the sequence, all its upstream nodes u are ranked before v. The backbone determination score is Score(e...). uv This is used to identify which stream is the main stream at the confluence point. The formula is: Score(e) uv )=a*(Φ u -Φ v )+b*UpSize(u); among them, Score(e uv ) is the confluence edge e uv The score for determining the main branch; 'a' is the weighting coefficient of the potential energy difference term; Φ u -Φ v denoted by , b represents the potential energy drop between upstream node u and confluence node v; b is the weighting coefficient of the upstream scale term; UpSize(u) represents the cumulative catchment area or cumulative river length upstream of node u. This formula combines the two characteristics of physical potential energy drop (rapid water flow) and catchment volume (large flow rate), which is more in line with hydrological principles than a single elevation determination.

[0091] Step 1101: If all branches flowing into the confluence node have the same river class, then set the river class of the downstream section of the confluence node to the same class plus one; if there are different river classes among the branches flowing into the confluence node, then set the river class of the downstream section of the confluence node to the highest river class among all the branches flowing into the confluence node.

[0092] To illustrate this propagation rule, the following numerical examples are provided: Assume a confluence node V has two upstream nodes U1 and U2. Case A (same-level confluence): The river segment containing U1 is at level 2, and the river segment containing U2 is at level 2. According to the rule, the downstream river segment level of V = 2 + 1 = 3. Case B (different-level confluence): The river segment containing U1 is at level 3 (main stream), and the river segment containing U2 is at level 2 (tributary). According to the rule, the downstream river segment level of V = max(3,2) = 3. The forward propagation algorithm based on topological order can complete the entire network classification in one traversal, with extremely high computational efficiency and rigorous logic, avoiding the risk of stack overflow caused by excessive recursion depth.

[0093] The grade calculation uses the following condition judgment process, executed in order of priority:

[0094] Condition A: If node v is a river source node (in-degree is 0), then R(v) = 1, and the processing is complete.

[0095] Condition B: If node v has a unique incoming edge (k=1), then R(v)=R(u1), that is, inheriting the upstream node level, and the processing is complete.

[0096] Condition C: If node v has multiple incoming edges (k≥2), then perform the following sub-condition:

[0097] Condition C1: If all incoming edges originate from the same level, i.e., R(u1) = R(u2) = ... = R(u... k )=R common Then R(v) = R common+1 This means that confluence of elements of the same level leads to a level upgrade;

[0098] Condition C2: If the levels of the incoming edge source nodes are not all the same, then R(v) = max{R(u1),R(u2),...,R(u... k That is, the level of the largest incoming edge is taken as the level of the current node.

[0099] The above conditions are evaluated in the order A→B→C→C1 / C2. If a preceding condition is met, subsequent conditions are not checked. Conditions C1 and C2 are mutually exclusive, determined by whether the source levels of the incoming edges are consistent.

[0100] Example 7 describes how to accurately spatially and logically associate discretely distributed water conservancy project location data with a constructed river network topology framework. Addressing the problem in existing technologies where locations and rivers often belong to different layers and lack logical connections, this example provides an automatic linking scheme based on geometric projection and attribute fusion.

[0101] Step 1201: Traverse each point in the key point data and calculate the vertical projection distance from the point to each river segment edge in the acyclic directed network structure; identify the river segment edge with the smallest vertical projection distance as the target river segment, and compare the vertical projection distance with the preset association tolerance threshold; if the vertical projection distance is less than the association tolerance threshold, establish the topological association relationship between the point and the target river segment, and map the attribute information of the point to the attribute table of the target river segment.

[0102] Specifically, for any key point P, its coordinates are (x0, y0). For any river segment E in the river network, let its starting point be A(x1, y1) and its ending point be B(x2, y2). Calculate the perpendicular distance d from point P to line AB. perp The calculation formula is: d perp =|(x2-x1)(y1-y0)-(x1-x0)(y2-y1)| / sqrt((x2-x1) 2 +(y2-y1) 2 ); where d perp (x0, y0) represents the vertical distance from the key point to the straight line of the river section; (x1, y1) represents the plane coordinates of the key point; (x2, y2) represents the plane coordinates of the starting point A of the river section; and (x2, y2) represents the plane coordinates of the ending point B of the river section.

[0103] It should be noted that simply calculating the straight-line distance may lead to incorrect associations, such as when a point is located on the extension line of a river segment. Therefore, this embodiment further employs the line segment constraint projection method. Specifically, it calculates the projection point P of point P onto line AB. proj Use the vector dot product to determine P proj Is it inside line segment AB? If it is inside the line segment, then the distance is d. perp If the point is outside the line segment, calculate the Euclidean distance from point P to endpoints A and B respectively, and take the smaller of the two as the final point-to-line distance.

[0104] After calculating the distances from point P to all surrounding candidate river segments, the system selects the edge E with the smallest distance. min As the target river segment, a preset correlation tolerance threshold T is introduced. tol For example, a distance of 1000 meters can be set in a plain with a network of waterways, while a distance of 500 meters can be set in a mountainous area. Only when the minimum distance is less than T... tol Only when the connection is established is the association considered successful. This threshold mechanism effectively prevents isolated points far from the river channel (such as mountaintop rain gauges) from being mistakenly and forcibly linked to the river.

[0105] After successful association, the system performs attribute mapping operations. At the database level, this manifests as adding the river code and river class to the location table, and entering the target river segment E. minThe corresponding attributes are then determined; simultaneously, a new field for a list of associated facilities is added to the river segment attribute table, incorporating the unique identifier of point P. Through a bidirectional indexing mechanism, a two-way topological query function is implemented, enabling both river-to-station and station-to-river queries.

[0106] Example 8, as Figure 2 As shown in the figure, the automatic identification system for river network hierarchy and river topology based on graph theory provided in this embodiment has a logical architecture that includes five core functional modules:

[0107] The data acquisition and preprocessing module is used to acquire vector river network data, key point data, and DEM elevation data of the target watershed, and to perform spatial coordinate unification and data cleaning. This module embeds a GIS format conversion engine (such as the Geographic Data Abstraction Library, i.e., the GDAL library) and a geometric cleaning algorithm library, which is responsible for converting multi-source heterogeneous data into standardized map data.

[0108] The graph model construction module is used to construct an initial river network graph model, including a set of nodes and edges, based on the connection relationships of line features in vector river network data. This module maintains an in-memory graph data structure and supports efficient CRUD operations on nodes and edges.

[0109] The flow potential energy calculation module is used to calculate the flow potential energy attributes of each node or edge in the network based on topographic evidence mapped from DEM elevation data and combined with the topological connectivity features of the initial river network model, thereby determining the flow direction relationships of the rivers. This module integrates sparse matrix operation libraries (such as Eigen or SciPy) to construct large-scale Laplacian matrices and solve linear equations. To improve computational efficiency, this module supports GPU-accelerated computation, especially when processing ultra-large-scale river networks with tens of millions of nodes.

[0110] The topology optimization module is used to identify and handle topological anomalies in the network based on flow relationships and graph theory connectivity analysis, generating acyclic directed network structures. This module implements the Tarjan strongly connected component algorithm and the confidence-based edge weight processing logic, and is a key component to ensure the robustness of the system.

[0111] The classification and output module is used to perform river classification based on the directed network structure and according to topology transfer rules, and to associate key points with the corresponding river network skeleton, outputting structured topology data. This module is responsible for serializing the calculation results into a user-readable format, such as GeoJSON or Shapefile, and generating a visualized topology diagram.

[0112] This embodiment also provides a computer device for running the above-described system. The device includes a memory, a processor, and a computer program stored in the memory and executable on the processor.

[0113] The processor can be a central processing unit (CPU), or other general-purpose processors, digital signal processors (DSPs), application-specific integrated circuits (ASICs), field-programmable gate arrays (FPGAs), or other programmable logic devices, discrete gate or transistor logic devices, or discrete hardware components. For preferred embodiments involving large-scale matrix solving, the processor preferably includes a high-performance computing unit supporting the Single Instruction Multiple Data (SIMD) instruction set, or is equipped with a dedicated graphics processing unit (GPU) to accelerate parallel computing.

[0114] The storage may include high-speed random access memory (RAM) and may also include non-volatile memory (NVM), such as at least one disk storage device. For database storage, a relational database supporting spatial indexes (such as PostgreSQL with the PostGIS spatial database extension plugin) is preferred for storing intermediate results.

[0115] The specific data table structure is designed as follows: NodeTable, containing fields: NodeID (node ​​number), Geometry (coordinates), Elevation (elevation), Potential (calculated potential energy value), and Type (endpoint / true / false). EdgeTable, containing fields: EdgeID (edge ​​number), StartNode (start point), EndNode (end point), Length (length), SlopeEvidence (slope evidence), Confidence (direction confidence), and FlowDir (flow direction indicator). RelationTable, containing fields: SegmentID (segment number), Level (river class), MainStemID (mainstream ID), and UpstreamSegments (list of upstream segments). AssociationTable, containing fields: PointID (point number), TargetSegmentID (associated segment ID), ProjDistance (projected distance), and RelPos (relative position percentage).

[0116] With the cooperation of the aforementioned hardware facilities and data structures, the technical solution provided in this application can operate stably and efficiently, realizing automated and high-precision identification of river network topology relationships.

[0117] A graph theory-based method for automatic identification of river network hierarchy and river topological relationships, the specific implementation steps of which are as follows: Figure 3 As shown in the accompanying drawings and embodiments, this application will be further described in detail below.

[0118] Taking a certain plain watershed as the study area, the watershed has existing linear river network data in Shp format (containing 2 unconnected water systems, with a total of 63 line elements) and point key location data (containing 1 hydrological station, 5 reservoir stations, and 19 rain gauge stations, with a total of 25 point elements). The DEM data accuracy is 30m, and the target coordinate system is WGS84.

[0119] Step S1, Data Preprocessing and Normalization:

[0120] Reading the river network and point shapefiles, no loop structures were detected, proceed to the next step;

[0121] Execution line file repair: There are 2617 points in the river network file. Among them, 41 groups of points were found to be within 10-4m apart. The coordinates of these points were corrected and unified to the coordinates of the first point in the group. After correction, there are a total of 2555 points.

[0122] All data are uniformly converted to the WGS84 coordinate system through affine transformation.

[0123] Step S2: Independent river system partitioning based on the disjoint-set data structure algorithm:

[0124] The disjoint-set data structure (DFS) algorithm is used to identify the connectivity relationships of line elements in the river network shapefile, dividing the river network into two independent, disconnected water systems, such as... Figure 4 As shown.

[0125] Step S3: Based on graph theory, an undirected vertex set containing 468 points and an undirected mapping relationship set containing 468 mapping relationships were established from the vector data of River System 1; an undirected vertex set containing 2088 points and an undirected mapping relationship set containing 2088 mapping relationships were established from the vector data of River System 2.

[0126] Step S4: Calculate the degree of each vertex in the undirected graph vertex set. The number of endpoints, pseudo nodes, and true nodes in the two river systems is shown in the table below.

[0127] Degree of vertices in an undirected graph

[0128] endpoints pseudo-node True Node River System 1 7 456 5 River System 2 18 2054 16

[0129] Step S5: Obtain the elevation values ​​of each endpoint from the DEM data, assign initial weights to each endpoint according to the elevation magnitude, and adjust the weights of each endpoint based on manual experience to identify the river source and estuary. Figure 5As shown, the DEM automatically identified endpoint 14 as having the lowest elevation value at 14m, which was designated as the river mouth. Endpoint 7 was automatically identified as having the highest elevation value at 748m, which was designated as the river source. However, based on manual experience, the actual river source is endpoint 4 (actual elevation 457m). Its elevation weight was manually adjusted to 857m. After adjustment, endpoint 4 has the highest elevation weight and is therefore the river source. Figure 6 As shown.

[0130] Endpoint DEM elevation weights and human experience elevation weights

[0131] Endpoint number X coordinate Y coordinate DEM elevation weighting Human experience elevation weight 1 12897880 3552934 357 357 2 12867170 3569953 376 376 3 12842421 3554657 125 125 4 12884449 3584896 457 857 5 12835520 3539724 88 88 6 12847131 3562227 84 84 7 12893774 3535706 748 748 8 12877157 3531540 117 117 9 12862883 3572508 160 160 10 12901051 3554396 249 249 11 12877803 3572897 558 558 12 12903444 3541745 567 567 13 12882371 3564919 304 304 14 12834269 3519433 14 14 15 12854371 3557300 72 72 16 12839290 3553589 248 248 17 12865115 3534117 230 230 18 12894588 3581168 436 436

[0132] Step S6: Using a recursive river network branching algorithm, reversely find and decompose all connected paths of the entire river network. Combined with an elevation-driven river network classification algorithm, forward classify the river levels to obtain the river classification results for river system 2, such as... Figure 7 As shown. Among them, there is 1 main stream (level 1 river), 12 level 2 rivers, 3 level 3 rivers, and 1 level 4 river.

[0133] Step S7: Identify and delete the true nodes in river system 2, reconstruct the node connection relationship, and obtain a simplified river network topology containing 18 endpoints, 16 true nodes, and 33 edges.

[0134] Step S8: Read the key point line shapefile layer and extract its station name and station type attributes. For key points with station types such as reservoir, gate, pumping station, hydrology, and water level, they are generally distributed on the river channel. Calculate the vertical distance from each key point to each edge using a vertical projection algorithm, setting a tolerance of 1000m. Associate each key point with the edge that has the smallest projected distance to that point and whose minimum distance is less than the tolerance, thus completing the automatic identification of the spatial location mapping of key points. Figure 8 As shown. In this embodiment, the point-based key location data (including 1 hydrological station, 5 reservoir stations, and 19 rain gauge stations, totaling 25 point elements) shows that 1 hydrological station and 2 reservoirs are associated with their corresponding river channels within the tolerance range (e.g., Figure 8 The Datong Reservoir, Huayuan Reservoir, and Xiheyi Hydrological Station shown are not included. The other three reservoir stations are not associated with the river system because their spatial distance from the river system is greater than the tolerance. The 19 rain gauge stations are not key locations distributed on the river channel and are not associated with them either.

[0135] Key locations and their relationship with the river channel

[0136] Station Name Website Minimum distance from the river channel Is it related to the river channel? Xihe Post Station Hydrological station 465.80m yes Datong Reservoir Station 11.96m yes garden Reservoir Station 401.72m yes Red Stone Bridge Reservoir Station 2724.26m no Yellow River Factory Reservoir Station 1981.49m no Lok Ma Bridge Reservoir Station 1904.62m no

[0137] Step S9: Result Output and Visualization. Output the river network topology table, edge topology table, and point topology attribute table in Excel and shapefile formats, with 100% data integrity. Generate a visualized network diagram, displaying node size by degree and coloring edges by level. River levels are categorized as follows: Figure 7 As shown.

[0138] To address the issue of flow direction conflicts caused by data noise interference in plain micro-topographic areas, the local elevation determination method is susceptible to local elevation noise. This solution employs a method of constructing and solving a system of linear equations for the global potential energy field. By transforming topographic slope evidence, distance attenuation characteristics, and boundary conditions into a weighted Laplace matrix and driving vector in the mathematical equations, a leap from local heuristic judgment to global optimal solution calculation is achieved. This ensures the physical consistency of the flow direction potential energy distribution across the entire network and effectively overcomes the interference of local DEM errors on flow direction identification.

[0139] To address the issue of computational interruption caused by the inability of existing recursive algorithms to resolve closed loops (strongly connected components), this solution introduces a directional confidence index and a minimum-cost loop-breaking strategy. By calculating the potential energy difference to determine the reliability of the flow direction of each edge, the solution intelligently identifies the weakest and least reliable edge in the loop structure and disconnects or flips it, successfully transforming the complex network of water systems into a standard directed acyclic structure, thus ensuring the stable operation of subsequent topology hierarchical algorithms.

[0140] To address the issue that the integration of key point data with river networks is limited to the geometric projection level and lacks logical constraints, this invention directly incorporates water conservancy project points such as reservoirs and pumping stations as mathematical soft constraints (potential energy anchor points) into the potential energy field equation. This allows the hydraulic characteristics of the points to reversely correct the flow direction calculation results of the river network, achieving deep integration of heterogeneous data at the physical and logical level.

Claims

1. A graph theory-based method for automatic identification of river network hierarchy and river topological relationships, characterized in that, include: Acquire vector river network data, key point data, and DEM elevation data for the target watershed, and perform spatial coordinate unification and data cleaning; Based on the connection relationship of the centerline elements in the vector river network data, an initial river network map model containing node sets and edge sets is constructed. Elevation values ​​of each node are extracted from DEM elevation data as topographic evidence. Combined with the topological connectivity features of the initial river network model, the flow potential energy attributes of each node or edge in the network are calculated, and the flow direction relationship of the river is determined based on the difference in flow potential energy attributes between adjacent nodes. Based on flow relationship and graph theory connectivity analysis, we identify and process topological anomalies in the network to generate acyclic directed network structures. Based on the directed network structure, the river level is classified according to the topology transfer rules, and key points are associated with the corresponding river network skeleton to output structured topology data. Constructing an initial river network graph model containing a set of nodes and a set of edges, specifically including: The disjoint-set data structure is used to process vector river network data, and an element grouping mechanism containing a parent node dictionary and a rank dictionary is constructed. Traverse the line features in the vector river network data. If two line features have a spatial intersection relationship, perform a merge operation to cluster line features with the same parent node into independent river systems. For each independent river system, extract the node coordinates and connection relationships to construct the node set and edge set of the initial river network map model; Calculating the directional potential energy properties of each node or edge in the network involves constructing and solving the global potential energy field equations, including: Based on the edge weight characteristics calculated from the length of each edge in the initial river network map model and the DEM elevation data, a weighted Laplace matrix L is constructed. w ; Based on the boundary conditions determined by key point data and DEM elevation data, a constraint matrix R and a boundary driving vector b are constructed. Construct a structure of the form (L) w The system of linear equations for the whole-world potential energy field is given by +λR)φ=b, where λ is the equilibrium coefficient and φ is the nodal potential energy vector to be solved. The sparse linear equations iterative solution method is used to solve the linear equations of the overall situation energy field, and the overall situation energy value of each node in the initial river network model is obtained. The overall situation energy value is then used as the flow potential energy attribute.

2. The method according to claim 1, characterized in that, Calculate the directional potential energy properties of each node or edge in the network, and determine the flow direction relationships of the river accordingly. Specifically, this includes: Calculate the degree of each node in the initial river network model, mark the node with a degree of 1 as an endpoint, mark the node with a degree of 2 as a pseudo node, and mark the node with a degree greater than 2 as a true node; Extract the elevation values ​​of the endpoints from the DEM elevation data and use the elevation values ​​as the initial flow direction potential energy attribute of the endpoints; The flow direction is determined based on the magnitude of the initial flow potential energy properties, with the endpoint having the highest flow potential energy property being identified as the source of the main stream and the endpoint having the lowest flow potential energy property being identified as the outlet of the basin. Based on the potential energy difference between the source of the main stream and the outlet of the basin, the flow direction of each side is determined to be from the high potential energy node to the low potential energy node along the connecting path in the initial river network model.

3. The method according to claim 2, characterized in that, River classification is performed according to topology transfer rules, specifically including: A recursive river network branching algorithm is adopted, starting from the watershed outlet, to reversely search and decompose the connected paths in the initial river network model until a true node or endpoint is encountered, thus forming a river segment path. Based on the principle of prioritizing the main stream, the river segments are assigned positive grade values ​​according to the flow potential energy attributes of the endpoints from high to low. The paths that flow into the existing graded river segments are marked as the next level tributaries.

4. The method according to claim 1, characterized in that, Construct the weighted Laplace matrix L w Specifically, it includes: For each edge in the initial river network model, the elevation values ​​of the nodes at both ends of the edge are obtained using DEM elevation data, and the ratio of the elevation difference of the edge to the Euclidean distance is calculated to obtain slope evidence. Based on slope evidence and the decay function of Euclidean distance, the connection weight of the edge is calculated, where the connection weight is negatively correlated with the Euclidean distance; Fill the weighted Laplacian matrix L according to the connection weights. w The diagonal elements are the sum of the connection weights of all connected edges of the node, and the off-diagonal elements are the negative values ​​of the connection weights of the corresponding connected edges.

5. The method according to claim 1, characterized in that, Constructing the constraint matrix R and the boundary driving vector b specifically includes: Candidate exit points are identified from the endpoints of the initial river network map model based on DEM elevation data, and potential energy anchoring constraints with a target potential energy value of 0 are set for the candidate exit points. Based on the type attributes of key point data, constraint points with upstream and downstream indication functions are selected, and the constraint points are mapped to the nodes of the initial river network model. The corresponding target potential energy values ​​are set according to the point type. Assign constraint reliability weights to each node that has a potential energy anchoring constraint or a target potential energy value set. Based on the constraint reliability weights, the diagonal elements of the constraint matrix R are constructed, and the boundary driving vector b is assembled by combining the target potential energy value with the slope characteristics of each side in the initial river network map model calculated based on DEM elevation data.

6. The method according to claim 1, characterized in that, The overall situation energy values ​​of each node in the initial river network model are obtained, including the calculation of directional confidence, specifically: For each edge in the initial river network model, calculate the difference in the total potential energy values ​​of the nodes at both ends of the edge to obtain the potential energy difference; If the absolute value of the potential energy difference is greater than the preset effective judgment threshold, then the edge is marked as a directed edge from the high potential energy node to the low potential energy node according to the direction of potential energy drop. Calculate the absolute value of the potential energy difference, define it as the direction confidence of the edge, and store the direction confidence as an attribute of the directed edge for subsequent topological anomaly processing.

7. The method according to claim 6, characterized in that, Identifying and processing topological anomalies in networks to generate acyclic directed network structures, specifically including: A directed graph model is constructed based on the flow direction relationship, and the closed loop structure existing in the directed graph model is identified by the strongly connected component algorithm. For each closed loop structure, traverse all directed edges contained within it and obtain the direction confidence of each directed edge; The directed edge with the lowest directional confidence is identified as the target edge for breaking the loop. The target edge is then deleted or its direction is reversed until the closed loop structure is eliminated, resulting in a loop-free directed network structure.

8. The method according to claim 7, characterized in that, River classification is performed according to topology transfer rules, specifically including: Perform topological sorting on acyclic directed network structures to generate topological sorting sequences with upstream and downstream dependencies; Traverse each node in the order of the topological sorting sequence. When traversing to a merging node with multiple merging edges, calculate the backbone determination score of each merging edge based on the difference in the overall situational energy value of the upstream nodes of each merging edge and the cumulative contribution of the upstream nodes. The incoming edge with the highest score in the main branch determination is identified as the main branch, and the remaining incoming edges are identified as tributary branches. Based on the logical relationship between the main branch and the tributary branches, the process is propagated downstream along the topological sorting sequence and the river level of each river segment is calculated.

Images