Urban inland inundation modeling method considering spatial-temporal heterogeneity

Abstract

The invention relates to an urban inland inundation modeling method considering spatio-temporal heterogeneity. The urban inland inundation modeling method comprises the following steps: constructing a rainfall time sequence and a multi-scale spatio-temporal grid with spatio-temporal heterogeneity; data preprocessing: dividing a research area into initial response units based on topographic features, land utilization types and drainage system distribution; measuring the waterlogging response process similarity among the response units by adopting a dynamic time warping method, and performing time sequence clustering on a division result according to the waterlogging response process similarity; further clustering and determining an optimal partition number by adopting an agglomerated hierarchical clustering method, and introducing a spatial adjacency constraint to ensure spatial connectivity; different hydrological parameters are given to different response units; and utilizing the multi-scale rainfall space-time grid data and the response unit division result to drive modeling of an urban inland inundation earth surface-pipe network coupling model. According to the method, the precision of urban waterlogging modeling is improved by constructing the multi-scale rainfall grid and dividing the urban waterlogging response units, and the modeling and simulation method of an existing urban waterlogging model is expanded.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of geographic information and hydrology, and particularly designs a city waterlogging modeling method considering spatio-temporal heterogeneity. BACKGROUND

[0002] Under the background of global climate change and rapid urbanization, extreme rainfall events occur frequently and intensively, which exceeds the design capacity of the traditional urban drainage system. At the same time, the rapid expansion of the hardening area of the ground surface leads to a large number of impervious surfaces replacing the natural underlying surface, which seriously weakens the natural infiltration and storage capacity of rainwater. The development trend of urban flood disasters is an increase in the frequency of occurrence, an expansion of the disaster range, and an intensification of the damage degree. As a typical "non-engineering" measure in the drainage and flood control system, city waterlogging modeling and simulation is a powerful support for flood disaster prediction and early warning.

[0003] The response unit division of the traditional city waterlogging modeling usually does not consider spatio-temporal heterogeneity and assumes that the hydrological parameters of the study area are homogenously distributed, which ignores the significant spatio-temporal heterogeneity in the city waterlogging process. It is difficult to fully reflect the influence of the spatio-temporal heterogeneity characteristics such as city local features and non-stationary rainfall on the city waterlogging process, which leads to the fact that the response unit division and parameter setting of the model cannot accurately reflect the spatio-temporal dynamic differences of the city internal topography, underlying surface conditions, drainage system, and human activities, thereby reducing the accuracy of the model in predicting the location, time, and intensity of city waterlogging and limiting its adaptability to complex scenarios.

[0004] Therefore, the present application proposes a city waterlogging modeling method considering spatio-temporal heterogeneity. Under the condition of considering spatio-temporal heterogeneity, various factors such as rainfall distribution, topographic features, land use types, and drainage system distribution are comprehensively considered to construct a surface-pipe network coupling model, realize high-precision simulation of the city waterlogging process, and provide scientific basis and technical support for city flood control and drainage decision-making. SUMMARY

[0005] The present application provides a city waterlogging modeling method considering spatio-temporal heterogeneity, which aims to introduce spatio-temporal heterogeneity into city waterlogging modeling, break through the limitations of the traditional city waterlogging modeling homogeneity assumption, realize the optimized adjustment of the model parameters according to the environmental conditions, thereby improving the precision and spatio-temporal adaptability of city waterlogging simulation, and providing more reliable scientific support for city flood control and disaster reduction decision-making.

[0006] The present application provides a city rainwater model modeling method based on city waterlogging response unit division, which is characterized by the following steps:

[0007] Step S1: constructing a multi-scale rainfall spatio-temporal grid with spatio-temporal heterogeneity;

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

[0009] Step S3: using dynamic time warping method to measure the similarity of waterlogging response process between response units, and clustering the division results according to the similarity;

[0010] Step S4: further clustering and determining the optimal number of partitions by using agglomerative hierarchical clustering method, and introducing spatial adjacency constraint to ensure spatial connectivity;

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

[0012] Step S6: using multi-scale rainfall spatio-temporal grid data and response unit division results to drive the surface-pipe network coupled model of urban waterlogging.

[0013] Further, the S1 comprises the following steps:

[0014] Step S11: based on the length of the time axis of the rainfall to be inverted and the time resolution, the time interval t0~t n is divided into n layers, each layer accurately corresponds to the spatial rainfall distribution of a specific time period, and each layer time slice satisfies:

[0015]

[0016] In the above formula, t0 represents the starting time of the rainfall process, t n represents the ending time of the rainfall process, t represents the time step of the main time slice, n represents the number of main time slices of the rainfall process, t i represents the starting time corresponding to the i-th main time slice;

[0017] Step S12: adding the constraint condition of the target time slice [t i， t i +t] when searching for rainfall observation points, screening the effective point rainfall data located in the modeling area Ω within the sub-time slice, and subdividing the main time slice into k sub-time slices according to the screening result, which satisfies:

[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 range of the modeling area, xi , y j represents the coordinates of the jth observation point, R j represents the rainfall observation value of the jth observation point, t j represents the observation time recorded at the jth observation point, k represents the number of subdivided sub-time slices within the current main time slice, S i represents the set of valid rainfall observation data within the i th main time slice and located in the modeling region Ω;

[0020] S13: Based on the rainfall point data obtained by filtering the time slice, the algorithm parameters are set according to the spatial scale and interpolation range, and the interpolation is performed to generate the rainfall grid data within the minimum time slice, and the rainfall distribution at the time index i is segmented and summed to accumulate the rainfall amount;

[0021] S14: The above process is repeatedly executed until all n layers of time slices are calculated, and finally the rainfall grid data set considering multiple time scales and multiple spatial resolutions is constructed by stacking along the time axis, which provides high-precision rainfall driving for subsequent urban flood models.

[0022] Further, the S2 comprises the following steps:

[0023] S21: Taking high-resolution DEM as the basic spatial data, the study area is divided into regular and uniform grid cells, and the cells have clear spatial topological relationship. The DEM is corrected and preprocessed in space to remove errors;

[0024] S22: The geographic position of each rainwater well node is expressed by vector structure, which is regarded as the initial control node of spatial confluence. Taking each rainwater well as the core, the spatial distance of the nearest area is determined by the Thiessen polygon method, and it is preliminarily classified as the response unit of the rainwater well;

[0025] S23: The land use types of the modeling region are spatially matched with the grid cells to determine the impervious rate of each cell, and the impervious rate C m of the response unit is calculated according to the following formula:

[0026]

[0027] In the above formula, C m is the comprehensive impervious rate of cell m, A m,l is the area of land use type l in cell m, α l is the impervious coefficient corresponding to land use type l, L represents the number of contained land use types, A m is the total area of cell m;

[0028] S24: Further introduce the direction similarity, distance similarity, land use similarity index, establish the joint similarity matrix S, use S to optimize the division result of Thiessen polygon, the direction similarity is mainly determined by calculating the included angle between the runoff directions of adjacent response units, the calculation formula is as follows:

[0029]

[0030] In the above formula, θ is the included angle between the centroid connecting line of response unit m and its adjacent response unit k and its own runoff direction, unit is radian; S θ is the direction similarity;

[0031] The distance similarity is used to measure the spatial proximity of adjacent response units, and the Euclidean distance is converted into the distance similarity, the calculation formula of Euclidean distance is as follows:

[0032]

[0033] 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; E mk is the Euclidean distance between (x m , y m ) and (x k , y k ), S d is the distance similarity;

[0034] The land use similarity is used to measure the similarity of the impervious rate of the underlying surface of adjacent response units m and k, and the calculation formula is as follows:

[0035]

[0036] In the above formula, C m is the comprehensive impervious rate of unit m, C k is the comprehensive impervious rate of adjacent unit, S l is the land use similarity, and the calculation formula of matrix element S mk is as follows:

[0037]

[0038] S25: Determine the runoff direction and runoff path of each unit through grid calculation, form the surface flow direction matrix, and provide basis for the fine division of subsequent response units;

[0039] S26: The obtained runoff path and flow information is used to correct the boundary of the vector response unit by using the grid surface runoff expression. If the grid runoff path direction shows that the grid catchment direction is not directed to the delineated 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 movement, thereby optimizing the accuracy of the response unit division.

[0040] Further, the S3 comprises the following steps:

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

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

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

[0044] In the above formula, u p and v p are the waterlogging response values at time p and q, respectively. The two time series are aligned, a cumulative distance matrix D p,q is constructed, and is calculated one by one through the dynamic programming recursion formula, and the calculation formula is 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 , and D p-1,q-1 are the minimum values of the upper, left, and upper left cumulative distances, D p,q is the cumulative distance matrix, d(u p , v q ) is the local distance, representing the difference between the response values at time p and time q, and the calculation formula of d(u p , v q ) is as follows:

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

[0048] u p , v p are the waterlogging response values at time p, q, respectively;

[0049] S33: Calculate the minimum cumulative distance D min (U, V) by constructing the cumulative distance matrix D p,q , eliminate the time lag of the water depth or runoff curve of different response units due to different terrain and drainage capacity, so that response units with similar waterlogging response can be accurately identified. The calculation formula of the minimum cumulative distance D min (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 end 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, D P,Q represents the right bottom element in the DTW cumulative matrix D p,q , which represents the minimum cumulative cost of the path from the starting point (1, 1) to (P, Q);

[0052] S34: Take the minimum cumulative distance D min (U, V) between all pairs of response units as the matrix elements to construct the waterlogging response distance matrix D and perform normalization processing to eliminate dimensional differences to obtain the waterlogging response similarity matrix S mk , S mk The calculation formula is as follows:

[0053]

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

[0055] S35: Based on the waterlogging response similarity matrix S mk , perform clustering operation on the response units. Use the similarity threshold μ to judge any pair of response units (m, k), if it satisfies:

[0056] S mk > μ

[0057] If two response units have highly similar flood response processes, they are merged into the same cluster, which eliminates the time lag of the water depth or runoff curve in different areas due to different topography and drainage capacity, and enables accurate identification of similar flood response areas.

[0058] Further, the S4 comprises the following steps:

[0059] S41: Based on the clustering result of S35, each initial cluster T g As the input of the agglomerative hierarchical clustering, the average linkage method is used to construct the flood response distance matrix between clusters based on the minimum cumulative DTW distance of the response units The calculation formula is as follows:

[0060]

[0061] In the above formula, T g , T r represent two clusters, respectively composed of multiple response units, is the minimum cumulative DTW distance between response units g, r calculated in S34, represents the minimum distance between clusters T g , T r

[0062] S42: To avoid the limitations brought by the artificially set zoning rules, agglomerative hierarchical clustering is used to adaptively adjust the zoning process: merging the cluster pair with the closest flood response distance, and recalculating the flood response distance and average silhouette coefficient between the new cluster after merging and other clusters, repeating this step until merging into a single cluster;

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

[0064]

[0065] In the above formula, arcmax represents the parameter value of the function maximization, and backtracks to the clustering situation when the number of partitions is h*;

[0066] ​S44: Spatial adjacency constraint is performed on the clustering result obtained in S43: if a cluster is not continuous in geographical space, it is split into multiple spatially connected sub-clusters to ensure that the final partitioning result meets the physical connectivity and structural integrity of the urban geographical space. Each cluster after splitting is the final response unit. Since the flood process is jointly constrained by topography, drainage network and land use, etc., the water flow movement has significant spatial dependence. Therefore, the urban flood spatio-temporal response unit should not only consider the temporal response similarity, but also ensure that adjacent areas are spatially connected, so as to avoid the generation of "discrete" partitions that do not conform to the actual hydrological and geographical rules.

[0067] Further, the S42 comprises the following steps:

[0068] S421: Find a pair of clusters T g , T r with the minimum internal flooding response distance between clusters, and merge them into a new cluster;

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

[0070] S423: At each merging, record the current clustering partition structure, and calculate the silhouette coefficient F(h) under different partition numbers in real time, and the calculation formula is as follows:

[0071]

[0072] In the above formula, for the ith cluster, a(i) is the average distance from it to other clusters, b(i) is the average normalized distance from it to the nearest neighboring cluster, S(i) is the silhouette coefficient, and the value range is -1 to 1. The closer 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, the S5 comprises 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 of each response unit, and calculate the width WIDTH = AREA / LENGTH based on the merged water flow path length LENGTH. The calculation formula of the response unit area AREA is as follows:

[0076]

[0077] In the above formula, cluster represents the initial unit set in the response unit, AREA mis the area of the initial response unit m, and the calculation formula of the average slope SLOPE of the response unit is as follows:

[0078]

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

[0080]

[0081] In the above formula, IMPERV m is the impermeable area ratio of the initial response unit m, and the calculation formula of the impermeable surface roughness NIMPERV of the response unit is as follows:

[0082]

[0083] In the above formula, NIMPERV m is the impermeable surface roughness of the initial response unit m, and P m,cluster is the area ratio of the initial response unit in the response unit, and the calculation formula of the permeable surface roughness NPERV of the response unit is as follows:

[0084]

[0085] In the above formula, NPERV m is the permeable surface roughness of the initial response unit m, and P m,cluster is the area ratio of the initial response unit in the response unit, and the calculation formula of the width WIDTH of the response unit is as follows:

[0086]

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

[0088] S52: Import the vector data of the response unit division result into Arcmap, create AREA, SLOPE, IMPERV, NIMPERV, NPERV, and WIDTH fields, and write the calculation into the attribute table to form differentiated parameters of the response unit that can be directly used for subsequent urban waterlogging modeling;

[0089] S53: Export the response unit vector data with differentiated parameter assignment, which is used for urban waterlogging modeling in S6. Compared with the prior art, this step provides more refined underlying surface parameters considering spatial heterogeneity for subsequent surface-pipe network coupling calculation by assigning accurate differentiated response unit parameters.

[0090] Further, the S6 comprises the following steps:

[0091] S61: clean up the underground pipe network data, eliminate the detection well and sewage inspection well without actual drainage function, and check the connectivity of the pipe network topology relationship;

[0092] S62: based on the constructed rainfall grid data, the spatial matching is performed to the delineated response unit, and the process of rainfall to runoff is simulated by combining the assigned differentiated parameters, so as to generate the runoff inflow of each response unit;

[0093] S63: the runoff inflow is loaded to the pipe network node, the one-dimensional pipe network model is used to simulate the drainage process, the inflow, backflow, full and overflow process is dynamically calculated by considering the geometric and hydraulic properties and spatial distribution of the pipe network data, and the spatial structure difference and time response characteristics of the pipe network system are embodied;

[0094] S64: the overflow result output by the one-dimensional pipe network model is taken as the inflow boundary of the two-dimensional surface model, the response unit is driven to perform the ponding calculation, the one-way coupling of the one-dimensional pipe network drainage result and the two-dimensional surface ponding process is realized, and the ponding depth distribution and waterlogging evolution process of different response units are finely simulated.

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

[0096] (1) The present application introduces the theory of spatio-temporal heterogeneity into the urban waterlogging modeling by combining the law of geography, analyzes the geographical modeling mechanism contained in the urban rainfall process and response unit division process under the waterlogging scenario, and perfects the application of GIS theory in urban waterlogging modeling;

[0097] (2) The present application overcomes the limitation that the traditional modeling is difficult to comprehensively reflect the influence of spatio-temporal heterogeneity characteristics such as urban local characteristics and rainfall non-stationarity on urban flood process, and introduces multi-scale rainfall spatio-temporal grid and differentiated hydrological response unit into urban waterlogging modeling, which can dynamically reflect the rainfall variation characteristics of different regions at different times, and effectively improves the accuracy and reliability of urban waterlogging process simulation;

[0098] (3) The present application simultaneously introduces multi-dimensional indexes such as terrain characteristics, impermeable rate, land use similarity, runoff path, spatial connectivity and drainage system distribution, establishes a spatial adjacency and time sequence clustering double-constraint mechanism, realizes the true reproduction of the complex surface and pipe network system in the city, and thus improves the adaptability of the model to different city forms. BRIEF DESCRIPTION OF DRAWINGS

[0099] In order to more clearly illustrate the technical solutions in the specific embodiments or prior art of the present application, the drawings needed in the specific embodiments or prior art description will be briefly introduced as follows.

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

[0101] Figure 2 The modeling area initial response unit division result map provided by the embodiment of the application;

[0102] Figure 3 The modeling area response unit clustering result map provided by the embodiment of the application;

[0103] Figure 4 The modeling area underground drainage pipe network processing result map provided by the embodiment of the application;

[0104] Figure 5 The inflow, inundated area and inundated volume change process map of the modeling area flood response process provided by the embodiment of the application. DETAILED DESCRIPTION

[0105] In order to make the objectives, technical solutions and advantages of the present application clearer, further detailed description will be made to the present application in combination with embodiments and drawings, and the schematic embodiments of the present application and the description thereof are only used to explain the present application, and do not limit the present application.

[0106] Taking a certain region as an example, the region is a certain urban area, the overall terrain presents a "peripheral high and central low" situation, and the region contains many squares, parks, hospitals, schools and other living places, and the living supporting facilities covered are complete, including many lakes, parks and other natural water storage facilities, and the ground hardening degree is relatively high, and the road surface is mostly cement asphalt material. The rivers around the region are all built with dams, and the pipe network is relatively independent, forming a relatively independent drainage unit.

[0107] According to Figure 1 , the present application example provides a city flood modeling method considering spatio-temporal heterogeneity, which specifically includes the following steps:

[0108] Step S1: constructing a multi-scale rainfall spatio-temporal grid with spatio-temporal heterogeneity;

[0109] Step S2: data preprocessing, dividing the research region into initial response units based on terrain features, land use types and drainage system distribution;

[0110] Step S3: using dynamic time warping method to measure the similarity of waterlogging response process between different response units, and further performing time series clustering based on the division;

[0111] Step S4: using agglomerative hierarchical clustering to adjust the partition process and calculate the optimal partition number, and introducing spatial adjacency constraint to ensure spatial connectivity;

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

[0113] Step S6: driving the surface-pipe network coupling model of urban waterlogging with the multi-scale spatiotemporal grid data of rainfall and the division results of the response unit.

[0114] Further, the S1 comprises the following steps:

[0115] Step S11: based on the time axis length and the time resolution of the rainfall to be inverted, the time interval t0~t n is divided into n layers according to the time step, and each layer accurately corresponds to the spatial rainfall distribution of a specific time period, and each layer time slice satisfies:

[0116]

[0117] In the above formula, t0 represents the starting time of the rainfall process, t n represents the ending time of the rainfall process, t represents the time step of the main time slice, n represents the number of main time slices of the rainfall process, t i represents the starting time corresponding to the i-th main time slice;

[0118] Step S12: additional constraint condition of target time slice [t i， t i +t] is added when searching for rainfall observation points, and effective point rainfall data located in the modeling area Ω in the sub-time slice is screened, and the main time slice is subdivided into k sub-time slices according to the screening result, and satisfies:

[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 range of the modeling area, x i ,y j represent the coordinates of the j-th observation point, R j represents the rainfall observation value of the j-th observation point, t j represents the observation time recorded by the j-th observation point, k represents the number of subdivided sub-time slices in the current main time slice, and S i represents the set of effective rainfall observation data in the i-th main time slice and located in the modeling area Ω;

[0121] S13: based on the additional screening target [t i , t iThe rainfall point data obtained by the time slice screening of t] is spatially interpolated to obtain a rainfall grid data set in a main time slice indexed as i, and the k rainfall grid data in the main time slice are summed to accumulate rainfall to obtain rainfall grid data of the main time slice indexed as i;

[0122] S14: The above process is repeatedly performed until all n layers of time slices are calculated, and finally, the time axis is superimposed to construct a rainfall grid data set considering multiple time scales and multiple spatial resolutions, thereby providing a high-precision rainfall driving for subsequent urban flood models.

[0123] Further, the S2 comprises the following steps:

[0124] S21: Taking high-resolution DEM as the basic spatial data, the research area is divided into regular and uniform grid cells, and the cells have clear spatial topological relationships. The DEM is terrain-corrected and spatially preprocessed to remove errors;

[0125] S22: The geographic position of each rainwater well node is expressed by a vector structure, and it is regarded as an initial control node of spatial confluence. The rainwater well is taken as the core, and the spatial distance of the nearest area is determined by the Thiessen polygon method, which is preliminarily classified as the response unit of the rainwater well;

[0126] S23: The land use type of the modeling area is spatially matched with the grid cell to determine the impervious rate of each cell, and the impervious rate C m of the response unit is calculated according to the following formula:

[0127]

[0128] In the above formula, C m is the comprehensive impervious rate of the cell m, A m,l is the area of the land use type l in the cell m, α l is the impervious coefficient corresponding to the land use type l, L represents the number of contained land use types, and A m is the total area of the cell m; the impervious coefficients of different land use types are shown in Table 1:

[0129] Table 1 Land use type parameter assignment table

[0130]

[0131] S24: The direction similarity, distance similarity, and land use similarity indexes are introduced to establish a joint similarity matrix S, and the division result of the Thiessen polygon is optimized by using S. The direction similarity is mainly determined by calculating the included angle between the runoff directions of adjacent response units to determine the spatial correlation, and the calculation formula is as follows:

[0132]

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

[0134] The distance similarity is used to measure the spatial proximity of adjacent response units. The Euclidean distance is converted into the distance similarity, and the calculation formula of the Euclidean distance is as follows:

[0135]

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

[0137] The land use similarity is used to measure the similarity of the impervious rate of the underlying surface of the adjacent response units m and k, and the calculation formula is as follows:

[0138]

[0139] In the above formula, C m is the comprehensive impervious rate of the unit m, C k is the comprehensive impervious rate of the adjacent unit, S l is the land use similarity, and the calculation formula of the matrix element S mk is as follows:

[0140]

[0141] S25: Determine the flow direction and flow path of each unit through grid calculation to form a surface flow direction matrix;

[0142] S26: Modify the boundary of the vector response unit by using the obtained runoff path and flow information of the grid surface runoff expression. If the grid runoff path direction shows that the water flow direction of a certain grid is not directed to the delineated 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 the surface water flow movement, so as to optimize the accuracy of the response unit division. The finally divided initial response unit is shown in FIG. 8. Figure 2

[0143] ​Further, the S3 comprises the following steps:

[0144] S31: Extracting the waterlogging response time series of each response calculation unit under typical rainfall events, including the water depth change series, the surface runoff change series, the flood peak flow change series and the waterlogging range expansion series, for characterizing the waterlogging response process of each unit;

[0145] S32: Calculating the time series similarity between any two hydrological response units using the dynamic time warping method, and the waterlogging response time series of response units m and k are respectively:

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

[0147] In the above formula, u p , v p are the waterlogging response values at time p and q, respectively, which are aligned with two time series, and the cumulative distance matrix D p,q is constructed and calculated one by one through the dynamic programming recursion formula, and the calculation formula is 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 are the minimum values of the upper, left and upper left cumulative distances, D p,q is the cumulative distance matrix, d(u p , v q ) is the local distance, which represents the difference between the response values at time p and time q, and the calculation formula of d(u p , v q ) is as follows:

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

[0151] In the above formula, u p , v p are the waterlogging response values at time p and q, respectively;

[0152] S33: Obtaining the minimum cumulative distance D min (U, V) through the constructed cumulative distance matrix D p,q, eliminate the different response unit of the water depth or runoff curve may due to the different time lag in the terrain, drainage capacity, make the response unit with similar waterlogging response can be accurately identified, the minimum cumulative distance D min The calculation formula of (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 beginning of sequence U to the end of sequence V, P and Q represent the length of time series X and Y respectively, that is, the number of elements in U and V, D P,Q represents the DTW cumulative matrix D p,q The right lower corner element in the matrix represents the minimum cumulative cost of the path from the starting point (1, 1) to (P, Q);

[0155] S34: Take the minimum cumulative distance D min (U, V) between all response units as matrix elements, and construct the waterlogging response distance matrix And carry out normalization processing to eliminate the dimensional difference, and obtain the waterlogging response similarity matrix S mk , S mk The calculation formula is as follows:

[0156]

[0157] In the above formula, max(D min ) is the maximum cumulative distance value between all response unit pairs, and S mk The closer to 1 indicates that the similarity of the waterlogging response process is higher;

[0158] S35: Based on the waterlogging response similarity matrix S mk , the clustering operation is carried out on the response unit, and the similarity threshold μ is used to judge any response unit pair (m, k), if it satisfies:

[0159] S mk > μ

[0160] It is considered that the two response units have highly similar waterlogging response process, and they are merged into the same cluster, and finally each merged cluster is regarded as the initial cluster.

[0161] Further, the S4 comprises the following steps:

[0162] S41: On the basis of the clustering result of S35, each initial cluster T g is taken as the input of the agglomerative hierarchical clustering, and the average connection method is adopted to construct the waterlogging response distance matrix between clusters based on the minimum cumulative DTW distance of the response unit The calculation formula is as follows:

[0163]

[0164] In the above formula, T g , T r represent two clusters, respectively composed of multiple response units, is the minimum cumulative DTW distance between response units g, r calculated in S34, represents the minimum distance between clusters T g , T r ;

[0165] S42: Adopting agglomerative hierarchical clustering, the partitioning process is adaptively adjusted: merging the pair of clusters with the closest response distance, and recalculating the response distance between the new cluster after merging and other clusters and the average profile coefficient, repeating this step until merging into a single cluster;

[0166] S43: Traversing the agglomerative clustering tree, selecting the number of clusters with the largest average profile coefficient as the optimal partition number h* on the time scale, the calculation formula is as follows:

[0167]

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

[0169] S44: Spatial adjacency constraint is performed on the clustering results obtained in S43: if a cluster is not continuous in geographical space, it is split into multiple connected clusters to ensure that the final partitioning result meets the physical connectivity of urban geographical space, and each cluster after splitting is the final response unit, and the final response unit clustering result is as shown in Figure 3 .

[0170] Further, the S42 comprises the following steps:

[0171] S421: Find a pair of clusters T g , T r with the minimum response distance between clusters, and merge them into a new cluster;

[0172] S422: After merging, the average connection method is used to update the response distance between the new cluster after merging and other clusters;

[0173] S423: At each merging, record the current clustering division structure, and calculate the profile coefficient F(h) under different partition numbers in real time, the calculation formula is as follows:

[0174]

[0175] In the above formula, for the βth cluster, z(β) is the average distance from it to other clusters, w(β) is the average normalized distance from it to the nearest neighboring cluster, F(β) is the silhouette coefficient, and the value ranges from -1 to 1, and the closer 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, the S5 comprises 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 of each response unit, and calculate the width WIDTH = AREA / LENGTH based on the length LENGTH of the merged water flow path. The calculation formula of the response unit area AREA is as follows:

[0179]

[0180] In the above formula, cluster represents the initial unit set in the response unit, AREA m is the area of the initial response unit m, and the calculation formula of the response unit average slope SLOPE is as follows:

[0181]

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

[0183]

[0184] In the above formula, IMPERV m is the impermeable area ratio of the initial response unit m, and the calculation formula of the response unit impermeable surface roughness NIMPERV is as follows:

[0185]

[0186] In the above formula, NIMPERV m is the impermeable surface roughness of the initial response unit m, and P m,cluster is the area ratio of the initial response unit in the response unit, and the calculation formula of the response unit permeable surface roughness NPERV is as follows:

[0187]

[0188] In the above formula, NPERV m is the permeable surface roughness of the initial response unit m, and Pm,cluster The initial response unit is the area ratio of the response unit, and the response unit width WIDTH is calculated as follows:

[0189]

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

[0191] S52: Import the vector data of the response unit division result into Arcmap, create AREA, SLOPE, IMPERV, NIMPERV, NPERV, WIDTH fields, and write the calculation to the corresponding attribute table, forming differentiated parameters of the response unit that can be directly used for subsequent urban waterlogging modeling;

[0192] S53: Export the response unit vector data with differentiated parameter assignment completed, which is used for urban waterlogging modeling in S6;

[0193] Further, the S6 includes the following steps:

[0194] S61: Clean the underground pipe network data, the present application adopts breadth-first algorithm to traverse the rainwater well and rainwater grate nodes, delete the exploration well, sewage well and isolated nodes not traversed, retain the main rainwater well nodes, and eliminate part of the secondary pipe network and other functional pipe networks. After data cleaning, finally 1600 rainwater well nodes, 1607 rainwater pipe sections and 7 drainage outlets are obtained. According to the node ID and pipe section inlet and outlet node ID attribute fields, the pipe network connectivity is checked by using ArcGIS, and part of the pipe section information with flow direction conflict is corrected. The final pipe network processing result is as shown in Figure 4 ;

[0195] S62: Based on the constructed rainfall grid data, match it to the delineated response unit in space, and simulate the process of rainfall to runoff in combination with the differentiated parameters, to generate the runoff inflow of each response unit;

[0196] S63: Load the runoff inflow to the pipe network node, adopt one-dimensional pipe network model to simulate the drainage process, consider the pipe network data geometry and hydraulic property, spatial distribution, dynamically calculate the inflow, backflow, full and overflow process, and embody the spatial structure difference and time response characteristics of the pipe network system;

[0197] S64: Take the overflow result output by the one-dimensional pipe network model as the inflow boundary of the two-dimensional surface model, drive the response unit to perform ponding calculation, realize one-way coupling of the one-dimensional pipe network drainage result and the two-dimensional surface ponding process, finely simulate the ponding depth distribution and waterlogging evolution process of different response units, and the inflow, flooded area and flooded volume change process of the flood response process is as shown in Figure 5 .

[0198] The above description of the embodiments is only used to help understand the method of the present application and its core idea; meanwhile, for those skilled in the art, according to the idea of the present application, there will be changes in the specific implementation and application range, and the above description should not be understood as a limitation on the present application.

Claims

1. A method for modeling urban flooding that takes into account temporal and spatial heterogeneity, characterized in that: The following steps are involved: Step S1: construct a multi-scale rainfall spatiotemporal grid with spatiotemporal heterogeneity; Step S2: Data preprocessing, dividing the study area into initial response units based on terrain characteristics, land use types and drainage system distribution; Step S3: Using the dynamic time warping method to measure the similarity of the waterlogging response process between each response unit, and based on this, perform time series clustering on the division results; Step S4: Use agglomerative hierarchical clustering method to further cluster and determine the optimal number of partitions, and introduce spatial adjacency constraints to ensure spatial connectivity; Step S5: assigning differentiated hydrological parameters to different response units; Step S6: Use the multi-scale rainfall spatiotemporal grid data and the response unit division results to drive the modeling of the surface-pipeline network coupling model of urban waterlogging.

2. The urban waterlogging modeling method according to claim 1, characterized in that: The 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 to t n The time slice is divided into n layers according to the time step. Each layer accurately corresponds to the spatial rainfall distribution of a specific time period. The time slice of each layer satisfies: In the above formula, t0 represents the starting time of the rainfall process, t n Indicates the end time of the rainfall process, t indicates the time step of the main time slice, n indicates the number of main time slices of the rainfall process, t i Indicates the start time corresponding to the i-th main time slice; Step S12: Add target time slice [t i， t i +t], filter the valid point rainfall data in the modeling area Ω within the sub-time slice, and subdivide the main time slice into k sub-time slices based on the screening results, specifically satisfying: S i ={(x j ,y j ,R j )|(x j ,y j )∈Ω,t j ∈[t i ,t i +t]} In the above formula, Ω represents the spatial range of the modeling area, x i ,y j represents the coordinates of the jth observation point, R j represents the rainfall observation value at the jth observation point, t j represents the observation time recorded at the jth observation point, k represents the number of sub-time slices subdivided in the current main time slice, S i represents the set of valid rainfall observation data within the i-th primary time slice and located in the modeling area Ω; Step S13: Perform spatial interpolation calculation on the rainfall point data obtained by filtering the time slices of the additional filtering target to generate a rainfall grid dataset within the main time slice indexed by i, and sum the rainfall grid data within the main time slice to accumulate the rainfall amount to obtain the rainfall grid data of the main time slice indexed by i; Step S14: The above process is repeated until all n layers of time slices are calculated. Finally, they are superimposed 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 for subsequent urban waterlogging models.

3. The urban flooding modeling method according to claim 1, characterized in that: The S2 comprises the following steps: Step S21: Using the 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 the cells, and the DEM is subjected to terrain correction and spatial preprocessing to remove errors; Step S22: Using a vector structure to express the geographical location of each rainwater well node, treating it as the initial control node of spatial convergence, using the Thiessen polygon method to determine the closest area in spatial distance with each rainwater well as the core, and preliminarily classifying it as the response unit of the rainwater well; Step S23: Spatially match the land use type of the modeling area with the grid cells, determine the imperviousness of each cell, and calculate the comprehensive imperviousness of the response cell C. m The calculation formula is as follows: In the above formula, C m is the comprehensive impermeability of unit m, A m,l is the area of ​​land use type l in unit m, α l is the impervious coefficient corresponding to land use type l, L represents the number of land use types included, A m is the total area of ​​unit m; Step S24: further introduce directional similarity, distance similarity, and land use similarity indicators to establish a joint similarity matrix S, and use S to optimize the division results of the Thiessen polygon. Directional similarity mainly determines spatial correlation by calculating the angle between the runoff directions of adjacent response units. The calculation formula is as follows: In the above formula, θ is the angle between the centroid of the response unit m and its adjacent response unit k and its own runoff direction, in radians; S θ is the direction similarity; Distance similarity is used to measure the spatial proximity of adjacent response units. The Euclidean distance is converted into distance similarity. The calculation formula of the Euclidean distance is as follows: In the above formula, (x m ,y m ) and (x k ,y k ) represent the geographical coordinates of the centroid of response units m and k respectively; E mk is (x m ,y m ) and (x k ,y k ) The Euclidean distance between two points, S d is the distance similarity; Land use similarity is used to measure the similarity of the underlying surface imperviousness of adjacent response units m, k. The calculation formula is as follows: In the above formula, C m is the comprehensive impermeability of unit m, C k is the comprehensive impermeability of adjacent units, S l is the land use similarity, the matrix element S mk The calculation formula is as follows: Step S25: Determine the flow direction and flow path of each unit through grid calculation to form a surface flow direction matrix, which provides a basis for the subsequent fine division of response units; Step S26: Using the runoff path and flow information obtained from the raster surface runoff expression, the boundaries of the vector response units are corrected. If the raster runoff path direction shows that the water flow direction of a certain raster is not pointing to the designated response unit but to other response units, the boundaries of the vector response units should be adjusted to ensure that the final spatial division conforms to the actual path of surface water flow movement, thereby optimizing the accuracy of the response unit division.

4. The urban flooding modeling method taking into account temporal and spatial heterogeneity according to claim 1 is characterized in that: The S3 comprises the following steps: Step S31: extracting the waterlogging response time series of each response calculation unit under a typical rainfall event, including the waterlogging 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: Dynamic time warping is used to calculate the temporal similarity between any two hydrological response units. The waterlogging response time series of response units m and k are: U=(u1,u2,…u b ),V=(v1,v,…v c ) In the above formula, u p 、v p They are the waterlogging response values ​​at time p and q respectively. The two are aligned to construct the cumulative distance matrix D p,q , and calculate them one by one through the dynamic programming recursive formula. The calculation formula is as follows: D p,q =d(x p ,y q )+min{D p-1,q ,D p,q-1 ,D p-1,q-1 } In the above formula, D p-1,q ,D p,q-1 ,D p-1,q-1 is the minimum value of the cumulative distance between the top, left, and top left, D p,q is the cumulative distance matrix, d(u p ,v q ) is the local distance, which represents the difference between the response values ​​at time p and time q, d(u p ,v q ) is calculated as follows: d(u p ,v q )=|u p -v q | In the above formula, u p 、v p are the waterlogging response values ​​at time p and q respectively; Step S33: Calculate the minimum cumulative distance D min (U, V), by constructing the cumulative distance matrix D p,q , eliminating the time lag of water depth or runoff curves of different response units due to differences in terrain and drainage capacity, so that response units with similar waterlogging responses can be accurately identified, and the minimum cumulative distance D min The calculation formula for (U,V) is as follows: D min (U,V)=D P,Q 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 end point of sequence V, P and Q represent the length of time series X and Y, that is, the number of elements in U and V, respectively. P,Q Denotes the DTW cumulative matrix D p,q The element in the lower right corner represents the minimum cumulative cost of the path from the starting point (1,1) to (P,Q); Step S34: The minimum cumulative distance D between all response units is calculated. min (U, V) is used as a matrix element to construct the waterlogging response distance matrix And normalization is performed to eliminate the dimension difference, and the waterlogging response similarity matrix S is obtained. mk , S mk The calculation formula is as follows: In the above formula, max(D min ) is the maximum cumulative distance value between all response unit pairs, S mk The closer it is to 1, the higher the similarity of the waterlogging response process; Step S35: Based on the waterlogging response similarity matrix S mk , cluster the response units and use the similarity threshold μ to judge any response unit pair (m, k). If it satisfies: S mk >m The two response units are considered to have highly similar waterlogging response processes and are merged into the same cluster. Each cluster that is finally merged is regarded as the initial cluster.

5. The urban flooding modeling method according to claim 1, characterized in that: The S4 comprises the following steps: Step S41: Each initial cluster T g As the input of agglomerative hierarchical clustering, the average connection method is used to construct the waterlogging response distance matrix between clusters based on the minimum cumulative DTW distance of the response unit. The calculation formula is as follows: In the above formula, T g 、T r Represents two clusters, each consisting of multiple response units. is the minimum cumulative DTW distance between the response units g and r calculated in S34, Represents cluster T g 、T r The minimum distance between Step S42: Adopting agglomerative hierarchical clustering to make the partitioning process adaptively adjusted: merge the cluster pairs with the closest waterlogging response distance, and recalculate the waterlogging response distance and average silhouette coefficient between the merged new cluster and other clusters. Repeat this step until a single cluster is merged. Step S43: traverse the agglomerative clustering tree and select the number of clusters with the largest average silhouette coefficient as the optimal partition number h* on the time scale. The calculation formula is as follows: In the above formula, arcmax represents the parameter value of the function that maximizes , and traces back 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 discontinuous in geographic space, split it into multiple connected clusters to ensure that the final partitioning result meets the physical connectivity of the urban geographic space. Each cluster after the splitting is a divided response unit.

6. The method of dynamically adjusting the partitioning results by combining the spatiotemporal clustering analysis with the partitioning algorithm according to claim 5 is characterized in that: The S42 includes the following steps: Step S421: Find a pair of clusters T with the smallest inter-cluster waterlogging response distance g 、T r Merge them into a new cluster; Step S422: after merging, the average connection method is used to update the waterlogging response distance between the merged new cluster and other clusters; Step S423: During each merging, the current clustering structure is recorded, and the silhouette coefficient F(h) under different numbers of partitions is calculated in real time. The calculation formula is as follows: In the above formula, for the βth cluster, z(β) is the average distance between it and other clusters, w(β) is the average normalized distance between it and the nearest neighboring 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 operations to gradually generate a complete agglomerative clustering tree from the initial N clusters to a single cluster.

7. The urban flooding modeling method taking into account temporal and spatial heterogeneity according to claim 1 is characterized in that: The S5 comprises the following steps: Step S51: Count the area AREA, average slope SLOPE, impervious area ratio IMPERV, impervious surface roughness NIMPERV, and pervious surface roughness NPERV of each response unit, and calculate the width WIDTH = AREA / LENGTH based on the combined water flow path length LENGTH. The calculation formula for the response unit area AREA is as follows: In the above formula, cluster represents the initial unit set in the response unit, AREA m is the area of ​​the initial response unit m, and the calculation formula of the average slope SLOPE of the response unit is as follows: In the above formula, SLOPE m is the slope of the initial response unit m, and the impervious area ratio IMPERV of the response unit is calculated as follows: In the above formula, IMPERV m is the impervious area ratio of the initial response unit m, and the calculation formula of the impervious surface roughness NIMPERV of the response unit is as follows: In the above formula, NIMPERV m is the impervious surface roughness of the initial response unit m, P m,cluster is the area ratio of the initial response unit to the response unit. The calculation formula of the permeable surface roughness NPERV of the response unit is as follows: In the above formula, NPERV m is the permeable surface roughness of the initial response unit m, P m,cluster The area ratio of the initial response unit to the response unit is calculated as follows: In the above formula, LENGTH is the water flow path length of the response unit; Step S52: Import the vector data of the response unit division result into Arcmap, create new AREA, SLOPE, IMPERV, NIMPERV, NPERV, and WIDTH fields, and write the calculated results into the attribute table to form the response unit differentiation parameters that can be directly used for subsequent urban waterlogging modeling; Step S53: exporting the response unit vector data after the differentiated parameter assignment is completed, for use in the urban waterlogging modeling in S6.

8. The urban flooding modeling method taking into account temporal and spatial heterogeneity according to claim 1 is characterized in that: The S6 comprises 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 delineated response units, and combine the assigned differentiation parameters to simulate the rainfall to runoff process and generate the runoff inflow of each response unit; Step S63: Load the runoff inflow to the pipe network nodes and use a one-dimensional pipe network model to simulate the drainage process. Considering the pipe network data geometry, hydraulic properties, and spatial distribution, the inflow, backflow, filling, and overflow processes are dynamically calculated to reflect the spatial structure differences and time response characteristics of the pipe network system. Step S64: Use the overflow result 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 of the one-dimensional pipe network drainage result and the two-dimensional surface water accumulation process, and accurately simulate the water accumulation depth distribution and waterlogging evolution process of different response units.

9. An electronic device comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein: When the processor executes the program, the list-driven intelligent data collection method according to any one of claims 1 to 7 is implemented.

10. A computer-readable storage medium having computer instructions stored thereon, characterized in that: When the computer instructions are executed by a processor, the list-driven intelligent data collection method according to any one of claims 1 to 7 is implemented.

Images