A two-dimensional waterlogging bidirectional coupling method based on GPU acceleration and multi-level physical restrictor

Abstract

The application discloses a one-dimensional and two-dimensional waterlogging two-way coupling method based on GPU acceleration and a multi-stage physical limiter, and comprises the following steps: obtaining drainage pipe network basic data of a target area and constructing one-dimensional drainage pipe network nodes and a space mapping dictionary; performing two-way coupling calculation on the drainage pipe network basic data; and repeatedly performing two-way coupling calculation to obtain corresponding waterlogging prediction data. The method provided by the application can solve the problems of numerical oscillation divergence and low CPU-GPU heterogeneous interaction efficiency caused by the traditional two-dimensional water dynamic model in explicit solving, so that the whole process of urban waterlogging overflow, spreading and water recession can be quickly deduced in a large range, high resolution and strict compliance with the law of conservation of energy.

Description

Technical Field

[0001] This invention belongs to the field of urban flood control and disaster reduction and smart water management technology, and in particular relates to a one- or two-dimensional bidirectional coupling method for urban flooding based on GPU acceleration and multi-level physical limiters. Background Technology

[0002] Existing technologies have the following limitations when performing pipeline-surface flood simulation and inundation early warning for large-scale, high-resolution urban complex terrain: Numerical oscillations and divergences in the explicit solution process of traditional numerical models: While the traditional full two-dimensional shallow water dynamic equations (SWEs) used for refined simulations are theoretically accurate, the massive matrix solutions in large-scale, high-resolution grid scenarios at the city level result in extremely time-consuming computations, making it difficult to meet the real-time simulation requirements of sudden rainstorms and flooding. To address this, simplified diffuse wave models or cellular automata (CA) models are often used for explicit time-step solutions to improve efficiency. However, in complex urban terrain (such as steep slopes) or extremely shallow water conditions, explicit solutions often lead to excessive water exchange between grids, causing severe "numerical oscillations" (i.e., abnormal rebounds of water flow between adjacent grids). This not only violates physical conservation laws but also easily leads to computational divergence and collapse of the entire model.

[0003] The existing anti-oscillation algorithms suffer from a mismatch in compatibility with the high-concurrency computing architecture of GPUs: To address the aforementioned numerical oscillation problem, existing high-order hydrodynamic models typically introduce complex approximate Riemannian solvers (such as the HLLC scheme) or employ extremely small adaptive / local time steps (LTS). These methods introduce a large amount of complex conditional branching logic and nonlinear equation iteration, completely destroying the parallelism of the algorithms and making them extremely difficult to adapt to the single-instruction multiple-data (SIMD) high-concurrency architecture of modern GPUs (Graphics Processing Units). This leads to the awkward technical bottleneck that "accurate calculations cannot be accelerated by GPUs, while calculations that can be accelerated by GPUs are extremely prone to oscillations."

[0004] Patent document CN119397862A discloses a method for simulating urban flooding based on CPU-GPU heterogeneous parallel computing. This method runs a one-dimensional pipe network model on the CPU and solves the two-dimensional shallow water equations on the GPU. Although heterogeneous computing is introduced, its core still relies on the traditional two-dimensional shallow water equations to describe the urban surface inundation process. When solving complex partial differential equations in parallel on the GPU, in order to maintain numerical stability, it is usually necessary to introduce complex Riemann solvers or extremely small time steps. This generates a large amount of conditional branching logic, which completely destroys the high concurrency advantage of the GPU's underlying single instruction multiple data (SIMD). If a simplified explicit solution is used under complex urban terrain with steep slopes or extremely shallow water conditions, it is extremely easy to cause excessive exchange of water between grids, leading to severe "numerical oscillations" or even model divergence and collapse.

[0005] Patent document CN121093812A discloses a threshold-optimized non-uniform grid method for simulating urban flooding based on surface features. This method uses machine learning to predict error thresholds and generates a dynamic non-uniform grid to reduce computational load. However, this scheme explicitly states that it uses an "integrated operation of a dynamic non-uniform grid and a surface-pipeline unidirectional coupling model." Unidirectional coupling only considers the overflow from the pipeline network to the surface due to overload, severely lacking the physical closed loop of "return flow" of surface water to underground pipelines during the later stages of heavy rainfall. In addition, existing explicit models based on grid state updates generally lack rigorous physical extreme value constraints. When dealing with extreme elevation differences beyond the normal range, they are prone to distortions that violate mass conservation and physical common sense, such as local grid "negative water depth" or "downstream water level abnormally exceeding upstream." Summary of the Invention

[0006] The purpose of this invention is to provide a one- or two-dimensional bidirectional coupling method for urban flooding based on GPU acceleration and multi-level physical limiters. This method can solve the problems of numerical oscillation and divergence that are easily caused by explicit solutions in traditional two-dimensional hydrodynamic models, as well as the low efficiency of CPU-GPU heterogeneous interaction. Thus, it can realize rapid simulation of the entire process of urban flooding "overflow-spreading-receding" on a large scale and at high resolution, while strictly following the physical conservation laws.

[0007] To achieve the objectives of this invention, the following technical solution is provided: a one- or two-dimensional flood bidirectional coupling method based on GPU acceleration and multi-level physical limiters, comprising the following steps: Step 1: Obtain the basic data of the drainage pipe network in the target area, as well as the corresponding high-resolution digital elevation model and rainfall and meteorological data, and establish a spatial mapping dictionary between the one-dimensional drainage pipe network nodes and the corresponding two-dimensional land surface in the basic data of the drainage pipe network. Step 2: Under a preset coupling time step, based on the high-resolution digital elevation model, spatial mapping dictionary, and rainfall and meteorological data, the rainfall and meteorological data are converted into surface runoff input using a hydrological runoff generation and confluence mechanism, and bidirectional coupling calculations are performed on the basic data of the drainage network. When a one-dimensional pipeline node overflows, the pipeline overflow volume is converted into an additional source term for the two-dimensional surface to synchronously update the two-dimensional surface water depth matrix at the current time step. When the one-dimensional pipe network is not overflowing and there is water accumulation on the ground surface, the backflow flow of the water flowing from the two-dimensional ground surface back to the underground pipe network is calculated based on the multi-level physical limiter, and used as the lateral inflow of the one-dimensional pipe network in the next time step to update the time series water level status. Step 3: Repeat the above bidirectional coupling calculation process in step 2 until the full-time simulation ends. Extract the two-dimensional surface water depth matrix and the time-series water level status of the target node at each time step according to the preset output frequency, and convert them into a spatial raster layer and node time-series sequence to obtain the corresponding urban flooding prediction data.

[0008] Specifically, the basic data of the drainage network includes land use roughness data of the underlying surface and one-dimensional network topology data.

[0009] Specifically, the construction process of the space mapping dictionary is as follows: Obtain the node coordinates of the one-dimensional drainage network nodes in the basic data of the drainage network; The node coordinates are converted into row and column indices in a raster matrix using a high-resolution digital elevation model. It then determines whether each index point in the raster matrix and its corresponding grid point in the two-dimensional map are within the valid computational domain. If they are within the valid computational domain, the row and column indices of the corresponding index point are recorded; otherwise, they are left blank, thus constructing a spatial mapping dictionary corresponding to the node coordinates.

[0010] Specifically, the process of transforming the overflow state into an additional source term is as follows: Based on each active node in the spatial mapping dictionary, determine whether the current overflow flow of the corresponding node in the one-dimensional drainage network simulation is greater than zero; If the value is greater than zero, then based on the current time step and the corresponding two-dimensional grid area, the equivalent value of converting the single-step overflow volume into the water depth increment is calculated, and the calculated equivalent value is accumulated into the water depth of the two-dimensional surface grid with the corresponding row and column index.

[0011] Specifically, the calculation process for the return flow rate of the wastewater is as follows: For each node in the spatial mapping dictionary, extract the current surface water depth of the corresponding two-dimensional grid. When the current surface water depth corresponding to a node is greater than the preset minimum calculated water depth, and the sum of the current grid absolute elevation and the surface water depth is greater than the wellhead elevation of the corresponding pipeline node, the theoretical return flow rate is calculated using the surface weir flow formula, and an upper limit constraint is applied in conjunction with the current total water storage of the grid to output the actual return flow rate.

[0012] Specifically, the two-dimensional surface water depth matrix is ​​updated through two-dimensional surface diffusion evolution, and the process is as follows: In the CUDA kernel function of the GPU, all two-dimensional effective grid cells are traversed concurrently, the water level difference and slope from the central grid cell to the corresponding 8 neighboring grid cells are calculated, the sum of the positive slopes and the maximum water level difference are calculated. If the sum of the positive slopes is greater than zero, the final outflow volume constrained by multi-level physical limiters is calculated to simulate the process of water spreading to the neighborhood, thereby obtaining the updated surface water depth at the current time step.

[0013] Specifically, the multi-level physical limiter includes a theoretical flow rate limit, a physical limit to prevent dewatering, and a water level reversal limit.

[0014] Specifically, the flooding prediction data includes the surface inundation range, water depth distribution map, and pipeline node status data.

[0015] Compared with the prior art, the beneficial effects of the present invention are as follows: By introducing a multi-level physical limiter consisting of "Manning velocity limit, anti-drainage physical limit, and water level reversal limit", this method completely eliminates water level oscillation and negative water depth phenomena between large slope or shallow water grids from the underlying physical mechanism. The extremely simple "minimum value discrimination" logic greatly improves the robustness of the algorithm, thereby breaking through the computational bottleneck of explicit models being prone to oscillation and divergence.

[0016] By establishing an efficient memory / video memory mapping mechanism between the CPU running pipeline and the GPU running on the ground, high-frequency and ultra-fast interaction of one- and two-dimensional model data is realized. With the support of oscillation-free algorithms, the GPU can perform massive grid state concurrent calculations under non-branching instruction flow, and solve video memory competition with the underlying atomic operations, realizing ultra-fast simulation of large-scale, high-resolution urban flooding spread process.

[0017] By introducing a drainage control boundary based on the surface weir flow formula, the physical process of surface water receding into the underground drainage network under the action of gravity during the rainstorm remission period was accurately simulated, which effectively improved the realism and early warning accuracy of long-term urban flooding full life cycle simulation. Attached Figure Description

[0018] Figure 1 This is an overall flowchart of a one- or two-dimensional bidirectional coupling method for waterlogging based on GPU acceleration and multi-level physical limiters provided in this embodiment; Figure 2 This is a schematic diagram of the one- or two-dimensional spatial mapping and initial water volume interaction provided in this embodiment; Figure 3 This is a schematic diagram illustrating the principle, control effect comparison, and GPU parallel acceleration architecture of the non-oscillating cellular automaton algorithm based on multi-level physical limiters provided in this embodiment. Figure 4 This is a verification diagram of the dynamic drainage mechanism of the surface-pipeline network and the two-way water volume interaction based on the surface weir flow formula provided in this embodiment. Figure 5 This is a time-series evolution diagram of the water depth submersion of the target street under heavy rain conditions provided in this embodiment. Detailed Implementation

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

[0020] like Figure 1 As shown in this embodiment, a one- or two-dimensional flood bidirectional coupling method based on GPU acceleration and multi-level physical limiters is provided, including the following steps: Step 1: Obtain the basic data of the drainage pipe network in the target area, as well as the corresponding high-resolution digital elevation model and rainfall and meteorological data, and establish a spatial mapping dictionary between the one-dimensional drainage pipe network nodes and the corresponding two-dimensional land surface in the basic data of the drainage pipe network.

[0021] Step 2: Under a preset coupling time step, based on the high-resolution digital elevation model, spatial mapping dictionary, and rainfall and meteorological data, the rainfall and meteorological data are converted into surface runoff input using a hydrological runoff generation and confluence mechanism, and bidirectional coupling calculations are performed on the basic data of the drainage network. When a one-dimensional pipeline node overflows, the pipeline overflow volume is converted into an additional source term for the two-dimensional surface to synchronously update the two-dimensional surface water depth matrix at the current time step. When the one-dimensional pipe network is not overflowing and there is water accumulation on the ground surface, the backflow of water flowing from the two-dimensional surface back to the underground pipe network is calculated based on the multi-level physical limiter, and used as the lateral inflow of the one-dimensional pipe network in the next time step to update the time series water level status.

[0022] Step 3: Repeat the above bidirectional coupling calculation process in Step 2 until the full-time simulation ends. Extract the two-dimensional surface water depth matrix and the time-series water level status of the target node at each time step according to the preset output frequency, so as to convert them into a spatial raster layer and node time-series sequence, thereby obtaining the corresponding urban flooding prediction data.

[0023] The solution provided in this embodiment is based on a CPU-GPU heterogeneous computing architecture to construct a two-way coupled model of urban flooding. It utilizes the affine transformation of a digital elevation model and a one-dimensional pipe network topology to establish an accurate spatial mapping mechanism. Secondly, a one- or two-dimensional water volume interaction module based on node states is constructed to receive overloaded overflow input from the underground pipe network and transform it into a source term for surface spread. Subsequently, an oscillation-free two-dimensional cellular automaton surface spread module is established. A multi-level physical limiter, consisting of the Manning theoretical limit, the anti-drainage physical limit, and the water level reversal limit, is introduced into the GPU parallel thread to strictly constrain the volume of water flowing out in a single step between grids. Synchronous updates of water depth are performed through atomic operations in the GPU memory, ensuring extreme value stability and conservation of mass and potential energy during surface flooding. Finally, a dynamic water return mechanism from the surface to the pipe network based on the weir flow formula is established, forming a complete closed-loop water exchange. This method effectively balances the physical accuracy, computational efficiency, and numerical stability of complex urban underlying surface simulation, significantly improving the ability to refine the full life-cycle simulation of large-scale urban flooding disasters.

[0024] More specifically, the construction process of the spatial mapping dictionary is as follows: read the node coordinate data in the configuration file of the one-dimensional pipeline network model; use the spatial affine inverse transformation matrix of the digital elevation model to convert the coordinates of each node into row and column indices in the two-dimensional raster matrix; determine whether the two-dimensional grid corresponding to the index is within the valid computational domain. If valid, pair the node ID with the corresponding row and column indices to generate the node spatial mapping dictionary.

[0025] The two-dimensional surface water depth matrix is ​​obtained by determining whether the current overflow flow rate of each active node in the spatial mapping dictionary is greater than zero in the one-dimensional pipe network simulation. If there is an overflow, the equivalent value of the single-step overflow volume converted into the water depth increment is calculated based on the current time step and the corresponding two-dimensional grid area, and then accumulated into the water depth of the two-dimensional surface grid with the corresponding row and column index.

[0026] The calculation process for the return flow of the drainage is as follows: For each node in the spatial mapping dictionary, the current surface water depth of its corresponding two-dimensional grid is extracted; when the surface water depth is greater than the preset minimum calculation water depth, and the sum of the current grid absolute elevation and the surface water depth is greater than the wellhead elevation of the corresponding pipeline node, the drainage calculation is initiated; the theoretical return flow is calculated using the surface weir flow formula, and an upper limit constraint is applied in conjunction with the current total water storage of the grid, with the actual drainage flow being the smaller of the two values; the theoretical return flow is the product of the weir flow coefficient, the weir width, and the surface water depth to the power of 1.5; the value of the upper limit constraint is the equivalent flow value obtained by multiplying the surface water depth and the grid area by the time step; the calculated actual drainage flow is converted into inflow boundary conditions and assigned to the corresponding one-dimensional pipeline node, and the corresponding water depth equivalent value is deducted from the two-dimensional surface grid.

[0027] The process of two-dimensional surface spread evolution using an oscillatory cellular automaton algorithm based on multi-level physical limiters is as follows: In the CUDA kernel function of the GPU, all effective two-dimensional grid cells are concurrently traversed, the water level difference and slope from the central grid cell to its 8 neighboring grid cells are calculated, the sum of the positive slopes and the maximum water level difference are calculated; if the sum of the positive slopes is greater than zero, the final outflow volume constrained by the multi-level physical limiters is calculated; the multi-level physical limiters are composed of the theoretical flow velocity limit, the anti-drainage physical limit and the water level reversal limit, and the final outflow volume is the minimum value of these three calculated volumes.

[0028] The calculation process for the final outflow volume under the constraints of multi-level physical limiters is as follows: First-level limit: Calculate the theoretical maximum flow rate under the Manning formula based on the water depth, Manning roughness, and average positive slope of the central grid, and multiply by the time step to obtain the theoretical outflow volume; Second-level limit: Based on the current volume of the central grid, set an upper limit for the anti-drainage volume, which is 50% of the product of the grid water depth and area, to ensure that mass conservation is not violated and to prevent negative water depths; Third-level limit: Based on the current maximum water level difference, set an upper limit for the anti-level reversal volume, which is 25% of the product of the maximum water level difference and the grid area, to prevent excessive single-step loss leading to downstream water levels surpassing upstream levels; Finally, the final outflow volume constrained by the multi-level physical limiters is determined as the minimum value among the theoretical outflow volume, the anti-drainage volume upper limit, and the anti-level reversal volume upper limit.

[0029] The parallel update process of two-dimensional surface water depth is as follows: Within the GPU computing framework, the final outflow volume is proportionally allocated to the diffusion flux weight of each outflow direction based on the proportion of the positive slope in the eight neighborhood directions; GPU memory atomic operations are called to perform the deduction of the water depth volume of the central grid and the increase of the water depth volume of the neighboring receiving grid, so as to avoid data contention for memory read and write during multi-threaded concurrent execution, and complete the evolution of two-dimensional surface water flow at the current time step.

[0030] To better illustrate the solution in this embodiment, the steps are explained by analyzing the flooding of a target street in a certain area during a rainstorm.

[0031] Step 1: Data acquisition and establishment of one- or two-dimensional spatial mapping.

[0032] First, basic data for the target area is acquired, including high-resolution digital elevation model (DEM) data, underlying land use roughness data, and a one-dimensional drainage network topology configuration file. Then, the node coordinate data in the one-dimensional drainage network configuration file is read to establish a spatial mapping relationship between the one-dimensional drainage network nodes and the two-dimensional surface grid. The specific mapping establishment process is as follows... Figure 2 As shown in (A): Using the spatial affine inverse transformation matrix of the digital elevation model, the coordinates of each one-dimensional pipeline node are ( , Convert the data to row and column indices in a two-dimensional raster matrix: ; in, For rotation parameters, These are the translation parameters.

[0033] Determine whether the two-dimensional grid corresponding to the index is within the valid computation domain. If valid, pair the node ID with the corresponding row and column index to generate a node space mapping dictionary, which serves as the index basis for subsequent one-dimensional and two-dimensional bidirectional water volume exchange.

[0034] Step 2: Within the preset coupling time step Next, initiate bidirectional coupling calculation between the one-dimensional pipeline network and the two-dimensional surface.

[0035] This process is completed through heterogeneous parallel collaboration between CPU and GPU. The specific iterative process includes: 2.1 Extracting one-dimensional overflow and converting it into two-dimensional surface water depth source terms (1D-2D). For each active node in the spatial mapping dictionary, the overflow status of the next-dimensional drainage network node in the current time step is extracted.

[0036] The physical processes of water interaction are as follows Figure 2 As shown in (B): Determine its current overflow flow rate. Is it greater than zero? If overflow exists, then based on the current time step. and the corresponding two-dimensional grid area Calculate the equivalent value of a single-step overflow volume converted into a water depth increment. And add it to the initial water depth of the two-dimensional surface grid corresponding to the row and column indices: ; 2.2 Two-dimensional performance of oscillatory cellular automata based on multi-level physical limiters.

[0037] After obtaining the updated two-dimensional surface water depth source term from step 2.1, the GPU-accelerated non-oscillating cellular automata (CA) spreading module is initiated. Its cellular network interaction principle, algorithm comparison, and hardware acceleration mechanism are as follows: Figure 3 As shown: First, establish cell interaction logic based on 8 neighborhoods (such as...). Figure 3 (As shown in (A)). Each two-dimensional grid is considered as an independent cell. Its properties include roughness. , water depth and absolute water level Within each computation step, the GPU concurrently traverses all valid cells and computes the central cell. and its 8 neighboring cells water level difference : When the water level in the central cell is higher than that in the neighboring cells and the slope is greater than zero, the water flow diffuses into the neighboring cells. To address the numerical oscillation problem that traditional explicit solution models are prone to under steep slopes or extremely shallow water conditions, a comparative analysis is shown below. Figure 3 (B) introduces a multi-stage physical limiter consisting of the theoretical flow velocity limit, the anti-drainage physical limit, and the water level reversal limit, which strictly constrains the final outflow volume in a single step. : ; in: Level 1 (Manning's theoretical limit): Based on the water depth of the central grid. Manning roughness and average positive slope Calculate the theoretical maximum volume: ; Level 2 (Anti-drainage limit): To prevent negative water depth, the upper limit is set to 50% of the current grid water volume. ; Level 3 (Water Level Reversal Limit): To prevent numerical fluctuations, the upper limit is set at 25% of the volume corresponding to the maximum water level drop. .

[0038] like Figure 3As shown in (C), the above algorithm logic is fully adapted to the parallel computing architecture of GPUs. In the CUDA kernel function, all threads execute unbranched instruction streams (SIMD), that is, uniformly execute the "minimum value" discrimination logic, thereby maximizing the utilization of computing power. When updating the water depth data in the video memory, atomic operations are used to perform the central grid volume reduction and the neighboring grid volume increase, effectively avoiding the video memory race condition during multi-threaded concurrent read and write, and realizing the stable and rapid evolution of the planar water flow.

[0039] 2.3 Calculate the return flow of water from the two-dimensional surface to the one-dimensional pipe network (2D - 1D).

[0040] After the two-dimensional performance at each time step is completed, the model establishes a dynamic drainage mechanism based on the surface weir flow formula to form a complete closed loop of water exchange. Its drainage triggering conditions, flow calculation logic, and water balance verification are as follows: Figure 4 As shown: Extracted updated two-dimensional surface water depth ,when The water depth is greater than the preset minimum calculation depth, and the current grid's absolute water level (elevation) is... + (The elevation of the wellhead is greater than that of the corresponding pipeline node) At that time, the drainage calculation is initiated. The actual drainage flow rate is calculated using the surface weir flow formula. Its value is determined by the following formula: in, The weir flow coefficient is... The effective weir width. For example... Figure 4 As shown in (B), the actual discharge flow rate was calculated. The inflow boundary conditions are transformed in real time and assigned to the corresponding one-dimensional pipe network nodes (such as the inflow term in the SWMM model), while the corresponding water depth equivalent value is subtracted from the corresponding two-dimensional surface grid. .

[0041] Based on a comparison of typical time-series timelines, it is evident that during periods of respite from heavy rainfall, surface water flows back into the underground pipe network due to gravity via this mechanism. As the simulation progresses, the surface water depth... Smooth decrease, while pipeline node flow (As the return receiver) the amount increases accordingly. This mechanism not only accurately simulates the real physical receding process, but also ensures the dynamic balance of water volume in the system during long-duration rainstorm simulations.

[0042] Step 3: Advance the time step and repeat the process in Step 2.

[0043] The bidirectional coupling calculation is repeated continuously until the set total simulation duration is reached. After reaching a specific output node, the full-time surface inundation range, water depth distribution map, and pipeline node status data are output, providing refined data support for urban flooding early warning, and ultimately forming a complete, detailed, and physically reliable dataset of the flood evolution process.

[0044] In this embodiment, the step of converting rainfall meteorological data into surface runoff input using hydrological runoff generation and confluence mechanisms can be specifically implemented based on a stormwater and flood management model (SWMM). Specifically, the surface runoff generation process uses the Horton infiltration model to calculate soil infiltration and net rainfall, while the pipe network confluence process uses the Dynamic Wave Evolution equation to calculate nodal head and pipe flow. Through this specific configuration, complex unsteady flow phenomena and waterlogging characteristics such as pipe network backflow and flooding can be accurately simulated.

[0045] Its simulation process is as follows Figure 5 As shown in (A), in the early stage of rainfall, the model accurately captures the initial overflow phenomenon at the rainwater pipe points caused by the one-dimensional pipe network reaching its limit.

[0046] At this time, the surface water is distributed in a point-source pattern, and the water depth is relatively shallow (mostly between 0.1 and 0.2 meters).

[0047] like Figure 5 As shown in (B), as rainfall enters its middle stage, surface water, driven by the CA algorithm based on multi-level physical limiters described in step 2.2, begins to rapidly evolve along the street slope and building edges towards surrounding low-lying areas. It can be seen that the waterlogged areas gradually connect from isolated points into a linear pattern, and the water depth at road intersections increases significantly. During this process, due to the constraints of the three-level limiters, the water flow boundary is smooth, and no numerical oscillations or abnormal rebounds occur.

[0048] like Figure 5 As shown in (C), the flooded area reaches its maximum when rainfall reaches its peak. The dark red area (water depth > 0.5 meters) clearly delineates the high-risk points for flooding in this area. At this time, the one-dimensional pipe network system and the two-dimensional surface system maintain a dynamic water balance through the drainage and return mechanism described in step 2.3.

[0049] In summary, by using a GPU parallel acceleration architecture, while maintaining centimeter-level high-resolution simulation, the computational efficiency is improved by tens of times compared to traditional methods. Furthermore, it fully simulates the entire life cycle physical process from "pipeline overflow - surface spread - gravity receding water", providing a reliable technical means for real-time and accurate early warning of urban flooding.

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

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

[0052] 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 one- or two-dimensional bidirectional coupling method for flooding based on GPU acceleration and multi-level physical limiters, characterized in that, include: Step 1: Obtain the basic data of the drainage pipe network in the target area, as well as the corresponding high-resolution digital elevation model and rainfall and meteorological data, and establish a spatial mapping dictionary between the one-dimensional drainage pipe network nodes and the corresponding two-dimensional land surface in the basic data of the drainage pipe network. Step 2: Under a preset coupling time step, based on the high-resolution digital elevation model, spatial mapping dictionary, and rainfall and meteorological data, the rainfall and meteorological data are converted into surface runoff input using a hydrological runoff generation and confluence mechanism, and bidirectional coupling calculations are performed on the basic data of the drainage network. When a one-dimensional pipeline node overflows, the pipeline overflow volume is converted into an additional source term for the two-dimensional surface to synchronously update the two-dimensional surface water depth matrix at the current time step. When the one-dimensional pipe network is not overflowing and there is water accumulation on the ground surface, the backflow flow of the water flowing from the two-dimensional ground surface back to the underground pipe network is calculated based on the multi-level physical limiter, and used as the lateral inflow of the one-dimensional pipe network in the next time step to update the time series water level status. Step 3: Repeat the above bidirectional coupling calculation process in Step 2 until the full-time simulation ends. Extract the two-dimensional surface water depth matrix and the time-series water level status of the target node at each time step according to the preset output frequency, so as to convert them into a spatial raster layer and node time-series sequence, thereby obtaining the corresponding urban flooding prediction data.

2. The one- or two-dimensional flood bidirectional coupling method based on GPU acceleration and multi-level physical limiters according to claim 1, characterized in that, The basic data of the drainage network includes land use roughness data of the underlying surface and one-dimensional network topology data.

3. The one- or two-dimensional flood bidirectional coupling method based on GPU acceleration and multi-level physical limiters according to claim 1, characterized in that, The construction process of the spatial mapping dictionary is as follows: Obtain the node coordinates of the one-dimensional drainage network nodes in the basic data of the drainage network; The node coordinates are converted into row and column indices in a raster matrix using a high-resolution digital elevation model. It then determines whether each index point in the raster matrix and its corresponding grid point in the two-dimensional map are within the valid computational domain. If they are within the valid computational domain, the row and column indices of the corresponding index point are recorded; otherwise, they are left blank, thus constructing a spatial mapping dictionary corresponding to the node coordinates.

4. The one- or two-dimensional flood bidirectional coupling method based on GPU acceleration and multi-level physical limiters according to claim 1, characterized in that, The process of transforming the overflow state into an additional source term is as follows: Based on each active node in the spatial mapping dictionary, determine whether the current overflow flow of the corresponding node in the one-dimensional drainage network simulation is greater than zero; If the value is greater than zero, then based on the current time step and the corresponding two-dimensional grid area, the equivalent value of converting the single-step overflow volume into the water depth increment is calculated, and the calculated equivalent value is accumulated into the water depth of the two-dimensional surface grid with the corresponding row and column index.

5. The one- or two-dimensional flood bidirectional coupling method based on GPU acceleration and multi-level physical limiters according to claim 1, characterized in that, The calculation process for the return flow rate of the wastewater is as follows: For each node in the spatial mapping dictionary, extract the current surface water depth of the corresponding two-dimensional grid. When the current surface water depth corresponding to a node is greater than the preset minimum calculated water depth, and the sum of the current grid absolute elevation and the surface water depth is greater than the wellhead elevation of the corresponding pipeline node, the theoretical return flow rate is calculated using the surface weir flow formula, and an upper limit constraint is applied in conjunction with the current total water storage of the grid to output the actual return flow rate.

6. The one- or two-dimensional flood bidirectional coupling method based on GPU acceleration and multi-level physical limiters according to claim 1, characterized in that, The two-dimensional surface water depth matrix is ​​updated through two-dimensional surface diffusion evolution, and the process is as follows: In the CUDA kernel function of the GPU, all two-dimensional effective grid cells are traversed concurrently, the water level difference and slope from the central grid cell to the corresponding 8 neighboring grid cells are calculated, the sum of the positive slopes and the maximum water level difference are calculated. If the sum of the positive slopes is greater than zero, the final outflow volume constrained by multi-level physical limiters is calculated to simulate the process of water spreading to the neighborhood, thereby obtaining the updated surface water depth at the current time step.

7. The one- or two-dimensional flood bidirectional coupling method based on GPU acceleration and multi-level physical limiters according to claim 1, characterized in that, The multi-level physical limiters include the theoretical flow rate limit, the anti-drainage physical limit, and the water level reversal limit.

8. The one- or two-dimensional flood bidirectional coupling method based on GPU acceleration and multi-level physical limiters according to claim 1, characterized in that, The flooding prediction data includes the surface inundation range, water depth distribution map, and pipeline node status data.

Images