Method for dividing calculation units of hydrological model in plain water network area

Abstract

This invention discloses a method for partitioning computational units in a hydrological model for a plain water network area. The method comprises: acquiring and preprocessing multi-source data to form an integrated dataset; extracting the river network skeleton and constructing the river network topology based on the integrated dataset; identifying water conservancy facilities and analyzing their connectivity to generate a regional connectivity feature map; and performing hydraulically driven adaptive computational unit partitioning based on the river network topology and the regional connectivity feature map. This method can effectively improve the accuracy of hydrological simulations in plain water networks and optimize computational efficiency.

Description

Technical Field

[0001] The invention relates to hydrological model technology, in particular to a method for dividing calculation units of a hydrological model in a plain water network area. Background Art

[0002] The division of calculation units in the hydrological model, as the premise and basis of hydrological simulation, directly affects the accuracy and efficiency of the model and is a key link in the construction of the hydrological model. It is especially important for plain water network areas with dense river networks, flat terrain, and numerous water conservancy projects to scientifically divide the calculation units.

[0003] Currently, commonly used computational unit division methods include DEM-based sub-basin division, regular grid division, and unstructured grid division. Traditional DEM-based algorithms such as D8 perform poorly in plain areas with low slopes; regular grid methods are computationally efficient but have difficulty expressing complex boundaries; and unstructured grids, while flexible and adaptable to terrain, are complex to construct and computationally burdensome. Existing technologies have proposed multi-factor overlay methods that comprehensively consider topography, river networks, and land use, as well as automatic computational unit division methods based on hydrological similarity. However, the applicability of these methods in plain water network areas remains limited.

[0004] The key problems existing in the existing technology in the division of computational units in plain water networks are mainly reflected in the following aspects: (1) It is difficult to maintain the topological relationship in densely interlaced river network areas. In particular, in the process of node simplification, traditional algorithms often destroy the key connection structure of the water system, resulting in errors in the judgment of water flow paths; (2) It is difficult to effectively express the dynamic connectivity under the regulation of water conservancy projects. The connectivity relationship of the water system under different operating conditions varies significantly. The existing static division framework cannot accurately reflect the changes in connectivity caused by artificial regulation; (3) In low-slope areas, the traditional terrain-based boundary judgment method is not effective. It is impossible to use hydraulic characteristics to optimize the boundaries of computational units in a targeted manner, resulting in reduced simulation accuracy. These problems seriously restrict the application effect and accuracy of hydrological models in plain water networks. Summary of the Invention

[0005] The purpose of the invention is to provide a method for dividing calculation units of a hydrological model in a plain water network area, in order to solve at least one technical problem existing in the prior art.

[0006] Technical solution, according to one aspect of the present application, provides a method for dividing calculation units of a hydrological model of a plain water network area, comprising the following steps:

[0007] Acquire and preprocess multi-source data, including digital elevation model data, remote sensing images, river network vector data, and water conservancy project facility data, to form an integrated data set under a unified reference system;

[0008] Based on the integrated dataset, the river network skeleton is extracted and the river network topology is constructed. The river network topology includes the connectivity and topological relationships of the river network.

[0009] Based on the integrated dataset and river network topology, water conservancy facilities are identified and their connectivity is analyzed to generate a regional connectivity feature map.

[0010] Based on the river network topology and regional connectivity characteristic map, the hydraulically driven calculation unit adaptive division is performed to form the plain water network hydrological calculation unit.

[0011] Beneficial effect: the present invention can effectively improve the accuracy of hydrological simulation in plain water network areas and optimize calculation efficiency. BRIEF DESCRIPTION OF THE DRAWINGS

[0012] Figure 1 A flowchart of the steps of a method for dividing calculation units of a plain water network area hydrological model provided in an embodiment of the present application.

[0013] Figure 2 A flowchart of the steps for extracting the river network skeleton and constructing the river network topology structure provided in an embodiment of the present application.

[0014] Figure 3 A flow chart of the steps for extracting stable water bodies provided in an embodiment of the present application.

[0015] Figure 4 A flowchart of the steps for calculating water flow characteristics and hydraulic characteristics provided in an embodiment of the present application. DETAILED DESCRIPTION

[0016] like Figure 1 As shown, according to one aspect of the present application, a method for dividing calculation units of a hydrological model of a plain water network area includes the following steps:

[0017] Acquire and preprocess multi-source data to obtain an integrated data set, including digital elevation model data, remote sensing images, river network vector data, time-series water body distribution information, and water conservancy project facility data;

[0018] Based on the integrated dataset, the river network skeleton is extracted and the river network topology is constructed;

[0019] Based on the integrated dataset and river network topology, water conservancy facilities are identified and their connectivity is analyzed to generate a regional connectivity feature map.

[0020] Based on the river network topology and regional connectivity characteristic map, the hydraulically driven calculation unit adaptive division is performed to form the plain water network hydrological calculation unit.

[0021] like Figure 2 As shown, according to one aspect of the present application, the steps of extracting the river network skeleton and constructing the river network topology structure include:

[0022] Analyze time-series water distribution information and extract stable water bodies;

[0023] Calculate water flow characteristics and hydraulic characteristics based on integrated data sets;

[0024] The enhanced Markov random field method is used to extract the river network skeleton by utilizing the flow characteristics and stable water bodies.

[0025] Based on the river network skeleton, identify and classify key river network nodes;

[0026] Based on key river network nodes and river network skeleton, the invariant topological feature preserving algorithm is used to construct the river network topology structure.

[0027] like Figure 3 As shown, according to one aspect of the present application, the steps of analyzing the time series water body distribution information and extracting the stable water body include:

[0028] Calculate the frequency of water bodies at each spatial location in the time series water body distribution information and generate a water body frequency grid;

[0029] Based on the water body frequency grid and time series water body distribution information, the Wavelet-SSA hybrid decomposition method is used to obtain the decomposed water body signal;

[0030] Based on the decomposed water body signal and DEM, the SSA parameters are optimized through regional segmentation strategy to obtain the optimized decomposition parameters;

[0031] By decomposing water body signals and optimizing decomposition parameters, the adaptive regional segmentation SSA method is used to reconstruct the signal and obtain the stable water body component and seasonal water body component.

[0032] The stable water body components and seasonal water body components are classified into stability grades and spatial consistency processing is performed to obtain a stability grade classification map, and the stable water bodies are determined based on it.

[0033] like Figure 4 As shown, according to one aspect of the present application, the steps of calculating the flow characteristics and hydraulic characteristics of the water body include:

[0034] Correct and fine-register remote sensing images to obtain a fine-registered image sequence;

[0035] The segmented robust optical flow algorithm is used to process the precisely registered image sequence and the temporal water distribution information to calculate the apparent displacement field.

[0036] Calculate the hydraulic gradient field based on DEM and time-series water distribution information;

[0037] Multi-scale wavelet transform is used to fuse the apparent displacement field and hydraulic gradient field to obtain multi-scale feature representation, which is then reconstructed into water flow characteristics and hydraulic characteristics through inverse wavelet transform and feature integration.

[0038] According to one aspect of the present application, the following steps are also included:

[0039] Based on the hydrological calculation units of the plain water network, a multidimensional water exchange relationship between the calculation units is constructed to form an adaptive exchange model.

[0040] The obtained computational unit set is combined with the adaptive exchange model to construct a complete hydrological model computation framework;

[0041] Calibrate and verify the parameters of the constructed hydrological model based on measured hydrological data, and evaluate the rationality of the calculation unit division;

[0042] Based on the verification results, necessary fine-tuning is performed on the computing unit boundaries to form the final computing unit division scheme.

[0043] According to one aspect of the present application, a method for dividing calculation units of a hydrological model in a plain water network area is also provided.

[0044] S1: Multi-source data acquisition and preprocessing.

[0045] S11: Acquisition and enhancement of high-precision digital elevation models (DEMs).

[0046] The original DEM data is read and processed using a combination of bidirectional filtering and multi-scale decomposition to obtain an enhanced DEM. Specifically, this involves using non-local means filtering to remove noise, then applying wavelet transform for multi-scale decomposition to enhance micro-topography features, and finally, using geomorphologically constrained reconstruction to obtain a high-precision DEM with enhanced micro-topography features.

[0047] S12: Multi-temporal remote sensing image acquisition and water body information extraction.

[0048] Satellite remote sensing images (including optical and SAR images) from multiple time phases were collected and the dynamic threshold water index (DTWI) method was used to extract temporal water distribution information. Specifically, radiometric and geometric corrections were performed on the remote sensing images from different seasons. The modified normalized difference water index (MNDWI) was then calculated. An adaptive threshold extraction algorithm was developed based on entropy feature analysis to ultimately generate spatial distribution maps of water bodies under different water conditions.

[0049] S13: River network vector data collection and structuring.

[0050] Obtain regional river network vector data and, through topological consistency checking and structural processing, generate structured river network data with hierarchical attributes. This includes checking the connectivity and topological consistency of the river network data, repairing breakpoints and hanging points, assigning Strahler-level attributes to river segments, and establishing a table of upstream and downstream connectivity relationships within the river network to form a complete river network data structure.

[0051] S14: Water conservancy project facility data acquisition and attribute assignment.

[0052] Collect data on regional water conservancy facilities (dams, culverts, pumping stations, etc.) and, through field research and historical operation records, establish a database of water conservancy project operational characteristics. Specifically, this includes collecting information on the spatial location, structural parameters, and operating rules of water conservancy projects. This database includes attributes such as scheduling rules, water flow capacity, and operating conditions. This provides data support for subsequent analysis of the impact of water conservancy projects on water flow.

[0053] S15: Data fusion and spatial consistency processing.

[0054] The multi-source heterogeneous data obtained above is coordinate-unified and spatially consistent to obtain an integrated dataset in a unified reference system. This includes: unifying the projected coordinate system, spatially registering the different data sources, addressing boundary matching issues for data at different scales, and establishing a spatial index to form a unified data access interface, providing a consistent data environment for subsequent analysis.

[0055] S2: River network skeleton extraction and topology construction based on hydraulic characteristics.

[0056] S21: Multi-temporal water body change analysis and stable water body extraction.

[0057] Analyze the temporal distribution of water bodies and extract stable and seasonal water bodies through time-frequency analysis. This involves calculating the frequency of water bodies in each pixel in the time series, applying singular spectrum analysis (SSA) to distinguish between stable and seasonal water bodies, and generating a stability classification map to provide a basis for river network skeleton extraction.

[0058] S22: Calculation of water flow characteristics and extraction of hydraulic features.

[0059] Based on multi-temporal remote sensing imagery and enhanced DEM, a method combining optical flow analysis and hydraulic gradient calculation is used to extract water flow and hydraulic characteristics. Specifically, this involves calculating the apparent displacement field of the water body using continuous temporal remote sensing imagery, combining it with the DEM to calculate the hydraulic gradient field. The resulting fusion generates a hydraulic characteristic field that describes the dynamic characteristics of the water flow, providing a basis for identifying major flow channels.

[0060] S23: River network skeleton extraction based on persistent homology.

[0061] Leveraging water flow characteristics and stable water bodies, and employing the theory of persistent homology, we extract a topologically stable river network skeleton. Specifically, we treat the hydraulic characteristic field as a high-dimensional manifold, construct its filter complex, calculate the persistence map to identify topologically significant features, and then filter through a persistence threshold to form a topologically stable river network skeleton, effectively preserving the loop and branching features in the network structure.

[0062] S24: River network skeleton extraction based on enhanced Markov random field.

[0063] Leveraging flow characteristics and stable water bodies, the Enhanced Markov Random Field (EMRF) method was employed to extract a topologically stable river network skeleton. Specifically, the method involved constructing the hydraulic characteristic field as a Markov random field with a nonlocal energy function, introducing dual prior constraints based on mobility and connectivity, designing an adaptive neighborhood structure to capture long-range dependencies, and solving the optimization problem using a graph cut algorithm. This method extracted a topologically stable river network skeleton, effectively preserving the loop and branching features within the network structure.

[0064] S25: River network node identification and classification.

[0065] Based on the river network skeleton, a hierarchical clustering method using both geometric and hydraulic characteristics was used to identify and classify key river network nodes. Specifically, this method involved extracting characteristic nodes such as river network intersections and bifurcations, calculating the hydraulic importance index of each node based on hydraulic characteristics, and applying spectral clustering to classify the nodes. Key nodes of different functional types were labeled to form a node classification table.

[0066] S26: Construction of river network topology based on invariant topological features.

[0067] Based on key river network nodes and the river network skeleton, an algorithm preserving invariant topological features is employed to construct the river network topology. Specifically, the network is represented as an attributed graph structure, topological invariants (such as the Betti number and the number of loops) are defined as topological features, these invariants are maintained during graph simplification, and key nodes are connected via minimum energy paths to form a river network model that preserves topological relationships.

[0068] S3: Identification and connectivity analysis of water conservancy project facilities.

[0069] S31: Spatial positioning of water conservancy projects and their association with river networks.

[0070] The facilities in the water conservancy project operation characteristics database are spatially associated with the river network topology to obtain river network-water conservancy project coupled data. Specifically, this involves mapping water conservancy project facilities onto the river network using spatial proximity analysis, establishing associations between water conservancy projects and river sections, and forming a coupled data structure that includes location information and connection relationships.

[0071] S32: Formal description of water conservancy project operation rules.

[0072] Based on the water conservancy project operation characteristic database, a formal expression of water conservancy project operation rules is constructed to form an operation rule knowledge base. Specifically, this includes: analyzing the water conservancy project scheduling rules, establishing a formal rule expression containing conditions, actions, and constraints, and designing an operation status prediction model based on rule reasoning to achieve a mathematical description of the water conservancy project operation behavior.

[0073] S33: Dynamic connectivity analysis of hydraulic projects based on operating conditions.

[0074] Combining river network-hydraulic project coupling data with an operational rule knowledge base, we analyze dynamic connectivity under different operating conditions and generate a connectivity state transition model. Specifically, we define a connectivity state space, establish a Petri net-based state transition model, analyze changes in river network connectivity under different operating conditions, and generate connectivity state transition rules.

[0075] S34: Connectivity uncertainty quantification and risk probability calculation.

[0076] Based on the connectivity state transition model, a Bayesian network approach was used to quantify connectivity uncertainty and obtain a connectivity probability distribution. Specifically, this approach involved identifying key uncertainties affecting connectivity, establishing a Bayesian network model to describe the conditional dependencies between these factors, and generating connectivity probability distributions under different conditions through Monte Carlo simulation, ultimately forming a risk probability map.

[0077] S35: Multi-scale connectivity index calculation.

[0078] Based on the connectivity probability distribution and river network topology, multi-scale connectivity indices are calculated to form a regional connectivity characteristic map. Specifically, this includes defining index calculation methods for local connectivity, regional connectivity, and global connectivity, taking into account the impact weights of water conservancy projects, and generating a characteristic map representing the spatial distribution of regional hydrological connectivity, providing a basis for subsequent calculation unit division.

[0079] S4: Hydraulics-driven adaptive partitioning of computational units.

[0080] S41: Preliminary identification of hydraulic response units.

[0081] Based on hydraulic characteristics and regional connectivity feature maps, a hydraulic response similarity clustering method was used to identify preliminary hydraulic response units. This involved calculating the hydraulic response characteristic vector for each grid cell within the region, applying a spectral clustering algorithm to analyze response similarity, and preliminarily dividing response units with similar hydraulic behavior, laying the foundation for finer division.

[0082] S42: Boundary optimization based on critical rheology point.

[0083] Analyze the boundary regions of the initial hydraulic response cells, identify critical flow change points based on a flow regime critical change point detection algorithm, and optimize the cell boundaries. This involves calculating the hydraulic gradient and flow regime change characteristics of the boundary region, identifying critical points where the flow regime changes significantly (such as the transition point from rapid to slow flow, diversion points, etc.), and using these points as control points for boundary optimization to generate cell boundaries that better match the hydraulic characteristics.

[0084] S43: Dynamic connectivity-driven unit adjustment.

[0085] Combined with the connectivity state transition model, the preliminary hydraulic response units are adjusted based on dynamic connectivity to form connectivity-optimized response units. This involves analyzing the connectivity changes between units under different operating conditions, identifying areas with significant connectivity changes, and adjusting unit boundaries to align with the connectivity change boundaries to ensure that the unit division reflects the dynamic connectivity characteristics.

[0086] S44: Cell boundary refinement under multi-objective constraints.

[0087] Based on the connectivity-optimized response unit, a multi-objective optimization model encompassing computational efficiency, physical consistency, and topological integrity was established to fine-tune unit boundaries. Specifically, this involved defining unit quality evaluation metrics, constructing a multi-objective optimization function, and employing a hybrid optimization algorithm based on gradient descent and simulated annealing to fine-tune unit boundaries, ensuring both physical significance and computational efficiency.

[0088] S45: Unit merging and segmentation based on similarity of hydrological processes.

[0089] After fine-tuning the unit boundaries, the hydrological process similarity of the units is analyzed and adaptively merged and split to form the final hydrological calculation unit of the plain water network. Specifically, this includes defining similarity indicators for multiple hydrological processes, including runoff generation, confluence, and river channel evolution, calculating the process similarity matrix between units, adaptively merging highly similar adjacent units based on similarity thresholds, and re-segmenting units with high internal heterogeneity, ultimately forming hydrological calculation units that balance physical significance and computational efficiency.

[0090] S5: Construction of multidimensional water exchange relationships between computing units.

[0091] S51: Analysis of hydraulic characteristics of unit interfaces.

[0092] Analyze the interfaces between hydrological calculation units in plain water networks and calculate the hydraulic characteristic parameters of these interfaces. This involves extracting the geometric characteristics of the interfaces between units, combining them with the hydraulic characteristic field to calculate the hydraulic gradient, flow direction probability, and exchange capacity of the interfaces, and establishing a parameter set describing the hydraulic characteristics of the interfaces, providing basic data for the subsequent construction of exchange relationships.

[0093] S52: Construction of switching networks based on spectral graph theory.

[0094] Based on hydraulic characteristic parameters and river network topology, spectral graph theory is applied to establish a computational unit exchange network. Specifically, computational units are represented as network nodes, and the interfaces between units are represented as edges. Edge weights are determined based on the hydraulic characteristic parameters. Laplace matrix spectral analysis is used to identify the key connections and main exchange paths in the exchange network, forming a network structure that describes the water exchange between units.

[0095] S53: Simplified representation of the multidimensional water exchange process.

[0096] Analyze and calculate the exchange process in the unit exchange network, develop a simplified expression model for multidimensional water exchange, and form an expression for the exchange relationship. Specifically, this involves decomposing the complex multidimensional exchange process into principal components, using orthogonal polynomial expansion techniques to express the spatiotemporal variation characteristics, screening the main influencing factors through sensitivity analysis, and establishing a computationally efficient simplified expression to improve computational efficiency while maintaining physical meaning.

[0097] S54: Uncertainty characterization and parameter calibration of exchange relations.

[0098] Based on the exchange relation expression, a Bayesian inference framework is constructed to quantify parameter uncertainty and form a calibrated exchange parameter set. Specifically, this involves defining the prior distribution of the exchange parameters, using observed data and watershed response characteristics to infer the posterior distribution of the parameters through the Markov Chain Monte Carlo method, and obtaining an exchange parameter set that represents uncertainty, thereby improving the robustness of the exchange calculation.

[0099] S55: Construction of adaptive exchange mechanism for dynamic weights.

[0100] Based on a calibrated exchange parameter set and connectivity state transition model, an adaptive exchange mechanism with dynamic weights is constructed to form an adaptive exchange model. Specifically, this involves designing an exchange weight function that automatically adjusts based on water conditions and water conservancy project operations, establishing multi-timescale weight update rules, and achieving dynamic adaptation of exchange relationships, thereby improving the model's adaptability to complex hydrological conditions.

[0101] According to one aspect of the present application, S21, multi-temporal water body change analysis and stable water body extraction, is specifically:

[0102] S211: Calculation of the frequency of water bodies in time series.

[0103] Read the time-series water body distribution information and construct a three-dimensional data cube (x, y, t), where x and y are spatial coordinates and t is the time dimension. The frequency of water bodies at each spatial location (x, y) is calculated to generate a water body frequency grid. Specifically, for each spatial location, the number of water body occurrences n in N time phases is counted, and the frequency f = n / N is calculated. This generates grid data with frequency values in the interval [0, 1].

[0104] S212: Water body time series signal decomposition based on Wavelet-SSA.

[0105] The time-series water distribution information and water frequency grid are read, and the time-series signal is processed using a Wavelet-SSA hybrid decomposition method to obtain a decomposed water signal. Specifically, the time-series water signal of each pixel is preprocessed using a wavelet transform to remove high-frequency noise. A trajectory matrix is then constructed and singular value decomposition is performed. Based on the energy distribution of the singular values, the signal is decomposed into trend, seasonal, and noise terms. The preprocessed and decomposed water signal components are then output.

[0106] S213: Regional adaptive SSA parameter optimization.

[0107] The decomposed water body signal and enhanced DEM are read in, and the SSA parameters are optimized through a regional segmentation strategy to obtain the optimized decomposition parameters. Specifically, the study area is divided into several sub-regions based on the DEM-derived terrain partitioning. For each sub-region, the optimal window length and SVD component selection threshold are determined through a multi-objective optimization method that minimizes reconstruction error and maximizes signal separation. Finally, a spatially distributed SSA parameter set is generated.

[0108] S214: Adaptive regional segmentation SSA signal reconstruction.

[0109] The decomposed water body signal is read and the decomposition parameters are optimized. The adaptive regional segmentation (SSA) method is used to reconstruct the signal to obtain stable water body components and seasonal water body components. Specifically, the decomposed water body signal is grouped and reconstructed using regional optimization parameters; trend terms are identified as stable water body components; seasonal terms are identified as seasonal water body components; regional boundary transition issues are addressed using wavelet boundary processing technology; and water body components with clear physical meaning are output.

[0110] S215: Stability level classification and spatial consistency optimization.

[0111] The stable and seasonal water components are read, and stability classification and spatial consistency processing are performed to produce a stability classification map. Specifically, this includes: designing a five-level classification system based on the stability index (permanent, highly stable, moderately stable, lowly stable, and temporary); applying a Markov random field model to optimize the spatial consistency of the classification; using a conditional random field algorithm to consider neighborhood information and terrain constraints; and outputting a water stability spatial distribution map to provide a basis for river network skeleton extraction.

[0112] According to one aspect of the present application, S22, water flow characteristic calculation and hydraulic characteristic extraction, is specifically:

[0113] S221: Multi-temporal remote sensing image preprocessing and registration.

[0114] Read multi-temporal remote sensing images, perform radiometric correction, geometric correction, and fine registration to obtain a precisely registered image sequence. This includes: using histogram matching to perform relative radiometric correction to eliminate atmospheric and sensor differences; using feature point matching and affine transformation for fine registration, keeping the error less than 0.5 pixels; and generating a spatially aligned multi-temporal image sequence.

[0115] S222: Segment-wise robust optical flow calculation.

[0116] The system reads the precisely registered image sequence and the time-series water distribution information, and uses the piecewise robust optical flow algorithm to calculate the apparent displacement field. Specifically, the method involves dividing the study area into multiple subregions based on the water distribution; setting different boundary constraints for each subregion; applying the improved Horn-Schunck algorithm, introducing L1-norm smoothing terms and data terms to improve robustness to discontinuities; using a multi-resolution strategy to optimize the flow field from coarse to fine scale; and outputting a two-dimensional vector field representing the apparent motion of the water.

[0117] S223: Hydraulic gradient calculation based on DEM.

[0118] The enhanced DEM and time-series water distribution information are read and the hydraulic gradient field based on the DEM is calculated to obtain the hydraulic gradient field. This includes: hydrological optimization of the DEM to remove depressions; calculation of eight-directional flow directions and cumulative flow; calculation of hydraulic gradients based on empirical hydraulic formulas; screening of effective areas based on water distribution information; and output of a gradient field representing potential flow paths.

[0119] S224: Multi-scale wavelet domain feature fusion.

[0120] The apparent displacement field and hydraulic gradient field are read and feature fused using a multi-scale wavelet transform to obtain a multi-scale feature representation. Specifically, the method involves performing a discrete wavelet transform on the two fields, decomposing them into sub-bands of different scales; fusing the corresponding sub-bands using an adaptive weighting method; determining weights based on local consistency and gradient strength; reconstructing the sub-bands; and outputting a multi-scale feature representation.

[0121] S225: Feature reconstruction and hydraulic characteristic field generation.

[0122] The multi-scale feature representation is read and, through inverse wavelet transform and feature integration, the final water flow and hydraulic characteristics are generated. This includes: reconstructing the fused feature field through weighted inverse wavelet transform; correcting the feature field by applying physical constraints (such as mass conservation and energy conservation); calculating the derived hydraulic characteristics (such as flow velocity, water depth, and Froude number); and outputting the hydraulic characteristic field that represents the dynamic characteristics of the water body.

[0123] According to one aspect of the present application, S24, river network skeleton extraction based on enhanced Markov random field, is specifically:

[0124] S241: Optimization and preprocessing of hydraulic characteristic fields.

[0125] The enhanced hydraulic feature field is obtained by reading the flow characteristics, hydraulic characteristics, and stable water bodies, optimizing and preprocessing them. This includes normalizing the hydraulic features and processing outliers; screening valid areas based on stable water body masks; applying anisotropic diffusion filtering to enhance linear features; calculating the main direction field of the hydraulic characteristics; and outputting the enhanced feature field for river network extraction.

[0126] S242: Non-local energy function construction.

[0127] The enhanced hydraulic characteristic field is read and an energy function with non-local interactions is constructed to obtain the non-local energy function. Specifically, the following steps are performed: defining the river network label set as {0,1}, representing both non-river and river networks; constructing an energy function that includes data terms, smoothing terms, and non-local terms; the data term is based on the hydraulic characteristic strength; the smoothing term uses the direction-aware Potts model; the non-local term considers interactions between distant pixels, where the interaction strength decays with distance but decays more slowly along the flow direction; and outputting an energy function that represents the probability of pixel classification.

[0128] S243: Dual prior constraint definition.

[0129] Read the stable water body and enhanced hydraulic characteristic field, define dual prior constraints based on mobility and connectivity, and obtain prior constraint conditions. Specifically, this includes defining a connectivity prior based on the stable water body to ensure the topological connectivity of the river network structure; defining a mobility prior based on the hydraulic characteristic flow direction to ensure that the river network extends along the main flow direction; constructing a mathematical expression for the constraint conditions; and outputting the prior constraints for the EMRF model.

[0130] S244: Adaptive Neighborhood Structure Design.

[0131] The enhanced hydraulic characteristic field and prior constraints are read in, and an adaptive neighborhood structure is designed to obtain an adaptive neighborhood system. Specifically, this includes: defining a direction-aware neighborhood based on local flow characteristics; designing a multi-scale neighborhood system, including local neighborhoods and non-local connections; adaptively adjusting neighborhood size and shape based on local characteristics; defining a spatially varying potential function for each pixel; and outputting the neighborhood structure for EMRF modeling.

[0132] S245: Graph cutting optimization solution and river network skeleton generation.

[0133] The non-local energy function, prior constraints, and adaptive neighborhood system are read in, and the optimization problem is solved using a graph cut algorithm to obtain the river network skeleton. Specifically, this involves constructing a weighted graph structure, where nodes represent pixels and edges represent interactions between pixels; transforming the energy minimization problem into a graph cut problem; applying the α-expansion algorithm to solve for the globally optimal label assignment; performing morphological refinement and topological restoration on the initial extraction results; and outputting the river network skeleton, which maintains the topological structure.

[0134] According to one aspect of the present application, S42, boundary optimization based on critical rheology point, is specifically:

[0135] S421: Multi-dimensional hydraulic characteristic calculation.

[0136] The initial hydraulic response unit, water flow characteristics, and hydraulic characteristics are read, and the multidimensional hydraulic characteristics of the boundary area are calculated to obtain the boundary hydraulic characteristic set. This includes: extracting the unit boundary extension area (5 pixels inside and outside each); calculating multidimensional hydraulic characteristics including water depth, flow velocity, water surface width, Froude number, hydraulic radius, etc.; normalizing the calculation results; and outputting the multidimensional hydraulic characteristic data set of the boundary area.

[0137] S422: Feature selection based on information entropy.

[0138] The boundary hydraulic feature set is read and information entropy analysis is used to select the most discriminative feature combination to obtain the optimal feature subset. This includes: calculating the information entropy of each feature and the mutual information between features; screening the feature subset based on the maximum relevance minimum redundancy (mRMR) criterion; verifying feature importance using recursive feature elimination; selecting the feature combination with the highest information gain; and outputting the optimal feature subset for flow regime change point detection.

[0139] S423: Change point detection with multi-feature fusion.

[0140] The optimal feature subset is read and multi-feature fusion technology is used to detect flow pattern change points to obtain a set of candidate flow pattern change points. This includes: constructing trajectory curves in feature space; applying geometric methods based on changes in curvature and tangent angles to detect change points; using the CUSUM algorithm to detect significant changes in statistical characteristics; combining DBSCAN-based spatial clustering to identify change areas; and outputting a set of candidate points representing significant flow pattern changes.

[0141] S424: Hydraulic significance verification and screening.

[0142] The candidate rheological point set and hydraulic characteristics are read, and hydraulic significance is verified and screened to obtain critical rheological points. This includes: verifying the physical rationality of the candidate points based on classical hydraulic theory (such as Froude number changes and hydraulic jump conditions); calculating the hydraulic importance index of each candidate point; setting an importance threshold to screen key rheological points; spatially grouping the rheological points and selecting representative points; and outputting critical rheological points with clear hydraulic significance.

[0143] S425: Boundary reconstruction based on rheological points.

[0144] Read critical rheological points and preliminary hydraulic response units, reconstruct the calculation unit boundaries, and obtain the optimized rheological point boundaries. This includes: converting critical rheological points into boundary control points; constructing smooth boundary curves using B-spline interpolation; resolving boundary conflicts and topological consistency issues; applying the Snake model to optimize boundary details; and outputting unit boundaries optimized based on hydraulic properties, providing a basis for subsequent unit adjustments.

[0145] According to one aspect of the present application, a simplified expression of S53, a multi-dimensional water exchange process, is specifically:

[0146] S531: Analysis of the physical mechanism of the exchange process.

[0147] The computational unit exchange network and hydraulic characteristic parameters are read to analyze the physical mechanisms of water exchange and derive a classification of exchange mechanisms. This includes: identifying the main exchange types (e.g., river channel exchange, overtopping exchange, groundwater exchange, etc.); analyzing the driving factors and governing equations for each type of exchange; establishing a correspondence between exchange types and physical parameters; and outputting a classification of the physical mechanisms of the exchange process.

[0148] S532: Exchange process decomposition and principal component extraction.

[0149] The exchange mechanism classification and hydraulic characteristic parameters are read, and principal component analysis is used to extract the main components of the exchange process, resulting in the exchange principal components. This includes: constructing a parameter matrix for the exchange process; applying singular value decomposition to extract the principal components; analyzing the physical meaning and contribution of each principal component; retaining the principal component that explains more than 90% of the variance; and outputting a simplified representation of the basic components.

[0150] S533: Orthogonal polynomial basis function construction.

[0151] Based on the principal components of the exchange and hydraulic characteristic parameters, an orthogonal polynomial basis function system suitable for expressing water exchange is constructed, resulting in an orthogonal basis function set. Specifically, this involves selecting an appropriate family of orthogonal polynomials (such as Legendre polynomials and Chebyshev polynomials) based on the physical characteristics of the exchange process; customizing basis functions for different types of exchange processes; ensuring the completeness and orthogonality of the basis function system; and outputting the basis function set used to express the exchange process.

[0152] S534: Polynomial coefficient optimization and expression generation.

[0153] Combining an orthogonal basis function set with commutative principal components, the polynomial coefficients are optimized to generate a simplified commutative relation expression. This includes: designing an objective function based on physical constraints; determining the polynomial coefficients using the least squares method; evaluating the accuracy of the expression through cross-validation; balancing accuracy and complexity to control the polynomial order; and outputting a physically meaningful simplified commutative relation expression.

[0154] S535: Expression efficiency verification and optimization.

[0155] Computational efficiency testing and optimization of commutative relation expressions are performed to obtain efficient commutative expressions. This includes: computational complexity analysis of different expressions; comparison of measured computation times; code optimization for frequently called parts; introduction of techniques such as lookup tables to accelerate computation; and output of optimized expressions that maintain accuracy while improving computational efficiency.

[0156] Example 1: This paper describes the process of extracting river network skeletons based on an enhanced Markov random field. This example applies the river network skeleton extraction method of the present invention to a typical subregion (approximately 800 square kilometers) of the Taihu Lake Basin plain water network, addressing the issue of preserving topological relationships when processing complex river network nodes.

[0157] First, we acquired multi-source data for the study area: high-resolution digital elevation model (DEM) data with a resolution of 5 meters; multi-temporal remote sensing imagery: 12 Sentinel-2 satellite images (10-meter resolution) from January to December 2022; existing river network vector data: 1:50,000 scale river network vectors; and water conservancy facility data: including the spatial location and attribute information of 42 sluice gates, 17 pumping stations, and 28 culverts within the study area. Through coordinate unification and spatial registration, we formed an integrated dataset in a unified reference system (CGCS2000). The DEM data was subjected to non-local means filtering for noise removal, and wavelet transforms were used for multi-scale decomposition to enhance microtopographic features, resulting in an enhanced DEM (E-DEM).

[0158] Based on multi-temporal Sentinel-2 remote sensing images, the modified normalized difference water index (MNDWI) was calculated as (G -SWIR) / (G + SWIR); where G is the reflectance of the green band and SWIR is the reflectance of the shortwave infrared band.

[0159] The water distribution of 12 images is extracted by dynamic threshold method to form the time series water distribution information W(x,y,t), where (x,y) is the spatial coordinate and t is the time series (1 to 12). The frequency of water presence is calculated for each spatial location: f(x,y) =Σ t=1 12 W(x,y,t) / 12; generates the water body frequency grid F.

[0160] The Wavelet-SSA hybrid decomposition method is used to process the time series water distribution information:

[0161] First, perform a discrete wavelet transform on W(x,y,t) to obtain the denoised water signal Wd(x,y,t). Then construct the trajectory matrix X. For each spatial location (x,y): X(x,y) = [Wd(x,y,1),...,Wd(x,y,L); Wd(x,y,2),...,Wd(x,y,L+1); ...; Wd(x,y,K),...,Wd(x,y,12)]; where L is the window length (4), and K = 12-L+1. Perform a singular value decomposition on X(x,y): X(x,y) = UΣV T ;

[0162] Based on the terrain zoning derived from E-DEM, the study area is divided into three sub-areas (flat area, slightly inclined area and transition area). For each sub-area, the optimal SSA parameters (window length and component selection threshold) are determined by minimizing the reconstruction error to form the optimized decomposition parameter set P.

[0163] Applying the corresponding optimization parameters to each subregion, the decomposition results are divided into a trend component Wt(x,y) (stable water bodies) and a seasonal component Ws(x,y) (seasonal water bodies). Wavelet boundary processing techniques are applied to address the transition between subregion boundaries and ensure the spatial continuity of the components.

[0164] Based on the stable water component Wt(x,y), the fuzzy C-means clustering algorithm is used to divide the water body into five stability levels (permanent, highly stable, moderately stable, low stable and temporary), forming a stability level classification map C(x,y).

[0165] The segmented robust optical flow algorithm is used to calculate the apparent displacement field of the water body for the 12-phase registered Sentinel-2 image. For the two adjacent images It and It+1, the optical flow field V(u,v) is calculated, where u and v are the displacement components in the x and y directions respectively:

[0166] grad I·V + dI / dt = 0;

[0167] Introduce the robustness objective function E(V) = ∫∫[ρd(grad I·V + dI / dt) + αρs(grad u) + αρs(grad v)]dxdy;

[0168] Where ρd and ρs are L1 norm loss functions, and α is the smoothing term weight (taken as 0.05).

[0169] The hydraulic gradient field G(x,y) = [dh(x,y) / dx, dh(x,y) / dy] is calculated based on E-DEM, where h(x,y) is the elevation value in E-DEM.

[0170] The apparent displacement field V and the hydraulic gradient field G are fused by multi-scale wavelet transform. First, discrete wavelet transform is applied to V and G, and decomposed into 4 scale sub-bands DWT(V) = {VA j , VD j ^h, VD j ^v, VD j ^d} (j = 1,2,3,4) DWT(G) = {GA j , GD j ^h, GD j ^v, GD j^d} (j = 1,2,3,4);

[0171] Where A denotes the approximation subband, D denotes the detail subband, and the superscripts h, v, and d denote the horizontal, vertical, and diagonal directions, respectively. ^h denotes h as a superscript, and ^ denotes a superscript. The following formulas are equivalent.

[0172] Adopt adaptive weighting method to fuse the corresponding sub-bands: FD j ^k = ω j ^k·VD j ^k + (1-ω j ^k)·GD j ^k (k = h,v,d; j = 1,2,3,4); where the fusion weight ω j ^k calculates ω based on local consistency and gradient strength j ^k = σ(VD j ^k) / [σ(VD j ^k) + σ(GD j ^k)]; σ(·) represents the local standard deviation.

[0173] The fused hydraulic characteristic field H(x,y) is reconstructed by inverse wavelet transform, which contains hydraulic information such as flow direction and flow velocity.

[0174] Based on the stable water distribution Wt(x,y) and the hydraulic characteristic field H(x,y), an enhanced Markov random field model is constructed to extract the river network skeleton. First, H(x,y) is normalized and outliers are processed to generate the enhanced hydraulic characteristic field He(x,y).

[0175] Define the river network label set S = {0,1}, representing non-river networks and river networks. Construct an energy function E(S) = Σi[Ed(si) + Σj∈Ni Es(si,sj) + Σk∈Ri\Ni Enl(si,sk)], which includes data terms, smoothing terms, and non-local terms. The data term Ed(si) = -log p(He(i)|si) represents the probability that pixel i is classified as si, based on the intensity of He(i). The smoothing term Es(si,sj) = λs·δ(si≠sj)·g(grad He(i),grad He(j)), where λs is the weight coefficient (taken as 0.8), δ(·) is the indicator function, and g(·) is the direction-aware Potts model. The non-local term Enl(si,sk) = λnl·δ(si≠sk)·wik, where λnl is the non-local weight (taken as 0.4), and wik is the non-local interaction strength. exp(-d(i,k) / σd)·exp(-θ(i,k) / σθ); d(i,k) is the Euclidean distance between pixels i and k, θ(i,k) is the flow direction angle difference between i and k, and σd and σθ are scale parameters.

[0176] Based on the stable water body Wt(x,y), the connectivity prior constraint Pc(S) = Σi,j[δ(si=1,sj=1)·c(i,j)] is defined; where c(i,j) represents the connectivity probability between i and j, based on the distribution of stable water bodies.

[0177] The flow prior constraint Pf(S) = Σi[δ(si=1)·f(i)] is defined based on the hydraulic characteristic flow direction; where f(i) represents the flow intensity of pixel i, based on the flow direction consistency in He(i).

[0178] Design an adaptive neighborhood structure N(i), including:

[0179] Local neighborhood Ni: Based on 8-neighborhood or 16-neighborhood, adaptively adjusted according to local flow direction;

[0180] Non-local connections Ri: distant pixels selected along the main flow direction, with a maximum distance of 50 pixels.

[0181] Finally, the posterior probability of the river network label field S is: P(S|He,Wt) ∝ exp(-E(S))·Pc(S)·Pf(S);

[0182] This optimization problem is solved using the α-expansion graph cut algorithm to obtain the preliminary extraction result S* of the river network skeleton. A topology repair algorithm is applied to deal with breakpoints and dangling points to ensure the connectivity of the river network, and finally the river network skeleton R is obtained.

[0183] Based on the river network skeleton R, a feature point detection algorithm based on both geometric and hydraulic characteristics was used to identify key river network nodes N = {n1, n2, ..., nm}, including intersections and bifurcations. The hydraulic importance index (HI) of each node was calculated as follows: HI(ni) = α·Dc(ni) + β·F(ni) + γ·A(ni), where Dc is the connectivity, F is the estimated mean flow, A is the estimated catchment area, and α, β, and γ are weighting coefficients (0.3, 0.4, and 0.3, respectively).

[0184] Based on the key nodes N and the river network skeleton R, the river network topology T is constructed using an invariant topology feature-preserving algorithm. The set of topological invariants is defined as: TV = {b0, b1, CL, BR}; where b0 is the number of connected components, b1 is the number of loops (the first Betti number), CL is the set of critical links, and BR is the branching ratio.

[0185] In the process of river network simplification, these topological invariants are kept unchanged, and the key nodes are connected through the minimum energy path algorithm to form a river network structure model T that maintains the topological relationship.

[0186] In this embodiment, the river network skeleton extracted using the enhanced Markov random field method has significant advantages in complex node processing compared with the traditional method: the river network connectivity accuracy is improved by 18.7%; the topology structure retention rate is improved by 22.3%; the node recognition accuracy in densely crossed areas reaches 92.5%, which is 15.6% higher than the traditional D8 algorithm.

[0187] Especially in the processing of key nodes such as river network bifurcations and intersections, this method can effectively maintain the accuracy of the water system topology relationship, laying a solid foundation for the subsequent calculation unit division.

[0188] Example 2: This describes the dynamic connectivity analysis process of water conservancy projects based on Petri nets. This example implements a Petri net-based dynamic connectivity analysis method for water conservancy projects in a plain river network area (approximately 500 square kilometers, including 28 sluices and 12 pumping stations) with typical water conservancy project control characteristics.

[0189] The water conservancy project facility data is spatially associated with the river network topology structure constructed in Example 1. The nearest distance method and buffer analysis (buffer distance is set to 30 meters) are used to establish the corresponding relationship between water conservancy projects and river sections, forming the river network-water conservancy project coupling data RH.

[0190] For each water conservancy project facility h, record its location coordinates, the river section r where it is located, and the adjacent upstream and downstream nodes nu and nd: RH(h) = {pos(h), r(h), nu(h), nd(h)};

[0191] Based on the water conservancy project operation characteristics database, a formal expression of water conservancy project operation rules is constructed. For each water conservancy project h, an operation rule set R(h) = {r1, r2, ..., rk} is defined; each rule ri consists of a condition set C, an action set A, and a constraint set CT: ri = {C, A, CT};

[0192] The condition set C includes water level conditions, time conditions, and scheduling instruction conditions: C = {WL(u) ⊙ wl_th, T∈ [t_start, t_end], CMD = cmd_type}; where WL(u) represents the upstream water level, ⊙ represents the comparison operator (>, <, =, ≠), wl_th is the water level threshold, T is the time, [t_start, t_end] is the applicable time period, CMD is the scheduling instruction, and cmd_type is the instruction type.

[0193] Action set A includes operations such as gate opening adjustment and pump station start and stop: A = {GO = go_val, PS = ps_status}; where GO represents the gate opening, go_val is the opening value (0-100%), PS represents the pump station status, and ps_status is the start and stop status (0 or 1).

[0194] The constraint set CT includes the operational constraints: CT = {GO_rate ≤ go_rate_max, WL(d) ≤ wl_max};

[0195] Where GO_rate is the gate opening change rate, go_rate_max is the maximum allowable change rate, WL(d) is the downstream water level, and wl_max is the maximum allowable water level. Through this formal expression, a complete operation rule knowledge base KB is established.

[0196] Based on the river network-water conservancy project coupling data RH and the operation rule knowledge base KB, a Petri net model is constructed to analyze dynamic connectivity.

[0197] First, the water conservancy project facilities are divided into four functional subsystems: control gates (18); diversion gates (10); drainage pumping stations (12); channel culverts (8);

[0198] A basic Petri net model PN(h) = {P, T, F, W, M0} is constructed for each subsystem; where P = {p1, p2, ..., pn} is the set of places, representing the connectivity state (such as "fully open", "half open", "closed", etc.); T = {t1, t2, ..., tm} is the set of transitions, representing the operating condition conversion conditions (such as "water level exceeds threshold", "time reaches scheduling moment", etc.); F contained in (P×T) ∪ (T×P) is the flow relation, representing the state transition path; W: F → N+ is the weight function, indicating the importance of the transition; M0: P → N is the initial identifier, indicating the initial connectivity state, represented by a token.

[0199] A hierarchical Petri net structure is used to organize the subsystem models. The top-level Petri net HP represents the overall connectivity: HP = {HSP, HT, HF, CC}; where HSP = {PN(h1), PN(h2), ..., PN(hk)} is the subnet set; HT is the global transition set; HF is the inter-layer connection relationship; CC is the color set, with different colors used to mark different types of connectivity relationships.

[0200] Through Petri net simulation analysis, a connectivity state transition model (CST) is generated. For each pair of adjacent computational units (i, j), a connectivity state set (CS(i, j) = {cs1, cs2, ..., csv}) is defined. Each connectivity state (csk) corresponds to a specific project operating condition, and transitions between these conditions are performed using a state transition function (ST): ST: CS(i, j) × CON → CS(i, j). CON is a set of conditions, including water level conditions, time conditions, and manual scheduling conditions.

[0201] Based on the connectivity state transition model (CST), a Bayesian network approach was used to quantify connectivity uncertainty. A Bayesian network (BN) = {V, E, CPT} was constructed. V = {v1, v2, ..., vl} is a node set consisting of water conservancy project status nodes, hydrological condition nodes, and connectivity status nodes. E contained in V × V is a directed edge set representing conditional dependencies. CPT is a conditional probability table that defines the probability distribution of each node under the conditions of its parent node.

[0202] For each pair of adjacent computational cells (i, j), the probability distribution of different connectivity states P(CS(i, j)) = {P(cs1), P(cs2), ..., P(csv)} is calculated. A connectivity probability distribution map is generated based on different hydrological scenarios and scheduling schemes using the Monte Carlo simulation method.

[0203] Based on the connectivity probability distribution and river network topology, multi-scale connectivity indices are calculated to form a regional connectivity characteristic map CI.

[0204] Three scales of connectivity indices are defined: local connectivity index LCI(i,j): direct connectivity between adjacent units i and j; regional connectivity index RCI(R): internal connectivity within region R; global connectivity index GCI: overall connectivity of the entire study area;

[0205] where the local connectivity index LCI(i,j) = Σk[P(csk) × W(csk)]; W(csk) is the weight of the connectivity state csk, which is determined based on the hydraulic flux capacity.

[0206] Regional connectivity index RCI(R) = [Σi,j∈R LCI(i,j)] / [n(n-1) / 2]; n is the number of computational units in region R.

[0207] The global connectivity index (GCI) is calculated as follows: λ1·DCI + λ2·FCI + λ3·TCI; DCI is the connectivity index during the flood season, FCI is the connectivity index during the dry season, TCI is the connectivity index during the transition period, and λ1, λ2, and λ3 are time period weight coefficients (0.4, 0.3, and 0.3, respectively).

[0208] This example uses a Petri net model to analyze the dynamic connectivity of water conservancy projects and form a regional connectivity feature map. Compared with traditional static connectivity analysis methods, this method can identify and quantify connectivity changes under nine typical operating conditions, achieve an 87.6% connectivity state prediction accuracy, improve adaptability to emergency scheduling scenarios by 23.4%, and provide a reliable basis for calculating the dynamic connection of units. Especially in gate-dense areas, this method can accurately reflect changes in water system connectivity under different operating conditions, significantly improving the accuracy of the hydrological model's response to artificial regulation in plain water networks.

[0209] Example 3 mainly describes the process of boundary optimization based on the detection of critical flow regime change points. This example implements a calculation unit boundary optimization method based on the detection of critical flow regime change points in a plain river network area (approximately 300 square kilometers) with gentle terrain (average slope <0.5‰).

[0210] Based on the hydraulic characteristic field H(x, y) in Example 1 and the regional connectivity characteristic map CI in Example 2, a hydraulic response similarity clustering method is used to identify preliminary hydraulic response units.

[0211] First, the hydraulic response characteristic vector HRV(x,y) = [h(x,y), v(x,y), d(x,y), f(x,y), c(x,y)] is calculated for each grid cell in the area; where h is the water level, v is the flow velocity, d is the water depth, f is the Froude number, and c is the connectivity index.

[0212] Normalize the eigenvector to obtain NHRV(x,y). Use the spectral clustering algorithm to perform response similarity analysis: construct a similarity matrix S, where the element sij = exp(-||NHRV(i)-NHRV(j)|| 2 / 2σ 2 ); calculate the Laplace matrix L = D - S, where D is a diagonal matrix and dii = Σj sij; solve the generalized eigenvalue problem Lx = λDx and obtain the first k eigenvectors; apply k-means clustering in the eigenspace to obtain the preliminary hydraulic response unit HRU0.

[0213] The boundary area of the preliminary hydraulic response unit HRU0 is analyzed, and based on the critical flow change point detection algorithm, the critical flow change point is identified and the unit boundary is optimized.

[0214] The extended area BE of the cell boundary (5 pixels inside and outside) is extracted, and the multidimensional hydraulic feature set BHF = {WD, V, WW, FR, HR} of the boundary area is calculated, where WD is the water depth, V is the flow velocity, WW is the water surface width, FR is the Froude number, and HR is the hydraulic radius.

[0215] The information entropy-based feature selection method is applied to calculate the information entropy H(X) of each feature and the mutual information I(X;Y) between features: H(X) = -Σx p(x)log p(x); I(X;Y) = Σx,yp(x,y)log[p(x,y) / (p(x)p(y))]; based on the maximum relevance minimum redundancy (mRMR) criterion, the most discriminative feature combination OFS = arg max[Σi I(Xi;C) - 1 / |S|·Σi,j∈SI(Xi;Xj)] is selected; where C is a categorical variable (indicating different flow states) and S is the selected feature set.

[0216] In this embodiment, the preferred feature subset finally selected is {FR, V, HR}.

[0217] Based on OFS, multi-feature fusion technology is used to detect flow pattern change points. The trajectory curve TR(t) = [FR(t), V(t), HR(t)] in the feature space is constructed along the boundary coordinate parameter t. The trajectory curvature κ(t) and the tangent angle change rate κ(t) are calculated as ||TR'(t)×TR''(t)|| / ||TR'(t)|| 3 τ(t) =<TR'(t),TR'(t+1)> / (||TR'(t)||·||TR'(t+1)||).

[0218] The CUSUM algorithm is used to detect significant changes in statistical characteristics: S(t) = max[0, S(t-1) + (X(t) - μ0- K)], where X(t) is the characteristic sequence, μ0 is the reference mean, and K is the sensitivity parameter (taken as 0.5).

[0219] When S(t) exceeds the threshold H (taken as 5.0), it is marked as a change point. Combined with DBSCAN-based spatial clustering (neighborhood radius ε = 3, minimum number of points MinPts = 4), the change area is identified and the candidate flow change point set CFP is formed.

[0220] Based on hydraulic principles, the physical rationality of candidate points was verified, and the hydraulic importance index (HII(p)) was calculated as follows: α1·κ(p) + α2·τ(p) + α3·ΔFR(p) + α4·CoV(p). ΔFR represents the Froude number jump, CoV represents the velocity coefficient change, and α1 to α4 represent weight coefficients (0.3, 0.25, 0.3, and 0.15, respectively). An importance threshold (HII > 0.6) was set to screen key rheological points (CPs). These were used as boundary control points, and the unit boundary was reconstructed using B-spline interpolation to obtain the optimized rheological point boundary (HRU1).

[0221] Combined with the connectivity state transition model CST in Example 2, the flow rheological point optimization boundary HRU1 is adjusted based on dynamic connectivity.

[0222] Analyze changes in inter-unit connectivity under different operating conditions and define a connectivity change index (CVI(i,j) = 1 - min_k,l[sim(CS_k(i,j), CS_l(i,j))], where sim is the connectivity similarity and CS_k represents the connectivity state under operating condition k. Identify regions with significant connectivity changes (CVI > 0.5) and adjust the unit boundaries to align with the connectivity change boundaries, forming a connectivity-optimized response unit (HRU2).

[0223] Based on the connectivity optimization response unit HRU2, a multi-objective optimization model is established to perform fine-tuning of the unit boundaries.

[0224] The element quality evaluation index Q is defined as Q = w1·QC + w2·QP + w3·QT; where QC is the computational efficiency index, QP is the physical consistency index, QT is the topological integrity index, and w1 to w3 are weight coefficients (0.3, 0.4, and 0.3, respectively).

[0225] Construct a multi-objective optimization function min[-Q(HRU)] st A(HRU) ≥ A_min C(HRU) ≤ C_max T(HRU) = T0; where A is the unit area constraint, C is the complexity constraint, and T is the topological invariant constraint.

[0226] A hybrid optimization algorithm (combining gradient descent and simulated annealing) is used to fine-tune the unit boundaries to form the boundary refinement unit HRU3.

[0227] The similarity of the hydrological processes of HRU3 is analyzed, and adaptive merging and segmentation are performed to form the final plain water network hydrological calculation unit HRU.

[0228] Define the hydrological process similarity index HPS(i,j) = β1·sim(R_i,R j ) + β2·sim(C_i,C j ) + β3·sim(F_i,F j ); where R represents the characteristics of the runoff generation process, C represents the characteristics of the confluence process, F represents the characteristics of the river channel evolution, sim is the similarity function, and β1 to β3 are weight coefficients (taken as 0.35, 0.35, and 0.3, respectively).

[0229] The process similarity matrix (SPM) between units is calculated, and highly similar adjacent units are adaptively merged based on the similarity threshold (HPS>0.75). Units with high internal heterogeneity (coefficient of variation CV>0.25) are re-segmented to ultimately form hydrological calculation units (HRUs) that balance physical significance and computational efficiency.

[0230] In this example, a boundary optimization method based on flow regime critical change point detection achieved a 32.4% improvement in boundary recognition accuracy in flat areas (slope <0.5‰) compared to traditional terrain-based unit division methods. The method also reduced the number of computational units by 17.3% while maintaining simulation accuracy. The Nash efficiency coefficient for flow simulation increased by 0.11, and the root mean square error (RMS) for water level simulation decreased by 0.08 meters. In particularly low-slope areas, this method adaptively adjusted computational unit boundaries based on flow characteristics, significantly improving the accuracy of hydrological computational unit division.

[0231] Example 4 describes the complete process of the method for dividing the hydrological model calculation unit in the plain water network area. This example is aimed at a plain water network area in the lower reaches of the Yangtze River (with an area of about 1,200 square kilometers), and applies the complete technical solution of the present invention to realize the division of the hydrological model calculation unit in the plain water network area. The characteristics of this area are: gentle terrain (average slope 0.3‰), dense river network (river network density 3.6km / km 2 ), and numerous water conservancy projects (including 63 sluices, 26 pumping stations, and 37 culverts).

[0232] S1: Multi-source data acquisition and preprocessing.

[0233] The original DEM data of the study area at a scale of 1:10000 (resolution 5 meters) were obtained and processed using the two-way filtering and multi-scale decomposition method:

[0234] Non-local means filtering was applied to remove DEM noise, with a filter window size of 7×7 pixels, a similarity window of 3×3 pixels, and a filter strength parameter of h=10. The filtered DEM was subjected to multi-scale wavelet transform decomposition using the "db4" wavelet basis function and four decomposition levels. The high-frequency components of the wavelet were adaptively enhanced with an enhancement coefficient related to the terrain relief: α(x,y) = 1 + γ·var(z), where var(z) is the local elevation variance and γ=5. The DEM was corrected using known river channel location information through geomorphic constraint reconstruction to form an enhanced DEM (E-DEM).

[0235] Twelve Sentinel-2 satellite images (10-meter resolution) and four Sentinel-1 SAR images (10-meter resolution) from January to December 2023 were collected to extract water body distribution information: the optical images were atmospherically corrected using the DOS (Dark Object Subtraction) method; the SAR images were scatter corrected and geometrically corrected, with the E-DEM as the reference image; the modified normalized difference water index (MNDWI) = (G-SWIR) / (G+SWIR) was calculated; and a water body extraction decision tree model based on entropy features was established in combination with the SAR backscatter coefficient. The OTSU adaptive threshold method was used to determine the optimal segmentation threshold and extract the spatial distribution of water bodies in different phases.

[0236] River network vector data at a scale of 1:50,000 was obtained and topological consistency checked and structured. The connectivity of the river network was checked, and 185 breakpoints and 76 hanging points were identified and repaired. River segments were assigned grade attributes based on the Strahler classification method, with a total of 6 grades. A table of upstream and downstream connectivity relationships of 458 river segments was established, recording the connectivity relationships of 458 river segments. Hydraulic characteristic parameters (river width, water depth, etc.) were assigned to each river segment.

[0237] Collect the spatial location and attribute information of water conservancy project facilities in the study area: obtain the GPS coordinates of 126 water conservancy projects through field surveys, with a positioning accuracy of 0.5 meters; collect the geometric parameters of water conservancy projects (such as gate width, height, pump station flow, etc.); organize the historical operation records from 2021 to 2023, and establish a scheduling rule database; construct an attribute table including project type, water flow capacity, and operating conditions.

[0238] The coordinates of multi-source heterogeneous data were unified and spatially consistent: the CGCS2000 coordinate system and Gauss-Krüger projection (central meridian 121°E) were uniformly adopted; different data sources were spatially aligned with 68 control points and an alignment error of <1.0 meter; the boundary matching problem of data of different scales was handled with a boundary smoothing transition algorithm; and a spatial index was established (using a quadtree structure) to form a unified data access interface.

[0239] S2: River network skeleton extraction and topology construction based on hydraulic characteristics.

[0240] Please see Example 1, which mainly includes: calculating the water body existence frequency f(x, y) of each pixel; using the Wavelet-SSA hybrid decomposition method to process the time series water body signal; optimizing the SSA parameters based on the terrain partitioning of E-DEM; using the adaptive regional segmentation SSA method to reconstruct stable water bodies and seasonal water bodies; and performing stability level classification to obtain a classification map C(x, y).

[0241] Perform fine registration processing on multi-temporal remote sensing images; use the segmented robust optical flow algorithm to calculate the apparent displacement field V; calculate the hydraulic gradient field G based on the E-DEM; perform feature fusion through multi-scale wavelet transform; and reconstruct the hydraulic characteristic field H(x,y).

[0242] Optimize the hydraulic characteristic field to obtain the enhanced hydraulic characteristic field He(x,y); construct the non-local energy function E(S); define dual prior constraints based on mobility and connectivity; design an adaptive neighborhood structure; solve the optimization problem through the graph cutting algorithm and extract the river network skeleton R.

[0243] Based on the river network skeleton R, key river network nodes are identified and classified: the Harris corner detection algorithm is used to preliminarily identify feature points and obtain the candidate node set CN; the morphological thinning algorithm is used to extract the topological skeleton Rt of the river network skeleton; the connectivity of each pixel is calculated based on the skeleton Rt, and the intersections and bifurcations with connectivity ≥ 3 are identified; the hydraulic importance index HI(ni) of the node is calculated: HI(ni) = 0.3·Dc(ni) + 0.4·F(ni) + 0.3·A(ni); where Dc is the connectivity, F is the estimated flow value, and A is the estimated catchment area value; according to the importance index and connection characteristics, the nodes are divided into four categories: main trunk intersections, main branch intersections, ordinary bifurcations, and endpoints.

[0244] Based on the key river network nodes and river network skeleton, the river network topology structure is constructed: the river network is represented as an attributed graph structure G = (V, E), where the node V represents the key point and the edge E represents the river section; the topological invariant set is defined as TV = {b0, b1,CL, BR}; where b0 is the number of connected components, b1 is the number of loops, CL is the key link set, and BR is the branching ratio; the topology preserving energy function is designed as follows: ET(G') = w1·|b0(G')-b0(G)| + w2·|b1(G')-b1(G)| + w3·dH(CL(G'),CL(G)) + w4·|BR(G')-BR(G)| where dH is the Hausdorff distance and w1 to w4 are weight coefficients (1.0, 0.8, 0.6, 0.4); minimize ET during the graph simplification process to ensure that the topological structure is maintained; calculate the minimum energy path between key nodes through the Dijkstra algorithm to form a complete river network topology T.

[0245] S3: Identification and connectivity analysis of water conservancy project facilities.

[0246] For details, see Example 2. Spatial association of water conservancy project facility data with river network topology: Use the nearest distance method to associate water conservancy project points with river network line elements; Set a buffer zone (30 meters) for spatial matching to resolve position deviation issues; Record the upstream and downstream river section IDs and associated nodes of each water conservancy project; Create a data table RH for association between water conservancy projects and river sections, containing 126 rows of records;

[0247] Analyze water conservancy project scheduling rules and extract conditions, actions, and constraints; construct a formal rule expression: ri = {C, A, CT}; define standardized operation rule templates for each type of water conservancy project; and establish a complete operation rule knowledge base KB.

[0248] The water conservancy project facilities are divided into functional subsystems; a basic Petri net model PN(h) = {P, T, F, W, M0} is constructed; a hierarchical Petri net structure is used to organize the subsystem models; a colored Petri net interaction protocol is designed; and a connectivity state transition model CST is generated.

[0249] Quantify connectivity uncertainty based on the connectivity state transition model: Identify key uncertainties affecting connectivity, including: water level observation error (standard deviation σw = 0.05m); gate control error (standard deviation σg = 0.02m); scheduling time uncertainty (standard deviation σt = 0.5h); pump station flow rate variation (coefficient of variation CV = 0.08);

[0250] A Bayesian network model BN = {V, E, CPT} was established to describe the conditional dependency between factors. A node conditional probability table was defined, such as the conditional probability of the gate state node GS: P(GS|WL,CMD) = P(GS|WL)·P(GS|CMD)·α; where α is the normalization factor. The posterior distribution of parameters was inferred using the Markov chain Monte Carlo (MCMC) method: the number of sampling times N = 10000; the burn-in period was set to 1000; the acceptance rate was controlled between 0.23 and 0.44; and the connectivity probability distribution P(CS(i,j)) under different water conditions was generated.

[0251] Based on the connectivity probability distribution and river network topology, the multi-scale connectivity index is calculated: the local connectivity index is defined as: LCI(i,j) = Σk[P(csk) × W(csk)]; where P(csk) is the connectivity state probability and W(csk) is the weight; the regional connectivity index is calculated as: RCI(R) = [Σi,j∈R LCI(i,j)] / [n(n-1) / 2]; n is the number of units in the region; the global connectivity index is calculated as: GCI = 0.4·DCI + 0.3·FCI + 0.3·TCI; DCI is the connectivity index during the flood season, FCI is the connectivity index during the dry season, and TCI is the connectivity index during the transition period; a regional connectivity feature map CI with a resolution of 50m is generated, which contains four layers (local, regional, global, and time-varying connectivity).

[0252] S4: Hydraulics-driven adaptive partitioning of computational units.

[0253] Based on the hydraulic characteristics and regional connectivity characteristic map, preliminary hydraulic response units are identified: the hydraulic response eigenvector of each grid cell is calculated: HRV(x,y) = [h(x,y), v(x,y), d(x,y), f(x,y), c(x,y)] where h is the water level, v is the flow velocity, d is the water depth, f is the Froude number, and c is the connectivity index; the eigenvector is normalized to obtain NHRV(x,y); the similarity matrix S is constructed, and the element sij = exp(-||NHRV(i)-NHRV(j)|| 2 / 2σ 2 ); calculate the normalized Laplace matrix L = I - D^(-1 / 2)·S·D^(-1 / 2); solve for the first k eigenvectors (k=18, determined by eigenvalue analysis); apply k-means clustering in the eigenspace to obtain the preliminary hydraulic response unit HRU0, which is divided into 63 units in total.

[0254] For details, please see Example 3, which calculates the multi-dimensional hydraulic characteristics of the boundary area; adopts a feature selection method based on information entropy; detects flow pattern change points through multi-feature fusion technology; performs hydraulic significance verification and screening; and reconstructs the calculation unit boundary.

[0255] Combined with the connectivity state transition model, the unit boundaries were adjusted: the connectivity changes between units under different working conditions were analyzed and the connectivity change index (CVI) was defined; the connectivity boundaries were determined using the watershed tracing algorithm; the boundary contours were smoothed using morphological opening and closing operations; the unit boundaries were adjusted where the connectivity changed significantly (CVI>0.5); and the connectivity optimization response unit HRU2 was formed, with the number of units adjusted to 57.

[0256] A multi-objective optimization model was established to fine-tune the unit boundaries: the unit quality evaluation index Q = 0.3·QC + 0.4·QP + 0.3·QT was defined; the computational efficiency index QC = f(n, A, Pe), where n is the number of units, A is the area, and Pe is the perimeter; the physical consistency index QP = g(grad H, V, CS), where grad H is the hydraulic gradient, V is the velocity field, and CS is the connectivity status; and the topological integrity index QT = h(Nc, Nr), where Nc is the number of truncated connected paths and Nr is the proportion of the retained river network skeleton.

[0257] Construct a multi-objective optimization function min[-Q(HRU)] with the following constraints: A(HRU) ≥ A_min, C(HRU) ≤ C_max, T(HRU) = T0;

[0258] A hybrid optimization algorithm was used: a simulated annealing algorithm was used to explore the global solution in the initial stage (initial temperature T0 = 100, cooling coefficient α = 0.95); a gradient descent method was used for local optimization in the refinement stage (learning rate η = 0.01, number of iterations M = 200); spline smoothing was performed on the boundary to ensure C1 continuity of the boundary curve; and a boundary refinement unit HRU3 was formed, with a total of 52 units.

[0259] Analyze the similarity of hydrological processes of the units and perform adaptive merging and segmentation: Define the hydrological process similarity index: HPS(i,j) = 0.35·sim(R_i,R j ) + 0.35·sim(C_i,C j ) + 0.3·sim(F_i,F j ); similarity of runoff process sim(R_i,R j ) uses the inverse of KL divergence; the similarity of the confluence process sim(C_i,C j ) uses the Euclidean distance of the unit line; the river evolution similarity sim(F_i,F j ) using the Manning coefficient and cross-sectional shape similarity;

[0260] The process similarity matrix (SPM) between units was calculated. A hierarchical clustering algorithm (Ward method) was used to adaptively merge highly similar adjacent units (HPS>0.75). The watershed algorithm was applied to re-segment units with high internal heterogeneity (CV>0.25). Finally, 46 hydrological calculation units (HRU) were formed, which balanced physical significance and computational efficiency.

[0261] S5: Construction of multidimensional water exchange relationships between computing units.

[0262] Analyze and calculate the interface between units and calculate hydraulic characteristic parameters:

[0263] Extract the geometric characteristics of the interface between units: interface length L(i,j); average width W(i,j); intersection angle θ(i,j); river network intersection number NC(i,j);

[0264] Combined with the hydraulic characteristic field to calculate the interface hydraulic parameters: hydraulic gradient grad H(i,j) = (H_i - H j ) / d(i,j); flow direction probability P_dir(i,j) = Φ(grad H(i,j) / σ), where Φ is the standard normal cumulative distribution function; exchange capacity EC(i,j) = K·A(i,j)·|grad H(i,j)|; K is the equivalent hydraulic conductivity, and A(i,j) is the effective area of the interface. A parameter set HPI(i,j) is established to describe the hydraulic characteristics of the interface.

[0265] Based on the hydraulic characteristic parameters and river network topology, a computational unit exchange network is established: the computational units are represented as network nodes, and the interfaces between units are represented as edges; the edge weight w(i,j) is defined as: w(i,j) = EC(i,j)·f(CI(i,j)) f(CI(i,j)) is an adjustment function based on the connectivity index; the weighted adjacency matrix W and degree matrix D are constructed; the Laplace matrix L = D – W is calculated; L is spectrally decomposed, and the eigenvalue distribution is analyzed to determine the key exchange paths; based on the second-order adjacency relationship, an extended exchange network EN including indirect exchange is established.

[0266] For details, please refer to Example 2, analyzing the physical mechanism of the exchange process; extracting the exchange principal components using principal component analysis; constructing an orthogonal polynomial basis function system; optimizing the polynomial coefficients; and performing computational efficiency testing and optimization.

[0267] First, the physical mechanisms of water exchange were analyzed in detail, identifying the main types of exchange: river exchange (65% of total exchange): direct water exchange between units connected by river networks; overtopping exchange (21% of total exchange): water exchange through overtopping areas during flood periods; and groundwater exchange (14% of total exchange): slow water exchange through aquifers.

[0268] For each exchange type, analyze its driving factors and governing equations. For example, the governing equation for river channel exchange is: Q_c(i,j) = f(ΔH, A_c, n, R, S); where ΔH is the water level difference, A_c is the cross-sectional area, n is the Manning coefficient, R is the hydraulic radius, and S is the slope. Using principal component analysis, construct the parameter matrix M for the exchange process: M = [m_11,m_12, ..., m_1p; m_21, m_22, ..., m_2p; ...; m_n1, m_n2, ..., m_np], where n is the number of exchange samples and p is the number of parameters. Perform singular value decomposition on M: M = USV^T. Analyze the singular value decay curve and retain the top k principal components (k = 4 in this example) that explain at least 90% of the variance.

[0269] According to the physical characteristics of the exchange process, the orthogonal polynomial family is selected to construct the basis function system: the Hermite polynomial is used for river channel exchange: H_n(x) = (-1)^n·e^(x 2 / 2)·d^n / dx^n(e^(-x 2 / 2)); overflow exchange uses Legendre polynomials: P_n(x) = 1 / (2^n·n!)·d^n / dx^n[(x 2-1)^n]; groundwater exchange uses the Laguerre polynomial: L_n(x) = e^x / n!·d^n / dx^n(x^n·e^(-x));

[0270] For different exchange types, construct multidimensional orthogonal basis functions: Ψ_ijk(x,y,z) = H_i(x)·P j (y)·L_k(z) where x, y, and z are normalized driving factors (such as water level difference, section parameters, and connectivity status).

[0271] Based on the principal components and orthogonal basis functions, a polynomial expression of the exchange process is constructed: Q(i,j) = Σ_l=1^La_l·Ψ_l(X) where X is the driving factor vector, a_l is the coefficient to be determined, and L is the number of polynomial terms (16-25, depending on complexity control).

[0272] The polynomial coefficients are determined using the least squares method with physical constraints: min||M·a - Y|| 2 + λ·||a|| 2 + μ·Σ_i C_i(a) where Y is the exchange flow, λ is the regularization parameter (taken as 0.01), μ is the constraint weight (taken as 1.0), and C_i is the physical constraint condition. The accuracy of the expression is evaluated by cross-validation (5-fold cross-validation, RMSE=0.037m 3 / s). Code optimization and lookup table technology are used to improve computing efficiency, build efficient exchange expressions, and increase computing speed by 42%.

[0273] Based on the exchange relation expression, the uncertainty of the parameters is quantified: the prior distribution of the exchange parameters is defined: the Manning coefficient of the river channel n ~ N(0.035, 0.005 2 ); overflow coefficient C ~ N(0.45, 0.07 2 ); Groundwater permeability coefficient K ~logN(-11.5, 0.5 2 ); using observed data and watershed response characteristics, a Bayesian inference framework was employed: P(θ|D) ∝ P(D|θ)·P(θ), where θ is the parameter vector and D is the observed data; the posterior distribution of the parameters was inferred via an adaptive Markov chain Monte Carlo method (AM-MCMC) with a chain length of 50,000 and a burn period of 10,000; adaptive step size adjustment with a target acceptance rate of 0.234; a Gelman-Rubin convergence diagnostic R < 1.1; a set of exchange parameters (EPS) representing uncertainty was generated; and 95% confidence intervals for the exchange flow were calculated to assess prediction uncertainty.

[0274] Based on the calibrated exchange parameter set and connectivity state transition model, an adaptive exchange mechanism is constructed:

[0275] Design an exchange weight function that automatically adjusts according to water conditions and project conditions: w(i,j,t) = w_base(i,j)·f_wl(WL(t))·f_cs(CS(t))·f_s(S(t)) ; w_base is the base weight; f_wl is the water level response function; f_cs is the connectivity response function; and f_s is the seasonal response function.

[0276] Establish weight update rules at multiple time scales: short-term response (hourly): based on real-time water level and gate opening; medium-term response (day): based on connectivity state transition; long-term response (monthly): based on seasonal pattern changes;

[0277] Develop a weighted smooth transition algorithm to avoid sudden changes in the exchange relationship; implement an adaptive exchange parameter table (AEPM) containing pre-calculated parameter values for different scenarios; and form a final adaptive exchange model (AEM) that can dynamically adjust the exchange relationship according to real-time conditions.

[0278] After completing the division of computing units and establishing exchange relationships, the system was verified and evaluated through two typical flood events in June 2022 and July 2023:

[0279] Evaluation of topological relationship preservation effect: the river network connectivity accuracy rate reached 94.2%; the topological structure preservation rate reached 96.8%; the key node correct identification rate reached 93.1%;

[0280] Dynamic connectivity expression of water conservancy projects: The accuracy of connectivity status prediction reached 89.4%; the recognition rate of connectivity changes for nine typical working conditions was 92.7%; and the adaptability to emergency dispatch scenarios was 26.2% higher than that of traditional methods.

[0281] The hydraulically driven cell division results in a 35.7% improvement in boundary recognition accuracy in flat areas (slope < 0.3‰). The number of computational cells was reduced from 83 in the traditional method to 46. Simulation accuracy was comparable to or better than that of the traditional method (Nash efficiency coefficient increased by 0.13).

[0282] The effect of multi-dimensional water exchange expression: the average relative error of exchange flow prediction is reduced by 46.3%; the calculation efficiency is improved by 42%; and the water balance error is controlled within ±2.5%.

[0283] Compared with traditional methods, the calculation unit division method of the hydrological model of the plain water network area in this embodiment is not only applicable to low slope areas, but also can accurately reflect the dynamic connectivity changes caused by the regulation of water conservancy projects. At the same time, it achieves a balance between computational efficiency and physical rationality, and provides more reliable technical support for the hydrological and hydrodynamic simulation of complex plain water network areas.

[0284] To address the issue of preserving topological relationships in densely intersecting river networks, this paper employs an enhanced Markov random field approach to extract the river network skeleton. Through the innovative application of constructing a nonlocal energy function, defining dual prior constraints, and designing an adaptive neighborhood structure, this method achieves precise extraction of the river network structure. In particular, by combining hydraulic characteristics with topological constraints, this approach ensures the accurate identification and retention of key nodes during river network simplification. Furthermore, by employing an algorithm that maintains invariant topological features to construct the river network topology, this approach effectively addresses the issue of preserving topological relationships in complex river network node processing and ensures accurate representation of flow paths.

[0285] To address the problem of expressing dynamic connectivity under water conservancy project regulation, this paper innovatively constructs a dynamic connectivity model based on Petri nets to analyze different operating conditions. This model divides water conservancy facilities into functional subsystems. Within the Petri net model, locations are defined to represent connectivity states, transitions to represent operating condition transition conditions, and tokens to indicate the current state. Connectivity uncertainty is then quantified using a Bayesian network approach. This method accurately describes the changes in connectivity states under different operating conditions of water conservancy projects, enabling conditional dynamic connections between computational units and effectively addressing the problem of expressing the impact of human regulation on water system connectivity.

[0286] To address the problem of inaccurate boundary determination in low-slope areas, the present invention proposes a method for identifying critical flow rheology points based on a flow regime critical change point detection algorithm. By calculating the multi-dimensional hydraulic characteristics of the boundary area of the preliminary hydraulic response unit, and applying feature selection and multi-feature fusion technology based on information entropy, the critical points where the flow regime changes significantly are accurately identified and used as control points for adjusting the unit boundary. This hydraulically driven boundary optimization method breaks through the limitations of traditional terrain-driven methods and can perform fine-grained boundary adjustments based on hydrodynamic characteristics in areas with flat terrain, significantly improving the physical rationality and simulation accuracy of the calculation unit division.

[0287] The preferred embodiments of the present invention are described in detail above. However, the present invention is not limited to the specific details in the above embodiments. Within the technical concept of the present invention, various equivalent transformations can be made to the technical solutions of the present invention, and these equivalent transformations all fall within the scope of protection of the present invention.

Claims

1. A method for dividing calculation units of a hydrological model in a plain water network area, characterized in that: The following steps are involved: Acquire and preprocess multi-source data to obtain an integrated dataset, including digital elevation model data, remote sensing images, river network vector data, time-series water body distribution information, and water conservancy project facility data; Based on the integrated dataset, the river network skeleton is extracted and the river network topology is constructed; Based on the integrated dataset and river network topology, water conservancy facilities are identified and their connectivity is analyzed to generate a regional connectivity feature map. Based on the river network topology and regional connectivity characteristic map, the hydraulic driven calculation unit is adaptively divided to form the plain water network hydrological calculation unit; The steps to extract the river network skeleton and construct the river network topology include: Analyze time-series water distribution information and extract stable water bodies; Calculate water flow characteristics and hydraulic characteristics based on integrated data sets; The enhanced Markov random field method is used to extract the river network skeleton by utilizing the flow characteristics and stable water bodies. Based on the river network skeleton, identify and classify key river network nodes; Based on the key river network nodes and river network skeleton, the invariant topological feature preservation algorithm is used to construct the river network topology structure; The steps of extracting the river network skeleton using the enhanced Markov random field method include: Optimize and pre-process the flow characteristics, hydraulic characteristics and stable water bodies to form an enhanced hydraulic characteristic field; Based on the enhanced hydraulic characteristic field, an energy function with non-local interaction is constructed to obtain a non-local energy function; Based on stabilizing water bodies and enhancing hydraulic characteristic fields, dual prior constraints of mobility and connectivity are defined to obtain prior constraint conditions. Based on the enhanced hydraulic characteristic field and prior constraints, an adaptive neighborhood structure is designed to obtain an adaptive neighborhood system. The non-local energy function, prior constraints and adaptive neighborhood system are used to solve the optimization problem through the graph cut algorithm to extract the river network skeleton.

2. The method according to claim 1, characterized in that The steps of analyzing the time series water distribution information and extracting stable water bodies include: Calculate the frequency of water bodies at each spatial location in the time series water body distribution information and generate a water body frequency grid; Based on the water body frequency grid and time series water body distribution information, the Wavelet-SSA hybrid decomposition method is used to obtain the decomposed water body signal; Based on the decomposed water body signal and digital elevation model, the SSA parameters are optimized through regional segmentation strategy to obtain the optimized decomposition parameters; By decomposing water body signals and optimizing decomposition parameters, the adaptive regional segmentation SSA method is used to reconstruct the signal and obtain the stable water body component and seasonal water body component. The stable water body components and seasonal water body components are classified into stability grades and spatial consistency processing is performed to obtain a stability grade classification map, and the stable water bodies are determined based on it.

3. The method according to claim 1, characterized in that The steps to calculate the flow characteristics and hydraulic characteristics of water bodies include: Correct and fine-register remote sensing images to obtain a fine-registered image sequence; The segmented robust optical flow algorithm is used to process the precisely registered image sequence and the temporal water distribution information to calculate the apparent displacement field. Calculate the hydraulic gradient field based on the digital elevation model and time-series water distribution information; Multi-scale wavelet transform is used to fuse the apparent displacement field and hydraulic gradient field to obtain multi-scale feature representation, which is then reconstructed into water flow characteristics and hydraulic characteristics through inverse wavelet transform and feature integration.

4. The method according to claim 1, characterized in that The steps to identify hydraulic facilities and analyze their connectivity include: Spatial association of water conservancy project facility data with river network topology to obtain river network-water conservancy project coupling data; Based on the water conservancy project operation characteristic data, a formal expression of water conservancy project operation rules is constructed to form an operation rule knowledge base; Combining river network-water conservancy project coupling data and operation rule knowledge base, a Petri net is constructed to analyze dynamic connectivity under different operating conditions and generate a connectivity state transition model; Based on the connectivity state transition model, the Bayesian network method is used to quantify connectivity uncertainty and obtain the connectivity probability distribution; Based on the connectivity probability distribution and river network topology, multi-scale connectivity indices are calculated to form a regional connectivity characteristic map.

5. The method according to claim 1, characterized in that The steps for adaptive partitioning of computational units driven by hydraulics include: Based on the hydraulic characteristics and regional connectivity characteristic map, the hydraulic response similarity clustering method was used to identify preliminary hydraulic response units; Analyze the boundary area of the preliminary hydraulic response unit, identify the critical flow change points based on the flow critical change point detection algorithm, and optimize the unit boundary; Combined with the connectivity state transition model, the preliminary hydraulic response unit is adjusted based on dynamic connectivity to form a connectivity optimized response unit; Based on the connectivity optimization response unit, a multi-objective optimization model including computational efficiency, physical consistency, and topological integrity is established to fine-tune the unit boundaries. After fine-tuning the unit boundaries, the similarity of the hydrological processes of the units is analyzed, and adaptive merging and segmentation are performed to form hydrological calculation units of the plain water network.

6. The method according to claim 1, characterized in that The steps to construct the nonlocal energy function include: An energy function consisting of data term, smoothing term and non-local term is constructed based on the enhanced hydraulic characteristic field, where the data term is calculated based on the hydraulic characteristic intensity, the smoothing term uses the direction-aware Potts model, and the non-local term considers the interaction between distant pixels. The energy function weight is set according to the physical constraints so that the nonlocal interaction strength decays with distance but the decay rate along the stream direction is less than a threshold; The energy function is normalized and converted into the expression of pixel classification probability to obtain the non-local energy function.

7. The method according to claim 2, characterized in that The steps of signal reconstruction using the adaptive segmentation SSA method include: The study area was divided into sub-areas based on the terrain zoning derived from the digital elevation model; For the decomposed water body signal of each sub-region, the optimized decomposition parameters are applied to perform singular value decomposition, and the decomposition results are separated into trend component, seasonal component and noise component; Reconstruct the trend component into a stable water component and the seasonal component into a seasonal water component; handle the transition problem of reconstruction results between adjacent sub-regions to ensure the spatial continuity of the components; The reconstruction results of all sub-regions are integrated to form a complete stable water body component and seasonal water body component.

8. The method according to claim 5, characterized in that The steps of identifying critical flow change points based on the flow critical change point detection algorithm include: Calculate the multi-dimensional hydraulic characteristics of the boundary area of the preliminary hydraulic response unit, including water depth, flow velocity, water surface width, Froude number and hydraulic radius, to form a boundary hydraulic characteristic set; The feature selection method based on information entropy is applied to the boundary hydraulic feature set to select the most discriminative feature combination and obtain the optimal feature subset. Based on the optimal feature subset, the flow pattern change points are detected by multi-feature fusion technology to obtain the candidate flow pattern point set; Verify and screen the hydraulic significance of the candidate rheological point set to determine the critical rheological points with clear hydraulic significance; Using the critical rheological point as the control point, the boundary of the computational unit is reconstructed to obtain the optimized boundary of the rheological point.

9. The method according to claim 4, characterized in that The steps for analyzing dynamic connectivity under different operating conditions based on Petri nets include: The water conservancy project facilities are divided into functional subsystems, and a basic Petri net model is constructed for each subsystem based on the river network-water conservancy project coupling data and the operation rule knowledge base; In the basic Petri net model, the definition library represents the connectivity state, the transition represents the working condition conversion condition, and the token represents the current state indication; Adopting hierarchical Petri net structure to organize each subsystem model and reduce the complexity of state space; Design an interaction protocol based on colored Petri nets, using different colors to mark different types of connectivity relationships between systems; Through model simulation and state reachability analysis, a complete connectivity state transition model is generated.

10. The method according to claim 1, characterized in that The method further includes constructing a multidimensional water exchange relationship between the hydrological calculation units of the plain water network to form an adaptive exchange model; wherein the step of constructing the multidimensional water exchange relationship between the calculation units includes: Analyze the interface between hydrological calculation units of plain water network and calculate the hydraulic characteristic parameters of the interface; Based on hydraulic characteristic parameters and river network topology, the computational unit exchange network is established using spectral graph theory. Analyze and calculate the exchange process in the unit exchange network, develop a simplified expression model for multi-dimensional water exchange, and form an exchange relationship expression; Based on the exchange relation expression, a Bayesian inference framework is constructed to quantify parameter uncertainty and form a calibrated exchange parameter set; Based on the calibrated exchange parameter set and connectivity state transition model, an adaptive exchange mechanism with dynamic weights is constructed to form an adaptive exchange model.

Images