River and lake water dynamic model reduced-order inversion initialization method for intelligent water conservancy

Abstract

The invention discloses a river and lake hydrodynamic model reduced-order inversion initialization method for intelligent water conservancy, and the method comprises the steps: dividing a computational domain into a plurality of sub-domains according to hydraulic characteristics, building a physical constraint equation represented by reduced-order description parameters, and determining an independent to-be-inverted parameter set through the analysis of the degree of freedom; estimating an initial state of a part of sub-domains by using observation information or a priori value at a starting moment; determining an inversion window according to the hydraulic disturbance propagation time; constructing an agent model, and performing joint inversion on the to-be-inverted parameters and the boundary process line parameters to minimize the cost function; and driving the proxy model to generate a global initial state by using the optimal parameter obtained by inversion, mapping the global initial state to a basic model grid, and taking the global initial state as an initial condition after free evolution and boundary smoothing processing. According to the method, actual measurement data are fused through order reduction-inversion, global initial field construction can be rapidly and accurately completed, and the requirement of an intelligent water conservancy system for cold start of a hydrodynamic model is met.

Description

Technical Field

[0001] This invention relates to hydrodynamic numerical simulation, and in particular to a method for order reduction and inversion initialization of river and lake hydrodynamic models for smart water conservancy. Background Technology

[0002] The intelligent water conservancy system, based on a digital twin watershed, requires hydrodynamic numerical models to start up and output reliable forecast results within a minute-level timescale. Model operation relies on an accurate initial field (the distribution of water level and velocity across the entire computational domain at startup); an inaccurate initial field will cause spurious fluctuations or even numerical divergence. In actual operation and maintenance, models often face cold starts due to system downtime, failures, or initial deployment, requiring rapid reconstruction of the entire domain's initial field based solely on limited observational data. For complex water bodies coupling plain river networks and lakes, the computational domain contains numerous one-dimensional river cross-sections and two-dimensional lake and beach grid cells, with state space dimensions reaching tens to hundreds of thousands. The construction of a cold-start initial field has become a key bottleneck restricting the continuity of intelligent water conservancy operations.

[0003] Currently, there are two main types of commonly used initialization methods. The first type is the long-term warm-up method, which starts with rough initial conditions and continuously drives the model operation with boundary conditions, relying on the dissipation mechanism of the hydrodynamic equations themselves to gradually reduce the initial error. This method often requires several hours or even longer warm-up times for plain river networks with lengths of tens to hundreds of kilometers, far exceeding the rapid start-up response requirements of smart water conservancy systems. At the same time, the warm-up method only relies on the indirect driving of boundary conditions and cannot actively integrate real-time observation information from stations within the computational domain. When there are uncertainties in the boundary conditions or complex hydraulic structure scheduling within the domain, the internal state after the warm-up may still deviate significantly from the actual hydrological conditions, greatly restricting the application of model results. The second type is the data assimilation method, including ensemble Kalman filtering and four-dimensional variational assimilation, which corrects the state estimate by incorporating observation data during model operation. However, data assimilation involves filtering or variational optimization within the complete state space of the basic model. The computational cost is closely related to the state dimension. Ensemble Kalman filtering requires maintaining a large set to characterize the error covariance, and four-dimensional variational methods require constructing an adjoint model and iterating multiple times. These methods are extremely difficult for the development and implementation of complex water areas coupled with river networks, lakes, and hydraulic structures, making it difficult to complete initialization within the time window required by smart water conservancy. The (basic) model used in existing technologies is a full-precision hydrodynamic numerical model for formal forecasting operations. It adopts the complete Saint-Venant equations or shallow water equations as the descriptive equations. Its definition file contains all the basic data required for model construction, such as computational grid, cross-sectional geometry, topographic elevation, roughness coefficient, hydraulic structure parameters, and boundary condition settings.

[0004] The common bottleneck of both methods lies in the fact that they operate within the full-dimensional state space of the basic hydrodynamic model. Current technology lacks a method for rapidly and accurately constructing a global initial field for complex river-lake water bodies. Summary of the Invention

[0005] Purpose of the invention: The purpose of this invention is to provide a method for order reduction and inversion initialization of river and lake hydrodynamic models for smart water conservancy that can quickly and accurately construct the initial field of a hydrodynamic simulation model of a subcritical flow river-lake composite water area, meeting the cold start requirements of smart water conservancy systems.

[0006] Technical Solution: The present invention provides a method for order reduction and inversion initialization of river and lake hydrodynamic models for smart water conservancy. The computational domain of the river and lake hydrodynamic model has a tree-like topology, and the Froude number of each open channel section is less than a set value under design conditions. The method includes the following steps:

[0007] S1. Divide the computational domain into multiple subdomains according to hydraulic characteristics, establish physical constraint equations characterized by reduced-order descriptive parameters for each subdomain, and determine the independent set of parameters to be inverted by performing degree-of-freedom analysis on the constraint equations.

[0008] S2. Using observation information or prior values ​​at the start time, estimate the initial state of some subdomains in the computational domain;

[0009] S3. Determine the time window for inversion based on the maximum propagation time of hydraulic disturbances in each subdomain to the nearest station or the open boundary with a measured process line;

[0010] S4. Construct a surrogate model with higher computational efficiency than the basic model. Perform joint inversion of the parameter set to be inverted and the boundary process line parameters within the time window to minimize the cost function. The independent parameters to be inverted when the termination condition is met are taken as the optimal parameters.

[0011] S5. Use the optimal parameters obtained by inversion to drive the surrogate model to generate the global initial state, and map it to the computation grid of the basic model. After the basic model undergoes free evolution and boundary condition smoothing at a set time, it is used as the initial condition of the basic model.

[0012] Step S1 compresses the high-dimensional state space of the river-lake composite water area (typically involving tens of thousands to hundreds of thousands of grid cells) into a small number of independent parameters to be inverted (usually only a few to dozens). As long as these independent parameters can be accurately calculated, the initial field of the entire area can be calculated, greatly reducing the computational burden of the initial field. Because some sub-domains (shoals) do not have independent parameters to be inverted, step S2 estimates the initial state of sub-domains such as shoals in advance, providing reasonable initial conditions for the first forward modeling of the surrogate model and avoiding numerical divergence or slow convergence due to incorrect initial wet / dryness judgments. Step S3 scientifically determines the length of historical data required for inversion based on the propagation time of hydraulic disturbances, avoiding computational waste caused by the use of irrelevant historical data and ensuring that the inversion window contains sufficient information to constrain the initial state. Step S4 employs a surrogate model for joint inversion. While ensuring computational efficiency, it actively integrates measured data to accurately solve for all independent parameters to be inverted. These parameters can then be used to solve for the remaining dependent parameters in the initial field, thus obtaining the global initial field. Compared to long-term preheating methods and data assimilation methods, this significantly reduces computational load and greatly improves efficiency. Step S5 substitutes the inverted parameter values ​​into the surrogate model to generate the global initial state, maps it back to the basic model mesh, and allows for free evolution. The dissipation mechanism of the basic model itself eliminates the initial imbalance between the surrogate model and the basic model caused by equation simplification and mesh differences. After free evolution, the boundary conditions are smoothed to achieve a smooth and shock-free transition from the initial state to the formal forecasting state. In summary, this method significantly reduces the computational load of initial field solution and can efficiently and accurately construct the global initial field by integrating measured data, effectively initializing the basic model so that the constructed global initial field can be used for subsequent forecasts.

[0013] Preferably, in step S1, the computational domain is divided into three subdomains: subcritical flow river section, lake, and beach. Among them, the lake subdomain must satisfy the requirement that the maximum propagation time of the shortest connection path along the water body between any two inlets and outlets does not exceed one-tenth of the time of change of the inflow characteristics. The beach is the connected area outside the main channel of the subcritical flow river section where the terrain elevation is not higher than the design flood level.

[0014] By explicitly dividing the computational domain into three subdomains with different hydraulic characteristics—subcritical flow sections, lakes, and floodplains—and imposing strict lumped parameterization conditions on the lake subdomain, we ensure that the lake surface can be approximated by a single water level parameter, thus minimizing the state space of the lake subdomain. Simultaneously, we define floodplains as connected areas outside the main channel of the river section where the terrain elevation does not exceed the design flood level. This definition not only delineates the spatial extent of the floodplain and its adjacency with the river section but also uses the design flood level as an elevation threshold, ensuring the floodplain's submergence under design conditions and its ability to switch between wet and dry states. This differentiated domain division method simplifies the description parameters of each subdomain while characterizing their physical constraints, resulting in a significantly smaller total number of description parameters in the entire computational domain compared to the original basic model. This lays the foundation for subsequent efficient and stable low-dimensional inversion.

[0015] Preferably, the reduced-order description parameters for the subcritical flow section are the instantaneous flow rate and water level at the downstream control section; the reduced-order description parameter for the lake sub-region is the lake surface water level; no reduced-order description parameters are introduced for the beach sub-region, and its state is determined by the water level of the adjacent sub-region.

[0016] To address the hydraulic characteristics of subcritical flow downstream controlling upstream, the downstream cross-section's flow rate and water level are used as reduced-order description parameters for the one-dimensional river segment. This allows for the unique determination of the water surface line for the entire river segment upstream using the energy equation, significantly reducing the number of description parameters for long river segments. For lakes and floodplains, based on the aforementioned domain conditions, a single water level is used as the description parameter, while no independent parameters are introduced. This reduces the number of parameters while ensuring the inherent physical constraints of each subdomain, thus guaranteeing the accuracy of subsequent inversion.

[0017] Preferably, the degree of freedom analysis in step S1 is as follows: establish a set of constraint equations between the reduced-order description parameters of each subdomain, normalize and scale the Jacobian matrix of the constraint equations and eliminate column pivots to determine the independent free variables as an independent set of parameters to be inverted; wherein, when choosing between parameters directly observed by the station and parameters not observed, the parameters directly observed by the station are selected as independent free variables.

[0018] By constructing a set of constraint equations among the parameters and performing numerical rank analysis, the number of truly independent free variables in the system can be accurately identified under the premise of satisfying all physical constraints (such as flow conservation and water level continuity). Prioritizing the retention of observable parameters as independent free variables allows subsequent inversion to utilize measured data more directly and effectively, thus simplifying the inversion problem into an optimization problem of solving a small number of key parameters, greatly improving computational efficiency and the stability of the inversion results.

[0019] As a preferred option, in step S3, the downstream propagation velocity of the subcritical flow section is calculated based on the motion wave velocity, and the upstream propagation velocity is calculated based on the inverse characteristic line velocity; the propagation time inside the lake is taken as zero; for each subdomain, the shortest time for the hydraulic disturbance inside to reach any station or the open boundary with a measured process line is calculated, and the maximum value of the whole domain is multiplied by a preset amplification factor as the inversion time window length, and the time window length does not exceed the preset maximum value.

[0020] To address the different propagation velocities of disturbances upstream and downstream in subcritical flows, the propagation time of hydraulic disturbances is estimated using both the motion wave velocity and the inverse characteristic line velocity, respectively. This allows for a more accurate determination of the propagation time. By calculating the time it takes for disturbances in each subdomain to reach the nearest effective monitoring target and using the maximum value across the entire domain to determine the inversion window, the length of the window is ensured to be sufficient for the hydraulic information of most key areas within the computational domain to propagate to the observation points. This guarantees the rationality of the independent inversion parameters obtained from the inversion.

[0021] Preferably, step S4 employs a two-layer solution strategy: the first layer uses the instantaneous observations at the start time to solve for the initial values ​​of the parameters to be inverted under the quasi-steady-state assumption; the second layer uses the solution from the first layer as the initial iteration starting point, and employs the Levenberg-Marquardt algorithm to jointly invert the parameters to be inverted and the boundary process line parameters within the time window, applying constraints such as non-negative flow rate and water level not lower than the corresponding subdomain topographic elevation during the iteration process.

[0022] The first layer utilizes instantaneous observations at startup to quickly provide an initial guess close to the optimal solution, preventing the second-layer time-varying inversion from getting trapped in local optima or iterative divergence. The second layer employs the Levenberg-Marquardt algorithm for joint inversion and ensures the physical authenticity of the inversion results by applying physical constraints, thus obtaining accurate optimal parameters while maintaining computational efficiency.

[0023] As a preferred option, the spatial grid of the surrogate model is independently sparsed relative to the basic model while maintaining topological consistency, and uses terrain data from the same source or with an elevation deviation not exceeding a preset threshold. The surrogate model uses the diffusion wave equation to describe the one-dimensional river segment and the two-dimensional beach, and uses the water balance ordinary differential equation to describe the lake. The flow velocity of the one-dimensional river segment and the two-dimensional beach is calculated by the Manning formula, and a minimum slope threshold is introduced into the Manning formula to avoid numerical singularities when the water surface tends to be horizontal.

[0024] The surrogate model significantly improves forward modeling speed by describing each subdomain using independent sparse spatial grids and computationally less computationally intensive diffusion wave and water balance equations, enabling efficient multiple iterative inversions within a time window. Simultaneously, sharing source topographic data with the basic model ensures consistency of key hydraulic features. Introducing a minimum slope threshold into the Manning formula resolves the numerical singularity problem caused by zero gradient when the water surface is horizontal, thus improving the numerical stability of the surrogate model.

[0025] Preferably, the expression for the cost function in step S4 is:

[0026]

[0027] in This refers to the serial number of the water level monitoring station. The time step number of the observed data exists. and These are the simulated water level and the measured water level at the k-th time step for the i-th water level station, respectively. This refers to the serial number of the flow measurement station. and The first Simulated and measured flow rates at the k-th time step of a flow meter station; For each parameter in the independent set of parameters to be inverted, For the j-th independent parameter to be inverted, Let be the prior value of the j-th independent parameter to be inverted. is the characteristic scale of the j-th independent parameter to be inverted; The open boundary number that needs to be inverted. For control point numbers, For the first The total number of control points for the open boundary; For the first The open boundary in the first Parametric values ​​at each control point; The weight of the fitting term for the i-th water level observation station is the reciprocal of the variance of the water level observation error at that station. For the first The weight of the fitted term for the flow observation station is the reciprocal of the variance of the flow observation error at that station; Let be the weight of the regularization term for the j-th independent parameter to be inverted, and let its value be such that the ratio of the regularization term to the sum of all observed fitted terms at the initial iteration is . ~ ; The weights for the boundary smoothing term are set such that the ratio of the smoothing term to the sum of all observed fitted terms is equal to the order of magnitude of the sum of the terms in the initial iteration. ~ .

[0028] This cost function incorporates water level and flow rate observation fitting terms to ensure that the inversion results fully utilize multi-source measured data. The parameter regularization term constrains the inversion parameters to avoid deviating too far from prior values, preventing overfitting to observational noise and thus improving the stability of the inversion results. The boundary condition time smoothing term ensures the temporal rationality of the inverted boundary process lines, avoiding drastic oscillations. Through adjustable weighting coefficients, this method can flexibly balance the importance of each term, thereby obtaining stable and reasonable inversion results.

[0029] As a preferred option, in step S2, the water level at the boundary of each beach area is estimated using water level observation or prior value at the start time, and compared with the lowest topographic elevation of the area. If the difference exceeds the tolerance threshold, it is determined to be wet; otherwise, it is determined to be dry. When a beach is adjacent to multiple water bodies, the highest water level is taken.

[0030] This prediction method is simple and efficient. It can quickly estimate the wet and dry state of the beach using limited information at the start time, providing key initial conditions for the first forward modeling calculation of the surrogate model.

[0031] Preferably, in step S5, the global initial state generated by the proxy model is mapped to the computational grid of the basic model according to the following rules: the water level of the subcritical flow section is directly mapped to the computational grid corresponding to the basic model; the lake water level is uniformly assigned to all computational grids within the lake area; the beach water level is mapped to the corresponding computational grid through interpolation; in the computational grid of the basic model, the water level of the part that exceeds the computational domain of the proxy model is taken as the water level of the nearest grid cell at the boundary of the computational domain of the proxy model; the value of the water level after mapping of each subdomain being lower than the terrain elevation is truncated to the terrain elevation; the cross-sectional average flow velocity of the subcritical flow section is mapped to the flow velocity component of the corresponding computational grid according to the river channel direction, and the initial flow velocity of the lake and beach is taken as zero;

[0032] The boundary condition for free evolution is that the value at the end of the inversion window remains unchanged. The termination sign is that the rate of change of water level in the whole area is lower than the preset threshold and the preset number of steps is continuously met after excluding units with water depth lower than the preset minimum water depth threshold, or the preset maximum time limit is reached.

[0033] After the free evolution ends, the boundary conditions smoothly transition to the real-time value according to the following cosine function within the preset transition period:

[0034]

[0035]

[0036] Where t is the current time, These are weighting coefficients. For transitional boundary conditions, These are the boundary condition values ​​at the end of the inversion window. These are real-time, time-varying boundary condition values. and These are the start and end times of the preset transition period, respectively.

[0037] Through clear and explicit mapping rules, the initial state generated by the surrogate model is accurately and completely converted into the initial field required by the basic model, ensuring the rationality of the spatial distribution. During the free evolution stage, fixed boundary conditions are used, and the dissipation mechanism of the basic model itself is utilized to eliminate the initial imbalances introduced by equation and mesh differences, allowing the model to quickly enter a stable state. Finally, a smooth transition in the form of a cosine function is employed to avoid numerical shocks caused by abrupt changes in boundary conditions, ensuring the stability of the model's subsequent operation.

[0038] Beneficial effects: By compressing the high-dimensional state space into low-dimensional independent parameters through a domain-based order reduction strategy, the computational complexity of the initialization problem is fundamentally reduced. By constructing a lightweight surrogate model and performing a two-layer joint inversion within the optimal time window, multi-source observation data are efficiently and accurately fused. Finally, through a reliable mapping, free evolution, and smooth transition process, a stable and accurate initial field is quickly generated for the basic model. This effectively solves key problems such as slow initial field construction speed, low accuracy, and difficulty in fusing observation data in cold-start scenarios for river-lake composite water bodies. Attached Figure Description

[0039] Figure 1 This is a flowchart of the method of the present invention;

[0040] Figure 2 This is a simplified topological diagram of rivers and lakes in an embodiment of the present invention;

[0041] Figure 3 This is an application effect diagram of an embodiment of the present invention. Detailed Implementation

[0042] As shown in the figure, the method for order reduction and inversion initialization of river and lake hydrodynamic models for smart water conservancy described in this invention has a tree-like topology in the computational domain of the river and lake hydrodynamic model. Under design conditions, the Froude number of each open channel section is less than a set value (usually 1.0). The method includes the following steps S1 to S5:

[0043] Step S1: Domain Reduction and Degrees of Freedom Analysis

[0044] The computational domain is divided into multiple subdomains based on hydraulic characteristics. Specifically, the computational domain is divided into three subdomains: one-dimensional subcritical flow river sections, lakes, and floodplains.

[0045] A one-dimensional subcritical flow section refers to an open channel section where the Froude number is less than 1.0 at all calculated cross-sections under design flow conditions.

[0046] A lake (including a reservoir) refers to a shoreline-enclosed body of water that meets the lumped parameterization condition. Under design flood levels, this enclosed body of water has no direct openings to external water bodies except for the inflow and outflow sections. The lumped parameterization condition is: the maximum propagation time along the shortest path between any two inflow and outflow points does not exceed one-tenth of the time of change in the inflow characteristic. The propagation time is calculated based on the shallow water gravity wave velocity. Conservative estimate The minimum still water depth along the propagation path under design conditions is given by , and g is the acceleration due to gravity. The time of change in the inflow characteristic is taken as the minimum time required for the flow in each inflow segment to rise from base flow to peak flow under design conditions. The above lumped parameterization conditions are determined once based on design conditions and serve as a static basis for subdomain classification, and do not change with actual operating conditions. Large water bodies that do not meet these conditions are not suitable for the lake lumped parameterization of this invention and should be excluded from the computational domain after setting open boundaries at their inflow and outflow gates.

[0047] The floodplain subdomain is defined as the connected area outside the main channel of the one-dimensional river segment, where the terrain elevation does not exceed the design flood level. The main channel extent is defined as the lateral range covered by the one-dimensional river segment calculation cross-section in the basic model definition file. Connectivity is determined according to the 8-directional connectivity criterion of the DEM raster (i.e., each cell in the DEM raster is considered connected to its 8 surrounding adjacent cells). The inundation state of the two-dimensional floodplain depends on the water level of the adjacent water body and does not introduce independent descriptive parameters.

[0048] Physical constraint equations characterized by reduced-order descriptive parameters are established for each subdomain.

[0049] For a one-dimensional subcritical flow river segment, under the assumptions of gradually varied flow and quasi-steady state, the instantaneous flow rate at the downstream control section of the river segment is used as the criterion. and downstream control section water level For describing parameters. The water surface line in the subcritical flow section is controlled by the downstream, given... and The water surface elevation and velocity at each cross-section can then be determined using the friction loss energy equation. The discretized form of the friction loss energy equation represents the energy conservation relationship between adjacent cross-sections, specifically:

[0050]

[0051] Among them, subscript and The upstream and downstream sections are numbered respectively. For terrain elevation, Because of the water depth, The cross-sectional average velocity is... This is the kinetic energy correction factor. It is the acceleration due to gravity. The average frictional slope between the two cross sections. The above water level lines are used to calculate the water level at the downstream control sections. Starting from the first point, the solution is obtained by proceeding upstream section by section.

[0052] Friction gradient is calculated using the Manning formula:

[0053]

[0054] in The roughness coefficient is Manning's coefficient. The hydraulic radius. The average frictional gradient between the two sections. The roughness coefficient is determined by averaging the frictional slopes of the upstream and downstream sections. It is taken from the basic model definition file.

[0055] At the inflow or outflow points and at the locations of hydraulic structures, nodes are inserted to divide the river segment into several sub-segments. The flow rate in each sub-segment is constant along the course and does not contain abruptly changing flow sections. The water surface line of the entire sub-segment can be controlled by the flow rate at the downstream control section. and water level The flow at the hydraulic structure itself is uniquely determined. Rapidly changing flows at the structure are handled by the flow formula, and the parameterized regions of each sub-section after decomposition do not include the rapidly changing flow sections of the structure itself. When the inflow or outflow originates outside the computational domain, its flow rate is taken from measured data or the basic model definition file, serving as a known external input; when it originates from other subdomains within the computational domain, it is reflected through the flow conservation constraints after decomposing the sub-sections.

[0056] The reduced-order descriptive parameter for the lake subdomain is the lake surface water level. The beach subdomain does not introduce a reduced-order description parameter; its state is determined by the water level of the adjacent subdomain.

[0057] Establish a set of constraint equations among the reduced-order description parameters of each subdomain, specifically including:

[0058] The first category concerns the flow relationships of hydraulic structures: formulas for weir flow and gate outflow establish functional relationships between upstream and downstream water levels and the flow rate. Each structure provides corresponding constraint equations (i.e., formulas for weir flow and gate outflow) under given operating conditions. Fixed geometric parameters and gate opening sequences are taken from the basic model definition file. The flow rate of the pumping station is directly given by the unit's operating status and treated as a known flow rate boundary condition.

[0059] The second type is the continuous water level condition at the coupling interface: when the river flows into the lake, the water level at the downstream control section is equal to the lake surface water level, that is... When a lake flows out, the water level at the upstream section of the outflowing river segment... satisfy ,and The flow rate can be controlled by the downstream section of this river. and water level The water level is uniquely determined by calculating cross-section by cross-section from downstream to upstream using the friction energy equation, and the water level remains continuous at the coupling interface. This equation, together with the energy equation along the river, constitutes the constraint equations; the water level in each river segment at the confluence node remains consistent. When the downstream control section of a river segment is located precisely at the open boundary of the computational domain, the water level at that section is included in the open boundary hydrograph parameter processing. For open boundaries with measured hydrographs, the water level at that section is taken as the measured value, which is a known quantity; for open boundaries without measured hydrographs, the water level at that section is used as the boundary control variable to be inverted in step S4, and the current iteration guess value is taken when deriving the dependent parameters in the constraint equations here. In all constraint equations involving open boundary water level parameters, the parameter is replaced with the above-mentioned known value or the current iteration guess value before participating in the constraint solution; if all unknown parameters in the equation after replacement have been determined, then the equation degenerates into an identity, and the equation is removed from the effective constraint equation set.

[0060] The third type is the continuous water level condition at the river segment split point: when a river segment is split into sub-segments due to lateral inflow, the water level at the downstream end of the upstream sub-segment at the split point is equal to the water level at the upstream end of the downstream sub-segment. If there is a head difference between upstream and downstream of a hydraulic structure, the continuous water level condition does not hold, and the third type of constraint does not apply to such split points; the third type of constraint only exists independently at split points with lateral inflow where there are no hydraulic structures.

[0061] The fourth type is the flow conservation condition: at the confluence node. ,in and These represent the flows at the inflow and outflow nodes, respectively. The flow conservation form at the split point is: the upstream sub-segment flow plus the net inflow from the lateral side at that point equals the downstream sub-segment flow, i.e. ,in, and These are the flow rates at the upstream and downstream control sections, respectively. This represents the net inflow from the side, with positive values ​​for inflow and negative values ​​for outflow; when there is no side inflow... At the upstream open boundary of the measured flow process line, the flow in that river segment equals the measured inflow at the corresponding time, constituting an additional constraint equation. The water balance constraint for each lake adopts... The steady-state form, where the net inflow source term of the lake surface is neglected. The unsteady-state effects of lake storage and The effects are handled by the water balance ordinary differential equation of the surrogate model in step S4. The initial deviation of the steady-state constraint is corrected by the time-varying inversion in step S4 as the initial guess value.

[0062] The above constraint equations are used to derive all dependent parameters at a given time given the values ​​of independent free variables, i.e., to construct an initial snapshot of a time section. This constraint system is applicable to any given time, and in step S4, the first layer is... The same applies when establishing parameter mappings. During the surrogate model integration in step S4, each parameter evolves independently by the describing equations, and the constraint equations in the steady-state form described above are no longer applied repeatedly.

[0063] Assume the entire domain contains The river section (including the sub-sections after division) and If there are lakes, then the total number of parameters in the subdomain is The water level at the downstream control section located at the open boundary is not included in the internal parameters but is classified as the boundary control variable in step S4; the flow parameter of the river segment located at the open boundary is still retained as an internal parameter. Let the total number of internal parameters after deducting the open boundary water level parameter be... Numerical rank analysis was performed on all valid constraint equations after removing the degenerate equations to determine the rank of their Jacobian matrix. The effective degrees of freedom of the system are The above. The internal parameters constitute the total parameter vector. ,in The parameters are independent free variables (i.e., independent parameters to be inverted), forming an independent free variable vector. ,the remaining The dependent variables are selected as follows: For all valid constraint equations, construct a Jacobian matrix at the prior values. First, normalize and scale each column according to the characteristic scale of the corresponding parameter (the characteristic scale for water level parameters is 1m, and the characteristic scale for flow parameters is the design flow of the corresponding river section). Then, perform column pivot elimination, selecting the parameters corresponding to the pivots as dependent variables, and the rest as independent free variables. Parameters directly observed by stations are preferentially retained as independent free variables without affecting the numerical stability of the elimination process.

[0064] Given After determining the values ​​of the independent free variables, the above nonlinear constraint equations need to be solved to derive all dependent parameters. The Newton-Raphson method is used for the solution. Initial guesses are taken from prior values, and the convergence criterion is that the residuals of each equation divided by the characteristic scale are less than [a certain value]. The water level residual is divided by 1m, and the flow residual is divided by the design flow rate of the corresponding river section. The maximum number of iterations does not exceed 30. If the process fails to converge, the iteration value with the smallest current residual is taken as the initial guess value for step S4, and its accuracy is corrected by the inversion in step S4.

[0065] Step S2: Initial State Estimation

[0066] Using observational information or prior values ​​at the start time, the initial inundation state of each beach area in the computational domain is estimated, providing the initial wet and dry conditions of the beach for the first forward modeling step S4. This estimation is based on... Based on the readily available observed water level, the error is automatically corrected by the wet / dry discrimination algorithm of the surrogate model in step S4 during the integration process.

[0067] Specifically, using Water levels at the boundaries of each floodplain area are estimated using either real-time water level observations or prior values. These levels are then compared with the lowest topographic elevation of the area. If the difference exceeds a tolerance threshold (e.g., 0.2m to 0.3m, but not less than the vertical accuracy of the DEM), the area is considered wet; otherwise, it is considered dry. The water level at the boundary of a floodplain area is calculated using descriptive parameters from the adjacent river segment or lake sub-region. When a floodplain is adjacent to multiple water bodies, the highest water level among them is used for determination.

[0068] Step S3: Determine the inversion time window

[0069] The time window for inversion is determined based on the maximum propagation time of hydraulic disturbances within each subdomain to the nearest station or the open boundary with a measured process line.

[0070] The downstream propagation velocity of a one-dimensional subcritical flow segment is estimated based on the motion wave velocity, which is determined by the channel cross-sectional geometry and hydraulic parameters. Calculate, where, For the speed of motion, The width of the water surface. For cross-sectional flow, The cross-sectional water level The velocity is obtained from Manning's formula and cross-sectional geometry; the upstream propagation velocity is estimated using the reverse characteristic line velocity, which equals the shallow water gravity wave velocity minus the average cross-sectional velocity, i.e. The cross-sectional parameters and flow velocities required for estimating propagation speed are preferentially taken from the previous model run archive; if unavailable, they are taken from the cross-sectional geometric data corresponding to historical average water levels and flow rates. In the tree-like topology, each node has a unique path to the nearest effective monitoring target (station or open boundary with a measured process line). When the path passes through a hydraulic structure, connectivity in that direction is determined based on the disturbance propagation characteristics of that structure under its current operating conditions; a weir in a free-flowing state blocks propagation in the upstream direction, and the propagation time in that direction is taken as infinite. Open boundaries without a measured process line are not considered as targets for propagation time.

[0071] Lake subdomains that meet the lumped parameterization conditions are treated as single nodes, with their internal propagation time set to zero. Beach subdomains do not introduce independent description parameters; their propagation time requirements are automatically included in adjacent river segments or lake subdomains.

[0072] For each sub-river segment and each lake sub-region, calculate the shortest time (not less than zero) for hydraulic disturbances within the sub-region to reach any monitoring station or the open boundary of a measured process line. Multiply the maximum value across the entire region by a preset amplification factor to obtain the inversion time window length T. Amplification factor of 1.2 to 1.5 is recommended to compensate for wave velocity estimation bias. The time window length should not exceed the preset maximum value (e.g., 12 hours). Sub-regions where disturbances cannot propagate to any effective monitoring target within T directly use prior values ​​as the initial state and do not participate in the inversion.

[0073] Step S4: Joint Inversion

[0074] A surrogate model with higher computational efficiency than the basic model is constructed. The surrogate model maintains consistency with the basic model in terms of water system topology and coupling interface type. The surrogate model uses the diffuse wave equation as a simplified descriptive equation, preserving mass conservation among subdomains. The surrogate model uses terrain data from the same source or with elevation deviations not exceeding a preset threshold, similar to the basic model. The spatial grid of the surrogate model can be the same as the basic model or independently sparsified relative to the basic model while maintaining topological consistency. The geometries of the sparsified sections should retain the main characteristics of the water level-area and water level-wetted perimeter relationships. The beach area is... The initial water depth at any given time is determined by the difference between the water level in the neighboring subdomain, which is determined by the independent free variables through constraint equations, and the elevation of the beach topography. When the difference is negative, the water depth is set to zero. During the inversion process, as the vector of independent free variables... The initial water depth of the beach is recalculated during iterative updates, overriding the prediction results from step S2. Areas identified as dry in step S2 are initialized to zero water depth in the first forward modeling iteration and participate in the calculation within the surrogate model. Their wet / dry status is dynamically updated by the model's wet / dry discrimination algorithm. The output time step of the surrogate model should not exceed the time interval of the observation data; the simulated values ​​are linearly interpolated to the observation time for alignment.

[0075] The surrogate model uses the diffusion wave equation to describe the one-dimensional river section and the two-dimensional floodplain, and the water balance ordinary differential equation to describe the lake. The flow velocity in the one-dimensional river section and the two-dimensional floodplain is calculated by the Manning formula, and a minimum slope threshold is introduced into the Manning formula to avoid numerical singularities when the water surface tends to be horizontal.

[0076] The specific describing equation is as follows:

[0077] The subcritical flow section is described by a one-dimensional diffusion wave equation, specifically:

[0078]

[0079]

[0080] in, The width of the water surface. Because of the water depth, For cross-sectional flow, The lateral inflow per unit river length This refers to the terrain elevation; The water surface gradient is defined as the negative gradient of the water level along the course of the river, x is the longitudinal coordinate of the one-dimensional subcritical flow segment along the main flow direction, and t is time.

[0081] Calculate according to the following formula

[0082]

[0083] in, V is the cross-sectional area of ​​the water passage; V is the average flow velocity of the cross-section.

[0084] V is calculated according to the following formula.

[0085]

[0086] in, The roughness coefficient is Manning's coefficient. For hydraulic radius, To avoid numerical singularities when the water surface tends to be horizontal, it is recommended to take [value missing]. ;

[0087] Lakes are described by water balance ordinary differential equations, specifically:

[0088]

[0089] in, For the lake at water level The water surface area at that location; The equivalent flow rate for net inflow to the lake surface; and These represent the flow rates at the inflow and outflow nodes of the lake, respectively.

[0090] The beach is described by a two-dimensional diffusion wave equation, specifically:

[0091]

[0092]

[0093] in, and They are respectively and The depth-average flow velocity in the direction of travel. and They are respectively and The water surface gradient in the direction, For surface source aggregation;

[0094] and The calculation formula is

[0095]

[0096] The above formula uses water depth Approximate hydraulic radius .

[0097] The control variables used in the surrogate model for inversion consist of two parts. The first part is the backtracking start time. independent free variable vector ,Include The parameters are 1,000 independent free variables, and their prior values ​​are selected according to the following priority: the archived state from the previous model run takes precedence, followed by the estimated value derived from the constraint equations in step S1 based on the measured values ​​at the stations, and finally the statistical mean of the historical observations of the parameter at the stations. The second part is the parameterized representation of the water level or flow rate process lines of each open boundary that needs to be inverted within the time window. Piecewise linear functions are used, and each open boundary is described by equally spaced control points, with the time interval between control points being [missing information]. to The open boundary at the upstream inlet is typically defined by a flow rate hydrograph, while the open boundary at the downstream outlet is typically defined by a water level hydrograph. Open boundary values ​​with measured hydrographs are fixed and not used as control variables. The first control point of the open boundary hydrograph to be inverted corresponds to... At each moment, its value is determined by the current independent free variable in each iteration. The full parameter state determined by the constraint equations is given and is not used as an independent control variable; the initial guess values ​​of the other control points are taken from the prior values, and if there is no prior information, the open boundary parameter values ​​corresponding to the solution of the first layer are taken.

[0098] During inversion, the subcritical flow segments and lakes are advanced using the same principal time step; the time step for the floodplain is determined by the stability conditions of the explicit scheme, and a sub-step approach is used within each principal time step if necessary. The lake water balance ordinary differential equation is explicitly coupled alternately with the adjacent flow segments. The coupling between the floodplain and adjacent flow segments or lakes is achieved by applying the water level of the adjacent water body as a boundary condition at the floodplain boundary. The inflow and outflow of the floodplain are cumulatively fed back to the lateral inflow terms of the adjacent flow segments or the water balance equation of the lake at the end of each principal time step.

[0099] The inversion step in this process employs a two-layer solution strategy. The first layer utilizes... The instantaneous water level or flow rate observations at each station at any given time do not involve time integration. Under the quasi-steady-state assumption, the constraint equations from step S1 are used in... Establish full parameter vector at all times The mapping between observations and independent variables; when the number of observations equals the number of independent variables, the solution is obtained directly; when the number of observations is greater than the number of independent variables, the objective is to minimize the weighted sum of squared residuals, with the weights being the reciprocals of the variances of the observation errors; when the number of observations is less than the number of independent variables, parameters lacking observation constraints are directly taken from prior values. The first-level solution is used as the initial guess. The first level in... The solution at time t is simultaneously used as vector of independent free variables at time points The initial guess value. The second layer uses the solution from the first layer as the initial iteration starting point, within the time window. The independent free variable vector and boundary process line parameters are jointly inverted. The Levenberg-Marquardt algorithm is used to minimize the cost function, and the Jacobian matrix is ​​approximated by forward numerical difference. Upper and lower bound constraints of the physical feasible region (non-negative flow rate and water level not lower than the corresponding subdomain topographic elevation) are imposed on the control variables during iterations, with out-of-bounds values ​​truncated to the constraint boundaries. The termination condition is that the relative decrease in the cost function compared to the previous iteration is less than... The maximum number of iterations does not exceed 50. The vector of independent free variables when the termination condition is met. As the optimal parameter.

[0100] The expression for the cost function is:

[0101]

[0102] in This refers to the serial number of the water level monitoring station. The time step number of the observed data exists. and These are the simulated water level and the measured water level at the k-th time step for the i-th water level station, respectively. This refers to the serial number of the flow measurement station. and The first Simulated and measured flow rates at the k-th time step of a flow meter station; For each parameter in the independent set of parameters to be inverted, For the j-th independent parameter to be inverted, Let be the prior value of the j-th independent parameter to be inverted. is the characteristic scale of the j-th independent parameter to be inverted; The open boundary number that needs to be inverted. For control point numbers, For the first The total number of control points for the open boundary; For the first The open boundary in the first Parametric values ​​at each control point; The weight of the fitting term for the i-th water level observation station is the reciprocal of the variance of the water level observation error at that station. For the first The weight of the fitted term for the flow observation station is the reciprocal of the variance of the flow observation error at that station; Let be the weight of the regularization term for the j-th independent parameter to be inverted, and let its value be such that the ratio of the regularization term to the sum of all observed fitted terms at the initial iteration is . ~ ; The weights for the boundary smoothing term are set such that the ratio of the smoothing term to the sum of all observed fitted terms is equal to the order of magnitude of the sum of the terms in the initial iteration. ~ .

[0103] Step S5: Initial Field Generation and Mapping

[0104] Drive the surrogate model using the optimal parameters obtained from the inversion. Points to ,extract The global water level and velocity field at time t is mapped onto the computational grid of the basic model. After free evolution and boundary condition smoothing by the basic model over a set time, these conditions are used as the initial conditions of the basic model.

[0105] The mapping rules are as follows: the water level of a one-dimensional subcritical flow segment is directly mapped to the computational grid corresponding to the basic model; the lake water level is uniformly assigned to all computational grids within the lake area; the floodplain water level is mapped to the corresponding computational grid through interpolation; for the portion of the computational grid of the basic model that extends beyond the computational domain of the surrogate model, the water level is taken as the water level of the nearest grid cell at the boundary of the surrogate model's computational domain; values ​​where the mapped water level of each subdomain is lower than the terrain elevation are truncated to the terrain elevation. The cross-sectional average velocity of the one-dimensional subcritical flow segment is mapped to the velocity component of the corresponding computational grid according to the river channel direction, and the initial velocity of the lake and floodplain mappings is zero.

[0106] The basic model is started with the initial field obtained from the mapping and undergoes free evolution. The boundary conditions for free evolution remain unchanged at the end time of the inversion window. The initial imbalance between the surrogate model and the basic model caused by equation simplification and mesh differences is eliminated by the dissipation mechanism of the basic model itself. The termination criteria for free evolution are: after excluding cells with water depths lower than the preset minimum water depth threshold, the rate of change of water level in the entire domain is lower than the preset threshold and continuously meets the preset number of steps, or the preset maximum time limit is reached. If the termination condition is not met even after reaching the maximum time limit, the state at that time is used as the initial condition to continue starting the basic model.

[0107] After the free evolution ends, the boundary conditions smoothly transition to the real-time value according to the following cosine function within the preset transition period:

[0108]

[0109]

[0110] Where t is the current time, These are weighting coefficients. For transitional boundary conditions, These are the boundary condition values ​​at the end of the inversion window. These are real-time, time-varying boundary condition values. and These are the start and end times of the preset transition period. A transition period of 5-20 minutes is recommended, after which the system will transition to normal operation.

[0111] To better illustrate this method, a specific example will be used below for further explanation.

[0112] Taking a typical plain river and lake system as an example, the complete execution process of the method of the present invention is explained.

[0113] Geometric description of the river system. The system comprises a central lake (approximately 5 km² in area, with an average depth of about 3 m), three inflow sections (denoted as R1, R2, and R3), and one outflow section (R4). R1 ​​is 20 km long, R2 is 15 km long, R3 is 10 km long, and R4 is 25 km long. Each section has an approximate trapezoidal cross-section, with a riverbed slope of approximately... The Manning roughness coefficient is 0.025~0.035. There is a beach area on the south side of the lake. A full-section overflow weir is located in the middle of the R2 river segment, dividing R2 into two sub-segments: R2a (upstream of the overflow weir, approximately 7.5 km long) and R2b (downstream of the overflow weir, approximately 7.5 km long). The downstream of R4 is the computational domain outlet (open boundary), with no measured flow rate profile. The upstream ends of R1, R2, and R3 are all computational domain inlets (open boundaries). Measured flow rate profiles are available at the upstream ends of R1 and R2, but not at the upstream end of R3. The station distribution is as follows: a flow station (W1, providing both flow rate and water level data) is located at the lake inlet downstream of R1; a water level station (W3) is located upstream of R3; and a water level station (WL) is located in the center of the lake, for a total of three stations. In this embodiment, the flow direction in each river segment does not reverse, and the lake satisfies the lumped parameterization conditions.

[0114] Step S1 is executed. All four river segments (R2 is divided into R2a and R2b via an overflow weir, totaling five sub-segments) satisfy the condition that the Froude number is less than 1.0. The downstream control section of R4 is located at the open boundary, and its water level... The flow rate of R4 is not included in the internal parameters and is given by the open boundary process line. Retain as an internal parameter.

[0115] The parameters for all sub-river sections and lakes are shown in Table 1:

[0116] Table 1. Parameters for each sub-river section and lake

[0117]

[0118] The total number of parameters for 5 sub-river sections and 1 lake is One; of which R4 The internal parameters are not included in the open boundary classification; the actual internal parameters are... indivual.

[0119] The constraint equations are listed one by one below, among which This represents the overflow weir flow relationship function, used to characterize the functional relationship between the upstream and downstream water levels and the flow rate over the weir. The specific form is selected according to the weir type and the outflow state, using the appropriate weir flow formula.

[0120] 1) R1 inflow water level is continuous (Type II constraint):

[0121] 2) Continuous R2b inflow level (Type II constraint):

[0122] 3) Continuous inflow water level of R3 (Type II constraint):

[0123] 4) R4 outflow water level is continuous (Type II constraint): (in (Current guess value for open boundary)

[0124] 5) R2 overflow weir flow relationship (first type of constraint):

[0125] 6) Flow rate conservation at the R2 overflow weir (Type IV constraint):

[0126] 7) Lake flow conservation (Fourth type constraint):

[0127] 8) Upstream open boundary flow constraint of R1 (fourth type of constraint): (in (For the upstream open boundary flow of R1)

[0128] 9) Upstream open boundary flow constraint of R2a (fourth type of constraint): (in (For the upstream of R2, open boundary flow)

[0129] In this embodiment, there is a head difference between the upstream and downstream of the overflow weir, so the third type of constraint is not applicable here. The above 9 constraint equations have a rank of 9 as verified by the Jacobian matrix, and the effective degrees of freedom are 9. .

[0130] Independent free variables were selected as Given Then, all dependent parameters can be determined sequentially using the above constraint equations.

[0131] Step S2 is executed. The area of ​​the beach on the south side of the lake is approximately 2 km², with the lowest topographic elevation being 12.3 m. The lake's water level was 12.8m at the time of observation, which is 0.5m different from the lowest topographic elevation. This exceeds the tolerance threshold and is therefore considered wet.

[0132] Step S3 is executed. In this embodiment, the overflow weir operates in a submerged outflow state, and hydraulic disturbances can cross the weir body in both directions. Effective monitoring targets are: station W1, upstream open boundary of R1, upstream open boundary of R2, station WL, and station W3; the downstream open boundary of R4 and the upstream open boundary of R3 have no measured process lines and are not considered effective monitoring targets. Through segment-by-segment estimation, the sub-segment and sub-domain with the longest effective propagation time is R4 (approximately 1.8 hours), multiplied by an amplification factor of 1.4, yields... Hour.

[0133] Step S4 is executed. Control variables include: One independent free variable at time 1 Five independent control points were used for the open boundary water level process of R4 downstream and five independent control points for the open boundary flow process of R3 upstream (each open boundary contains six equally spaced control points; the first control point is determined by the constraint equation under the current iteration and is not considered an independent control variable), totaling 11 independent control variables. The time interval between control points was set to 0.5 hours.

[0134] The surrogate model adopts the same diffusion wave equation as the technical solution, and the spatial grid and time step are set according to stability requirements. Prior values ​​for lake water level: m. Actual flow rate measured at each boundary at different times: m³ / s, m³ / s.

[0135] First layer of utilization Instantaneous observations from the stations at specific times (water level W1 12.80m and flow rate 38m³ / s, water level W3 13.1m, water level WL 12.8m) are used to solve for the first layer of solutions under the quasi-steady-state assumption, with the objective of minimizing the weighted sum of squared residuals: m. All dependent parameters are determined by the constraint equations. The first-level solution also serves as... Initial guess value at time.

[0136] The second layer starts with the solution from the first layer, and uses the forward modeling results within a 2.5-hour window of the surrogate model to jointly invert with the continuous observation data from the station, and the optimal solution after convergence is: m. Root mean square error within the inversion window for each station: flow rate at station W1 is approximately 1.0 m³ / s, water level at station W3 is approximately 0.08 m, and water level at station WL is approximately 0.03 m.

[0137] Step S5 is executed. The surrogate model is integrated for 2.5 hours with optimal control variables, and the results are extracted. The global water level and velocity field at time t is used to generate the initial conditions for the basic model mesh according to the mapping rules. The basic model is then allowed to evolve freely, and the termination condition is met after about 200 seconds. Subsequently, the boundary conditions are smoothly transitioned to real-time time-varying values ​​within a 7-minute transition period, and the initialization is complete.

[0138] Comparison with the preheating method. To verify the effectiveness of the method of this invention, a long-term preheating method was used as a comparison. The preheating method uses prior values ​​as initial conditions; for open boundaries with measured process lines, measured values ​​are applied; for open boundaries without measured process lines, prior values ​​are taken, driving the basic model from... Run to The preheating method does not fuse observation data within the inversion window. At the end of the same 2.5-hour window, the root mean square error of the preheating method at each station is: flow rate at station W1 approximately 3.5 m³ / s, water level at station W3 approximately 0.21 m, and water level at station WL approximately 0.05 m, all greater than the corresponding values ​​of the method in this invention. This is because the preheating method relies solely on boundary condition driving and model dissipation to eliminate initial biases, and the initial error has not been sufficiently decayed within the limited preheating time; while this invention directly fuses observation data through order reduction-inversion to determine the initial state, and the accuracy of the initial field is not limited by the preheating time.

Claims

1. A method for order reduction and inversion initialization of a river and lake hydrodynamic model for smart water conservancy, wherein the computational domain of the river and lake hydrodynamic model has a tree-like topology, and the Froude number of each open channel section is less than a set value under design conditions, characterized in that... Includes the following steps: S1. Divide the computational domain into multiple subdomains according to hydraulic characteristics, establish physical constraint equations characterized by reduced-order descriptive parameters for each subdomain, and determine the independent set of parameters to be inverted by performing degree-of-freedom analysis on the constraint equations. S2. Using observation information or prior values ​​at the start time, estimate the initial state of some subdomains in the computational domain; S3. Determine the time window for inversion based on the maximum propagation time of hydraulic disturbances in each subdomain to the nearest station or the open boundary with a measured process line; S4. Construct a surrogate model with higher computational efficiency than the basic model. Perform joint inversion of the parameter set to be inverted and the boundary process line parameters within the time window to minimize the cost function. The independent parameters to be inverted when the termination condition is met are taken as the optimal parameters. S5. Use the optimal parameters obtained by inversion to drive the surrogate model to generate the global initial state, and map it to the computation grid of the basic model. After the basic model undergoes free evolution and boundary condition smoothing at a set time, it is used as the initial condition of the basic model.

2. The method according to claim 1, characterized in that: In step S1, the computational domain is divided into three subdomains: subcritical flow river section, lake, and beach. Among them, the lake subdomain must satisfy the requirement that the maximum propagation time of the shortest connection path along the water body between any two inlets and outlets does not exceed one-tenth of the time of change of the inflow characteristics. The beach refers to the connected area outside the main channel of the subcritical flow river section where the terrain elevation is not higher than the design flood level.

3. The method according to claim 2, characterized in that: The reduced-order descriptive parameters for the subcritical flow section are the instantaneous flow rate and water level at the downstream control section; the reduced-order descriptive parameter for the lake sub-region is the lake surface water level; no reduced-order descriptive parameters are introduced for the beach sub-region, and its state is determined by the water level of the adjacent sub-region.

4. The method according to claim 1, characterized in that: The degree of freedom analysis in step S1 is as follows: establish a set of constraint equations between the reduced-order description parameters of each subdomain, normalize and scale the Jacobian matrix of the constraint equations and eliminate column pivots to determine the independent free variables as an independent set of parameters to be inverted; wherein, when choosing between parameters directly observed by the station and parameters without observation, the parameters directly observed by the station are selected as independent free variables.

5. The method according to claim 2, characterized in that: In step S3, the downstream propagation velocity of the subcritical flow section is calculated based on the motion wave velocity, and the upstream propagation velocity is calculated based on the inverse characteristic line velocity; the propagation time inside the lake is taken as zero; for each subdomain, the shortest time for the hydraulic disturbance inside to reach any station or the open boundary with a measured process line is calculated, and the maximum value of the whole domain is multiplied by the preset amplification factor as the inversion time window length, and the time window length does not exceed the preset maximum value.

6. The method according to claim 1, characterized in that, Step S4 employs a two-layer solution strategy: the first layer utilizes the instantaneous observations at the start time to solve for the initial values ​​of the parameters to be inverted under the quasi-steady-state assumption; The second layer uses the solution from the first layer as the initial iteration starting point. Within the time window, the Levenberg-Marquardt algorithm is used to jointly invert the parameters to be inverted and the boundary process line parameters. During the iteration process, constraints are applied to ensure that the flow rate is non-negative and the water level is not lower than the corresponding subdomain topographic elevation.

7. The method according to claim 2, characterized in that: The spatial grid of the surrogate model is independently sparsed relative to the basic model while maintaining topological consistency, and uses terrain data from the same source or with an elevation deviation not exceeding a preset threshold. The surrogate model uses the diffusion wave equation to describe the one-dimensional river segment and the two-dimensional beach, and the water balance ordinary differential equation to describe the lake. The flow velocity of the one-dimensional river segment and the two-dimensional beach is calculated by the Manning formula, and a minimum slope threshold is introduced into the Manning formula to avoid numerical singularities when the water surface tends to be horizontal.

8. The method according to claim 1, characterized in that: The expression for the cost function in step S4 is: ， in This refers to the serial number of the water level monitoring station. The time step number of the observed data exists. and These are the simulated water level and the measured water level at the k-th time step for the i-th water level station, respectively. This refers to the serial number of the flow measurement station. and The first Simulated and measured flow rates at the k-th time step of a flow meter station; For each parameter in the independent set of parameters to be inverted, For the j-th independent parameter to be inverted, Let be the prior value of the j-th independent parameter to be inverted. is the characteristic scale of the j-th independent parameter to be inverted; The open boundary number that needs to be inverted. For control point numbers, For the first The total number of control points for the open boundary; For the first The open boundary in the first Parametric values ​​at each control point; The weight of the fitting term for the i-th water level observation station is the reciprocal of the variance of the water level observation error at that station. For the first The weight of the fitted term for the flow observation station is the reciprocal of the variance of the flow observation error at that station; Let be the weight of the regularization term for the j-th independent parameter to be inverted, and let its value be such that the ratio of the regularization term to the sum of all observed fitted terms at the initial iteration is . ~ ; The weights for the boundary smoothing term are set such that the ratio of the smoothing term to the sum of all observed fitted terms is equal to the order of magnitude of the sum of the terms in the initial iteration. ~ .

9. The method according to claim 2, characterized in that: In step S2, the water level at the boundary of each beach area is estimated using water level observation or prior value at the start time. It is compared with the lowest topographic elevation of the area. If the difference exceeds the tolerance threshold, it is judged as wet; otherwise, it is judged as dry. When a beach is adjacent to multiple water bodies, the highest water level is taken.

10. The method according to claim 2, characterized in that, In step S5, the global initial state generated by the proxy model is mapped to the computational grid of the basic model according to the following rules: the water level of the subcritical flow section is directly mapped to the computational grid corresponding to the basic model; the lake water level is uniformly assigned to all computational grids within the lake area; the beach water level is mapped to the corresponding computational grid through interpolation; for the portion of the computational grid of the basic model that exceeds the computational domain of the proxy model, the water level is taken as the water level of the nearest grid cell at the boundary of the computational domain of the proxy model; the value of the water level after mapping in each subdomain that is lower than the terrain elevation is truncated to the terrain elevation; the cross-sectional average flow velocity of the subcritical flow section is mapped to the flow velocity component of the corresponding computational grid according to the river channel direction, and the initial flow velocity of the lake and beach is zero. The boundary condition for free evolution is that the value at the end of the inversion window remains unchanged. The termination sign is that the rate of change of water level in the whole area is lower than the preset threshold and the preset number of steps is continuously met after excluding units with water depth lower than the preset minimum water depth threshold, or the preset maximum time limit is reached. After the free evolution ends, the boundary conditions smoothly transition to the real-time value according to the following cosine function within the preset transition period: ， ， Where t is the current time, These are weighting coefficients. For transitional boundary conditions, These are the boundary condition values ​​at the end of the inversion window. These are real-time, time-varying boundary condition values. and These are the start and end times of the preset transition period, respectively.

Images