Water network dike-flood simulation and disaster loss assessment method

Abstract

The application discloses a kind of water network dike overtopping flood simulation and disaster loss assessment method, belong to the technical field of flood control and disaster mitigation.The method comprises: based on hydrological geographic data, water power analysis is carried out, and potential overtopping area is determined;For potential overtopping area, based on cumulative erosion energy criterion simulation breach evolution, generate breach flow process;According to the breach flow process, simulate flood evolution, generate spatialized disaster intensity field;Combined with disaster intensity field and social and economic data, direct and indirect economic losses are assessed, and comprehensive disaster loss results are generated.The application introduces the physical criterion based on cumulative erosion energy, uses comprehensive intensity index, three-dimensional vulnerability surface and applies graph neural network to simulate cascading failure, solves the problems of missing breach trigger mechanism, single disaster intensity measurement and static disaster loss assessment.

Description

TECHNICAL FIELD

[0001] The application belongs to the technical field of flood control and disaster reduction, and particularly relates to a water network dike overflowing inundation simulation and disaster loss evaluation method. BACKGROUND

[0002] In plain areas with dense river networks and gentle terrain, the safety and stability of dikes, as the first barrier against floods, are of great importance. Once a dike is overtopped by a flood exceeding the standard, it is prone to breach in a short time, leading to rapid spread of the flood in the river network system and large-scale inundation disaster. Therefore, a method capable of accurately simulating the whole process of water network dike overflowing inundation and scientifically evaluating direct and indirect disaster losses is urgently needed to improve the scientificity, foresight and effectiveness of flood control and disaster reduction decision-making, protect the sustainable development of the society and economy, and provide technical support for planning and design of flood control projects, development and optimization of emergency plans, and fine management of disaster risks.

[0003] At present, for the simulation and evaluation of flood inundation, the existing technical solutions include, in terms of hydrodynamic simulation, river network models based on one-dimensional Saint-Venant equation set (such as HEC-RAS) and flood routing models based on two-dimensional shallow water equation (such as TELEMAC and FLO-2D); both of which can effectively numerically solve the water level evolution of the river under specified boundary conditions and the overland flow process of the flood in two dimensions. In the field of disaster loss evaluation, the mainstream method is to establish a vulnerability curve between the disaster-causing factor (usually the inundation water depth) and the loss rate of the disaster-bearing body; for example, by statistical analysis of historical disaster data, a loss rate function of different types of buildings under different inundation water depths is constructed, which is superimposed with asset spatial distribution data to calculate direct economic losses. For indirect losses, some studies have begun to try to use input-output models or computable general equilibrium models. In addition, to deal with the uncertainty in the model, methods such as Monte Carlo simulation are also introduced in the risk assessment process.

[0004] Under the above working conditions, the existing solutions still have problems such as lack of physical criteria, incomplete measurement of disaster-causing intensity, and insufficient consideration of dynamic social factors in disaster loss evaluation. This is because the existing models have limited ability to finely depict the physical process when simulating the complete chain from overtopping to disaster loss, and have failed to achieve a true reflection of multi-system coupling. Therefore, further research is needed to solve the above problems existing in the prior art. SUMMARY

[0005] The application aims to provide a water network dike overflowing inundation simulation and disaster loss evaluation method to solve the above problems existing in the prior art.

[0006] Technical solution: According to one aspect of the application, the water network dike overflowing inundation simulation and disaster loss evaluation method comprises:

[0007] Hydrodynamic analysis based on hydrogeographic data is used to identify potential overtopping areas;

[0008] For potential overtopping regions, the breach evolution is simulated based on the cumulative scour energy criterion to generate the breach flow process;

[0009] Based on the breach flow process, flood evolution is simulated to generate a spatialized disaster intensity field;

[0010] By combining disaster intensity field and socioeconomic data, direct and indirect economic losses are assessed to generate a comprehensive disaster loss result.

[0011] According to another aspect of this application, the method for simulating flooding and assessing damage caused by water network overtopping can also be:

[0012] Hydrodynamic analysis based on hydrogeographic data is used to identify potential overtopping areas;

[0013] For potential overtopping regions, the overtopping overflow process is simulated and the breach evolution is simulated based on the cumulative scour energy criterion to generate the breach flow process;

[0014] Based on the breach flow process, flood evolution is simulated to generate a spatialized disaster intensity field;

[0015] By combining disaster intensity field and socioeconomic data, direct and indirect economic losses are assessed to generate a comprehensive disaster loss result.

[0016] According to another aspect of this application, simulating breach evolution based on the cumulative scour energy criterion includes:

[0017] For each potential overtopping area, determine the instantaneous scouring power of the overtopping flow;

[0018] The instantaneous scouring power is integrated over time to obtain the cumulative scouring energy.

[0019] When the accumulated scour energy exceeds the scour resistance energy threshold set based on the physical properties of the dike material, breach evolution is triggered.

[0020] The time-varying process of geometric morphology after the breach is simulated and coupled with the storage state variables in the river network channel. Based on this (the time-varying process coupled with the storage state variables in the river network channel), the breach flow process is calculated and generated.

[0021] According to another aspect of this application, determining the instantaneous scouring power of the overhead water flow includes:

[0022] Analyze the hydraulic conditions of the potential overtopping area to obtain the unit width flow rate of the overtopping flow;

[0023] Extracting local topographic slopes of potential overtopping areas from hydrogeographic data;

[0024] The Froude number is calculated based on the unit width of the overpass flow, and the kinetic energy correction factor is generated accordingly.

[0025] The instantaneous scouring power is determined by taking into account the unit width flow rate of the overpass water, the local topographic slope, and the kinetic energy correction factor.

[0026] According to another aspect of this application, setting an impact energy threshold includes:

[0027] The physical properties of embankment materials in potential overtopping areas were extracted from embankment monitoring data and embankment deformation radar remote sensing results based on INSAR. These physical properties include soil cohesion and internal friction angle.

[0028] Calculate the critical shear stress of the embankment material based on the soil cohesion and internal friction angle;

[0029] The critical shear stress is converted into an impact energy threshold.

[0030] According to another aspect of this application, after triggering the breach evolution, it further includes:

[0031] The amount by which the accumulated scouring energy exceeds the scouring energy threshold is calculated to obtain the excess scouring energy.

[0032] Based on the super-threshold scour energy, the initial breach geometric parameters are determined;

[0033] Based on the initial breach geometry parameters, the time-varying process of the geometric morphology after the breach is triggered is simulated.

[0034] According to another aspect of this application, a spatialized disaster intensity field is generated, including:

[0035] Based on the breach flow process and hydrogeographic data, the spatiotemporal distribution of inundation depth and velocity field is calculated;

[0036] Based on the spatiotemporal distribution of inundation depth and flow velocity field, a comprehensive intensity index and impact duration are determined for any calculation unit within the inundation area;

[0037] The comprehensive intensity index and impact duration of each calculation unit constitute a spatialized comprehensive intensity index field and an impact duration field, which together form a spatialized disaster intensity field.

[0038] The flooded area is the spatial range determined based on the flooded water depth being greater than a set water depth threshold.

[0039] According to another aspect of this application, determining a comprehensive intensity index for any calculation unit within a flooded area includes:

[0040] Based on the spatiotemporal distribution of submerged water depth and flow velocity field, the water depth time series and flow velocity time series of this computing unit are analyzed.

[0041] Based on hydrogeographic data, the local topographic slope of this calculation unit is extracted;

[0042] For each flooding moment of the computing unit, its water depth, flow velocity, and local topographic slope are coupled to obtain the instantaneous impact power density;

[0043] The comprehensive strength index of the computational unit is obtained by integrating the instantaneous impact power density along the time dimension.

[0044] According to another aspect of this application, before coupling its water depth, flow velocity, and local topographic slope to determine the instantaneous impact power density, the method further includes:

[0045] Identify the direction of water flow evolution based on flow velocity time sequence, and identify the slope direction based on local topographic slope;

[0046] Based on the vector relationship between the direction of water flow evolution and the direction of slope, the directional coupling coefficient is determined;

[0047] The local terrain slope is modulated using the directional coupling coefficient to generate an effective slope.

[0048] The coupling operation specifically involves coupling water depth, flow velocity, and effective slope.

[0049] According to another aspect of this application, the assessment of direct economic losses includes:

[0050] For any calculation unit within the inundation area, extract its comprehensive intensity index and impact duration from the spatialized disaster intensity field;

[0051] Obtain early warning timeframes from socioeconomic data;

[0052] By applying a pre-constructed three-dimensional vulnerability surface, the loss rate of the calculation unit is queried and determined based on the comprehensive strength index, impact duration, and early warning lead time.

[0053] By combining the spatial distribution data of assets and the loss rate contained in the socio-economic data, the direct economic loss of this calculation unit is calculated.

[0054] According to another aspect of this application, after determining the loss rate of the computing unit, the method further includes:

[0055] Assess the prediction uncertainty of the three-dimensional vulnerability surface at the query point and generate a confidence interval for the loss rate;

[0056] Specifically, assessing the uncertainty of prediction includes: calculating the Jacobian matrix of the three-dimensional vulnerability surface at the query point with respect to its input dimension; combining the Jacobian matrix with the uncertainty covariance of the input data to calculate the predicted variance of the loss rate, and generating a confidence interval based on this (prediction variance).

[0057] Beneficial effects: This invention fills the gap in physical criteria for breaching by using the cumulative scour energy criterion, improves the measurement of disaster intensity by using a spatialized disaster intensity field, and enriches disaster damage assessment by combining socio-economic data, thereby achieving accurate simulation of the entire chain from flooding to disaster damage, and providing reliable support for disaster response decision-making in water network areas. Attached Figure Description

[0058] Figure 1 A flowchart of the water network flooding simulation and disaster damage assessment method provided in the embodiments of this application.

[0059] Figure 2 A flowchart illustrating the simulation of breach evolution based on the cumulative scour energy criterion provided in this application embodiment.

[0060] Figure 3 A flowchart for determining the instantaneous scouring power of the overhead water flow, provided for an embodiment of this application.

[0061] Figure 4 This is a flowchart illustrating the setting of the impact energy threshold provided in an embodiment of this application.

[0062] Figure 5 This is a flowchart for determining the comprehensive intensity index of any calculation unit within the flooded area, provided as an embodiment of this application. Detailed Implementation

[0063] 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. Obviously, the described embodiments are merely 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.

[0064] 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 the embodiments of the invention described herein can be implemented in orders other than those illustrated or described herein.

[0065] It should be noted that, for the purpose of clearly demonstrating the steps of this application, each step has been numbered in the specification. These numbers are for ease of explanation only and do not limit the execution order of the steps. In actual operation, depending on the technical requirements of the specific implementation scenario, the steps may be executed in a different order than that shown in the specification, and in some cases, parallel processing between steps can also be achieved.

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

[0067] To determine and simulate the phase transition process from overtopping weir flow to dynamic breach evolution, existing models mostly rely on manually setting the time and location of the breach or using simplified hydraulic thresholds (such as fixed flow velocity), resulting in significant deviations from the actual physical process. Weir failure is not an instantaneous event, but rather a process of continuous work done by water flow on the soil, energy accumulation, and ultimately, material fatigue failure. Current technologies lack physical criteria to characterize this cumulative damage process, failing to clearly define the time and intensity thresholds for breach triggering under overtopping conditions. This simplification of the crucial flood-generating process directly impacts the accuracy of subsequent flood simulations.

[0068] Furthermore, existing disaster loss assessments rely too heavily on a single metric for measuring disaster intensity, failing to comprehensively reflect the actual destructive power of floods. The widely used inundation depth-loss rate curve neglects the powerful impact and scouring capabilities of high-velocity water flows, as well as the cumulative destructive effects over time. For example, the destructive power of still water at a depth of 0.5 meters differs drastically from that of a flash flood at a depth of 0.5 meters but with a flow velocity of 3 m / s. Current technologies lack a comprehensive intensity index that integrates multiple physical dimensions such as water depth (potential energy), flow velocity (kinetic energy), duration of action (cumulative effect), and terrain slope (exacerbating scouring), which can easily lead to distorted disaster loss assessment results in complex terrain and high-velocity scenarios.

[0069] Based on this, existing technologies, when correlating physical damage with loss assessment, do not fully consider the moderating effect of dynamic factors such as human emergency response on losses. For example, for the same physical damage, if there is sufficient warning time and personnel and materials can be evacuated in advance, the resulting loss will be less than in a sudden incident. Most existing models are static and fail to intrinsically couple the warning lead time into the loss assessment model. This makes it impossible to assess the actual loss reduction effect of the warning system and weakens the practical guiding value of the assessment results.

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

[0071] Example 1: A method for simulating flooding and assessing damage caused by water network overtopping is provided. This method can be applied to urban or regional flood warning and disaster assessment systems. Its operating environment can be deployed on a server cluster consisting of multiple high-performance computers, integrating geographic information system software, hydrodynamic simulation software, and database management system.

[0072] In the context of this invention, to maintain consistency in description, some technical terms and data items of the embodiments are defined as follows: Z _levee h is the elevation of the top of the embankment._levee Here, h is the dam height; h0 is the initial water level; v0 is the initial flow velocity; q and q(t) are the unit width flow rates; and ρ is the density of water, usually taken as 1000 kg / m³. 3 q' represents the unit width flow rate of the overpass, i.e., the volumetric flow rate of water passing through a unit width; ▽z represents the terrain slope; P(τ) is the instantaneous scour power, which is the instantaneous destructive capacity of the water flow acting on the surface of the levee at a certain moment τ; τ is the integral variable (representing any moment, used to traverse the entire process from the initial moment 0 to the current moment t); B is the breach width; B(t) is the time-varying breach width; B _max τ is the maximum width. _b The bottom shear stress is calculated based on the Manning formula; τ _c τ is the critical shear stress. _b_turb The corrected bottom shear stress; h refers to the breach depth, h(t) is the time-varying breach depth; H _up The upstream water level; h _b It is the scour increment relative to the initial breach depth; h _b0 B is the initial breach depth. _b0 The initial breach width; m _b0 Z represents the initial breach slope ratio (dimensionless); m represents the breach slope ratio; m(t) is the time-varying breach slope ratio; _b0 Z represents the initial breach bottom elevation. _b The bottom elevation of the breach; Z _b (t) represents the time-varying breach bottom elevation; Q _b (t) represents the time-varying breach flow rate; V _total (t) represents the time-varying cumulative flow volume; V is the water volume; IM _i T is the quantized value of the adjusted intensity field of the i-th sample in the historical disaster damage case database; _i W represents the impact duration of the i-th sample in the historical disaster damage case database. _i Let S be the advance warning amount for the i-th sample in the historical disaster damage case database; in the graph network, i is the master node number, j is the neighbor node number of i, and S is the early warning amount for the i-th sample. _i S represents the current state of node i; _i S' represents the current state of node i; _i S represents the updated state of node i; _j This represents the current state of node j; t _predict Indicates the predicted flood peak time; t _now Represents the current time; the computational grid (i, j) indicates that i is the row index of the flooded region grid and j is the column index, representing a specific spatial unit; the Logistic function P=1 / (1+exp(-α(XX)) 50In the logistic curve, the slope parameter (sensitivity coefficient) of the α function is used to characterize the sensitivity of the input factor X to the probability P, i.e., the steepness of the logistic curve; X is a quantitative indicator directly related to the target event (such as facility damage, regional disaster, disease occurrence, etc.), and its value is usually a continuous variable (such as depth, time, intensity, etc.), and it shows a positive correlation with P (i.e., the larger X is, the closer P is to 1; the smaller X is, the closer P is to 0); X 50 For the critical value parameter of the function; {IM _k T _k W _k L _k In the dataset, k is the index number of the sample in the training dataset.

[0073] In this embodiment, one possible implementation is as follows:

[0074] Step S1.1: Conduct hydrodynamic analysis based on hydrogeographic data to determine potential overtopping areas.

[0075] Hydrological and geographic data may include, but is not limited to: high-resolution digital elevation model (DEM) data, such as data with a spatial resolution of 5 to 30 meters; vector data of river networks, including river centerlines and cross-sectional dimensions; levee engineering data, including the spatial location of levees along rivers and a database of levee crest elevations; and underlying surface characteristic data, such as land use type maps and soil type distribution maps. Socioeconomic data may include: spatial distribution data of assets, such as grid-divided value density maps of housing, farmland, and industrial and commercial assets; and critical infrastructure data, including the type, location, capacity, and interdependencies of substations, water plants, hospitals, transportation hubs, etc.

[0076] Optionally, this step can be achieved by constructing a hydrodynamic network model that couples one-dimensional river channels with zero-dimensional storage units. Specifically, based on river network vector data and cross-sectional measurement data, a river node-segment adjacency matrix is ​​constructed using a topology cleaning algorithm. Spatial proximity analysis is then used to determine the connection relationships between zero-dimensional storage units such as lakes and polder areas and one-dimensional river channels, forming a complete river network topology map.

[0077] Furthermore, based on this topology, the Saint-Venant equations are spatially discretized using the Preissmann four-point implicit scheme. Real-time monitored hydrological sequences (such as upstream water level and flow rate) are configured as time-varying boundary conditions for the model using spatial interpolation methods such as Kriging interpolation. Solving this discretized nonlinear equation set yields the time-series distribution of water level and flow rate across all computational nodes in the network. Optionally, the Newton-Raphson iterative method can be used for solving the equations, supplemented by Aitken acceleration techniques to improve convergence speed.

[0078] Based on this, potential overtopping areas are identified, including: calculating the superelevation margin between the dike elevation and the water level; monitoring the rate of change of the superelevation margin; setting tiered early warning systems based on the superelevation margin values; and predicting the future overtopping time based on the superelevation margin and its rate of change. Specifically, this involves real-time comparison of the calculated water level H(i,t) at each node with the dike crest elevation Z in the dike elevation database. _levee (i) Perform comparison and calculate the super-high margin ΔH=Z _levee -H. When the super-high margin ΔH is less than the preset danger threshold (e.g., 0.3 meters) or less than or equal to 0, the cross section is identified as a potential overtopping area or a dangerous overtopping cross section. Through this identification process, the overtopping risk is transformed from a static defense standard to a dynamic process early warning, providing a target area for subsequent breach simulation.

[0079] Alternatively, another implementation scheme is provided below: Based on river network topology data and real-time hydrological monitoring sequences, a one-dimensional river channel-zero-dimensional storage unit coupled hydrodynamic network model is constructed, and the temporal distribution of water levels at all network nodes is obtained by solving the Saint-Venant equations; the calculated node water levels are compared in real time with the levee elevation database to screen out a set of potential overtopping sections with water levels close to or exceeding the levee crest elevation, and the initial overtopping unit width discharge of each dangerous section is initially estimated using the broad-crested weir formula.

[0080] Alternatively, another implementation scheme is provided below:

[0081] The original river centerline vector data and river cross-section measurement data are read, and suspended nodes and duplicate river segments are removed using a topology cleaning algorithm. A river node-segment adjacency matrix is ​​constructed to generate a standardized river network topology map G(V, E). The boundary polygons of storage units such as lakes and polder areas are converted into zero-dimensional storage nodes. Spatial proximity analysis is used to determine their lateral connectivity with the one-dimensional river, forming an extended topology matrix G', which contains complete connectivity information between river channels and between river channels and storage units.

[0082] Based on the extended topological matrix G', the Saint-Venant equations are spatially discretized using the Preissmann four-point implicit scheme, and the real-time water level monitoring data Z is then processed. _obs (t) and flow monitoring data Q _obs (t) is mapped to the model boundary nodes through Kriging interpolation to generate the time-varying boundary condition set BC(t). For internal nodes, the initial water level field is set according to historical water level statistics, and the initial velocity field is estimated by Manning's formula and the river roughness parameter library to complete the construction of the initial state vector X0=[h0, v0].

[0083] An implicit finite difference solver was run, and the hydrodynamic simulation was performed in 15-minute time steps to obtain the time series water level H(i,t) and cross-sectional flow Q(i,t) for all network nodes. The calculated water level H(i,t) was compared in real time with the levee elevation database Z. _levee (i) Perform comparison and calculate the super-high margin ΔH=Z _levee -H, when ΔH < 0.3m, the section is marked as a warning section, and when ΔH ≤ 0, it is marked as a dangerous section with roof fall. _risk For the overlying section, the broad-crested weir formula q=Cd*sqrt(2g)H is adopted. 1.5 Calculate the initial overflow unit width flow q0(i), where the overflow coefficient Cd is dynamically adjusted according to the roughness of the embankment crest.

[0084] Step S1.2: For the potential overtopping region, simulate the breach evolution based on the cumulative scour energy criterion to generate the breach flow process. In some embodiments, for the potential overtopping region, the overtopping overflow process can also be simulated and the breach evolution can be simulated based on the cumulative scour energy criterion to generate the breach flow process.

[0085] In some embodiments, after identifying the potential overtopping region, the system will activate a physics-driven overtopping-breach phase transition simulation module for that region. Since breaching is not an instantaneous event, but rather the result of continuous work done by the water flow on the embankment, accumulating energy to the material's scour resistance limit, the module calculates the instantaneous scour power of the overtopping flow in real time and integrates it over time to obtain the accumulated scour work (i.e., accumulated scour energy). When this accumulated value exceeds the critical scour resistance threshold set based on the physical properties of the embankment soil material (such as cohesion and internal friction angle), the model will trigger a phase transition, switching from a stable weir flow mode to a dynamic breach evolution mode. In the breach evolution mode, the model will simulate the expansion process of the breach geometry (such as width and depth) over time, and based on this time-varying geometry, calculate the flow rate through the breach in real time, generating a complete breach flow process line Q. _breach (t). This process line accurately describes the total amount and peak value of floodwater released during the breach and is a key input for subsequent flood simulations.

[0086] In other embodiments, for the identified potential overtopping cross-section and its initial overtopping unit width flow, the unit width power P=ρgq'|▽z| of the overtopping flow is calculated. When the cumulative power ∫Pdt exceeds the scour resistance threshold of the embankment material, a phase transition from stable weir flow to time-varying breach is triggered. The breach propagation rate is dynamically updated based on the nonlinear relationship between local shear power and soil erosion resistance, generating a time-varying breach geometric parameter sequence (bottom width, depth, slope). The breach flow process curve is calculated in real time using the improved Breach formula, outputting the spatially distributed cumulative scour power field for subsequent damage assessment.

[0087] In some other embodiments, for the dangerous overtopping section D _risk The initial overpass unit width flow rate q0, combined with the local terrain slope ▽z, is extracted from the DEM. The instantaneous unit width scour power P(t) is calculated as: P(t) = ρg * q0(t) * |▽z| * (1 + 0.5Fr) 2 The Froude number Fr reflects the amplification effect of water kinetic energy. A power accumulator integrator is established to calculate the cumulative scouring power W in real time. _cum (t)=∫0 t P(τ)dτ, and the critical scour resistance threshold W in the embankment material database. _cr A comparison was performed, among which W _cr =f(soil cohesion c, internal friction angle φ, compaction degree γ).

[0088] When the cumulative flushing power W _cum Exceeding the critical threshold W _cr When the phase transition flag Flag_breach is set to 1, the system automatically switches from weir flow mode to breach evolution mode. At the moment of transition, the excess power ΔW = W... _cum -W _cr Calculate the initial breach depth h _b0 =k1*(ΔW / W _cr ) 0.6 Initial bottom width B _b0 =k2*h _b0 *(1+tanβ), where k1 and k2 are material correlation coefficients, and β is the angle of repose. These initial breach geometry parameters serve as the starting point for time-varying evolution and are input into the next calculation module.

[0089] Establish the system of differential equations for breach propagation: Base width evolution dB / dt = α*τ _b 1.5 *(1-B / B _max ), deep evolution dh / dt=β*(τ) _b -τ _c )*exp(-γ*t), where τ _b The bottom shear stress is determined by the instantaneous breach flow rate Q. _b Obtained through inverse calculation. The adaptive Runge-Kutta method is used to solve the equation system, generating a time-varying breach geometry sequence [B(t), h(t), Z]. _b (t)]. Based on the updated breach geometry, the improved weir-orifice flow combination formula Q is used. _b =f(H _up B, h, Z _b Real-time calculation of the breach flow rate process line Q _breach (t) is fed back as a lateral source term to the river network hydrodynamic model to ensure mass conservation.

[0090] Monitor the breach propagation rates dB / dt and dh / dt, and when both are simultaneously less than the threshold ε (10 -4 When the speed is (m / s) and the duration exceeds 30 minutes, the breach is considered to have reached a quasi-steady state; the final breach size [B] is recorded. _final h _final and peak flow Q _peak Output the cumulative flow volume V throughout the entire evolution process. _total This is used for subsequent calculations of the flooding range.

[0091] Step 1.3: Simulate flood evolution based on breach flow process to generate a spatialized disaster intensity field.

[0092] In this step, the breach flow process line Q generated in the previous step is... _breach (t) serves as a point source or line source boundary condition, driving the two-dimensional shallow water equation solver. This solver simulates the evolution of a flood on a two-dimensional surface after leaving a breach, using high-resolution DEM topographic data, and calculates the spatiotemporal distribution of water depth h(x, y, t) and flow velocity v(x, y, t) for each computational grid (e.g., 30 m × 30 m) within the inundation area.

[0093] Optionally, the solver can employ an HLLC approximate Riemann solver based on the finite volume method to handle shock waves and discontinuities. Further, this invention uses the comprehensive intensity index IM as the disaster intensity index; it represents the total impact energy borne per unit area and can be obtained by integrating the instantaneous impact power density (coupled with water depth, flow velocity, and local topographic slope) of each grid over the entire inundation duration. Based on this, the IM values ​​of all computational grids collectively constitute a spatialized comprehensive intensity index field IM(x, y); the effective impact duration T(x, y) of each grid can also be statistically analyzed. Both together form a spatialized disaster intensity field, providing a basis for measuring hazard-causing factors in subsequent disaster damage assessment.

[0094] In another possible implementation, the breach flow process line is used as the source term to drive the two-dimensional shallow water equations, simulating flood evolution on a high-resolution DEM topography to obtain the spatiotemporal distributions of the inundation depth field h(x, y, t) and the velocity field v(x, y, t). Further integrating water depth, velocity, and topographic slope information, a comprehensive intensity index IM=∫ρgh*v*|▽z|dt is constructed, which physically characterizes the cumulative impact work per unit area. Combining this with a land use type layer, exposure adjustment coefficients for different land types are calculated differentially, generating a spatialized intensity index field IM(x, y) and the corresponding impact duration distribution T(x, y) considering land cover characteristics.

[0095] In another possible implementation, the breach flow process line Q will be... _breach(t) is converted to point source boundary conditions for a two-dimensional model, and a two-dimensional shallow water equation solver is configured on a 30m resolution DEM grid. Manning roughness is set differently based on land use type maps: 0.025 for water bodies, 0.035 for farmland, and 0.05-0.08 for built-up areas. An adaptive time step Δt = CFL * min(Δx / sqrt(gh)) is used to ensure numerical stability, where the CFL number is set to 0.7.

[0096] A two-dimensional solver based on the finite volume method is used, and an HLLC approximate Riemann solver is employed to handle discontinuities, obtaining the transient water depth field h(x, y, t) and the velocity vector field v(x, y, t). For the wet-dry boundary, a thin-layer water depth technique (h... _min =0.001m) to avoid numerical oscillations. The water depth and velocity output at each calculation time step are stored as a spatiotemporal data cube DC[h, v, t], with a time resolution of 15 minutes and a spatial resolution of 30m.

[0097] The water depth time series h(t) and flow velocity time series v(t) of each grid are extracted from the spatiotemporal data cube DC. Combined with the slope field |▽z| of the DEM, the instantaneous impact power density p(x, y, t) = ρg*h*v*|▽z|*(1+v) is calculated. 2 / 2g). Integrating over the time dimension yields the cumulative intensity index IM(x,y)=∫p(t)dt, which represents the total impact energy per unit area. The effective impact duration T for each grid is then calculated. _iMpact (x, y) = Σ(Δt|h>0.1m), generation intensity-duration binary [IM, T] _iMpact ].

[0098] The original intensity index (IM) is overlaid with a high-resolution land use layer and adjusted according to the erosion resistance characteristics of different land features: IM for farmland areas. _adj =1.2*IM (crop vulnerability), building area IM _adj =0.8*IM (structural impact resistance), road IM _adj =0.9*IM. Generate the adjusted intensity field IM. _adjusted (x, y). Based on statistical quantiles, the intensity is divided into 5 levels: very low [0, P20), low [P20, P40), medium [P40, P60), high [P60, P80), and very high [P80, P100]. Output the intensity level distribution map.

[0099] Step 1.4: Combine disaster intensity field and socioeconomic data to assess direct and indirect economic losses and generate comprehensive disaster loss results.

[0100] For example, after obtaining the spatialized disaster intensity field, it is converted into economic loss through a vulnerability mapping module. For direct economic loss, this invention introduces a three-dimensional vulnerability surface V(IM, T, W). This surface takes the comprehensive intensity index IM, impact duration T, and warning lead time W as inputs, and queries to obtain the loss rate for each computational grid. Multiplying this loss rate by the pre-gridized asset spatial distribution data yields the spatialized direct economic loss L. _direct (x, y). For indirect economic losses, the physical damage intensity of the location of critical infrastructure (such as substations and water plants) is mapped to its functional failure probability. Next, these initial failure probabilities are input into a cascading failure model based on a graph neural network (GNN). This model simulates the cascading failure process of systems such as power, water supply, and transportation by message passing on the infrastructure dependency graph, predicting the cumulative indirect losses L. _indirect (t). Furthermore, by integrating direct and indirect losses and considering uncertainty analysis (e.g., through Monte Carlo simulation), a comprehensive disaster loss assessment report is generated that includes the spatial distribution, temporal evolution, and confidence intervals of the losses.

[0101] Another example utilizes the spatialized intensity index field IM, impact duration T, and early warning lead time W. The loss rate of each grid is queried using a pre-constructed three-dimensional vulnerability surface V(IM, T, W). Combined with spatial distribution data of socio-economic assets, the gridded direct economic loss L is calculated. _direct (x, y). Furthermore, the physical damage level of critical infrastructure nodes is mapped to the probability of functional failure, input into a cascade model based on a graph neural network, and a message-passing mechanism is used to simulate the cascading failure process of electricity → water supply → transportation → healthcare, predicting the time-varying cascade failure sequence and the indirect loss accumulation curve L. _indirect (t). Integrate direct and indirect losses to generate a comprehensive disaster loss assessment report that includes spatial distribution, temporal evolution, and uncertainty intervals.

[0102] Another example involves constructing a vulnerability query function V = f(IM, T, W), where IM is the adjusted intensity field, T is the impact duration, and W is the warning lead time (obtained from the emergency response system). The vulnerability surface is interpolated using radial basis functions: V(IM, T, W) = Σ _i λ _i *φ(||[IM,T,W]-[IM _i T _i W _i ]||), where φ is the Gaussian kernel function, and the parameters are calibrated using a historical disaster loss case database. The obtained loss rate V(x, y) is compared with the asset value density map A(x, y) (yuan / m²). 2 Multiply by , and calculate the direct loss L of the gridding. _direct(x, y) = V*A, and aggregated to obtain the categories of direct losses: agricultural losses, housing losses, infrastructure losses, etc.

[0103] Extracting the local intensity value (IM) of key facility locations (substations, water plants, hospitals, etc.) _facility Through the facility's specific vulnerability function P _fail =1 / (1+exp(-α(IM-IM 50 )))Calculate the failure probability, where IM 50 The strength threshold corresponding to 50% failure is given, and α is the curve steepness parameter. The continuous probability is discretized into failure states: P < 0.3 (normal), 0.3 ≤ P < 0.7 (damaged), P ≥ 0.7 (failed), generating the facility state vector S. _facility .

[0104] Constructing the infrastructure dependency graph G _infra Nodes represent facilities, and edges represent dependencies (e.g., hospital → electricity, water supply). The facility state vector S... _facility As the initial state, run the graph neural network cascade model: each node has a state S' _i =max(S _i , Σ _j w _ij *S _j *δ _ij ), where w _ij For weights, δ _ij The propagation coefficients are used. Iterative updates are performed until convergence, outputting the cascading failure time sequence {S(t1), S(t2), ...}, recording key time points: the first large-scale power outage t... _power Water supply interruption _water wait.

[0105] Based on the cascading failure sequence, an input-output model is used to calculate the supply chain disruption loss: L _indirect =A*(I-(ID)*A) -1 *Y, where A is the input coefficient matrix, D is the diagonal matrix of capacity reduction due to facility failure, and Y is the final demand vector. The cumulative indirect loss time curve L is calculated. _indirect (t). Integration of direct losses L _direct and indirect losses L _indirect Considering the recovery time factor, calculate the total economic loss L. _total =L _direct +∫L _indirect (t)*e -λt dt. By propagating parameter uncertainty through Monte Carlo simulation (n=1000), a 90% confidence interval [L] for the loss is generated. _low L _highThe final output includes a comprehensive disaster loss assessment report that incorporates spatial distribution, temporal evolution, and classification statistics.

[0106] By calculating the ultra-high margin change rate dΔH / dt in real time and performing linear extrapolation t _predict =t _now The method, using -ΔH / (dΔH / dt), upgrades traditional static water level monitoring to dynamic trend prediction, enabling early warnings of overtopping risks 2-4 hours in advance. This method incorporates an automatic warning escalation mechanism triggered by a rate-of-change threshold (-0.1 m / h). When water levels rise sharply, the system automatically upgrades the warning level from yellow to orange or red, avoiding delays caused by manual judgment. In practical applications, this dynamic warning mechanism extends the average evacuation preparation time in downstream hazardous areas, improving evacuation effectiveness, especially for elderly people with mobility issues, demonstrating the impact of timely warnings on disaster reduction.

[0107] Example 2 provides an optional scheme for simulating breach evolution based on the cumulative scour energy criterion, especially the implementation process of the overtopping-breach phase transition simulation, as follows:

[0108] Step 2.1: For each potential overtopping area, determine the instantaneous scour power of the overtopping flow. This determination includes: analyzing the hydraulic conditions of the potential overtopping area to obtain the overtopping flow unit width discharge; extracting the local topographic slope of the potential overtopping area from hydrogeographic data; calculating the Froude number based on the overtopping flow unit width discharge, and generating a kinetic energy correction factor based on this Froude number; and comprehensively considering the overtopping flow unit width discharge, local topographic slope, and kinetic energy correction factor to determine the instantaneous scour power.

[0109] Alternatively, the instantaneous unit width scouring power P(t) can be calculated using the following formula: P(t) = ρ * g * q(t) * |▽z| _smooth *κ(Fr(t)); where P(t) is the instantaneous unit width scouring power at time t, in watts per meter (W / m); ρ is the density of water, usually taken as 1000 kg / m³. 3 g is the acceleration due to gravity, usually taken as 9.8 m / s². 2 q(t) represents the instantaneous overhead unit width flow rate at time t, in meters per second (m). 2 / s, which can be expressed by the improved broad-crested weir formula q(t)=Cd*(2g). 0.5 *h _over (t) 1.5 The calculation yields h. _over (t) represents the head of water flowing over the flood at time t, and Cd is the flow coefficient. Cd can be dynamically adjusted according to the degree of flooding to improve accuracy; |▽z| _smoothThe local terrain slope is smoothed and is a dimensionless pure number. It can be calculated from the DEM data within a 3×3 or 5×5 window centered on the potential overtopping area using the Horn algorithm, and then Gaussian smoothed to eliminate the influence of outliers. κ(Fr(t)) is the kinetic energy amplification factor based on the Froude number at time t, also a dimensionless pure number. It is used to correct the sharp increase in kinetic energy when the flow transitions from subcritical to supercritical, and can be expressed as κ(Fr) = 1 + 0.5 * Fr(t). 2 The function *tanh(Fr(t)-1) is used for calculation; Fr(t) is the Froude number at time t, Fr(t)=v _over (t) / (g*h _flow (t)) 0.5 , where v _over (t) and h _flow (t) represents the average velocity and water depth of the overpass flow, respectively.

[0110] In this embodiment, when using this formula to calculate the instantaneous scour power, the unit width flow rate q(t) characterizing the scale of the water flow, the topographic dynamics |▽z| that causes scour, and the Froude number Fr(t) characterizing the internal energy state of the water flow are coupled to form a physical quantity with clear physical meaning and high integration, which can be used to characterize the destructive power of water flow on dikes.

[0111] Step 2.2: Integrate the instantaneous scouring power over time to obtain the cumulative scouring energy; and set the anti-scouring energy threshold.

[0112] For example, the cumulative flushing energy W _cum (t) can be obtained by integrating the instantaneous scouring power P(t) from 0 to time t, for example, using the trapezoidal integration method W. _cum (t)=W _cum (t-Δt)+0.5*[P(t)+P(t-Δt)]*Δt.

[0113] Setting a scour resistance energy threshold specifically includes: extracting the physical properties of the dike material in the potential overtopping area from dike monitoring data and INSAR-based dike deformation radar remote sensing results; the physical properties include soil cohesion and internal friction angle; calculating the critical shear stress of the dike material based on the soil cohesion and internal friction angle; and converting the critical shear stress into a scour resistance energy threshold. Alternatively, setting a scour resistance energy threshold may include: extracting the physical properties of the dike material in the potential overtopping area from hydrogeological data, including at least soil cohesion and internal friction angle; calculating the critical shear stress of the dike material based on these (soil cohesion and internal friction angle); and converting the critical shear stress into a scour resistance energy threshold.

[0114] Optionally, soil physical parameters for the corresponding levee section can be extracted from the levee database, such as cohesion c (in kPa) and internal friction angle φ (in °); the critical shear stress τ can then be calculated. _c =c+σ _n *tan(φ), where σ _n The effective normal stress can be estimated based on the soil's dry density and moisture content; this critical shear stress is then determined by W. _cr =τ _c *L _ref *T _ref Converting / (ρ*g) into the critical impact resistance threshold W _cr (i.e., impact energy threshold), where L _ref and T _ref These are empirical reference scour length and reference scour duration, for example, 10 meters and 300 seconds respectively.

[0115] The criterion for breach occurrence is changed from an instantaneous mechanical threshold (critical shear stress) to a process-based energy threshold (critical scour resistance energy), which is related to the fatigue and cumulative damage process of the dike material under continuous scour, and is used to determine the occurrence of breach.

[0116] Another example is the calculation and cumulative monitoring of overhead flow power, which can also be achieved through the following method:

[0117] From the dangerous section D of the flood peak _risk Extract real-time water level H(t) and dike crest elevation Z _levee Calculate the head of water above the top, h _over (t)=H(t)-Z _levee When h _over When the value is greater than 0, the instantaneous unit width flow rate q(t) = Cd*(2g) is calculated using the improved broad-crested weir formula. 0.5 *h _over 1.5 *ψ(θ), where the flow coefficient Cd=0.385*(1+h _over / P) 0.5 The dynamic adjustment is based on the degree of inundation, where P is the weir height, and the angle correction coefficient ψ(θ)=cos(θ / 2) takes into account the influence of the water inflow angle θ. The calculated time-varying single-width flow sequence q(t) is stored at 1-minute intervals for power calculation.

[0118] Using the location coordinates of the breach (x) _b y _bCentered on a 30m DEM, a 3×3 window elevation matrix is ​​extracted from the data. The Horn algorithm is used to calculate the slope components: Яz / Яx = [(z3-z1) + 2(z6-z4) + (z9-z7)] / (8*Δx), Яz / Яy = [(z1-z7) + 2(z2-z8) + (z3-z9)] / (8*Δy); where Я is the partial derivative. The local slope value |▽z| = sqrt((Яz / Яx)) is generated. 2 +(Яz / Яy) 2 ); Calculate the prevailing scour direction angle α = arctan(Яz / Яy, Яz / Яx) to determine the dominant scour direction; smooth the slope |▽z| _smooth =0.7*|▽z|+0.3*|▽z| _avg To avoid the influence of local outliers, where |▽z| _avg The average slope is for a 5×5 window.

[0119] Based on the instantaneous unit width flow rate q(t) and the estimated average depth of the overpass flow h _flow =(q / (1.5*sqrt(g*h _over ))) 2 / 3 Calculate the overpass velocity v _over (t)=q(t) / h _flow ; Calculate the Froude number Fr(t) = v _over / sqrt(g*h _flow The kinetic energy characteristics of water flow are characterized by κ(Fr), and the kinetic energy amplification factor κ(Fr) = 1 + 0.5 * Fr is constructed. 2 *tanh(Fr-1), where the tanh function is less modified in subcritical flow (Fr<1) and amplified in supercritical flow (Fr>1); Fr(t) and κ(Fr) are stored as time-varying energy correction sequences.

[0120] Combined unit width flow rate q(t), smooth slope |▽z| _smooth Given the kinetic energy amplification factor κ(Fr), calculate the instantaneous unit width scour power P(t) = ρ*g*q(t)*|▽z| _smooth *κ(Fr(t)), where ρ=1000kg / m 3 Establish a cumulative power integrator and use the trapezoidal integration method to calculate the cumulative scouring power W. _cum (t)=W _cum (t-Δt)+0.5*[P(t)+P(t-Δt)]*Δt, time step Δt=60 seconds; calculate the power change rate dP / dt=[P(t)-P(t-Δt)] / Δt, when dP / dt>0.1kW / m / s, mark it as the moment of rapid power increase, and generate the power time series triplet [P(t), W _cum[(t), dP / dt], are used for phase transition criteria.

[0121] Step 2.3: When the accumulated scour energy exceeds the scour resistance energy threshold, breach evolution is triggered. After breach evolution is triggered, the following steps are also included: calculating the amount by which the accumulated scour energy exceeds the scour resistance energy threshold to obtain the super-threshold scour energy; determining the initial breach geometric parameters based on the super-threshold scour energy; and simulating the time-varying process of the geometric morphology after breach triggering, starting from the initial breach geometric parameters. Alternatively, when the accumulated scour energy exceeds the scour resistance energy threshold, the breach evolution program is initiated. After breach evolution is initiated, the following steps are also included: calculating the portion of the accumulated scour energy exceeding the scour resistance energy threshold to obtain the super-threshold scour energy; determining the initial breach geometric parameters based on the super-threshold scour energy; and simulating the change in the geometric morphology over time after breach triggering based on the initial breach geometric parameters.

[0122] Furthermore, when W _cum (t)>W _cr At that time, the system triggers a breach evolution. Immediately calculate the excess scour energy ΔW = W. _cum (t)-W _cr ; where ΔW is related to the destructive energy exceeding the levee's erosion resistance limit. Based on this threshold energy ΔW, the initial breach geometry parameters (such as the initial depth h) are determined. _b0 and initial bottom width B _b0 ) can be determined by the following or similar relation: h _b0 =k _h *(P _breach / P _ref ) 0.6 *h _over B _b0 =k _b *h _b0 *(2+tan(φ)); where h _b0 B is the initial breach depth. _b0 P is the initial breach base width; _breach The scouring power at the moment of triggering the breach; h _over To trigger the instantaneous overhead; k _h and k _b It is related to the relative damage degree D _r =ΔW / W _cr Related dimensionless coefficients, such as k _h =0.3+0.5*tanh(D _r ); P _ref It is a reference power, for example, 100kW / m; φ is the internal friction angle of the soil.

[0123] In this step, a physical relationship is established between the initial morphology of the breach and the actual degree of overload damage to the levee (represented by ΔW), which reduces the reliance on empirical assumptions about the initial breach size and makes the initial conditions of the entire breach evolution simulation correspond to the actual overload damage state of the levee.

[0124] Another example is that the implementation can be as follows: extract the soil physical parameters from the embankment material database: cohesion c (kPa), internal friction angle φ (°), and dry density ρ. _d (g / cm 3 ), moisture content w (%). Calculate the critical shear stress τ. _c =c+σ _n *tanφ, where the effective stress σ _n =ρ _d *g*h _levee *(1-0.4w) Considering the effect of moisture content; converting the critical shear stress into the critical scour work threshold W. _cr =τ _c *L _ref *T _ref / (ρ*g), where the reference length L _ref =10m (typical initial breach length), reference time T _ref =300s (impact resistance duration); Introducing a safety factor SF = 0.7~0.9 (depending on the dike grade), the adjusted threshold W is obtained. _cr_adj =W _cr *SF.

[0125] Real-time comparison of cumulative flushing power (W) _cum (t) and adjustment threshold W _cr_adj Construct the main criterion J1=(W _cum >W _cr_adj Monitor the power change rate dP / dt, and construct an auxiliary criterion J2 = (dP / dt > 0.2 kW / m / s for 30 seconds) to characterize scour acceleration; check the local slope |▽z|, and construct a terrain criterion J3 = (|▽z| > 0.02) to ensure sufficient scour force; when J1*J3 = 1 or J2*J3 = 1, trigger the breach transition flag Flag_breach = 1, and record the transition time t. _breach and conversion instantaneous power P _breach .

[0126] Calculate the energy exceeding the threshold ΔW=W _cum (t _breach )-W _cr_adj This characterizes the degree to which the impact resistance is exceeded; it defines the relative damage degree D. _r =ΔW / W _cr_adj , range [0, ∞); establish damage-geometric mapping relationship: when D _r∈[0, 0.5] represents mild damage, D _r ∈[0.5, 1.5] represents moderate damage, D _r A value >1.5 indicates severe damage; the initial scour depth coefficient k is calculated based on the damage level. _h =0.3+0.5*tanh(D _r ) and lateral expansion coefficient k _b =0.5+0.8*tanh(0.5*D _r Output damage parameter set [D] _r k _h k _b ].

[0127] Based on damage parameter [k] _h k _b and instantaneous power P _breach Calculate the initial breach depth h _b0 =k _h *(P _breach / P _ref ) 0.6 *h _over , where P _ref =100kW / m is the reference power; initial bottom width B _b0 =k _b *h _b0 *[2+tan(φ)], where φ is the internal friction angle of the soil, which determines the natural collapse angle; initial slope ratio m _b0 =1 / tan(45°+φ / 2), following Coulomb's earth pressure theory; Breach bottom elevation Z _b0 =Z _levee -h _b0 The initial four parameters of the breach [B] will be used. _b0 h _b0 m _b0 Z _b0 Pack and output.

[0128] At the instant of phase transition triggering, the flow mode is switched from weir flow to breach flow; the initial breach flow rate Q at the instant of transition is calculated. _b0 =(2 / 3)*Cd*B _b0 *sqrt(2g)*[(H _up -Z _b0 ) 1.5 -(H _up -Z _levee ) 1.5 To ensure flow continuity, update the upstream boundary conditions, replacing the original weir outflow term with the time-varying breach outflow term. Record the mode transition log [t]. _breach Flag_breach, Q _b0It sends a breach evolution initiation signal to downstream modules.

[0129] Step 2.4: Simulate the time-varying process of the geometric morphology after the breach is triggered, and calculate the breach discharge process based on this time-varying process. Alternatively, this step can be implemented as follows: Simulate the time-varying process of the geometric morphology after the breach is triggered, coupled with the storage state variables within the river network channel, and calculate the breach discharge process based on this (time-varying process coupled with storage state variables within the river network channel). After determining the initial breach geometric parameters, the model enters the time-varying evolution stage. The simulated time-varying process of the geometric morphology is driven by the bottom shear stress; the bottom shear stress is obtained by applying a turbulence correction to the shear stress calculated based on the Manning formula.

[0130] Specifically, a system of differential equations can be established regarding the evolution of the breach width B(t) and depth h(t), such as dB / dt=f(τ). _b_turb ,...) and dh / dt=g(τ _b_turb The driving force of the equation is the bottom shear stress τ. _b_turb This shear stress can be controlled by τ. _b_turb =τ _b *(1+0.3*(v _b 2 / (g*R _h )-1) + The calculation yields τ, where τ is the τ value. _b It is the shear stress calculated based on the Manning formula, v _b and R _h These are the average flow velocity and hydraulic radius of the breach cross-section, (...) + The positive part of the function is represented by the equation; by solving this system of differential equations using adaptive Runge-Kutta numerical methods, the time-varying sequence of the breach geometry can be obtained. Correspondingly, at each time step, the breach flow rate Q is calculated in real-time using the weir-orifice flow combination formula based on the updated breach geometry. _b (t), generating the complete breach flow process.

[0131] Another possible implementation is: based on the current breach flow Q _b Given the t and the breach geometry [B(t), h(t), m(t)], calculate the flow area A. _b =(B+m*h)*h and wetted perimeter P _w =B+2h*sqrt(1+m 2 ), thus obtaining the hydraulic radius R _h =A _b / P _w The average flow velocity v was calculated using the Manning formula. _b =Q _b / A _b Calculate the bottom shear stress τ_b (t)=ρ*g*n 2 *v _b 2 / R _h 1 / 3 Where n is the erosion roughness (initially 0.035, decreasing to 0.025 as erosion smoothing occurs); turbulence correction τ is introduced. _b_turb =τ _b *(1+0.3*(v _b 2 / g / R _h -1) + Output the corrected shear stress sequence τ _b_turb (t).

[0132] Construct the vertical scour differential equation dh _b / dt=β _v *(τ _b_turb -τ _c ) 1.5 *exp(-γ _v *t / T _scale )*f _sat (h _b ), where the scouring coefficient β _v =5×10 -4 m*s -1 *Pa -1.5 Adjust the time decay coefficient γ according to the soil type. _v =0.5 reflects the gradual compaction of the soil, characteristic time T _scale =3600s, saturation function f _sat (h _b )=(1-h _b / h _max ) 0.5 Maximum drawing depth h _max =0.8*H _channel The fourth-order Runge-Kutta method is used to solve the problem, with a time step Δt. _r K=min(60s, 0.1*h) _b / |dh _b / dt|) Adaptive adjustment to generate the time history of the ulcer depth h _b (t) and bottom elevation time history Z _b (t)=Z _b0 -h _b (t).

[0133] Establish the lateral expansion equation dB / dt=β _h *τ _b_turb 1.5 *(1-B / B _max )*g_shape(B / h), where the transverse scour coefficient β_h =3×10 -4 m*s -1 *Pa -1.5 Maximum width B _max =min(0.3*L _reach (100m) subject to the length L of the river section _reach The shape function is defined as g_shape(B / h) = 1 + 2 * exp(-((B / h - 5) / 3) 2 Promotes the width-to-depth ratio to tend towards a stable value of 5; slope evolution dm / dt=-0.01*(mm) _stable )*(τ _b_turb / τ _c -1) + This causes the slope to gradually stabilize with a slope ratio of m. _stable =2.0; Parallel solution generates the breach width time history B(t) and slope time history m(t).

[0134] Based on the updated breach geometry [B(t), h(t), m(t), Z] _b (t)] and upstream water level H _up (t), the breach flow rate is calculated using the broad-crested weir-orifice combination formula: when H _up >Z _b At +1.5h (submerged outflow), Q _b (t)=C _sub *A _b *sqrt(2g*(H _up -Z _b -0.5h), submersion coefficient C _sub =0.65; when H _up ≤Z _b At +1.5h (free outflow), Q _b (t)=(2 / 3)*C _free *B _eff *sqrt(2g)*(H _up -Z _b ) 1.5 Free flow coefficient C _free =0.45, effective width B _eff =B+2m*h / 3; Perform mass conservation check: |Q _in -Q _out -dV / dt|<0.01*Q _in Where V is the control volume, Q _in Q represents the total flow rate into the control body. _out The total flow rate out of the control body (including the initially calculated breach flow rate Q) _b ), dV / dt is the rate of change of the water volume within the control body; ensure the breach flow rate Q after verification._b_verified (t) satisfies continuity.

[0135] The verified breach flow rate Q _b_verified (t) is fed back as a lateral outflow boundary condition to the one-dimensional river network model, updating the momentum equation source term S of node i. _i =-Q _b_verified / A _i Update the mass equation: ЯA / Яt + ЯQ / Яx = -q _b , where q _b =Q _b_verified / Δx represents the flow rate per unit length of the breach; record the cumulative flow volume V. _total (t)=∫Q _b dt and peak flow Q _peak =max(Q _b Generate a complete data packet for the vulnerability evolution [t, B(t), h(t), Z]. _b (t), Q _b (t), V _total (t)], which is transmitted in real time to the downstream flooding simulation module.

[0136] In a specific embodiment, the breach flow rate Q is established. _b A two-way real-time coupling mechanism with the river network hydrodynamic model is introduced, along with a mass conservation check |Q _in -Q _out -dV / dt|<0.01*Q _in This mechanism ensures the conservation of water volume throughout the river network system during the breach evolution process. The breach discharge is fed back as a lateral source term to the one-dimensional river network model in real time, and changes in river network water level dynamically affect the breach discharge calculation, forming a coupled feedback loop. This is achieved through the bottom shear stress τ. _b The synchronous iteration of the driven breach geometric evolution equations dB / dt, dh / dt, and flow updates eliminates the water imbalance problem caused by the disconnect between traditional offline breach simulation and river network simulation, reduces the system mass conservation error, and improves the numerical stability and physical reliability of long-term simulations. In other words, this step is used for the technical implementation of water quantity, mass conservation, numerical stability, and physical correlation in river network system simulation.

[0137] Example 3 provides a possible implementation scheme for generating a spatialized disaster intensity field. After obtaining the breach flow process, a comprehensive intensity field that can fully reflect the disaster-causing capacity is generated by using a two-dimensional flood evolution simulation and a comprehensive intensity index coupled with multiple physical quantities.

[0138] Step 3.1: Based on the breach flow process and hydrogeographic data, calculate the spatiotemporal distribution of inundation depth and velocity field; based on the spatiotemporal distribution of inundation depth and velocity field, determine the comprehensive intensity index and impact duration for any calculation unit within the inundation area; the comprehensive intensity index and impact duration of each calculation unit constitute the spatialized comprehensive intensity index field and impact duration field, respectively, which together form the spatialized disaster intensity field; wherein, the inundation area is the spatial range determined based on the inundation depth being greater than the set depth threshold.

[0139] In other words, based on the breach flow process and hydrogeographic data, the spatiotemporal distribution calculations of inundation depth and velocity field are performed. Furthermore, the definition criteria for the inundation area are clarified as follows: the inundation depth is greater than a set depth threshold. Based on the calculated spatiotemporal distribution of inundation depth and velocity field, the comprehensive intensity index and impact duration are calculated for any calculation unit within the area. On this basis, the comprehensive intensity index and impact duration of all calculation units are integrated in spatial dimensions to generate a spatialized comprehensive intensity index field and a spatialized impact duration field, which together constitute a spatialized disaster intensity field.

[0140] Specifically, this step uses the breach flow process line Q _breach (t) is the starting point, which is used as the upstream boundary condition for the two-dimensional hydrodynamic model. Hydrogeographic data here primarily refers to high-resolution digital elevation models (DEMs), such as 30-meter resolution raster data. The model solves a set of two-dimensional shallow-water equations within the two-dimensional computational domain defined by the DEM to calculate flood evolution.

[0141] In this embodiment, the solution of the two-dimensional shallow water equations can be achieved using a numerical scheme based on the finite volume method. Specifically, when processing the interface flux of each computational grid cell, the HLLC approximate Riemann solver is preferably used, as it can effectively handle complex flow regimes such as shock waves and wet-dry boundaries that may occur during flood evolution. Furthermore, the time step Δt can be adaptively adjusted, for example, by setting Δt = CFL * min(Δx / sqrt(g * h)) based on the CFL condition, where CFL is the Courant number, which can be between 0.5 and 0.8, Δx is the grid size, and h is the water depth. For the handling of wet-dry boundaries, in other words, the delineation of the flooded and unflooded areas, a very small water depth threshold (e.g., h) can be introduced. _min =0.001 meters) to distinguish between wet and dry grids (unsubmerged and submerged areas); flow velocities in newly submerged grids are limited to avoid non-physical numerical oscillations. From this, the water depth time series h(x,y,t) and flow velocity vector time series v(x,y,t) on all computational grids (x,y) within the submerged area can be obtained. These data form the basis for subsequent calculations of the comprehensive intensity index.

[0142] Step 3.2: Determine the comprehensive intensity index for any computational unit within the inundation area, including: analyzing the water depth time series and flow velocity time series of the computational unit from the spatiotemporal distribution of the inundation water depth and flow velocity fields; extracting the local topographic slope of the computational unit from hydrogeographic data; coupling the water depth, flow velocity and local topographic slope for each inundation moment of the computational unit to determine the instantaneous impact power density; integrating the instantaneous impact power density along the time dimension to obtain the comprehensive intensity index of the computational unit.

[0143] Before coupling water depth, flow velocity, and local topographic slope to determine the instantaneous impact power density, the method further includes: identifying the water flow direction based on the flow velocity time series and identifying the slope direction based on the local topographic slope; determining the directional coupling coefficient based on the vector relationship between the water flow direction and the slope direction; modulating the local topographic slope using the directional coupling coefficient to generate an effective shaped slope; the coupling operation specifically involves coupling water depth, flow velocity, and effective shaped slope. In other embodiments, before coupling water depth, flow velocity, and local topographic slope to determine the instantaneous impact power density, the method further includes: identifying the water flow evolution direction based on the flow velocity time series and identifying the slope direction based on the local topographic slope; determining the directional coupling coefficient based on the vector relationship between the water flow evolution direction and the slope direction; modulating the local topographic slope using the directional coupling coefficient to generate an effective shaped slope; the coupling operation specifically involves coupling water depth, flow velocity, and effective shaped slope.

[0144] In this step, a comprehensive strength index IM will be constructed for each computation grid (i, j). _ij From the spatiotemporal distribution data generated in the previous step, the water depth time series h of this computing unit is extracted. _ij (t) and the time series of the velocity vector v _ij (t); Extract the local terrain slope vector ▽z of this cell from the DEM data. _ij Next, for each flooding time t, the instantaneous impact power density p is calculated. _ij (t), its formula can be described as: p _ij (t)=ρ*g*h _ij (t)*|v _ij (t)|*|▽z _eff_ij |*(1+|v _ij (t)| 2 / (2*g*h _ij (t)));where p _ij (t) represents the instantaneous impact power density of the calculation unit (i, j) at time t, in W / m³. 2 ρ and g are the density of water and the acceleration due to gravity, respectively; h _ij (t) and |v _ij (t)| represents the water depth and flow velocity of the unit at time t; |▽z _eff_ij|To effectively shape the slope; the term in parentheses (1+Fr) 2 / 2) is a correction for the impact kinetic energy of the water flow, where Fr is the Froude number.

[0145] Effectively form slope | ▽z _eff_ij The calculation method for | is: |▽z _eff_ij |=|▽z _ij |*max(C _min 1+C _weight *cos(θ _ij )); where, |▽z _ij | represents the original terrain slope of this unit; θ _ij Let v be the direction vector of water flow at time t. _ij (t) and the topographic slope direction vector ▽z _ij The included angle between them; cos(θ) _ij This is the directional coupling coefficient, with a value between [-1, 1]. It is positive when the water flows downhill and negative when it flows uphill; C _min and C _weight The adjustment coefficients, for example, can be 0.5 and 0.5 respectively, to retain the scouring effect of the foundation on the reverse slope and to adjust the enhancement effect on the downslope.

[0146] It can be seen that by introducing an effective slope design, the calculation of instantaneous impact power density is correlated with the interaction between water flow and topography; downslope flow (cos(θ) _ij The flow rate (cos(θ)) > 0 will have a stronger destructive force due to the accelerated conversion of gravitational potential energy, while the flow rate against the slope (cos(θ)) will have a stronger destructive force. _ij The destructive power of water flowing along the slope of the terrain (i.e., cos(θ) < 0) will be correspondingly reduced. _ij When θ > 0, the efficiency of converting the gravitational potential energy of the water flow into kinetic energy increases; when the water flows against the slope of the terrain (i.e., cos(θ) > 0), the efficiency of converting the gravitational potential energy of the water flow into kinetic energy increases; _ij When )<0), the efficiency of converting the gravitational potential energy of water flow into kinetic energy decreases, resulting in differences in the destructive effects of water flow in the two types of flow states.

[0147] Based on this, the instantaneous impact power density is integrated along the entire inundation time dimension to obtain the comprehensive intensity index IM of the calculation unit. _ij IM _ij =∫p _ij This integral (t)dt can be implemented using numerical integration methods, such as Simpson's integral. _ij The cumulative impact energy borne per unit area is measured in joules per meter. 2 (J / m 2Its calculation integrates four dimensions: inundation depth (potential energy), flow velocity (kinetic energy), topography (exacerbating effect), and duration of effect (cumulative effect), and is used to measure disaster intensity.

[0148] Step 3.3, after obtaining the comprehensive intensity index, also includes: identifying the land use type corresponding to the calculation unit; adjusting the comprehensive intensity index according to the land use type using a preset exposure adjustment coefficient to generate the adjusted intensity index.

[0149] After calculating the original comprehensive intensity index field IM(x, y), it is further spatially overlaid with the land use type layer. For each calculation unit (i, j), a preset exposure adjustment coefficient k is applied according to its corresponding land use type (such as farmland, built-up area, forest land, water body, etc.). _lu Its IM _ij The value is adjusted to generate the adjusted intensity index IM. _adj_ij =k _lu *IM _ij .

[0150] For example, for farmland areas where crops are susceptible to erosion and flooding damage, the adjustment coefficient k _lu The adjustment factor can be set to a value greater than 1, such as 1.2, to amplify the disaster intensity; for building areas with certain structural resistance, the adjustment factor can be set to a value less than 1, such as 0.8; for infrastructure such as roads, it can be set to 0.9. The adjustment factor can be calibrated based on historical disaster data or expert experience.

[0151] In this step, the physical disaster intensity index is transformed into a disaster risk intensity index, taking into account the differences in exposure and vulnerability of different disaster-bearing bodies (ground features). This allows the intensity field to reflect the potential degree of damage to the socio-economic system, providing input for subsequent loss assessment.

[0152] Example 4 provides an optional implementation scheme for a direct economic loss assessment method based on a three-dimensional vulnerability surface. It describes how to use a disaster intensity field, combined with a vulnerability model that includes a warning time dimension, to assess direct economic losses and quantify the uncertainty of the assessment results.

[0153] Step 4.1: The pre-constructed three-dimensional vulnerability surface specifically includes: obtaining a historical disaster case library, which contains the comprehensive strength index, impact duration, early warning period and corresponding actual loss rate of historical cases; training a radial basis function network based on the historical disaster case library to establish a nonlinear mapping relationship between the comprehensive strength index, impact duration, early warning period and loss rate; the trained radial basis function network constitutes the three-dimensional vulnerability surface.

[0154] In this embodiment, the construction of the three-dimensional vulnerability surface is an offline preprocessing process. Accordingly, detailed case data of historical flood disasters are collected and organized, and information on four key dimensions is extracted from the affected unit of each case: the comprehensive intensity index IM (which can be obtained by post-processing historical floods using the method of this invention), the impact duration T, the warning lead time W issued by the emergency response system, and the corresponding actual property loss rate L.

[0155] After obtaining the training dataset {IM _k T _k W _k L _k After that, the input data can be preprocessed. For example, since the range and distribution of data in each dimension vary greatly, non-linear standardization can be performed, such as IM. _norm =(IM-IM _mean ) / IM _std ;T _norm =log(T / T _median );W _norm =1-exp(-W / W _char ), where W _char The warning time can be set as a feature, such as 2 hours. For sparse data areas, data augmentation methods such as Bootstrap resampling can also be used.

[0156] Furthermore, a radial basis function (RBF) network is used to fit this set of high-dimensional nonlinear relationships. The RBF network can be in the form of V(x) = Σ _i λ _i *φ(||xc _i || / r _i ), where x is generated by IM _norm T _norm W _norm The resulting three-dimensional input vector, V(x), represents the predicted loss rate, φ is typically a Gaussian kernel function, and c... _i r _i and λ _i These are the network centers, radial parameters, and weights, respectively. The network training process includes: determining the RBF centers using clustering algorithms such as K-means. _i The location of the center; the radial parameter r of each center is optimized through methods such as cross-validation. _i The weight λ is determined by solving a system of linear equations. _i The trained RBF network, which numerically defines a mapping function from (IM, T, W) to the loss rate L, constitutes a three-dimensional vulnerability surface.

[0157] Alternatively, in addition to the RBF network, other machine learning models capable of handling nonlinear regression problems can be used to construct this surface, such as Support Vector Regression (SVR), Gaussian Process Regression (GPR), or Deep Neural Networks (DNN). Introducing the early warning period W as an independent dimension couples the impact of social emergency response capabilities into the loss assessment process, moving beyond mere physical process simulation to support the correlation between the assessment and real-world scenarios.

[0158] Step 4.2, assessing direct economic losses includes: for any calculation unit within the inundation area, extracting its comprehensive intensity index and impact duration from the spatialized disaster intensity field; obtaining the early warning lead time contained in the socio-economic data; applying a pre-constructed three-dimensional vulnerability surface, querying and determining the loss rate of the calculation unit based on the comprehensive intensity index, impact duration, and early warning lead time; and calculating the direct economic loss of the calculation unit by combining the asset spatial distribution data and loss rate contained in the socio-economic data.

[0159] During the real-time assessment phase, for any calculation unit (i, j) within the flooded area, its adjusted comprehensive intensity index IM is extracted. _adj_ij and impact duration T _ij Accordingly, the system early warning lead time W for this event is obtained from the emergency management system. _system These three values ​​(IM) _adj_ij T _ij W _system Construct the query vector x _query Perform the same standardization process as during training. Input the standardized query vector into the trained RBF network to obtain the point estimation loss rate V of that computational unit. _ij .

[0160] Based on this, the loss rate is overlaid with a layer of asset spatial distribution data from the socioeconomic data. Specifically, L _direct_ij =V _ij *A _ij A _ij This represents the asset value (e.g., yuan / square meter) within the calculation unit (i, j). The total direct economic loss L is obtained by summing the direct losses of all flooded units. _direct =ΣL _direct_ij .

[0161] Step 4.3, after determining the loss rate of the calculation unit, also includes: assessing the prediction uncertainty of the three-dimensional vulnerability surface at the query point and generating a confidence interval for the loss rate; assessing the prediction uncertainty specifically includes: calculating the Jacobian matrix of the three-dimensional vulnerability surface at the query point with respect to its input dimension; combining the Jacobian matrix with the uncertainty covariance of the input data to calculate the prediction variance of the loss rate, and generating a confidence interval based on the prediction variance.

[0162] In some embodiments, the prediction uncertainty of the loss rate is further quantified. Specifically, after obtaining the point-estimated loss rate V... _ij Simultaneously, the RBF network is calculated at query point x. _query The partial derivatives of the loss rate with respect to its three input dimensions (IM, T, W) form the Jacobian matrix J = [ЯV / ЯIM, ЯV / ЯT, ЯV / ЯW]; this matrix characterizes the sensitivity of the loss rate to changes in each input factor. In other words, this matrix is ​​used to correlate the loss rate with the response to changes in each input factor.

[0163] Assuming the input data itself has uncertainty, this uncertainty can be represented by the covariance matrix Σ. _x To describe (if the inputs are independent, then it is a diagonal matrix); according to error propagation theory, the prediction variance σ of the output loss rate _v 2 It can be done through σ _v 2 =J*Σ _x *J T To approximate the calculation. After obtaining the prediction variance σ... _v 2 Then, confidence intervals for the loss rate can be constructed. For example, a 90% confidence interval can be represented as [V...]. _ij -1.645*σ _v V _ij +1.645*σ _v Multiplying the upper and lower limits of this confidence interval by the asset value respectively yields the confidence interval for direct economic loss [L]. _lower L _upper ].

[0164] In this embodiment, a three-dimensional vulnerability surface V(IM, T, W) with a warning time dimension is introduced. This surface is constructed using the comprehensive intensity index IM, impact duration T, and warning lead time W as inputs, addressing the issues of static disaster loss assessment models and their failure to couple the dynamic social factors such as human emergency response. Through a radial basis function network, the model learns and fits complex nonlinear relationships in historical disaster data. The comprehensive intensity index IM, impact duration T, and warning lead time W are specific data directly related to the business rules of emergency management processes. When assessing losses, the model uses both physical disaster intensity (IM, T) and social disaster mitigation capacity (W) as variables determining the loss rate. In practical applications, the assessment results can quantify the disaster mitigation benefits brought by the warning system. For example, the model can compare the economic losses under two scenarios with warning lead times of 6 hours and 1 hour, showing that even with the same physical impact, the economic losses will differ. This method enables managers not only to know where disasters will occur, but also to assess the extent to which early warnings can reduce losses. It provides quantifiable decision-making basis for the construction investment of early warning systems and the optimization of emergency plans, and integrates physical process simulation with social management practices.

[0165] Example 5 provides an optional technical solution for assessing cascading failures and indirect losses of infrastructure based on graph neural networks. It describes how, after assessing direct economic losses, indirect economic losses can be further assessed by simulating the chain reaction of critical infrastructure networks.

[0166] Step 5.1: For each critical infrastructure node identified in the socio-economic data, determine its initial failure state based on the disaster intensity corresponding to the node extracted from the spatialized disaster intensity field; construct an infrastructure dependency graph representing the interdependencies between critical infrastructure nodes.

[0167] Specifically, critical infrastructure nodes may include, but are not limited to, substations, water plants, communication base stations, transportation hubs, and hospitals. The geographical coordinates of these nodes are overlaid with the spatialized disaster intensity field generated in Example 3 to extract the disaster intensity information for each facility node f, such as the adjusted comprehensive intensity index IM. _f Maximum submersion depth h _f and maximum flow velocity v _f .

[0168] In this embodiment, instead of relying solely on a single strength index, a multi-factor failure probability comprehensive evaluation model is used to determine the initial failure probability P. _fail For example: P _fail =1-(1-P _iM )*(1-P _h ), where P _iMThe failure probability, P, is caused by the comprehensive strength index IM. _h The failure probability is caused by the submerged water depth h, and both can be expressed using a vulnerability function specific to this type of facility (such as the Logistic function P=1 / (1+exp(-α(XX)). 50 The calculation yielded this result. This multi-factor model can more comprehensively reflect the damage mechanism of the facility. Accordingly, the successive failure probabilities P are calculated. _fail Discretize into initial failure state S _f (0), for example, when P _fail When P < 0.3, the state is normal (S=0); when 0.3 ≤ P _fail When P < 0.7, the state is damaged (S = 0.5); when P _fail When the value is ≥0.7, the state is invalid (S=1).

[0169] Furthermore, construct the infrastructure dependency graph G. _infra This graph is a multi-layered network, where nodes represent critical infrastructure and edges represent dependencies between them. The edge weight w _ij The dependency from node j to node i can be set according to the dependency type, for example: strong physical dependency (such as a water plant depending on a specific substation for power supply), with a weight w. _phy It can be set to 0.9; the weight w for geographical proximity dependence (e.g., two facilities within 1 kilometer may affect each other due to road interruption) _geo It can be set to 0.3; the weight w represents information or functional dependencies (such as hospitals relying on communication networks). _inf It can be set to 0.6.

[0170] Step 5.2: Starting from the initial failure state, simulate the failure propagation process on the infrastructure dependency graph to generate a cascading failure sequence. The simulation of the failure propagation process is specifically implemented by applying a graph neural network model. This implementation includes: on the topology of the infrastructure dependency graph, based on the initial failure state, performing iterative message passing operations to update the state of each critical infrastructure node; the message passing operations include: aggregating the state information received by any node from its neighboring nodes to form an aggregated message, and updating its state at the next moment based on the aggregated message and the node's current state; the cascading failure sequence consists of the state evolution history of the nodes during the iteration process.

[0171] In this step, the initial failure state vector S(0) is used as input in the dependency graph G. _infra Run a pre-trained GNN model on it to simulate the propagation of failure states.

[0172] Specifically, the iterative process of GNN is as follows: At each time step t, for any node i in the network, its state S at the next time step... _i(t+1) is calculated using the following or a similar message passing and state update mechanism:

[0173] Node i receives its current failure state S from all its dependent neighbor nodes j (j→i). _j (t), and based on the dependency weight w _ij Perform weighted aggregation to form aggregated message m _i (t)=Σ _j w _ij *S _j (t). The new state S of node i _i (t+1) is determined by its current state S _i (t) and the received aggregation message m _i (t) are jointly determined. The specific update function can be S. _i (t+1)=σ(W _1 *S _i (t)+W _2 *m _i (t)+b), where W _1 W _2 b are the parameters learned by the GNN network, and σ is the activation function (such as ReLU or Sigmoid) to keep its state values ​​within a reasonable range.

[0174] Furthermore, to make the simulation more realistic, more complex update logic can be introduced, such as: S _i (t+1)=min(1,S) _i (t)+Σ _j w _ij *T _ij *S _j (t)), where T _ij This is the propagation coefficient, which can be related to dependency type or distance. Correspondingly, a recovery mechanism can be introduced; for example, if node i is in a damaged state, but all its upstream dependent nodes have recovered, then its state S... _i (t) has a certain probability of decreasing.

[0175] The above iterative process will continue until the state vector S(t) of the entire network converges (i.e., ||S(t+1)-S(t)||2 is less than a very small threshold ε) or the preset maximum number of iterations is reached. The state evolution history matrix S[i,t] of all nodes during the entire iteration process constitutes the cascading failure sequence. Using a GNN for simulation can automatically learn and capture high-order, nonlinear dependencies in the network topology, discover potential failure propagation paths that traditional rule-based models may overlook, and improve the fidelity of the simulation.

[0176] Step 5.3: Calculate indirect economic losses based on the cascading failure sequence.

[0177] After obtaining the cascading failure sequence S[i, t], it is transformed into economic loss. Specifically, the functional loss of each infrastructure node i can be associated with its industry sector. At each time t, the functional level of node i can be represented as 1-S _i (t). This can be mapped to the capacity decline coefficient of the corresponding industry sector, forming a time-varying diagonal matrix D(t) of capacity decline. Further, this matrix D(t) is substituted into the national economic input-output model to calculate the indirect economic losses caused by supply chain disruptions. For example, the indirect loss L... _indirect (t) can be expressed by formula L _indirect (t)=A*(I-(ID(t))*A) -1 *Y is calculated, where A is the direct consumption coefficient matrix in the input-output table, Y is the final demand vector, and I is the identity matrix. For L... _indirect (t) By integrating over the entire impact period and considering the time value of money, the total indirect economic loss can be obtained.

[0178] Run a graph neural network model on a pre-built infrastructure dependency graph, using a message passing mechanism S. _i (t+1)=min(1,S) _i (t)+Σ _j w _ij *T _ij *S _j This method simulates the propagation of failure states, enabling dynamic prediction from single-point failure to network cascading collapse. By distinguishing between three types of relationships—physical dependence (weight 0.9), information dependence (weight 0.6), and geographical dependence (weight 0.3 / d)—it accurately characterizes the cascading failure path from electricity to water supply to transportation to healthcare, extending traditional static loss assessment to dynamic failure sequence prediction that includes a time dimension. In urban flooding scenarios, this model can predict the operational status of the water supply system within 2-6 hours after a substation failure and the service capacity of the healthcare system within 12 hours, providing quantitative evidence for the time-series allocation of emergency resources and improving the completeness of indirect loss assessment. In other words, this model provides data support for the time-series allocation of emergency resources and provides correlation parameters for indirect loss assessment.

[0179] Example 6 provides a specific numerical case for calculating breach triggering and disaster intensity. Exemplarily, it includes:

[0180] Assume the physical properties of the soil at a certain embankment cross section are as follows: cohesion c = 15 kPa, internal friction angle φ = 25°, and effective unit weight γ' = 18 kN / m. 3 During a flood event, the overtopping depth h at this section was monitored. _over(t) varies with time as summarized later. The terrain slope |▽z| behind this section is 0.03. Calculate the impact energy threshold W. _cr Assuming the effective normal stress σ at the top of the dike _n ≈γ'*h _levee / 2≈18kN / m 3 *5m / 2 = 45kPa. Critical shear stress τ _c =c+σ _n *tan(φ) = 15kPa + 45kPa *tan(25°) ≈ 35.98kPa. Impact energy threshold W _cr =τ _c *L _ref *T _ref / (ρ*g)≈35980Pa*10m*300s / (1000kg / m 3 *9.8m / s 2 )≈11014kJ / m. Calculate the cumulative scouring energy W. _cum (t), calculate each physical quantity according to the time step, and summarize as follows:

[0181] Time t (minutes), water depth at the top of the roof h _over (m), unit width flow rate q(t)(m) 2 / s), unit width power P(t)(kW / m), cumulative unit width energy W _cum (t)(kJ / m)}={(0,0.20,0.137,7.0,0),(10,0.30,0.253,12.9,5970(Calculation formula: (7.0+12.9) / 2*600),(20,0.45,0.458,23.4,16860(Calculation formula: 5970+(12.9+23.4) / 2*600),(30,0.50,0.540,27.6,32160(Calculation formula: 16860+(23.4+27.6) / 2*600))}; where q(t) is expressed as q=0.385*(2g) 0.5 *h _over 1.5 Calculate; v _over The estimate is (2g*h) _over ) 0.5 Fr is estimated as v _over / (g*h _over ) 0.5 =2 0.5 ≈1.414 (supercritical flow); κ(Fr)≈1+0.5*1.414 2 *tanh(0.414)≈1.39; P(t)=ρ*g*q*|▽z|*κ(Fr); W _cum(t) is accumulated through trapezoidal integration. At t=20 minutes, the accumulated scouring energy W _cum (16860kJ / m) for the first time exceeded the impact energy threshold W. _cr (11014 kJ / m). At this point, the breach is triggered. Calculate the initial breach geometry parameters: the superthreshold energy ΔW = W. _cum (20min)-W _cr =16860-11014=5846kJ / m. Relative damage degree D _r =ΔW / W _cr =5846 / 11014≈0.53. Depth coefficient k _h =0.3 + 0.5 * tanh(0.53) ≈ 0.543. Trigger instantaneous power P _breach =23.4kW / m; Triggering instantaneous head h _over =0.45m. Initial breach depth h _b0 =k _h *(P _breach / P _ref ) 0.6 *h _over ≈0.543*(23.4 / 100) 0.6 *0.45≈0.11m.

[0182] After the breach, the floodwaters flowed downstream. The topographic slope of a downstream computational cell (i, j) is |▽z| = 0.02. The simulated hydrodynamic time series for this cell is as follows: {time t (minutes), h(t)(m), v(t)(m / s), p(t)(W / m)} 2 ), cumulative impact energy IM (kJ / m 2 )}={(0,0.0,0.0,0,0,(5,0.5,1.5,291.5,43.7(10,1.2,2.5,1342.3,288.8(15,0.8,1.8,611.7,581.9(20,0.2,0.5,29.8,678.0(25,0.0,0.0,0,682.5)}.

[0183] Assuming the angle θ between the water flow direction and the terrain slope direction is constant at 30 degrees, the directional coupling coefficient cos(θ) ≈ 0.866. Calculate the effective terrain slope |▽z. _eff ||▽z _eff |=|▽z|*max(0.5, 1+0.5*cos(θ))=0.02*max(0.5, 1+0.5*0.866)≈0.0287. Calculate the comprehensive strength index IM. _ij The instantaneous impact power density p(t) is calculated based on the time step: p(t) = ρ * g * h * v * |▽z_eff |*(1+v 2 Integrate / (2gh)).

[0184] Based on the above, the final comprehensive strength index IM of this calculation unit _ij The value is 682.5 kJ / m 2 .

[0185] This invention introduces a comprehensive intensity index (IM), addressing the problem that the simplistic measurement of disaster intensity fails to fully reflect the actual destructive power of floods. Specifically, for any computational unit within the inundation zone, the model, at each time step, couples water depth h, flow velocity v, and effective slope |▽z. _eff |, calculate the instantaneous impact power density p(t); where, the effective slope |▽z _eff The calculation also incorporates the coupling coefficient cos(θ) between the flow direction and the slope direction, enabling it to distinguish between different physical effects such as downslope scouring and uplift. Next, the instantaneous power density over the entire inundation period is integrated to obtain the comprehensive intensity index IM. In specific disaster assessment scenarios, for two types of floods with the same water depth but significantly different flow velocities (such as a breach breach rapids area and a polder calm water area), the IM index can distinguish their differences in destructive force, avoiding the potential underestimation of risk in high-velocity areas. This transforms disaster risk assessment from a vague concept of water depth into an energy metric, providing input for subsequent vulnerability analysis and loss assessment, and is particularly suitable for disaster-bearing structures such as buildings and crops that are highly sensitive to the impact of water flow.

[0186] Example 7 provides optional optimization methods for the stability and efficiency of hydrodynamic simulation and data processing, which can be used to improve computational efficiency, numerical stability and data management capabilities.

[0187] As a preferred implementation of the hydrodynamic analysis steps, to shorten the computation time while ensuring accuracy, the solution process of the nonlinear equation system can be optimized. Specifically, when using the Newton-Raphson method to iteratively solve X... k+1 =X k -J -1 *F(X k When ), AitkenΔ can be introduced. 2 Acceleration methods are used to improve convergence speed. This is achieved by performing three consecutive iterations X... k X k+1 X k+2 Based on this, construct an accelerated solution X. _acc X _acc =X k -(ΔX k ) 2 / (Δ 2 Xk ); where ΔX k =X k+1 -X k It is a first-order difference; Δ 2 X k =X k+2 -2*X k+1 +X k This is a second-order difference. By periodically applying Aitken acceleration (e.g., after every 3-5 regular iterations), the convergence criterion ||F(X)||<10 can be reached faster with fewer iterations. -6 .

[0188] Furthermore, in river sections where the flow state is close to the critical flow (i.e., Froude number Fr≈1), the Newton-Raphson iteration is prone to numerical oscillations or even non-convergence. Therefore, a relaxation factor ω can be introduced to modify the iteration step size. That is, X k+1 =X k -ω*J -1 *F(X k ); where ω is the relaxation factor, whose value ranges between (0, 1). For example, it can be dynamically adjusted according to the real-time calculated Fr number. When Fr is close to 1, ω takes a smaller value, such as 0.7, to suppress oscillations and ensure the stability of the solution. Through the above combination optimization, the hydrodynamic simulation module of this invention can achieve fast and robust solutions under complex water networks and extreme hydraulic conditions, providing input for subsequent overtopping judgment and breach simulation.

[0189] As a preferred implementation for generating water depth and velocity fields, this embodiment provides a spatiotemporal data cube management strategy to address the storage and query performance issues arising from the massive spatiotemporal data (up to TB levels) generated by two-dimensional flood evolution simulation. Specifically, in terms of data storage, a hybrid compression scheme combining incremental and wavelet methods is adopted. Specifically, instead of storing the complete spatial field data for every time step, the complete two-dimensional water depth field h(x, y, t) is stored only at key time nodes (e.g., every hour on the hour). _key ) and velocity field v(x, y, t) _key These keyframe data are compressed using two-dimensional discrete wavelet transform, achieving a high compression ratio by retaining over 95% of the energy coefficients. At intermediate moments between two key nodes, only the difference between the current and previous moments is stored, i.e., the incremental data Δh(t) = h(t) - h(t-1) and Δv(t) = v(t) - v(t-1). Since incremental data typically exhibits high sparsity, lossless compression algorithms such as run-length encoding or Huffman coding can be used for efficient compression. Correspondingly, this hybrid compression strategy can reduce the storage requirements of the original data by orders of magnitude (compression ratio up to 10:1) while ensuring that the reconstruction error remains within an acceptable range.

[0190] In terms of data querying, a multi-level index structure is constructed to support fast access to compressed data. Specifically, in the time dimension, a B index structure can be built. + Tree indexes map timestamps to their physical locations within storage files, supporting fast time-slice queries. Spatially, a quadtree index can be built for each keyframe, mapping spatial coordinate regions to corresponding data blocks, supporting efficient spatial range queries. When querying data at any spatiotemporal point (x, y, t), the system uses a tree index... + The tree locates the nearest keyframe t before t. _key Extract the data block containing (x, y) from the quadtree, decompress it, and then sequentially read and accumulate it from t. _key By using all the incremental data up to time t, the required data can be quickly reconstructed.

[0191] Based on this, the present invention can handle big data-related problems in complex simulations, enabling efficient execution of operations that require frequent access to spatiotemporal data, such as the calculation of the comprehensive intensity index (IM) and the visualization analysis of the flooding process. This strategy can be used to improve the operational efficiency and response speed of the entire system.

[0192] Example 8: By using an optimized judgment model for assessing breach triggering and facility failure, the model's decisions are made more closely in line with physical and social realities by incorporating more multi-dimensional information.

[0193] As a preferred implementation of the breach triggering judgment step, this embodiment provides a composite triggering logic that integrates energy accumulation, scour intensity, and terrain conditions to replace the single accumulated scour energy threshold criterion.

[0194] Specifically, the mechanism monitors three independent criteria in parallel: the primary criterion for energy accumulation J1: J1=(W _cum >W _cr_adj );W _cr_adj This is the adjusted critical scour resistance threshold. This criterion is mainly used to identify fatigue failure of dikes caused by long-duration, moderate-intensity scour. Scour Intensity Auxiliary Criterion J2: This criterion is used to capture situations where scour intensity increases sharply in a short period of time; its expression is J2=(dP / dt>P). _rate_threshold ); where dP / dt is the rate of change of instantaneous scouring power, P _rate_threshold A preset power surge threshold, such as 0.2 kW / m / s, is set, and this state must persist for a certain period (e.g., 30 seconds) to exclude fluctuations in power. This criterion is particularly effective in identifying rapid damage caused by sudden, high-intensity water flow resulting from upstream dam breaches or extreme rainfall. Topographic Dynamics Criterion J3: This criterion is used when a breach occurs in a terrain location with sufficient scour dynamics, and its expression is J3=(|▽z|>z_slope_threshold ); where |▽z| is the local terrain slope, z _slope_threshold The minimum slope threshold is set, for example, 0.02. This criterion can effectively prevent false breaches caused by calculation errors in almost level levee crest areas.

[0195] The final breach phase transition trigger signal, Flag_breach, is generated by logically combining these three criteria: Flag_breach = (J1 AND J3) OR (J2 AND J3). In other words, when the levee is in an area with erosion-prone terrain, and has experienced sufficient energy accumulation or a severe erosion impact, the system will trigger breach evolution. This multi-criteria fusion mechanism is used to cover different breach occurrence scenarios, improving the accuracy and robustness of model predictions; in other words, it supports the model's prediction of breach evolution.

[0196] As a preferred method for determining the initial failure state of critical infrastructure, this embodiment provides a comprehensive failure probability assessment model that couples multiple physical disaster-causing factors, replacing assessment methods that rely solely on a single disaster intensity index. Accordingly, for a specific infrastructure node (such as a substation), its comprehensive failure probability P... _fail The marginal failure probability can be calculated by combining different catastrophic factors. The preferred combination method is based on the inverse probability multiplication formula of the independent event assumption: P _fail =1-(1-P _iM )*(1-P _h )*(1-P _v ); where P _fail P represents the overall failure probability of the facility. _iM P represents the failure probability caused by the comprehensive strength index IM (primarily reflecting impact and erosion). _h The maximum submerged water depth h _max (Primarily reflects the probability of failure caused by immersion and equipment submersion); P _v For the maximum flow velocity v _max (Mainly reflects the probability of failure caused by dynamic pressure and drag force)

[0197] Each marginal failure probability P _k (k=IM, h, v) can all be determined using a vulnerability function specific to this facility type. For example, they can all be in the form of a Logistic function: P _k =1 / (1+exp(-α _k *(X _k -X _k , 50 )));wherein, X _k It corresponds to the magnitude of the hazard factor (IM, h) _max or v_max ); α _k and X _k，50 This refers to the vulnerability parameter (curve steepness and threshold corresponding to a 50% failure rate) of this type of facility to this type of disaster-causing factor. This parameter can be calibrated through statistical analysis of historical disaster data or structural reliability simulation. For example, for outdoor open substations, they may be very sensitive to the flooding depth h (h... 50 The water depth is relatively low, while the sensitivity to flow velocity v is relatively low. For bridge piers, however, they are very sensitive to IM and v, while their sensitivity to h is reflected in whether the water depth exceeds the bridge deck.

[0198] This multi-factor coupled assessment model can characterize the differentiated failure mechanisms of different types of infrastructure. The model incorporates the different characteristics of the destructive effects of different disaster factors on different facilities, which are used to provide the initial failure state for cascading failure simulation, thus improving the accuracy of the entire indirect loss assessment chain.

[0199] This invention employs a combination of instantaneous scouring power P(t) and cumulative scouring energy W. _cum And the erosion energy threshold W of dike materials _cr The proposed breach triggering criterion addresses the issues of missing physical criteria for breach occurrence and reliance on manually set or simplified hydraulic thresholds. Specifically, this method calculates the scouring power P(t) of the water flow on the levee at each moment by coupling the overpass unit width flow rate q(t), the local topographic slope |▽z|, and a kinetic energy correction factor based on the Froude number Fr; the cumulative scouring energy W is obtained by integrating this power over time. _cum This energy quantifies the cumulative damage to the dike caused by the water flow over the entire overtopping period; this cumulative energy is then compared with the erosion resistance energy threshold W, which is calibrated based on physical parameters such as soil cohesion c and internal friction angle φ. _cr Real-time comparisons are performed. When the accumulated energy exceeds a threshold, the system automatically triggers a phase transition from stable weir flow to dynamic breach. This algorithm transforms the levee failure process into a calculable and predictable physical model. In flood emergency scenarios in plain river network areas, the simulation system no longer needs to rely on expert experience to guess which levee section might breach at what time. Instead, it can identify dangerous levee sections online and predict their breach timing based on real-time evolving hydraulic loads, improving the automation level of flood forecasting and the accuracy of early warnings, thus buying time for the deployment of emergency resources and the evacuation of downstream personnel. In other words, the simulation system can be used to support automated operations and early warning-related calculations in flood forecasting.

[0200] Example 9 provides an end-to-end system integration method from physical simulation to comprehensive risk assessment, demonstrating the integration of various technical modules into an automated, end-to-end comprehensive risk assessment system; it describes the data flow, module collaboration, and final output of the methods in the aforementioned examples at the system level.

[0201] In this embodiment, the entire evaluation system can be divided into two stages: an offline preparation stage and an online real-time evaluation stage. Specifically:

[0202] The offline preparation phase includes a series of preparatory tasks to build the system's decision-making framework and data foundation before a flood event. Specifically, the system acquires and processes a historical disaster damage case database. This process includes preprocessing the input features (comprehensive intensity index IM, impact duration T, and warning lead time W), such as nonlinear standardization (logarithmic transformation, exponential transformation) and data augmentation through Bootstrap resampling. Based on this, hyperparameters are optimized using methods such as K-means clustering and cross-validation to train and solidify the RBF network, which constitutes the three-dimensional vulnerability surface knowledge base.

[0203] Accordingly, the system needs to pre-train the GNN model. This requires utilizing historical infrastructure cascading failure data or simulation data generated by expert systems to allow the GNN to learn the complex nonlinear patterns of failure propagation between different types of infrastructure. The learning results are reflected in the weight parameter matrix between network layers.

[0204] Furthermore, the system will analyze socioeconomic data, extract node attributes (type, capacity, location, etc.) of critical infrastructure, and construct a weighted, multi-layered infrastructure dependency graph G based on multi-dimensional dependencies including physical, geographical, and informational factors. _infra This map serves as the foundational topology for online inference by the GNN model. Based on this, the system integrates and preprocesses all necessary geographic and engineering data, including but not limited to DEM, river network topology, levee elevation, soil parameters, land use types, and asset spatial distribution, establishing an efficient database index to support rapid invocation during the online phase.

[0205] In the online real-time assessment phase, specifically: when a real flood event occurs, the system enters online operation mode, with the data flow and calculation chain as follows: The system receives real-time hydrological monitoring data (such as time-series data from upstream water level stations) as input. The hydrodynamic model is activated, efficiently solving for the water level and flow distribution across the entire river network, monitoring the super-high margin of each levee section in real time, and dynamically identifying potential overtopping areas. For identified overtopping areas, the overtopping-breach phase transition simulation module is triggered. This module, based on the cumulative scour energy criterion, determines whether a breach has occurred. Once triggered, it accurately simulates the dynamic evolution of the breach, generating the breach-closing flow process line Q. _breach (t). The breach flow process line Q _breach(t) serves as the source term input for the two-dimensional flood evolution model. The model simulates the evolution of floods on complex terrain, calculating and generating a spatialized comprehensive intensity index field IM(x,y) and an impact duration field T(x,y).

[0206] Furthermore, the intensity field IM(x, y) and duration field T(x, y), along with the early warning lead time W obtained from the emergency response system, are passed to the three-dimensional vulnerability assessment module. This module invokes a pre-trained RBF network to quickly retrieve the spatialized loss rate field and calculates the direct economic loss L by combining it with asset data. _direct The initial state vector is then fed into the GNN cascaded failure simulation module. The GNN performs rapid inference on a pre-constructed dependency graph, outputting the cascaded failure sequence S[i,t]. This sequence is then fed into the input-output model to calculate the time evolution curve L of the indirect economic loss. _indirect (t).

[0207] Accordingly, the system integrates all assessment results to generate a multi-dimensional comprehensive disaster loss assessment report. This report not only includes the total direct and indirect economic losses, but also provides: spatial visualization maps, such as spatial distribution maps of direct economic losses and comprehensive disaster intensity level maps, intuitively showing the areas affected by the disaster; temporal evolution analysis, such as curves showing the accumulation of indirect economic losses over time, revealing the persistence and dynamic changes in the disaster's impact; categorized statistical information, with loss statistics categorized by surface type (agriculture, industry, residential areas) and infrastructure type (electricity, transportation, water supply), providing a basis for targeted rescue and recovery; and uncertainty quantification. For all key assessment results (such as total loss), the system preferably uses Monte Carlo simulation to propagate the uncertainty across the entire chain from hydrological inputs and model parameters to economic model parameters, presenting it in the form of expected values ​​and confidence intervals (e.g., 90% confidence intervals), providing decision-makers with a comprehensive understanding of the risk.

[0208] According to one aspect of this application, the time-varying evolution equation process of the breach geometry specifically involves the following: after the breach is triggered and its initial geometry is determined, its subsequent lateral, longitudinal, and slope evolution is driven by a set of coupled differential equations. This set of equations can preferably take the following form: Bottom vertical scour evolution equation: dh / dt=β _v *(τ _b_turb -τ _c ) 1.5 *exp(-γ _v *t / T _scale )*f _sat (h _b ); where dh / dt is the rate of change of the ulcer depth h with time t; β_v This is the vertical scour factor, for example, it can be taken as 5 × 10⁻⁶. -4 m*s -1 *Pa -1.5 Its value can be determined according to the type of soil in the dike; τ _b_turb The bottom shear stress is corrected for turbulence; τ _c γ is the critical shear stress of the soil; _v T is the time decay coefficient, which can be taken as 0.5 to reflect the effect of soil gradually compacting as erosion deepens; _scale For the characteristic time scale, e.g., 3600 seconds; f _sat (h _b ) is a saturation function, for example (1-h _b / h _max ) 0.5 Used to limit the maximum scour depth h _max (e.g. h) _max It can be taken as 80% of the depth of the main channel of the river.

[0209] Lateral expansion evolution equation: dB / dt=β _h *τ _b_turb 1.5 *(1-B / B _max )*g_shape(B / h); where dB / dt is the rate of change of the ulcer base width B with time t; β _h This is the lateral scour coefficient, for example, it can be taken as 3×10. -4 m*s -1 *Pa -1 *5;(1-B / B _max ) is a saturation term used to limit the maximum ulcer width B. _max (e.g., 100 meters); g_shape(B / h) is the shape function, for example, 1+2*exp(-((B / h-5) / 3) 2 ), used to promote the width-to-depth ratio of the ulcer to a stable value (such as 5 here).

[0210] Slope evolution equation: dm / dt = -0.01*(mm) _stable )*max(0, τ _b_turb / τ _c -1); where dm / dt is the rate of change of the slope ratio m with time t; m _stable Let be the stable slope ratio of the soil resting underwater, for example, 2.0; the max(0,...) term indicates that the slope will only adjust to a stable state when the actual shear stress exceeds a critical value. The driving force in the above equations is the base shear stress τ. _b_turb The calculation process is as follows: Calculate the flow area A based on the current breach geometry parameters [B(t), h(t), m(t)]. _b=(B+m*h)*h and wetted perimeter P _w =B+2h*sqrt(1+m 2 ), thus obtaining the hydraulic radius R _h =A _b / P _w Furthermore, based on the current breach flow rate Q... _b Calculate the average flow velocity v _b =Q _b / A _b Next, the foundation shear stress τ is calculated using Manning's formula. _b =ρ*g*n 2 *v _b 2 / R _h 1 / 3 Based on this, turbulence correction is applied to obtain τ. _b_turb =τ _b *(1+0.3*max(0, v) _b 2 / (g*R _h )-1)).

[0211] According to one aspect of this application, the calculation process for time-varying breach flow under different scenarios is as follows: at each time step of breach evolution, based on the updated breach geometric parameters [B(t), h(t), m(t), Z...]... _b (t)] and upstream water level H _up (t), distinguishing between submerged outflow and free outflow for calculation: when H _up >Z _b At +1.5h, it is determined to be a submerged outflow, with a breach flow rate Q. _b (t) Calculated using the orifice flow formula: Q _b (t)=C _sub *A _b *sqrt(2g*(H _up -Z _b -0.5h)); where C _sub The submerged outflow coefficient can be taken as 0.65; A _b Z represents the cross-sectional area of ​​the water passage; _b H is the elevation of the bottom of the breach; g is the acceleration due to gravity. _up ≤Z _b At +1.5h, it is determined to be free outflow, and the breach flow rate Q _b (t) Calculated using the broad-crested weir formula: Q _b (t)=(2 / 3)*C _free *B _eff *sqrt(2g)*(H _up -Z _b ) 1.5 Among them, C_free The free outflow coefficient can be taken as 0.45; B _eff For effective overcurrent width, B _eff =B+2m*h / 3.

[0212] According to one aspect of this application, the mass conservation verification mechanism for breach evolution specifically involves: at the end of each calculation time step, performing a water balance verification on the control volume containing the breached river segment. The verification formula is: |Q _in -Q _out -dV / dt|<ε_mass*Q _in Where ε_mass is a very small tolerance, for example, 0.01. If this inequality does not hold, it indicates that mass is not conserved. In this case, the breach flow rate needs to be checked, and the checked breach flow rate Q is generated. _b_verified The preferred verification method is to adjust proportionally: Q _b_verified =Q _b *(Q _in -(Q _out -Q _b )-dV / dt) / Q _b Q _b_verified The final breach flow at that time step is fed back to the river network model.

[0213] According to another aspect of this application, the stability judgment and termination condition process for breach evolution specifically involves: during the breach evolution simulation, the system monitors the expansion rate of its geometric shape in real time, i.e., dB / dt and dh / dt. When the values ​​of dB / dt and dh / dt are both less than a preset minimum threshold ε for a continuous period of time (e.g., exceeding 30 minutes),... _rate (For example, 10) -4 When the flow rate reaches m / s, the system determines that the breach morphology has reached a quasi-steady state. At this point, the breach evolution calculation will be automatically terminated, maintaining the final breach morphology and flow rate.

[0214] According to one aspect of this application, the complete process for determining the initial breach geometry parameters specifically involves: at the initial depth h _b0 and bottom width B _b0 Based on the calculations, the initial slope ratio m is further determined. _b0 and the initial breach bottom elevation Z _b0 This forms the initial morphological parameters. The initial slope ratio m _b0 Based on Coulomb's earth pressure theory, it can be related to the internal friction angle φ of the soil, and its calculation formula is: m _b0 =1 / tan(45°+φ / 2). Initial breach bottom elevation Z _b0 : From the top elevation Z of the levee _levee and initial scour depth h _b0 Z can be obtained directly from calculation._b0 =Z _levee -h _b0 .

[0215] According to another aspect of this application, the process of monitoring ultra-high margin and predicting overtopping time specifically involves: in hydrodynamic analysis, the ultra-high margin ΔH=Z is... _levee The monitoring of -H is deepened. Accordingly, the real-time rate of change of the super-high margin is calculated: dΔH / dt=[ΔH(t)-ΔH(t-Δt)] / Δt; where Δt is the calculation time step. The system can set a warning threshold for the rate of change; for example, when dΔH / dt is less than -0.1 m / h, it indicates a rapid rise in water level, and the warning level can be automatically upgraded by one level. Based on this, the future flooding moment can be predicted by linear extrapolation: t _predict =t _now -ΔH / (dΔH / dt); where t _now This is the current time. When the predicted flood peak time t... _predict When the emergency response preparation time is less than 2 hours, the system can issue the highest level of emergency warning.

[0216] According to another aspect of this application, the concept of terrain gradient directional coupling coefficient is specifically: the directional coupling coefficient cos(θ) when calculating the comprehensive intensity index IM. _ij The calculation method is: cos(θ) _ij )=(v _ij •▽z _ij ) / (|v _ij |•|▽z _ij |); where v _ij The direction vector of water flow [vx] _ij vy _ij ];▽z _ij The terrain slope direction vector is [Яz / Яx, Яz / Яy]. _ij • indicates the vector dot product.

[0217] According to another aspect of this application, the classification and weighting process for infrastructure dependency networks specifically involves: when constructing an infrastructure dependency graph, the dependency edges between nodes can be subdivided into various types and assigned different weights to reflect the differences in dependency relationships. Preferred classification methods include: physical dependency, which refers to the fact that the normal operation of a facility must depend on the physical output of another facility, such as a water plant depending on the power supply from a substation. This type of dependency is the strongest, with a weight w. _phy This can be set to a higher value, such as 0.9. Information dependency refers to a facility needing to obtain information or data from another facility to function properly, such as a traffic management system relying on a communication network. This type of dependency is of secondary strength, with a weight w. _infThis can be set to a moderate value, such as 0.6. Geographical dependence refers to the potential for facilities in similar geographical locations to influence each other due to shared access, secondary disasters, etc. The strength of this type of dependence is inversely proportional to distance, and its weight w _geo It can be set to C / d _ij Where C is a constant (e.g., 0.3), d _ij This refers to the distance between facilities.

[0218] According to another aspect of this application, the concept of GNN propagation rules specifically refers to: updating S in the GNN state. _i (t+1)=min(1,S) _i (t)+Σ _j w _ij *T _ij *S _j In (t)), the propagation coefficient T _ij Used to adjust the efficiency of failure state propagation from node j to node i. T _ij It can be a function related to multiple factors, and a preferred implementation is: T _ij =β*exp(-d _ij / d_0)*type_match(i,j); where β is the basic propagation rate, such as 0.4; exp(-d _ij The term / d_0) represents the propagation effect as a function of distance d. _ij Attenuation, d_0 is the feature distance; type_match(i,j) is the type matching coefficient. Failure propagation may be faster between facilities of the same type (coefficient>1), while it is slower between facilities of different types (coefficient<1).

[0219] According to another aspect of this application, the specific optimization process of the RBF network parameters is as follows: during the training of the RBF network V(x)=Σ _i λ _i *φ(||xc _i || / r _i When ), its parameter c _i (center) and r _i Determining the radial direction is crucial for model performance. One method is to apply K-means clustering to the input features x of all training samples, and use the centroids of the resulting M clusters as the M centers c of the RBF network. _i Accordingly, for each center c _i Its radial parameter r _i The distance from the center to its k-th nearest neighbor center can be calculated, where the value of k can be optimized using methods such as leave-one-out cross-validation to minimize the prediction error.

[0220] This invention introduces a physical criterion based on accumulated scour energy to achieve automatic phase transition simulation from stable overtopping to dynamic breach; by constructing a comprehensive intensity index that integrates water depth, flow velocity, slope, and duration, it improves the scientific nature of disaster intensity measurement; and by constructing a three-dimensional vulnerability surface that includes the early warning time dimension and applying graph neural networks to simulate cascade failure, it enhances the realism and accuracy of disaster damage assessment.

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

Claims

1. A method for simulating and evaluating flood damage of water network dike overflow, characterized in that, The method comprises the following steps: performing water dynamics analysis based on hydrological geographical data to determine potential overtopping areas; for the potential overtopping areas, performing overtopping flow process simulation and simulating breach evolution based on cumulative scour energy criterion to generate a breach flow process; simulating flood evolution according to the breach flow process to generate a spatialized disaster intensity field; combining the disaster intensity field and social and economic data to assess direct and indirect economic losses and generate a comprehensive disaster loss result; wherein simulating the breach evolution based on the cumulative scour energy criterion comprises: determining the instantaneous scour power of overtopping flow for each potential overtopping area; time-integrating the instantaneous scour power to obtain the cumulative scour energy; triggering the breach evolution when the cumulative scour energy exceeds the anti-scour energy threshold value set based on the physical properties of the embankment material; simulating the time-varying process of the geometric morphology after the breach is triggered, and coupling it with the in-channel storage state quantity to calculate and generate the breach flow process; wherein generating the spatialized disaster intensity field comprises: calculating the time and space distribution of the submerged water depth and flow velocity field according to the breach flow process and the hydrological geographical data; determining the comprehensive intensity index and impact duration for any calculation unit in the submerged area based on the time and space distribution of the submerged water depth and flow velocity field; constructing the spatialized comprehensive intensity index field and the impact duration field from the comprehensive intensity index and impact duration of each calculation unit, respectively, to form the spatialized disaster intensity field; wherein the submerged area is a spatial range determined based on the submerged water depth being greater than a set water depth threshold value; wherein determining the comprehensive intensity index for any calculation unit in the submerged area comprises: analyzing the water depth time sequence and flow velocity time sequence of the calculation unit based on the time and space distribution of the submerged water depth and flow velocity field; extracting the local terrain slope of the calculation unit according to the hydrological geographical data; coupling the water depth, flow velocity and local terrain slope of the calculation unit for each submerged time to obtain the instantaneous impact power density; integrating the instantaneous impact power density along the time dimension to obtain the comprehensive intensity index of the calculation unit.

2. The method of claim 1, wherein, determining the instantaneous scour power of overtopping flow comprises: analyzing the hydraulic conditions of the potential overtopping area to obtain the unit width flow of overtopping flow; extracting the local terrain slope of the potential overtopping area from the hydrological geographical data; calculating the Froude number according to the unit width flow of overtopping flow to generate a kinetic energy correction factor accordingly; determining the instantaneous scour power by integrating the unit width flow of overtopping flow, the local terrain slope and the kinetic energy correction factor.

3. The method of claim 1, wherein, setting the anti-scour energy threshold value comprises: extracting the physical properties of the embankment material of the potential overtopping area from the embankment monitoring data and the INSAR-based embankment deformation radar remote sensing results, the physical properties including the soil cohesion and internal friction angle; calculating the critical shear stress of the embankment material according to the soil cohesion and internal friction angle; converting the critical shear stress into the anti-scour energy threshold value.

4. The method of claim 1, wherein, after triggering the breach evolution, further comprising: calculating the amount of the cumulative scour energy exceeding the anti-scour energy threshold value to obtain the super-threshold scour energy; defining the initial breach geometric parameters based on the super-threshold scour energy; simulating the time-varying process of the geometric morphology after the breach is triggered based on the initial breach geometric parameters.

5. The method of claim 1, wherein, Before coupling the water depth, the flow velocity and the local terrain slope to determine the instantaneous impact power density, further comprising: identifying the flow evolution direction according to the flow velocity time series and the slope direction according to the local terrain slope; based on the vector relationship between the flow evolution direction and the slope direction, determining the directional coupling coefficient; modulating the local terrain slope by using the directional coupling coefficient to generate the effective terrain slope; coupling operation, specifically coupling the water depth, the flow velocity and the effective terrain slope.

6. The method of claim 1, wherein, evaluating the direct economic loss, including: for any calculation unit in the flooded area, extracting its comprehensive intensity index and impact duration from the spatialized disaster intensity field; obtaining the warning lead time contained in the socio-economic data; applying the pre-constructed three-dimensional vulnerability surface, according to the comprehensive intensity index, the impact duration and the warning lead time, querying and determining the loss rate of the calculation unit; combining the asset spatial distribution data contained in the socio-economic data and the loss rate, calculating the direct economic loss of the calculation unit.

7. The method of claim 6, wherein, After determining the loss rate of the calculation unit, further comprising: evaluating the prediction uncertainty of the three-dimensional vulnerability surface at the query point to generate the confidence interval of the loss rate; wherein evaluating the prediction uncertainty specifically includes: calculating the Jacobian matrix of the three-dimensional vulnerability surface at the query point with respect to its input dimensions; combining the Jacobian matrix and the uncertainty covariance of the input data to calculate the prediction variance of the loss rate, and accordingly generating the confidence interval.

Images