Plain water network area hydrological model calculation unit division method

Abstract

The invention discloses a method for dividing calculation units of a hydrological model in a plain water network area. The method comprises the following steps: acquiring and preprocessing multi-source data to form an integrated data set; extracting a river network skeleton based on the integrated data set and constructing a river network topological structure; identifying hydraulic engineering facilities and analyzing the connectivity of the hydraulic engineering facilities, and generating a regional connectivity feature map; and carrying out hydraulically-driven calculation unit adaptive division based on the river network topological structure and the regional connectivity feature map. According to the method, the plain water network area hydrological simulation precision can be effectively improved, and the calculation efficiency is optimized.

Description

Technical Field

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

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

[0003] Currently, commonly used calculation unit division methods mainly include sub-basin division based on DEM, regular grid division, and unstructured grid division, etc. Traditional algorithms such as D8 based on DEM perform poorly in plain areas with low slope gradients; the regular grid method has high calculation efficiency but is difficult to represent complex boundaries; although unstructured grids can flexibly adapt to terrain, the construction process is complex and the calculation burden is large. Existing technologies have proposed a multi-factor superposition method that comprehensively considers terrain, river networks, and land use, as well as an automatic calculation unit division method based on hydrological similarity, but the applicability of these methods in plain water network areas is still limited.

[0004] The key problems existing in the division of calculation units in plain water network areas in the prior art are mainly reflected in: (1) It is difficult to maintain the topological relationship in areas with dense and intertwined river networks. Especially during the node simplification process, traditional algorithms often damage the key connection structures of the water system, resulting in incorrect judgment of water flow paths; (2) It is difficult to effectively represent the dynamic connectivity under the regulation of water conservancy projects. The water system connectivity relationships vary significantly under different operating conditions, and the existing static division framework cannot accurately reflect the connectivity changes caused by artificial regulation; (3) In areas with low slope gradients, the traditional terrain-based boundary determination method has poor effects and cannot use hydraulic characteristics to optimize the calculation unit boundaries specifically, resulting in a decrease in simulation accuracy. These problems seriously restrict the application effect and accuracy of hydrological models in plain water network areas. Summary of the Invention

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

[0006] Technical solution: 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 provided, including the following steps:

[0007] Obtain 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, extract the river network skeleton and construct the river network topological structure, which includes the connectivity and topological relationships of the river network;

[0009] Based on the integrated dataset and the river network topological structure, identify hydraulic engineering facilities and analyze their connectivity to generate a regional connectivity feature map;

[0010] Based on the river network topological structure and the regional connectivity feature map, perform hydraulics-driven adaptive partitioning of computational units to form plain water network hydrological computational units.

[0011] Advantageous effects: The present invention can effectively improve the accuracy of hydrological simulation in the plain water network area and optimize the calculation efficiency. Description of the Drawings

[0012] Figure 1 It is a step flow chart of a method for dividing hydrological model computational units in a plain water network area provided by an embodiment of the present application.

[0013] Figure 2 It is a step flow chart of extracting the river network skeleton and constructing the river network topological structure provided by an embodiment of the present application.

[0014] Figure 3 It is a step flow chart of extracting stable water bodies provided by an embodiment of the present application.

[0015] Figure 4 It is a step flow chart of calculating water body flow characteristics and hydraulic characteristics provided by an embodiment of the present application. Detailed Embodiments

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

[0017] Obtain and preprocess multi-source data to obtain an integrated dataset, including digital elevation model data, remote sensing images, river network vector data, temporal water body distribution information, and hydraulic engineering facility data;

[0018] Based on the integrated dataset, extract the river network skeleton and construct the river network topological structure;

[0019] Based on the integrated dataset and the river network topological structure, identify hydraulic engineering facilities and analyze their connectivity to generate a regional connectivity feature map;

[0020] Based on the river network topological structure and the regional connectivity feature map, perform hydraulics-driven adaptive partitioning of computational units to form plain water network hydrological computational units.

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

[0022] Analyze the temporal water body distribution information and extract stable water bodies;

[0023] Based on the integrated dataset, calculate the water body flow characteristics and hydraulic characteristics;

[0024] Utilize the water body flow characteristics and stable water bodies, and adopt the enhanced Markov random field method to extract the river network skeleton;

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

[0026] Based on the key river network nodes and the river network skeleton, adopt the invariant topological feature preservation algorithm to construct the river network topological structure.

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

[0028] Calculate the water body presence frequency at each spatial position in the temporal water body distribution information to generate a water body frequency grid;

[0029] Based on the water body frequency grid and the temporal water body distribution information, process using the Wavelet-SSA hybrid decomposition method to obtain the decomposed water body signal;

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

[0031] Utilize the decomposed water body signal and the optimized decomposition parameters, and adopt the adaptive regional segmentation SSA method for signal reconstruction to obtain the stable water body component and the seasonal water body component;

[0032] Perform stability level classification and spatial consistency processing on the stable water body component and the seasonal water body component to obtain the stability level classification map, and determine the stable water bodies accordingly.

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

[0034] Correct and finely register the remote sensing images to obtain a sequence of finely registered images;

[0035] Utilize the segmented robust optical flow algorithm to process the sequence of finely registered images and the temporal water body distribution information to calculate the apparent displacement field;

[0036] Based on DEM and the temporal water body distribution information, calculate the hydraulic gradient field;

[0037] The apparent displacement field and the hydraulic gradient field are subjected to feature fusion using multi-scale wavelet transform to obtain a multi-scale feature representation, and through inverse wavelet transform and feature integration, it is reconstructed into water flow characteristics and hydraulic characteristics.

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

[0039] Based on the hydrological calculation units in the plain water network, a multi-dimensional water volume exchange relationship between the calculation units is constructed to form an adaptive exchange model.

[0040] The obtained set of calculation units is combined with the adaptive exchange model to construct a complete hydrological model calculation framework;

[0041] Based on the measured hydrological data, the parameters of the constructed hydrological model are calibrated and verified to evaluate the rationality of the calculation unit division;

[0042] According to the verification result, the boundary of the calculation unit is fine-tuned as necessary to form the final calculation unit division scheme.

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

[0044] S1: Acquisition and preprocessing of multi-source data.

[0045] S11: Acquisition and enhancement of a high-precision digital elevation model (DEM).

[0046] Read the original DEM data, and process it using a method combining bidirectional filtering and multi-scale decomposition to obtain an enhanced DEM. Specifically, it includes: using non-local means filtering to remove noise, then applying wavelet transform for multi-scale decomposition to strengthen the micro-topography features, and finally obtaining a high-precision DEM with enhanced micro-topography features through geomorphic constraint reconstruction.

[0047] S12: Acquisition of multi-temporal remote sensing images and extraction of water body information.

[0048] Collect satellite remote sensing images of multiple time phases (including optical and SAR images), and extract the temporal water body distribution information through the dynamic threshold water body index (DTWI) method. Specifically, it includes: performing radiometric correction and geometric correction on the remote sensing images of different seasons, then calculating the modified normalized water body index (MNDWI), and establishing an adaptive threshold extraction algorithm by combining entropy feature analysis, and finally obtaining the water body spatial distribution map under different water conditions.

[0049] S13: Collection and structuring of river network vector data.

[0050] Obtain the regional river network vector data, and through topological consistency check and structured processing, form structured river network data with hierarchical attributes. Specifically, it 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 an upstream and downstream connection relationship table for the river network to form a complete river network data structure.

[0051] S14: Acquisition and attribute assignment of water conservancy project facility data.

[0052] Collect water conservancy project facility data (such as sluices, culverts, pumping stations, etc.) in the region, and establish a database of the operating characteristics of water conservancy projects through on-site research and historical operation records. Specifically, it includes: collecting information such as the spatial location, structural parameters, and operation rules of water conservancy projects, constructing a database containing attributes such as scheduling rules, water passing capacity, and working conditions, and providing data support for subsequent analysis of the impact of water conservancy projects on water flow.

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

[0054] Perform coordinate unification and spatial consistency processing on the above-mentioned multi-source heterogeneous data to obtain an integrated data set under a unified reference system. Specifically, it includes: unifying the projection coordinate system, performing spatial registration on different data sources, processing the boundary matching problem of data with different scales, and establishing a spatial index to form a unified data access interface to provide a consistent data environment for subsequent analysis.

[0055] S2: Extraction of river network skeleton and topological construction based on hydraulic characteristics.

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

[0057] Analyze the temporal distribution information of water bodies, and extract stable water bodies and seasonal water bodies through time-frequency characteristic analysis. Specifically, it includes: calculating the water body presence frequency of each pixel in the time series, applying singular spectrum analysis (SSA) to distinguish stable water bodies and seasonal water bodies, and generating a stability level classification map to provide a basis for river network skeleton extraction.

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

[0059] Based on multi-temporal remote sensing images, combined with enhanced DEM, adopt a method combining optical flow field analysis and hydraulic gradient calculation to extract water body flow characteristics and hydraulic characteristics. Specifically, it includes: calculating the apparent displacement field of water bodies through consecutive temporal remote sensing images, calculating the hydraulic gradient field in combination with DEM, and fusing to obtain a hydraulic characteristic field describing the dynamic characteristics of water flow to provide a basis for identifying the main water flow channels.

[0060] S23: Extraction of river network skeleton based on persistent homology.

[0061] Utilize the water body flow characteristics and stable water bodies, and adopt the persistent homology theory to extract the river network skeleton with topological stability. Specifically, it includes: regarding the hydraulic characteristic field as a high-dimensional manifold, constructing its filtration complex, calculating the persistence diagram to identify topologically significant features, and screening through the persistence threshold to form a river network skeleton with topological stability, effectively retaining the circular and branching features in the network structure.

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

[0063] Utilize the water body flow characteristics and stable water bodies, and adopt the Enhanced Markov Random Field (EMRF) method to extract the river network skeleton with topological stability. Specifically, it includes: constructing the hydraulic characteristic field as a Markov random field with a non-local energy function, introducing dual prior constraints based on fluidity and connectivity, designing an adaptive neighborhood structure to capture long-range dependencies, solving the optimization problem through the graph cut algorithm, extracting the river network skeleton that maintains the topological structure, and effectively retaining the circular and branching features in the network structure.

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

[0065] Based on the river network skeleton, adopt a hierarchical clustering method with dual geometric and hydraulic characteristics to identify and classify key river network nodes. Specifically, it includes: extracting characteristic nodes such as river network intersection points and bifurcation points, calculating the hydraulic importance index of each node in combination with hydraulic characteristics, classifying the nodes using spectral clustering, marking key nodes of different functional types, and forming a node classification table.

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

[0067] Based on the key river network nodes and the river network skeleton, adopt an invariant topological feature preservation algorithm to construct the river network topological structure. Specifically, it includes: representing the river network as a graph structure with attributes, defining topological invariants (such as Betti numbers, loop numbers, etc.) as topological features, maintaining these invariants during the graph simplification process, and connecting key nodes through the minimum energy path to form a river network structure model that maintains 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 the river network.

[0070] Spatially associate the facilities in the water conservancy project operation characteristics database with the river network topology to obtain coupled river network - water conservancy project data. Specifically, it includes: mapping water conservancy project facilities onto the river network using spatial proximity analysis, establishing the association relationship between water conservancy projects and river reaches, and forming a coupled data structure containing location information and connection relationships.

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

[0072] Based on the water conservancy project operation characteristics database, construct a formal expression of the operation rules of water conservancy projects to form an operation rule knowledge base. Specifically, it includes: analyzing the scheduling rules of water conservancy projects, establishing a formal rule expression containing conditions, actions, and constraints, designing an operation state prediction model based on rule reasoning, and realizing the mathematical description of the operation behavior of water conservancy projects.

[0073] S33: Dynamic connectivity analysis of water conservancy projects based on operation conditions.

[0074] Combining the coupled river network - water conservancy project data and the operation rule knowledge base, analyze the dynamic connectivity under different operation conditions to generate a connectivity state transition model. Specifically, it includes: defining the connectivity state space, establishing a state transition model based on Petri nets, analyzing the connectivity changes of the river network under different condition scenarios, and generating connectivity state transition rules.

[0075] S34: Quantification of connectivity uncertainty and calculation of risk probability.

[0076] Based on the connectivity state transition model, use the Bayesian network method to quantify the connectivity uncertainty and obtain the connectivity probability distribution. Specifically, it includes: identifying the key uncertain factors affecting connectivity, establishing a Bayesian network model to describe the conditional dependence relationship between factors, generating the connectivity probability distribution under different conditions through Monte Carlo simulation, and forming a risk probability map.

[0077] S35: Calculation of multi - scale connectivity indices.

[0078] Based on the connectivity probability distribution and the river network topology, calculate the multi - scale connectivity indices to form a regional connectivity characteristic map. Specifically, it includes: defining the index calculation methods for three scales of local connectivity, regional connectivity, and global connectivity, considering the influence weight of water conservancy projects, and generating a characteristic map representing the spatial distribution of regional hydrological connectivity to provide a basis for subsequent calculation unit division.

[0079] S4: Hydraulics - driven adaptive division of calculation units.

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

[0081] Based on the hydraulic characteristics and regional connectivity characteristic maps, the hydraulic response similarity clustering method is adopted to identify the preliminary hydraulic response units. Specifically, it includes: calculating the hydraulic response feature vectors of each grid cell in the region, applying the spectral clustering algorithm for response similarity analysis, and preliminarily dividing the response units with similar hydraulic behaviors, laying a foundation for fine division.

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

[0083] Analyze the boundary regions of the preliminary hydraulic response units, and based on the flow regime critical change point detection algorithm, identify the critical rheological points and optimize the unit boundaries. Specifically, it includes: calculating the hydraulic gradients and flow regime change characteristics of the boundary regions, identifying the critical points where significant changes in the flow regime occur (such as the transition points from rapid flow to slow flow, diversion points, etc.), using these points as the control points for boundary optimization, and generating unit boundaries that are more in line with the hydraulic characteristics.

[0084] S43: Unit adjustment driven by dynamic connectivity.

[0085] Combined with the connectivity state transition model, perform dynamic connectivity-based adjustment on the preliminary hydraulic response units to form connectivity-optimized response units. Specifically, it includes: analyzing the changes in connectivity between units under different working conditions, identifying the regions with significant connectivity changes, adjusting the unit boundaries to be consistent with the connectivity change boundaries, and ensuring that the unit division can reflect the dynamic connectivity characteristics.

[0086] S44: Refinement of unit boundaries under multi-objective constraints.

[0087] Based on the connectivity-optimized response units, establish a multi-objective optimization model including computational efficiency, physical consistency, and topological integrity, and perform refined adjustment of the unit boundaries. Specifically, it includes: defining the unit quality evaluation index, constructing the multi-objective optimization function, and adopting a hybrid optimization algorithm combining gradient descent and simulated annealing to finely adjust the unit boundaries, taking into account both computational efficiency while ensuring physical meaning.

[0088] S45: Unit merging and splitting based on hydrological process similarity.

[0089] Analyze the hydrological process similarity of the units after the refined adjustment of the unit boundaries, perform adaptive merging and splitting to form the final plain water network hydrological calculation units. Specifically, it includes: defining the similarity index including multiple hydrological processes such as runoff generation, confluence, and river channel evolution, calculating the process similarity matrix between units, adaptively merging adjacent units with high similarity according to the similarity threshold, and re-splitting the units with high internal heterogeneity, finally forming hydrological calculation units that balance physical meaning and computational efficiency.

[0090] S5: Construction of multi-dimensional water volume exchange relationships between units.

[0091] S51: Analysis of the hydraulic characteristics at the unit interface.

[0092] Analyze the interfaces between the hydrological calculation units in the plain water network, and calculate the hydraulic characteristic parameters of the interfaces. Specifically, it includes: extracting the geometric characteristics of the interfaces between units, calculating the hydraulic gradient, flow direction probability, and exchange capacity of the interfaces in combination with the hydraulic characteristic field, establishing a parameter set describing the hydraulic characteristics of the interfaces, and providing basic data for the subsequent construction of exchange relationships.

[0093] S52: Construction of an exchange network based on spectral graph theory.

[0094] Based on the hydraulic characteristic parameters and the river network topological structure, apply spectral graph theory to establish an exchange network for calculation units. Specifically, it includes: representing the calculation units as network nodes, representing the interfaces between units as edges, determining the edge weights based on the hydraulic characteristic parameters, and determining the key connections and main exchange paths of the exchange network through Laplacian matrix spectral analysis, forming a network structure describing the water volume exchange between units.

[0095] S53: Simplified expression of the multi-dimensional water volume exchange process.

[0096] Analyze the exchange process in the exchange network of calculation units, develop a simplified expression model for multi-dimensional water volume exchange, and form an exchange relationship expression. Specifically, it includes: decomposing the complex multi-dimensional exchange process into principal components, using orthogonal polynomial expansion technology to express the spatio-temporal variation characteristics, screening the main influencing factors through sensitivity analysis, establishing a computationally efficient simplified expression, and improving the computational efficiency while maintaining the physical meaning.

[0097] S54: Characterization of the uncertainty of the exchange relationship and parameter calibration.

[0098] Based on the exchange relationship expression, construct a Bayesian inference framework to quantify the parameter uncertainty and form a calibrated set of exchange parameters. Specifically, it includes: defining the prior distribution of the exchange parameters, using the observed data and the characteristics of the basin response, inferring the posterior distribution of the parameters through the Markov chain Monte Carlo method, obtaining a set of exchange parameters characterizing the uncertainty, and improving the robustness of the exchange calculation.

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

[0100] Based on the calibrated set of exchange parameters and the connectivity state transition model, construct an adaptive exchange mechanism with dynamic weights and form an adaptive exchange model. Specifically, it includes: designing an exchange weight function that automatically adjusts according to the water regime conditions and the operating conditions of water conservancy projects, establishing a weight update rule at multiple time scales, realizing the dynamic adaptability of the exchange relationship, and improving the adaptability of the model 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, specifically:

[0102] S211: Temporal water body presence frequency calculation.

[0103] Read the temporal water body distribution information, construct a three-dimensional data cube (x, y, t), where x and y are spatial coordinates and t is the time dimension. Calculate the water body presence frequency for each spatial position (x, y) to obtain the water body frequency grid. Specifically include: For each spatial position, count the number of times n that water bodies appear in N time phases, calculate the presence frequency f = n / N, and generate grid data with frequency values in the range [0, 1].

[0104] S212: Decomposition of water body temporal signal based on Wavelet-SSA.

[0105] Read the temporal water body distribution information and the water body frequency grid, and use the Wavelet-SSA hybrid decomposition method to process the temporal signal to obtain the decomposed water body signal. Specifically include: First, perform wavelet transform preprocessing on the temporal water body signal of each pixel to remove high-frequency noise; then construct a trajectory matrix and perform singular value decomposition; according to the energy distribution of the singular values, decompose the signal into a trend term, a seasonal term, and a noise term; output the preprocessed and decomposed water body signal components.

[0106] S213: Regional adaptive SSA parameter optimization.

[0107] Read the decomposed water body signal and the enhanced DEM, and optimize the SSA parameters through a regional segmentation strategy to obtain the optimized decomposition parameters. Specifically include: Based on the terrain partition derived from the DEM, divide the study area into several sub-regions; for each sub-region, determine the optimal window length and the SVD component selection threshold through multi-objective optimization of minimizing the reconstruction error and maximizing the signal separation degree; generate a spatially distributed SSA parameter set.

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

[0109] Read the decomposed water body signal and the optimized decomposition parameters, and use the adaptive regional segmentation SSA method for signal reconstruction to obtain the stable water body component and the seasonal water body component. Specifically include: Use the regional optimization parameters to perform grouped reconstruction on the decomposed water body signal; the trend term is identified as the stable water body component; the seasonal term is identified as the seasonal water body component; process the regional boundary transition problem through wavelet boundary processing technology; output the water body components with clear physical meanings.

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

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

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

[0113] S221: Preprocessing and registration of multi-temporal remote sensing images.

[0114] Read multi-temporal remote sensing images, perform radiometric correction, geometric correction, and fine registration processing to obtain a sequence of finely registered images. Specifically, it includes: using the histogram matching method for relative radiometric correction to eliminate atmospheric and sensor differences; using feature point matching and affine transformation for fine registration to control the error to be less than 0.5 pixels; generating a sequence of multi-temporal images with spatial alignment.

[0115] S222: Calculation of piecewise robust optical flow.

[0116] Read the sequence of finely registered images and the temporal water body distribution information, and use the piecewise robust optical flow algorithm to calculate the apparent displacement field to obtain the apparent displacement field. Specifically, it includes: dividing the study area into multiple sub-regions according to the water body distribution; for each sub-region, setting different boundary constraint conditions; applying the improved Horn-Schunck algorithm, introducing the L1 norm smoothing term and data term to improve the robustness to discontinuities; using a multi-resolution strategy to optimize the flow field from coarse to fine; outputting a two-dimensional vector field representing the apparent movement of the water body.

[0117] S223: Calculation of hydraulic gradient based on DEM.

[0118] Read the enhanced DEM and the temporal water body distribution information, calculate the hydraulic gradient field based on DEM to obtain the hydraulic gradient field. Specifically, it includes: performing hydrological optimization processing on the DEM to remove depressions; calculating the eight-direction flow direction and cumulative flow; calculating the hydraulic gradient based on empirical hydraulic formulas; screening effective regions in combination with the water body distribution information; outputting a gradient field representing potential water flow paths.

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

[0120] Read the apparent displacement field and hydraulic gradient field, and use multi-scale wavelet transform to achieve feature fusion to obtain multi-scale feature representation. Specifically, it includes: performing discrete wavelet transform on the two fields and decomposing them into sub-bands of different scales; fusing the corresponding sub-bands using an adaptive weighting method; determining the weights based on local consistency and gradient intensity; performing sub-band reconstruction; and outputting the multi-scale feature representation.

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

[0122] Read the multi-scale feature representation, and generate the final water body flow features and hydraulic features through inverse wavelet transform and feature integration. Specifically, it 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 features (such as flow velocity, water depth, Froude number, etc.); and outputting the hydraulic feature field characterizing the dynamic characteristics of the water body.

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

[0124] S241: Optimization and preprocessing of hydraulic feature field.

[0125] Read the water body flow features, hydraulic features and stable water bodies, perform optimization and preprocessing to obtain an enhanced hydraulic feature field. Specifically, it includes: normalizing and dealing with outliers of hydraulic features; screening effective regions based on the stable water body mask; applying anisotropic diffusion filtering to enhance linear features; calculating the main direction field of hydraulic features; and outputting the enhanced feature field for river network extraction.

[0126] S242: Construction of non-local energy function.

[0127] Read the enhanced hydraulic feature field, construct an energy function with non-local interactions to obtain the non-local energy function. Specifically, it includes: defining the river network label set as {0, 1}, representing non-river network and river network; constructing an energy function including a data term, a smoothing term and a non-local term; the data term is based on the hydraulic feature intensity; the smoothing term uses a direction-aware Potts model; the non-local term considers the interaction between pixels at a long distance, where the interaction intensity decays with distance but decays more slowly along the flow direction; and outputting the energy function characterizing the pixel classification probability.

[0128] S243: Definition of double prior constraints.

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

[0130] S244: Adaptive neighborhood structure design.

[0131] Read the enhanced hydraulic feature field and prior constraint conditions, design an adaptive neighborhood structure, and obtain an adaptive neighborhood system. Specifically, it includes: defining a direction-aware neighborhood based on local flow features; designing a multi-scale neighborhood system, including local neighborhoods and non-local connections; adaptively adjusting the neighborhood size and shape according to local features; defining a spatially varying potential function for each pixel; and outputting the neighborhood structure for EMRF modeling.

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

[0133] Read the non-local energy function, prior constraint conditions, and adaptive neighborhood system, solve the optimization problem through the graph cut algorithm, and obtain the river network skeleton. Specifically, it includes: constructing a weighted graph structure, where nodes represent pixels and edges represent the interaction between pixels; transforming the energy minimization problem into a graph cut problem; applying the α-expansion algorithm to solve the global optimal label assignment; performing morphological thinning and topological repair on the preliminary extraction results; and outputting the river network skeleton that maintains the topological structure.

[0134] According to one aspect of the present application, S42. Boundary optimization based on the critical rheological point is specifically as follows:

[0135] S421: Multidimensional hydraulic feature calculation.

[0136] Read the preliminary hydraulic response units, water body flow characteristics, and hydraulic features, calculate the multidimensional hydraulic features of the boundary region, and obtain the boundary hydraulic feature set. Specifically, it includes: extracting the unit boundary expansion region (5 pixels inside and outside); calculating multidimensional hydraulic features including water depth, flow velocity, water surface width, Froude number, hydraulic radius, etc.; normalizing the calculation results; and outputting the multidimensional hydraulic feature dataset of the boundary region.

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

[0138] Read the boundary hydraulic feature set, select the most discriminative feature combination using information entropy analysis, and obtain the optimal feature subset. Specifically, it includes: calculating the information entropy of each feature and the mutual information between features; screening the feature subset based on the maximum correlation minimum redundancy (mRMR) criterion; verifying the feature importance using the recursive feature elimination method; selecting the feature combination with the maximum information gain; and outputting the best feature subset for detecting the flow regime change points.

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

[0140] Read the optimal feature subset, detect the flow regime change points through multi-feature fusion technology, and obtain the candidate flow regime change point set. Specifically, it includes: constructing the trajectory curve in the feature space; applying the geometric method based on the change of curvature and tangent angle to detect the change points; using the CUSUM algorithm to detect the significant change of statistical characteristics; combining the spatial clustering based on DBSCAN to identify the change region; and outputting the candidate point set representing the significant change of the flow regime.

[0141] S424: Hydraulic significance verification and screening.

[0142] Read the candidate flow regime change point set and the hydraulic characteristics, conduct hydraulic significance verification and screening, and obtain the critical flow regime change points. Specifically, it includes: verifying the physical rationality of the candidate points based on classical hydraulic theories (such as the change of Froude number, hydraulic jump conditions, etc.); calculating the hydraulic importance index of each candidate point; setting the importance threshold to screen the key flow regime change points; performing spatial grouping and representative point selection for the flow regime change points; and outputting the critical flow regime change points with clear hydraulic significance.

[0143] S425: Boundary reconstruction based on the flow regime change points.

[0144] Read the critical flow regime change points and the preliminary hydraulic response units, reconstruct the calculation unit boundary, and obtain the optimized boundary of the flow regime change points. Specifically, it includes: converting the critical flow regime change points into boundary control points; constructing a smooth boundary curve using B-spline interpolation; solving the boundary conflict and topological consistency problems; applying the Snake model for boundary detail optimization; and outputting the unit boundary optimized based on the hydraulic characteristics to provide a basis for subsequent unit adjustment.

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

[0146] S531: Physical mechanism analysis of the exchange process.

[0147] Read the computational unit exchange network and hydraulic characteristic parameters, analyze the physical mechanism of water volume exchange, and obtain the classification of exchange mechanisms. Specifically, it includes: identifying the main exchange types (such as river channel exchange, overbank exchange, groundwater exchange, etc.); analyzing the driving factors and control equations of various exchanges; establishing the corresponding relationship between exchange types and physical parameters; and outputting the classification of the physical mechanism of the exchange process.

[0148] S532: Decomposition of the exchange process and extraction of principal components.

[0149] Read the classification of exchange mechanisms and hydraulic characteristic parameters, adopt the principal component analysis method, extract the main components of the exchange process, and obtain the exchange principal components. Specifically, it includes: constructing the parameter matrix of the exchange process; applying singular value decomposition to extract the principal components; analyzing the physical meaning and contribution degree of each principal component; retaining the principal components that explain more than 90% of the variance; and outputting the basic components with simplified expressions.

[0150] S533: Construction of orthogonal polynomial basis functions.

[0151] Based on the exchange principal components and hydraulic characteristic parameters, construct an orthogonal polynomial basis function system suitable for expressing water volume exchange, and obtain the orthogonal basis function set. Specifically, it includes: selecting a suitable orthogonal polynomial family according to the physical characteristics of the exchange process (such as Legendre polynomials, Chebyshev polynomials, etc.); customizing the basis functions for different types of exchange processes; ensuring the completeness and orthogonality of the basis function system; and outputting the basis function set for expressing the exchange process.

[0152] S534: Optimization of polynomial coefficients and generation of expressions.

[0153] Combine the orthogonal basis function set and the exchange principal components, optimize the polynomial coefficients, and generate a simplified exchange relationship expression. Specifically, it includes: designing an objective function based on physical constraints; using the least squares method to determine the polynomial coefficients; evaluating the accuracy of the expression through cross-validation; balancing accuracy and complexity, and controlling the polynomial order; and outputting a simplified exchange relationship expression with physical meaning.

[0154] S535: Verification and optimization of expression efficiency.

[0155] Conduct computational efficiency tests and optimizations on the exchange relationship expressions to obtain efficient exchange expressions. Specifically, it includes: analyzing the computational complexity of different expressions; comparing the measured computational times; optimizing the code for the frequently called parts; introducing techniques such as lookup tables to accelerate the calculation; and outputting optimized expressions that improve computational efficiency while maintaining accuracy.

[0156] Example 1: It mainly describes the process of extracting the river network skeleton based on the enhanced Markov random field. In this example, for a typical sub-region (with an area of about 800 square kilometers) in the plain water network area of the Taihu Lake Basin, the river network skeleton extraction method of the present invention is applied to solve the problem of maintaining topological relationships in the processing of complex river network nodes.

[0157] First, obtain multi-source data of the study area: High-resolution digital elevation model (DEM) data: DEM data with a resolution of 5 meters; Multi-temporal remote sensing images: 12 Sentinel-2 satellite images from January to December 2022 (resolution 10 meters); Existing river network vector data: River network vector map at a scale of 1:50,000; Water conservancy project facility data: Spatial location and attribute information of 42 sluices, 17 pumping stations, and 28 culverts within the study area; Through coordinate unification and spatial registration, an integrated dataset in the unified reference system (CGCS2000 coordinate system) is formed. Among them, non-local means filtering is applied to the DEM data for noise removal, and multi-scale decomposition is performed through wavelet transform to enhance micro-topographic features, obtaining an enhanced DEM (E-DEM).

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

[0159] Extract the water body distribution of the 12 images through the dynamic threshold method to form the temporal water body distribution information W(x,y,t), where (x,y) is the spatial coordinate and t is the time series (1 to 12). Calculate the water body presence frequency for each spatial position: f(x,y) = Σ t=1 12 W(x,y,t) / 12; Generate the water body frequency raster F.

[0160] Adopt the Wavelet-SSA hybrid decomposition method to process the temporal water body distribution information:

[0161] First, perform discrete wavelet transform on W(x,y,t) to obtain the denoised water body signal Wd(x,y,t). Then construct the trajectory matrix X, for each spatial position (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 (take 4) and K = 12 - L + 1. Perform singular value decomposition on X(x,y): X(x,y) = UΣV T ;

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

[0163] Apply the corresponding optimized parameters to each sub-region, and divide the decomposition results into a trend component Wt(x, y) (stable water body) and a seasonal component Ws(x, y) (seasonal water body). Apply wavelet boundary processing technology to handle the sub-region boundary transition problem to ensure the spatial continuity of the components.

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

[0165] For the 12 registered Sentinel-2 images, use the piecewise robust optical flow algorithm to calculate the apparent displacement field of the water bodies. For two adjacent images It and It+1, calculate the optical flow field V(u, v), where u and v are the displacement components in the x and y directions respectively:

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

[0167] Introduce a robust 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 weight of the smoothing term (take 0.05).

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

[0170] Fuse the features of the apparent displacement field V and the hydraulic gradient field G through multi-scale wavelet transform. First, apply the discrete wavelet transform to V and G, and decompose them into sub-bands of 4 scales 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] Among them, A represents the approximation sub-band, D represents the detail sub-band, and the superscripts h, v, and d represent the horizontal, vertical, and diagonal directions respectively. ^h means h is the superscript, ^ represents the superscript, and the same expression is used in the following formulas.

[0172] The corresponding sub-bands are fused using an adaptive weighted method: 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 is calculated based on local consistency and gradient intensity ω 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 through inverse wavelet transform, which contains hydraulic information such as flow direction and flow velocity.

[0174] Based on the stable water body 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 outlier processed to generate the enhanced hydraulic characteristic field He(x, y).

[0175] Define the river network label set S = {0, 1}, representing non-river network and river network. Construct an energy function E(S) that includes a data term, a smooth term, and a non-local term: E(S) = Σi[Ed(si) + Σj∈Ni Es(si,sj) + Σk∈Ri\Ni Enl(si,sk)]; where the data term Ed(si) = -log p(He(i)|si), representing the probability that pixel i is classified as si, based on the intensity of He(i); the smooth term Es(si,sj) = λs·δ(si≠sj)·g(grad He(i),grad He(j)), where λs is the weight coefficient (taking 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 (taking 0.4), and wik is the non-local interaction intensity wik = exp(-d(i,k) / σd)·exp(-θ(i,k) / σθ); d(i,k) is the Euclidean distance between pixel i and k, θ(i,k) is the difference in flow direction between i and k, and σd and σθ are scale parameters.

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

[0177] Define the mobility prior constraint Pf(S) = Σi[δ(si=1)·f(i)] based on the flow direction of hydraulic characteristics; where f(i) represents the mobility 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 the 8-neighborhood or 16-neighborhood, adaptively adjusted according to the local flow direction;

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

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

[0182] Solve this optimization problem through the α-expansion graph cut algorithm to obtain the preliminary extraction result S* of the river network skeleton. Apply the topological repair algorithm to process breakpoints and hanging points to ensure the connectivity of the river network, and finally obtain the river network skeleton R.

[0183] Based on the river network skeleton R, a feature point detection algorithm with dual geometric and hydraulic characteristics is used to identify key river network nodes N = {n1, n2, ..., nm}, including intersection points, bifurcation points, etc. Calculate the hydraulic importance index HI of each node: HI(ni) = α·Dc(ni) + β·F(ni) + γ·A(ni); where Dc is the connectivity, F is the estimated average flow, A is the estimated catchment area, and α, β, γ are weight coefficients (taking 0.3, 0.4, 0.3 respectively).

[0184] Based on the key nodes N and the river network skeleton R, an invariant topological feature preservation algorithm is used to construct the river network topological structure T. Define the set of topological invariants: 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 key links, and BR is the bifurcation ratio.

[0185] During the river network simplification process, 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, compared with the traditional method, the river network skeleton extracted by the enhanced Markov random field method has significant advantages in dealing with complex nodes: the accuracy of river network connectivity has increased by 18.7%; the topological structure preservation rate has increased by 22.3%; the accuracy of node recognition in the densely crossed area 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 bifurcation points and intersections, this method can effectively maintain the accuracy of the water system topological relationship, laying a solid foundation for the subsequent calculation unit division.

[0188] Embodiment 2: Mainly describes the process of dynamic connectivity analysis of water conservancy projects based on Petri nets. This embodiment is aimed at a plain river network area with typical water conservancy project control characteristics (about 500 square kilometers, including 28 sluices and 12 pumping stations), and implements a method for dynamic connectivity analysis of water conservancy projects based on Petri nets.

[0189] Spatially associate the water conservancy project facility data with the river network topological structure constructed in Embodiment 1, and use the nearest distance method and buffer analysis (the buffer distance is set to 30 meters) to establish the corresponding relationship between water conservancy projects and river reaches, forming the river network - water conservancy project coupling data RH.

[0190] For each water conservancy project facility h, record its position coordinates, the river reach 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 operation characteristic database of water conservancy projects, construct a formal expression of the operation rules of water conservancy projects. For each water conservancy project h, define the operation rule set R(h) = {r1, r2, ..., rk}; 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 a comparison operator (>、<、=、≠), wl_th is the water level threshold, T is the time, [t_start, t_end] is the applicable period, CMD is the scheduling instruction, and cmd_type is the instruction type.

[0193] The action set A includes operations such as gate opening adjustment and pump station start / 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 / stop status (0 or 1).

[0194] The constraint set CT includes operation constraint conditions: 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, establish a complete operation rule knowledge base KB.

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

[0197] First, divide the water conservancy project facilities into 4 types of functional subsystems: control gates (18); diversion gates (10); drainage pump stations (12); channel culverts (8);

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

[0199] Organize the subsystem models using a hierarchical Petri net structure. 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 set of subnets; HT is the global transition set; HF is the inter-layer connection relation; CC is the color set, using different colors to mark different types of connectivity relations.

[0200] Through Petri net simulation analysis, generate a connectivity state transition model CST. For each pair of adjacent computing units (i, j), define the connectivity state set CS(i, j) = {cs1, cs2, ..., csv}; where each connectivity state csk corresponds to a specific engineering operating condition, and the condition conversion between operating conditions is performed through the state transition function ST: ST: CS(i, j) × CON → CS(i, j); CON is the set of conditions, including water level conditions, time conditions, and manual scheduling conditions.

[0201] Based on the connectivity state transition model CST, use the Bayesian network method to quantify the connectivity uncertainty. Build a Bayesian network BN = {V, E, CPT}; where V = {v1, v2, ..., vl} is the set of nodes, including water conservancy project state nodes, hydrological condition nodes, and connectivity state nodes; E is included in V × V as the set of directed edges, representing the conditional dependence relationship; CPT is the conditional probability table, defining the probability distribution of each node under the values of its parent nodes.

[0202] For each pair of adjacent computing units (i, j), calculate the probability distribution of different connectivity states P(CS(i, j)) = {P(cs1), P(cs2), ..., P(csv)}; through the Monte Carlo simulation method, generate a connectivity probability distribution map based on different hydrological scenarios and scheduling schemes.

[0203] Based on the connectivity probability distribution and river network topological structure, calculate the connectivity indices at multiple scales to form the regional connectivity characteristic map CI.

[0204] Define the connectivity indices at three scales: Local Connectivity Index LCI(i,j): The direct connectivity between adjacent cells i and j; Regional Connectivity Index RCI(R): The internal connectivity within region R; Global Connectivity Index GCI: The overall connectivity of the entire study area.

[0205] Among them, the Local Connectivity Index LCI(i,j) = Σk[P(csk) × W(csk)]; W(csk) is the weight of the connectivity state csk, determined based on the hydraulic flux capacity.

[0206] The Regional Connectivity Index RCI(R) = [Σi,j∈R LCI(i,j)] / [n(n - 1) / 2]; n is the number of calculation cells within region R.

[0207] The Global Connectivity Index GCI = λ1·DCI + λ2·FCI + λ3·TCI; where DCI is the connectivity index during the flood period, FCI is the connectivity index during the dry period, TCI is the connectivity index during the transition period, and λ1, λ2, λ3 are the time period weight coefficients (taking 0.4, 0.3, 0.3 respectively).

[0208] In this embodiment, the dynamic connectivity of water conservancy projects is analyzed through the Petri net model to form the regional connectivity characteristic map. Compared with the traditional static connectivity analysis method: it can identify and quantify the connectivity changes under 9 typical working conditions; the prediction accuracy of the connectivity state reaches 87.6%; the adaptability to sudden scheduling scenarios is improved by 23.4%; it provides a reliable basis for the dynamic connection of calculation cells; especially in the area with dense gates, this method can accurately reflect the changes in the water system connectivity relationship under different operating conditions, significantly improving the accuracy of the hydrological model's response to artificial regulation in the plain water network area.

[0209] Example 3 mainly describes the process of boundary optimization based on the detection of critical change points of flow patterns. This example is aimed at the plain river network area (about 300 square kilometers) with gentle terrain (average slope <0.5‰), and implements the calculation unit boundary optimization method based on the detection of critical change points of flow patterns.

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

[0211] First, calculate the hydraulic response feature vector HRV(x, y) = [h(x, y), v(x, y), d(x, y), f(x, y), c(x, y)] for each grid cell within the region; 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 feature vector to obtain NHRV(x, y). Conduct response similarity analysis using the spectral clustering algorithm: construct a similarity matrix S, with elements sij = exp(-||NHRV(i)-NHRV(j)|| 2 / 2σ 2 ); calculate the Laplacian matrix L = D - S, where D is a diagonal matrix with dii = Σj sij; solve the generalized eigenvalue problem Lx = λDx to obtain the first k eigenvectors; apply k-means clustering in the feature space to obtain the preliminary hydraulic response unit HRU0.

[0213] Analyze the boundary region of the preliminary hydraulic response unit HRU0, and based on the flow state critical change point detection algorithm, identify the critical rheological points and optimize the unit boundary.

[0214] Extract the unit boundary expansion region BE (5 pixels inside and outside each), and calculate the multi-dimensional hydraulic feature set BHF = {WD, V, WW, FR, HR} for the boundary region; 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] Apply the feature selection method based on information entropy to calculate the information entropy H(X) and the mutual information I(X; Y) between features: H(X) = -Σx p(x)log p(x); I(X; Y) = Σx,y p(x,y)log[p(x,y) / (p(x)p(y))]; based on the maximum correlation and minimum redundancy (mRMR) criterion, select the most discriminative feature combination OFS = arg max[Σi I(Xi; C) - 1 / |S|·Σi,j∈S I(Xi; Xj)]; where C is the classification variable (representing different flow states) and S is the selected feature set.

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

[0217] Based on OFS, detect the flow regime change points through multi-feature fusion technology. Construct the trajectory curve TR(t) = [FR(t), V(t), HR(t)] in the feature space along the boundary coordinate parameter t. Calculate the trajectory curvature κ(t) and the tangent angle change rate κ(t)= ||TR'(t)×TR''(t)|| / ||TR'(t)|| 3 τ(t) = <TR'(t),TR'(t+1)> / (||TR'(t)||·||TR'(t+1)||).

[0218] Apply the CUSUM algorithm to detect significant changes in statistical characteristics S(t) = max[0, S(t-1) + (X(t) - μ0- K)]. Where X(t) is the feature sequence, μ0 is the reference mean, and K is the sensitivity parameter (take 0.5).

[0219] When S(t) exceeds the threshold H (take 5.0), it is marked as a change point. Combine spatial clustering based on DBSCAN (neighborhood radius ε = 3, minimum number of points MinPts = 4) to identify the change region and form the candidate rheological change point set CFP.

[0220] Verify the physical rationality of the candidate points based on the hydraulic principle. Calculate the hydraulic importance index HII(p) = α1·κ(p) + α2·τ(p) + α3·ΔFR(p) + α4·CoV(p). Where, ΔFR is the Froude number jump value, CoV is the velocity coefficient change, and α1 to α4 are weight coefficients (take 0.3, 0.25, 0.3, 0.15 respectively). Set the importance threshold (HII>0.6) to screen the key rheological change points CP, use them as boundary control points, and reconstruct the unit boundary through B-spline interpolation to obtain the optimized boundary of rheological change points HRU1.

[0221] Combine the connectivity state transition model CST in Example 2 to make dynamic connectivity-based adjustments to the optimized boundary of rheological change points HRU1.

[0222] Analyze the connectivity changes between units under different working conditions. Define the connectivity change index CVI(i,j) = 1 - min_k,l[sim(CS_k(i,j), CS_l(i,j))]; where, sim is the connectivity state similarity, and CS_k represents the connectivity state under working condition k. Identify the regions with significant connectivity changes (CVI>0.5), adjust the unit boundary to be consistent with the connectivity change boundary, and form the connectivity optimized response unit HRU2.

[0223] Based on the connectivity optimized response unit HRU2, establish a multi-objective optimization model for fine adjustment of the unit boundary.

[0224] Define the unit quality evaluation index \(Q = w_1\cdot Q_C + w_2\cdot Q_P + w_3\cdot Q_T\); where \(Q_C\) is the computational efficiency index, \(Q_P\) is the physical consistency index, \(Q_T\) is the topological integrity index, and \(w_1\) to \(w_3\) are weight coefficients (taking 0.3, 0.4, and 0.3 respectively).

[0225] Construct a multi-objective optimization function \(\min[-Q(HRU)]\ s.t.\ A(HRU)\geq A_{min}\ C(HRU)\leq C_{max}\ T(HRU) = T_0\); where \(A\) is the unit area constraint, \(C\) is the complexity constraint, and \(T\) is the topological invariant constraint.

[0226] Adopt a hybrid optimization algorithm (combining gradient descent and simulated annealing) to finely adjust the unit boundary to form a boundary-refined unit HRU3.

[0227] Analyze the hydrological process similarity of HRU3, and perform adaptive merging and splitting to form the final plain water network hydrological calculation unit HRU.

[0228] Define the hydrological process similarity index \(HPS(i,j)=\beta_1\cdot sim(R_i,R\) j )+\beta_2\cdot sim(C_i,C\) j )+\beta_3\cdot 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 process, \(sim\) is the similarity function, and \(\beta_1\) to \(\beta_3\) are weight coefficients (taking 0.35, 0.35, and 0.3 respectively).

[0229] Calculate the process similarity matrix SPM between units, and adaptively merge adjacent units with high similarity according to the similarity threshold (\(HPS > 0.75\)). Re-segment the units with high internal heterogeneity (coefficient of variation \(CV>0.25\)) to finally form a hydrological calculation unit HRU that balances physical meaning and computational efficiency.

[0230] In this embodiment, the boundary optimization method based on the detection of the critical change point of the flow regime, compared with the traditional terrain-based unit division method: the boundary recognition accuracy in the flat area (slope < 0.5‰) is increased by 32.4%; the number of calculation units is reduced by 17.3%, while maintaining the simulation accuracy; the Nash efficiency coefficient of the flow simulation is increased by 0.11; the root mean square error of the water level simulation is reduced by 0.08 m. Especially in the low slope drop area, this method can adaptively adjust the calculation unit boundary according to the water flow characteristics, significantly improving the accuracy of the hydrological calculation unit division.

[0231] Example 4 describes the complete process of the hydrological model calculation unit division method in the plain water network area. In this example, for a plain water network area in the lower reaches of the Yangtze River (with an area of about 1,200 square kilometers), the complete technical solution of the present invention is applied to realize the division of the hydrological model calculation units 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.6 km / km 2 ), and numerous water conservancy projects (including 63 sluices, 26 pumping stations, and 37 culverts).

[0232] S1: Acquisition and preprocessing of multi-source data.

[0233] Obtain the original DEM data (resolution 5 meters) of the study area at a scale of 1:10,000, and process it using the bidirectional filtering and multi-scale decomposition method:

[0234] Apply non-local mean filtering to remove DEM noise, set the filtering window size to 7×7 pixels, the similar window to 3×3 pixels, and the filtering intensity parameter h = 10; perform wavelet transform multi-scale decomposition on the filtered DEM, use the "db4" wavelet basis function, and the decomposition level is 4; perform adaptive enhancement on the wavelet high-frequency components, and the enhancement coefficient is related to the terrain undulation degree: α(x,y) = 1 + γ·var(z), where var(z) is the local elevation variance and γ = 5; through geomorphic constraint reconstruction, use the known river channel position information to correct the DEM to form an enhanced DEM (E-DEM);

[0235] Collect 12 Sentinel-2 satellite images (resolution 10 meters) from January to December 2023 and 4 Sentinel-1 SAR images (resolution 10 meters), and extract the water body distribution information: perform atmospheric correction on the optical images using the DOS (DarkObject Subtraction) method; perform scattering correction and geometric correction on the SAR images, with the reference image being the E-DEM; calculate the modified normalized difference water index MNDWI = (G - SWIR) / (G + SWIR); combine the SAR backscattering coefficient to establish a water body extraction decision tree model based on entropy features; use the OTSU adaptive threshold method to determine the optimal segmentation threshold to extract the water body spatial distribution of different time phases.

[0236] Obtain the river network vector data at a scale of 1:50,000, and perform topological consistency check and structuring processing: check the river network connectivity, identify and repair 185 breakpoints and 76 hanging points; assign level attributes to the river reaches based on the Strahler grading method, and a total of 6 grades are divided; establish a connection relationship table between the upstream and downstream of the river reaches to record the connection relationships of 458 river reaches; assign hydraulic characteristic parameters (river width, water depth, etc.) to each river reach.

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

[0238] Perform coordinate unification and spatial consistency processing on multi-source heterogeneous data: Uniformly adopt the CGCS2000 coordinate system and Gauss-Kruger projection (central meridian 121°E); perform spatial registration on different data sources, with 68 control points and a registration error < 1.0 meter; handle the boundary matching problem of data with different scales and adopt a boundary smooth transition algorithm; establish a spatial index (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] For details, please refer to Embodiment 1, which mainly includes: calculating the water body presence frequency f(x, y) of each pixel; processing the time-series water body signal using the Wavelet-SSA hybrid decomposition method; optimizing the SSA parameters based on the terrain partition of E-DEM; reconstructing the stable water body and seasonal water body using the adaptive regional segmentation SSA method; performing stability level classification to obtain the classification map C(x, y).

[0241] Perform fine registration processing on multi-temporal remote sensing images; calculate the apparent displacement field V using the segmented robust optical flow algorithm; calculate the hydraulic gradient field G based on E-DEM; perform feature fusion through multi-scale wavelet transform; 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 the dual prior constraints based on fluidity and connectivity; design an adaptive neighborhood structure; solve the optimization problem through the graph cut algorithm to extract the river network skeleton R.

[0243] Based on the river network skeleton R, identify and classify key river network nodes: Apply the Harris corner detection algorithm to initially identify feature points and obtain the candidate node set CN; Use the morphological thinning algorithm to extract the topological skeleton Rt of the river network skeleton; Calculate the connectivity of each pixel based on the skeleton Rt, and identify the intersection points and bifurcation points with a connectivity ≥ 3; Calculate the hydraulic importance index HI(ni) of the node: 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; According to the importance index and connection characteristics, classify the nodes into 4 categories: main trunk intersection points, main branch confluence points, ordinary bifurcation points, and end points.

[0244] Based on the key river network nodes and the river network skeleton, construct the river network topological structure: Represent the river network as a graph structure G = (V, E) with attributes, where the nodes V represent key points and the edges E represent river segments; Define the set of topological invariants: TV = {b0, b1, CL, BR}; where b0 is the number of connected components, b1 is the number of loops, CL is the set of key links, and BR is the branch ratio; Design a topological preservation energy function: 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 (taking 1.0, 0.8, 0.6, 0.4 respectively); Minimize ET during the graph simplification process to ensure the preservation of the topological structure; Calculate the minimum energy path between key nodes through the Dijkstra algorithm to form the complete river network topological structure T.

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

[0246] Specifically, see Embodiment 2. Spatially associate the water conservancy project facility data with the river network topological structure: Associate the water conservancy project points with the river network line elements using the nearest distance method; Set a buffer zone (30 meters) for spatial matching to solve the position deviation problem; Record the upstream and downstream river segment IDs and associated nodes of each water conservancy project; Establish an association data table RH between water conservancy projects and river segments, containing 126 rows of records;

[0247] Analyze the water conservancy project operation rules, extract conditions, actions, and constraints; Construct a formal rule expression: ri = {C, A, CT}; Define a standardized operation rule template for each type of water conservancy project; Establish a complete operation rule knowledge base KB.

[0248] Divide the water conservancy project facilities into functional subsystems; construct the basic Petri net model PN(h) = {P, T, F, W, M0}; organize the subsystem models using a hierarchical Petri net structure; design a colored Petri net interaction protocol; generate a connectivity state transition model CST.

[0249] Based on the connectivity state transition model, quantify the connectivity uncertainty: identify the key uncertain factors 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 variation (coefficient of variation CV = 0.08);

[0250] Establish a Bayesian network model BN = {V, E, CPT} to describe the conditional dependence relationship between factors; define the node conditional probability table, 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; perform parameter posterior distribution inference through the Markov chain Monte Carlo (MCMC) method: the number of samples N = 10000; set the burn-in period to 1000; control the acceptance rate between 0.23 - 0.44; generate the connectivity probability distribution P(CS(i,j)) under different water conditions.

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

[0252] S4: Adaptive division of computational units driven by hydraulics.

[0253] Based on the hydraulic characteristics and regional connectivity characteristic maps, identify the preliminary hydraulic response units: Calculate the hydraulic response feature vector for each grid cell: 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; Normalize the feature vector to obtain NHRV(x,y); Construct the similarity matrix S, with elements sij = exp(-||NHRV(i)-NHRV(j)|| 2 / 2σ 2 ); Calculate the normalized Laplacian 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 feature space to obtain the preliminary hydraulic response units HRU0, which are divided into 63 units in total.

[0254] For details, please refer to Example 3. Calculate the multi-dimensional hydraulic characteristics of the boundary region; Adopt a feature selection method based on information entropy; Detect the flow state change points through multi-feature fusion technology; Conduct hydraulic significance verification and screening; Reconstruct the calculation unit boundary.

[0255] Combine the connectivity state transition model to adjust the unit boundary: Analyze the connectivity changes between units under different working conditions, and define the connectivity change index CVI; Apply the watershed tracing algorithm to determine the connectivity boundary; Use morphological opening and closing operations to smooth the boundary contour; Adjust the unit boundary at significant connectivity change points (CVI>0.5); Form the connectivity-optimized response units HRU2, and the number of units is adjusted to 57.

[0256] Establish a multi-objective optimization model to conduct refined adjustment of the unit boundary: Define the unit quality evaluation index Q = 0.3·QC + 0.4·QP + 0.3·QT; Calculate the 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 state; The topological integrity index QT = h(Nc, Nr), where Nc is the number of truncated connectivity paths and Nr is the proportion of the retained river network skeleton;

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

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

[0259] Analyze the similarity of the hydrological processes of the units, and perform adaptive merging and splitting: 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 ) ; The similarity of the runoff generation process sim(R_i,R j ) uses the reciprocal of the KL divergence; the similarity of the confluence process sim(C_i,C j ) uses the Euclidean distance of the unit hydrograph; the similarity of the river channel evolution sim(F_i,F j ) uses the Manning coefficient and the cross-section shape similarity.

[0260] Calculate the process similarity matrix SPM between the units; use the hierarchical clustering algorithm (Ward method) to adaptively merge adjacent units with high similarity (HPS>0.75); apply the watershed algorithm to re-segment the units with high internal heterogeneity (CV>0.25); finally form a hydrological calculation unit HRU that balances physical meaning and computational efficiency, with 46 units.

[0261] S5: Calculate and construct the multi-dimensional water volume exchange relationship between the units.

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

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

[0264] Combine the hydraulic characteristic field to calculate the hydraulic parameters of the interface: 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; establish a parameter set HPI(i,j) that describes the hydraulic characteristics of the interface.

[0265] Based on the hydraulic characteristic parameters and river network topological structure, 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)), where f(CI(i,j)) is an adjustment function based on the connectivity index; a weighted adjacency matrix W and a degree matrix D are constructed; the Laplacian 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 exchanges is established.

[0266] For details, please refer to Example 2 to analyze the physical mechanism of the exchange process; use principal component analysis to extract the exchange principal components; construct an orthogonal polynomial basis function system; optimize the polynomial coefficients; conduct computational efficiency tests and optimizations.

[0267] First, the physical mechanism of water volume exchange is analyzed in detail to identify the main exchange types: River channel exchange (accounting for 65% of the total exchange volume): direct water volume exchange between units connected by the river network; Overflow exchange (accounting for 21% of the total exchange volume): water volume exchange through the overflow area during the flood period; Groundwater exchange (accounting for 14% of the total exchange volume): slow water volume exchange through the aquifer.

[0268] For each exchange type, its driving factors and control equations are analyzed. For example, the control equation for river channel exchange: 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 of the water passage, n is the Manning coefficient, R is the hydraulic radius, and S is the slope. Using the principal component analysis method, a parameter matrix M of the exchange process is constructed: 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. The singular value decomposition of M is performed: M = USV^T. Analyze the singular value decay curve and retain the first k principal components (k = 4 in this example) that explain more than 90% of the variance.

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

[0270] For different exchange types, construct multi-dimensional 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, cross-section parameters, connectivity status).

[0271] Based on the principal components and orthogonal basis functions, construct a polynomial expression for the exchange process: 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 (taking 16 - 25, controlled according to complexity).

[0272] Use the least squares method with physical constraints to determine the polynomial coefficients: min||M·a - Y|| 2 + λ·||a|| 2 + μ·Σ_i C_i(a) where Y is the exchange flow rate, λ is the regularization parameter (taking 0.01), μ is the constraint weight (taking 1.0), and C_i is the physical constraint condition. Evaluate the accuracy of the expression through cross-validation (5-fold cross-validation, RMSE = 0.037m 3 / s). Use code optimization and lookup table techniques to improve the calculation efficiency, construct an efficient exchange expression, and the calculation speed is increased by 42%.

[0273] Based on the exchange relationship expression, quantify the parameter uncertainty: Define the prior distribution of the exchange parameters: Manning coefficient n of the river channel ~ 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 the observed data and the characteristics of the basin response, adopt the Bayesian inference framework: P(θ|D) ∝ P(D|θ)·P(θ) where θ is the parameter vector and D is the observed data; Perform posterior distribution inference of the parameters through the adaptive Markov chain Monte Carlo method (AM-MCMC): the chain length is 50000, and the burn-in period is 10000; Adaptive step size adjustment, with a target acceptance rate of 0.234; Gelman-Rubin convergence diagnosis R < 1.1; Generate an exchange parameter set EPS representing uncertainty; Calculate the 95% confidence interval of the exchange flow rate and evaluate the prediction uncertainty.

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

[0275] Design an exchange weight function that automatically adjusts according to water conditions and engineering operating 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 basic weight; f_wl is the water level response function; f_cs is the connectivity state response function; f_s is the season response function;

[0276] Establish weight update rules for multiple time scales: Short-term response (hourly level): Based on real-time water level and gate opening; Medium-term response (daily level): Based on connectivity state conversion; Long-term response (monthly level): Based on seasonal pattern changes;

[0277] Develop a weight smooth transition algorithm to avoid sudden changes in exchange relationships; Implement an Adaptive Exchange Parameter Table (AEPM) that contains pre-computed parameter values for different scenarios; Form the final Adaptive Exchange Model AEM that can dynamically adjust exchange relationships according to real-time conditions;

[0278] After completing the calculation unit division and the construction of exchange relationships, verify and evaluate through two typical flood events in June 2022 and July 2023:

[0279] Evaluate the effect of maintaining topological relationships: The accuracy rate of river network connectivity reaches 94.2%; The retention rate of topological structure reaches 96.8%; The correct recognition rate of key nodes is 93.1%;

[0280] Expression effect of dynamic connectivity of water conservancy projects: The prediction accuracy rate of connectivity state reaches 89.4%; The recognition rate of connectivity changes for 9 typical operating conditions is 92.7%; The adaptability to sudden scheduling scenarios is improved by 26.2% compared with traditional methods;

[0281] Effect of unit division driven by hydraulics: The accuracy rate of boundary recognition is increased by 35.7% in flat areas (slope < 0.3‰); The number of calculation units is reduced from 83 in the traditional method to 46; The simulation accuracy is equivalent to or better than that of the traditional method (the Nash efficiency coefficient is increased by 0.13);

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

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

[0284] Regarding the problem of maintaining topological relationships in areas with dense and intertwined river networks, the present invention uses an enhanced Markov random field method to extract the river network skeleton. Through the innovative applications of constructing a non-local energy function, defining dual prior constraints, and designing an adaptive neighborhood structure, the accurate extraction of the river network structure is achieved. In particular, by combining hydraulic characteristics with topological constraints, the accurate identification and retention of key nodes during the river network simplification process are ensured. At the same time, an invariant topological feature preservation algorithm is used to construct the river network topological structure, effectively solving the problem of maintaining topological relationships in the processing of complex river network nodes and ensuring the accurate expression of water flow paths.

[0285] Regarding the problem of expressing dynamic connectivity under the regulation of water conservancy projects, the present invention innovatively constructs a dynamic connectivity model for different operating conditions based on Petri nets. The water conservancy project facilities are divided into functional subsystems. In the Petri net model, places are defined to represent connectivity states, transitions are defined to represent condition conversion conditions, and tokens are defined to represent current state indicators. And the uncertainty of connectivity is quantified through the Bayesian network method. This method can accurately describe the changes in connectivity states under different operating conditions of water conservancy projects, realize the conditional dynamic connection between computing units, and effectively solve the problem of expressing the impact of manual regulation on water system connectivity.

[0286] Regarding the problem of inaccurate boundary determination in low slope regions, the present invention proposes a method for identifying critical rheological points based on a flow regime critical change point detection algorithm. By calculating the multi-dimensional hydraulic characteristics of the boundary region of the preliminary hydraulic response unit and applying feature selection and multi-feature fusion techniques based on information entropy, the critical points where significant changes in the flow regime occur are accurately identified and used as control points for adjusting the unit boundary. This hydraulic-driven boundary optimization method breaks through the limitations of traditional terrain-driven methods and can finely adjust the boundary according to hydrodynamic characteristics in flat terrain areas, significantly improving the physical rationality and simulation accuracy of computing unit division.

[0287] The preferred embodiments of the present invention have been described in detail above. However, the present invention is not limited to the specific details in the above embodiments. Within the scope of 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 protection scope 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 integrated data sets, including digital elevation model data, remote sensing images, river network vector data, time-series water 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 hydraulically driven calculation unit adaptive division is carried out to form the plain water network hydrological calculation unit.

2. The method according to claim 1, characterized in that The steps to extract the river network skeleton and construct the river network topology structure include: Analyze the 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 key river network nodes and river network skeleton, the river network topology structure is constructed by adopting the invariant topological feature preserving algorithm.

3. The method according to claim 2, 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 to 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 DEM, the SSA parameters are optimized through the 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 component and the seasonal water body component are classified into stability grades and processed for spatial consistency to obtain a stability grade classification map, and the stable water body is determined based on it.

4. The method according to claim 2, 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 time-series water distribution information to calculate the apparent displacement field. Calculate the hydraulic gradient field based on DEM 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.

5. The method according to claim 2, characterized in that 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 fluidity and connectivity are defined to obtain prior constraints. 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 and extract the river network skeleton through the graph cut algorithm.

6. The method according to claim 1, characterized in that The steps to identify hydraulic facilities and analyze their connectivity include: Spatially associate the water conservancy project facility data with the 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 the water conservancy project operation rules is constructed to form an operation rule knowledge base; Combining the river network-hydraulic project coupling data and the operation rule knowledge base, the dynamic connectivity under different operating conditions is analyzed based on the Petri net construction to 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 indexes are calculated to form a regional connectivity characteristic map.

7. The method according to claim 1, characterized in that The steps for adaptive partitioning of hydraulically driven computational units include: Based on hydraulic characteristics and regional connectivity feature maps, 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 point 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 boundary. 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.

8. The method according to claim 5, characterized in that The steps to construct a non-local energy function include: Based on the enhanced hydraulic characteristic field, an energy function including data term, smoothing term and non-local term is constructed, wherein 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 weights are set according to the physical constraints so that the nonlocal interaction strength decays with distance but the decay speed 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.

9. The method according to claim 3, characterized in that The steps of signal reconstruction using the adaptive regional segmentation SSA method include: Based on the DEM-derived terrain zoning, the study area was divided into sub-areas; 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-areas are integrated to form complete stable water body components and seasonal water body components.

10. The method according to claim 7, characterized in that The steps of identifying the critical flow change point 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 candidate rheological point set for hydraulic significance, and 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.

11. The method according to claim 6, characterized in that The steps of 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 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; Adopt 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.

12. The method according to claim 1, characterized in that The method also includes constructing a multidimensional water exchange relationship between the calculation units based on the plain water network hydrological calculation units 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 expression for the exchange relationship; 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