A high-resolution flood simulation method coupling topological optimization of motion wave and cellular automaton

Abstract

The application discloses a high-resolution flood simulation method coupling topological optimization motion wave and cellular automata, comprising the following steps: obtaining the basin data of a target region, and re-projecting the basin data to a high-resolution grid coordinate system to construct an upstream and downstream relationship grid graph; traversing the upstream and downstream relationship grid graph, taking the rainfall runoff obtained at a current time as a lateral inflow, combining a water flow state at a previous time step and an inflow state of a corresponding grid in an upstream to construct a one-dimensional motion wave equation, and solving the one-dimensional motion wave equation to obtain a river water level at the current time; taking the river water level as an initial condition, starting a two-dimensional diffusion process on the basis of the upstream and downstream relationship grid graph, and iteratively obtaining a submerged range and a water depth distribution graph at the current time step; and repeating the above steps to obtain the submerged range and the water depth distribution graph at all time steps. The method provided by the application can realize fast deduction of a large-scale, high-resolution and physically strictly conservative flood evolution process.

Description

Technical Field

[0001] This invention belongs to the field of smart water conservancy and flood control safety technology, and in particular relates to a high-resolution flood simulation method that couples topology optimization moving waves and cellular automata. Background Technology

[0002] Existing technologies have the following limitations when performing flood simulation and inundation early warning for large-scale, high-resolution complex terrain: The contradiction between the computational efficiency of traditional physical models and the detailed representation of flooding: Fully two-dimensional hydrodynamic models used for fine simulation (such as models that solve the complete Saint-Venant equations) are theoretically accurate and can describe the extent of surface inundation and water depth distribution.

[0003] However, such models require solving massive mathematical matrices, resulting in extremely high computational costs and time consumption when dealing with large-scale, high-resolution grid scenarios at the city or watershed level. This makes it difficult to meet the real-time simulation and early warning requirements for the sudden and intense nature of extreme heavy rainfall. While some simplified one-dimensional routing models are faster, they often fail to describe the specific inundation range, water depth distribution, and spatial dynamics of flood evolution.

[0004] Limitations of numerical models in water conservation: Existing distributed physical models or simplified hydrodynamic models, when simulating complex flood evolution processes, are often limited by numerical calculation formats, making it difficult to strictly maintain the water balance within the system. Some simplified numerical schemes may lead to non-physical increases or decreases in the total water volume when facing atypical terrain or extreme flow velocity scenarios. Furthermore, simulation methods lacking strict physical mass conservation constraints are prone to producing anomalies that violate the principles of mass conservation or basic hydraulic principles when dealing with extreme climate events beyond the scope of observational data, thus reducing the reliability of early warning results.

[0005] The lack of micro-connectivity in topographic representation: Flood simulations heavily rely on digital elevation models (DEMs) to characterize micro-topographic connectivity. However, the accuracy of existing publicly available DEM data often falls short of the requirements for high-resolution simulations, especially in narrow valleys or complex river networks. This can easily lead to the formation of "pseudo-depressions" or broken "pseudo-flow paths," preventing upstream water from accurately converging or discharging, severely impacting the accuracy of confluence calculations and inundation simulations. Furthermore, the massive data volume generated by large-scale high-resolution simulations is constrained by traditional storage and transmission methods, hindering the rapid release of simulation results.

[0006] Patent document CN121684661A discloses a method, device, and equipment for calculating the drainage and flood control capacity of urban low-lying roads. The method includes: acquiring and quantifying road characteristics and constructing a road characteristic parameter library; then associating rainfall and road water accumulation data with the parameter library to form a road water accumulation scenario parameter library; subsequently selecting characteristic parameters from the scenario library and training a prediction model with the submerged area or water depth as the target; and finally using the actual characteristic parameters of the target road to predict its water accumulation index through the model.

[0007] Patent document CN121436692A discloses a method for real-time flood control scheduling risk dynamic early warning and emergency plan formulation in a watershed, including: acquiring real-time full-element information tensor and engineering flood control design parameter set; calculating the real-time hierarchical risk entropy time series based on safety margin distribution; inputting real-time data and pre-constructed entropy flux Markov diagram structure data into a trained entropy-aware dual-flow prediction model for inference to obtain a future risk prediction set. Summary of the Invention

[0008] The purpose of this invention is to provide a high-resolution flood simulation method that couples topology-optimized moving waves with cellular automata. This method simplifies the problem that hydrodynamic models cannot strictly maintain water balance, thereby enabling rapid deduction of large-scale, high-resolution, and physically strictly conserved flood evolution processes.

[0009] To achieve the objectives of this invention, the following technical solution is provided: a high-resolution flood simulation method coupling topology-optimized moving waves and cellular automata, comprising: Step 1: Obtain watershed data for the target area and reproject the watershed data onto a high-resolution grid coordinate system to construct a grid diagram of upstream and downstream relationships based on flow direction. The watershed data includes a high-resolution digital elevation model, flow direction data, and catchment area data. Step 2: Traverse the upstream and downstream relationship grid at a preset time step, take the rainfall runoff collected at the current time as the lateral inflow, and construct a one-dimensional kinematic wave equation by combining the water flow state of the previous time step and the inflow state of the corresponding upstream grid, and solve the one-dimensional kinematic wave equation to obtain the river level at the current time step. Using the calculated river water level as the initial condition, a two-dimensional diffusion process is initiated based on the upstream and downstream relationship grid diagram to iteratively obtain the inundation range and water depth distribution map at the current time step. Step 3: Repeat Step 2 to obtain the inundation range and water depth distribution maps for all time steps.

[0010] This invention constructs a river network topology based on high-resolution topographic data, uses the runoff area to topologically sort the computational grid to determine an efficient computational time series, and establishes an empirical mechanism for estimating river width parameters. Secondly, a one-dimensional kinematic wave confluence module based on the topological sequence is constructed to receive meteorological runoff input. An implicit finite difference algorithm combined with a Newton-Raphson iterative algorithm is used to quickly solve for the flow evolution and initial water level within the river network channel. Subsequently, a conserved two-dimensional cellular automaton diffusion module is established, driven by the one-dimensional calculated water level. Damping coefficients and half-height difference flow limiting thresholds are introduced to calculate the diffusion flux between grids, and a synchronous state update mechanism ensures strict mass conservation during surface flooding. Finally, data quantization and efficient compression strategies are employed to optimize the storage of large-scale, high-resolution, long-term simulation results. This method effectively balances simulation accuracy, computational efficiency, and numerical stability, significantly improving the ability to perform refined flood extrapolation in large-scale, complex river network areas.

[0011] Specifically, the construction process of the upstream and downstream relationship mesh diagram is as follows: The upstream and downstream connections between grids are established based on the raster flow direction data, and all effective grid cells are sorted in ascending topology according to the runoff area to construct the corresponding grid diagram. The river width of each grid cell is calculated using an empirical formula based on the catchment area, and the river slope and Manning roughness coefficient field of the entire region are initialized.

[0012] Specifically, the river width is calculated based on the catchment area corresponding to the grid cell using an empirical power-law formula.

[0013] Specifically, the ascending topological sorting process is as follows: The runoff area data in two-dimensional raster form is flattened and converted into a one-dimensional array; An indirect sorting algorithm is applied to the one-dimensional array to arrange the flow area values ​​in ascending order, so as to obtain an integer index array for filling the grid cells.

[0014] Specifically, implicit finite difference and Newton-Raphson iteration are used to solve the one-dimensional motion wave equation.

[0015] Specifically, the solution process is as follows: For each grid cell in the topology sequence, the actual flow length parameter is represented by a flow direction code; A nonlinear equation for a one-dimensional moving wave is constructed based on the mass conservation equation and the Manning momentum equation. The nonlinear equation is solved iteratively using the Newton-Raphson method. The iteration is repeated until the difference between two adjacent calculation results is less than the preset convergence threshold, so as to obtain the grid outflow rate at the current moment. Finally, the calculated grid outflow is added to the inflow of the downstream receiving unit of the grid at the current moment to complete the water transfer calculation for a single grid.

[0016] Specifically, the two-dimensional diffuse diffusion process is as follows: Determine whether the current water depth of each target grid cell within the target area is greater than zero; Traverse the N neighboring grid cells of the target grid cell. If the water level of the neighboring grid cell is lower than that of the target grid cell, then determine that the corresponding neighboring grid cell is a potential water-receiving cell. Calculate the water level difference between the target grid cell and the potential water-receiving cell, and sum all water level differences to obtain the total water level difference. If the total water level difference is greater than zero, then the flow from the target grid to the [missing grid name] grid is calculated based on the principle of proportional distribution of elevation differences. Diffusion flux of a potential water-receiving unit; Subtract all calculated diffusion fluxes from the current water depth of the target grid cell and increase the water depth of each potential receiving cell accordingly.

[0017] Specifically, the calculation process for the lateral inflow is as follows: At the start of each simulation time step, acquire gridded surface runoff depth rate data covering the entire region at the current moment; For each specific grid cell, the total runoff volume generated within the time step is calculated by combining the corresponding grid resolution size and the current time step. The total runoff volume is converted into the lateral inflow per unit river length in the one-dimensional kinematic wave equation.

[0018] Compared with the prior art, the beneficial effects of the present invention are as follows: By introducing a topological sorting mechanism based on runoff area, the shortcomings of low computational efficiency in traditional full two-dimensional models are solved. Furthermore, by constructing a conservation-type state synchronization update mechanism, the problem of the simplified hydrodynamic model being unable to strictly maintain water balance is solved. Ultimately, a rapid simulation of the large-scale, high-resolution, and physically strictly conserved flood evolution process is achieved. Attached Figure Description

[0019] Figure 1 This is a flowchart of the high-resolution flood simulation method using coupled topology optimization motion waves and cellular automata provided in this embodiment; Figure 2 This is a schematic diagram illustrating the principle of grid topology sorting and one-dimensional motion wave confluence provided in this embodiment; Figure 3 This is a schematic diagram of the conservation two-dimensional cellular automaton diffusion and water exchange mechanism provided in this embodiment; Figure 4This is a result diagram of a high-resolution (90-meter) flood simulation example of the target watershed provided in this embodiment; Figure 5 This figure shows the results of long-term daily runoff simulation verification and accuracy evaluation at a typical hydrological station, as provided in this embodiment. Detailed Implementation

[0020] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. The components of the embodiments of the present invention described and shown in the accompanying drawings can generally be arranged and designed in various different configurations. Therefore, the following detailed description of the embodiments of the present invention provided in the accompanying drawings is not intended to limit the scope of the claimed invention, but merely to illustrate selected embodiments of the invention. All other embodiments obtained by those skilled in the art based on the embodiments of the present invention without inventive effort are within the scope of protection of the present invention.

[0021] like Figure 1 The figure shows the high-resolution flood simulation method provided in this embodiment, and its specific steps are as follows: Step 1: First, obtain the high-resolution digital elevation model (DEM), flow direction data (FDR), and catchment area data (UPA) of the target simulation area (large-scale urban watershed).

[0022] All data are uniformly reprojected and aligned to the same high-resolution grid coordinate system. Upstream and downstream connections between grids are established based on raster flow direction data, and all effective grid cells are sorted in ascending topology according to the catchment area (UPA) to construct an upstream-downstream relationship grid diagram.

[0023] Meanwhile, the river width of each grid cell is calculated using an empirical formula based on the catchment area, and the river slope and Manning roughness coefficient field of the entire region are initialized.

[0024] The river width for each grid cell is calculated using an empirical power-law formula based on the catchment area (UPA) of that grid cell. The specific calculation formula is as follows:

[0025] in, For grid cells Estimated river width (unit: meters). This is the cumulative catchment area value obtained by the grid cell during the preprocessing stage (the original unit is square kilometers or grid number, which is converted to square meters here).

[0026] Based on the catchment area (UPA), all valid grid cells are sorted in ascending topology to construct an upstream-downstream relationship grid diagram, as follows: The two-dimensional raster-style runoff area data (UPA) is flattened and converted into a one-dimensional array. Then, an indirect sorting algorithm is applied to this one-dimensional array to arrange the runoff area values ​​in ascending order.

[0027] The output of the sorting operation is an integer index array that stores the original position index of the grid cell before flattening. These indices are strictly arranged in order from the smallest confluence area (corresponding to the upstream source area) to the largest confluence area (corresponding to the downstream outlet area), thus ensuring that the calculation process strictly follows the physical causal relationship of water flow from upstream to downstream. This avoids complex iterative solutions or large matrix operations, and significantly improves computational efficiency.

[0028] Furthermore, a minimum width threshold constraint mechanism is applied to the calculated river width: if a certain grid cell calculates the river width according to the above formula... If the width is less than the preset minimum physical width threshold (0.5 meters), the river width of the grid cell will be forcibly assigned to this minimum value to ensure the physical rationality and numerical stability of subsequent one-dimensional moving wave hydraulic calculations in extremely small watersheds or source regions.

[0029] Step 2: Within each simulation time step, based on the upstream and downstream relationship mesh diagram obtained in the above steps, traverse all mesh cells strictly in order from upstream to downstream. The specific process is as follows:

[0030] For the currently traversed grid, rainfall runoff from meteorological data is received as lateral inflow. Combined with the flow state of the previous time step and the inflow transmitted from the upstream grid, the one-dimensional motion wave equation is solved using the finite difference method and the Newton-Raphson iterative algorithm to calculate the outflow, channel storage, and channel water level of the current grid.

[0031] Among them, rainfall runoff received from meteorological data is used as lateral inflow. The specific calculation and processing process is as follows: at the beginning of each simulation time step, gridded surface runoff depth rate data covering the entire area at the current moment is acquired. (The original unit is millimeters per hour). First, convert this runoff depth rate to velocity form in standard SI units (meters per second): ; Subsequently, for each specific grid cell, combined with its grid resolution size ( (unit: meters) and current time step ( (Unit: seconds), calculate the total runoff volume generated within this time step ( (Unit: cubic meters) ; Finally, this total runoff volume is converted into a lateral inflow rate per unit river length suitable for the one-dimensional kinematic wave equation. (Unit: square meters per second or cubic meters per second / meter) ; The result That is, the lateral inflow term at the current moment ( Substitute the equation into the one-dimensional motion wave finite difference equation described in step (2) to solve it.

[0032] The river water level is derived from the calculated river storage capacity based on the assumption of a wide and shallow rectangular river channel. The specific steps are as follows: The river channel storage capacity obtained by solving the equation under the one-dimensional moving wave approximation ( ) is considered equivalent to the cross-sectional area of ​​the water passage ( ).

[0033] Assuming the cross-sectional shape of the river channel within each grid cell is a wide and shallow rectangle, its bottom width is the river width of that grid cell calculated in step (1). ).

[0034] Based on the formula for the area of ​​a rectangular cross-section, the average water depth within this grid cell is calculated using the following relationship ( (Unit: meters) ; During the calculation process, if the grid river width If the value is zero, then the water depth is set to zero.

[0035] Finally, the calculated average water depth is added to the ground elevation (DEM) of the grid cell to obtain the river level of the grid, which is used as the initial driving head condition for the two-dimensional diffusion process in subsequent steps.

[0036] Simultaneously, the calculated outflow is directly added to the inflow of the downstream grid at the next time step, completing the unidirectional transfer of water volume. The specific process is as follows:

[0037] Using the calculated river level as the initial condition, a two-dimensional diffusion process is initiated based on the upstream and downstream relationship grid diagram to iteratively obtain the inundation range and water depth distribution map at the current time step.

[0038] For each grid cell containing water depth or received water volume, search for valid grids within its 8-neighborhood and calculate the water level difference between the central grid and its neighboring grids.

[0039] The diffusion flux is calculated based on the water level difference, and a damping coefficient and a half-height difference flow limit threshold are introduced to prevent numerical oscillations.

[0040] A strict water conservation exchange operation is performed, that is, while deducting the outflow of water in the central grid, an equal amount of water is added in the corresponding neighboring grid. Through multiple iterations, the two-dimensional spread and depression filling process of water flow on the surface is realized, and finally the flood range and water depth distribution map of the current time step are generated.

[0041] Step 3: Within the set time window, repeat the process in Step 2 until the simulation of all preset time steps is completed.

[0042] Upon reaching a specific output node, the high-resolution flow field of the entire watershed will be displayed. The two-dimensional flood depth field (m) and the two-dimensional inundation depth field are stored and output in geospatial raster (such as GeoTIFF) format, ultimately forming a complete, detailed and physically reliable dataset of the flood evolution process.

[0043] like Figure 2 The diagram shows the specific implementation logic of grid topology sorting and one-dimensional motion wave confluence calculation based on the catchment area (UPA) provided in this embodiment.

[0044] Traditional two-dimensional hydrodynamic models typically require the construction and solution of a massive global sparse matrix when solving the complete Saint-Venant equations. This leads to an exponential increase in computation time when dealing with large-scale, high-resolution complex watersheds, making it difficult to meet the timeliness requirements of flood warnings.

[0045] like Figure 2 As shown in (a), this represents the initial state of the physical terrain grid, including weights for flow direction (FDR) and catchment area (UPA). During model preprocessing, the unidirectional flow direction of each grid cell is first extracted based on the high-resolution DEM, and the upstream catchment area of ​​each grid cell is calculated accordingly. In the original two-dimensional physical space, due to the crisscrossing river network, the convergence of water flows and the intricate upstream-downstream dependencies, directly performing hydrodynamic iterations on the two-dimensional grid easily leads to frequent boundary condition exchanges and non-convergence of nonlinear equations.

[0046] like Figure 2 As shown in (b), this is a diagram of the dimensionality reduction calculation sequence after UPA ascending sort and the one-dimensional motion wave confluence principle. In this embodiment, the two-dimensional grid cell data is flattened and strictly sorted in ascending order from the minimum to the maximum value according to the confluence area contained therein, thereby constructing a brand-new upstream and downstream relationship grid diagram.

[0047] This UPA-based ascending index essentially reconstructs a physically directed acyclic graph in computer memory that absolutely follows the causal relationship of "water flowing downhill".

[0048] The diagram clearly distinguishes between the computer's underlying execution logic and the logic for transmitting physical water volume.

[0049] Within each specific simulation time step, the computation program completely abandons complex spatial two-dimensional neighborhood searches or large matrix inversion operations. The program only needs to perform a single traversal along the dashed path, according to the topologically sorted one-dimensional array indices, strictly from the upstream to the downstream.

[0050] For the currently traversed grid, the system invokes the implicit finite difference and Newton-Raphson iterative algorithm to quickly solve the one-dimensional nonlinear equation of the moving wave, obtaining the current outflow and river level of the grid. The specific execution logic employs a numerical solution method combining implicit finite difference and Newton-Raphson iteration. For each grid cell in the topology sequence, the actual flow length is first determined based on its D8 flow direction code. If the flow direction is due east, due south, due west, or due north, then Equal to grid resolution If the flow direction is northeast, southeast, southwest, or northwest, then .

[0051] Subsequently, a nonlinear implicit difference scheme for a one-dimensional moving wave is constructed based on the mass conservation equation and the Manning momentum equation. The grid outflow rate at the current moment is defined as... Construct the following residual function : ; in, For time step; The experience index is 0.6. To comprehensively reflect the cross-sectional shape and roughness of the river channel, the approximate relationship for its calculation is as follows: = ,in This is the Manning coefficient. For wet period, The slope of the river channel. The set of known items, which includes the known state from the previous time step and the inflow from the current time step, is calculated using the following formula: ; in, For the grid outflow of the previous time step, This represents the cumulative inflow from upstream at the current moment. and These represent the lateral inflow rates at the previous and current time points, respectively. The Newton-Raphson method is used to iteratively solve the nonlinear equations. In the first In each iteration, the flow estimate is updated using the residual function and its derivative with respect to the flow rate. : ; Repeat the iteration until the difference between two consecutive calculation results is less than a preset convergence threshold. After the iteration converges, the final result is obtained. As the grid outflow at the current moment The current channel storage (cross-sectional area) is updated based on the water level-discharge relationship under the approximation of kinematic waves. : ; Finally, the calculated results The inflow rate is added to the current inflow rate of the downstream receiving unit of the grid to complete the water transfer calculation for a single grid.

[0052] At the same time, once the computer has traversed along the dotted line and completed the calculation of the outflow of a certain grid, it strictly follows the... Figure 2 The predetermined flow direction code described in (a) will directly and accurately deliver and accumulate the calculated water volume to the absolute downstream grid corresponding to it in the original two-dimensional physical space, as the known inflow boundary condition of the downstream grid.

[0053] Through the decoupled design of the one-dimensional dimensionality reduction traversal calculation with the two-dimensional physical precision delivery with the orange line, this embodiment ensures that when the program traverses and calculates any mid-to-downstream grid, the runoff of all upstream grids must have been calculated and delivered in advance.

[0054] The above mechanism not only perfectly matches the physical causal law of rainfall confluence propulsion, but also enables extremely high computational efficiency and absolute numerical stability when dealing with tens of millions of high-resolution grids.

[0055] like Figure 3 The diagram shows a specific implementation of the conservation-type two-dimensional cellular automaton diffusion and water exchange mechanism for the surface flooding process in this embodiment.

[0056] When the water level in the one-dimensional channel exceeds the riverbank elevation, or when the surface grid receives over-infiltration / saturation runoff, the following steps are performed: For each target grid cell within the simulation region First, determine its current water depth. Is it greater than zero? If so, calculate its current water level elevation. ,in This represents the ground elevation of the grid.

[0057] Then, the 8-neighborhood grids of the target grid cell are traversed. For each valid neighboring grid... Calculate its water level elevation If the water level in the neighboring area is lower than the target grid water level (i.e. If the neighboring area is identified as a potential water-receiving unit, the water level difference between the two units is calculated. = The total water level difference of all potential water-receiving units is calculated. .

[0058] If there is a potential water-receiving unit ( Then, based on the principle of proportional allocation of elevation differences, the flow from the target grid to the first grid is calculated. Diffusion flux of each water-receiving neighborhood grid : ; in, A preset damping coefficient (ranging from 0.2 to 0.3) is used to control the diffusion rate to maintain numerical stability. Simultaneously, a half-height difference flow-limiting constraint condition is introduced to force... This is to prevent non-physical oscillations in the water flow between adjacent grids.

[0059] Finally, a water depth update operation with strict water conservation is performed: all calculated outflow fluxes are immediately deducted from the current water depth of the target grid cell, and the water depth of each receiving neighboring grid is increased accordingly, ensuring that the total water volume of the entire computational domain remains absolutely balanced during the water exchange process in a single computation step.

[0060] like Figure 3 (a) and Figure 3 As shown in (b) in the figure, it is a schematic diagram of the local hydraulic gradient search and diffusion judgment process based on 8 neighborhoods.

[0061] Based on the target topographic map, local high-resolution grid subdivision is performed to construct sub-catchments and refined grids. Within each two-dimensional evolution time step, the model uses any grid with water depth as the central cell and performs a traversal search within its Moore Neighborhood (i.e., the surrounding 8 neighboring grids). The system first extracts the absolute water level elevation of the central cell and each neighboring cell, which is the sum of the topographic elevation and the surface water depth.

[0062] Subsequently, strictly based on the physical principle that water flows downhill, the model only selects adjacent grids with an absolute water level lower than the central cell as potential water-receiving grids. This local hydraulic gradient search mechanism establishes the basic addressing logic for the natural diffusion of water flow into the surrounding depressions.

[0063] After determining the direction of spread, such as Figure 3 As shown in (d), this is a schematic diagram of the synchronous update mechanism of half-height difference flow restriction and anti-oscillation constraint and strict water conservation state based on physical profile. When calculating the diffusion flux transferred from the central grid to the receiving grid, if the calculation is directly based on empirical formulas or Manning's formula, the calculated transfer volume is very likely to be too large when the single-step time step is long or the terrain difference is extremely large. This will cause the water level of the receiving grid to exceed the water level of the central grid in the next moment, thus causing the "non-physical oscillation" phenomenon of water flow turbulence between the two grids in subsequent calculations.

[0064] like Figure 3 As shown in (c), in this embodiment, by setting a "flush threshold line", it is mandatory that within a single time step, the amount of water depth change exchanged between two grids must not exceed half of the initial water level elevation difference between the two (i.e., half-elevation difference flow restriction).

[0065] This constraint, as the highest priority cutoff condition, ensures a smooth and gradual diffusion of water flow.

[0066] In addition, such as Figure 3 The data flow architecture is shown in (b) above. To prevent data read / write conflicts during sequential grid traversal (i.e., the result of the previous grid calculation contaminates the initial state of the current grid), a state replica temporary storage mechanism is also provided. Each water exchange calculation is performed on the static read-only matrix at the current moment, while the calculated water deduction at the outflow end and the accumulated water at the inflow end are synchronously updated and written to the temporary replica matrix at the next moment. This "in-site read, off-site synchronous addition / subtraction" design achieves the strictest mass conservation at the computer level, ensuring the absolute balance of the system's global water volume, specifically including: Before each diffusion time step begins, a complete copy of the current water depth distribution data is first constructed and used as a temporary data matrix to record the state at the next time step.

[0067] During the process of traversing all grid cells to calculate the diffusion flux, the diffusion flux value from the target grid to a certain neighboring grid is calculated based on the water depth data at the current moment. Subsequently, instead of directly modifying the water depth data at the current moment, a deduction and accumulation operation is performed synchronously on the temporary data matrix: that is, in the temporary data matrix, the diffusion flux is deducted from the water depth value of the target grid cell, and the same diffusion flux is simultaneously added to the water depth value of the corresponding water-receiving neighboring grid cell.

[0068] This synchronous increase / decrease operation ensures that the spatial movement of any water volume within the entire computational domain does not lead to an increase or decrease in the total water volume of the system, thus satisfying strict mass conservation constraints. After completing the traversal and flux exchange calculation of all grid cells in the current time step, the updated temporary data matrix is ​​used to replace the current water depth data as the basis for the next round of diffusion iteration calculation, thereby completing a strictly conserved water depth update process.

[0069] like Figure 4 As shown, this embodiment demonstrates the comprehensive implementation effect in a large-scale complex watershed. The results directly verify the superior computational performance and physical reliability of this invention while balancing large-scale spatial dimensions and high-resolution details. For example... Figure 4 (a) shows the spatial distribution of one-dimensional river confluence in a typical local area. Thanks to the aforementioned grid topology sorting and one-dimensional moving wave confluence calculation scheme, the model accurately and rapidly depicts the evolution of water flow from tiny capillary tributaries to the main stream, and the logarithmic span of the flow is continuously and reasonably represented spatially. Meanwhile, as... Figure 4 (b) shows the spatial distribution of two-dimensional diffuse water depth on the surface in the corresponding local typical area.

[0070] The aforementioned conservation cellular automata can be used to reconstruct the extent and depth of localized flooding caused by overflowing river water or torrential rain on the ground surface with high fidelity.

[0071] Meanwhile, by coupling the efficient drainage of one-dimensional river networks with the fine spread of two-dimensional land surfaces in both spatial scale and physical process, the technical bottleneck of limited computing power in traditional full two-dimensional models in ultra-high resolution large-basin flood simulation is broken, providing a technical solution with extremely high engineering application value for urban flood control emergency response, urban waterlogging risk assessment and water conservancy infrastructure planning.

[0072] In addition, such as Figure 5 The figure shown is a graph illustrating the results of long-term daily runoff simulation verification and accuracy evaluation at a typical hydrological station.

[0073] pass Figure 5 The text visually demonstrates the dynamic comparison between the measured flow rate of the hydrological station and the corrected simulated flow rate output by the scheme provided in this embodiment over several consecutive years from 2015 to 2021.

[0074] It can be seen that the curve of the corrected simulated flow rate and the curve of the measured flow rate at the hydrological station maintain a high degree of consistency in their overall evolution trend. Facing the sudden change in flood peaks driven by extreme continuous precipitation during the high-water season, and the long-term baseflow recession during the dry season, the method in this embodiment demonstrates extremely high dynamic response capture capabilities. It is worth noting that this embodiment introduces an automatic phase alignment and nonlinear peak-shaving correction mechanism in the output stage, effectively compensating for the problems of flood peak sharpening and overestimation caused by the lack of dynamic regulation of floodplains in traditional one-dimensional physical models under extreme flood conditions, making the simulation results highly consistent with real physical laws.

[0075] In order to objectively quantify the simulation accuracy of the method in this embodiment, Figure 5 The study introduced the Nash efficiency coefficient (NSE) and Kling-Gupta efficiency coefficient (KGE), both recognized by the hydrological community, as core evaluation indicators.

[0076] The NSE value of this long-series daily runoff simulation reached 0.72, and the KGE value reached 0.80. These specific quantitative indicators fully demonstrate that, compared with traditional simplified models, the method of this invention achieves extremely high accuracy standards while maintaining strict water conservation.

[0077] In summary, this scheme overcomes the technical bottleneck of limited computing power in large-scale, high-resolution complex watersheds by traditional two-dimensional hydrodynamic models. Under the premise of ensuring strict water conservation and numerical stability, it provides a brand-new technical approach for urban flood prevention, watershed flood control emergency response, and disaster prevention and mitigation planning.

[0078] Furthermore, the terms "upper," "lower," "inner," "outer," "front," and "rear" are used for descriptive purposes only and should not be construed as indicating or implying relative importance. Unless otherwise specifically stated, the relative steps, numerical expressions, and values ​​of the components and steps set forth in these embodiments do not limit the scope of the invention.

[0079] Of course, the above description is only a specific embodiment of the present invention and is not intended to limit the scope of the present invention. All equivalent changes or modifications made to the structure, features and principles described in the claims of the present invention should be included in the scope of the claims of the present invention.

[0080] Finally, it should be noted that the above-described embodiments are merely specific implementations of the present invention, used to illustrate the technical solutions of the present invention, and not to limit it. The scope of protection of the present invention is not limited thereto. Although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that any person skilled in the art can still modify or easily conceive of changes to the technical solutions described in the foregoing embodiments within the technical scope disclosed in the present invention, or make equivalent substitutions for some of the technical features; and these modifications, changes, or substitutions do not cause the essence of the corresponding technical solutions to deviate from the spirit and scope of the technical solutions of the embodiments of the present invention, and should all be covered within the scope of protection of the present invention. Therefore, the scope of protection of the present invention should be determined by the scope of the claims.

Claims

1. A high-resolution flood simulation method coupling topology optimization motion waves and cellular automata, characterized in that, include: Step 1: Obtain watershed data for the target area and reproject the watershed data onto a high-resolution grid coordinate system to construct a grid diagram of upstream and downstream relationships based on flow direction. The watershed data includes a high-resolution digital elevation model, flow direction data, and catchment area data. Step 2: Traverse the upstream and downstream relationship grid at a preset time step, take the rainfall runoff collected at the current time as the lateral inflow, and construct a one-dimensional kinematic wave equation by combining the water flow state of the previous time step and the inflow state of the corresponding upstream grid, and solve the one-dimensional kinematic wave equation to obtain the river level at the current time step. Using the calculated river water level as the initial condition, a two-dimensional diffusion process is initiated based on the upstream and downstream relationship grid diagram to iteratively obtain the inundation range and water depth distribution map at the current time step. Step 3: Repeat Step 2 to obtain the inundation range and water depth distribution maps for all time steps.

2. The high-resolution flood simulation method based on coupled topology optimization of moving waves and cellular automata according to claim 1, characterized in that, The process of constructing the upstream and downstream relationship mesh diagram is as follows: The upstream and downstream connections between grids are established based on the raster flow direction data, and all effective grid cells are sorted in ascending topology according to the runoff area to construct the corresponding grid diagram. The river width of each grid cell is calculated using an empirical formula based on the catchment area, and the river slope and Manning roughness coefficient field of the entire region are initialized.

3. The high-resolution flood simulation method based on coupled topology optimization of moving waves and cellular automata according to claim 2, characterized in that, The river width is calculated based on the catchment area corresponding to the grid cell using an empirical power-law formula.

4. The high-resolution flood simulation method based on coupled topology optimization of moving waves and cellular automata according to claim 2, characterized in that, The process of the ascending topological sort is as follows: The runoff area data in two-dimensional raster form is flattened and converted into a one-dimensional array; An indirect sorting algorithm is applied to the one-dimensional array to arrange the flow area values ​​in ascending order, so as to obtain an integer index array for filling the grid cells.

5. The high-resolution flood simulation method based on coupled topology optimization of moving waves and cellular automata according to claim 1, characterized in that, The one-dimensional motion wave equation is solved using implicit finite difference and Newton-Raphson iteration.

6. The high-resolution flood simulation method based on coupled topology optimization of moving waves and cellular automata according to claim 5, characterized in that, The solution process is as follows: For each grid cell in the topology sequence, the actual flow length parameter is represented by a flow direction code; A nonlinear equation for a one-dimensional moving wave is constructed based on the mass conservation equation and the Manning momentum equation. The nonlinear equation is solved iteratively using the Newton-Raphson method. The iteration is repeated until the difference between two adjacent calculation results is less than the preset convergence threshold, so as to obtain the grid outflow rate at the current moment. Finally, the calculated grid outflow is added to the inflow of the downstream receiving unit of the grid at the current moment to complete the water transfer calculation for a single grid.

7. The high-resolution flood simulation method based on coupled topology optimization of moving waves and cellular automata according to claim 1, characterized in that, The two-dimensional diffuse process is as follows: Determine whether the current water depth of each target grid cell within the target area is greater than zero; Traverse the N neighboring grid cells of the target grid cell. If the water level of the neighboring grid cell is lower than that of the target grid cell, then determine that the corresponding neighboring grid cell is a potential water-receiving cell. Calculate the water level difference between the target grid cell and the potential water-receiving cell, and sum all water level differences to obtain the total water level difference. If the total water level difference is greater than zero, then the flow from the target grid to the [missing grid name] grid is calculated based on the principle of proportional distribution of elevation differences. Diffusion flux of a potential water-receiving unit; Subtract all calculated diffusion fluxes from the current water depth of the target grid cell and increase the water depth of each potential receiving cell accordingly.

8. The high-resolution flood simulation method based on coupled topology optimization of moving waves and cellular automata according to claim 1, characterized in that, The calculation process for the lateral inflow is as follows: At the start of each simulation time step, acquire gridded surface runoff depth rate data covering the entire region at the current moment; For each specific grid cell, the total runoff volume generated within the time step is calculated by combining the corresponding grid resolution size and the current time step. The total runoff volume is converted into the lateral inflow per unit river length in the one-dimensional kinematic wave equation.

Images