Flood conduction prediction method and system coupling three defense lines and multi-source information

Abstract

The application discloses a flood conduction prediction method and system coupled with three-line multi-source information. The method comprises the following steps: acquiring meteorological forecast rainfall, rainfall station measured rainfall and hydrological station measured water level and flow data; performing space-time scale normalization processing on the multi-source rainfall, and constructing a corrected forecast rainfall, and splicing and constructing a mixed rainfall driving sequence according to a time partition structure; extracting a trunk line in a river channel to perform feature point identification and constraint point supplementing, generating cross section data to construct a hydrodynamic terrain basic grid; inputting the mixed rainfall driving sequence into a cascade coupling model based on the terrain grid configuration to deduce, and updating the initial state of the hydrodynamic module through a spatialization gain mechanism by using the measured data; advancing to a preset period to perform rolling prediction at a prediction time, and outputting a water level and flow space-time prediction result of a target river section in a prediction period. The application realizes the dual benefits of prolonging the effective prediction period and improving the prediction accuracy, and can provide effective support for fine flood control decision-making.

Description

Technical Field

[0001] This invention belongs to the field of hydrology, water resources and flood forecasting technology, and more specifically, relates to a flood transmission forecasting method and system that couples multi-source information from three lines of defense. Background Technology

[0002] Timely and continuous flood forecasting is of great technical importance in ensuring flood control safety along the coast. Traditional forecasting models mostly rely on lumped hydrological simulation methods, whose calculation accuracy is limited by the density and timeliness of surface rain gauges.

[0003] Currently, the field of hydrological forecasting often employs a single-station monitoring-single-point forecasting model, which uses observational data from discrete ground-based rain gauge stations to drive hydrological models and extrapolate the flow at a single cross-section over time. However, in practical applications, this model faces an inherent contradiction between predictable lead time and forecast accuracy: if relying solely on measured rainfall, the forecast lead time is often relatively short due to the rigid constraints of the basin's runoff generation and confluence time; if macro-meteorological forecasts are introduced to drive the forecast, the forecast results often diverge rapidly due to systematic spatial underestimation or overestimation biases in meteorological forecasts. Furthermore, the traditional single-point forecasting framework cannot solve for the continuous flow movement along the river channel, resulting in severe spatial blind spots in riparian inundation risk areas not covered by measured hydrological stations.

[0004] Therefore, traditional flood forecasting methods struggle to balance long-term forecast accuracy with spatiotemporal dimensions when dealing with river scenarios lacking measured survey data and with sparse hydrological station networks. Furthermore, they remain insufficient in their ability to jointly assimilate multi-source meteorological and hydrological information and perform continuous extrapolation across the entire basin. Thus, it is necessary to research a method that can effectively integrate monitoring and forecasting information across spatiotemporal scales and improve the precision of river flood calculations. Summary of the Invention

[0005] In view of this, this application provides a flood propagation forecasting method and system that couples multi-source information from three lines of defense, in order to solve at least one of the above-mentioned problems existing in the prior art.

[0006] According to one aspect of this application, a flood propagation forecasting method coupling multi-source information from three lines of defense is provided, comprising:

[0007] Acquire meteorological forecast rainfall data, actual rainfall data from rain gauges, and actual water level and flow data from hydrological stations;

[0008] The meteorological forecast rainfall data and the actual rainfall data measured by rain gauges are normalized in terms of time and space. Based on the rainfall deviation characteristics in historical periods, a revised forecast rainfall is constructed. The actual rainfall data measured by rain gauges, the revised forecast rainfall, and the meteorological forecast rainfall data are spliced ​​together according to the time partition structure to construct a hybrid rainfall driving sequence.

[0009] The target river channel midline is extracted using pre-acquired digital elevation data for feature point identification and constraint point supplementation. Based on the supplemented set of control points, river channel cross-section data is generated, and a hydrodynamic topographic basic grid is constructed.

[0010] The mixed rainfall-driven sequence is input into a cascaded coupled model based on hydrodynamic topography grid configuration for propagation and simulation. During the simulation, measured water level and flow data from hydrological stations are used to update the initial state of the hydrodynamic module in the cascaded coupled model through a spatialized gain mechanism.

[0011] The forecast time is advanced by a preset period and the process of constructing and extrapolating the mixed rainfall-driven sequence is executed in a rolling manner, outputting the spatiotemporal forecast results of water level and flow rate of the target river segment within the forecast period.

[0012] According to another aspect of this application, a flood propagation forecasting system coupling multi-source information from three lines of defense is also provided, comprising:

[0013] At least one processor; and a memory communicatively connected to the at least one processor; wherein the memory stores instructions executable by the processor to implement a flood propagation forecasting method that couples three lines of defense and multi-source information, as proposed in this application.

[0014] Based on the technical solution of this application, the present invention achieves the dual benefits of extending the effective forecast period and improving forecast accuracy, which can provide effective support for refined flood control decision-making. Attached Figure Description

[0015] The accompanying drawings, which are included to provide a further understanding of this application and form part of this application, illustrate exemplary embodiments and are used to explain this application, but do not constitute an undue limitation of this application. In the drawings:

[0016] Figure 1 This is a flowchart illustrating a flood propagation forecasting method that couples multi-source information from three lines of defense, as provided in this application.

[0017] Figure 2 This is a flowchart illustrating the spatiotemporal scale normalization process for meteorological forecast rainfall data and actual rainfall data from rain gauges, as provided in this application.

[0018] Figure 3 This is a schematic diagram of the process for constructing revised rainfall forecasts based on rainfall deviation characteristics within historical time periods, as provided in this application.

[0019] Figure 4 This application provides a schematic diagram of the process of extracting the channel line of the target river using pre-acquired digital elevation data, identifying and constraining control points, and generating river cross-section data based on the supplemented set of control points.

[0020] Figure 5 This is a schematic diagram of the process for updating the initial state of the hydrodynamic module in the cascaded coupling model through a spatialized gain mechanism, as provided in this application. Detailed Implementation

[0021] To enable those skilled in the art to better understand the present invention, the technical solutions of the present invention will be clearly and completely described below with reference to the accompanying drawings of the embodiments of the present invention. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort should fall within the scope of protection of the present invention.

[0022] Example 1

[0023] This embodiment provides a flood propagation forecasting method that couples multi-source information from three lines of defense, referring to... Figure 1 This method may include, but is not limited to, the following steps:

[0024] Step 101: Obtain meteorological forecast rainfall data, actual rainfall data from rain gauge stations, and actual water level and flow data from hydrological stations.

[0025] In practical implementation, a three-tiered defense system for rainfall and water level monitoring and forecasting can be established, serving as the data input source to continuously and in real-time access multi-source hydrological and meteorological elements. Meteorological forecast rainfall data can be represented as a three-dimensional grid matrix with spatiotemporal coordinates, provided by the first line of defense—meteorological satellites and rainfall radars. Its effective forecast period can reach 24 to 72 hours, extending the system's information horizon to the long-term forecast axis, providing a broad-area areal rainfall monitoring and advanced forecasting background. Rain gauge measured rainfall data can be represented as a one-dimensional time-series vector of discrete stations, provided by the second line of defense—the surface rain gauge network. It possesses high point rainfall observation accuracy and is used to accurately capture real-time rainfall within the watershed. Hydrological station measured water level and flow data are provided by the third line of defense—the hydrological station network—specifically including one-dimensional time-series vectors of measured water levels and flows at key control sections.

[0026] The three lines of defense complement each other, possessing different accuracy and visibility characteristics across different forecast periods and spatial scales. Together, they form a multi-source information input system that spans forecast periods and accuracy levels, providing complete and closed-loop underlying data support for subsequent propagation forecasts.

[0027] Step 102: Perform spatiotemporal scale normalization on the meteorological forecast rainfall data and the actual rainfall data measured by the rain gauge, and construct a corrected forecast rainfall based on the rainfall deviation characteristics in historical periods. Then, splice the actual rainfall data measured by the rain gauge, the corrected forecast rainfall, and the meteorological forecast rainfall data according to the time partition structure to construct a hybrid rainfall driving sequence.

[0028] To eliminate inconsistencies in spatiotemporal resolution of multi-source rainfall information and extend the forecast period, this step performs data field assimilation and driven reconstruction. After unifying data with different temporal resolutions to the same temporal resolution, inverse distance weighted interpolation is used to transform discrete rain gauge station measured data into a unified spatial grid, thereby achieving deep fusion of multi-source data on the same spatiotemporal reference.

[0029] Based on this, by tracing back to historical periods before the current forecast time, and statistically extracting the deviation trend between measured rainfall and forecast rainfall within the historical window, the systemic underestimation or overestimation bias of rainfall forecasts is captured, and the forecast rainfall in the near future is assimilated and corrected in real time to obtain the revised forecast rainfall.

[0030] Furthermore, the driving field is hard-stitched according to the time partition structure. That is, high-precision rain gauge measured rainfall data is used in the historical period, assimilated and debiased predicted rainfall is used in the near future period, and the original meteorological forecast rainfall data is retained in the far future period.

[0031] The aforementioned three-stage driving structure ensures that the data source with the highest information validity is selected within any forecast period, which helps to alleviate the problem that the effective forecast period is constrained by the time of production and convergence in traditional forecast models, and suppresses the phenomenon that forecast errors diverge as the forecast period lengthens.

[0032] Step 103: Use the pre-acquired digital elevation data to extract the channel line of the target river for feature point identification and constraint point supplementation. Generate river cross-section data based on the supplemented control point set and construct the hydrodynamic topography basic grid.

[0033] To establish a topographical basis for hydrodynamic simulation in river basins lacking measured cross-sectional survey data, this step adaptively reconstructs the river channel geometry and topology based on publicly available satellite digital elevation data. After extracting the midcourse line of the target river channel, feature point identification and thinning based on morphological adaptive changes are performed to preserve the control points at bends that reflect sharp bends and abrupt changes in flow direction. However, directly using the thinned nodes to build the hydrodynamic mesh would result in sparse mesh nodes in the straight river sections between two control points. Consequently, when solving the Saint-Venant equations, the excessively large mesh scale would violate the Courant stability condition for numerical computation, leading to non-convergence or divergence in the calculation.

[0034] To address this, this step introduces a mandatory spacing threshold constraint to supplement nodes. Equidistant control points are automatically inserted between straight river sections where the actual distance exceeds the preset threshold, and these are merged with the original bend control points to form a supplemented control point set. Centered on this supplemented control point set, cross-sectional sampling lines are generated along the normal direction, and digital elevation data are extracted. This constructs a hydrodynamic topographical grid that accurately captures the flow characteristics of bends while ensuring stable numerical calculations for straight sections, providing a topographic data foundation for hydrodynamic modeling in areas lacking data.

[0035] Step 104: Input the mixed rainfall-driven sequence into the cascaded coupling model based on the hydrodynamic topography grid configuration for propagation and simulation. During the simulation, use the measured water level and flow data from the hydrological station to update the initial state of the hydrodynamic module in the cascaded coupling model through a spatialized gain mechanism.

[0036] After reconstructing the terrain mesh and driving sequence, a cascaded coupling model is used to achieve a complete physical mapping from watershed surface rainfall to water level and flow rate processes at any cross-section throughout the watershed. The cascaded coupling model seamlessly cascades runoff generation and concentration calculations, reservoir and dam scheduling rules, and the solution of the one-dimensional Saint-Venant equations. The mixed rainfall driving sequence, used as boundary condition input, is passed through the cascaded models at each level to drive the flood forward propagation within the river network.

[0037] To mitigate the systematic and time-amplified propagation errors accumulated in historical runoff generation and hydrodynamic calculations due to uncertainties in various physical parameters, a state constraint mechanism is triggered upon receiving measured water level and flow data from hydrological stations. This mechanism calculates the state deviation between the measured and model-calculated values ​​and, utilizing a spatialized gain mechanism that considers flood wave velocity, distributes the error correction at the measuring point upstream and downstream along the river channel. This ensures that the overall flow field state variables of the hydrodynamic module are reasonably corrected at the current moment, improving the physical accuracy of the initial state in subsequent forward forecasts and indirectly enhancing the accuracy of flood projections over long-term prediction periods.

[0038] Step 105: Advance the forecast time with a preset period (i.e., preset rolling period) and execute the mixed rainfall-driven sequence construction and propagation process in a rolling manner, and output the spatiotemporal forecast results of water level and flow rate of the target river section within the forecast period.

[0039] As time progresses and new observational data is continuously incorporated, the current forecast time is automatically advanced at a preset fixed time cycle, entering the next forecast cycle. In each new rolling cycle, the observed area extends into the future, the assimilation window shifts to the right accordingly, the bias correction parameters are refreshed using newly acquired measured rainfall data, and the initial flow field state of the cascaded coupled model is also updated in real time using the latest acquired measured water level and flow rate data from hydrological stations. Based on the refreshed driving field and initial state, the model again extrapolates the complete forecast period.

[0040] With this cyclical rolling mechanism, under the condition that effective measured data can be obtained in each round and there are no abnormal jumps in the forecast input, due to the continuous injection of high-precision measured information, the forecast uncertainty of the same target flood peak time shows an overall decreasing trend throughout the entire life cycle. This achieves the effect of continuous convergence of forecast accuracy with rolling progress, thereby continuously outputting spatiotemporal forecast results of water level and flow with a high Nash efficiency coefficient to the flood control decision terminal, covering any section of the entire river.

[0041] According to another aspect of this application, a flood propagation forecasting system coupling multi-source information from three lines of defense is also provided, comprising:

[0042] At least one processor; and a memory communicatively connected to the at least one processor; wherein the memory stores instructions executable by the processor to implement a flood propagation forecasting method that couples three lines of defense and multi-source information, as proposed in this application.

[0043] Example 2

[0044] Based on the above embodiment 1, this embodiment further describes the spatiotemporal scale normalization process in detail.

[0045] In one possible implementation, refer to Figure 2 Spatiotemporal scale normalization of meteorological forecast rainfall data and actual rainfall data from rain gauges can include, but is not limited to, the following steps:

[0046] Step 201: Unify the meteorological forecast rainfall data and the actual rainfall data measured by the rain gauge to the same time resolution; obtain the time-normalized meteorological forecast rainfall data and the time-normalized actual rainfall data measured by the rain gauge.

[0047] Due to inherent differences in the acquisition frequency of multi-source hydrological observation data—for example, meteorological satellite data may have a long update cycle, while ground rain gauge data is collected at a relatively high frequency—time base calibration is required before cross-source fusion. The precipitation products from each source are decomposed or aggregated using a pre-defined hydrodynamic model calculation step size. When the resolution of the original input data is lower than the target calculation step size, a uniform distribution method is used to smoothly decompose the precipitation into various sub-time steps. When the resolution of the original input data is higher than the target calculation step size, cumulative summation is performed according to time windows to ensure that the time base of all data streams remains consistent, outputting a normalized precipitation sequence with equal time steps.

[0048] Step 202: Use spatial interpolation to interpolate the measured rainfall data from the rain gauge station to the spatial grid corresponding to the meteorological forecast rainfall data to generate a gridded measured rainfall reference field.

[0049] Specifically, spatial interpolation is used to interpolate the time-normalized rain gauge measured rainfall data to the spatial grid corresponding to the time-normalized meteorological forecast rainfall data, thereby generating a gridded measured rainfall reference field.

[0050] After aligning the time reference, it is necessary to address the issue that ground rainfall stations are spatially discrete and lack grid topological characteristics. The two-dimensional spatial grid resolution parameters contained in the meteorological forecast rainfall data are extracted as a global projection reference. For each grid node to be calculated in this reference coordinate system, in this embodiment, the inverse distance weighting method is used as a specific implementation of the spatial interpolation method to perform spatial interpolation calculations on the surrounding discrete measured station data, outputting an intermediate layer with continuous surface rainfall distribution characteristics. Kriging interpolation or other spatial interpolation methods can also be used.

[0051] As an optional specific implementation method, the spatial grid resolution corresponding to the meteorological forecast rainfall data can be set to 1km*1km. The station's measured values ​​are mapped to the grid system by the inverse distance weighting method to obtain a gridded measured rainfall reference field with the same resolution as the forecast rainfall, providing a structural alignment basis for subsequent isotopic grid algebra operations.

[0052] Step 203: In the area covered by the measured rainfall data from the rain gauge station, the (time-normalized) meteorological forecast rainfall data is locally corrected using a gridded measured rainfall reference field. In areas without measured coverage, the rain-measuring radar monitoring portion of the (time-normalized) meteorological forecast rainfall data is retained. In areas without rain-measuring radar monitoring, the satellite forecast portion of the (time-normalized) meteorological forecast rainfall data is retained, forming a multi-level nested spatial fusion rainfall field.

[0053] To provide a continuous rainfall-driven background across the entire spatial domain while ensuring data reliability, a conditional fusion operation with decreasing priority is performed. Spatially, effective coverage areas are divided with each rain gauge station as the center and a preset distance as the radius. Within the coverage area, the confidence level of the measured data is at the highest level, and the gridded measured rainfall reference field is directly used to replace or numerically smooth the corresponding meteorological forecast grid values ​​within that area. For areas without measured coverage beyond this radius limit, rainfall radar monitoring grid data is used for interpolation; if the geographic grid has neither rain gauge stations nor is it in a radar beam detection blind zone, the meteorological satellite forecast data with the widest coverage is used as a fallback rainfall background field. Relying on a hierarchical spatial nesting mechanism, a multi-level nested spatially fused rainfall field without data gaps is output.

[0054] In some alternative implementations, the preset distance for delineating the coverage area can be set to 30km to establish a trust boundary for the measured data and prevent false rainfall information from being introduced by spatial interpolation calculations over excessively long distances.

[0055] Example 3

[0056] Building upon Example 2 above, this example further details the rainfall bias assimilation correction process based on an adaptive window. This example focuses on addressing the feature lag caused by a fixed historical time window during the rainfall forecast assimilation process, as well as the numerical divergence problem resulting from a single proportional correction.

[0057] In one possible implementation, a revised rainfall forecast is constructed based on rainfall deviation characteristics within historical time periods, referencing... Figure 3 This may include, but is not limited to, the following steps:

[0058] Step 301: At the current forecast time, based on the trend of the deviation between the measured rainfall data from the rain gauge and the forecast rainfall data, determine the adaptive time window for the historical period. This specifically includes the following process:

[0059] Calculate the spatial average of the absolute deviation of all spatial grid points at each time point, and construct the surface average absolute deviation sequence.

[0060] Extract the one-lag autocorrelation coefficient of the surface mean absolute deviation sequence, and use it (the one-lag autocorrelation coefficient) to perform dynamic interpolation between the preset lower limit and upper limit of the window to determine the length of the adaptive time window.

[0061] As each rolling forecast cycle progresses to the current forecast time, the statistical deviation characteristics between the forecasted precipitation field and surface observation data over a past period are first analyzed. To enable the system to dynamically adjust the information extraction range based on the suddenness and persistence of precipitation events, this method does not use a fixed historical backtracking time, but instead utilizes the error sequence characteristics within historical periods to dynamically divide the time boundaries.

[0062] Specifically, a preset maximum backtracking period is loaded, and the spatial average of the absolute deviation of all spatial grid points at each time point is calculated to construct a surface average absolute deviation sequence.

[0063] Furthermore, to quantitatively assess the persistence or variability of the bias over time, the lagged one-step autocorrelation coefficient of the surface-mean absolute bias sequence is extracted. The formula for calculating the lagged one-step autocorrelation coefficient is as follows:

[0064] ;

[0065] in, The autocorrelation coefficient is the one-step lag, where t0 is the current forecast time, and T is the time lag.max Let D(t) be the preset maximum backtracking time, D(t) be the spatial average of the absolute deviation of all spatial grid points at time t, and D(t-1) be the spatial average of the absolute deviation of all spatial grid points at time t-1. This represents the mean of the surface-averaged absolute deviation sequence over the period of maximum backtracking. From time t0-T max The summation sign for time traversal from +1 to time t0. From time t0-T max The summation symbol for time traversal at time t0.

[0066] As an example, the maximum backtracking time T max The time window can be set to 12 hours. After extracting the one-step autocorrelation coefficient, it is dynamically interpolated between the preset lower and upper window limits to determine the length of the adaptive time window. The specific calculation is as follows:

[0067] T w (t0)=T w,min +(T w,max -T w,min )*max(0,ρ1);

[0068] Among them, T w (t0) represents the length of the adaptive time window determined at the current forecast time t0, where T w,min T is the preset lower limit of the window. w,max ρ1 is the preset upper limit of the window, max is the maximum value operator, and ρ1 is the autocorrelation coefficient with one step lag.

[0069] In this embodiment, the lower limit of the window T w,min Set to 2h, window limit T w,max The time window is set to 12 hours. When the deviation exhibits high persistence (i.e., the correlation coefficient approaches 1), the calculated window length shrinks towards the upper limit to fully utilize the stable deviation trend information over a long period. When the deviation exhibits rapid alternation or abrupt changes (i.e., the correlation coefficient approaches 0 or is negative), the calculated window length is shortened to the lower limit to avoid introducing outdated historical rainfall information. This outputs the adaptive time window for the current moment.

[0070] Step 302: Within the adaptive time window, extract the relative deviation factor and absolute deviation of each spatial grid point at each time step.

[0071] After determining the length of the adaptive time window, a spatially consistent historical time period boundary is defined up to the current time. For each discrete time step and each corresponding grid coordinate point within this time period, statistical errors at two scales are extracted simultaneously. The absolute deviation is calculated using the difference between the measured and predicted values. To eliminate the masking effect of varying rainfall baselines on the magnitude of the deviation, a relative deviation factor for each grid point is calculated concurrently.

[0072] α x,y,t =(P obs (x,y,t)-P fcst (x,y,t)) / (P fcst (x,y,t)+ε);

[0073] β x,y,t =P obs (x,y,t)-P fcst (x,y,t);

[0074] Where, α x,y,t β is the relative deviation factor of the grid point at coordinate position (x,y) at time t. x,y,t P is the absolute deviation of a grid point at coordinate position (x, y) at time t. obs (x,y,t) represents the measured rainfall value at time t for a grid point with coordinates (x,y). fcst (x,y,t) represents the forecast rainfall value at time t for the grid point at coordinate position (x,y), and ε is a preset small positive number used to prevent division by zero errors.

[0075] In this embodiment, the tiny positive number ε is set to 10. -6 This is used to ensure that the denominator of the relative deviation factor is mathematically valid at grid points where the predicted rainfall is 0. By traversing all grids and time steps within the adaptive time window, the relative and absolute deviation field space matrices are obtained.

[0076] Step 303 involves using a time decay mechanism to perform a weighted average of the extracted relative deviation factor and absolute deviation, resulting in a weighted relative deviation and a weighted absolute deviation. This specifically includes the following process:

[0077] The attenuation coefficient is determined based on the length of the adaptive time window, and a time weighting function that decays exponentially with increasing time distance is constructed.

[0078] By using a time weighting function to assign higher weights to deviation samples that are closer to the current forecast time, the relative deviation factor and the absolute deviation are weighted and summed to obtain the weighted relative deviation and the weighted absolute deviation.

[0079] To enhance the sensitivity of the rainfall deviation assimilation process to recent error evolution, this step introduces time weights that decay exponentially with increasing time span to the historical samples within the window. A time decay weight function is constructed based on the length of the adaptive time window. The time weight function is calculated as follows:

[0080] ω t =exp(-φ(t0)*(t0-t)) / Σ k exp(-φ(t0)*(t0-t k ));

[0081] Where, ω t Let be the decay weight value corresponding to the historical time t, exp be the exponential function with the natural constant as the base, t0 be the current forecast time, and t be the current historical time being calculated. k For any historical moment traversed within the adaptive time window, Σ k The sign is the summation of all historical moments within the adaptive time window, and φ(t0) is the attenuation coefficient corresponding to the current forecast time t0.

[0082] It should be noted that, in order to establish a rigid physical constraint between the decay rate and the window range, the decay coefficient φ(t0) is determined by the window half-life, and is calculated as follows:

[0083] φ(t0)=ln2 / T half ;

[0084] Where ln2 is the value of the natural logarithm, T half The half-life length is a preset value. As a preferred option, the half-life length is constrained to be half the length of the adaptive time window.

[0085] Using the constructed time weighting function, a weighted summation is performed on the relative deviation factor sequence and the absolute deviation sequence within the adaptive time window. By multiplying the deviation at each historical moment by its corresponding weight value and summing the results, the weighted relative deviation and weighted absolute deviation, representing the recent average error trend, are output respectively.

[0086] Step 304: Based on the weighted relative bias and weighted absolute bias, perform a hybrid correction operation on the forecast rainfall within the future time window, and output the corrected forecast rainfall. This specifically includes the following process:

[0087] Step 3041: For the forecast rainfall amount of the forecast rainfall field within the future time window, multiply the weighted relative deviation by the forecast rainfall amount to obtain the relative correction amount;

[0088] Step 3042: Compare the absolute value of the relative correction with the absolute value of the weighted absolute deviation; if the absolute value of the relative correction is less than the absolute value of the weighted absolute deviation, use the relative correction as the target correction; otherwise, use the absolute value of the weighted absolute deviation after taking the sign of the relative correction as the target correction.

[0089] Step 3043: The target correction amount is superimposed on the forecast rainfall amount and non-negative truncation is performed to obtain the corrected forecast rainfall.

[0090] After obtaining the bias characteristics of the weighted average, it is used as a correction operator for the short-term future period after the current time. In order to prevent numerical divergence or abnormal amplification of the calculation results due to a single proportional correction when encountering sudden heavy rain or sporadic light rainfall in conventional rainfall assimilation methods, this step implements a nonlinear hybrid truncation correction mechanism.

[0091] For each forecast grid point in the forecast rainfall field within the future time window, the weighted relative deviation of that point is first multiplied by the original forecast rainfall to calculate the relative correction. The relative correction is then compared with the corresponding weighted absolute deviation of that point.

[0092] In the first branch condition, if the calculated relative correction is less than the weighted absolute deviation, it indicates that the relative error term has not experienced scale anomaly, and the relative correction is determined as the target correction.

[0093] In the second branch case, if the calculated relative correction is greater than or equal to the weighted absolute deviation, it indicates that the original forecast base may be severely distorted on the predetermined grid. Continuing to use relative proportional multiplication will cause the correction value to be amplified exponentially. Therefore, the absolute deviation-dominated logic is forcibly switched to determine the weighted absolute deviation as the target correction.

[0094] As a supplementary explanation for the second branch scenario, suppose at another spatial grid point, the normalized value of the original forecast rainfall is 10.0, the weighted relative bias is -0.8, and the weighted absolute bias is -3.0. The relative correction is calculated as (-0.8) × 10.0 = -8.0, and the absolute value is taken as 8.0. Comparing 8.0 with the absolute value of the weighted absolute bias, 3.0, since 8.0 > 3.0, the absolute bias-dominated logic is triggered, and the target correction is -3.0. The corrected rainfall is 10.0 + (-3.0) = 7.0, and after non-negative truncation, 7.0 is output. This mechanism also prevents excessive reduction of the product term in the overestimation scenario.

[0095] After determining the target correction amount, it is superimposed onto the original forecast rainfall at that point. To ensure that the corrected physical quantities conform to the fundamental law of nonnegativity in hydrological evolution, the superposition result is truncated using a maximum value operator. The formula for the above hybrid correction operation is expressed as follows:

[0096] P corr (x,y,t)=max(P fcst (x,y,t)+ΔP(x,y,t),0);

[0097] Among them, P corr (x,y,t) represents the revised forecast rainfall value corresponding to the grid point at location (x,y) at future time t, P fcst (x,y,t) represents the original forecast rainfall value for this point, ΔP(x,y,t) represents the target correction value for this point, and max is the maximum value operator.

[0098] To illustrate the logic of the hybrid correction operation, a simplified example with normalization is provided. Assume that at a spatial grid point in a future time period, the normalized value of the original forecast rainfall is 1.0. If the weighted relative deviation for this point, calculated through previous steps, is 9.0, and the weighted absolute deviation is 5.0, the relative correction is calculated as 9.0 * 1.0 = 9.0. This relative correction of 9.0 is compared with the weighted absolute deviation of 5.0. Since 9.0 > 5.0, the absolute deviation-dominated logic is triggered, the system cuts off the relative correction path, and directly limits the target correction to the weighted absolute deviation of 5.0. This target correction is then added to the original forecast rainfall, resulting in 1.0 + 5.0 = 6.0. After non-negative truncation, the corrected forecast rainfall is output as 6.0. This mechanism maintains dynamic error correction while preventing the divergence of the product term.

[0099] Example 4

[0100] Based on the above embodiment 1, the process of extracting the midline of the target river channel using pre-acquired digital elevation data, identifying feature points and constraining point supplementation, and generating river channel cross-section data based on the supplemented set of control points is further described in detail.

[0101] In one possible implementation, refer to Figure 4 This process may include, but is not limited to, the following steps:

[0102] Step 401: Extract the midline of the target river channel based on digital elevation data.

[0103] When constructing the geometric boundaries of the hydrodynamic model, the acquired digital elevation data (DEM) provides a topographic reference surface for the entire watershed. Flow direction tracing and runoff accumulation calculations are performed on this DEM data, extracting grid paths where the runoff accumulation exceeds a preset area threshold. This generates a two-dimensional vector line representing the main flow channel, i.e., the midchannel line of the target river. This line segment serves as the longitudinal distance reference axis for one-dimensional hydrodynamic calculations in subsequent steps.

[0104] Step 402: Perform recursive thinning processing based on thinning threshold on the middle channel of the target river to identify the set of bend control points with obvious morphological changes.

[0105] To eliminate redundant spatial nodes while preserving the river channel's geometric features, a recursive vector simplification logic, such as the classic Douglas-Peucker algorithm, is used to thin the nodes along the channel's midline. Specifically, a baseline is formed by connecting the first and last nodes, and the vertical distance from each intermediate node to this baseline is calculated. When the maximum vertical distance exceeds a preset thinning threshold, the node corresponding to this maximum distance is retained, and the curve is divided into two segments using this point as the boundary. This recursive judgment continues until the distance from all unretained nodes to their corresponding baselines is less than or equal to the thinning threshold. The retained nodes constitute the set of bend control points. These nodes spatially correspond to sharp bends or abrupt changes in the river's direction, reflecting the local hydraulic resistance caused by changes in flow direction.

[0106] In some alternative implementations, the thinning threshold is determined by a preset ratio based on the minimum curvature of the channel in the target river channel to preserve key control points in the bend. As an example, the thinning threshold is set to be less than 10% of the minimum curvature of the channel, thereby maintaining a high density of geometric control points in bend areas with significant topographic relief.

[0107] Step 403: Calculate the actual distance between adjacent control points in the bend control point set. For straight river sections where the actual distance is greater than the preset spacing threshold, add equidistant intermediate control points between the corresponding adjacent control points and merge them to form a supplemented control point set.

[0108] After the preceding recursive thinning process, the removal of a large number of nodes in the gently sloping straight river section leads to excessively large physical spacing between adjacent control points. If such sparse nodes are directly used to generate the computational mesh, the model is prone to divergence during hydrodynamic differential solutions due to the spatial step size exceeding the stability boundary. To mitigate this computational risk, the computational unit performs a spacing scan on the set of bend control points.

[0109] In one specific implementation, equidistant intermediate control points are added between adjacent control points. This includes: dividing the actual distance by a preset spacing threshold and rounding up to determine the number of equal divisions; and determining the number of intermediate control points to be added and their spatial coordinates based on the number of equal divisions. The calculation process is as follows:

[0110] N d =ceil(d / d thresh );

[0111] Where, N d The number of equal parts is denoted by 'ceil', which is the floor function, and 'd' is the actual distance between adjacent control points. thresh This is the preset spacing threshold.

[0112] After obtaining the number of equal parts, calculate the number of intermediate control points required, which is N. d -1. Using linear interpolation, generate a corresponding number of intermediate control point coordinates between adjacent control points, ensuring that the distance between any two adjacent nodes after insertion is less than or equal to a preset distance threshold. Merge and sort the original turning control points with the newly generated intermediate control points, and output the supplemented control point set. As a specific parameter setting example, the preset distance threshold can be set to 500m.

[0113] Step 404: Using each control point in the supplemented control point set as the center, generate cross-sectional sampling lines along a direction perpendicular to the channel line of the target river.

[0114] After obtaining a spatially distributed control point topology, the supplemented control point set is traversed. For each control point, its tangent vector on the middle channel of the target river is calculated, and its orthogonal normal vector is obtained. Using this control point as the reference coordinate, bidirectional extensions are made along the orthogonal normal vector to both banks of the river, generating a cross-sectional sampling line with a line segment shape. When performing this extension operation, the total width parameter of the sampling line is set to be greater than or equal to the maximum flood discharge width corresponding to the historical highest water level of the target river section, thereby ensuring that the generated cross-section can cover the inundation range under extreme flood conditions.

[0115] Step 405: Extract elevation data from digital elevation data along the cross-section sampling line according to the preset sampling interval to generate river cross-section data that meets the topology requirements of hydrodynamic simulation.

[0116] After establishing the lateral extraction path, discrete sampling point arrays are generated on each cross-sectional sampling line according to the preset sampling interval. For the two-dimensional plane coordinates of each sampling point, a spatial query is initiated to the digital elevation data layer to extract the corresponding surface elevation value.

[0117] In some optional implementations, the preset sampling interval is set to half the spatial resolution of the digital elevation data to meet the Nyquist sampling theorem's requirement for capturing the frequency of terrain undulations. After extracting the elevation data, a terrain feature point recognition algorithm is executed on each cross-section to search for the lowest elevation point as the riverbed coordinates, and to search for the highest local elevation points on both the left and right banks as the riverbank top boundaries. After filtering out invalid elevation points outside the boundaries, the hydraulic connection topology between each cross-section is established, structured, and exported, thereby generating river cross-section data that supports the solution of the one-dimensional Saint-Venant equations.

[0118] Example 5

[0119] Based on the above embodiment 4, the process of inputting the mixed rainfall-driven sequence into the cascaded coupling model based on the hydrodynamic topographic grid configuration for propagation and deduction is further described in detail.

[0120] In one possible implementation, the process may include, but is not limited to, the following steps:

[0121] Step 501: Input the mixed rainfall-driven sequence into the runoff generation and confluence model for simulation, calculate and output the boundary conditions of each tributary in the target watershed, as well as the reservoir inflow data.

[0122] The cascaded coupling model, as the core of the entire forecasting system, is responsible for mapping rainfall physical quantities to river hydrodynamic boundary parameters. The reconstructed mixed rainfall driving sequence is input into the front-end runoff generation and confluence module to drive the rainfall-runoff conversion process. This embodiment provides several parallel alternatives for the specific structure of the runoff generation and confluence model.

[0123] Optionally, the runoff generation and confluence model includes any one of the following: a watershed hydrological model based on physical mechanisms, a lumped runoff generation and confluence model, or a runoff generation and confluence inference model based on machine learning.

[0124] In the first alternative, the runoff generation and confluence model employs a watershed hydrological model based on physical mechanisms, specifically a distributed or semi-distributed runoff generation and confluence model. This approach discretizes the target watershed into multiple hydrological computational units with independent physical properties. For each computational unit, based on the corresponding land use type, soil permeability, and vegetation interception parameters, the net rainfall of that unit is calculated. The surface and subsurface runoffs are then collected into the tributary channels using kinematic wave or unit hydrograph methods, outputting the lateral boundary conditions for each tributary.

[0125] In the second alternative, a lumped runoff model is adopted, specifically the Xin'anjiang model. This scheme treats the predetermined tributary catchment area as a whole. The total water yield of the basin is calculated based on the principle of full-storage runoff, and the total water yield is separated into surface runoff, interflow, and groundwater runoff using a three-source division mechanism. Finally, the lateral flow process into the main channel is obtained through linear reservoir thrust flow calculation.

[0126] In the third alternative, the runoff generation and confluence model employs a machine learning-based runoff generation and confluence inference model. It utilizes a pre-trained deep neural network to obtain the temporal characteristics of the mixed rainfall-driven sequence and, based on nonlinear mapping weights learned from historical data, directly outputs the flow process sequence of each tributary.

[0127] In one specific implementation, the deep neural network can employ architectures capable of extracting temporal features, such as long short-term memory networks, gated recurrent units, or temporal convolutional networks. It uses a gridded rainfall sequence within a time window as the input feature and the flow process at the corresponding tributary control section as the output label. Offline training is performed by minimizing the mean squared error loss function. Those skilled in the art can use conventional methods such as backpropagation and stochastic gradient descent to complete the model training, and the hyperparameters during training can be determined using conventional hyperparameter tuning methods.

[0128] Furthermore, to ensure the operational capability of the aforementioned runoff generation and confluence model with physical mechanisms in river regions lacking long-term hydrological observation records, an offline parameter calibration procedure was implemented before real-time simulations. Based on the area and bottom slope topography of the target watershed, initial parameter values ​​were extracted from similar watersheds with abundant hydrological data to complete regional transplantation. Hydrological parameters of similar watersheds can be obtained from publicly available hydrological yearbooks, water resources survey data, or academic literature. Those skilled in the art can select the reference watershed with the highest similarity based on physical characteristics such as watershed area, slope, and soil type.

[0129] When a limited number of measured flood data sessions are received, the Hybrid Complex Evolutionary Algorithm (SCE-UA) global optimization algorithm is used to iteratively fine-tune sensitive parameters such as free water storage capacity and groundwater retreat coefficient, outputting a converged parameter system. This pre-set process ensures the numerical reliability of the runoff generation and confluence model in data-scarce scenarios.

[0130] Step 502: Input the reservoir inflow data output from the runoff generation and confluence model into the reservoir scheduling model for calculation, and calculate and output the upper boundary conditions.

[0131] For main stream sections within a basin containing large-scale water conservancy projects, upstream flow is heavily influenced by artificial scheduling rules. Inflow data from the runoff generation and collection model is extracted and imported into the reservoir scheduling model. Pre-configured reservoir level-capacity relationship curves, floodgate operation procedures, and the current reservoir starting water level are loaded. The outflow at each time step is calculated by solving the water balance equation. This outflow sequence serves as the starting input for hydrodynamic extrapolation, outputting the upper boundary conditions. If the target basin has no reservoir nodes, the runoff generation and collection output from the upstream catchment area of ​​the main stream is directly used as the upper boundary conditions.

[0132] Step 503: Input the upper boundary conditions and the lateral boundary conditions of each tributary into the hydrodynamic simulation model, and use the hydrodynamic topographic base grid as the topographic boundary to continuously simulate and calculate the propagation process of flood in the river channel in order to solve the flow state along the river.

[0133] After obtaining the boundary driving forces and topographic basis, the numerical solution for river flood evolution is initiated. The upper boundary conditions are configured at the top node of the main channel, and the boundary conditions of each tributary are configured at the corresponding confluence nodes according to spatial topological coordinates. Downstream boundary conditions are configured at the end node of the main channel. The downstream boundary conditions can be determined using a water level-discharge relationship curve or a normal water depth assumption. Those skilled in the art can select an appropriate downstream boundary type based on the actual hydraulic conditions at the outlet of the target river segment. A hydrodynamic topographic base grid is loaded as the spatial computational domain. Regarding the selection of the hydrodynamic simulation model, this embodiment also provides several parallel alternatives to adapt to different computing power conditions and flood control accuracy requirements.

[0134] Optionally, the hydrodynamic simulation model includes any one of the following: a one-dimensional hydrodynamic simulation model, a one- or two-dimensional coupled hydrodynamic simulation model, or the Muskingen hydrological evolution model.

[0135] In the first application scenario, a one-dimensional hydrodynamic simulation model is adopted. A numerical model is established based on the Saint-Venant equations governing fluid motion, and the equations are discretized and solved on a hydrodynamic topographic grid using an implicit finite difference scheme. This process iteratively calculates the mass and momentum conservation of the water flow at each time step, outputting the water level and flow rate distributed longitudinally along the river channel.

[0136] In the second application scenario, a coupled one-dimensional and two-dimensional hydrodynamic simulation model is adopted. In this scheme, the main channel area is solved using one-dimensional equations, while the floodplains on both sides of the channel are solved using two-dimensional shallow water equations. When the one-dimensional channel water level exceeds the levee crest elevation, a two-way exchange of water mass and momentum occurs between the one-dimensional and two-dimensional grids through the side weir flow formula, thereby refining the spatial inundation process after the flood overflows the levee.

[0137] In the third application scenario, the hydrodynamic simulation model adopts the Muskingen hydrological evolution model. When computational resources are limited or high-precision topographic data is lacking, the flow rate for each river segment is calculated directly by time-series recursion based on the continuity equation and the channel storage equation, using preset channel storage coefficients and specific gravity coefficients. This approach sacrifices rigorous solution to the hydrodynamic equations in exchange for increased computational speed.

[0138] Through the aforementioned parallel optional mechanism, the hydrodynamic simulation model continuously outputs the flow state along the entire river section under different boundary constraints.

[0139] Example 6

[0140] Based on Embodiment 5 above, this embodiment further details the process of updating the initial state of the hydrodynamic module in the cascaded coupled model using measured water level and flow data from hydrological stations through a spatialized gain mechanism. This embodiment focuses on solving the numerical oscillation and overcorrection problems caused by overlapping spatial decay of state variables during weighted correction in densely populated river sections with multiple monitoring stations.

[0141] In one possible implementation, refer to Figure 5 This process may include, but is not limited to, the following steps:

[0142] Step 601: When the measured water level and flow rate data of the hydrological station are received, compare the measured water level and flow rate data of the hydrological station with the free calculation state of the corresponding section of the cascaded coupling model, and extract the state deviation.

[0143] When the cascaded coupled model reaches the predetermined update time, the assimilation update module detects the latest transmission status of the third-line hydrological station network. If newly input measured water level and flow rate data from a hydrological station are obtained, the cross-sectional node in the cascaded coupled model with the same geographical location as the hydrological station is located, and the water level value generated by the model's free calculation without external constraints is loaded. The measured water level value in the hydrological station's measured water level and flow rate data is compared with the water level calculated by the model in its free calculation state, and the two are subtracted to extract the water level state deviation, which serves as the basic error driving source for subsequent spatialized gain adjustment.

[0144] Step 602: Based on a preset spatial correlation scale, construct a spatialized gain function that decreases exponentially with increasing spatial distance.

[0145] To reasonably interpret and propagate the measured correction amount at a single hydrological station to surrounding blind sections lacking monitoring capabilities, the assimilation and update module establishes spatial correlations. Based on the boundary influence law, the diffusion influence weight of this correction amount exhibits an exponential decay as the physical distance between the target section and the measured hydrological station increases.

[0146] The attenuation rate of the spatialized gain function is determined by the spatial correlation scale, which is defined as the product of the flood physical wave velocity and the forecast rolling period (i.e., the preset rolling period). The flood physical wave velocity is estimated using empirical formulas based on the hydraulic radius of the representative cross-section in the hydrodynamic topography grid and the riverbed slope corresponding to the grid. Within a rolling period, the physical distance that the flood wave can propagate determines the spatial validity and timeliness of the current measured observation information. The spatial correlation scale is calculated as follows:

[0147] L s =c w *T roll ;

[0148] Among them, L s For spatially relevant scales, c w To estimate the obtained physical wave velocity of the flood, T roll This is a preset rolling period. Wherein, the physical wave velocity of the flood, c... w Further estimation was performed using Manning's empirical formula for flow velocity, as follows:

[0149] c w =(5 / 3)*(1 / n)*(R h 2 / 3 )*(S0 1 / 2 );

[0150] Among them, c w R represents the physical wave velocity of the flood, n is the Manning roughness coefficient preset for the target river section, and R is the water velocity. h S0 represents the hydraulic radius of the cross-section, and S0 is the average bottom slope value of the target river segment. This spatial correlation scale, L, with its clear hydraulic background, is utilized. s A basic spatial distribution gain model was constructed.

[0151] Step 603: Use the spatialized gain function to propagate the state deviation upstream and downstream along the river channel to generate a spatially distributed constraint correction.

[0152] For river locations where the influence areas of multiple hydrological stations overlap, a normalized compression operation of summing the distance exponents of multiple stations is introduced into the denominator of the spatialized gain function to suppress excessive state correction in the overlapping areas.

[0153] The spatial gain function constructed in this step, which includes multi-station overlap normalization and deduplication logic, is shown in the following formula:

[0154] ;

[0155] Among them, K j(s) represents the spatialized gain function value of the j-th hydrological station at any cross-sectional position s on the river channel, K0 is the preset local gain coefficient, exp is an exponential function with the natural constant as the base, and |ss| obs,j | represents the distance s along the path between any cross-sectional location s and the j-th hydrological station. obs,j The absolute distance between them, L s The calculated spatial correlation scale is given by max, which is the maximum value operator. The summation symbol is |ss| for iterating over all M hydrological stations within the basin. obs,l | represents the distance s between any cross-sectional location s and the l-th hydrological station. obs,l The absolute distance between them.

[0156] It should be noted that the local gain coefficient K0 reflects the ratio of the model calculation variance to the observation variance. Since the observation accuracy of water level stations in small and medium-sized rivers is typically on the order of centimeters, while the water level simulation error of a one-dimensional hydrodynamic model under data-free conditions is typically on the order of decimeters to meters, the observation error variance is much smaller than the model error variance. Therefore, in some possible implementations, the value of the local gain coefficient K0 is directly simplified to 1. In this case, the gain mapping under single-point water level observation degenerates into a unit mapping. This simplification reduces the computational cost of online covariance calculation.

[0157] By comparing the state deviation of each station with the corresponding spatialized gain function value K j (s) performs multiplication calculations to generate constraint correction quantities in the two-dimensional or one-dimensional longitudinal space of the entire basin, which are the result of the combined effects of various hydrological stations and are normalized and deduplicated by the denominator. The maximum value operator in the denominator can ensure that the complete spatial natural exponential decay pattern is maintained in the single-station control area, while in the multi-station overlapping and dense control area, the overall value is adaptively normalized and compressed by summing multiple denominators, which suppresses numerical oscillations caused by excessive correction of state variables.

[0158] To facilitate a more intuitive understanding of the working mechanism of the spatialized gain function, a typical scenario of a dual-station controlled river section is presented below. Assume the total length of the target river section is 20km, with hydrological station A (located 5km upstream) and hydrological station B (located 15km upstream) deployed along the river. The flood physical wave velocity is estimated to be approximately 2m / s using the Manning formula, and the preset rolling period is 6h. Therefore, the spatial correlation scale Ls = 2 × 6 × 3600 / 1000 = 43.2km. At a location 10km in the middle of the river, this point is 5km from station A and 5km from station B, with symmetrical influence between the two stations. If the state deviation extracted by station A is +0.5m and that of station B is -0.3m, the denominator ≈ 1.782 > 1, and the gain of each station is compressed to K. A ≈0.50, K B≈0.50, the constraint correction at this point = 0.50 × (+0.5) + 0.50 × (-0.3) = +0.10m. At station A's own location, the denominator = max(1, 1 + exp(-10000 / 43200)) = max(1, 1.793) = 1.793, K A =1 / 1.793≈0.56. The above results show that the bi-station correction in the overlapping area is adaptively compressed, avoiding the oscillations that may be caused by the simple linear superposition of the positive and negative corrections of the two stations.

[0159] Step 604: The constraint correction is superimposed on the current free calculation state of the cascaded coupled model to obtain the updated initial state.

[0160] After obtaining the constraint correction amount distributed throughout the space, the correction amount matrix is ​​directly summed and superimposed onto the current flow field water level state vector of the hydrodynamic module in the cascaded coupling model to correct the accumulated error of historical calculations and output the updated initial state.

[0161] As another complex alternative to this step, under conditions where high-frequency synchronous observation of water level and flow rate is available in the basin, the above-mentioned single-point update logic based on physical wave velocity can be upgraded to an ensemble Kalman filter data assimilation framework. By introducing disturbance terms into the runoff generation and hydrodynamic models to construct a set of sampled model states, the empirical covariance matrix of the sample set is used to dynamically replace the above-mentioned fixed spatialized gain function, thereby achieving adaptive constraint updates of multiple state variables.

[0162] Example 7

[0163] Based on the method concept of Embodiment 1 above, this paper further introduces the process of constructing and transmitting mixed rainfall-driven sequences with a preset period for advancing forecast times and rolling execution, as well as the engineering fault-tolerant operation mechanism in the scenario of missing hydrological station data.

[0164] An alternative approach is to advance the forecast time with a preset period and execute the construction and propagation process of the mixed rainfall-driven sequence in a rolling manner.

[0165] The entire forecasting process is a dynamically advancing closed-loop structure in terms of time sequence. A preset fixed period (i.e., a preset rolling period) is loaded as the step size for time progression. Whenever the actual time passes through this fixed period, the current forecast time is moved forward by one step size.

[0166] In practice, the preset fixed period can be set to 6 hours, but the specific time can be determined based on the actual runoff time in the watershed. As the current forecast time advances, the historical time period boundary is simultaneously shifted 6 hours into the future. Correspondingly, the adaptive time window required for rainfall assimilation also slides to the right, reloading the newly added measured rainfall data from rain gauges to refresh the one-step lag autocorrelation coefficient and time decay weight, and reconstructing the mixed rainfall driving sequence. After receiving the refreshed driving field, the cascaded coupled model initiates a new round of forward extrapolation, outputting the spatiotemporal processes of the cross-section along the route within the future forecast period.

[0167] An alternative approach is to pause the spatialized gain mechanism update operation if the measured water level and flow data of the hydrological station cannot be successfully obtained at the current forecast time or if the obtained data is abnormally missing during the rolling execution of the propagation simulation. Instead, the free calculation state of the cascaded coupling model at the previous forecast time is used as the initial state for continuous simulation to maintain the temporal continuity of the forecast process.

[0168] Hydrological stations in small and medium-sized river basins are typically located in the field. During extreme weather events such as torrential rains and floods, they may experience missing data or faulty points due to communication interruptions, power outages, or sensor drift. This step provides a fault-tolerant mechanism by establishing conditional flow branching to ensure the system's continuous forecasting capability under adverse conditions.

[0169] In the first operating condition branch, if the data access unit successfully obtains the measured water level or flow data of the hydrological station that meets the quality control specifications at the current forecast time, the spatialized gain algorithm described in Example 6 is triggered to perform the assimilation correction operation of the initial state of the flow field.

[0170] In the second operating condition branch, if anomalies are detected in the measured data of the third-line hydrological station network at the current forecast time, such as missing measurements, unresponsiveness due to timeouts, or values ​​exceeding the physical reasonableness threshold, the condition judgment mechanism intercepts the state update command. The calculation of the spatialized gain mechanism is automatically paused to prevent numerical oscillations caused by loading bad data. At this time, the system switches to the free calculation fault-tolerant path, directly extracting the calculated flow field state of the cascaded coupled model at the previous forecast time for the current time. This calculated state is used as a seamless initial condition, substituted into the control equations of the hydrodynamic module, directly driving the model forward. This mechanism reduces the risk of interruption in the flood forecast of the entire basin due to the failure of individual station sensors, ensuring the robustness of long-duration rolling simulations.

[0171] Example 8

[0172] Based on Embodiment 1 above, the application of flood control early warning and uncertainty assessment after outputting the spatiotemporal forecast results of water level and flow rate for the target river section within the forecast period is further described in detail. In one possible implementation, this step includes the following steps:

[0173] Step S1: Based on the spatiotemporal forecast results of water level and flow rate, extract the flood arrival time characteristics, peak water level characteristics, and peak flow rate characteristics of the preset target inundation risk point sections.

[0174] After the cascaded coupled model completes one round of simulation, the output basic data is a two-dimensional spatiotemporal matrix covering all computational grids across the entire river section. To transform the underlying continuous hydrodynamic state variables into decision indicators directly supporting flood control operations, the grid index of the target inundation risk point section is mapped and located based on pre-entered geospatial coordinates. For this predetermined grid, a numerical extreme value search operation is performed on the time axis of the forecast period. By traversing the water level time series of this section, the maximum value and its corresponding timestamp are extracted and determined as the peak flood level characteristic and flood arrival time characteristic, respectively. Similarly, the maximum value is searched in the corresponding flow time series to determine the peak flow characteristic. This extraction mechanism discretizes the continuous flow field simulation process into specific risk parameters.

[0175] As an optional implementation method, if there are multiple protection objects of different levels along the target river section, multiple sets of coordinate indices can be loaded in parallel, and the arrival time and extreme value characteristic sequences of each section along the river can be output in batches.

[0176] Step S2 involves performing ensemble statistical analysis on the rolling peak flow characteristics over multiple forecast periods, extracting the convergence interval of the peak flow characteristics as the rolling progresses, and combining the flood arrival time characteristics and peak water level characteristics to generate and output flood warning information containing forecast uncertainties.

[0177] Considering that the output of a single forecast period is constrained by the uncertainty of input boundary errors and initial states, a time-based ensemble statistical evaluation mechanism is established. As the forecast time rolls forward at fixed intervals, for the same future potential flood peak occurrence time, the flood peak flow characteristic values ​​generated in each rolling simulation are extracted and recorded, thereby constructing a forecast ensemble sample sequence that continuously expands with the increase of rolling batches.

[0178] After obtaining the forecast set sample sequence, the mean, variance, and probability density distribution of the sequence are calculated. Based on the set confidence parameters, the upper and lower limits of the values ​​under this distribution are extracted and used as the convergence interval for the peak flow characteristics as the forecast progresses. In the actual rolling assimilation process, as the forecast time continuously approaches the actual physical time of the flood peak, high-precision measured rainfall from rain gauges and observation data from hydrological stations are continuously absorbed, driving the error inward correction, resulting in a gradually narrowing width of the extracted convergence interval. Based on the extracted extreme value features and the convergence interval, flood warning information is generated and sent to the terminal interface. The above output method, which includes the probability distribution range, changes the traditional single-point value forecast format and provides a flood control scheduling basis with uncertainty quantification reference.

[0179] Example 9

[0180] Based on Embodiment 1 above, the mechanism for evaluating the effectiveness of system forecast information and determining the maximum effective lead time is further described in detail. In one possible implementation, the steps may include, but are not limited to, the following:

[0181] Step S01: For each forecast period, calculate the reference variance when there is no forecast information and the forecast error variance of each data source of the three lines of defense based on the pre-stored historical observation statistics. Then, evaluate the effectiveness of the forecast information of meteorological forecast rainfall data, rain gauge measured rainfall data and hydrological station measured water level and flow data.

[0182] When evaluating the data quality and forecast boundaries of flood control systems, dimensionless information validity is introduced as a quantitative evaluation index. For the target hydrological variable, historical long-series observation data are acquired, and the climatological natural statistical variance under conditions without forecast information access is calculated and determined as the reference variance. This parameter is equivalent to the statistical variance generated when only the historical statistical mean is used as the baseline forecast result. For a set lead time axis, the discrete distribution characteristics between the simulated output and the measured true value corresponding to each flood control line data source are statistically analyzed to generate the forecast error variance corresponding to each data source.

[0183] After obtaining the above variance parameters, the effectiveness of the forecast information of a single data source at the forecast period is calculated according to the preset mapping relationship, as follows:

[0184] e i (τ)=σ ref 2 / (σ ref 2 +σ i 2 (τ));

[0185] Among them, e i (τ) represents the effectiveness of the forecast information from the i-th data source at the lead time τ, σref 2 σ is the reference variance calculated based on historical data. i 2 (τ) represents the forecast error variance of the i-th data source at the forecast period τ.

[0186] This calculation process maps the error characteristics of various heterogeneous data sources to a set dimensionless validity range, providing a numerical basis for subsequent multi-source information integration and evaluation.

[0187] Step S02: Calculate the global forecast validity of the composite system based on the forecast information validity of each data source, compare it (global forecast validity) with the preset forecast validity threshold, determine and output the current maximum effective forecast period of the system.

[0188] Because the main effective time intervals of the acquired meteorological forecast rainfall data, rain gauge measured rainfall data, and hydrological station measured water level and flow data differ along the time axis, and because the physical mechanisms generating errors exhibit decoupled properties, the data processing module treats them as independent observation sources at the statistical level. Based on the optimal linear unbiased estimation theoretical framework, the forecast error variance of each data source is extracted and nonlinear algebraic operations are performed to calculate the overall forecast error variance of the composite system driven by fused information. The overall forecast error variance is calculated as follows:

[0189] σ sys 2 (τ)=1 / (1 / σ1 2 (τ)+1 / σ2 2 (τ)+1 / σ3 2 (τ));

[0190] Where, σ sys 2 (τ) represents the overall forecast error variance of the composite system at the lead time τ, σ1 2 (τ) represents the forecast error variance corresponding to the meteorological forecast rainfall data, σ² 2 (τ) represents the forecast error variance corresponding to the measured rainfall data from the rain gauge, σ3 2 (τ) represents the forecast error variance corresponding to the measured water level and flow rate data of the hydrological station.

[0191] After obtaining the overall forecast error variance, it is substituted into the effectiveness evaluation framework to generate the global forecast effectiveness, which characterizes the overall performance of the system.

[0192] e sys (τ)=σ ref 2 / (σ ref 2 +σ sys 2 (τ));

[0193] Among them, e sys (τ) represents the global forecast effectiveness of the composite system at the lead time τ, σ ref 2 For reference variance, σ sys 2 (τ) represents the calculated variance of the overall forecast error.

[0194] After generating a global forecast validity sequence covering the forecast period time axis, a pre-set forecast validity threshold is obtained. A comparison mechanism performs a point-by-point numerical scan along the future direction of the time axis, searching for and recording the furthest time span node that makes the global forecast validity greater than or equal to the forecast validity threshold. This time span node is extracted and established as the system's current maximum effective forecast period. This mechanism, from a numerical analysis perspective, demonstrates and quantifies that a multi-source information coupling mechanism can overcome the theoretical limits of a single-source forecast architecture in terms of time scale.

[0195] In a specific verification scenario, historical rainstorm and flood events in a small to medium-sized river basin were selected as test samples. The method of this invention was compared with traditional forecasting methods driven by only a single data source. The Nash efficiency coefficient of the method of this invention for peak flow forecasting within a 6-24 hour lead time is significantly improved compared to schemes driven solely by measured rainfall from rain gauges. Simultaneously, compared to schemes driven directly by uncorrected meteorological forecast rainfall, the error divergence is effectively suppressed in lead times exceeding 12 hours. Those skilled in the art will understand that the specific accuracy improvement may vary depending on the characteristics of the basin, the type of rainfall event, and the configuration of model parameters.

[0196] This application employs an adaptive time window based on the lag autocorrelation coefficient of the error sequence and a hybrid truncation mechanism of relative and absolute deviations, which effectively suppresses the abnormal amplification of the product term when the base of sudden rainfall is extremely small or extremely large. It suppresses the systematic error of the rainfall driving field from a divergent and distorted state to a stable and convergent state, thus extending the effective prediction period while ensuring the accuracy of the driving information.

[0197] A morphologically adaptive cross-section extraction method based on thinning threshold and forced spacing constraint of straight river sections is adopted. This method not only preserves the geometric features of sharp bends but also fills the node gaps in straight river sections, avoiding hydrodynamic solution divergence caused by excessively large spatial grid step size. This provides a spatial topological foundation with high numerical computational stability for watersheds without data.

[0198] By introducing the physical wave velocity of flood as the spatial correlation scale and combining it with the spatial gain update mechanism of normalizing the denominator of the multi-station distance exponent, the natural spatial exponential decay of the correction amount is maintained in the single-station control area, while the excessive superimposed correction weight is adaptively compressed in the multi-station dense area. This effectively suppresses numerical overflow and oscillation in the state assimilation process and realizes the smooth reconstruction of the initial state of the flow field in the entire basin.

[0199] By adopting a state substitution fault tolerance mechanism based on the integrity condition judgment of hydrological data, the risk of invalid bad points contaminating the flow field is reduced, and the continuity and overall robustness of the system flood evolution sequence are ensured under extremely harsh operating conditions.

[0200] It should be noted that the various specific technical features described in the above embodiments can be combined in any suitable manner without contradiction. To avoid unnecessary repetition, the present invention will not describe the various possible combinations separately.

Claims

1. A flood propagation forecasting method coupling multi-source information from three lines of defense, characterized in that, include: Acquire meteorological forecast rainfall data, actual rainfall data from rain gauges, and actual water level and flow data from hydrological stations; The meteorological forecast rainfall data and the actual rainfall data measured by rain gauges are normalized in terms of time and space. Based on the rainfall deviation characteristics in historical periods, a revised forecast rainfall is constructed. The actual rainfall data measured by rain gauges, the revised forecast rainfall, and the meteorological forecast rainfall data are spliced ​​together according to the time partition structure to construct a hybrid rainfall driving sequence. The target river channel midline is extracted using pre-acquired digital elevation data for feature point identification and constraint point supplementation. Based on the supplemented set of control points, river channel cross-section data is generated, and a hydrodynamic topographic basic grid is constructed. The mixed rainfall-driven sequence is input into a cascaded coupled model based on hydrodynamic topography grid configuration for propagation and simulation. During the simulation, measured water level and flow data from hydrological stations are used to update the initial state of the hydrodynamic module in the cascaded coupled model through a spatialized gain mechanism. The forecast time is advanced by a preset period and the process of constructing and extrapolating the mixed rainfall-driven sequence is executed in a rolling manner, outputting the spatiotemporal forecast results of water level and flow rate of the target river segment within the forecast period; Using measured water level and flow rate data from hydrological stations, the initial state of the hydrodynamic module in the cascaded coupled model is updated through a spatialized gain mechanism, including: When the measured water level and flow rate data of the hydrological station are received, the state deviation is extracted by comparing it with the free calculation state of the corresponding section of the cascaded coupling model. Construct a spatialized gain function that decays exponentially with increasing spatial distance; The spatialized gain function is used to propagate the state deviation along the river channel upstream and downstream, generating a spatially distributed constraint correction. The constraint correction is superimposed on the current free calculation state of the cascaded coupled model to obtain the updated initial state; The spatialized gain function is expressed as follows: ； Among them, K j (s) represents the spatialized gain function value of the j-th hydrological station at any cross-sectional position s on the river channel, K0 is the preset local gain coefficient, exp is an exponential function with the natural constant as the base, and |ss| obs,j | represents the distance s along the path between any cross-sectional location s and the j-th hydrological station. obs,j The absolute distance between them, L s The calculated spatial correlation scale is given by max, which is the maximum value operator. The summation symbol is |ss| for iterating over all M hydrological stations within the basin. obs,l | represents the distance s between any cross-sectional location s and the l-th hydrological station. obs,l The absolute distance between them.

2. The method according to claim 1, characterized in that, Based on historical rainfall deviation characteristics, revised rainfall forecasts are constructed, including: At the current forecast time, based on the trend of the deviation between the measured rainfall data from the rain gauge and the forecast rainfall data, an adaptive time window for backtracking historical periods is determined; Within the adaptive time window, the relative deviation factor and absolute deviation of each spatial grid point are extracted time by time. The extracted relative deviation factor and absolute deviation are weighted and averaged using a time decay mechanism to obtain the weighted relative deviation and weighted absolute deviation. Based on weighted relative bias and weighted absolute bias, a hybrid correction operation is performed on the forecast rainfall within the future time window to output the corrected forecast rainfall.

3. The method according to claim 2, characterized in that: Based on the trend of the deviation between measured rainfall data from rain gauges and forecasted rainfall data, an adaptive time window for backtracking historical periods is determined, including: Calculate the spatial average of the absolute deviation of all spatial grid points at each time point, and construct the surface average absolute deviation sequence. The one-step lag autocorrelation coefficient of the surface mean absolute deviation sequence is extracted, and its value is used to dynamically interpolate between the preset lower and upper limits of the window to determine the length of the adaptive time window. The extracted relative deviation factor and absolute deviation are weighted and averaged using a time decay mechanism, including: Construct a time weighting function that decays exponentially with increasing time distance; By using a time weighting function to assign higher weights to deviation samples that are closer to the current forecast time, the relative deviation factor and the absolute deviation are weighted and summed to obtain the weighted relative deviation and the weighted absolute deviation.

4. The method according to claim 2, characterized in that, Based on weighted relative bias and weighted absolute bias, a hybrid correction operation is performed on the forecast rainfall within the future time window, outputting the corrected forecast rainfall, including: For the predicted rainfall amount of the forecasted rainfall field within the future time window, the weighted relative deviation is multiplied by the predicted rainfall amount to obtain the relative correction amount; Compare the absolute value of the relative correction with the absolute value of the weighted absolute deviation; When the absolute value of the relative correction is less than the absolute value of the weighted absolute deviation, the relative correction is used as the target correction; otherwise, the absolute value of the weighted absolute deviation after taking the sign of the relative correction is used as the target correction. The target correction is superimposed on the forecast rainfall and then subjected to non-negative truncation to obtain the corrected forecast rainfall.

5. The method according to claim 1, characterized in that, The mixed rainfall-driven sequence is input into a cascaded coupled model based on hydrodynamic topography grid configuration for propagation and derivation, including: The mixed rainfall-driven sequence is input into the runoff generation and confluence model for simulation, and the boundary conditions of each tributary side and the inflow data of the reservoir within the target watershed are calculated and output. The inflow data output from the runoff generation and confluence model is input into the reservoir scheduling model for calculation, and the upper boundary conditions are calculated and output. The upper boundary conditions and the lateral boundary conditions of each tributary are input into the hydrodynamic simulation model, and the hydrodynamic topographic base grid is used as the topographic boundary to continuously simulate and calculate the propagation process of flood in the river channel in order to solve the flow state along the route.

6. The method according to claim 1, characterized in that, During the rolling execution of the propagation simulation, if the measured water level and flow data of the hydrological station cannot be successfully obtained at the current forecast time or if the obtained data is abnormally missing, the update operation of the spatialized gain mechanism will be paused, and the free calculation state of the cascaded coupling model at the previous forecast time will be used as the initial state for continuous simulation to maintain the temporal continuity of the forecast process.

7. The method according to claim 1, characterized in that, After outputting the spatiotemporal forecast results of water level and flow rate for the target river segment within the forecast period, it also includes: Based on the spatiotemporal forecast results of water level and flow rate, the characteristics of flood arrival time, peak water level and peak flow rate of the target inundation risk point section are extracted; By performing ensemble statistical analysis on the rolling peak flow characteristics over multiple forecast periods, the convergence interval of the peak flow characteristics as the rolling progress is extracted, and flood control early warning information containing forecast uncertainties is generated and output.

8. The method according to claim 1, characterized in that, Also includes: For each forecast period, the reference variance and forecast error variance of each data source in the three lines of defense are calculated when there is no forecast information, and the effectiveness of forecast information based on meteorological forecast rainfall data, rain gauge measured rainfall data, and hydrological station measured water level and flow data is evaluated. The global forecast validity of the composite system is calculated based on the forecast validity of each data source, and compared with a pre-set forecast validity threshold to determine the current maximum effective forecast period of the system.

9. A flood propagation forecasting system coupling multi-source information from three lines of defense, characterized in that, include: At least one processor; as well as, A memory that is communicatively connected to at least one processor; wherein, The memory stores instructions executable by a processor to implement the method according to any one of claims 1 to 8.

Images