Optimization Methods for Storage and Transportation of Flood and Drought Disaster Prevention Materials

Abstract

This invention discloses an optimization method for the storage and transportation of flood and drought disaster prevention materials, comprising: acquiring real-time hydrological data of the monitoring area and calculating the river-lake hydrodynamic coupling index reflecting the interaction between main streams and tributaries; performing time-series analysis on the index to identify the current river-lake hydrodynamic state and assessing disaster risk accordingly; constructing a dynamic water level field by combining hydrological connectivity analysis and extracting the dynamic inundation range; constructing a time-varying transport network accessibility matrix based on the topological relationship between the inundation range and road network elevation; and generating material scheduling and route planning schemes by comprehensively considering the material demand correction coefficients under different hydrodynamic states and the dynamic blockage of the road network. This invention introduces physical mechanism-driven state recognition and connectivity interpolation, solving the problems of inaccurate prediction of complex water system disasters and scheduling deviating from physical reality, and realizing accurate emergency response and dynamic closed-loop.

Description

Technical Field

[0001] This invention relates to a method for optimizing the storage and transportation of materials, and more particularly to a method for optimizing the storage and transportation of materials for flood and drought disaster prevention. Background Technology

[0002] In large river systems such as the middle and lower reaches of the Yangtze River, where rivers and lakes are interconnected, flood and drought prevention relies on the assessment of water conditions and the timely allocation of resources. Accurately identifying the hydrodynamic interactions (such as backflow and flooding) between the main stream and connected lakes, and assessing the risk of inundation and the need for resources accordingly, helps to ensure flood control and drought relief safety.

[0003] Currently, disaster relief material dispatching methods mainly rely on real-time water level data from hydrological monitoring stations and static path planning based on Geographic Information Systems (GIS). Existing technical solutions typically trigger early warnings directly by detecting water level exceeding limits at monitoring stations, and use inverse distance weighting (IDW) based on Euclidean distance or the Thiessen polygon method to estimate water level distribution in areas without data. In calculating material demand, allocation is often based on a fixed proportion of the population within an administrative region. In terms of transportation dispatching, vehicle assignment is mostly done using shortest path algorithms based on fixed road network topology.

[0004] However, existing technologies suffer from insufficient accuracy and feasibility when handling complex river and lake interaction scenarios, including simplistic water level assessment leading to misjudgments of disaster status, inaccurate water level reconstruction due to a lack of consideration for engineering barriers, and failure to integrate the dynamic impact of floods into emergency dispatching. Therefore, further research and innovation are needed to address these issues in existing technologies. Summary of the Invention

[0005] Purpose of the invention: This application provides an optimized method for the storage and transportation of flood and drought disaster prevention materials to solve the aforementioned problems.

[0006] Technical solution: According to one aspect of this application, an optimized method for the storage and transportation of flood and drought disaster prevention materials includes:

[0007] Acquire real-time hydrological data and basic geographic data of the monitoring area, and calculate the river-lake hydrodynamic coupling index reflecting the interaction between the main stream and tributaries based on the real-time hydrological data;

[0008] Time-series characteristic analysis of the river-lake hydrodynamic coupling index is performed to identify the current river-lake hydrodynamic state in the monitoring area;

[0009] Disaster risk classification is based on the hydrodynamic state of rivers and lakes, dynamic water level field is constructed by combining hydrological connectivity analysis, and dynamic inundation range and spatiotemporal distribution of material demand in disaster-stricken areas are determined based on the dynamic water level field.

[0010] Based on the spatial topological relationship between the dynamic inundation range and the road elevation information in the basic geographic data, a transportation network accessibility matrix reflecting the real-time traffic capacity of the road network is constructed.

[0011] Using the spatiotemporal distribution of material demand as input and the accessibility matrix of the transportation network as constraint, a material scheduling and route planning scheme is generated.

[0012] Beneficial effects: This invention introduces physical mechanism-driven state recognition and connectivity interpolation, solving the problems of inaccurate prediction of complex water system disasters, distorted water level reconstruction, and scheduling deviating from physical reality, thus achieving accurate emergency response and dynamic closed-loop. The related technical effects will be described in detail below with reference to specific embodiments. Attached Figure Description

[0013] Figure 1 A flowchart illustrating an optimized method for the storage and transportation of flood and drought disaster prevention materials, provided as an embodiment of this application.

[0014] Figure 2 A flowchart for calculating the river-lake hydrodynamic coupling index, which reflects the interaction between tributaries, is provided for an embodiment of this application.

[0015] Figure 3 This is a flowchart illustrating disaster risk classification based on the hydrodynamic state of rivers and lakes, provided as an embodiment of this application.

[0016] Figure 4 A flowchart for constructing a dynamic water level field by combining hydrological connectivity analysis, provided for embodiments of this application.

[0017] Figure 5 A flowchart illustrating the introduction of a correction mechanism based on the hydrodynamic state of rivers and lakes, provided for embodiments of this application. Detailed Implementation

[0018] 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.

[0019] It should be noted that the terms "first," "second," etc., in the specification and accompanying drawings of this invention are used to distinguish similar objects and are not necessarily used to describe a specific order or sequence. It should be understood that such data can be interchanged where appropriate so that embodiments of the invention described herein can be implemented in orders other than those illustrated or described herein. Furthermore, the terms "including" and "having," and any variations thereof, are intended to cover non-exclusive inclusion; for example, a process, method, system, product, or apparatus that includes a series of steps or units is not necessarily limited to those steps or units explicitly listed, but may include other steps or units not explicitly listed or inherent to such processes, methods, products, or apparatus.

[0020] To address the aforementioned issues, the applicant conducted in-depth searches and analyses, and discovered:

[0021] Correspondingly, relying solely on a single water level threshold cannot distinguish between fundamentally different physical states such as high-water-level outflow and high-water-level backflow, leading to misjudgments of the nature and evolution trend of disasters. Furthermore, interpolation methods based on Euclidean distance ignore the obstructive effect of water conservancy projects such as dikes and sluices on water flow paths, often resulting in distorted reconstruction of the water level field.

[0022] Based on this, static road network planning does not take into account dynamic road blockages caused by flooding, and the forecast of material demand fails to reflect the different impacts of different water flow patterns, such as still water soaking and rapid current impact, on material consumption, making it difficult to achieve accurate and feasible emergency dispatch.

[0023] To solve these problems, combined with Figures 1 to 5 The present invention will be specifically described through the following embodiments.

[0024] On the one hand, an exemplary scheme for optimizing the storage and transportation of flood and drought disaster prevention materials is provided. This not only demonstrates the steps from multi-source data acquisition to scheduling decision generation, but also illustrates the data flow and logical dependencies between each stage, verifying the technical feasibility of the scheme in complex flood and drought disaster scenarios. Accordingly, this embodiment can be implemented through the following steps:

[0025] Step 101: Obtain real-time hydrological data and basic geographic data of the monitoring area, and calculate the river-lake hydrodynamic coupling index reflecting the interaction between the main stream and tributaries based on the real-time hydrological data.

[0026] The acquired real-time hydrological data specifically includes physical quantities such as water level, flow rate, and hydraulic distance between key control stations. This data typically originates from hydrological monitoring stations or remote sensing inversion systems deployed within the watershed. Basic geographic data includes high-precision digital elevation models (DEMs), road network topology, and historical hydrological data. Furthermore, this step introduces the river-lake hydrodynamic coupling index, a dimensionless physical index that does not process water level or flow rate data through linear superposition but rather uses a specific physical mechanism to nonlinearly couple factors reflecting interactions such as backflow and flooding.

[0027] Alternatively, by weighting and integrating the normalized hydraulic gradient component and flow throughput component, the potential energy factor and kinetic energy factor reflecting the interaction forces such as backflow and inflow can be comprehensively characterized into a unified index, namely the river-lake hydrodynamic coupling index.

[0028] Specifically, calculating the river-lake hydrodynamic coupling index involves normalizing and weighting the hydraulic gradient and flow ratio. This index compresses the fluid dynamic state into measurable time-series values, providing a physical basis for subsequent state identification. For example, in the interaction between the Yangtze River and its connected lakes (such as Poyang Lake), this index can capture the critical point where natural flow transitions into backflow.

[0029] Step 102 involves performing time-series characteristic analysis on the river-lake hydrodynamic coupling index to identify the current river-lake hydrodynamic state in the monitoring area. Alternatively, the river-lake hydrodynamic state can also be referred to as the hydrodynamic condition.

[0030] Furthermore, the hydrodynamic state of rivers and lakes refers to the discrete physical patterns classified based on the temporal evolution characteristics of the coupling index, such as a strong outflow state, a weak backflow state, or a backflow state. The identification process uses time-series analysis algorithms, such as cluster analysis or probabilistic models, to capture the changing trends and transformation patterns of the index over time. Identification methods based on time-series characteristics can distinguish between instantaneous numerical fluctuations and actual physical state transformations, avoiding misjudgments caused by data noise.

[0031] For example, when the coupling index exhibits a specific pattern over a period of time, such as sustained high-level oscillation or a sharp decline, the system identifies it as the corresponding river and lake hydrodynamic state. The identified state not only reflects the current water conditions but also implies future evolution trends, providing a label for subsequent risk assessment and demand forecasting.

[0032] Step 103: Based on the hydrodynamic state of the river and lake, conduct disaster risk classification, construct a dynamic water level field by combining hydrological connectivity analysis, and determine the dynamic inundation range and the spatiotemporal distribution of material needs in the disaster-stricken area based on the dynamic water level field.

[0033] In other words, disaster risk classification is performed based on the hydrodynamic state of rivers and lakes to obtain the disaster risk level; a dynamic water level field is constructed based on the water connectivity information in the basic geographic data and combined with hydrological connectivity analysis; the dynamic inundation range and inundation depth are determined according to the dynamic water level field and the digital elevation model in the basic geographic data; and the spatiotemporal distribution of material demand in the disaster-stricken area is determined based on the dynamic inundation range, inundation depth and disaster risk level.

[0034] Accordingly, based on the identified hydrodynamic state of the rivers and lakes, the corresponding risk assessment model is invoked to determine the current disaster risk level. Under different conditions, such as backflow and leakage, the same water level may correspond to drastically different risk levels. Furthermore, hydrological connectivity analysis is introduced to estimate the water level in unmonitored areas based on the actual physical path of the water flow, i.e., the connectivity distance.

[0035] Accordingly, the dynamic inundation range can be extracted based on the overlay analysis of the reconstructed dynamic water level field and the digital elevation model (DEM). Furthermore, by combining the population distribution, asset density, and risk level of the affected area, the material needs of each affected location can be calculated. As shown above, material needs exhibit spatiotemporal distribution characteristics; that is, the demand for materials (such as flood control sandbags, drainage pumps, or daily necessities) at different locations at different times is dynamically changing and directly influenced by the hydrodynamic conditions of rivers and lakes (such as whether backflow prevention is required).

[0036] Step 104: Based on the spatial topological relationship between the dynamic flooding range and the road elevation information in the basic geographic data, construct a transportation network accessibility matrix that reflects the real-time traffic capacity of the road network.

[0037] In this step, the accessibility of the transportation network is defined as a variable constrained by the dynamic inundation range. Specifically, the system maps the predicted dynamic inundation range onto the road network and compares the ground elevation of road segments with the predicted water level. When the water depth of a road segment exceeds the safe passage threshold for vehicles, that road segment is marked as blocked in the topology network, and its passage cost in the adjacency matrix is ​​set to infinity. Based on this, the constructed transportation network accessibility matrix is ​​a dynamically updated matrix over time, reflecting the road network survival status under disaster scenarios.

[0038] It should be understood that the dynamic update mechanism ensures that the subsequently generated route planning schemes will not guide transport vehicles into dangerous areas that have been flooded or are about to be flooded, thus guaranteeing the feasibility of scheduling.

[0039] Step 105: Using the spatiotemporal distribution of material demand as input and the transportation network accessibility matrix as constraint, generate a material scheduling and route planning scheme.

[0040] The system aims to meet the spatiotemporally distributed material needs, such as maximizing coverage and minimizing response time, while strictly adhering to the physical constraints defined by the reachability matrix. The generated solutions specifically include: the material outflow volume for each reserve warehouse, the travel route sequence of transport vehicles, and the delivery time windows for each disaster-stricken area.

[0041] For example, in one application scenario, when the system identifies backflow and predicts an impending closure of a critical road, the generated scheduling plan prioritizes the deployment of drainage equipment and automatically plans a detour around the low-lying section, ensuring that supplies reach the high-risk area before the road is blocked. It should be understood that this process achieves an end-to-end closed loop from water situation monitoring to logistics response, improving the scientific rigor and timeliness of emergency dispatch.

[0042] On the other hand, the optional implementation method for the physical definition and calculation of the river-lake hydrodynamic coupling index solves the technical problem that traditional methods, relying solely on a single water level or flow rate index, cannot accurately describe the nonlinear interactions (such as backflow and flooding) between the main stream and tributaries (or connected lakes), thus providing a physical benchmark for subsequent state identification. Furthermore, this embodiment specifically includes:

[0043] Step 201: Based on the water levels of the main stream control station, the tributary control station, and the hydraulic distance between the two stations in the real-time hydrological data, calculate the normalized hydraulic gradient component, wherein the normalization process is performed using the reference hydraulic gradient.

[0044] In this application, the hydraulic gradient component is used to quantify the potential energy gradient driving water flow. Accordingly, normalization can be performed to eliminate the dimensional effects caused by varying monitoring station spacing, making the indicator comparable across different watersheds. The specific calculation formula is as follows:

[0045] G _norm (t)=(H _main (t)-H _sub (t)) / (ΔL×i _ref );

[0046] Among them, G _norm (t) represents the normalized hydraulic gradient component at time t, which is dimensionless; H _main (t) represents the water level at a main stream control station (e.g., Hankou Station on the Yangtze River) at time t, in meters; H _sub (t) represents the water level at time t at the outlet control station of a tributary or lake (e.g., Xingzi station of Poyang Lake), in meters; ΔL represents the hydraulic distance between the main stream control station and the tributary control station, in kilometers; i _ref This is the preset reference hydraulic gradient, in meters per kilometer.

[0047] Specifically, the hydraulic gradient is typically taken as the natural water level gradient of the river section during the multi-year average dry season, with a value ranging from 0.01 to 0.05 meters per kilometer. The commonly used unit for hydraulic gradient, meters per kilometer, is numerically equal to a few thousandths (‰). For example, 0.02 meters per kilometer means that the water level drops by 0.02 meters per kilometer, which can also be written as 0.02‰. When G... _norm When (t) is positive and large, it indicates that the water level of the main stream is higher than that of the tributaries, and there is a tendency for backflow or even reverse flow; when G _norm When (t) is negative, it indicates that the tributary is draining smoothly. This is further explained by introducing ΔL and i. _ref This formula transforms the water level difference into a multiple relationship relative to the characteristic gradient, characterizing the driving force of gravitational potential energy on the flow.

[0048] Step 202: Calculate the normalized flow throughput component based on the outlet cross-sectional flow and reference flow in the real-time hydrological data.

[0049] The throughput component reflects the intensity and direction of water transport. Accordingly, using a reference flow rate for dimensionless calculation eliminates the influence of differences in absolute flow rate magnitudes, such as the significant difference between flow rates during dry and flood seasons. The calculation formula is as follows:

[0050] Q _norm (t)=Q _out (t) / Q _ref ;

[0051] Among them, Q _norm (t) represents the normalized throughput component at time t, which is dimensionless; Q _out (t) represents the measured flow rate at the outlet section of the tributary or lake at time t, expressed in cubic meters per second. Generally, outflow (towards the main stream) is considered positive, and backflow (towards the tributary) is considered negative. _ref This is the preset reference flow rate, in cubic meters per second.

[0052] For example, the reference flow rate can be set as the characteristic discharge rate of the outlet section at the warning water level or the multi-year average maximum flow rate. Through this process, Q... _norm The absolute value of (t) intuitively reflects the proportion of the current flow rate relative to the full load state of the cross section, and can distinguish the essential difference in hydrodynamic energy between small flow backflow and large flow backflow.

[0053] Step 203: Use weighting coefficients to perform a weighted summation of the hydraulic gradient component and the flow throughput component to obtain the river-lake hydrodynamic coupling index.

[0054] In this step, the river-lake hydrodynamic coupling index is a comprehensive representation of the potential energy term (hydraulic gradient) and the kinetic energy term (flow rate), and the calculation formula is as follows:

[0055] Φ(t)=α×G_norm (t)+β×Q _norm (t);

[0056] Where Φ(t) is the river-lake hydrodynamic coupling index at time t; α is the weighting coefficient of the hydraulic gradient component; β is the weighting coefficient of the flow throughput component; and α+β=1. Specifically, the values ​​of the weighting coefficients α and β depend on the hydrological characteristics of the monitoring area. For river-type tributaries dominated by gravity flow, hydraulic gradient plays a dominant role, and α can be taken as a larger value, such as 0.6-0.8; for the outlet of large lakes that are greatly affected by wind, waves, or inertia, the flow term better reflects the instantaneous state, and β can be appropriately increased.

[0057] For example, suppose the measured data at a certain moment is as follows: Main stream water level H _main =20.0 meters, tributary water level H _sub =19.8 meters, the distance between the two stations ΔL=50 kilometers, reference slope i _ref =0.02 m / km; Outlet flow rate Q _out = -500 cubic meters per second, a negative value indicates backflow, reference flow rate Q _ref =10000 cubic meters / second. Assume weights α=0.6, β=0.4.

[0058] Accordingly, the hydraulic gradient component G is calculated. _norm =(H _main -H _sub ) / (ΔL×i _ref ) = (20.0 - 19.8) / 1.0 = 0.2; where, H _main For the water level at the main stream control station, H _sub The water level at the tributary control station is ΔL, where ΔL is the hydraulic distance between the two stations, and i _ref The reference hydraulic gradient is used. This result indicates that the currently measured hydraulic gradient is 0.2 times the reference gradient.

[0059] Furthermore, the flow throughput component Q is calculated. _norm =-500 / 10000=-0.05. Calculate the coupling index Φ=0.6×0.2+0.4×(-0.05)=0.12-0.02=0.1. The result Φ=0.1 is positive, indicating that although backflow has occurred and the flow rate is negative, the overall hydrodynamic situation is still dominated by the backflow of the main stream (positive slope), which is the initial stage of weak backflow.

[0060] Alternatively, the weighted summation process can be replaced with a nonlinear fusion function, such as introducing a sign-preserving squared term in the component calculation, or dynamically adjusting the weighting coefficients α and β according to the season. For example, during the main flood season, the weight of α can be dynamically increased to give more emphasis to the impact of the flood warning water level difference.

[0061] According to one aspect of this application, the reference hydraulic gradient can be determined using a historical statistical method, with the specific steps as follows:

[0062] Collect daily water level monitoring data of the target river section over the past N years, and select water level records for the normal water period each year, usually from March to May or from October to November, where N is not less than 10.

[0063] For each day's data, calculate the measured hydraulic gradient i between the main stream control station and the tributary control station. _actual (t)=(H _main (t)-H _sub (t)) / ΔL; where, i _actual (t) represents the measured hydraulic gradient on day t, in meters per kilometer.

[0064] The arithmetic mean of the measured hydraulic gradients during all normal water periods is taken as the reference hydraulic gradient i. _ref =(1 / N _total )×Σ t=1 N_total i _actual (t);

[0065] Among them, i _ref For reference hydraulic gradient, N _total Let be the total number of days in the normal water period sample, and Σ represent the summation over all samples. For the confluence of rivers and lakes in the middle and lower reaches of the Yangtze River, i _ref The empirical range for this value is 0.01-0.05 m / km. Reference flow rate Q _ref The characteristic flow method can be used to determine it, and the specific calculation formula is as follows:

[0066] Q _ref =η×Q _warning ;

[0067] Where η is the flow margin coefficient, with a value ranging from 0.8 to 1.2, and Q... _warning This refers to the design discharge capacity at the outlet section corresponding to the warning water level, expressed in cubic meters per second, which can be obtained from the design documents issued by the water resources department. When design data is unavailable, Q... _ref Alternatively, a quantile of the measured flow series over many years at this cross-section can be used, such as the 75th or 90th quantile, i.e.:

[0068] Q _ref =Percentile(Q _historical ,p);

[0069] Among them, Q _historical Here, p is the historical flow sequence, p is the quantile parameter, usually taken as 0.75 or 0.90, and Percentile is the quantile calculation function.

[0070] Optionally, the weighting coefficients α and β are determined using a multiple linear regression calibration method. This method utilizes historical disaster event data and uses regression analysis to find the weight combination that best correlates the coupling index with the actual severity of the disaster. The specific steps are as follows:

[0071] Collect data on K typical flood and drought disaster events in the target area's history. Each event record includes: the hydraulic gradient component G at the time of the event. _norm (k) Flow throughput component Q _norm (k), and the disaster loss index L(k) of the event, which can be the normalized value of economic loss or affected population. Using the disaster loss index L as the dependent variable, and G as the variable... _norm and Q _norm Using the independent variable, establish a linear regression model, that is:

[0072] L(k) = α' × G _norm (k)+β'×Q _norm (k)+ε(k);

[0073] Where L(k) is the disaster loss index of the kth event, α' and β' are the regression coefficients to be estimated, and ε(k) is the residual term.

[0074] Furthermore, the regression coefficients α' and β' are solved using the least squares method, and then normalized to obtain the final weights, i.e.:

[0075] α = |α'| / (|α'| + |β'|);

[0076] β = |β'| / (|α'| + |β'|);

[0077] Where α is the weight of the hydraulic gradient component, β is the weight of the flow throughput component, and |…| represents taking the absolute value. After normalization, α+β=1 is satisfied.

[0078] When historical disaster samples are insufficient, initial values ​​can be set using expert experience, with α=0.6 and β=0.4 recommended. During system operation, the weights are dynamically adjusted based on newly accumulated case data.

[0079] On the other hand, this paper provides a specific implementation process for an iterative temporal clustering algorithm with state transition constraints, which solves the technical problem that traditional static clustering algorithms such as K-means ignore the temporal continuity and dynamic evolution of hydrological processes, easily leading to frequent jumps in state identification results between adjacent time steps. By introducing a state transition penalty term and dynamic time warping (DTW) distance, this scheme can be used to identify continuous hydrodynamic states.

[0080] Step 301: Call the clustering objective function, which includes a clustering compactness error term and a state transition penalty term. The state transition penalty term is used to constrain the state abrupt changes between adjacent time steps, and the clustering compactness error term is used to measure the degree of aggregation of the time series samples of the river-lake hydrodynamic coupling index within each cluster.

[0081] Accordingly, the clustering objective function can be used to measure the quality of state recognition. It not only focuses on the distance between sample points and cluster centers, i.e., compactness, but also incorporates a time smoothing regularization term. The specific mathematical expression is as follows:

[0082] J=Σ t=1 T (Dist(Φ _t μ _S_t )+λ _trans ×Cost(S _t-1 S _t ));

[0083] Where J is the total objective function value; T is the total length of the time series; Φ _t S is the exponential vector of the river-lake hydrodynamic coupling at time t; _t The state label assigned at time t, for example, k=1, 2, ..., K; S _t-1 The corresponding state label assigned at the previous time step; μ _S_t State S _t The corresponding cluster centers; Dist(...) is the distance metric function between the observations and the cluster centers; λ _trans S is the state transition penalty coefficient, used to adjust the strength of timing smoothing; Cost(...) is the state transition cost function, S _t-1 S _t The Σ value is positive (e.g., 1) when the inequalities are not equal, and 0 otherwise. In other words, Σ t=1 T (Dist(Φ _t μ _S_t This corresponds to the cluster compactness error term.

[0084] Furthermore, when calculating the distance Dist(...), this embodiment may optionally use Dynamic Time Warping (DTW) distance, the specific formula of which is as follows:

[0085] Dist _DTW (A, B) = min _π (Σ (i，j)∈π d(A _i B _j ));

[0086] Where A and B are two time series segments to be compared; π is the alignment path; d(A _i B _jLet (i, j) represent the basis distance between pairs of points in the sequence, such as the Euclidean distance. Alternatively, (i, j) can be interpreted as an index pair in π, where i represents the element index of time series A, j represents the element index of time series B, and min represents the minimum value function. Through DTW (Digital Time-Displacement Wrap), the algorithm can identify hydrological processes with similar shapes but misaligned time periods as belonging to the same state, rather than misclassifying them as outliers.

[0087] Furthermore, the state transition penalty coefficient is used to balance the relationship between cluster compactness and temporal smoothness, and its value directly affects the stability of the state recognition results. λ _trans A value that is too small will cause frequent state transitions, λ _trans If the value is too large, the algorithm may tend to classify the entire sequence into a single state. A specific method for determining its value can be a grid search method based on a validation set, as follows:

[0088] The historical hydrological sequence was divided into a training set (70%) and a validation set (30%), where the true state labels of the validation set were pre-labeled by water conservancy experts based on historical records;

[0089] Set λ _trans The candidate value set Λ = {0.1, 0.5, 1.0, 2.0, 5.0, 10.0}.

[0090] For each candidate value λ _trans (i) Run an iterative temporal clustering algorithm on the training set to obtain the clustering model parameters; perform state prediction on the validation set and calculate the consistency index between the predicted labels and the true labels, such as adjusting the RAND index ARI or normalized mutual information NMI, i.e.:

[0091] Score(i) = ARI(S _pred (i), S _true );

[0092] Where Score(i) is the evaluation score corresponding to the i-th candidate value, S _pred (i) is to use candidate value λ _trans (i) The obtained predicted state sequence, S _true Given the true state sequence, ARI is the adjusted Rand index calculation function. The candidate value that maximizes the score is selected as the state transition penalty coefficient, as follows:

[0093] λ _trans_opt =argmax λ_trans(i)∈Λ Score(i);

[0094] Where, λ _trans_opt λ represents the optimal state transition penalty coefficient, and argmax represents the parameter value that maximizes the objective function. Based on practical application experience in the Yangtze River Basin, λ _transThe value range is 1.0-5.0. When the monitoring data is noisy, it is recommended to use a larger value, such as 5.0, to enhance the smoothing effect; when a rapid response to sudden changes in state is required, it is recommended to use a smaller value, such as 1.0.

[0095] Step 302: The clustering objective function is solved by an alternating optimization strategy. When the cluster centers are fixed, the dynamic programming algorithm is used to optimize the state label sequence to minimize the cumulative cost. When the state label sequence is fixed, the cluster centers of each state are updated according to the current classification samples.

[0096] Since the objective function contains discrete state sequence variables, it cannot be solved by direct differentiation. Therefore, an alternating optimization strategy is adopted to decompose the problem into two subproblems.

[0097] When the cluster center μ _k When the state is fixed, find the optimal state sequence S={S _1 S _T This is equivalent to finding the shortest path across the time grid. It can be solved using dynamic programming, with the specific recurrence relation as follows:

[0098] V _t (k)=Dist(Φ _t μ _k )+min _j (V _t-1 (j)+λ _trans ×Cost(j, k));

[0099] Among them, V _t (k) represents the minimum cumulative objective function value at time t when the state is k, where j is the state category index, and V _t-1 (j) represents the minimum cumulative objective function value at time t-1 when the state is j; t-1 is the time preceding t. The forward recursion calculates all V... _t (k) records the path pointer. By backtracking from the endpoint, a globally optimal state label sequence can be obtained.

[0100] When the state sequence is fixed, the problem degenerates into a standard center update problem. For each state class k, all sample points Φ labeled as that state are collected. _t Calculate their arithmetic mean, or the centroid in DTW space, as the new cluster center μ. _k '.

[0101] μ _k '=(1 / N _k )×Σ _t∈{τ|S_τ=k} Φ _t ;

[0102] Where, N _kS is the total number of samples assigned to state k. _τ =k represents the state label at time τ as k.

[0103] Step 303: When the iteration meets the convergence condition, the output state label sequence is taken as the hydrodynamic state of the river and lake. In other words, the convergence condition is that the objective function value decreases by less than a preset threshold in two consecutive iterations, or the state label sequence no longer changes.

[0104] The above sub-steps will be executed alternately in a loop. Each iteration will cause the overall objective function to decrease monotonically. When the decrease in the objective function value between two consecutive iterations is less than a preset threshold, such as 1 × 10⁻⁶, the objective function will be considered complete. -4 The algorithm is considered converged when the state label sequence no longer changes.

[0105] For example, suppose a coupling index sequence Φ=[0.1, 0.12, 0.8, 0.85, 0.1] of length T=5 is monitored.

[0106] Initially, if only K-means is used, the two middle high-value points may be clustered into state B, and the first and last low-value points into state A, with the sequence being AABBA.

[0107] However, when introducing the state transition cost Cost(A, B) = 10 and λ _trans When the value is large, if the duration of the high value in the middle is short, only two time points, the increase in the penalty term resulting from converting it to state B may exceed the reduction in compactness error. The algorithm may tend to smooth the entire sequence into a single state A, or identify it as a special short-duration pulse state C. It should be understood that this mechanism filters out instantaneous fluctuations, making the output hydrodynamic state of rivers and lakes stable and of disaster prevention guidance significance.

[0108] In some alternative implementations, the state transition cost Cost(j,k) can be set as an asymmetric matrix. For example, the cost of transitioning from leakage to backflow can be set lower than that of the reverse transition to reflect the need for early warning of disaster prevention in response to severe operating conditions.

[0109] According to one aspect of this application, the iterative temporal clustering algorithm requires pre-specifying the number of cluster states K, i.e., the total number of categories of river and lake hydrodynamic states. The value of K is determined using an automatic selection method based on silhouette coefficients, with the following steps:

[0110] The candidate range of K can be set to K. _min To K _max K is usually chosen. _min =2,K _max =6;K _min K _max Let it be its minimum and maximum values.

[0111] For each candidate value K(j), run the iterative temporal clustering algorithm until convergence to obtain the state label sequence S(j) and the cluster center set μ(j); where j is the index.

[0112] Accordingly, the average silhouette coefficient of this clustering result is calculated as Silhouette(j) = (1 / T) × Σ t=1 T ((b(t)-a(t)) / max(a(t),b(t)));

[0113] Where Silhouette(j) is the average silhouette coefficient corresponding to the candidate value K(j), T is the total length of the time series, and a(t) is the sample Φ. _t The average distance to other samples within the same cluster, b(t) is the sample Φ. _t The average distance to all samples in the nearest neighbor cluster, where max represents the larger of the two values.

[0114] Furthermore, the candidate value that maximizes the average profile coefficient is selected as the state number, i.e.:

[0115] K _opt =argmax K(j)∈{K_min，...，K_max} Silhouette(j);

[0116] Among them, K _opt This represents the optimal number of cluster states.

[0117] According to hydrological knowledge, the hydrodynamic state of rivers and lakes can generally be divided into 3 to 5 categories, including: strong outflow state, weak outflow / balance state, weak backwater state, strong backwater state, and backflow state. Therefore, in practical applications, it is recommended to set the search range of K to 3 to 5.

[0118] As one possible implementation method, this paper describes the specific implementation process of state identification and risk classification based on Hidden Markov Models. It addresses the technical problems of simple clustering methods struggling to handle nonlinear transition probabilities between states and the inability to directly utilize historical prior knowledge. By employing a probabilistic graphical model, it achieves the coupling of state identification and risk assessment. Specifically, it includes:

[0119] Step 401 calls the state transition matrix describing the probability of transitions between the hydrodynamic states of the river and lake, and the observation probability distribution describing the probability of observing a specific coupling index in each state; where the state transition matrix and the observation probability distribution both belong to the Hidden Markov Model.

[0120] Alternatively, one can invoke a Hidden Markov Model and extract its state transition matrix and observation probability distribution.

[0121] The Hidden Markov Model (HMM) is defined by the triple λ = (A, B, π). The state transition matrix A describes the system's transition from physical state S. _i Evolved to another state S _j The inherent laws governing hydrological processes play a crucial role because certain state transitions are unlikely or rare, such as a direct jump from extremely low water levels to a major flood, while other transitions are highly probable, such as a flood peak following a rising water level. The observation probability distribution B describes the likelihood of observing a specific river-lake hydrodynamic coupling index under a given state, and its specific formula is expressed as follows:

[0122] A _ij =P(S _t =j|S _t-1 =i);

[0123] Among them, A _ij Let be the state transition probability, satisfying Σ _j A _ij =1, P(…|…) corresponds to the conditional probability operator, or in other words, the state label S assigned in the previous time step. _t-1 When i, the state label S assigned at time t _t The probability of being j.

[0124] B _j (Φ _t )=P(Φ _t |S _t =j);

[0125] Among them, B _j (Φ _t Let N(μ) be the observation probability density function, which is usually assumed to follow a Gaussian distribution. _j , Σ _j Above, μ _j Let Σ represent the mean vector in the j-th state. _j Let represent the covariance matrix of the j-th state, and N represent the Gaussian distribution (normal distribution).

[0126] Step 402: Use the Baum-Welch algorithm to estimate the parameters of the state transition matrix and the observation probability distribution based on historical observation data.

[0127] The model parameters are not manually set, but are automatically learned from a large number of historical hydrological sequences through machine learning algorithms. The Baum-Welch algorithm is an iterative algorithm based on the expectation-maximization (EM) principle. In the E-step (expectation step), based on the current model parameters λ, the posterior probability of each state at each time step is calculated, i.e., the forward-backward variables. In the M-step (maximization step), the parameters A and B are updated using the posterior probabilities to maximize the likelihood probability of the historical observation sequence.

[0128] For example, by inputting the coupling index sequence of water level and flow of the Yangtze River and Poyang Lake over the past 20 years, the algorithm can automatically discover statistically significant state patterns, such as the period from April to June being the period of rising water and backwater, and the period from July to September being the period of high water level backflow, and solidify them into a transition matrix A.

[0129] In other words, the above state patterns are not random fluctuations or accidental phenomena, but rather objective laws that repeatedly appear in long-term observation data; solidifying them into a state transition matrix A helps to achieve accurate modeling and reliable prediction of the evolution laws of hydrological systems.

[0130] Step 403: Based on the current river-lake hydrodynamic coupling index sequence, the Viterbi algorithm is used to decode the most likely state sequence as the river-lake hydrodynamic state.

[0131] Once the model parameters λ are determined, the Viterbi algorithm is used to find an optimal path for real-time monitoring data, i.e., a globally optimal sequence of state labels, maximizing the probability of generating an observation sequence along that path. Similar to dynamic programming, but based on probability product rather than distance accumulation.

[0132] Step 404: Invoke the conditional probability risk assessment model based on the hydrodynamic state of rivers and lakes. The model is configured with independent risk scoring functions for different hydrodynamic states of rivers and lakes.

[0133] Furthermore, a state-dependent risk model was established. For each identified state k, such as a backflow state, a dedicated scoring function f was trained. _k (X), the specific conditional probability formula is as follows:

[0134] P(R=r|S=k,X)=exp(w _k ×X+b _k ) / Σ r' exp(w _k_r' ×X+b _k_r' );

[0135] Where P(R=r|S=k,X) is the probability of disaster risk level r given the current state k and feature vector X (including water level, rainfall, etc.); w _k and b _k These are the weight vector and bias term trained for state k, respectively. Alternatively, r' can be described as the summation index (iterating through all possible values ​​of R); w _k_r' b is the weight vector of X in the k-th state when the corresponding result value is r'; _k_r' This is the bias term for the result value r' in the k-th state.

[0136] Step 405: Using the identified current hydrodynamic state of the river and lake as a condition variable, select the corresponding risk scoring function.

[0137] Once the current status is output, the system immediately activates the corresponding scoring function. For example, if the current status is identified as backflow in the main stream, the system will call a model specifically trained for backflow risk, focusing on features such as tributary levee pressure and backwater duration; if it is identified as local rainstorm overflow, it will call a model focusing on rainfall intensity and drainage capacity. This mechanism improves the targeting of risk assessment.

[0138] Step 406: Input the real-time hydrological data into the selected risk scoring function to calculate the conditional probability distribution of the current disaster risk level, and determine the disaster risk level based on this (conditional probability distribution).

[0139] By substituting real-time hydrological data into the selected model, the probability distribution vector of risk levels is calculated, such as [low risk: 0.1, medium risk: 0.2, high risk: 0.7]. The system selects the level with the highest probability as the final output.

[0140] Optionally, a Long Short-Term Memory (LSTM) network can be combined to predict future states. Specifically, the current and historical HMM state sequences are input into the LSTM network to predict the probability distribution of the state at the next moment, and a risk scoring function is called in advance for pre-assessment. The HMM+LSTM combined architecture realizes a complete chain from current state identification to future trend prediction.

[0141] In this scheme, a river-lake hydrodynamic coupling index was constructed, which includes normalized hydraulic gradient and flow throughput components. By performing HMM-based or iterative temporal clustering analysis on this index, different hydrodynamic states such as backflow, backwater, and leakage can be distinguished from a physical mechanism perspective. This solves the problem that the nature of a disaster (such as whether it is internal flooding or external flooding) cannot be determined solely by water level, and provides a physical benchmark for subsequent decision-making.

[0142] As another possible implementation method, an exemplary scheme combining dynamic water level field construction and inundation extraction based on hydrological connectivity is provided. This solves the technical problems that in complex river networks or areas with dike barriers, Euclidean distance-based interpolation methods, such as inverse distance weighted (IDW), are prone to generating erroneous water level fields by traversing terrain obstacles, and that elevation thresholding methods are prone to misclassifying isolated low-lying areas as inundation zones. Accordingly, this embodiment includes:

[0143] Step 501: Define the hydrological connectivity distance between any two points within the monitoring area as the shortest hydraulic path length along the connected waterway between them. When there is no connected waterway path between the two points, the hydrological connectivity distance is infinite.

[0144] In this step, the hydrological connectivity distance d_H The shortest path distance (SPD) is the sole criterion for measuring the degree of correlation between two points in a hydraulic sense. Unlike straight-line distance, it is defined based on the physical fact that water flow can follow paths such as channels, ditches, or floodplains. The calculation process is typically performed on a rasterized or gridded water map using Dijkstra's algorithm (single-source shortest path algorithm) or the Faster Movement Method (FMM). The specific definition is as follows:

[0145] d _H (A, B) = min _path ∫ A B cost(x)ds;

[0146] Where path represents all possible paths connecting points A and B; ds is the infinitesimal length of the path; cost(x) is the cost of traveling to point x on the path. If x is in water or a floodable area, cost(x) = 1; if it is in a water-blocking structure (such as a dike), cost(x) = ∞. If A and B are blocked, the integral result is infinity.

[0147] For example, suppose point A is the center of the river channel and point B is a depression outside the levee. Although the straight-line distance between A and B is only 50 meters, the water flow cannot reach them directly due to the insurmountable levee in between. Therefore, d _H (A, B) = ∞. The above definition ensures that the water level inside the dike will not affect the water level outside the dike through interpolation.

[0148] Step 502: Invoke the spatiotemporal joint mutation function that includes spatial and temporal dimensions, where the independent variable of the spatial dimension is the hydrological connectivity distance.

[0149] Accordingly, spatiotemporal kriging is used for interpolation, and the variogram γ(h) is used. _s h _t The construction of ). Spatial lag h _s It is no longer Euclidean distance, but the aforementioned hydrological connectivity distance d _H A product-sum model is used to describe spatiotemporal correlation:

[0150] γ _st (h _s h _t )=(k _s ×γ _s (h _s )+1)×(k _t ×γ _t (h _t )+1)-1;

[0151] Where, γ _st (h _s h _t ) represents the spatiotemporal variability function, h_s h _t These correspond to spatial lag and time lag, respectively, k _s and k _t Represents the spatiotemporal coupling coefficient, γ _t (h _t ) represents the pure time variogram, γ _s (h _s Since is a purely spatial variogram, a spherical model can be used, i.e.:

[0152] γ _s (h _s )=C _0 +C _s ×(1.5×(h _s / a _s )-0.5×(h _s / a _s ) 3 ), when h _s ≤a _s ;

[0153] Among them, a _s For variable range, representing the maximum connectivity distance of hydrological correlation; C _0 For the nugget effect; C _s The arch height is determined by h. _s Replace with d _H This function can express the physical property that the correlation is zero even if the distance is close but the waterway is blocked.

[0154] Furthermore, the spatiotemporal joint variogram includes spatial variogram parameters, temporal variogram parameters, and the spatiotemporal coupling coefficient k. _s and k _t The spatial variogram parameters and time variogram parameters can include nugget value, dome height, and range. These parameters are determined using an empirical variogram fitting method. The spatial variogram parameters can be determined using the following method:

[0155] The empirical semivariance between all pairs of water level monitoring stations was calculated, as detailed below:

[0156] γ _emp (h)=(1 / (2×N(h)))×Σ d_H(i，j)∈h (H _i -H _j ) 2 ;

[0157] Where, γ _emp (h) represents the empirical semivariance of hydrological connectivity at a distance of h, N(h) represents the number of station pairs whose distances fall within the interval h, and H _i and H _jThe water level observations at stations i and j are respectively, and d _H (i, j) represents the hydrological connectivity distance between stations i and j.

[0158] Furthermore, the empirical semivariance values ​​are grouped by distance and plotted as scatter plots. A spherical model is then fitted using the nonlinear least squares method, i.e.:

[0159] γ _s (h)=C _0 +C _s ×(1.5×(h / a _s )-0.5×(h / a _s ) 3 When h≤a _s ;

[0160] γ _s (h)=C _0 +C _s When h > a _s ;

[0161] Where, γ _s (h) represents the theoretical space variability function value, C _0 To reflect measurement error and microscale variation, C _s To reflect the intensity of spatial structural variation in arch height, a _s The range represents the maximum distance that reflects spatial correlation. The goodness of fit is determined using the coefficient of determination R0. 2 Evaluation, requiring R 2 Not less than 0.8.

[0162] Based on this, the parameters of the time variation function are determined as follows:

[0163] Similarly, the empirical semivariance between different times at the same station is calculated as follows:

[0164] γ _emp (τ)=(1 / (2×N(τ)))×Σ |t_i-t_j|∈τ (H(t _i )-H(t _j )) 2 ;

[0165] Where, γ _emp (τ) represents the empirical semivariance with a time lag of τ, and N(τ) represents the number of sample pairs with a time lag falling within the interval τ. _i t _j For different observation times, H(t) _i H(t) _j (This corresponds to the water level observation values ​​at different observation times.)

[0166] Accordingly, an exponential model can be used to fit the time variability function γ. _t (τ)=C _0t +C _t ×(1-exp(-3×τ / a _t ));

[0167] Where, γ _t (τ) is the theoretical time variation function value, C _0t C represents the nugget value over time. _t For the camber in the time dimension, a _t The time range is a characteristic scale reflecting the time autocorrelation of water level, and the unit is hours.

[0168] Furthermore, the determination of the spatiotemporal coupling coefficient is as follows:

[0169] The spatiotemporal coupling coefficient can be determined using cross-validation. Specifically, 10% of the monitoring station data is retained as the validation set, and data is processed at different k... _s and k _t For spatiotemporal kriging interpolation under a combination of parameters, the parameter combination that minimizes the root mean square error (RMSE) of the validation set is selected. The formula can be described as follows:

[0170] (k _s_opt k _t_opt )=argmin _k_s，k_t RMSE _cv (k _s k _t );

[0171] Where, k _s_opt and k _t_opt These are the optimal spatial and temporal coupling coefficients, RMSE. _cv For cross-validation root mean square error, argmin is the minimum operator, and the subscript... _k_s，k_t Indicates the spatial correlation coefficient k _s Time correlation coefficient k _t To find the best option.

[0172] Based on practical application experience in the middle and lower reaches of the Yangtze River, the parameter value range is set as follows: Spatial range a _s The distance is 20-100 kilometers, and the time range is a. _t The time range is 6-48 hours, and the spatiotemporal coupling coefficient k _s and k _t All range from 0.5 to 2.0.

[0173] Step 503: Calculate the interpolation weights between known water level monitoring points and the location to be estimated using the spatiotemporal joint variogram function, and solve the dynamic water level field of the entire domain using the spatiotemporal kriging interpolation algorithm.

[0174] Furthermore, the point to be estimated, x _0 At time t _0 water level H * (x _0 , t _0 The following is obtained by linearly weighted summation of known monitoring point data:

[0175] H * (x _0 , t _0 )=Σ i=1 N λ _i ×H(x _i , t _i );

[0176] Where i is the sample index, N is the number of samples, and H(x) _i , t _i Let x be the measured water level value of the i-th sample point. _i Let t be the spatial coordinates of the i-th sample point. _i Let λ be the observation time node of the i-th sample point, and λ be the weight. _i The coefficient matrix of the Kriging equations is obtained by solving the system of equations, which is derived from the variogram γ. _st (d _H The interpolation is filled with Δt. This ensures that the interpolation result is statistically unbiased and has the minimum variance.

[0177] Step 504: Extract the dynamic flooding range using the seed point diffusion algorithm. The algorithm uses the main water body as the seed point to perform neighborhood search, and only includes grids with elevations lower than the predicted water level and spatially connected to the flooded area into the flooding range, excluding isolated depressions.

[0178] After obtaining the continuous dynamic water level field H(x, y) across the entire region, the inundation range needs to be extracted using a digital elevation model Z(x, y); where x is the spatial abscissa and y is the spatial ordinate. To avoid misclassifying pseudo-inundation areas with low elevations but no water source (such as dry basins with high surrounding areas and low central areas) as disaster areas, a seed point diffusion process can be performed, with the specific logic as follows:

[0179] Mark all permanent water areas of the main stream, tributaries, and lakes as submerged seed points and add them to queue Q.

[0180] When queue Q is not empty, remove a point P. Check the 8-neighborhood grid N of P. _i .

[0181] For each N _i If the elevation of that point is Z(N) _i ) < Predicted water level H(N) _i And N _i If it has not yet been labeled, then N_i Mark it as submerged and add it to queue Q.

[0182] Until the queue is emptied, all marked grids constitute the dynamic flooding range M(t).

[0183] Step 505: Overlay the dynamic flooding range with the pre-stored population and asset distribution data to calculate the disaster exposure degree; construct a water depth loss function based on the difference between the predicted water level and the ground elevation, and calculate the basic demand for various materials in combination with the disaster exposure degree.

[0184] The demand for materials depends not only on the inundated area but also on the inundation depth and the value of the disaster-bearing structures. Accordingly, the inundation depth h(x, y) = H(x, y) - Z(x, y) is calculated for each grid. Next, a depth loss function D(h) is defined, for example, using a normalized saturated growth function: D(h) = min(1.0, h / h) _max ); where h _max The maximum water depth that would result in loss, such as 3.0 meters.

[0185] Based on this, the disaster exposure level E is calculated. _total And transformed into basic requirement Q _base ,Right now:

[0186] E _total =Σ (x，y)∈M(t) (Pop(x,y)×w) _pop +Asset(x, y)×w _asset )×D(h(x,y));

[0187] Q _base =E _total ×q _unit ;

[0188] Where Pop(x, y) is the grid population; Asset(x, y) is the asset value; q _unit w is the standard for material consumption per unit of exposure. _pop w _asset For the corresponding weights.

[0189] The above demonstrates the mapping from one-dimensional water level monitoring to three-dimensional flood risk, providing a physical basis for subsequent material calculations.

[0190] This scheme introduces hydrological connectivity distance to replace the traditional Euclidean distance, and constructs a spatiotemporal joint variogram based on this distance for Kriging interpolation. Combined with a seed point diffusion inundation extraction algorithm, it solves the problems of erroneous water level penetration in dike-blocked areas and misclassification of isolated low-lying areas as inundation zones, thus enabling the reconstruction of dynamic inundation ranges under complex terrain.

[0191] As another possible implementation method, this paper describes an alternative implementation of physical state-driven material demand forecasting to address the technical problem that traditional material demand forecasting methods rely solely on historical statistical patterns, such as population proportions, and neglect the impact of different hydrodynamic disaster forms, such as slow soaking and rapid impact, on the differentiated types and quantities of materials. Specifically, this solution can be implemented by performing the following steps:

[0192] Step 601: Preset the physical mapping relationship between the hydrodynamic state of rivers and lakes and the types and quantities of materials, introduce a physical state correction coefficient, and dynamically adjust the basic material requirements under different states.

[0193] Because the consumption rate of a specific material varies drastically under different hydrodynamic conditions at the same flood depth, a system maintenance state-material correction table is used to define correction coefficients. The physical state correction coefficient ξ(S, k) is a function of state S and material type k. The corrected demand Q is... _adj The calculation is as follows:

[0194] Q _adj (k, t) = Q _base (k, t) × ξ(S) _t ,k);

[0195] Among them, S _t The hydrodynamic state of the river and lake identified at time t.

[0196] For example, when sandbags are being discharged, the water flow rises gently, primarily for raising the dikes. In this case, ξ(discharge, sandbag) might be set to 1.0 as a baseline value.

[0197] During backflow, the river water carries high energy and impacts the tributaries, making it prone to piping or breaches, requiring urgent sealing. At this time, ξ (backflow, sandbags) may reach as high as 2.5, and the demand surges by 1.5 times.

[0198] For drainage equipment, in the external discharge state, gravity drainage is the main method, ξ(external discharge, drainage) = 0.5. In the backflow state, the water level of the outer river is higher than that of the inner river, so forced drainage is the only option, ξ(backflow, drainage) = 3.0.

[0199] Through this mechanism, the system can sense the physical properties of the water situation and achieve accurate control of demand.

[0200] On the other hand, the demand adjustment coefficient ξ(S, k) represents the ratio of the demand for the k-th type of material to the baseline demand under the hydrodynamic state S of the river and lake. This coefficient is determined using a historical case regression method, with the following steps:

[0201] Collect material consumption data for M historical flood and drought disaster events in the target area. Each event record includes: the dominant river and lake hydrodynamic state S (m) during the event, and the basic demand Q calculated based on the inundated area and population. _base (m, k), and the actual amount of resources consumed, Q. _actual (m, k), where m represents the event index.

[0202] For each state S and each type of material k, select all event samples where the dominant state is S, and calculate the ratio ξ of actual consumption to basic demand. _sample (m, S, k) = Q _actual (m, k) / Q _base (m, k);

[0203] Above, ξ _sample (m, S, k) is the sample correction coefficient for the k-th type of material under state S in the m-th event.

[0204] The median of the correction coefficients for all samples under the same state S is taken as the final correction coefficient ξ(S, k) for that state-material combination. _sample (m, S, k)|S(m=S});

[0205] Median indicates the median calculation, using the median instead of the mean to reduce the impact of outliers.

[0206] When historical samples are insufficient, the following empirical recommendations can be used as a reference:

[0207] {Lake and river hydrodynamic conditions, flood control sandbags ξ, drainage equipment ξ, life-saving equipment ξ, living supplies ξ;

[0208] Strong leakage state, 0.8, 0.5, 0.6, 1.0;

[0209] Weak leakage / balanced state, 1.0, 1.0, 1.0, 1.0;

[0210] Weak support level, 1.5, 1.5, 1.2, 1.0;

[0211] Strong support state, 2.0, 2.5, 1.5, 1.2;

[0212] Backflow state, 2.5, 3.0, 2.0, 1.5}.

[0213] The above recommended values ​​are based on flood control experience in the middle and lower reaches of the Yangtze River. In practical applications, they should be adjusted according to local conditions.

[0214] Step 602: Based on the known disaster monitoring point data, calculate the composite distance between the point to be predicted and the monitoring point. The composite distance is composed of a weighted average of the geographic spatial distance and the disaster characteristic distance.

[0215] Furthermore, to extrapolate detailed disaster conditions in areas without monitoring stations, such as house damage rates and injury rates, spatial interpolation is required. However, the principle of geographical proximity fails here; for example, two villages may be geographically adjacent, but one may be in high ground and the other in low-lying area, resulting in different disaster conditions. Therefore, a composite distance d is introduced. _composite The calculation formula is as follows:

[0216] d _composite (i, j) = sqrt(w _geo ×d _geo (i, j) 2 +w _feat ×d _feat (i, j) 2 );

[0217] Where i is a known point, j is a point to be predicted; d _geo (i, j) represents the geographic Euclidean distance (or road network distance) between the two points; d _feat (i, j) represents the distance to the disaster characteristics; w _geo and w _feat These are the normalized weighting coefficients.

[0218] Or rather, d _composite (i, j) = sqrt(w _H ×d _H (i, j) 2 +w _feat ×d _feat (i, j) 2 );

[0219] Where i is a known point, j is a point to be predicted; d _H (i, j) represents the hydrological connectivity distance between two points; d _feat (i, j) represents the distance to the disaster characteristics; w _H and w _feat These are the normalized weighting coefficients.

[0220] Step 603: The disaster characteristic distance is calculated based on the difference in water level, difference in terrain elevation, and difference in risk level between the two points.

[0221] Where, d _feat This is used to measure the similarity between two points in terms of their disaster environment. If two points have similar environments, even if they are geographically far apart, their resource demand patterns may be highly correlated. The specific formula can be expanded as follows:

[0222] d_feat (i, j) = |H _i -H _j |×δ _h +|Z _i -Z _j |×δ _z +I(R _i ≠R _j )×δ _r ;

[0223] Where H is the water level; Z is the ground elevation; δ _h δ _z Here, I is the dimension conversion coefficient; I(...) is the characteristic function; R is the risk level, and if the risk levels are different, a fixed penalty distance δ is added. _r .

[0224] Introducing d _feat Subsequently, the interpolation algorithm tends to use sample data from distant reference points with the same elevation, water level, and risk to predict the current point's demand, rather than blindly referencing geographically closest but environmentally different points.

[0225] In other words, obtain the hydrological connectivity distance between two points;

[0226] Based on the differences in disaster risk level and inundation depth between two points, the similarity distance of disaster situations is calculated.

[0227] By weighted and fused hydrological connectivity distance and disaster similarity distance, a composite demand distance is obtained to measure the degree of correlation between material demand between two points.

[0228] Step 604: Using the inverse distance weighted algorithm, the reciprocal of the composite distance is used as the weight to interpolate and calculate the material demand of the point to be predicted.

[0229] Accordingly, the final material requirement Q _final (j) The improved inverse distance weighted (IDW) calculation can be expressed as:

[0230] Q _final (j)=(Σ i=1 N (Q _i / d _composite (i, j) p )) / (Σ i=1 N (1 / d _composite (i, j) p ));

[0231] Among them, Q _iLet be the corrected material demand for the known monitoring point i; p is the power exponent, which can be 2. Since the denominator uses composite distance, it applies to areas with large environmental differences, i.e., d... _feat With large sample points, the weight of the sample points will decrease sharply, ensuring that the prediction results conform to the physical laws of flood occurrence.

[0232] In this embodiment, some methods can be further described as follows:

[0233] Call the association mapping table between different river and lake hydrodynamic states and different types of material demand. The association mapping table stores the corresponding demand correction coefficients.

[0234] Based on the identified current hydrodynamic state of the river and lake, the corresponding demand correction coefficient is matched from the association mapping table;

[0235] The demand for basic materials calculated based on disaster risk classification is adjusted using a demand correction coefficient to obtain the corrected spatiotemporal distribution of material demand.

[0236] For example, machine learning models such as Extreme Gradient Boosting (XGBoost) can also be introduced as residual correctors. That is, the physical model is first used to calculate the base values, and then XGBoost is used to predict the residual ΔQ based on historical error data, with the final result being Q. _final +ΔQ. The above, a dual-drive model combining physics and data, preserves physical interpretability while leveraging the fitting capabilities of data mining.

[0237] On the one hand, it provides a specific implementation method for closed-loop scheduling decision based on dynamic flooding and reinforcement learning, which solves the problem that static path planning cannot cope with road interruption caused by disaster evolution, and the technical problem that traditional optimization algorithms are difficult to respond to large-scale sudden demand in milliseconds, thus constructing a closed loop from environmental perception to intelligent decision-making.

[0238] Step 701: In some embodiments, road intersections within the monitoring area are used as nodes and road segments are used as edges to construct a road network topology map; for each edge, its ground elevation is compared with the predicted water level within the dynamic flooding range, and when the predicted water level is higher than the ground elevation and exceeds the safe passage threshold, the edge is determined to be blocked.

[0239] In other embodiments, a transport network accessibility matrix reflecting the real-time traffic capacity of the road network is constructed, specifically including:

[0240] The road network is abstracted into a graph structure consisting of nodes and road segments, and the ground elevation data along each road segment is extracted.

[0241] Map the water level information corresponding to the dynamic flooding range to the road section and compare the ground elevation of the road section with the predicted water level.

[0242] When the ground elevation of a certain road segment is lower than the predicted water level or the difference is less than the safety margin, the road segment is determined to be in a state of flooding and blockage, and the passage cost value of the corresponding road segment in the transport network accessibility matrix is ​​set to infinity.

[0243] Accordingly, a basic map G=(V, E) is constructed based on GIS data, where V corresponds to vertices and E corresponds to edges. For each edge e... _ij ∈E, obtain the elevation Z of the lowest point of its road surface. _min (e _ij Based on the dynamic water level field H(x, y, t) at time t, calculate the water depth h of this road section. _depth The specific blocking determination logic is as follows:

[0244] If H(e) _ij ,t)-Z _min (e _ij )>h _safe Then the state (e) _ij ,t)=BLOCKED;

[0245] Otherwise, State(e) _ij ,t)=AVAILABLE;

[0246] Among them, H(e) _ij (t) represents road segment e at time t. _ij The predicted water level at the location, State(e) _ij ,t) represents the condition of road segment e at time t. _ij The passage status, h _safe This is a preset safe wading depth threshold. For example, for a typical transport truck, h _safe The maximum depth can be set to 0.3 meters; for special wading vehicles, it can be set to 0.8 meters. That is, under the same flooding scenario, the available road network is different for different vehicle types.

[0247] Alternatively, the specific blocking determination logic can be described as follows:

[0248] If H(e) _ij ,t)-Z _min (e _ij )>-h _safe Then the state (e) _ij ,t)=BLOCKED;

[0249] Otherwise, State(e) _ij ,t)=AVAILABLE.

[0250] The above are the safe wading depth thresholds h. _safeThe value used to determine whether a road segment is passable depends on the type and performance of the transport vehicles. Optionally, based on the vehicle configuration of the emergency transport fleet, a vehicle type-wading capability comparison can be established as follows:

[0251] Vehicle type, typical vehicle, safe wading depth (h) _safe ;

[0252] Ordinary transport vehicles, light trucks, and vans, 0.25 meters;

[0253] Medium-sized transport vehicles, medium-sized trucks and buses, 0.35 meters;

[0254] Heavy transport vehicles, heavy trucks and trailers, 0.45 meters;

[0255] Specialized wading vehicles, military transport vehicles, amphibious vehicles, 0.80 meters;

[0256] Inflatable boats / hovercrafts, water transport vehicles, are not restricted.

[0257] When constructing the accessibility matrix of the transportation network, the road segment traffic status is calculated separately for different vehicle types. For vehicle type c, road segment e... _ij The formula for determining the passage status is:

[0258] State(e _ij ,t,c)=BLOCKED, when H(e _ij ,t)-Z _min (e _ij )>h _safe (c);

[0259] State(e _ij ,t,c)=AVAILABLE, when H(e _ij ,t)-Z _min (e _ij )≤h _safe (c);

[0260] Among them, State(e _ij (t, c) represents vehicle type c at time t for road segment e. _ij The passage status, H(e) _ij (t) represents road segment e at time t. _ij The predicted water level at location Z _min (e _ij ) is road segment e _ij The lowest ground elevation, h _safe (c) represents the safe wading depth threshold for vehicle model c.

[0261] In subsequent route planning, the system automatically selects the corresponding accessibility matrix based on the assigned vehicle type, which helps ensure the safety and feasibility of the planned route for that vehicle type. Furthermore, considering the impact of flow velocity on wading safety, when the area is experiencing backflow or strong backwater conditions, the water flow velocity is relatively high, and the safe wading depth threshold should be reduced by 20% to 30% from the above-mentioned threshold.

[0262] h _safe_adj (c, S) = h _safe (c)×(1-δ _v (S));

[0263] Among them, h _safe_adj (c, S) is the adjustment threshold considering the hydrodynamic state, δ _v (S) is the flow velocity reduction factor, which is 0.2-0.3 under strong jacking or backflow conditions, and 0 under other conditions.

[0264] Step 702: Construct the transport network accessibility matrix using the Big M method, setting the passage cost of blocked edges to infinity, and the passage cost of unblocked edges to vehicle travel time or distance.

[0265] Alternatively, construct an adjacency matrix that varies with time t, i.e., the reachability matrix of the transportation network. Matrix element a _ij (t) represents the travel cost from node i to j. If State(e _ij ,t)==BLOCKED,a _ij (t) = M, where M is a large positive number, such as 10. 9 If State(e) _ij ,t)==AVAILABLE,a _ij (t)=Length(e _ij ) / Speed(e _ij ,t);

[0266] Where Length(...) and Speed(...) correspond to the side length and speed calculation functions, respectively, and Speed(e _ij The speed (t) can be reduced based on rainfall intensity, such as halving the vehicle speed during heavy rain. In this way, the originally connected graph structure is dynamically reconstructed mathematically, forcing subsequent pathfinding algorithms to automatically avoid high-risk areas.

[0267] Step 703: Predict future road traffic conditions based on spatiotemporal graph neural networks to capture the spatiotemporal dependence of road network traffic and flooding risk.

[0268] In other embodiments, invoking the transport network reachability matrix reflecting the real-time traffic capacity of the road network also includes future state prediction based on a spatiotemporal graph neural network, specifically including:

[0269] A spatiotemporal graph neural network model containing graph convolutional layers and gated recurrent units is invoked, where the graph convolutional layers are used to capture the spatial topological correlation of the road network, and the gated recurrent units are used to capture the temporal evolution characteristics of the flooding state.

[0270] Input the historical road segment traffic status sequence and the corresponding water level sequence into the spatiotemporal graph neural network model to predict the probability of the traffic status of each road segment in the future period.

[0271] The transport network accessibility matrix for future time periods is updated based on the predicted traffic state probabilities, supporting cross-time period route planning.

[0272] Optionally, a spatiotemporal graph neural network (ST-GNN) model is introduced to predict the road network state at the next T time steps. The input feature matrix X of the model... _t Includes: current road network traffic, water depth, and rainfall. Graph Convolutional (GCN) operation is used to capture spatial dependencies, yielding a spatial feature output matrix H. _spatial =ReLU(D -1 / 2 AD -1 / 2 X _t W); where ReLU is the modified linear activation function, D is the degree matrix, A is the adjacency matrix, and W is the weight matrix. Temporal convolutional operations (TCN, Temporal Convolutional Network) or LSTM units are used to capture temporal evolution.

[0273] Based on this, the system outputs the passage speed and the probability of blockage at future time t+k. This allows the scheduling system to predict that a bridge may be blocked half an hour later and plan detour routes in advance.

[0274] The spatiotemporal graph neural network model is used to predict road traffic conditions in future time periods. Its network structure adopts an encoder-decoder architecture, with the following specific configuration:

[0275] Input feature matrix X∈R N×T _in ×F Where N is the number of road network nodes (i.e., the number of road segments), and T _in The historical time window length can be 12 time steps, each step is 15 minutes, for a total of 3 hours. F is the feature dimension, including: road section water depth, road section flow rate, road section speed, rainfall intensity, and upstream water level, for a total of 5 dimensions.

[0276] Optionally, a two-layer graph convolutional network (GCN) is used to capture spatial correlations, namely:

[0277] H _1 =ReLU(D _norm ×A×X×W _1 );

[0278] H _2 =ReLU(D _norm ×A×H _1 ×W _2 );

[0279] Among them, H _1 and H _2 These are the outputs of the first and second layer graph convolutions, respectively. ReLU is the corrected linear activation function, and D... _norm Let W be the normalized form of the degree matrix, where A is the road network adjacency matrix, and W is the degree matrix. _1 ∈R F×64 W is the first layer weight matrix. _2 ∈R 64 ×32 This is the weight matrix for the second layer.

[0280] Optionally, a two-level gated recurrent unit (GRU) is used to capture the temporal evolution characteristics, which can be specifically expressed as the following formula:

[0281] h _t =GRU(H _2 (t), h _t-1 );

[0282] Among them, h _t Let h be the hidden state vector at time t. _t-1 The hidden state vector is the one from the previous time step. GRU is a gated recurrent unit operation, and the hidden layer dimension is set to 64.

[0283] Optionally, a fully connected layer is used to map the hidden state to the future T. _out The probability of passage status at each time step can be expressed by the formula:

[0284] Y _pred =Sigmoid(FC(h _T ));

[0285] Among them, Y _pred ∈R N×T_out The predicted output represents the probability of blockage for each road segment at each future time step, T. _out To predict the time window length, four time steps, totaling one hour, can be used. The Sigmoid activation function compresses the output to between 0 and 1. The FC (fully connected) layer... _T This represents the GRU hidden state vector at the last time step T.

[0286] Furthermore, the loss function is a binary cross-entropy loss: Loss = -(1 / (N×T)). _out ))×Σ _i，t (Y _true (i, t) × log(Y) _pred(i,t))+(1-Y _true (i,t))×log(1-Y _pred (i, t)));

[0287] Among them, Y _true This is a label indicating the actual passage status; 0 means passage is permitted, and 1 means passage is blocked.

[0288] Based on this, the optimizer uses the Adam algorithm, with an initial learning rate of 0.001, a batch size of 32, and 100 training rounds. An early stopping strategy is used to prevent overfitting, and training stops when the validation set loss does not decrease for 10 consecutive rounds.

[0289] Step 704: Construct a multi-objective material scheduling optimization model with the objective functions of minimizing total delivery time, minimizing material shortage, and maximizing the coverage rate of high-priority areas, and using the transportation network accessibility matrix as a hard constraint.

[0290] In other embodiments, generating a material scheduling and route planning scheme specifically includes:

[0291] A multi-objective optimization model is constructed. The objective functions of the multi-objective optimization model include minimizing the weighted transportation time, minimizing the total transportation cost, and maximizing the coverage of material demand.

[0292] In the multi-objective optimization model, dynamic accessibility constraints are set. The constraints limit the establishment of transportation path associations only when the travel cost between the warehouse and the disaster point in the transportation network accessibility matrix is ​​less than infinity.

[0293] The multi-objective optimization model was solved using a non-dominated sorting genetic algorithm, and the Pareto optimal solution set was obtained as a material scheduling and route planning scheme.

[0294] Furthermore, the optimization model is used to find the optimal set of vehicle paths. The objective function F consists of the following three parts:

[0295] Total travel time objective function F _time =Σ(travel time + node waiting time) → min;

[0296] Objective function for material shortage F _shortage =Σ(Demand - Actual Delivery) → min;

[0297] High-risk area coverage objective function F _coverage =Σ(High-risk area weight × whether delivery was successful) → max;

[0298] Where →min indicates that the objective function takes its minimum value.

[0299] In the constraints, x _ijk ×a_ij (t)<M;→max indicates that the objective function takes the maximum value.

[0300] Among them, x _ijk is a 0-1 variable, indicating whether vehicle k passes through section i-j. This constraint ensures that the vehicle never passes through the blocked section, because once it does, the total cost will become infinite through M and will be excluded by the optimizer.

[0301] Step 705, use a dynamic decision-making algorithm based on reinforcement learning to solve the optimization model. Define the state space to include the current material inventory, the demand at each demand point, and the road network state, and the action space is the material allocation plan and the vehicle path adjustment instruction.

[0302] In some other embodiments, generating a material scheduling and path planning scheme further includes real-time dynamic decision-making based on deep reinforcement learning:

[0303] Model the material storage and transportation process as a Markov decision process. Define the state space to include the current river hydrodynamic state, the disaster risk level, the warehouse inventory, the material demand, and the road network traffic state;

[0304] Define the reward function as the sum of the positive incentive of the material demand coverage rate and the negative penalties of the transportation time, transportation cost, and unmet demand;

[0305] Use the Soft Actor-Critic (SAC) algorithm to train the policy network and output dynamic material scheduling actions or replenishment actions according to the real-time observed state space.

[0306] Given the high dynamics of the disaster environment, static planning may fail. The SAC algorithm can be used for real-time decision-making. Further, the Markov decision process MDP five-tuple is defined as follows:

[0307] State = [vector of inventory at each warehouse, vector of demand at each disaster point, vector of road network blockage status, vehicle location].

[0308] Action = [next node of the vehicle, vehicle loading capacity].

[0309] Reward = -(delivery time) + λ _1 × (quantity of successfully delivered materials) - λ _2 × (penalty for entering the dangerous area);

[0310] Among them, the penalty for entering the dangerous area is set relatively high, and λ _1 and λ _2 correspond to the weight coefficients.

[0311] On this basis, after offline training, the policy network can output optimal actions in milliseconds according to the real-time state.

[0312] On the other hand, the SAC algorithm is used to train the dynamic decision-making strategy network for resource scheduling, and it can be configured with the following hyperparameters:

[0313] Optionally, the policy network employs a three-layer fully connected network with hidden layer dimensions of 256 and 256, respectively, using ReLU as the activation function. The output layer uses a dual-head structure to output the action mean and log-standard deviation, respectively. The value network uses a dual-Q network architecture to alleviate the overestimation problem. Each Q network is a three-layer fully connected network with hidden layer dimensions of 256 and 256, respectively.

[0314] Optionally, the learning rate can be configured as follows:

[0315] lr _actor =3×10 -4 ;lr _critic =3×10 -4 ;lr _α =3×10 -4 ;

[0316] Among them, lr _actor For the policy network learning rate, lr _critic For the value network learning rate, lr _α The learning rate is automatically adjusted for the entropy coefficient.

[0317] Other hyperparameters can be γ=0.99; τ=0.005; where γ is the discount factor, used to balance the weights of immediate and future rewards, τ is the target network soft update coefficient, used to control the speed at which the target network parameters approach the current network parameters; α _init =0.2; where α _init The initial value of the entropy coefficient is used to control the degree of exploration of the strategy. This coefficient is dynamically updated during training through an automatic adjustment mechanism.

[0318] buffer _size =1000000; batch _size =256; warmup _steps =10000;

[0319] Among them, buffer _size For the experience replay buffer capacity, batch _size The batch size for sampling during each update, warmup _steps This is the minimum number of samples that need to be collected before starting a policy update.

[0320] Optionally, the state vector s _t Its dimensions and composition are as follows:

[0321] s _t =[S _hydro R_level I _warehouse D _demand A _road P _vehicle ];

[0322] Among them, S _hydro R is a one-hot encoded vector of the hydrodynamic state of the river and lake, with the dimension being the number of state categories. _level This is a risk level vector for each disaster site, with the dimension being the number of disaster sites, N. _d I _warehouse This represents the current inventory vector for each warehouse, with dimension N. _w ×M (number of warehouses multiplied by number of material types), D _demand This represents the current demand vector for each disaster site, with dimension N. _d ×M, A _road This is a road network traffic state vector, where 0 indicates traffic and 1 indicates blockage. The dimension is the number of road segments, P. _vehicle Encode the current location of each transport vehicle, with the dimension being the number of vehicles N. _v .

[0323] Action vector a _t Hybrid discrete-continuous space design:

[0324] a _t =[a _dispatch a _route ];

[0325] Among them, a _dispatch For resource allocation actions, this represents a matrix of resource quantities distributed from various warehouses to various disaster sites, with dimension N. _w ×N _d ×M, using a continuous action space, a _route The path selection action represents the selection of the next hop node for each vehicle, with dimension N. _v It employs a discrete action space.

[0326] Optionally, the reward function r _t It consists of multiple weighted sub-items:

[0327] r _t =w _1 ×r _coverage +w _2 ×r _time +w _3 ×r _cost +w _4 ×r _safety ;

[0328] Where, r _t Let w be the total reward value at time t. _1 to w _4These are the weighting coefficients for each item. The calculation formulas for each item are as follows:

[0329] r _coverage =Σ i=1 N_d (min(Q _delivered (i), Q _demand (i)) / Q _demand (i));

[0330] r _time =-Σ j=1 N_v T _travel (j) / T _max ;

[0331] r _cost =-C _total / C _budget ;

[0332] r _safety =-λ _danger ×N _danger ;

[0333] Where, r _coverage As a reward for demand coverage, Q _delivered (i) represents the amount of supplies that have been delivered to disaster site i, Q. _demand (i) represents the demand at disaster point i; r _time As a penalty for delivery time, T _travel (j) represents the cumulative travel time of vehicle j, T _max The maximum allowable response time; r _cost C is penalized for transportation costs. _total To accumulate transportation costs, C _budget The upper limit of the budget; r _safety As a safety penalty, N _danger The number of vehicles entering dangerous road sections, i.e., those with water depth exceeding 80% of the threshold, λ _danger This is the penalty coefficient, with a value of 10.

[0334] Alternatively, the weight configuration can be: w _1 =1.0, w _2 =0.3, w _3 =0.2, w _4 =5.0. The safety penalty weight is set relatively high, causing the strategy to prioritize avoiding dangerous road sections.

[0335] Step 706: Iteratively optimize the scheduling strategy based on the feedback from the reward function, which comprehensively considers both the timeliness of material delivery and the success rate of risk avoidance.

[0336] By continuously interacting with the simulated environment (a digital twin built based on historical flood data), the agent learns complex strategies. For example, the strategy of preemptive resource accumulation: when it predicts that a certain area will become isolated, i.e., when ST-GNN predicts a decline in connectivity, the agent will rush to transport and stockpile supplies at any cost before the road is cut off, even if the current demand has not yet surged.

[0337] In this scheme, the dynamic inundation range is mapped to the road network, constructing a real-time updated transportation network accessibility matrix, and introducing a physical state correction coefficient to dynamically adjust material demand. This not only prevents vehicles from entering inundated and blocked road sections, but also calculates the demand for materials such as sandbags based on the differentiated water flow impact intensity, thereby ensuring the physical feasibility and targeted nature of the scheduling scheme.

[0338] On the other hand, an exemplary scheme for model training and parameter acquisition is provided, solving the problems of how to initialize the algorithm and how to inject prior knowledge. Furthermore, this embodiment includes:

[0339] Step 801: Construct a training dataset using historical hydrological monitoring data and historical disaster records, and conduct offline training on the hidden Markov model, risk scoring function, and dynamic decision-making algorithm.

[0340] Alternatively, the dataset can be constructed using the following process:

[0341] Collect time-series data on water level, flow rate, and rainfall for the target watershed over the past 10-20 years, as well as corresponding historical disaster reports, including inundation range maps, disaster statistics tables, and material dispatch logs.

[0342] Missing hydrological data were filled using linear interpolation.

[0343] Water conservancy experts were organized to review typical historical periods and manually label real-state conditions such as leakage, backflow, and backwater as ground truth values ​​for HMM and clustering algorithms, which were then used to verify or supervise fine-tuning.

[0344] Based on historical disaster losses (economic losses, affected population), each event is classified into risk levels I-IV according to national standards.

[0345] Step 802, optionally, for the risk scoring function, a gradient boosting decision tree algorithm is used for training, and optimal hyperparameters are determined through cross-validation.

[0346] For each HMM state k, train an independent gradient boosting decision tree (GBDT, such as XGBoost or LightGBM) model f. _k (X).

[0347] The input features X include the current water level change rate and the cumulative rainfall over the past 24 hours, upstream flow forecast, and soil moisture content. The target variable Y represents the historical actual risk level (0-3). The loss function is multi-class Log-Loss (log loss). Furthermore, a grid search is used to determine parameters such as tree depth and learning rate; 5-fold cross-validation is typically employed to prevent overfitting.

[0348] Step 803: For the dynamic decision-making algorithm, construct a digital twin simulation environment containing historical disaster scenarios, and iterate the strategy through multiple rounds of scenario simulation.

[0349] In this context, reinforcement learning agents cannot be directly tested and failed in the real world; instead, they can be trained in simulators. Accordingly, based on DEM and hydrodynamic models, such as the two-dimensional water area simulation system Mike21 or the river analysis system HEC-RAS, the evolution of historical floods can be reconstructed, generating dynamic transport network accessibility matrices and demand.

[0350] Further, initialize the agent policy network.

[0351] At each simulation step, such as 15 minutes, the agent observes the state of the simulated environment and outputs a scheduling action.

[0352] The simulator uses a physics engine to calculate vehicle arrival times and provides rewards accordingly.

[0353] Update the parameters of the policy network using SAC or the Proximal Policy Optimization (PPO) algorithm.

[0354] Based on this, when the average material shortage rate of the agent in the test set scenario is less than 5% and no transportation vehicle is damaged, training is stopped and the model parameters are saved for loading by the online system.

[0355] On the other hand, the determination of the parameters of the variogram (gold value, arch height, range) can also be achieved by fitting the empirical variogram using the least squares method.

[0356] Specifically, calculate the semivariance γ(i,j) of all monitoring point pairs (i,j) = 0.5 × (H _i -H _j ) 2 It is plotted on the hydrological connectivity distance d _H On a graph where (i, j) is the horizontal axis, the optimal parameter combination for the spherical model is determined using curve fitting tools, quantifying the spatial heterogeneity characteristics unique to this watershed.

[0357] Through the aforementioned offline training process, this system has accumulated historical hydrological wisdom about the basin, enabling it to make judgments that conform to physical laws and historical experience during online operation.

[0358] According to one aspect of this application, prior to executing the method, a model is also pre-built, specifically:

[0359] Collect long-term historical hydrological monitoring data and train hidden Markov models or optimize the parameters of iterative time-series clustering algorithms;

[0360] Collect historical disaster data and corresponding hydrological data, and train the risk scoring function and material demand correction coefficient;

[0361] Collect historical records of road submersion and train a spatiotemporal neural network model.

[0362] Collect historical scheduling and exercise data, and train the strategy network by interacting with the simulation environment.

[0363] On the other hand, assume that the monitoring area is the confluence of the Yangtze River and a certain connected lake, including 1 main stream control station (station A), 1 lake outlet control station (station B), 3 material reserve warehouses (W1, W2, W3), 5 potential disaster points (D1-D5), and a road network consisting of 20 nodes (N1-N20) and 35 road sections.

[0364] At time t _0 The system obtains the following real-time hydrological data: water level H at station A. _main =22.5 meters; B station water level H _sub =22.8 meters; distance between the two stations ΔL=45 kilometers; outlet flow rate Q of station B _out = -800 cubic meters per second (negative value indicates backflow).

[0365] Furthermore, the system pre-stores the following parameters: i _ref =0.025 m / km, Q _ref =12000 cubic meters / second, α=0.6, β=0.4.

[0366] Calculate the hydraulic gradient component G _norm =(22.5-22.8) / (45×0.025)=-0.3 / 1.125=-0.267; where, the negative value indicates that the lake water level is higher than the main stream, and there is a reverse slope.

[0367] Calculate the flow rate and throughput component Q _norm =-800 / 12000=-0.067; where the negative value indicates that the water flow direction is backflow.

[0368] Calculate the river-lake hydrodynamic coupling index Φ(t) _0 ) = 0.6 × (-0.267) + 0.4 × (-0.067) = -0.160 - 0.027 = -0.187; where, Φ(t) _0 ) for t_0 The coupling index at any given time, with a negative value and a large absolute value, indicates a clear backflow situation.

[0369] Accordingly, the system obtains the coupling index sequence Φ from the past 24 hours. _history =[-0.05, -0.08, -0.12, -0.15, -0.18, -0.187, ...], with one sampling point every 2 hours, for a total of 12 points. The values ​​of the first 6 time points are shown here to illustrate the sequence characteristics.

[0370] The sequence is input into a trained Hidden Markov Model (HMM), and the Viterbi algorithm is used to decode it to obtain the state sequence. Based on the decoding result, the current time t... _0 It has been identified as a backflow state, with state number S=5.

[0371] Furthermore, the system calls the risk scoring function f_5(X) trained for the backflow state, and inputs the current hydrological feature vector X=[H _main H _sub Q _out [Cumulative rainfall in the past 6 hours, soil moisture content] = [22.5, 22.8, -800, 45, 0.85];

[0372] The model outputs a risk level probability distribution of [Level I: 0.05, Level II: 0.15, Level III: 0.55, Level IV: 0.25]. The system determines the current risk level to be Level III, indicating a relatively high risk.

[0373] Based on dynamic water level field interpolation and DEM overlay analysis, the system extracts the dynamic inundation range, involving disaster-affected points D2, D3, and D4. The basic infrastructure demand is calculated based on the inundation depth and population distribution. Then, correction coefficients for backflow conditions are applied: flood control sandbags ξ=2.5, and drainage equipment ξ=3.0, i.e.:

[0374] Point D2: Initial requirement: 2000 sandbags, revised to 5000 sandbags; Drainage pumps: 5 units, revised to 15 units.

[0375] Point D3: Initial requirement: 3000 sandbags, revised to 7500 sandbags; Drainage pumps: 8 units, revised to 24 units.

[0376] Point D4: Initial requirement: 1500 sandbags, revised to 3750 sandbags; Drainage pumps: 3 units, revised to 9 units.

[0377] Based on this, the system maps the predicted water level to the road network to determine the traffic status of road segments. The results show that road segment e... _7 The connecting nodes N3 and N8 have a water depth of 0.45 meters, exceeding the threshold of 0.25 meters for ordinary transport vehicles, and are marked as BLOCKED; section e _15Its connecting nodes N10 and N14 have a water depth of 0.55 meters and are marked as BLOCKED.

[0378] The updated reachability matrix shows that the shortest path from warehouse W1 to disaster point D3 is blocked and requires a detour; the path from warehouse W2 to disaster point D4 remains connected.

[0379] The ST-GNN model predicts that within the next hour, road segment e _12 There is a 75% probability that the floodwaters will be blocked. The current water depth is 0.20 meters, and it is predicted to rise to 0.40 meters.

[0380] Optionally, the policy network outputs a scheduling action based on the current state vector, as follows:

[0381] 6,000 sandbags and 20 drainage pumps were transferred from W1 and transported to D2 and D3 by vehicles V1 and V2;

[0382] 4,000 sandbags and 10 drainage pumps were transferred from W2 and transported to D4 by vehicle V3;

[0383] 3,000 sandbags were pre-positioned from W3 to W1 to replenish inventory.

[0384] Optionally, the output path planning results are as follows:

[0385] V1 path, W1→N1→N2→N5→N9→D2, avoids the blocked section e. _7 ;

[0386] V2 path, W1→N1→N4→N6→N11→D3, avoids the blocked section e. _7 and predicting high-risk road sections e _12 ;

[0387] V3 path, W2→N15→N14→D4, direct route.

[0388] Accordingly, the estimated total response time is 45 minutes for V1, 62 minutes for V2, and 28 minutes for V3. The coverage of material needs meets the emergency response requirements.

[0389] Through the above process, after identifying the backflow situation, the system increased the demand for flood control sandbags and drainage equipment, while avoiding already blocked and soon-to-be-blocked road sections, generating a physically feasible and targeted dispatching plan. The calculation process is time-efficient, meeting the timeliness requirements of emergency response.

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

Claims

1. An optimized method for the storage and transportation of flood and drought disaster prevention materials, characterized in that, include: Acquire real-time hydrological data and basic geographic data of the monitoring area, and calculate the river-lake hydrodynamic coupling index reflecting the interaction between the main stream and tributaries based on the real-time hydrological data; Time-series characteristic analysis of the river-lake hydrodynamic coupling index is performed to identify the current river-lake hydrodynamic state in the monitoring area; Disaster risk classification is based on the hydrodynamic state of rivers and lakes, dynamic water level field is constructed by combining hydrological connectivity analysis, and dynamic inundation range and spatiotemporal distribution of material demand in disaster-stricken areas are determined based on the dynamic water level field. Based on the spatial topological relationship between the dynamic inundation range and the road elevation information in the basic geographic data, a transportation network accessibility matrix reflecting the real-time traffic capacity of the road network is constructed. Using the spatiotemporal distribution of material demand as input and the accessibility matrix of the transportation network as constraints, a material scheduling and route planning scheme is generated. The calculation of the river-lake hydrodynamic coupling index, which reflects the interaction between the main stream and tributaries, specifically includes: calculating the normalized hydraulic gradient component based on the water levels of the main stream control station, the water levels of the tributary control station, and the hydraulic distance between the two stations in real-time hydrological data, wherein the normalization process uses a reference hydraulic gradient; calculating the normalized flow throughput component based on the outlet cross-section flow and reference flow in real-time hydrological data; and weighting the hydraulic gradient component and the flow throughput component using weighting coefficients to obtain the river-lake hydrodynamic coupling index. The process of constructing a dynamic water level field by combining hydrological connectivity analysis includes: defining the hydrological connectivity distance between any two points within the monitoring area as the shortest hydraulic path length along the connected waterway between them; when there is no connected waterway path between the two points, the hydrological connectivity distance is infinite; calling a spatiotemporal joint variogram function containing spatial and temporal dimensions, where the independent variable of the spatial dimension is the hydrological connectivity distance; using the spatiotemporal joint variogram function to calculate the interpolation weights between known water level monitoring points and the location to be estimated; and solving the dynamic water level field of the entire area using a spatiotemporal kriging interpolation algorithm. The construction of a transport network accessibility matrix that reflects the real-time traffic capacity of the road network includes: abstracting the road network into a graph structure composed of nodes and road segments, and extracting ground elevation data along each road segment; mapping the water level information corresponding to the dynamic flooding range to the road segment, and comparing the ground elevation of the road segment with the predicted water level; when the ground elevation of a road segment is lower than the predicted water level or the difference is less than the safety margin, the road segment is determined to be in a flooded and blocked state, and the passage cost value of the corresponding road segment in the transport network accessibility matrix is ​​set to infinity.

2. The method according to claim 1, characterized in that, Temporal feature analysis of the river-lake hydrodynamic coupling index is performed to identify the current river-lake hydrodynamic state in the monitoring area. This is specifically achieved through an iterative temporal clustering algorithm with state transition constraints, including: Call the clustering objective function which includes a clustering compactness error term and a state transition penalty term. The state transition penalty term is used to constrain the state abrupt changes between adjacent time steps, and the clustering compactness error term is used to measure the degree of aggregation of the time series samples of the river-lake hydrodynamic coupling index within each cluster. An alternating optimization strategy is adopted to solve the clustering objective function. When the cluster centers are fixed, a dynamic programming algorithm is used to optimize the state label sequence to minimize the cumulative cost. When the state label sequence is fixed, the cluster centers of each state are updated according to the current classification samples. When the iteration meets the convergence condition, the output state label sequence is taken as the hydrodynamic state of the river and lake.

3. The method according to claim 1, characterized in that, A time-series characteristic analysis of the river-lake hydrodynamic coupling index was performed to identify the current river-lake hydrodynamic state in the monitoring area. This was specifically achieved using a hidden Markov model, including: The system calls a state transition matrix that describes the probability of transitions between hydrodynamic states in rivers and lakes, and an observation probability distribution that describes the probability of observing a specific coupling index in each state; both the state transition matrix and the observation probability distribution belong to a hidden Markov model. The parameters of the state transition matrix and the observation probability distribution are estimated based on historical observation data using the Baum-Welch algorithm. Based on the current river-lake hydrodynamic coupling index sequence, the Viterbi algorithm is used to decode the most likely state sequence, which is then used as the river-lake hydrodynamic state.

4. The method according to claim 1, characterized in that, Disaster risk classification based on the hydrodynamic state of rivers and lakes specifically includes: The conditional probability risk assessment model based on the hydrodynamic state of rivers and lakes is invoked. The model is configured with independent risk scoring functions for different hydrodynamic states of rivers and lakes. The current hydrodynamic state of the river and lake is identified as a condition variable, and the corresponding risk scoring function is selected. Input real-time hydrological data into the selected risk scoring function to calculate the conditional probability distribution of the current disaster risk level, and determine the disaster risk level accordingly.

5. The method according to claim 1, characterized in that, Determining the spatiotemporal distribution of material needs in the disaster-stricken area includes introducing a correction mechanism based on the hydrodynamic state of rivers and lakes, specifically: Call the association mapping table between different river and lake hydrodynamic states and different types of material demand. The association mapping table stores the corresponding demand correction coefficients. Based on the identified current hydrodynamic state of the river and lake, the corresponding demand correction coefficient is matched from the association mapping table; The demand for basic materials calculated based on disaster risk classification is adjusted using a demand correction coefficient to obtain the corrected spatiotemporal distribution of material demand.

6. The method according to claim 1, characterized in that, The spatiotemporal distribution of material demand is a continuous field constructed through spatial interpolation. The interpolation process uses composite demand distance, which is calculated as follows: Obtain the hydrological connectivity distance between two points; Based on the differences in disaster risk level and inundation depth between two points, the similarity distance of disaster situations is calculated. By weighted and fused hydrological connectivity distance and disaster similarity distance, a composite demand distance is obtained to measure the degree of correlation of material demand between two points.

7. The method according to claim 1, characterized in that, Constructing a transportation network accessibility matrix that reflects the real-time traffic capacity of the road network also includes future state prediction based on spatiotemporal graph neural networks, specifically including: A spatiotemporal graph neural network model containing graph convolutional layers and gated recurrent units is invoked, where the graph convolutional layers are used to capture the spatial topological correlation of the road network, and the gated recurrent units are used to capture the temporal evolution characteristics of the flooding state. Input the historical road segment traffic status sequence and the corresponding water level sequence into the spatiotemporal graph neural network model to predict the probability of the traffic status of each road segment in the future period. The transport network accessibility matrix for future time periods is updated based on the predicted traffic state probabilities, supporting cross-time period route planning.

Images