Urban waterlogging modeling method considering spatio-temporal heterogeneity

Abstract

This invention relates to a method for modeling urban flooding that takes into account spatiotemporal heterogeneity, comprising the following steps: constructing a spatiotemporally heterogeneous rainfall time series and multi-scale spatiotemporal grid; data preprocessing, dividing the study area into initial response units based on topographic features, land use types, and drainage system distribution; using a dynamic time warping method to measure the similarity of flooding response processes among response units, and performing time series clustering on the division results accordingly; using a cohesive hierarchical clustering method to further cluster and determine the optimal number of partitions, and introducing spatial adjacency constraints to ensure spatial connectivity; assigning differentiated hydrological parameters to different response units; and using multi-scale rainfall spatiotemporal grid data and response unit division results to drive the surface-pipeline coupling modeling of urban flooding. This invention improves the accuracy of urban flooding modeling by constructing a multi-scale rainfall grid and dividing urban flooding response units, thus expanding the existing modeling and simulation methods for urban flooding.

Description

Technical Field

[0001] This invention relates to the fields of geographic information technology and hydrology, and in particular to a method for modeling urban flooding that takes into account spatiotemporal heterogeneity. Background Technology

[0002] Against the backdrop of global climate change and rapid urbanization, extreme rainfall events are becoming more frequent and intense, exceeding the design capacity of traditional urban drainage systems. Simultaneously, the rapid expansion of hardened surfaces has led to the replacement of natural underlying surfaces with large amounts of impermeable surfaces, severely weakening the natural infiltration and retention capacity of rainwater. As a result, urban flooding disasters are showing a trend of increasing frequency, expanding disaster area, and intensified damage. As a typical "non-engineering" measure in drainage and flood control systems, urban flood modeling and simulation provide strong support for flood disaster prediction and early warning.

[0003] Traditional urban flood modeling typically does not consider spatiotemporal heterogeneity in its response unit division and assumes a homogeneous distribution of hydrological parameters in the study area. This ignores the significant spatiotemporal heterogeneity in urban flood processes and makes it difficult to fully reflect the impact of spatiotemporal heterogeneity characteristics such as local urban features and rainfall non-stationarity on urban flood processes. Consequently, the model's response unit division and parameter settings cannot accurately reflect the spatiotemporal dynamic differences in urban topography, underlying surface conditions, drainage systems, and human activities. This reduces the accuracy of the model in predicting the location, time, and intensity of urban floods and limits its adaptability to complex scenarios.

[0004] Therefore, this invention proposes an urban flooding modeling method that takes into account spatiotemporal heterogeneity. Under the condition of taking into account spatiotemporal heterogeneity, it comprehensively considers multiple factors such as rainfall distribution, topographic features, land use type, and drainage system distribution to construct a surface-pipeline coupled model, thereby achieving high-precision simulation of urban flooding processes and providing scientific basis and technical support for urban flood control and drainage decisions. Summary of the Invention

[0005] This invention addresses the shortcomings of existing urban flooding modeling methods by providing a method that considers spatiotemporal heterogeneity. The aim is to introduce spatiotemporal heterogeneity into urban flooding modeling, overcome the limitations of the homogeneous assumption in traditional urban flooding modeling, and enable the optimization and adjustment of model parameters according to environmental conditions. This improves the accuracy and spatiotemporal adaptability of urban flood simulation, providing more reliable scientific support for urban flood control and disaster reduction decisions.

[0006] This invention provides a method for modeling urban stormwater models based on urban flood response unit division, characterized by the following steps:

[0007] Step S1: Construct a multi-scale spatiotemporal grid of rainfall with spatiotemporal heterogeneity;

[0008] Step S2: Data preprocessing, dividing the study area into initial response units based on topographic features, land use types, and drainage system distribution;

[0009] Step S3: Use the dynamic time warping method to measure the similarity of the waterlogging response process among the response units, and perform time series clustering on the partitioning results accordingly;

[0010] Step S4: Further clustering is performed using agglomerative hierarchical clustering to determine the optimal number of partitions, and spatial adjacency constraints are introduced to ensure spatial connectivity;

[0011] Step S5: Assign differentiated hydrological parameters to different response units;

[0012] Step S6: Use multi-scale rainfall spatiotemporal grid data and response unit partitioning results to drive the surface-pipeline coupling model of urban flooding.

[0013] Further, S1 includes the following steps:

[0014] Step S11: Based on the time axis length and time resolution of the rainfall to be inverted, divide the time interval t0~t n Divide the time slices into n equal layers based on time steps. Each layer precisely corresponds to the spatial rainfall distribution over a specific time period. Each time slice satisfies the following:

[0015]

[0016] In the above formula, t0 represents the start time of the rainfall process, t n The time of the end of the rainfall event is indicated by t, where t represents the time step of the main time slice, and n represents the number of main time slices in the rainfall event. i This represents the start time corresponding to the i-th main time slice;

[0017] Step S12: When retrieving rainfall observation points, attach the target time slice [t] i， t i The constraint condition +t] is used to filter valid point rainfall data located in the modeling region Ω within this sub-time slice. Based on the filtering results, the main time slice is subdivided into k sub-time slices, specifically satisfying:

[0018] S i ={(x j ,y j ,R j )|(x j ,y j )∈Ω,t j ∈[t i ,t i +t]}

[0019] In the above formula, Ω represents the spatial extent of the modeling region, and xi ,y j Let R represent the coordinates of the j-th observation point. j Let t represent the rainfall observation value at the j-th observation point. j S represents the observation time recorded at the j-th observation point, k represents the number of sub-time slices within the current main time slice, and S represents the sub-time slice. i This represents the set of valid rainfall observation data within the i-th main time slice and located in the modeling region Ω;

[0020] S13: Spatial interpolation calculation is performed on the point-like rainfall data obtained by filtering time segments. The algorithm parameters are set according to the spatial scale and interpolation range, and the interpolation is performed to generate rainfall raster data within the smallest time segment. The rainfall distribution at time index i is summed in segments to accumulate the rainfall.

[0021] S14: Repeat the above process until all n time slices are calculated. Finally, stack them along the time axis to build a rainfall grid dataset that takes into account multiple time scales and multiple spatial resolutions, providing high-precision rainfall driving force for subsequent urban flood models.

[0022] Further, S2 includes the following steps:

[0023] S21: Using high-resolution DEM as the basic spatial data, the study area is divided into regular and uniform grid cells with clear spatial topological relationships between cells. The DEM is then subjected to terrain correction and spatial preprocessing to remove errors.

[0024] S22: The geographic location of each rainwater well node is expressed using a vector structure, which is regarded as the initial control node of spatial confluence. With each rainwater well as the core, the Thiessen polygon method is used to determine the area with the closest spatial distance, and it is initially classified as the response unit of the rainwater well.

[0025] S23: Spatial matching of land use types in the modeling area with grid cells to determine the impermeability of each cell, and the overall impermeability C of the response cell. m The calculation formula is as follows:

[0026]

[0027] In the above formula, C m Let A be the overall impermeability of unit m. m,l Let α be the area of ​​land use type l within unit m. l Let A be the impermeability coefficient corresponding to land use type l, where L represents the number of land use types included. m The total area of ​​unit m;

[0028] S24: Further, directional similarity, distance similarity, and land use similarity indices are introduced to establish a joint similarity matrix S. S is used to optimize the Thiessen polygon division results. Directional similarity is mainly determined by calculating the angle between the runoff directions of adjacent response units to determine spatial correlation. The calculation formula is as follows:

[0029]

[0030] In the above formula, θ is the angle between the line connecting the centroids of response unit m and its adjacent response unit k and its own runoff direction, in radians; S θ For directional similarity;

[0031] Distance similarity is used to measure the spatial proximity of adjacent response units. It converts Euclidean distance into distance similarity. The formula for calculating Euclidean distance is as follows:

[0032]

[0033] In the above formula, (x m y m ) and (x k y k E represents the centroid geographic coordinates of response units m and k, respectively; mk For (x) m y m ) and (x k y k The Euclidean distance between two points, S d For distance similarity;

[0034] Land use similarity is used to measure the similarity of the impermeability of the underlying surface of adjacent response units m and k. The calculation formula is as follows:

[0035]

[0036] In the above formula, C m C represents the overall impermeability of unit m. k S represents the combined impermeability of adjacent units. l For land use similarity, matrix element S mk The calculation formula is as follows:

[0037]

[0038] S25: The flow direction and flow path of each unit are determined by grid calculation to form a surface flow direction matrix, which provides a basis for the fine division of subsequent response units;

[0039] S26: Use the runoff path and flow information obtained from the grid surface runoff representation to correct the boundaries of the vector response cells. If the grid runoff path direction shows that the water flow direction of a certain grid does not point to the designated response cell but to other response cells, the boundaries of the vector response cells should be adjusted to ensure that the final spatial division conforms to the actual path of surface water flow, thereby optimizing the accuracy of the response cell division.

[0040] Further, S3 includes the following steps:

[0041] S31: Extract the waterlogging response time series of each response calculation unit under a typical rainfall event, including the water depth change series, surface runoff change series, peak flow change series, and waterlogging range expansion series, to characterize the waterlogging response process of each unit;

[0042] S32: The time series similarity between any two hydrological response units is calculated using the dynamic time warping method. The waterlogging response time series of response units m and k are as follows:

[0043] U = (u1, u2, ... u) b V = (v1, v, ... v) c )

[0044] In the above formula, u p v p Let p and q be the flood response values ​​at times p and q, respectively. These are two time series that are aligned, and a cumulative distance matrix D is constructed from them. p,q And calculate them one by one using a dynamic programming recursive formula, as follows:

[0045] D p,q =d(x p ,y q )+min{D p-1,q D p,q-1 D p-1,q-1}

[0046] In the above formula, D p-1,q D p,q-1 D p-1,q-1 D is the minimum cumulative distance between the top, left, and top-left sides. p,q Let d(u) be the cumulative distance matrix. p ,v q ) represents the local distance, indicating the difference between the response values ​​at time p and time q, d(u p ,v q The calculation formula for ) is as follows:

[0047] d(v p ,v q )=|v p -vq |

[0048] In the above formula, u p v p These are the flood response values ​​at times p and q, respectively;

[0049] S33: Find the minimum cumulative distance D min (U,V), through the constructed cumulative distance matrix D p,q This eliminates the time lag that may occur in the water depth or runoff curves of different response units due to differences in terrain and drainage capacity, enabling accurate identification of response units with similar flooding responses, with a minimum cumulative distance D. min The formula for calculating (U,V) is as follows:

[0050] D min (U,V)=D P,Q

[0051] In the above formula, D min (U,V) represents the total cost of the optimal path from the starting point of sequence U to the ending point of sequence V, P and Q represent the lengths of time series X and Y, i.e., the number of elements in U and V, respectively, and D P,Q D represents the cumulative matrix D of DTW p,q The bottom right element represents the minimum cumulative cost of the path from the starting point (1,1) to (P,Q);

[0052] S34: The minimum cumulative distance D between all pairs of response units min Using (U,V) as matrix elements, construct the waterlogging response distance matrix. The normalization process was then performed to eliminate dimensional differences, resulting in the waterlogging response similarity matrix S. mk S mk The calculation formula is as follows:

[0053]

[0054] In the above formula, max(D) min S represents the maximum cumulative distance between all pairs of response units. mk The closer the value is to 1, the higher the similarity of the flood response process;

[0055] S35: Based on the similarity matrix S of urban flooding response mk Clustering is performed on the response units, and a similarity threshold μ is used to judge any pair of response units (m,k). If the following conditions are met:

[0056] S mk >μ

[0057] The two response units are considered to have highly similar waterlogging response processes and are merged into the same cluster. This operation eliminates the time lag that may occur in the water depth or runoff curves of different areas due to differences in terrain and drainage capacity, so that similar flood response areas can be accurately identified.

[0058] Further, S4 includes the following steps:

[0059] S41: Based on the clustering results in S35, each initial cluster T... g As input for agglomerative hierarchical clustering, an average connectivity approach is adopted, and a flood response distance matrix between clusters is constructed based on the minimum cumulative DTW distance of response units. The calculation formula is as follows:

[0060]

[0061] In the above formula, T g T r This represents two clusters, each composed of multiple response units. This refers to the minimum cumulative DTW distance between response units g and r calculated in S34. Indicates cluster T g T r The minimum distance between them;

[0062] S42: To avoid the limitations of manually set partitioning rules, agglomerative hierarchical clustering is adopted to make the partitioning process adaptively adjust: merge the cluster pairs with the closest waterlogging response distance, and recalculate the waterlogging response distance and average profile coefficient between the merged new cluster and other clusters. Repeat this step until it is merged into a single cluster.

[0063] S43: Traverse the agglomerative clustering tree and select the cluster with the largest average silhouette coefficient as the optimal number of partitions h* on the time scale. The calculation formula is as follows:

[0064]

[0065] In the above formula, arcmax represents the parameter value that is maximized by the function, backtracking to the clustering situation when the number of partitions is h*;

[0066] S44: Spatial adjacency constraints are applied to the clustering results obtained in S43: If a cluster is not geographically discontinuous, it is split into multiple spatially connected subclusters to ensure that the final partitioning results satisfy the physical connectivity and structural integrity of the urban geographic space. Each cluster after the splitting is the partitioned response unit. Since flood processes are constrained by multiple factors such as topography, drainage networks, and land use, water flow has significant spatial dependence. Therefore, urban flood spatiotemporal response units should not only consider temporal response similarity but also ensure that adjacent areas remain spatially connected, thereby avoiding the generation of "discrete" partitions that do not conform to actual hydrogeographical laws.

[0067] Further, S42 includes the following steps:

[0068] S421: Find the pair of clusters T with the smallest inter-cluster waterlogging response distance. g T r Merge them into a new cluster;

[0069] S422: After merging, the average connectivity method is used to update the flood response distance between the merged cluster and other clusters;

[0070] S423: During each merge, record the current clustering structure and calculate the silhouette coefficient F(h) for different numbers of partitions in real time. The calculation formula is as follows:

[0071]

[0072] In the above formula, for the i-th cluster, let a(i) be its average distance to other clusters, b(i) be its average normalized distance to the nearest neighbor cluster, and S(i) be the silhouette coefficient, which ranges from -1 to 1. The closer it is to 1, the better the clustering effect.

[0073] S424: Repeat the above operation to gradually generate a complete agglomerative clustering tree from the initial N clusters to a single cluster.

[0074] Further, S5 includes the following steps:

[0075] S51: Calculate the area (AREA), average slope (SLOPE), impermeable area ratio (IMPERV), impermeable surface roughness (NIMPERV), and permeable surface roughness (NPERV) for each response unit. Then, calculate the width (WIDTH) based on the merged flow path length (LENGTH) = AREA / LENGTH. The formula for calculating the area (AREA) of the response unit is as follows:

[0076]

[0077] In the above formula, cluster represents the initial set of units within the response unit, and AREA mGiven the area of ​​the initial response element m, the formula for calculating the average slope (SLOPE) of the response element is as follows:

[0078]

[0079] In the above formula, SLOPE m Given the slope of the initial response unit m, the formula for calculating the impermeable area ratio IMPERV of the response unit is as follows:

[0080]

[0081] In the above formula, IMPERV m Given the proportion of the impermeable area of ​​the initial response element m, the formula for calculating the impermeable surface roughness NIMPERV of the response element is as follows:

[0082]

[0083] In the above formula, NIMPERV m Let P be the surface roughness of the impermeable surface of the initial response element m. m,cluster The initial response element represents the area proportion of the response element within its respective response element. The formula for calculating the permeable surface roughness (NPERV) of the response element is as follows:

[0084]

[0085] In the above formula, NPERV m Let P be the surface roughness of the permeable surface of the initial response element m. m,cluster The initial response element represents the area of ​​the response element within its respective response element. The formula for calculating the response element width (WIDTH) is as follows:

[0086]

[0087] In the above formula, LENGTH is the length of the water flow path of the response unit;

[0088] S52: Import the vector data of the response unit division results into Arcmap, create AREA, SLOPE, IMPERV, NIMPERV, NPERV, and WIDTH fields, and write the calculated values ​​into the attribute table to form response unit differentiation parameters that can be directly used for subsequent urban flooding modeling.

[0089] S53: Export the response unit vector data with completed differential parameter assignment for urban flooding modeling in S6. Compared with existing technologies, this step provides more refined underlying surface parameters that take into account spatial heterogeneity through precise assignment of differential response unit parameters, which is essential for subsequent surface-pipeline coupling calculations.

[0090] Further, S6 includes the following steps:

[0091] S61: Clean the underground pipe network data, remove detection wells and sewage inspection wells that have no actual drainage function, and check the connectivity of the pipe network topology.

[0092] S62: Based on the constructed rainfall grid data, spatially match it to the defined response units, and combine it with the assigned differential parameters to simulate the process from rainfall to runoff, generating the runoff inflow for each response unit;

[0093] S63: Load the inflow into the network nodes, use a one-dimensional network model to simulate the drainage process, consider the geometric and hydraulic properties and spatial distribution of the network data, dynamically calculate the inflow, backflow, filling and overflow processes, and reflect the spatial structural differences and time response characteristics of the network system.

[0094] S64: The overflow results output by the one-dimensional pipe network model are used as the inflow boundary of the two-dimensional surface model to drive the response unit to perform water accumulation calculation, thereby realizing the one-way coupling between the one-dimensional pipe network drainage results and the two-dimensional surface water accumulation process, and accurately simulating the water accumulation depth distribution and waterlogging evolution process of different response units.

[0095] Compared with the prior art, the advantages of the present invention are as follows:

[0096] (1) This invention combines geographical laws and introduces the theory of spatiotemporal heterogeneity into urban flooding modeling, analyzes the geographical modeling mechanism contained in the urban rainfall process and response unit division process under flooding scenario, and improves the application of GIS theory in urban flooding modeling.

[0097] (2) This invention overcomes the limitations of traditional modeling, which makes it difficult to fully reflect the spatiotemporal heterogeneity of urban characteristics and rainfall nonstationarity on the impact of urban flooding processes. It incorporates multi-scale spatiotemporal grids of rainfall and differentiated hydrological response units into urban flooding modeling, which can dynamically reflect the rainfall change characteristics of different regions at different times, effectively improving the accuracy and reliability of urban flooding process simulation.

[0098] (3) This invention introduces multi-dimensional indicators such as topographic features, impermeability, land use similarity, runoff path, spatial connectivity and drainage system distribution, and establishes a dual constraint mechanism of spatial adjacency and time series clustering to realize the realistic reproduction of the complex surface and pipe network system inside the city, thereby improving the model's adaptability to different urban forms. Attached Figure Description

[0099] To more clearly illustrate the specific embodiments of the present invention or the technical solutions in the prior art, the accompanying drawings used in the description of the specific embodiments or the prior art will be briefly introduced below.

[0100] Figure 1 Method framework diagram of the present invention;

[0101] Figure 2 The initial response unit partitioning result diagram of the modeling region provided in this embodiment of the invention;

[0102] Figure 3 Clustering results of modeling region response units provided in this embodiment of the invention;

[0103] Figure 4 The resulting diagram of the underground drainage network in the modeling area provided in this embodiment of the invention;

[0104] Figure 5 The embodiments of the present invention provide a diagram showing the changes in inflow, inundated area, and inundated volume of the modeled regional flood response process. Detailed Implementation

[0105] To make the objectives, technical solutions, and advantages of the present invention clearer, the present invention will be further described in detail below with reference to the embodiments and accompanying drawings. The illustrative embodiments and descriptions of the present invention are only used to explain the present invention and are not intended to limit the present invention.

[0106] Taking a certain region as an example, this area is a certain urban area. The overall terrain presents a pattern of "high surrounding areas and low center." The region includes numerous squares, parks, hospitals, schools, and other living spaces, with comprehensive supporting facilities. It also includes several lakes and parks, natural water storage facilities, and a high degree of ground hardening, with most roads made of cement and asphalt. The rivers surrounding the region are all dammed, and the pipeline network is relatively independent, forming a relatively independent drainage unit.

[0107] according to Figure 1 As shown, this invention provides a method for urban flood modeling that takes into account spatiotemporal heterogeneity, specifically including the following steps:

[0108] Step S1: Construct a multi-scale spatiotemporal grid of rainfall with spatiotemporal heterogeneity;

[0109] Step S2: Data preprocessing, dividing the study area into initial response units based on topographic features, land use types, and drainage system distribution;

[0110] Step S3: Use the dynamic time warping method to measure the similarity of the waterlogging response process among different response units, and further perform time series clustering based on the partitioning;

[0111] Step S4: Adopt agglomerative hierarchical clustering to adjust the partitioning process and calculate the optimal number of partitions, and introduce spatial adjacency constraints to ensure spatial connectivity;

[0112] Step S5: Assign differentiated hydrological parameters to different response units;

[0113] Step S6: Use multi-scale rainfall spatiotemporal grid data and response unit partitioning results to drive the surface-pipeline coupling model of urban flooding.

[0114] Further, S1 includes the following steps:

[0115] Step S11: Based on the time axis length and time resolution of the rainfall to be inverted, divide the time interval t0~t n Divide the time slices into n equal layers based on time steps. Each layer precisely corresponds to the spatial rainfall distribution over a specific time period. Each time slice satisfies the following:

[0116]

[0117] In the above formula, t0 represents the start time of the rainfall process, t n The time of the end of the rainfall event is indicated by t, where t represents the time step of the main time slice, and n represents the number of main time slices in the rainfall event. i This represents the start time corresponding to the i-th main time slice;

[0118] Step S12: When retrieving rainfall observation points, attach the target time slice [t] i， t i The constraint condition +t] is used to filter valid point rainfall data located in the modeling region Ω within this sub-time slice. Based on the filtering results, the main time slice is subdivided into k sub-time slices, specifically satisfying:

[0119] S i ={(x j ,y j ,R j )|(x j ,y j )∈Ω,t j ∈[t i ,t i +t]}

[0120] In the above formula, Ω represents the spatial extent of the modeling region, and x i ,y j Let R represent the coordinates of the j-th observation point. j Let t represent the rainfall observation value at the j-th observation point. j S represents the observation time recorded at the j-th observation point, k represents the number of sub-time slices within the current main time slice, and S represents the sub-time slice. i This represents the set of valid rainfall observation data within the i-th main time slice and located in the modeling region Ω;

[0121] S13: Based on additional screening objectives [t] i , t iSpatial interpolation is performed on the point rainfall data obtained by filtering the time segments of +t] to obtain the rainfall grid dataset within the main time segment with index i. The rainfall grid data within the main time segment are summed to accumulate the rainfall, thus obtaining the rainfall grid data of the main time segment with index i.

[0122] S14: Repeat the above process until all n time slices are calculated. Finally, stack them along the time axis to build a rainfall grid dataset that takes into account multiple time scales and multiple spatial resolutions, providing high-precision rainfall driving force for subsequent urban flood models.

[0123] Further, S2 includes the following steps:

[0124] S21: Using high-resolution DEM as the basic spatial data, the study area is divided into regular and uniform grid cells with clear spatial topological relationships between cells. The DEM is then subjected to terrain correction and spatial preprocessing to remove errors.

[0125] S22: The geographic location of each rainwater well node is expressed using a vector structure, which is regarded as the initial control node of spatial confluence. With each rainwater well as the core, the Thiessen polygon method is used to determine the area with the closest spatial distance, and it is initially classified as the response unit of the rainwater well.

[0126] S23: Spatial matching of land use types in the modeling area with grid cells to determine the impermeability of each cell, and the overall impermeability C of the response cell. m The calculation formula is as follows:

[0127]

[0128] In the above formula, C m Let A be the overall impermeability of unit m. m,l Let α be the area of ​​land use type l within unit m. l Let A be the impermeability coefficient corresponding to land use type l, where L represents the number of land use types included. m The total area of ​​unit m; the impermeability coefficients for different land use types are shown in Table 1:

[0129] Table 1. Land Use Type Parameter Assignment Table

[0130]

[0131] S24: Introducing directional similarity, distance similarity, and land use similarity indices, a joint similarity matrix S is established. S is used to optimize the Thiessen polygon division results. Directional similarity is mainly determined by calculating the angle between the runoff directions of adjacent response units to determine spatial correlation. The calculation formula is as follows:

[0132]

[0133] In the above formula, θ is the angle between the line connecting the centroids of response unit m and its adjacent response unit k and its own runoff direction, in radians; S θ For directional similarity;

[0134] Distance similarity is used to measure the spatial proximity of adjacent response units. It converts Euclidean distance into distance similarity. The formula for calculating Euclidean distance is as follows:

[0135]

[0136] In the above formula, (x m y m ) and (x k y k E represents the centroid geographic coordinates of response units m and k, respectively; mk For (x) m y m ) and (x k y k The Euclidean distance between two points, S d For distance similarity;

[0137] Land use similarity is used to measure the similarity of the impermeability of the underlying surface of adjacent response units m and k. The calculation formula is as follows:

[0138]

[0139] In the above formula, C m C represents the overall impermeability of unit m. k S represents the combined impermeability of adjacent units. l For land use similarity, matrix element S mk The calculation formula is as follows:

[0140]

[0141] S25: The flow direction and flow path of each cell are determined by grid calculation to form a surface flow direction matrix;

[0142] S26: Using the runoff path and flow information obtained from the raster surface runoff representation, the boundaries of the vector response cells are corrected. If the raster runoff path direction shows that the water catchment direction of a certain raster does not point to the designated response cell but to other response cells, the boundaries of the vector response cells should be adjusted to ensure that the final spatial division conforms to the actual path of surface water flow, thereby optimizing the accuracy of the response cell division. The final initial response cells are as follows: Figure 2 As shown.

[0143] Further, S3 includes the following steps:

[0144] S31: Extract the waterlogging response time series of each response calculation unit under a typical rainfall event, including the water depth change series, surface runoff change series, peak flow change series, and waterlogging range expansion series, to characterize the waterlogging response process of each unit;

[0145] S32: The time series similarity between any two hydrological response units is calculated using the dynamic time warping method. The waterlogging response time series of response units m and k are as follows:

[0146] U = (u1, u2, ... u) b V = (v1, v, ... v) c )

[0147] In the above formula, u p v p Let p and q be the flood response values ​​at times p and q, respectively. These are two time series that are aligned, and a cumulative distance matrix D is constructed from them. p,q And calculate them one by one using a dynamic programming recursive formula, as follows:

[0148] D p,q =d(x p ,y q )+min{D p-1,q D p,q-1 D p-1,q-1}

[0149] In the above formula, D p-1,q D p,q-1 D p-1,q-1 D is the minimum cumulative distance between the top, left, and top-left sides. p,q Let d(u) be the cumulative distance matrix. p ,v q ) represents the local distance, indicating the difference between the response values ​​at time p and time q, d(u p ,v q The calculation formula for ) is as follows:

[0150] d(u p ,v q )=|u p -v q |

[0151] In the above formula, u p v p These are the flood response values ​​at times p and q, respectively;

[0152] S33: Find the minimum cumulative distance D min (U,V), through the constructed cumulative distance matrix D p,qThis eliminates the time lag that may occur in the water depth or runoff curves of different response units due to differences in terrain and drainage capacity, enabling accurate identification of response units with similar flooding responses, with a minimum cumulative distance D. min The formula for calculating (U,V) is as follows:

[0153] D min (U,V)=D P,Q

[0154] In the above formula, D min (U,V) represents the total cost of the optimal path from the starting point of sequence U to the ending point of sequence V, P and Q represent the lengths of time series X and Y, i.e., the number of elements in U and V, respectively, and D P,Q D represents the cumulative matrix D of DTW p,q The bottom right element represents the minimum cumulative cost of the path from the starting point (1,1) to (P,Q);

[0155] S34: The minimum cumulative distance D between all pairs of response units min Using (U,V) as matrix elements, construct the waterlogging response distance matrix. The normalization process was then performed to eliminate dimensional differences, resulting in the waterlogging response similarity matrix S. mk S mk The calculation formula is as follows:

[0156]

[0157] In the above formula, max(D) min S represents the maximum cumulative distance between all pairs of response units. mk The closer the value is to 1, the higher the similarity of the flood response process;

[0158] S35: Based on the similarity matrix S of urban flooding response mk Clustering is performed on the response units, and a similarity threshold μ is used to judge any pair of response units (m,k). If the following conditions are met:

[0159] S mk >μ

[0160] If two response units are considered to have highly similar waterlogging response processes, they are merged into the same cluster, and each of the final merged clusters is regarded as the initial cluster.

[0161] Further, S4 includes the following steps:

[0162] S41: Based on the clustering results in S35, each initial cluster T... g As input for agglomerative hierarchical clustering, an average connectivity approach is adopted, and a flood response distance matrix between clusters is constructed based on the minimum cumulative DTW distance of response units. The calculation formula is as follows:

[0163]

[0164] In the above formula, T g T r This represents two clusters, each composed of multiple response units. This refers to the minimum cumulative DTW distance between response units g and r calculated in S34. Indicates cluster T g T r The minimum distance between them;

[0165] S42: Adopt agglomerative hierarchical clustering to make the partitioning process adaptively adjust: merge the cluster pairs with the closest waterlogging response distance, and recalculate the waterlogging response distance and average profile coefficient between the merged new cluster and other clusters. Repeat this step until they are merged into a single cluster.

[0166] S43: Traverse the agglomerative clustering tree and select the cluster with the largest average silhouette coefficient as the optimal number of partitions h* on the time scale. The calculation formula is as follows:

[0167]

[0168] In the above formula, arcmax represents the parameter value that is maximized by the function, backtracking to the clustering situation when the number of partitions is h*;

[0169] S44: Apply spatial adjacency constraints to the clustering results obtained in S43: If a cluster is not geographically contiguous, split it into multiple connected clusters to ensure that the final partitioning results satisfy the physical connectivity of the urban geographic space. Each cluster after splitting is the partitioned response unit. The final response unit clustering results are as follows: Figure 3 As shown.

[0170] Further, S42 includes the following steps:

[0171] S421: Find the pair of clusters T with the smallest inter-cluster waterlogging response distance. g T r Merge them into a new cluster;

[0172] S422: After merging, the average connectivity method is used to update the flood response distance between the merged cluster and other clusters;

[0173] S423: During each merge, record the current clustering structure and calculate the silhouette coefficient F(h) for different numbers of partitions in real time. The calculation formula is as follows:

[0174]

[0175] In the above formula, for the β-th cluster, z(β) is its average distance to other clusters, w(β) is its average normalized distance to the nearest neighbor cluster, and F(β) is the silhouette coefficient, which ranges from -1 to 1. The closer it is to 1, the better the clustering effect.

[0176] S424: Repeat the above operation to gradually generate a complete agglomerative clustering tree from the initial N clusters to a single cluster.

[0177] Further, S5 includes the following steps:

[0178] S51: Step S51: Calculate the area (AREA), average slope (SLOPE), impermeable area ratio (IMPERV), impermeable surface roughness (NIMPERV), and permeable surface roughness (NPERV) for each response unit. Then, calculate the width (WIDTH) based on the merged flow path length (LENGTH) = AREA / LENGTH. The formula for calculating the area (AREA) of the response unit is as follows:

[0179]

[0180] In the above formula, cluster represents the initial set of units within the response unit, and AREA m Given the area of ​​the initial response element m, the formula for calculating the average slope (SLOPE) of the response element is as follows:

[0181]

[0182] In the above formula, SLOPE m Given the slope of the initial response unit m, the formula for calculating the impermeable area ratio IMPERV of the response unit is as follows:

[0183]

[0184] In the above formula, IMPERV m Given the proportion of the impermeable area of ​​the initial response element m, the formula for calculating the impermeable surface roughness NIMPERV of the response element is as follows:

[0185]

[0186] In the above formula, NIMPERV m Let P be the surface roughness of the impermeable surface of the initial response element m. m,cluster The initial response element represents the area proportion of the response element within its respective response element. The formula for calculating the permeable surface roughness (NPERV) of the response element is as follows:

[0187]

[0188] In the above formula, NPERV m Let P be the surface roughness of the permeable surface of the initial response element m.m,cluster The initial response element represents the area of ​​the response element within its respective response element. The formula for calculating the response element width (WIDTH) is as follows:

[0189]

[0190] In the above formula, LENGTH is the length of the water flow path of the response unit;

[0191] S52: Import the vector data of the response unit division results into Arcmap, create AREA, SLOPE, IMPERV, NIMPERV, NPERV, and WIDTH fields, and write the calculated values ​​into the attribute table to form response unit differentiation parameters that can be directly used for subsequent urban flooding modeling.

[0192] S53: Export the response unit vector data with completed differential parameter assignment for use in urban flooding modeling in S6;

[0193] Further, S6 includes the following steps:

[0194] S61: The underground pipe network data is cleaned. This invention uses a breadth-first search algorithm to traverse rainwater wells and rainwater grate nodes, deleting detection wells, sewage wells, and isolated nodes that have not been traversed, retaining the main rainwater well nodes, and removing some secondary pipe networks and other functional pipe networks. After data cleaning, 1600 rainwater well nodes, 1607 rainwater pipe segments, and 7 drainage outlets are finally obtained. Based on the node ID and the inlet / outlet node ID attribute fields of the pipe segments, ArcGIS is used to check the pipe network connectivity, and some pipe segment information with conflicting flow directions is corrected. The final pipe network processing result is as follows: Figure 4 As shown;

[0195] S62: Based on the constructed rainfall grid data, spatially match it to the defined response units, and combine it with the assigned differential parameters to simulate the process from rainfall to runoff, generating the runoff inflow for each response unit;

[0196] S63: Load the inflow into the network nodes, use a one-dimensional network model to simulate the drainage process, consider the geometric and hydraulic properties and spatial distribution of the network data, dynamically calculate the inflow, backflow, filling and overflow processes, and reflect the spatial structural differences and time response characteristics of the network system.

[0197] S64: The overflow results output from the one-dimensional pipe network model are used as the inflow boundary of the two-dimensional surface model to drive the response unit to perform water accumulation calculations, realizing unidirectional coupling between the one-dimensional pipe network drainage results and the two-dimensional surface water accumulation process. This allows for precise simulation of the water accumulation depth distribution and waterlogging evolution process in different response units. The changes in inflow, inundated area, and inundated volume during the flood response process are shown in the following figures. Figure 5 As shown.

[0198] The above description of the embodiments is only for the purpose of helping to understand the method and core idea of ​​the present invention; at the same time, for those skilled in the art, there will be changes in the specific implementation and application scope based on the idea of ​​the present invention. Therefore, the content of this specification should not be construed as a limitation of the present invention.

Claims

1. A method for modeling urban waterlogging considering spatio-temporal heterogeneity, characterized in that, Includes the following steps: Step S1: Construct a multi-scale spatiotemporal grid of rainfall with spatiotemporal heterogeneity; Step S2: Data preprocessing, dividing the study area into initial response units based on topographic features, land use types, and drainage system distribution; Step S3: Use the dynamic time warping method to measure the similarity of the waterlogging response process among the response units, and perform time series clustering on the partitioning results accordingly; Step S4: Further clustering is performed using agglomerative hierarchical clustering to determine the optimal number of partitions, and spatial adjacency constraints are introduced to ensure spatial connectivity; Step S5: Assign differentiated hydrological parameters to different response units; Step S6: Use multi-scale rainfall spatiotemporal grid data and response unit partitioning results to drive the surface-pipeline coupling model of urban flooding; Step S1 includes the following steps: Step S11: based on the time axis length and time resolution of the rainfall to be inverted, the time interval t0~t n Divided into n layers according to time steps, each layer accurately corresponds to the spatial rainfall distribution of a specific time period, and each layer time slice satisfies: In the above formula, t0 represents the start time of the rainfall process, t n The time of the end of the rainfall event is indicated by t, where t represents the time step of the main time slice, and n represents the number of main time slices in the rainfall event. i This represents the start time corresponding to the i-th main time slice; Step S12: additional constraint condition of target time slice [t i， t i +t] is added when searching for rainfall observation points, effective point rainfall data located in the modeling area Ω within the target time slice is screened, and the main time slice is subdivided into k sub-time slices according to the screening result, which specifically meets: In the above formula, Ω represents the spatial extent of the modeling region, and x i ,y j Let R represent the coordinates of the j-th observation point. j Let t represent the rainfall observation value at the j-th observation point. j S represents the observation time recorded at the j-th observation point, k represents the number of sub-time slices within the current main time slice, and S represents the observation time recorded at the j-th observation point. i This represents the set of valid rainfall observation data within the i-th main time slice and located in the modeling region Ω; Step S13: Perform spatial interpolation calculation on the point rainfall data obtained by filtering the time segments based on the additional filtering target to generate a rainfall grid dataset within the main time segment with index i, and sum the rainfall grid data within the main time segment to accumulate the rainfall amount, thereby obtaining the rainfall grid data of the main time segment with index i. Step S14: Repeat the above process until all n time slices are calculated. Finally, stack them along the time axis to construct a rainfall grid dataset that takes into account multiple time scales and multiple spatial resolutions, providing high-precision rainfall driving force for subsequent urban flooding models.

2. The urban flooding modeling method considering spatiotemporal heterogeneity according to claim 1, characterized in that, S2 includes the following steps: Step S21: Using high-resolution DEM as the basic spatial data, the study area is divided into regular and uniform grid cells with clear spatial topological relationships between cells. The DEM is then subjected to terrain correction and spatial preprocessing to remove errors. Step S22: Use a vector structure to express the geographical location of each rainwater well node, and regard it as the initial control node of spatial confluence. With each rainwater well as the core, use the Thiessen polygon method to determine the area with the closest spatial distance, and initially classify it as the response unit of the rainwater well. Step S23: Spatial match the land use type of the modeling area with the grid cells to determine the impermeability of each cell, and the overall impermeability C of the response cell. m The calculation formula is as follows: In the above formula, C m Let A be the overall impermeability of unit m. m,l Let m be the area of ​​land use type l within unit m. Let A be the impermeability coefficient corresponding to land use type l, where L represents the number of land use types included. m The total area of ​​unit m; Step S24: Further introduce directional similarity, distance similarity, and land use similarity indices to establish a joint similarity matrix S. Use S to optimize the Thiessen polygon division results. Directional similarity is mainly determined by calculating the angle between the runoff directions of adjacent response units to determine spatial correlation. The calculation formula is as follows: In the above formula, is the angle between the line connecting the response unit m and its adjacent response unit k and its own direction of runoff, in radians; is the direction similarity; Distance similarity is used to measure the spatial proximity of adjacent response units. It converts Euclidean distance into distance similarity. The formula for calculating Euclidean distance is as follows: In the above formula, (x m y m ) and (x k y k () represent the centroid geographic coordinates of response units m and k, respectively; For (x) m y m ) and (x k y k The Euclidean distance between two points. For distance similarity; Land use similarity is used to measure the similarity of the impermeability of the underlying surface of adjacent response units m and k. The calculation formula is as follows: In the above formula, C m C represents the overall impermeability of unit m. k The combined impermeability of adjacent units, For land use similarity, matrix element S mk The calculation formula is as follows: Step S25: Determine the flow direction and flow path of each cell through grid calculation to form a surface flow direction matrix, which provides a basis for the subsequent fine division of response cells; Step S26: Using the runoff path and flow information obtained from the grid surface runoff representation, the boundary of the vector response unit is corrected. If the direction of the grid runoff path shows that the water flow direction of a certain grid does not point to the designated response unit but to other response units, the boundary of the vector response unit should be adjusted to ensure that the final spatial division conforms to the real path of surface water flow, thereby optimizing the accuracy of the response unit division.

3. The urban flooding modeling method considering spatiotemporal heterogeneity according to claim 1, characterized in that, S3 includes the following steps: Step S31: Extract the waterlogging response time series of each response calculation unit under a typical rainfall event, including the water depth change series, surface runoff change series, peak flow change series, and waterlogging range expansion series, to characterize the waterlogging response process of each unit; S32: The time series similarity between any two hydrological response units is calculated using the dynamic time warping method. The waterlogging response time series of response units m and k are as follows: In the formula, u p , v p are the waterlogging response values at time p, q respectively, and both are aligned with two time series, a cumulative distance matrix D p,q is constructed, and is calculated one by one through dynamic programming recursive formula, and the calculation formula is as follows: In the above formula, This represents the minimum cumulative distance to the top, left, and top-left. Let d(u) be the cumulative distance matrix. p ,v q ) represents the local distance, indicating the difference between the response values ​​at time p and time q, d(u p ,v q The calculation formula for ) is as follows: In the above formula, u p , v p are the waterlogging response values at time p, q, respectively. Step S33: Calculate the minimum cumulative distance D min (U,V), through the constructed cumulative distance matrix D p,q This eliminates the time lag that may occur in the water depth or runoff curves of different response units due to differences in terrain and drainage capacity, enabling accurate identification of response units with similar flooding responses, with a minimum cumulative distance D. min The formula for calculating (U,V) is as follows: In the above formula, D min (U,V) represents the total cost of the optimal path from the starting point of sequence U to the ending point of sequence V, P and Q represent the lengths of time series X and Y, i.e., the number of elements in U and V, respectively, and D P,Q D represents the cumulative matrix D of DTW p,q The bottom right element represents the minimum cumulative cost of the path from the starting point (1,1) to (P,Q); Step S34: Calculate the minimum accumulated distance D between all response units two by two min (U, V) as matrix elements, construct the waterlogging response distance matrix , and normalize to eliminate dimensional differences to obtain the waterlogging response similarity matrix S mk , S mk The calculation formula is as follows: In the above equation, max(D min ) is the maximum accumulated distance value between all pairs of response units, S mk The closer to 1 indicates the higher similarity of the waterlogging response process; Step S35: Based on the waterlogging response similarity matrix S mk The clustering operation is performed on the response units, and the similarity threshold μ is used to judge any response unit pair (m, k). If the following condition is met: If two response units are considered to have highly similar waterlogging response processes, they are merged into the same cluster, and each of the final merged clusters is regarded as the initial cluster.

4. The urban waterlogging modeling method of claim 1, wherein, S4 includes the following steps: Step S41: divide each initial cluster T g As the input of the agglomerative hierarchical clustering, the average linkage method is adopted to construct the matrix of the response distance between clusters based on the minimum cumulative DTW distance of the response units The calculation formula is as follows: In the above formula, T g T r This represents two clusters, each composed of multiple response units. This refers to the minimum cumulative DTW distance between response units g and r calculated in S34. Represents cluster T g T r The minimum distance between them; Step S42: Adopt agglomerative hierarchical clustering to make the partitioning process adaptively adjust: merge the cluster pairs with the closest waterlogging response distance, and recalculate the waterlogging response distance and average profile coefficient between the merged new cluster and other clusters. Repeat this step until they are merged into a single cluster. Step S43: Traverse the agglomerative clustering tree and select the cluster with the largest average silhouette coefficient as the optimal number of partitions h* on the time scale. The calculation formula is as follows: In the above formula, arcmax represents the parameter value that is maximized by the function, backtracking to the clustering situation when the number of partitions is h*; Step S44: Apply spatial adjacency constraints to the clustering results obtained in S43: If a cluster is not continuous in geographic space, split it into multiple connected clusters to ensure that the final partitioning result satisfies the physical connectivity of the urban geographic space. Each cluster after the split is the partitioned response unit.

5. The urban waterlogging modeling method of claim 4, wherein, S42 includes the following steps: Step S421: find a pair of clusters T with the minimum inter-cluster waterlogging response distance g , T r merge them into a new cluster; Step S422: After merging, update the flood response distance between the merged cluster and other clusters using the average connectivity method; Step S423: During each merge, record the current clustering structure and calculate the silhouette coefficient F(h) for different numbers of partitions in real time. The calculation formula is as follows: In the above formula, for the β-th cluster, z(β) is its average distance to other clusters, w(β) is its average normalized distance to the nearest neighbor cluster, and F(β) is the silhouette coefficient, which ranges from -1 to 1. The closer it is to 1, the better the clustering effect. Step S424: Repeat the above operation to gradually generate a complete agglomerative clustering tree from the initial N clusters to a single cluster.

6. The urban waterlogging modeling method of claim 1, wherein, The S5 step includes the following steps: Step S51: Calculate the area (AREA), average slope (SLOPE), impermeable area ratio (IMPERV), impermeable surface roughness (NIMPERV), and permeable surface roughness (NPERV) for each response unit. Then, calculate the width (WIDTH) based on the merged flow path length (LENGTH) = AREA / LENGTH. The formula for calculating the area (AREA) of the response unit is as follows: In the above formula, cluster represents the initial set of units within the response unit, and AREA m Given the area of ​​the initial response element m, the formula for calculating the average slope (SLOPE) of the response element is as follows: In the above formula, SLOPE m Given the slope of the initial response unit m, the formula for calculating the impermeable area ratio IMPERV of the response unit is as follows: In the above formula, IMPERV m The proportion of the impermeable area of the initial response unit m is calculated by the following formula: In the above formula, NIMPERV m is the impervious surface roughness of the initial response unit m m,cluster is the area proportion of the initial response unit in the response unit where it is located, and the response unit impervious surface roughness NPERV is calculated according to the following formula: In the above formula, NPERV m is the water permeable surface roughness of the initial response unit m m,cluster is the area proportion of the initial response unit in the response unit in which it is located, and the response unit width WIDTH is calculated according to the following formula: In the above formula, LENGTH is the length of the water flow path of the response unit; Step S52: Import the vector data of the response unit partitioning results into Arcmap, create fields AREA, SLOPE, IMPERV, NIMPERV, NPERV, and WIDTH, and write the calculated values ​​into the attribute table to form response unit differentiation parameters that can be directly used for subsequent urban flooding modeling. Step S53: Export the response unit vector data with completed differential parameter assignment for use in urban flooding modeling in S6.

7. The urban flooding modeling method considering spatiotemporal heterogeneity according to claim 1, characterized in that, S6 includes the following steps: Step S61: Clean the underground pipe network data, remove detection wells and sewage inspection wells that have no actual drainage function, and check the connectivity of the pipe network topology. Step S62: Based on the constructed rainfall grid data, spatially match it to the defined response units, and combine it with the assigned differential parameters to simulate the process from rainfall to runoff, generating the runoff inflow for each response unit; Step S63: Load the runoff into the pipe network nodes, use a one-dimensional pipe network model to simulate the drainage process, consider the geometric and hydraulic properties and spatial distribution of the pipe network data, dynamically calculate the inflow, backflow, filling and overflow processes, and reflect the spatial structural differences and time response characteristics of the pipe network system. Step S64: Use the overflow results output by the one-dimensional pipe network model as the inflow boundary of the two-dimensional surface model to drive the response unit to perform water accumulation calculation, realize the one-way coupling between the one-dimensional pipe network drainage results and the two-dimensional surface water accumulation process, and finely simulate the water accumulation depth distribution and waterlogging evolution process of different response units.

8. An electronic device comprising a memory, a processor, and a computer program stored on the memory and executable on the processor, characterized in that, When the processor executes the program, it implements the urban flooding modeling method that takes into account spatiotemporal heterogeneity as described in any one of claims 1 to 7.

9. A computer readable storage medium having stored thereon computer instructions, wherein, When the computer instruction is executed by the processor, it implements an urban flooding modeling method that takes into account spatiotemporal heterogeneity as described in any one of claims 1-7.

Images