Ground surface accumulated water rasterization calculation method for pipe network data missing region

Abstract

The invention discloses an earth surface accumulated water rasterization calculation method for a pipe network data missing region, which considers an earth surface runoff and confluence process, and adopts a rainfall intensity grid and drainage capacity calculated by design rainstorm intensity in a pipe duct design return period to participate in calculation of earth surface direct accumulated water quantity. According to the method, the calculation simulation of the ponding amount under the condition of no pipe network data is realized, and the earth surface ponding rebalance calculation method is realized by applying a D8 algorithm and calculating the confluence cumulant: the dynamic rebalance process of the ponding overflow process is considered and calculated, and the ponding simulation result precision is improved.

Description

Technical Field

[0001] The invention relates to the field of rainfall processing, and in particular to a rasterized calculation method for surface water accumulation in areas where pipe network data is missing. Background Art

[0002] Against the backdrop of global climate change, my country's extreme rainfall events are increasing. my country is located in a monsoon region and is strongly affected by the monsoon, with concentrated rainstorms and floods and severe flood disasters. With the continuous development of urbanization, urban flood disasters have become one of the main natural disasters affecting the economic and social development of urban areas. In the study of urban flood disasters, the simulation of waterlogging processes in most areas requires the participation of pipeline network data, but urban pipeline network data is often difficult to obtain. For example, it is lost due to age, or its original pipeline network data is not saved at the beginning; even if there is a pipeline network design drawing, the actual layout will be inconsistent with the design drawing due to the involvement of many departments and coordination difficulties during the construction period; and pipeline network data in some areas cannot be made public due to confidentiality issues.

[0003] In summary, a rasterized calculation method for surface water in areas where pipe network data is missing is needed to solve the shortcomings of the existing technology. Summary of the invention

[0004] In view of the deficiencies in the prior art, the present invention provides a rasterized calculation method for surface waterlogging in areas where pipe network data is missing, aiming to solve the above problems.

[0005] To achieve the above object, the present invention provides the following technical solution: a rasterized calculation method for surface waterlogging in areas where pipe network data is missing, comprising the following steps: Step S1: Obtain rainfall data of the weather station that changes over time and calculate rainfall intensity data; Step S2: Calculate the amount of water accumulation, calculate the rainfall spatial distribution grid and the drainage capacity grid, and calculate the direct surface water accumulation based on the rainfall intensity data and the drainage capacity; Step S3: Calculate the depression depth spatial grid by using hydrological analysis and regional analysis to obtain the depression depth; Step S4: Determine the relationship between the amount of direct surface water accumulation and the depth of the depression based on the results of steps S2 and S3; Step S5: Rebalance the surface water volume, use the D8 algorithm to calculate the direction of the overflow water flow, calculate the cumulative amount of confluence, and perform rebalancing.

[0006] Optionally, the rainfall intensity data is calculated in step S1 in the following manner: I=P t / t, where I is the rainfall intensity, P t is the accumulated rainfall in time t, where t is the rainfall time.

[0007] Optionally, the rainfall spatial distribution grid in step S2 is calculated by: The inverse distance weighted interpolation method is used to calculate the spatial distribution grid of rainfall intensity. The calculation formula is as follows: , , where is the rainfall intensity value of the interpolation point S; N is the number of meteorological stations used to interpolate at S, For the i The rainfall intensity at the ith weather station at λ; i is the weight of the ith weather station; d i is the distance from the interpolation point to the i-th weather station; b is the weight index.

[0008] Optionally, the drainage capacity grid is calculated by the following steps: Step A1: Determine the drainage capacity of the region and use rainfall data to analyze the relationship between rainfall intensity I, rainfall duration t and return period P; Step A2: Calculation of direct surface water accumulation, the calculation formula is: △I=I in -I out , where △I is the direct surface water accumulation per hour, I in is the spatial grid value of rainfall intensity, I out It is the drainage capacity spatial grid value.

[0009] Optionally, the relationship formula in step A1 is: , where I is rainfall intensity, t is rainfall duration, P is return period, A1 is the rainfall force of the design rainfall with a return period of 1 year, C is the rainfall force variation parameter, b is the rainfall duration correction parameter, and n is the rainfall attenuation index.

[0010] Optionally, the depression depth spatial grid calculation in step S3 is performed in the following manner: Step B1: Depression extraction, obtaining boundary data of the depression area to be extracted and digital elevation model data to analyze the water flow direction, and obtaining water flow direction raster data; Step B2: Identify and extract depression data, traverse the water flow grid data in a set order, and identify the depression grid data; Step B3: Overlay analysis is performed on the water flow direction grid data and the depression grid data to obtain the depression contribution area; Step B4: Overlay analysis of each depression contribution area with the digital elevation model data to calculate the depression depth.

[0011] Optionally, the water flow direction is analyzed in the following manner: The digital elevation model is clipped according to the boundary data to obtain the DEM data of the specific area, which is then rasterized. A starting grid unit in the depression detection area after removing the boundary is selected, and a matrix is ​​constructed with the grid as the center. The D8 algorithm is used to analyze the water flow direction.

[0012] Optionally, the D8 algorithm is used to analyze the water flow direction as follows: The direction of water flow is determined by the steepest slope method, and the distance difference between the central grid and the adjacent grid is calculated. The calculation formula is: , Where Drop is the maximum distance weighted drop between the central grid and its adjacent grid cells, H is the elevation value of the DEM data center grid cell, Hi is the elevation value of the adjacent grid cell, and i is the flow direction code relative to the central grid cell.

[0013] Optionally, the depression depth is calculated in the following manner: Traverse the digital elevation model data grid cells in each depression contribution area and obtain its lowest elevation grid value Sink min , and then traverse the digital elevation model data grid cells of each depression contribution area boundary, extract the minimum elevation value of the regional boundary grid cell, that is, get the lowest elevation value of the depression contribution area boundary Sink max , Depth of Sink Depth =Sink max —Sink min .

[0014] Optionally, the determination in step S4 is performed in the following manner: If the direct surface water accumulation is less than or equal to the depression depth, the required surface water accumulation depth is equal to the direct surface water accumulation. Otherwise, the surface water overflows and a rebalance calculation of the surface water accumulation is required.

[0015] Beneficial effects of the present invention: 1. In the present invention, the spatial distribution grid of rainfall intensity is obtained by combining the rainfall data of the meteorological station with the inverse distance weighted interpolation algorithm, and the urban drainage capacity grid is determined according to the outdoor drainage standard. The difference between the two is used to obtain the direct surface water accumulation. The GIS analysis methods such as hydrological analysis and regional analysis are used to extract depressions, calculate the depression contribution area and depression depth, and for overflow water whose water accumulation is greater than the depression depth, the D8 algorithm is used to calculate the water flow direction, and based on this, the confluence accumulation is calculated to achieve surface water rebalancing and obtain the final water accumulation depth result; 2. In the present invention, the surface runoff generation and convergence process is taken into account, and the drainage capacity calculated by using the rainfall intensity grid and the design rainstorm intensity of the pipe design recurrence period is used to participate in the calculation of the direct surface water accumulation, thus realizing the calculation simulation of the water accumulation without pipe network data; 3. In the present invention, the D8 algorithm and the calculation of the accumulated runoff are used to realize the surface water rebalancing calculation method: this method takes into account and calculates the dynamic rebalancing process of the water overflow process, thereby improving the accuracy of the water accumulation simulation results. BRIEF DESCRIPTION OF THE DRAWINGS

[0016] Figure 1 The present invention is a schematic flow chart of a method.

[0017] Figure 2 This is a schematic diagram of the four-category urban area division according to the present invention.

[0018] Figure 3 This is a depression depth spatial grid calculation graph of the present invention.

[0019] Figure 4 This is a schematic diagram of water flow direction coding according to the present invention.

[0020] Figure 5 This is a schematic diagram of depression data extraction according to the present invention.

[0021] Figure 6 This is a schematic diagram of a depression contribution area of ​​the present invention.

[0022] Figure 7 It is a schematic diagram of water accumulation, overflow and confluence accumulation according to the present invention. DETAILED DESCRIPTION

[0023] In order to more clearly illustrate the embodiments of the invention or the technical solutions in the prior art, the drawings required for use in the embodiments or the description of the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying any creative work.

[0024] like Figures 1 to 7 As shown in the figure, a rasterized calculation method for surface water in areas where pipe network data is missing includes the following contents: 1. Calculation of direct surface water accumulation 1.1 Calculating the spatial distribution grid of rainfall 1.1.1 Obtain rainfall data from weather stations over time and calculate rainfall intensity Rainfall intensity (i.e., rain strength) usually refers to the amount of rainfall per unit time. After obtaining the rainfall data for a specific period of time, the rainfall intensity is calculated using the following formula: I=P t / t, where I is the rainfall intensity in millimeters per hour (mm / h); Pt is the accumulated rainfall in a period of time t in millimeters (mm); and t is the rainfall duration in hours (h).

[0025] 1.1.2 Use the inverse distance weighted interpolation method to calculate the spatial distribution grid of rainfall intensity.

[0026] The inverse distance weighted interpolation method (IDW) is a commonly used spatial interpolation method. The weight is determined according to the distance to the interpolation point for interpolation. The weight is inversely proportional to the distance. The farther away from the interpolation point, the smaller the weight is. The calculation formula is as follows: , , In the formula, is the rainfall intensity value of the interpolation point S; N is the number of meteorological stations used to interpolate at S, For the i The rainfall intensity at the ith weather station at λ; i is the weight of the ith weather station; d i is the distance from the interpolation point to the i-th meteorological station; b is the weight index, and as the commonly used inverse distance method, b=2 here.

[0027] Based on the rainfall intensity I of each meteorological station calculated in 1.1.1, the grid data value I is further interpolated according to the inverse distance weighted interpolation formula. * That is the local rainfall intensity spatial distribution grid value I in .

[0028] 1.2 Calculation of drainage capacity spatial grid 1.2.1 Classification of four types of urban areas According to the "Outdoor Drainage Design Standard" and "Urban Area Determination Regulations", the urban area is divided into four types: non-central urban area, central urban area, important areas of central urban area and central urban area flood-prone points (underpasses, underground passages and sunken squares, etc.), such as Figure 2 shown.

[0029] 1.2.2 Determining the drainage capacity of different types of urban areas At present, my country has accumulated complete automatic rainfall record data, and the method of mathematical statistics can be used to calculate and determine the rainstorm intensity formula. my country's "Outdoor Drainage Design Standard" stipulates that areas with more than 20 years of automatic rainfall records should use the annual maximum value method to compile rainfall intensity formulas, and areas with less than 20 years of automatic rainfall records can use the annual multiple sample method. The longer the statistical data, the more the compiled rainstorm intensity formula can reflect the local rainfall intensity law. The relationship between the rainstorm intensity I, rainfall duration t and recurrence period P is obtained by analyzing the rainfall data, and the values ​​of each parameter in the designed rainstorm intensity formula are calculated based on it.

[0030] , where I is rainfall intensity, t is rainfall duration, P is return period, A1 is the rainfall force of the design rainfall with a return period of 1 year, C is the rainfall force variation parameter, b is the rainfall duration correction parameter, and n is the rainfall attenuation index. A1, C, b, and n are all local parameters determined using statistical methods.

[0031] According to the "Outdoor Drainage Design Standard", the design return period of rainwater pipes is shown in Table 1:

[0032] Among them, towns with dense populations, prone to waterlogging and good economic conditions should adopt the prescribed upper limit of the design recurrence period. Therefore, taking large cities as an example, the recurrence period of non-central urban areas is once in 3 years, the recurrence period of central urban areas is once in 5 years, the recurrence period of important areas in central urban areas is once in 10 years, and the recurrence period of flood-prone points in central urban areas such as tunnels, underground passages and sunken squares is once in 30 years. According to the recurrence period, the design rainfall intensity i is calculated by the local rainstorm intensity formula, which is the drainage capacity value I of each grid of different urban types. out .

[0033] 1.3 Calculation of direct surface water accumulation Without considering the influence of terrain, the calculation method of direct surface water accumulation is: △I=I in -I out , where △I is the direct surface water accumulation per hour, I in is the spatial grid value of rainfall intensity, I out It is the drainage capacity spatial grid value.

[0034] That is, the hourly direct surface water accumulation is the difference between the spatial grid value of rainfall intensity and the spatial grid value of drainage capacity.

[0035] 2. Depression depth spatial grid calculation The calculation process of the depression depth spatial grid is as follows: Figure 3 shown.

[0036] 2.1 Depression extraction 2.1.1 Analysis of water flow direction Obtain the boundary data and digital elevation model (DEM) data of the depression area to be extracted. Clip the digital elevation model (DEM) data according to the boundary data to obtain the digital elevation model (DEM) data of the area, and rasterize it. Use the local digital elevation model (DEM) data to remove the boundary grid cells as the depression detection area, and select the grid cell in the upper left corner of the detection area as the detected grid cell. In the 3×3 grid matrix centered on the detected grid cell, the D8 algorithm is used to calculate the flow direction of water on the DEM surface. The D8 algorithm assumes that the water flow in a single grid can only flow into the 8 adjacent grid cells, and the steepest slope method is used to determine the direction of the water flow, that is, on the 3×3 grid matrix, calculate the distance weighted drop between the central grid and its adjacent grid cells, that is, the drop of the grid center point divided by the distance between the grid center points. The calculation formula is as follows: , In the formula, Drop is the maximum distance weighted drop between the central grid and its adjacent grid cells, H is the elevation value of the DEM data center grid cell, Hi is the elevation value of the adjacent grid cell, and i is the flow direction code relative to the central grid cell, such as Figure 4 As shown. d is the distance between the center point of the central grid unit and the center point of the adjacent grid unit. When i=1,4,16,64, d is the size of the grid unit; when i=2,8,32,128, d is times the grid cell size.

[0037] If there is a maximum distance weight drop, the grid with the largest distance weight drop is taken as the outflow grid of the center grid, and the vector direction between the outflow grid and the detected grid unit is the flow direction. The outflow grid is encoded according to the flow direction, such as Figure 4 As shown. The coded value of the outflow grid is assigned to the detected grid cell. If the distance weight difference between the detected grid cell and the adjacent pixel is 0, the value of the detected grid cell is 0.

[0038] Traverse the digital elevation model grid cells of the detection area from left to right and from top to bottom, assign values ​​to each grid cell according to the above steps, and obtain the flow direction grid data displayed in grayscale color according to the flow direction value.

[0039] 2.1.2 Depression identification and extraction Traverse the flow direction grid data from left to right and from top to bottom, and extract the grid cells with a value of 0. The extracted grid values ​​are judged, and if one of the following three conditions is met, it is determined to be a depression grid: ① The values ​​of the detected grid cell and its 8 adjacent grid cells are all 0; ② The positions of the 8 adjacent grid cells relative to the detected grid cell are opposite to the flow direction, that is, the water flow direction points to the detected grid cell; ③ The positions of the 8 adjacent grid cells relative to the detected grid cell are opposite to the flow direction and the values ​​of the detected grid cell and its 8 adjacent grid cells are all 0. In this way, the identified depression grids that meet the conditions are extracted to obtain the depression grid data. An example of extracting depression raster data is as follows: Figure 5 As shown, the depression grid cells are determined and extracted according to the flow direction grid cell values ​​and the water flow direction.

[0040] 2.2 Calculating the Contribution Area of ​​Depression The obtained flow direction grid data and depression grid data are overlaid and analyzed. For each independent depression in the depression grid data, an attribute value different from other depressions is assigned to the depression grid. That is, assuming that there are five independent depressions in the depression grid data, the attribute value "1" is assigned to depression 1, the attribute value "2" is assigned to depression 2, and so on.

[0041] The depression raster data is used as seed data, and the regional growing algorithm is used for the flow direction raster data. The algorithm requires that "the position of the adjacent raster cells relative to the detected raster cells is opposite to the flow direction", and all raster cells that meet this requirement are assigned the same attribute value as the seed data. The calculation ends when a raster cell that does not meet this condition is encountered. Since the attribute values ​​of each independent depression used as seed data are different, the flow direction data is formed into multiple blocks with different attributes through the regional growing algorithm. Each block is the contribution area of ​​the corresponding independent depression. The collection of all blocks is the depression contribution area raster data. An example of depression contribution area extraction is shown below. Figure 6 As shown: 2.3 Calculating the Depth of Depression For each depression contribution area, it is overlaid with the local digital elevation model (DEM) data for analysis, and the digital elevation model data grid cells in each depression contribution area are traversed to obtain its lowest elevation grid value Sink min Then, the digital elevation model data grid cells of each depression contribution area boundary are traversed to extract the minimum elevation value of the regional boundary grid cell, that is, the lowest elevation value of the depression contribution area boundary Sink is obtained. max Finally, the depth of the depression can be calculated by the following formula: Sink Depth =Sink max —Sink min

[0042] That is, the depth of the depression is the difference between the lowest elevation value of the exit of the depression contribution area boundary and the lowest elevation grid value in the area.

[0043] Determine the relationship between the direct water accumulation volume △I and the depth of the depression Depth If △I is less than or equal to Depth, then the final required surface water depth result I Depth =△I; if △I is greater than Depth, it means that the surface water has overflowed and a rebalance calculation of the surface water volume is required.

[0044] 4. Rebalance the surface water volume.

[0045] 4.1.1 Use the D8 algorithm to calculate the direction of overflow water flow.

[0046] 4.1.2 Calculation of cumulative flow like Figure 7 As shown in the figure, each grid of the DEM digital elevation model has a unit of water volume, that is, it is assumed that each grid has an equal and symbolic amount of water at the beginning. Usually this unit can be arbitrary, such as "1" or other convenient values ​​for calculation. According to the natural law that natural water flows from high to low, the water flow direction digital matrix of the regional terrain calculated based on the D8 algorithm is used to obtain the water volume flowing through each point. The difference between the amount of accumulated water and the depth of the depression is taken as the overflow volume. The multiple unit water volumes collected by the grid unit are taken as weight parameters, that is, the proportion of the overflow volume flowing to the grid unit, multiplied by the total overflow volume of the original grid unit, which is the cumulative flow of the area, that is, the final water depth result.

[0047] Urban stormwater runoff model: A mathematical model that is established based on the laws of urban rainfall runoff and using modern hydrological and hydrodynamic knowledge to simulate stormwater runoff and drainage systems, and further analyze water flow levels and flooding conditions.

[0048] Rainfall intensity: the amount of rainfall per unit time.

[0049] Raster Data: It is a format used to represent spatial data. It divides the earth's surface into a regular grid (i.e., grid cells). Each grid cell (also called pixel or picture element) contains a specific value, which can represent the surface's elevation, temperature, humidity, vegetation coverage and other properties.

[0050] Inverse distance weighted interpolation (IDW): Also known as inverse distance weighted interpolation, it is a weighted average interpolation method. In this method, the point to be interpolated is assigned the product of a specific data point value and its weight, and the weight is inversely proportional to the distance between the known data point and the point to be interpolated.

[0051] Recurrence period: The average interval between the occurrence of rainfall greater than or equal to a certain rainstorm intensity during the statistical period of rainfall record data in a certain era, which is the inverse of the frequency of the rainstorm. It can be generally understood as how many years such a large amount of rainfall occurs, which is the standard for road drainage design.

[0052] Digital Elevation Model (DEM): It is a digital expression of the terrain surface morphology by using limited terrain elevation data to realize digital simulation of the terrain. It is a kind of raster data, and each raster unit has a corresponding ground elevation value.

[0053] D8 algorithm: Based on the digital elevation model (DEM), it determines the direction of water flow by calculating the maximum distance weighted drop between the central grid and the neighboring grid. It assumes that when rain falls on a grid cell in the terrain, the water flow of the grid will flow to the grid cell with the smallest elevation among the surrounding 8 adjacent grid cells.

[0054] Depression contribution area: In the DEM data, the area where water flows to a specific depression is the contribution area of ​​the depression.

[0055] Overlay analysis: It is a commonly used spatial analysis method in geographic information systems. It involves geometric overlay operations on two or more spatial data layers to obtain comprehensive spatial information of these data.

[0056] Region growing calculation: an image segmentation technology whose core idea is to start from a set of "seed points" and add adjacent pixels with similar attributes to the corresponding region according to some similarity criteria (such as grayscale difference, color difference, etc. between pixels) until there are no pixels that meet the conditions to be added.

[0057] Cumulative runoff: In the digital elevation model (DEM) represented by a regular grid, there is a unit of water in each grid cell. According to the natural law that water flows from high to low, the amount of water flowing through each point is calculated based on the water flow direction data of the regional terrain, which is the cumulative runoff of the area.

[0058] The present invention obtains the spatial distribution grid of rainfall intensity through the rainfall data of the meteorological station in combination with the inverse distance weighted interpolation algorithm, and determines the urban drainage capacity grid according to the outdoor drainage standard. The difference between the two is used to obtain the direct surface water accumulation. The GIS analysis methods such as hydrological analysis and regional analysis are used to extract depressions, calculate the depression contribution area and depression depth, and for overflow water whose water accumulation is greater than the depression depth, the D8 algorithm is used to calculate the water flow direction, and based on this, the confluence accumulation is calculated to achieve surface water rebalancing and obtain the final water accumulation depth result; This method takes into account the surface runoff generation and convergence process, and uses the drainage capacity calculated by the rainfall intensity grid and the design rainstorm intensity of the pipe design return period to participate in the calculation of the direct surface water accumulation, thus realizing the calculation and simulation of water accumulation without pipe network data. The surface water rebalancing calculation method is realized by using the D8 algorithm and calculating the accumulated runoff: This method takes into account and calculates the dynamic rebalancing process of the water overflow process, thus improving the accuracy of the water accumulation simulation results.

[0059] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modification, equivalent substitution or improvement made within the spirit and principle of the present invention should be included in the protection scope of the present invention.

Claims

1. A rasterized calculation method for surface water in areas where pipe network data is missing, characterized in that: The following steps are involved: Step S1: Obtain rainfall data of the weather station that changes over time and calculate rainfall intensity data; Step S2: Calculate the amount of water accumulation, calculate the rainfall spatial distribution grid and the drainage capacity grid, and calculate the direct surface water accumulation based on the rainfall intensity data and the drainage capacity; Step S3: Calculate the depression depth spatial grid by using hydrological analysis and regional analysis to obtain the depression depth; Step S4: Determine the relationship between the amount of direct surface water accumulation and the depth of the depression based on the results of steps S2 and S3; Step S5: Rebalance the surface water volume, use the D8 algorithm to calculate the direction of the overflow water flow, calculate the cumulative amount of confluence, and perform rebalancing.

2. The method for calculating surface water gridding in areas where pipe network data is missing according to claim 1 is characterized in that: The rainfall intensity data is calculated in step S1 in the following manner: I=P t / t, where I is the rainfall intensity, P t is the accumulated rainfall in time t, where t is the rainfall time.

3. The method for calculating surface water gridding in areas where pipe network data is missing according to claim 1 is characterized in that: The rainfall spatial distribution grid in step S2 is calculated in the following way: The inverse distance weighted interpolation method is used to calculate the spatial distribution grid of rainfall intensity. The calculation formula is as follows: ， ， In the formula, is the rainfall intensity value of the interpolation point S; N is the number of meteorological stations used to interpolate at S, For the i The rainfall intensity at the ith weather station at λ; i is the weight of the ith weather station; d i is the distance from the interpolation point to the i-th weather station; b is the weight index.

4. The method for calculating surface water gridding in areas where pipe network data is missing according to claim 1 is characterized in that: The drainage capacity grid is calculated by the following steps: Step A1: Determine the drainage capacity of the region and use rainfall data to analyze the relationship between rainfall intensity I, rainfall duration t and return period P; Step A2: Calculation of direct surface water accumulation, the calculation formula is: △I=I in -I out , where △I is the direct surface water accumulation per hour, I in is the spatial grid value of rainfall intensity, I out It is the drainage capacity spatial grid value.

5. The method for calculating surface water gridding in areas where pipe network data is missing according to claim 4 is characterized in that: The relationship formula in step A1 is: , where I is rainfall intensity, t is rainfall duration, P is return period, A1 is the rainfall force of the design rainfall with a return period of 1 year, C is the rainfall force variation parameter, b is the rainfall duration correction parameter, and n is the rainfall attenuation index.

6. The method for calculating surface water gridding in areas where pipe network data is missing according to claim 1 is characterized in that: The depth spatial grid of the depression in step S3 is calculated in the following way: Step B1: Depression extraction, obtaining boundary data of the depression area to be extracted and digital elevation model data to analyze the water flow direction, and obtaining water flow direction raster data; Step B2: Identify and extract depression data, traverse the water flow grid data in a set order, and identify the depression grid data; Step B3: Overlay analysis is performed on the water flow direction grid data and the depression grid data to obtain the depression contribution area; Step B4: Overlay analysis of each depression contribution area with the digital elevation model data to calculate the depression depth.

7. The method for calculating surface water gridding in areas where pipe network data is missing according to claim 6 is characterized in that: The analysis of water flow direction is done in the following way: The digital elevation model is clipped according to the boundary data to obtain the DEM data of the specific area, which is then rasterized. A starting grid unit in the depression detection area after removing the boundary is selected, and a matrix is ​​constructed with the grid as the center. The D8 algorithm is used to analyze the water flow direction.

8. The method for calculating surface water gridding in areas where pipe network data is missing according to claim 7 is characterized in that: The specific analysis of water flow direction using the D8 algorithm is as follows: The direction of water flow is determined by the steepest slope method, and the distance difference between the central grid and the adjacent grid is calculated. The calculation formula is: ， In the formula, Drop is the maximum distance weighted drop between the center grid and its adjacent grid cells, H is the elevation value of the DEM data center grid cell, and H i is the elevation value of the adjacent grid cell, and i is the flow direction code relative to the central grid cell.

9. The method for calculating surface waterlogging in areas where pipe network data is missing according to claim 7 is characterized in that: The depression depth is calculated in the following way: Traverse the digital elevation model data grid cells in each depression contribution area and obtain its lowest elevation grid value Sink min , and then traverse the digital elevation model data grid cells of each depression contribution area boundary, extract the minimum elevation value of the regional boundary grid cell, that is, get the lowest elevation value of the depression contribution area boundary Sink max , Depth of Sink Depth =Sink max —Sink min .

10. The method for calculating surface water gridding in areas where pipe network data is missing according to claim 1, characterized in that: The step S4 is judged in the following manner: If the direct surface water accumulation is less than or equal to the depression depth, the required surface water accumulation depth is equal to the direct surface water accumulation. Otherwise, the surface water overflows and a rebalance calculation of the surface water accumulation is required.

Images