A method for tracing sources of urban rainstorm waterlogging multi-source pollution based on bidirectional dynamic coupling

Abstract

This invention discloses a method for tracing the source of multi-source pollution in urban stormwater flooding based on bidirectional dynamic coupling. First, multi-source data is acquired. Then, the computational domain of the two-dimensional shallow water equation solver ANUGA is constructed and initialized, i.e., the surface hydrodynamic model is initialized. Next, a bidirectional dynamic coupling model of the surface and the pipe network is constructed. Subsequently, water quality transport analysis is performed, and by constructing a pollution source fingerprint database, qualitative identification and quantitative source tracing analysis of pollution sources are completed based on the multi-source water quality tracing model. This invention achieves integrated simulation and analysis of the entire chain of surface hydrodynamics, drainage network hydraulics, water quality transport, and pollution source tracing, solving problems such as model decoupling, missing water quality dimensions, imperfect tracing methods, and lack of a systematic tracing framework. It meets the engineering requirements for high-precision, real-time wastewater source tracing in urban stormwater flooding scenarios.

Description

Technical Field

[0001] This invention relates to the field of data processing technology for urban rainstorm and flooding analysis, specifically a method for tracing the source of multi-source pollution from urban rainstorm and flooding based on bidirectional dynamic coupling. Background Technology

[0002] As the impermeability of urban surfaces continues to rise and the frequency and intensity of extreme rainfall events increase significantly, urban surface runoff and drainage pipe network systems are often overloaded simultaneously under heavy rainfall conditions, resulting in a large amount of sewage overflowing from the pipe network to the surface, forming a complex urban flooding disaster with synergistic harm of "water quality and water quantity".

[0003] In this scenario, the overflowing wastewater contains high levels of chemical oxygen demand (COD) and ammonia nitrogen. ), nitrate nitrogen ( Pollutants such as total phosphorus (TP) and suspended solids (TSS) can easily cause black and odorous water bodies in cities, secondary pollution of bottom sediments, and public health risks. How to quickly and accurately track the source location, source type, and contribution ratio of sewage overflows during or after urban flooding events has become a key technical requirement for urban water environment management, emergency response, and post-event assessment.

[0004] After decades of development, the existing pollution source tracing technology system has formed five mainstream technical routes, represented by receptor model method, isotope tracing method, multivariate statistical method, numerical model method and machine learning method. The above five mainstream technical routes have the following main defects: (1) Model decoupling, lacking a two-way dynamic coupling mechanism between surface and pipe network. The surface model and drainage pipe network model run independently, which cannot reflect the substantial contribution of pipe network overflow to surface water dynamics in real time, nor can it capture the backwater effect of surface water accumulation on pipe network, resulting in distorted simulation results; (2) Lack of water quality dimension. The existing coupling framework only simulates water volume and lacks the ability to simulate the multi-component pollutants (COD, etc.) carried in the overflow water body. , (2) The ability to simulate the transport of pollutants (TP, TSS, etc.) is insufficient to quantitatively assess the impact of urban flooding on urban water environment quality; (3) The pollutant source tracing methods are imperfect: existing pollution source tracing technologies are mostly used for point source emission assessment, lacking real-time source tracing methods for dynamic overflow of multiple outlets in the urban flooding scenario, resulting in low computational efficiency and inability to support the rapid source tracing needs under large-scale parallel high-resolution simulation; (4) There is a lack of a systematic pollution fingerprint database and source tracing analysis framework, making it impossible to distinguish the relative contribution ratio of pollutants from different sources such as domestic sewage, industrial wastewater, and agricultural runoff. Summary of the Invention

[0005] The technical problem to be solved by this invention is to provide a method for tracing the source of multi-source pollution in urban rainstorms and floods based on bidirectional dynamic coupling. This method realizes the integrated simulation and analysis of the entire chain of surface water dynamics, drainage network hydraulics, water quality transport and pollution source tracing, and solves the problems of model decoupling, missing water quality dimensions, imperfect source tracing methods and lack of a systematic source tracing framework. It meets the engineering requirements of high-precision and real-time sewage source tracing in urban rainstorm and flood scenarios.

[0006] The technical solution of this invention is as follows: A method for tracing the sources of multi-source pollution from urban stormwater flooding based on bidirectional dynamic coupling includes the following steps: (1) Acquisition of multi-source data: Collect DEM data and drainage network node data of the area to be studied, and construct rainfall driving sequences; (2) Initialize the surface hydrodynamic model: Construct and initialize the computational domain of the two-dimensional shallow water equation solver ANUGA; (3) Construct a two-way dynamic coupling model between the surface and the pipeline network: The overflow of the SWMM node is transmitted to the surface hydrodynamic model, and the surface hydrodynamic model feeds back the water surface elevation increment to the storm flood management model SWMM, forming a two-way dynamic coupling, and thus constructing a two-way dynamic coupling model between the surface and the pipeline network. (4) Conduct water transport analysis: Initialize the concentration of various pollutants in the study area on the surface hydrodynamic model, and use the water multi-component transport module to dynamically update the pollutant concentration; (5) Pollution source tracing analysis: Construct a pollution source fingerprint database and a multi-pollution source water quality tracing model based on conservative tracers to conduct qualitative identification and quantitative tracing analysis of pollution sources.

[0007] The DEM data is in GeoTIFF format; The drainage network node data is obtained by transforming the original coordinates of the drainage network nodes from the WGS-84 geographic coordinate system to the projected coordinate system corresponding to the area under study. The transformation process is shown in the following formula (1): (1); In equation (1), Represents longitude in the WGS-84 geographic coordinate system. Represents latitude in the WGS-84 geographic coordinate system. This represents the X-axis coordinate value in the UTM projected coordinate system. This represents the Y-axis coordinate value in the UTM projected coordinate system. This represents the transformation function from latitude and longitude to UTM projected X-axis coordinates. This represents the transformation function from latitude and longitude to the UTM projected Y-axis coordinates; The construction of the rainfall driving sequence first involves collecting rainfall observation data from meteorological stations or simulating future precipitation data using atmospheric numerical models. Then, the rainfall observation data or future precipitation data are arranged in time sequence and the units are converted. The specific unit conversion formula is shown in the following formula (2): (2); In equation (2), Represents observed rainfall data or future precipitation data, in units of... ; Represents rainfall observation data or future precipitation data after unit conversion, in units of ; When the duration of rainfall observation data or future precipitation data is shorter than the total simulation duration, the rainfall observation data or future precipitation data for the excess simulation time will be set to zero to ensure the integrity of the rainfall driving sequence.

[0008] The initialization of the surface hydrodynamic model specifically includes the following steps: S21. Within the spatial coverage of the DEM data, the finite volume method is used to generate an unstructured triangular computational grid in the area to be studied, and the computational domain of the two-dimensional shallow water equation solver ANUGA is constructed. The maximum triangle area in the unstructured triangular computational grid is constrained by the following formula (3). (3); In equation (3), This represents the basic resolution, measured in meters (m). Equation (3) represents the area of ​​the largest triangle; Equation (3) represents that the area of ​​a single triangle grid does not exceed the area of ​​the square corresponding to the edge of the basic resolution. S22. Set the computational domain parameters, specifically by using the flow algorithm to calculate the water flow wave velocity and constructing the CFL stability conditions, as shown in the following formula (4): (4); In equation (4), Represents the number of CFLs; The water flow wave velocity calculated by the representative flow algorithm; The time step representing the surface hydrodynamic model; This represents the size of the triangular grid, i.e., the basic resolution. Equation (4) represents the stability of the numerical calculation of the surface hydrodynamic model by controlling the ratio of water flow wave velocity, time step and triangular mesh size; S23. The DEM data is assigned to the centroid of all triangular meshes using spatial interpolation, and a smoothing coefficient is set to suppress interpolation noise. The initial water surface elevation of the surface hydrodynamic model is calculated by the following formula (5): (5); In equation (5), Represents the initial water surface elevation; This represents the bed elevation after spatial interpolation of the DEM data; This represents the initial water depth under thin water film conditions.

[0009] The construction of the bidirectional dynamic coupling model of the land surface and the pipeline network specifically includes the following steps: S31. In each coupling time step, the surface water depth at all SWMM node coordinates is calculated in batches using spatial interpolation, and the surface water depth at all SWMM node coordinates is assigned to the corresponding SWMM node. The formula for calculating the surface water depth at the SWMM node coordinates is shown in the following formula (6): (6); In equation (6), Represents SWMM node coordinates The surface water depth at that location is the node in the SWMM stormwater and flood management model. represent The water surface elevation at the location is calculated using spatial interpolation based on the water surface elevation output from the surface hydrodynamic model. represent The surface elevation at that location; This represents taking the maximum value to ensure that the surface water depth is not negative. S32. For the SWMM nodes that overflow in the storm flood management model, calculate the total overflow volume according to the coupling time step. With each SWMM node that overflows as the center and a circular area with a set radius as the overflow influence range, uniformly superimpose the total overflow volume onto the water surface elevation of the corresponding triangular grid within the overflow influence range. This transmits the location of the SWMM node that overflows and the total overflow volume to the surface hydrodynamic model, thus constructing a two-way dynamic coupling model between the surface and the pipe network. The total overflow volume The calculation formula is shown in the following formula (7): (7); In equation (7), This represents the total overflow volume generated by a single SWMM node overflow within the coupling time step; Represents the instantaneous overflow of a single SWMM node (m³ / s); Represents the coupling time step; The water surface elevation increment caused by overflow in the surface hydrodynamic model The calculation formula is shown in the following formula (8): (8); In equation (8), The area representing the extent of the overflow's influence, i.e., the area of ​​a circular region with a set radius; Steps S33, S31, and S32 are executed cyclically at each coupling time step. That is, at each coupling time step, the total overflow volume of the stormwater flood management model (SWMM) is transferred to the surface hydrodynamic model, and then the surface hydrodynamic model feeds back the water surface elevation increment to the stormwater flood management model (SWMM), forming a two-way dynamic coupling.

[0010] The water transport analysis specifically includes the following steps: S41, for each pollutant Create a concentration array with the same length as the total number of triangular meshes in the surface hydrodynamic model. , Represents chemical oxygen demand. Represents ammonium ions. Represents nitrate ions. Represents total phosphorus, Represents suspended matter; wherein, the initial concentration of each pollutant is determined by the background concentration of the pollutant, as shown in the following formula (9): (9); In equation (9), This represents the pollutants in the i-th triangular grid. The initial concentration at the initial time t=0; Represents the background concentration of pollutants; This represents the total number of triangular meshes in the surface hydrodynamic model; S42. Construct the kinetic equations for each pollutant, as shown in equations (10) to (14) below: (10); In equation (10), represent Rate of change of concentration over time; represent At ℃, The first-order decay rate constant of the concentration; Represents the temperature coefficient. Represents the actual water temperature; Represents actual water temperature Real-time concentration; The negative sign represents Concentration decreases over time; (11); In equation (11), represent Rate of change of concentration over time; represent At ℃, The baseline nitration rate constant; Represents the temperature coefficient; Represents actual water temperature Real-time concentration; The negative sign represents Concentration decreases over time; (12); In equation (12), represent Rate of change of concentration over time; represent The baseline nitration rate constant at ℃; Represents the temperature coefficient; Represents actual water temperature Real-time concentration; represent The baseline denitrification rate constant at ℃; Represents actual water temperature Real-time concentration; (13); In equation (13), represent Rate of change of concentration over time; represent At ℃, The baseline deposition rate constant; Represents actual water temperature Real-time concentration; (14); In equation (14), represent Rate of change of concentration over time; represent At ℃, The baseline deposition rate constant; Represents actual water temperature Real-time concentration; S43. The conservative tracer is injected into the surface water dynamic model. The conservative tracer channels of various pollution sources are completely independent. The kinetic equation of the conservative tracer is shown in the following equation (15): (15); In equation (15), This represents the rate of change of the conservative tracer concentration over time. This represents the conservative tracer concentration at the actual water temperature. Represents the surface velocity vector; This represents the divergence operator.

[0011] The pollution source tracing analysis specifically includes the following steps: S51. Construct a pollution source fingerprint database, which includes pollutant fingerprint concentrations corresponding to various pollution source types. S52, Adopt a pollution source-centered approach... To influence the radius, a distance-weighted strategy is used to superimpose the fingerprint concentration of the pollution source onto the triangular grid within the influence area, specifically: First, calculate the distance weighting coefficients of the distance-weighted strategy. , see the following formula (16) for details: (16); In equation (16), This represents the Euclidean distance between the centroid of the i-th triangular grid and the center of the pollution source; The radius of influence corresponding to the area of ​​influence; represent No more than the radius of influence For a triangular mesh, equation (16) holds true; Then, the fingerprint concentration of the pollution source is superimposed onto the triangular grid within the affected area to update the pollutant concentration at the triangular grid, as shown in the following formula (17): (17); In equation (17), This represents the pollutant concentration at the i-th triangular grid. The fingerprint concentration represents the concentration of each pollutant at the pollution source. The assignment operator is used for assignment updates. Equation (17) represents that the pollutant concentration of the i-th triangular grid at the current time step is calculated by adding the pollutant concentration of the i-th triangular grid at the previous time step to the fingerprint concentration of each pollutant at the pollution source multiplied by the distance weighting coefficient. Since the conservative tracer channels of each pollution source are completely independent, the conservative tracer concentration corresponding to the pollution source at the i-th triangular grid is calculated and updated simultaneously. The specific update of the conservative tracer concentration is shown in the following formula (18): (18); In equation (18), The conservative tracer concentration at the i-th triangular grid is represented by Equation (18). The conservative tracer concentration at the current time step of the i-th triangular grid is updated by adding the conservative tracer concentration at the previous time step of the i-th triangular grid to the distance weighting coefficient multiplied by 1000. 1000 is the initial concentration of the conservative tracer, which is dimensionless and is the baseline value for conservative tracer deployment. S53. The multi-source water quality tracing model calculates the Pearson correlation coefficient between the concentration of each pollutant and the concentration of the conservative tracer at each triangular grid within the affected area based on the time-series data of the monitoring points. , see the following formula (19) for details: (19); In equation (19), The pollutant concentrations at different time steps are obtained by updating the values ​​using equation (17); This represents the average concentration of pollutants at all time steps within the monitoring period. The conservative tracer concentrations at different time steps are obtained by updating the values ​​using equation (18); This represents the average conservative tracer concentration at all time steps within the monitoring period; At a certain triangular grid location, the Pearson correlation coefficient calculated from the conservative tracer concentration and pollutant concentration of a pollution source is... A value close to 1 indicates a high degree of synchronization between the changes in the conservative tracer and the pollutant, suggesting that a certain pollution source is the main source of pollution at this triangular grid location; when the Pearson correlation coefficient is close to 1... A value close to 0 or negative indicates that the conservative tracer and pollutant do not change synchronously, thus excluding a certain pollution source as the main pollution source at this triangular grid. S54. After determining the main pollution source at a certain triangular grid within the affected area using the multi-source water quality tracing model, a linear regression is performed on the scatter plot data of the conservative tracer and pollutants to establish a quantitative tracing relationship, as shown in the following formula (20): (20); In equation (20), The pollutant concentrations at different time steps are obtained by updating the values ​​using equation (17). The conservative tracer concentrations at different time steps are obtained by updating the values ​​using equation (18); Represents the regression coefficient; The intercept represents the background concentration of pollutants unaffected by the pollution source. Represents the coefficient of determination. The closer the value is to 1, the stronger the relationship between pollutants and conservative tracers, and the more accurate the prediction results of the multi-source water quality tracing model. The sum of squared errors between the predicted pollutant concentration and the actual pollutant concentration predicted by formula (20); This represents the sum of squares of the fluctuations in the actual pollutant concentration relative to its mean. The linear relationship between different pollutants and conservative tracers is calculated by equation (20). By substituting the concentration of the conservative tracer at the current time step into the linear relationship between different pollutants, the concentration of different pollutants at the current time step can be obtained. That is, the pollution load of the pollution source corresponding to the triangular grid can be calculated by the signal intensity of the conservative tracer.

[0012] It includes three operating modes: simulation mode, analysis mode, and combined mode. The simulation mode is performed in the MPI parallel environment, which divides the computational domain of the surface hydrodynamic model ANUGA into multiple subdomains. Different processes calculate the hydrodynamic state, water quality transport state and conservative tracer state in each subdomain, and exchange data at the subdomain boundaries. This completes the acquisition of multi-source data, bidirectional dynamic coupling of surface and pipeline networks, water quality transport analysis, pollution source injection and conservative tracer migration process, and outputs the spatial distribution results and time series data of pollutants and conservative tracers in the study area. The analysis mode is based on the time series data of pollutants and conservative tracers in the study area, calculates Pearson correlation coefficient and quantitative source tracing relationship, completes pollution source tracing correlation analysis, automatically generates a visual source tracing analysis report, and rapidly iterates and adjusts regression coefficients based on the same batch of simulation data. The combined mode is a combination of simulation mode and analysis mode. First, the simulation mode generates time-series data, and then the analysis mode completes the pollution source tracing correlation analysis and visualization output.

[0013] The MPI parallel speedup of the aforementioned MPI parallel environment satisfies Amdahl's law, and the MPI parallel speedup... It satisfies the following equation (21): (twenty one); In equation (21), This represents the proportion of the serial portion. This represents the number of parallel processes.

[0014] Advantages of this invention: (1) The present invention adopts bidirectional dynamic coupling of ANUGA and SWMM. In each coupling time step, the top support effect of surface water depth on drainage network and the feedback of network overflow on surface water volume are updated synchronously, which significantly improves the accuracy of urban rainstorm flood simulation and solves the simulation distortion problem caused by model decoupling.

[0015] (2) This invention performs water quality transport analysis on the computational domain of the surface hydrodynamic model, shares the grid data structure and time stepping mechanism with the surface hydrodynamic model, and does not require additional spatial interpolation or format conversion. It achieves strict synchronization of water quantity and water quality, reduces computational complexity and data transmission error, and fills the gap in the existing coupled framework that lacks water quality dimension simulation.

[0016] (3) The present invention adopts a strategy of combining a pollutant fingerprint database with multiple independent conservative tracer channels. Each independent pollution source is configured with a dedicated conservative tracer channel. Through Pearson correlation analysis, the dominant pollution source can be quantitatively identified and the contribution rate estimated, thus solving the problem of source tracing and unmixing in the case of multiple superimposed sources.

[0017] (4) Based on MPI large-scale parallel computing, the present invention distributes the ANUGA computing domain to multiple parallel processes for independent computing, making full use of multi-core computing resources, enabling high-resolution urban rainstorm and flood simulation to be completed within an engineering-acceptable time scale, and providing computing support for the transformation of real-time emergency response systems.

[0018] (5) This invention completes the entire source tracing analysis from data acquisition to a visualized source tracing analysis report, and the source tracing analysis stage can be executed independently without relying on the hydrodynamic simulation environment, thus reducing the threshold for engineering applications.

[0019] In summary, this invention has achieved substantial improvements in five dimensions: hydrodynamic coupling accuracy, water quality integration, source tracing quantification, computational efficiency, and engineering practicality. Attached Figure Description

[0020] Figure 1 This is a flowchart of the present invention.

[0021] Figure 2 This is a time-series graph of the pollutant concentrations of the present invention.

[0022] Figure 3 This is a time-series diagram of the conservative tracer concentration of the present invention.

[0023] Figure 4 This is a bar graph showing the Pearson correlation coefficients calculated from the conservative tracer concentration and various pollutant concentrations of a pollution source according to the present invention.

[0024] Figure 5 This is a scatter plot showing the correlation between the conservative tracer of this invention and the pollutant COD.

[0025] Figure 6 Time-series diagrams of water depth at each monitoring point in this embodiment of the invention.

[0026] Figure 7 This is a bar graph showing the Pearson correlation coefficients between the COD concentration at each monitoring point and the conservative tracer concentrations corresponding to the two pollution sources in this embodiment of the invention.

[0027] Figure 8 This is a time series diagram of the conservative tracer concentrations corresponding to the two pollution sources at MP3 in an embodiment of the present invention.

[0028] Figure 9This is a matrix heatmap showing the Pearson correlation coefficients between the COD concentration at each monitoring point and the concentrations of the conservative tracers corresponding to the two pollution sources in this embodiment of the invention. Detailed Implementation

[0029] 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. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0030] See Figure 1 A method for tracing the sources of multi-source pollution from urban stormwater flooding based on bidirectional dynamic coupling, specifically including the following steps: (1) Acquisition of multi-source data: Collect DEM data and drainage network node data of the area to be studied, and construct rainfall driving sequences; DEM data is in GeoTIFF format to directly obtain pixel elevation and spatial reference information; The drainage network node data is obtained by transforming the original coordinates of the drainage network nodes from the WGS-84 geographic coordinate system to the projected coordinate system corresponding to the area under study. The transformation process is shown in the following formula (1): (1); In equation (1), Represents longitude in the WGS-84 geographic coordinate system. Represents latitude in the WGS-84 geographic coordinate system. This represents the X-axis coordinate value in the UTM projected coordinate system. This represents the Y-axis coordinate value in the UTM projected coordinate system. This represents the transformation function from latitude and longitude to UTM projected X-axis coordinates. This represents the transformation function from latitude and longitude to the UTM projected Y-axis coordinates; The construction of the rainfall driving sequence first involves collecting rainfall observation data from meteorological stations or simulating future precipitation data using atmospheric numerical models. Then, the rainfall observation data or future precipitation data are arranged in time sequence and the units are converted. The specific unit conversion formula is shown in the following formula (2): (2); In equation (2), Represents observed rainfall data or future precipitation data, in units of... ; Represents rainfall observation data or future precipitation data after unit conversion, in units of ; When the duration of rainfall observation data or future precipitation data is shorter than the total simulation duration, the rainfall observation data or future precipitation data for the simulation time that exceeds the duration will be set to zero to ensure the integrity of the rainfall driving sequence. (2) Initialize the surface hydrodynamic model: Construct and initialize the computational domain of the two-dimensional shallow water equation solver ANUGA, which includes the following steps: S21. Within the spatial coverage of the DEM data, the finite volume method is used to generate an unstructured triangular computational grid in the area to be studied, and the computational domain of the two-dimensional shallow water equation solver ANUGA is constructed. The maximum triangle area in the unstructured triangular computational grid is constrained by the following formula (3). (3); In equation (3), This represents the basic resolution, measured in meters (m). Equation (3) represents the area of ​​the largest triangle; Equation (3) represents that the area of ​​a single triangle grid does not exceed the area of ​​the square corresponding to the edge of the basic resolution. S22. Set the computational domain parameters, specifically by using the flow algorithm to calculate the water flow wave velocity and constructing the CFL stability conditions, as shown in the following formula (4): (4); In equation (4), Represents the number of CFLs; The water flow wave velocity calculated by the representative flow algorithm; The time step representing the surface hydrodynamic model; This represents the size of the triangular grid, i.e., the basic resolution. Equation (4) represents the stability of the numerical calculation of the surface hydrodynamic model by controlling the ratio of water flow wave velocity, time step and triangular mesh size; S23. Assign the DEM data to the centroids of all triangular meshes using spatial interpolation, and set the smoothing coefficient. To suppress interpolation noise, the initial water surface elevation of the surface hydrodynamic model is calculated by the following equation (5): (5); In equation (5), Represents the initial water surface elevation; This represents the bed elevation after spatial interpolation of the DEM data; Represents the initial water depth under thin water film conditions; (3) Construct a two-way dynamic coupling model of the surface and the pipeline network, which specifically includes the following steps: S31. In each coupling time step, the surface water depth at all SWMM node coordinates is calculated in batches using spatial interpolation, and the surface water depth at all SWMM node coordinates is assigned to the corresponding SWMM node. The formula for calculating the surface water depth at the SWMM node coordinates is shown in the following formula (6): (6); In equation (6), Represents SWMM node coordinates The surface water depth at that location is the node in the SWMM stormwater and flood management model. represent The water surface elevation at the location is calculated using spatial interpolation based on the water surface elevation output from the surface hydrodynamic model. represent The surface elevation at that location; This represents taking the maximum value to ensure that the surface water depth is not negative. S32. For the SWMM nodes that overflow in the storm flood management model, calculate the total overflow volume according to the coupling time step. With each SWMM node that overflows as the center and a circular area with a set radius (default 2 m) as the overflow influence range, uniformly superimpose the total overflow volume onto the water surface elevation of the corresponding triangular grid within the overflow influence range. This transmits the location of the SWMM node that overflows and the total overflow volume to the surface hydrodynamic model, thus constructing a two-way dynamic coupling model between the surface and the pipe network. Total overflow volume The calculation formula is shown in the following formula (7): (7); In equation (7), This represents the total overflow volume generated by a single SWMM node overflow within the coupling time step; Represents the instantaneous overflow of a single SWMM node (m³ / s); Represents the coupling time step; Increment in water surface elevation caused by overflow in surface hydrodynamic models The calculation formula is shown in the following formula (8): (8); In equation (8), The area representing the extent of the overflow's influence, i.e., the area of ​​a circular region with a set radius; Steps S33, S31, and S32 are executed cyclically at each coupling time step. That is, at each coupling time step, the total overflow volume of the storm flood management model SWMM is transferred to the surface hydrodynamic model, and then the surface hydrodynamic model feeds back the water surface elevation increment to the storm flood management model SWMM, forming a two-way dynamic coupling. (4) Conduct water transport analysis: Initialize the concentrations of various pollutants in the study area on the surface hydrodynamic model, and use the water multi-component transport module to dynamically update the pollutant concentrations. This includes the following steps: S41, for each pollutant Create a concentration array with the same length as the total number of triangular meshes in the surface hydrodynamic model. , Represents chemical oxygen demand. Represents ammonium ions. Represents nitrate ions. Represents total phosphorus, Represents suspended matter; wherein, the initial concentration of each pollutant is determined by the background concentration of the pollutant, as shown in the following formula (9): (9); In equation (9), This represents the pollutants in the i-th triangular grid. The initial concentration at the initial time t=0; Represents the background concentration of pollutants; This represents the total number of triangular meshes in the surface hydrodynamic model; S42. Construct the kinetic equations for each pollutant, as shown in equations (10) to (14) below: (10); In equation (10), represent Rate of change of concentration over time; ,represent At ℃, The first-order decay rate constant of the concentration; , representing the temperature coefficient Represents the actual water temperature; Represents actual water temperature Real-time concentration; The negative sign represents Concentration decreases over time; (11); In equation (11), represent Rate of change of concentration over time; ,represent At ℃, The baseline nitration rate constant; ; Represents actual water temperature Real-time concentration; The negative sign represents Concentration decreases over time; (12); In equation (12), represent Rate of change of concentration over time; ; Represents the temperature coefficient; Represents actual water temperature Real-time concentration; ,represent The baseline denitrification rate constant at ℃; Represents actual water temperature Real-time concentration; (13); In equation (13), represent Rate of change of concentration over time; ,represent At ℃, The baseline deposition rate constant; Represents actual water temperature Real-time concentration; (14); In equation (14), represent Rate of change of concentration over time; ,represent At ℃, The baseline deposition rate constant; Represents actual water temperature Real-time concentration; S43. The conservative tracer is injected into the surface water dynamic model. The conservative tracer does not participate in any chemical or biochemical reaction. Its governing equation only contains the convective transport term and satisfies the mass conservation law. The conservative tracer channels of various pollution sources are completely independent. The kinetic equation of the conservative tracer is shown in the following equation (15): (15); In equation (15), This represents the rate of change of the conservative tracer concentration over time. This represents the conservative tracer concentration at the actual water temperature. Represents the surface velocity vector; Represents the divergence operator; (5) Pollution source tracing analysis: Construct a pollution source fingerprint database and a multi-pollution source water quality tracing model based on conservative tracers to perform qualitative identification and quantitative tracing analysis of pollution sources. The specific steps include: S51. Construct a pollution source fingerprint database, which includes pollutant fingerprint concentrations corresponding to various pollution source types, as shown in Table 1 below. Table 1

[0031] S52, Adopt a pollution source-centered approach... To influence the radius, a distance-weighted strategy is used to superimpose the fingerprint concentration of the pollution source onto the triangular grid within the influence area, specifically: First, calculate the distance weighting coefficients of the distance-weighted strategy. , see the following formula (16) for details: (16); In equation (16), This represents the Euclidean distance between the centroid of the i-th triangular grid and the center of the pollution source; The radius of influence corresponding to the area of ​​influence; represent No more than the radius of influence For a triangular mesh, equation (16) holds true; Then, the fingerprint concentration of the pollution source is superimposed onto the triangular grid within the affected area to update the pollutant concentration at the triangular grid, as shown in the following formula (17): (17); In equation (17), This represents the pollutant concentration at the i-th triangular grid. The fingerprint concentration represents the concentration of each pollutant at the pollution source. The assignment operator is used for assignment updates. Equation (17) represents that the pollutant concentration of the i-th triangular grid at the current time step is calculated by adding the pollutant concentration of the i-th triangular grid at the previous time step to the fingerprint concentration of each pollutant at the pollution source multiplied by the distance weighting coefficient. Since the conservative tracer channels of each pollution source are completely independent, the conservative tracer concentration corresponding to the pollution source at the i-th triangular grid is calculated and updated simultaneously. The specific update of the conservative tracer concentration is shown in the following formula (18): (18); In equation (18), The conservative tracer concentration at the i-th triangular grid is represented by Equation (18). The conservative tracer concentration at the current time step of the i-th triangular grid is updated by adding the conservative tracer concentration at the previous time step of the i-th triangular grid to the distance weighting coefficient multiplied by 1000. 1000 is the initial concentration of the conservative tracer, which is dimensionless and is the baseline value for conservative tracer deployment. S53. The multi-source water quality tracing model calculates the Pearson correlation coefficient between the concentration of each pollutant and the concentration of the conservative tracer at each triangular grid within the affected area based on the time-series data of the monitoring points. , see the following formula (19) for details: (19); In equation (19), The pollutant concentrations at different time steps are obtained by updating the values ​​using equation (17); This represents the average concentration of pollutants at all time steps within the monitoring period. The conservative tracer concentrations at different time steps are obtained by updating the values ​​using equation (18); This represents the average conservative tracer concentration at all time steps within the monitoring period; pollutant COD and The time series plot of concentration is shown in Figure 2 The time series plot of conservative tracer concentration is shown in [reference needed]. Figure 3 The bar graph showing the Pearson correlation coefficients calculated from the conservative tracer concentration and the concentrations of various pollutants for a certain pollution source is shown below. Figure 4 ; At a certain triangular grid location, the Pearson correlation coefficient calculated from the conservative tracer concentration and pollutant concentration of a pollution source is... A value close to 1 indicates a high degree of synchronization between the changes in the conservative tracer and the pollutant, suggesting that a certain pollution source is the main source of pollution at this triangular grid location; when the Pearson correlation coefficient is close to 1... A value close to 0 or negative indicates that the conservative tracer and pollutant do not change synchronously, thus excluding a certain pollution source as the main pollution source at this triangular grid. S54. After identifying the main pollution sources at a triangular grid location within the affected area using the multi-source water quality tracing model, linear regression is performed on the scatter plot data of conservative tracers and pollutants (see the scatter plot of the correlation between conservative tracers and pollutant COD). Figure 5 Establish a quantitative traceability relationship, as shown in the following formula (20): (20); In equation (20), The pollutant concentrations at different time steps are obtained by updating the values ​​using equation (17). The conservative tracer concentrations at different time steps are obtained by updating the values ​​using equation (18); Represents the regression coefficient; The intercept represents the background concentration of pollutants unaffected by the pollution source. Represents the coefficient of determination. The closer the value is to 1, the stronger the relationship between pollutants and conservative tracers, and the more accurate the prediction results of the multi-source water quality tracing model. The sum of squared errors between the predicted pollutant concentration and the actual pollutant concentration predicted by formula (20); This represents the sum of squares of the fluctuations in the actual pollutant concentration relative to its mean. The linear relationship between different pollutants and conservative tracers is calculated by equation (20). By substituting the concentration of the conservative tracer at the current time step into the linear relationship between different pollutants, the concentration of different pollutants at the current time step can be obtained. That is, the pollution load of the pollution source corresponding to the triangular grid can be calculated by the signal intensity of the conservative tracer.

[0032] The validation metrics for multi-source water quality tracing models include mean absolute error (MAE) and root mean square error (RMSE), which are used to quantitatively assess the model's accuracy.

[0033] A method for tracing the source of multi-source pollution from urban rainstorms and flooding based on bidirectional dynamic coupling includes three operating modes: simulation mode, analysis mode, and joint mode. The simulation mode operates in an MPI parallel environment, dividing the computational domain of the ANUGA surface hydrodynamic model into multiple subdomains. Different processes calculate the hydrodynamic state, water transport state, and conservative tracer state within each subdomain, exchanging data at the subdomain boundaries. This process acquires multi-source data, dynamically couples surface and pipe networks, analyzes water transport, simulates pollutant injection and conservative tracer migration, and outputs spatial distribution results and time-series data of pollutants and conservative tracers within the study area. The MPI parallel speedup in this environment satisfies Amdahl's law. It satisfies the following equation (21): (twenty one); In equation (21), This represents the proportion of the serial portion. Represents the number of parallel processes; The analysis mode is based on the time series data of pollutants and conservative tracers in the study area, calculates the Pearson correlation coefficient and quantitative source tracing relationship, completes the pollution source tracing correlation analysis, automatically generates a visual source tracing analysis report, and rapidly iterates and adjusts the regression coefficients based on the same batch of simulation data; The combined mode is a combination of simulation mode and analysis mode. First, the simulation mode generates time-series data, and then the analysis mode completes the pollution source tracing correlation analysis and visualization output.

[0034] Example (I) Acquisition of multi-source data:

[0035] This embodiment uses a drainage area in a central urban area as the study area. The study area is located in the UTM Zone 50N projected coordinate system (EPSG:32650) and includes a mixed underlying surface of residential areas, industrial areas, and urban green spaces, with an impermeability of approximately 62%. The digital elevation model (DEM) data is sourced from the Geospatial Data Cloud (https: / / www.gscloud.cn), in GeoTIFF format, with a spatial resolution of 30 m, a raster size of 555×474 pixels, an effective pixel ratio of 92.3%, and a terrain elevation range of 4.00–75.00 m, covering the entire study area.

[0036] The drainage network node data has been simplified, containing 425 drainage network nodes and several pipe segments. The SWMM input file (.inp format) is automatically generated by the system program. The original coordinates of the drainage network nodes have been batch converted from WGS-84 (EPSG:4326) to UTM Zone 50N (EPSG:32650) using a coordinate transformation tool and stored separately for subsequent construction of a two-way dynamic coupling model of the surface and the network.

[0037] The rainfall observation data for the rainfall-driven sequence were continuous measured rainfall time series (1 h time resolution, mm / h) from 06:00 on June 19, 2025 to 21:00 on June 21, 2025, with a total rainfall duration of 63 hours. The rainfall process exhibited a multi-peak pattern: the rainfall intensity was relatively weak in the early stage (approximately 5 mm accumulated in the first 12 hours), followed by a period of heavy rainfall. The rainfall intensity increased to 2.5 mm / h at t≈40 h, reached a peak of 5 mm / h at t≈48 h, and then reached another peak of 5 mm / h at t≈54 h. The cumulative rainfall during the heavy rainfall phase was approximately 42 mm, which was the main driving force for surface water accumulation and pipe network overflow.

[0038] (II) Parameter Configuration: The basic resolution is set using the DEM-scale polygon as the boundary of the ANUGA computational domain. The minimum triangle angle is 28°, generating an unstructured triangular computational mesh. The flow algorithm is set to DE0, the CFL number is 0.7, the initial water depth is 0.05 m, and the Manning roughness coefficient n = 0.1; all four boundaries are set as reflective boundaries.

[0039] Coupling time step The simulation output time step yieldstep = 900 s, water quality output interval 300 s, and MPI parallel process number N = 20. The pollutant types are set to COD, , Five substances were tested: COD, TP, TSS, etc., with background concentrations of: COD = 10.0 mg / L, etc. = 0.1 mg / L = 0.5 mg / L, TP = 0.05 mg / L, TSS = 5.0 mg / L. Reaction kinetic parameters: COD decay rate The baseline nitration rate constant The baseline denitrification rate constant Reference deposition rate constant Reference deposition rate constant Temperature coefficient .

[0040] (III) Pollution Sources and Monitoring Point Setup: Based on the DEM low-lying area analysis (excluding edge pixels of 10 pixels and pixels with elevation ≤ 0.1 m, the entire area was divided into 3×3 partitions and the lowest point of each partition was extracted), and combined with the topographic slope and drainage direction, two pipeline overflow pollution sources, source_1 and source_2, were set in the low-lying area in the middle of the study area (elevation 13 m). Three source tracing monitoring points, MP1, MP2, and MP3, were set up due south of the two sources (radius of overflow influence range = 30 m). All coordinates are UTM Zone 50N. Both pipeline overflow pollution sources are located in the central area of ​​the ANUGA computational domain, more than 4500 m from the nearest boundary.

[0041] Source_1 (domestic sewage): coordinates (519189, 3511327), triggered at t = 144000 s (40 h, at the start of heavy rainfall), discharge flow rate 0.08 m³ / s, duration 7200 s (2 h), total discharge volume 576 m³; pollutant fingerprint concentration: COD = 300 mg / L. = 25 mg / L, TP = 4.0 mg / L, TSS = 150 mg / L; equipped with an independent conserved tracer channel Tracer_S1 (initial concentration 1000).

[0042] Source_2 (industrial wastewater): Coordinates (520389, 3511327), located 1200 m east of source_1, triggered at t = 172800 s (48 h, peak rainfall 5 mm / h), discharge flow rate 0.04 m³ / s, duration 3600 s (1 h), total discharge volume 144 m³; pollutant fingerprint concentration: COD = 500 mg / L. = 10 mg / L, TP = 1.0 mg / L, TSS = 100 mg / L; equipped with an independent conserved tracer channel Tracer_S2 (initial concentration 1000).

[0043] The monitoring points are designed with a symmetrical layout, making full use of the low-lying catchment area to the south of the two water sources: MP1 (519189, 3510527) is located 800 m south of source_1 (estimated propagation time 1.9 h) and approximately 1442 m (3.3 h) from source_2. MP2 (519789, 3510727) is located due south between the two sources, equidistant from source_1 and source_2 (approximately 849 m and 2.0 h respectively), and is designed as a mixed signal verification point. MP3 (520389, 3510527) is located 800 m (1.9 h) south of source_2 and about 1442 m (3.3 h) from source_1.

[0044] The monitoring point positions of MP1 and MP3 are completely symmetrical along the midpoint of the line connecting source_1 and source_2.

[0045] (IV) Simulation Operation and Result Analysis: After launching a full simulation in the MPI parallel environment and completing a 63-hour full-domain two-way dynamic coupling and water transport simulation, the water depth time series of three monitoring points was extracted from the simulation output file, and a visual source tracing analysis report was automatically generated, including water depth time series plots for each monitoring point (see...). Figure 6 Bar chart of Pearson correlation coefficient (see) Figure 7 Time series plot of conservative tracer concentration at MP3 (see) Figure 8 A matrix plot of the correlation coefficient between Pearson and Pearson (see...) Figure 9 ).

[0046] See Figure 6The water depth time series at the three monitoring points were all derived from the simulation output of the surface hydrodynamic model, verifying the end-to-end connectivity of the water depth extraction process. The water depth at MP1 initially decreased to approximately 0.002 m due to the near-zero drainage slope at the beginning of rainfall (t ≈ 5–25 h), then rose again to 0.035 m following heavy rainfall at t = 40 h, exhibiting a dynamic response of initial decline followed by accumulation, reflecting the influence of micro-topographic undulations in the actual terrain. The water depth at MP2 remained between 0.03 and 0.04 m, indicating moderate topographic drainage conditions. MP3, located in a low-lying catchment area, saw its water depth continuously accumulate from an initial 0.06 m with rainfall, reaching a peak of 0.138 m at t = 60.3 h, showing a continuous upward trend throughout, reflecting the runoff and water accumulation characteristics under actual topographic conditions.

[0047] The pollutant concentration time series was generated by overlaying a Gaussian diffusion model with a pollutant fingerprint database, with the pollution source triggering time corresponding to the period of heavy rainfall. The peak COD concentration at MP1 was 94.7 mg / L, occurring at t = 42.7 h (approximately 2.7 h after source_1 triggering, which is basically consistent with the 1.9 h estimate corresponding to the 800 m propagation distance); the peak COD concentration at MP2 was 91.0 mg / L, occurring at t = 43.0 h; and the peak COD concentration at MP3 was 94.2 mg / L, occurring at t = 50.2 h (approximately 2.2 h after source_2 triggering, which is consistent with the 1.9 h propagation time for 800 m). The peak COD concentrations at the three monitoring points were similar in magnitude (91–95 mg / L), and the differences in the peak times reflected the propagation time from each point to different dominant pollution sources.

[0048] See Figures 7-9 Pearson correlation analysis was performed on the time series of the conservative tracer at each monitoring point. The results showed a significant spatial mirror-symmetric structure, which was completely consistent with the theoretical expectation: the Pearson correlation coefficient of MP1 with the conservative tracer Tracer_S1 and COD was... Linear regression Pearson correlation coefficient of MP1 with the conservative tracer Tracer_S2 and COD Linear regression Based on the Pearson correlation coefficient Conclusion: The COD pollution at MP1 mainly originates from source_1 (domestic sewage); MP3 correlation coefficient between the conservative tracer Tracer_S1 and COD. Linear regression Pearson correlation coefficient of MP3 with the conservative tracer Tracer_S2 and COD Linear regression Based on the Pearson correlation coefficient The conclusion is that the COD pollution at MP3 mainly comes from source_2 (industrial wastewater).

[0049] MP2 showed a Pearson correlation coefficient between the conservative tracer Tracer_S1 and COD. Linear regression Pearson correlation coefficient of MP2 with the conservative tracer Tracer_S2 and COD Linear regression Source_1 is slightly dominant, which is consistent with the setting that the emission of source_1 (0.08 m³ / s) is greater than that of source_2 (0.04 m³ / s).

[0050] Figure 9 The matrix heatmaps exhibit a clear diagonal high correlation structure, which intuitively verifies the spatial differentiation ability of multiple independent conservative tracer channels for multiple pollution sources.

[0051] Figure 8 The time series plot of the conservative tracer concentration at MP3 is shown: the peak 123 of the conservative tracer Tracer_S1 appears at t = 44.0 h, and the peak 321 of the conservative tracer Tracer_S2 appears at t = 50.3 h. The time difference between the two peaks is about 6.3 h, which is clearly distinguishable, verifying the ability of the present invention to identify two independent pollution source signals that arrive at the same monitoring point.

[0052] (V) Performance Analysis Summary: This embodiment, based on real DEM (555×474 pixels, 30 m resolution) data of a certain urban area and a 63-hour measured Meiyu (plum rain) precipitation process, completed a two-way dynamic coupling and multi-component water quality transport simulation in a parallel environment of the MPI 20 process. The implementation results verified the following key characteristics: (1) Real hydrodynamic background: The water depth time series of the three monitoring points are all derived from the simulation output of the ANUGA computational domain. At MP3, water continued to accumulate to 0.138 m, while at MP1, a real micro-topographic response was shown to first recede and then accumulate, which verified the effective coupling between the water depth extraction process and the surface hydrodynamic field.

[0053] (2) Spatial source differentiation capability: The Pearson correlation coefficient difference between the two conservative tracers at MP1 and MP3 is Δr = 0.628 and Δr = 0.391, respectively, which is significantly higher than the random threshold. The two independent conservative tracer channels successfully achieved quantitative differentiation of two pollution sources with a distance of 1200 m. The source attribution conclusion is completely consistent with the spatial location of the pollution source.

[0054] (3) Time-series distinguishability: The peak time difference between the two conservative tracers at MP3 is 6.3 h, which is clearly distinguishable in the 63 h simulated time series, verifying the ability to identify signals from multiple pollution sources that arrive one after another.

[0055] (4) Full-chain process verification: From DEM data reading, coordinate transformation, bidirectional dynamic coupling simulation, water depth extraction, test data generation to the automatic output of views in the visualization traceability analysis report, the entire chain runs smoothly within the acceptable calculation time of the project, providing a feasible technical verification for the transformation of the method of this invention into an actual engineering emergency response system.

[0056] Although embodiments of the invention have been shown and described, it will be understood by those skilled in the art that various changes, modifications, substitutions and alterations can be made to these embodiments without departing from the principles and spirit of the invention, the scope of which is defined by the appended claims and their equivalents.

Claims

1. A method for tracing the sources of multi-source pollution from urban stormwater flooding based on bidirectional dynamic coupling, characterized in that: Specifically, it includes the following steps: (1) Acquisition of multi-source data: Collect DEM data and drainage network node data of the area to be studied, and construct rainfall driving sequences; (2) Initialize the surface hydrodynamic model: Construct and initialize the computational domain of the two-dimensional shallow water equation solver ANUGA; (3) Construct a two-way dynamic coupling model between the surface and the pipeline network: The overflow of the SWMM node is transmitted to the surface hydrodynamic model, and the surface hydrodynamic model feeds back the water surface elevation increment to the storm flood management model SWMM, forming a two-way dynamic coupling, and thus constructing a two-way dynamic coupling model between the surface and the pipeline network. (4) Conduct water transport analysis: Initialize the concentration of various pollutants in the study area on the surface hydrodynamic model, and use the water multi-component transport module to dynamically update the pollutant concentration; (5) Pollution source tracing analysis: Construct a pollution source fingerprint database and a multi-pollution source water quality tracing model based on conservative tracers to conduct qualitative identification and quantitative tracing analysis of pollution sources.

2. The method for tracing the source of multi-source pollution from urban stormwater flooding based on bidirectional dynamic coupling as described in claim 1, characterized in that: The DEM data is in GeoTIFF format; The drainage network node data is obtained by transforming the original coordinates of the drainage network nodes from the WGS-84 geographic coordinate system to the projected coordinate system corresponding to the area under study. The transformation process is shown in the following formula (1): （1）； In equation (1), Represents longitude in the WGS-84 geographic coordinate system. Represents latitude in the WGS-84 geographic coordinate system. This represents the X-axis coordinate value in the UTM projected coordinate system. This represents the Y-axis coordinate value in the UTM projected coordinate system. This represents the transformation function from latitude and longitude to UTM projected X-axis coordinates. This represents the transformation function from latitude and longitude to the UTM projected Y-axis coordinates; The construction of the rainfall driving sequence first involves collecting rainfall observation data from meteorological stations or simulating future precipitation data using atmospheric numerical models. Then, the rainfall observation data or future precipitation data are arranged in time sequence and the units are converted. The specific unit conversion formula is shown in the following formula (2): （2）； In equation (2), Represents observed rainfall data or future precipitation data, in units of... ; Represents rainfall observation data or future precipitation data after unit conversion, in units of ; When the duration of rainfall observation data or future precipitation data is shorter than the total simulation duration, the rainfall observation data or future precipitation data for the excess simulation time will be set to zero to ensure the integrity of the rainfall driving sequence.

3. The method for tracing the source of multi-source pollution from urban stormwater flooding based on bidirectional dynamic coupling as described in claim 1, characterized in that: The initialization of the surface hydrodynamic model specifically includes the following steps: S21. Within the spatial coverage of the DEM data, the finite volume method is used to generate an unstructured triangular computational grid in the area to be studied, and the computational domain of the two-dimensional shallow water equation solver ANUGA is constructed. The maximum triangle area in the unstructured triangular computational grid is constrained by the following formula (3). （3）； In equation (3), This represents the basic resolution, measured in meters (m). Equation (3) represents the area of ​​the largest triangle; Equation (3) represents that the area of ​​a single triangle grid does not exceed the area of ​​the square corresponding to the edge of the basic resolution. S22. Set the computational domain parameters, specifically by using the flow algorithm to calculate the water flow wave velocity and constructing the CFL stability conditions, as shown in the following formula (4): （4）； In equation (4), Represents the number of CFLs; The water flow wave velocity calculated by the representative flow algorithm; The time step representing the surface hydrodynamic model; This represents the size of the triangular grid, i.e., the basic resolution. Equation (4) represents the stability of the numerical calculation of the surface hydrodynamic model by controlling the ratio of water flow wave velocity, time step and triangular mesh size; S23. Assign the DEM data to the centroid of all triangular meshes using spatial interpolation, and set a smoothing coefficient to suppress interpolation noise. The initial water surface elevation of the surface hydrodynamic model is calculated by the following formula (5): （5）； In equation (5), Represents the initial water surface elevation; This represents the bed elevation after spatial interpolation of the DEM data; This represents the initial water depth under thin water film conditions.

4. The method for tracing the source of multi-source pollution from urban stormwater flooding based on bidirectional dynamic coupling as described in claim 3, characterized in that: The construction of the bidirectional dynamic coupling model of the land surface and the pipeline network specifically includes the following steps: S31. In each coupling time step, the surface water depth at all SWMM node coordinates is calculated in batches using spatial interpolation, and the surface water depth at all SWMM node coordinates is assigned to the corresponding SWMM node. The formula for calculating the surface water depth at the SWMM node coordinates is shown in the following formula (6): （6）； In equation (6), Represents SWMM node coordinates The surface water depth at that location is the node in the SWMM stormwater and flood management model. represent The water surface elevation at the location is calculated using spatial interpolation based on the water surface elevation output from the surface hydrodynamic model. represent The surface elevation at that location; This represents taking the maximum value to ensure that the surface water depth is not negative. S32. For the SWMM nodes that overflow in the storm flood management model, calculate the total overflow volume according to the coupling time step. With each SWMM node that overflows as the center and a circular area with a set radius as the overflow influence range, uniformly superimpose the total overflow volume onto the water surface elevation of the corresponding triangular grid within the overflow influence range. This transmits the location of the SWMM node that overflows and the total overflow volume to the surface hydrodynamic model, thus constructing a two-way dynamic coupling model between the surface and the pipe network. The total overflow volume The calculation formula is shown in the following formula (7): （7）； In equation (7), This represents the total overflow volume generated by a single SWMM node overflow within the coupling time step; Represents the instantaneous overflow of a single SWMM node (m³ / s); Represents the coupling time step; The water surface elevation increment caused by overflow in the surface hydrodynamic model The calculation formula is shown in the following formula (8): （8）； In equation (8), The area representing the extent of the overflow's influence, i.e., the area of ​​a circular region with a set radius; Steps S33, S31, and S32 are executed cyclically at each coupling time step. That is, at each coupling time step, the total overflow volume of the stormwater flood management model (SWMM) is transferred to the surface hydrodynamic model, and then the surface hydrodynamic model feeds back the water surface elevation increment to the stormwater flood management model (SWMM), forming a two-way dynamic coupling.

5. The method for tracing the source of multi-source pollution from urban stormwater flooding based on bidirectional dynamic coupling as described in claim 3, characterized in that: The water transport analysis specifically includes the following steps: S41, for each pollutant Create a concentration array with the same length as the total number of triangular meshes in the surface hydrodynamic model. , Represents chemical oxygen demand. Represents ammonium ions. Represents nitrate ions. Represents total phosphorus, Represents suspended matter; wherein, the initial concentration of each pollutant is determined by the background concentration of the pollutant, as shown in the following formula (9): （9）； In equation (9), This represents the pollutants in the i-th triangular grid. The initial concentration at the initial time t=0; Represents the background concentration of pollutants; This represents the total number of triangular meshes in the surface hydrodynamic model; S42. Construct the kinetic equations for each pollutant, as shown in equations (10) to (14) below: （10）； In equation (10), represent Rate of change of concentration over time; represent At ℃, The first-order decay rate constant of the concentration; Represents the temperature coefficient. Represents the actual water temperature; Represents actual water temperature Real-time concentration; The negative sign represents Concentration decreases over time; （11）； In equation (11), represent Rate of change of concentration over time; represent At ℃, The baseline nitration rate constant; Represents the temperature coefficient; Represents actual water temperature Real-time concentration; The negative sign represents Concentration decreases over time; （12）； In equation (12), represent Rate of change of concentration over time; represent The baseline nitration rate constant at ℃; Represents the temperature coefficient; Represents actual water temperature Real-time concentration; represent The baseline denitrification rate constant at ℃; Represents actual water temperature Real-time concentration; （13）； In equation (13), represent Rate of change of concentration over time; represent At ℃, The baseline deposition rate constant; Represents actual water temperature Real-time concentration; （14）； In equation (14), represent Rate of change of concentration over time; represent At ℃, The baseline deposition rate constant; Represents actual water temperature Real-time concentration; S43. The conservative tracer is injected into the surface water dynamic model. The conservative tracer channels of various pollution sources are completely independent. The kinetic equation of the conservative tracer is shown in the following equation (15): （15）； In equation (15), This represents the rate of change of the conservative tracer concentration over time. This represents the conservative tracer concentration at the actual water temperature. Represents the surface velocity vector; This represents the divergence operator.

6. The method for tracing the source of multi-source pollution from urban stormwater flooding based on bidirectional dynamic coupling as described in claim 5, characterized in that: The pollution source tracing analysis specifically includes the following steps: S51. Construct a pollution source fingerprint database, which includes pollutant fingerprint concentrations corresponding to various pollution source types. S52, Adopt a pollution source-centered approach... To influence the radius, a distance-weighted strategy is used to superimpose the fingerprint concentration of the pollution source onto the triangular grid within the influence area, specifically: First, calculate the distance weighting coefficients of the distance-weighted strategy. , see the following formula (16) for details: （16）； In equation (16), This represents the Euclidean distance between the centroid of the i-th triangular grid and the center of the pollution source; The radius of influence corresponding to the area of ​​influence; represent No more than the radius of influence For a triangular mesh, equation (16) holds true; Then, the fingerprint concentration of the pollution source is superimposed onto the triangular grid within the affected area to update the pollutant concentration at the triangular grid, as shown in the following formula (17): （17）； In equation (17), This represents the pollutant concentration at the i-th triangular grid. The fingerprint concentration represents the concentration of each pollutant at the pollution source. The assignment operator is used for assignment updates. Equation (17) represents that the pollutant concentration of the i-th triangular grid at the current time step is calculated by adding the pollutant concentration of the i-th triangular grid at the previous time step to the fingerprint concentration of each pollutant at the pollution source multiplied by the distance weighting coefficient. Since the conservative tracer channels of each pollution source are completely independent, the conservative tracer concentration corresponding to the pollution source at the i-th triangular grid is calculated and updated simultaneously. The specific update of the conservative tracer concentration is shown in the following formula (18): （18）； In equation (18), The conservative tracer concentration at the i-th triangular grid is represented by Equation (18). The conservative tracer concentration at the current time step of the i-th triangular grid is updated by adding the conservative tracer concentration at the previous time step of the i-th triangular grid to the distance weighting coefficient multiplied by 1000. 1000 is the initial concentration of the conservative tracer, which is dimensionless and is the baseline value for conservative tracer deployment. S53. The multi-source water quality tracing model calculates the Pearson correlation coefficient between the concentration of each pollutant and the concentration of the conservative tracer at each triangular grid within the affected area based on the time-series data of the monitoring points. , see the following formula (19) for details: （19）； In equation (19), The pollutant concentrations at different time steps are obtained by updating the values ​​using equation (17); This represents the average concentration of pollutants at all time steps within the monitoring period. The conservative tracer concentrations at different time steps are obtained by updating the values ​​using equation (18); This represents the average conservative tracer concentration at all time steps within the monitoring period; At a certain triangular grid location, the Pearson correlation coefficient calculated from the conservative tracer concentration and pollutant concentration of a pollution source is... A value close to 1 indicates a high degree of synchronization between the changes in the conservative tracer and the pollutant, suggesting that a certain pollution source is the main source of pollution at this triangular grid location; when the Pearson correlation coefficient is close to 1... A value close to 0 or negative indicates that the conservative tracer and pollutant do not change synchronously, thus excluding a certain pollution source as the main pollution source at this triangular grid. S54. After determining the main pollution source at a certain triangular grid within the affected area using the multi-source water quality tracing model, a linear regression is performed on the scatter plot data of the conservative tracer and pollutants to establish a quantitative tracing relationship, as shown in the following formula (20): （20）； In equation (20), The pollutant concentrations at different time steps are obtained by updating the values ​​using equation (17). The conservative tracer concentrations at different time steps are obtained by updating the values ​​using equation (18); Represents the regression coefficient; The intercept represents the background concentration of pollutants unaffected by the pollution source. Represents the coefficient of determination. The closer the value is to 1, the stronger the relationship between the pollutant and the conservative tracer, and the more accurate the pollutant concentration predicted by the multi-source water quality tracing model. The sum of squared errors between the predicted pollutant concentration and the actual pollutant concentration predicted by formula (20); This represents the sum of squares of the fluctuations in the actual pollutant concentration relative to its mean. The linear relationship between different pollutants and conservative tracers is calculated by equation (20). By substituting the concentration of the conservative tracer at the current time step into the linear relationship between different pollutants, the concentration of different pollutants at the current time step can be obtained. That is, the pollution load of the pollution source corresponding to the triangular grid can be calculated by the signal intensity of the conservative tracer.

7. The method for tracing the source of multi-source pollution from urban stormwater flooding based on bidirectional dynamic coupling as described in claim 6, characterized in that: It includes three operating modes: simulation mode, analysis mode, and combined mode. The simulation mode is performed in the MPI parallel environment, which divides the computational domain of the surface hydrodynamic model ANUGA into multiple subdomains. Different processes calculate the hydrodynamic state, water quality transport state and conservative tracer state in each subdomain, and exchange data at the subdomain boundaries. This completes the acquisition of multi-source data, bidirectional dynamic coupling of surface and pipeline networks, water quality transport analysis, pollution source injection and conservative tracer migration process, and outputs the spatial distribution results and time series data of pollutants and conservative tracers in the study area. The analysis mode is based on the time series data of pollutants and conservative tracers in the study area, calculates Pearson correlation coefficient and quantitative source tracing relationship, completes pollution source tracing correlation analysis, automatically generates a visual source tracing analysis report, and rapidly iterates and adjusts regression coefficients based on the same batch of simulation data. The combined mode is a combination of simulation mode and analysis mode. First, the simulation mode generates time-series data, and then the analysis mode completes the pollution source tracing correlation analysis and visualization output.

8. The method for tracing the source of multi-source pollution from urban stormwater flooding based on bidirectional dynamic coupling as described in claim 7, characterized in that: The MPI parallel speedup of the MPI parallel environment satisfies Amdahl's law, and the MPI parallel speedup... It satisfies the following equation (21): （21）； In equation (21), This represents the proportion of the serial portion. This represents the number of parallel processes.

Images