Urban multi-granularity catchment area division method based on high-precision three-dimensional data

Abstract

The invention relates to an urban multi-granularity catchment area division method based on high-precision three-dimensional data. The method comprises the following steps: acquiring urban high-precision three-dimensional data by using a laser point cloud technology; classifying the point cloud data by adopting a progressive triangulation network filtering method, and extracting vector boundary lines of building, road and flyover elements; hydrological analysis is carried out based on a high-precision DEM, a catchment area is extracted and combined with a boundary line, and a macroscopic boundary layer is constructed; irregular grids are dynamically divided by adopting a quadtree grid subdivision and merging mechanism, and a detail grid layer is constructed; dividing the detail grid layer into a runoff generation grid layer, a confluence grid layer and a runoff generation and confluence grid layer; the macroscopic boundary layer and the detail grid layer are integrated, hydrological parameters are given, and urban multi-granularity three-dimensional catchment area division is formed; and efficient storage and rapid retrieval are realized through object-oriented spatial data organization and GeoJSON storage. According to the method, the dividing fineness of the catchment area and the accuracy of dynamic simulation in the complex urban environment are improved.

Description

TECHNICAL FIELD

[0001] The present application belongs to the technical field of geographic information system and urban hydrology, and particularly relates to a method for dividing urban multi-granularity catchment area based on high-precision three-dimensional data. BACKGROUND

[0002] With the acceleration of urbanization, the underlying surface environment of urban three-dimensional space is increasingly complex, and buildings, roads, overpasses and ordinary ground jointly constitute a three-dimensional catchment network, which significantly changes the original hydrological cycle process of the city and deeply affects the surface permeability, water flow path and water distribution.

[0003] However, the traditional two-dimensional catchment division method ignores the complexity of the urban three-dimensional space, cannot accurately reflect the vertical distribution of hydrological geographical elements such as buildings and overpasses in the city, lacks the expression of the structural complexity of geographical space, and cannot accurately depict the slight ups and downs of the city terrain and the internal complex structure of the hydrological geographical elements; the traditional method is based on a relatively fixed catchment division method and a simplified hydrological parameter setting, and it is difficult to reflect the significant differences in elevation and water permeability in different regions of the city; the traditional method lacks the vertical layered water transfer mechanism of rainfall between different spatial levels, and it is difficult to accurately simulate the effects of hydrological geographical elements on water flow, such as resistance and guidance, aggregation and dispersion.

[0004] Therefore, the present application provides a method for dividing urban multi-granularity catchment area based on high-precision three-dimensional data, which integrates the fine depiction of the complex underlying surface structure of the city, the four-tree grid subdivision and merging mechanism, the rainfall intensity four-dimensional matrix driven runoff calculation and the water exchange mechanism between irregular grids, takes the hydrological geographical elements as the hydrological calculation unit, and divides the multi-granularity three-dimensional catchment area considering the macro boundary and micro details, effectively improving the fineness of the catchment division and the accuracy of the dynamic simulation in the complex urban environment, and having strong adaptability and expansibility. SUMMARY

[0005] In view of the above technical problems, the present application provides a method for dividing urban multi-granularity catchment area based on high-precision three-dimensional data, which aims to make full use of high-precision three-dimensional data to divide three-dimensional catchment area that can express the characteristics of urban three-dimensional space, multi-scale hydrological response mechanism and distributed heterogeneity elements, and provide a new method for improving the fineness of catchment division and the accuracy of dynamic simulation in complex urban environment.

[0006] In order to achieve the above application purpose, the technical scheme adopted by the present application is as follows:

[0007] A method for dividing urban multi-granularity catchment area based on high-precision three-dimensional data, the method comprising the following steps,

[0008] S1: Obtain high-precision three-dimensional data of a city by using a laser point cloud technology;

[0009] S2: Classify the point cloud data by using a progressive triangulation network filtering method, and extract the boundary lines of building, road and overpass elements having a hydrological regulation function in a vector form;

[0010] S3: Perform hydrological analysis based on high-precision DEM to extract a catchment area, and perform spatial merging with the vector boundary lines of S2 to construct a macro boundary layer;

[0011] S4: Divide irregular grids dynamically by using a quadtree grid subdivision and merging mechanism to construct a detail grid layer;

[0012] S5: Further divide the detail grid layer constructed in S4 into runoff, confluence and runoff-confluence grid layers based on the three-dimensional vertical layering characteristics of the city, and perform runoff-confluence calculation in combination with rainfall input driven by a four-dimensional matrix of rainfall intensity;

[0013] S6: Integrate the macro boundary layer of S3 and the detail grid layer of S5 in combination with the geometric morphological characteristics and hydrological processes of building, road, overpass and ordinary ground elements, assign corresponding hydrological parameters, and form a city multi-granularity three-dimensional catchment area;

[0014] S7: Realize efficient storage and rapid retrieval of the city multi-granularity three-dimensional catchment area by using an object-oriented spatial data organization and GeoJSON storage method;

[0015] Further, S2 includes the following steps:

[0016] S21: Divide the point cloud data into a plurality of grid blocks, select the lowest point in each grid block as an initial ground seed point to establish an initial ground triangulation network, and the specific formula for dividing the grid block is:

[0017]

[0018] In the above formula, row is the number of rows of the grid block, col is the number of columns of the grid block, l is the length of the point cloud, D is the diagonal size of the largest ground object, and p is the point cloud density;

[0019] S22: Determine whether a non-seed point P in the grid block is a ground point, if the two determination parameters of the included angle θ of the point P with the triangle and the distance d of the point P to the triangle plane are less than the set iteration angle and distance threshold, then the point P is determined to be a ground point, and the ground triangulation network is reconstructed, and the specific calculation formula is:

[0020]

[0021] In the above formula, θ is the included angle of the point P with the triangle, V iand V j are two vertices of the triangular face, d is the distance from point P to the triangular face, V k are vertices of the triangular face, P is the point cloud, and T is the triangular face.

[0022] S23: Filtering the ground points and non-ground points by iteratively processing S22 until all point cloud data in the grid block is processed, and finally classifying the building point cloud, road point cloud and overpass point cloud in the point cloud data;

[0023] S24: Converting the point cloud data classified in S23 into a raster surface, and extracting the boundary lines of the building, road and overpass hydrographic features by vector conversion.

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

[0025] S31: Based on the laser point cloud data, obtaining the ordinary ground point cloud after removing the buildings, roads and overpasses, and generating high-precision DEM by IDW inverse distance weighted interpolation;

[0026] S32: Adjusting the elevation value of the grid in the depression area to the same as the lowest elevation grid around by identifying the depression area in the DEM where the surrounding pixel elevation is higher than the center pixel, and filling the depression;

[0027] S33: Using D8 algorithm to calculate the flow direction of the DEM after filling the depression according to the elevation difference between adjacent grids;

[0028] S34: Starting from each grid cell of the DEM, calculating the number of cells flowing into the grid along the flow direction, marking the grid with a runoff accumulation amount reaching a set threshold as part of the river, and identifying the river network from the DEM to generate the river network;

[0029] S35: Based on the extracted river network, determining the outlet of each catchment area, and according to the flow direction data of each grid cell, tracking the upstream grid cells and assigning each grid cell to a unique catchment area, thereby defining the exact range of each catchment area;

[0030] S36: Spatially merging the catchment area extracted in S35 with the vector boundary line of the hydrographic feature extracted in S2 to construct a macro boundary layer.

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

[0032] S41: The missing part of the secondary layer point cloud of the double-layer roof building is filled by interpolation, and the terrain data under the overpass is generated by interpolating the ordinary ground point cloud. After obtaining the complete building, road, overpass and ordinary ground point cloud data, the inverse distance weighted interpolation method is used to convert to regular grid to generate DEM;

[0033] S42: Based on DEM, the terrain slope and curvature information are extracted, and the grid subdivision index is calculated, and the specific calculation formula is:

[0034]

[0035] In the above formula, C t is the grid subdivision index, S s is the slope of the current grid, C c is the curvature, S max and C max are the maximum values of the slope and curvature in the region respectively;

[0036] S43: According to the DEM, an initial quadtree structure is constructed, each node corresponds to a rectangular grid, and the initial resolution of the grid is set to I0, and the elevation value is obtained by bilinear interpolation on the DEM, thereby forming a regular initial grid structure, and the initial quadtree depth H0 is calculated according to the following formula:

[0037]

[0038] In the above formula, max(grid x ,grid y ) is the maximum pixel number of DEM data in the horizontal and vertical directions, H0 is the initial quadtree depth, and I0 is the initial resolution;

[0039] S44: Calculate the average grid subdivision index of each quadtree node coverage area, if it is greater than the grid subdivision threshold corresponding to the initial resolution, then the current node is subdivided into 4 subgrids, and the depth is increased by 1; otherwise, the grid remains independent, and the quadtree grid is subdivided until the subdivision termination condition is met or the maximum depth is reached. The elevation of the subdivided subgrid is also assigned by bilinear interpolation;

[0040] S45: Based on the subdivision in S44, calculate the grid merging index of each quadtree node, if it is less than the merging threshold, then the current node is merged with its parent node or adjacent low merging index node, until all mergable grid units are completed. The specific calculation formula is:

[0041]

[0042] E t = |S slope -Smean |

[0043] In the formula, S slope is the slope gradient in the current quadtree grid cell, S mean is the average slope gradient of the current grid and its adjacent grids, w is the elevation, a and b are the horizontal and vertical coordinates on the two-dimensional plane, respectively, and E t is the grid merging index.

[0044] S46: Extract all the leaf nodes from the quadtree structure after subdivision and merging, which together constitute a non-overlapping irregular grid set, thereby constructing a detailed grid layer.

[0045] Further, the S5 comprises the following steps:

[0046] S51: The irregular grid cells corresponding to the area participating in runoff calculation of the direct rainfall are divided into runoff grid layers, the runoff calculation adopts the Horton model, the infiltration rate of each runoff grid cell is calculated according to the rainfall duration and the set Horton model parameters, and the runoff is calculated in combination with the four-dimensional rainfall intensity matrix;

[0047] Further, the S51 comprises the following steps:

[0048] S511: For the rainfall stratification, runoff asynchronous superposition and vertical dynamic coupling phenomena caused by the three-dimensional structure of multi-layer roofs of buildings and spaces above and below interchanges, a four-dimensional rainfall intensity matrix including four collaborative dimensions of plane dimension, vertical dimension, time dimension and feature dimension is introduced;

[0049] S512: According to the time span of rainfall data, the rainfall data is uniformly segmented according to the set time step, each segment corresponds to the rainfall field at a specific time, and if there is no rainfall observation data at a certain time, the corresponding rainfall intensity is set to 0;

[0050] S513: Convert the rainfall data at each time into the rainfall grid of the study area at the corresponding time, assign a single rainfall value to small-scale areas, and use spline interpolation for large-scale areas;

[0051] S514: For the multi-level rainfall structure of different hydrological geographical elements, dynamically adjust the rainfall grid value of each rainfall level at each time, except for the part of the secondary roof and the space under the bridge that is blocked by the top roof, which adjusts the corresponding grid rainfall intensity to 0, and other elements do not need to be adjusted, finally the rainfall grids of different rainfall levels are spliced into complete three-dimensional rainfall distribution;

[0052] S515: Repeat S512 and S513 to obtain rainfall grids at n time points, and superimpose all rainfall grids at different time points to obtain a rainfall intensity four-dimensional matrix of the study area, and finally store the rainfall intensity four-dimensional matrix in NetCDF format.

[0053] S52: The irregular grid cells corresponding to the areas that do not directly receive precipitation but receive water input from the upper runoff-producing grid layer or other surrounding land surfaces are divided into a confluence grid layer, and the interlayer confluence in the surface layer and the interlayer confluence are calculated based on the water transfer in the vertical three-dimensional structure;

[0054] Further, the S52 comprises the following steps:

[0055] S521: Intra-layer grid confluence calculation in surface confluence, based on the irregular structure in the detailed grid layer, a two-dimensional shallow water equation is used to calculate the exchange of water flow between adjacent grids according to the water level difference between adjacent grids, and the specific calculation formula is:

[0056]

[0057]

[0058] In the above formula, t is time, q is flow vector, h is water depth, u and v are flow velocities in the horizontal and vertical directions, respectively, z is river bed elevation, u and v are flow velocities in the horizontal and vertical directions, respectively, n is Manning's coefficient, f and g are water flux vectors in the horizontal and vertical directions, respectively, S is the source term, including R is the source term vector of runoff, S b is the source term vector of bottom slope, S f is the source term vector of friction, and n is the Manning roughness coefficient.

[0059] S522: Interlayer grid confluence calculation in surface confluence, the common edge of the grid edge overlap in the vertical direction between the runoff-producing grid layer and the confluence grid layer is used to exchange water between layers, and the specific calculation formula is:

[0060] h n = h x,y + Δt(F x-1 / 2,y -F x+1 / 2,y +F x,y-1 / 2 -F x,y+1 / 2 ) / A x,y

[0061] In the above formula, F x-1 / 2,y is the water flux at the left boundary of the grid horizontal axis, F x,y+1 / 2Water flux at the boundary in the grid longitudinal axis direction, At is the time step, Δi and Δj are the size of the grid cell in the horizontal and vertical directions, respectively, (uh) and (vh) are the flow in the horizontal and vertical directions, respectively, (x, y) is the row and column index of the grid cell, A x,y is the (x, y) grid cell area, h n is the water depth after the current grid water exchange;

[0062] S53: The irregular grid cell corresponding to the area with runoff and confluence functions is divided into a runoff and confluence grid layer. Some areas directly receive precipitation to participate in runoff calculation and also receive water input from the upper runoff grid layer or the surrounding grids in the same layer. Through the coupling calculation of runoff and confluence, the total water quantity and flow characteristics of the runoff and confluence grid layer are obtained.

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

[0064] S61: The building element is accurately defined by the macro boundary layer to outline the building profile. The single-layer roof modifies the elevation of the grid by the runoff grid layer. The top roof of the double-layer roof uses the runoff grid layer to finely depict the building roof divide line to simulate the flow, distribution, water accumulation and drainage path of rainwater on the top roof. The secondary roof of the double-layer roof that is shielded by the top roof receives water from the top roof or the surrounding grids in the same layer by the confluence grid layer. The unshielded part of the secondary roof considers the combined action of direct rainfall input and top roof confluence by the runoff and confluence grid layer;

[0065] S62: The road element is described by the macro boundary layer to describe the road edge profile and is divided from the surrounding ordinary ground. The elevation fluctuation inside the road is finely depicted by the runoff grid layer to simulate the real water flow evolution path;

[0066] S63: The overpass element is described by the macro boundary layer to describe the bridge edge profile. The bridge surface uses the runoff grid layer to accurately depict the slope change details. The bridge underpass receives the water vertically transmitted from the upper bridge surface runoff grid layer and the lateral flow from the surrounding ground by the confluence grid layer;

[0067] S64: The ordinary ground element is described by the macro boundary layer to describe the ground edge profile. The terrain change is captured by the runoff grid layer;

[0068] S65: Based on the geometric shape characteristics of each element, the macro boundary layer boundary control coefficient, impermeable percentage and depression storage parameters are assigned;

[0069] S66: Based on the micro hydrological process of each element, the infiltration coefficient, surface friction coefficient and evaporation coefficient parameters of the detailed grid layer are assigned;

[0070] S67: Integrating the macro boundary layer and the detailed grid layer data of all hydrological geographical elements to form a multi-granularity three-dimensional catchment division.

[0071] Further, the S7 comprises the following steps:

[0072] S71: Building an object-oriented spatial data organization method, regarding each multi-granularity three-dimensional catchment as an object with a unique identifier, and realizing the index mapping and synchronous expression of the macro boundary layer and the detailed grid layer in space and attributes;

[0073] S72: Storing the multi-granularity catchment data in a physical storage in the form of GeoJSON.

[0074] An electronic device comprises a memory, a processor, and a computer program stored on the memory and executable on the processor, and the processor realizes the method for dividing a city multi-granularity catchment based on high-precision three-dimensional data when executing the program.

[0075] A computer readable storage medium has computer instructions stored thereon, and the computer instructions realize the method for dividing a city multi-granularity catchment based on high-precision three-dimensional data when executed by a processor.

[0076] Compared with the prior art, the present application has the following advantages and beneficial effects:

[0077] 1. The multi-granularity three-dimensional catchment division method fusing the macro boundary layer and the detailed grid layer is first proposed, the overall control ability and the local fine expression ability of the complex city space are considered, and the hydrological simulation is accurately expressed by assigning different parameters to each hydrological unit.

[0078] 2. The limitation of the traditional two-dimensional expression method is broken through, the laser point cloud data and the progressive triangular net filtering method are used, the three-dimensional expression of the hydrological geographical elements such as buildings, roads, overpasses and ordinary ground in the whole space is realized, and the distribution characteristics of the city space in the vertical direction are accurately described.

[0079] 3. For the vertical layering characteristics of the complex underlying surface of the city, three types of functional units of the runoff grid layer, the confluence grid layer and the runoff and confluence grid layer are constructed, the details of the hydrological process in different regions are reflected by dividing the irregular grid, the influence of the three-dimensional hydrological geographical elements such as buildings, roads, overpasses and ordinary ground on the hydrological process is combined, the multi-granularity spatial modeling of the hydrological process is realized, and the fine water flow path simulation is supported.

[0080] 4. A water exchange path in the vertical direction is constructed by identifying the common edges between the different level grid units of the runoff grid layer and the confluence grid layer.

[0081] 5. To overcome the rainfall input error caused by traditional rainfall homogenization, a four-dimensional rainfall intensity matrix is ​​proposed. This matrix comprehensively considers the plane, time, vertical, and characteristic dimensions, effectively simulating the process of rainfall reception, diversion, and vertical transmission in the three-dimensional structures of multi-story roofs and overpasses.

[0082] 6. Propose an object-oriented spatial data modeling method and physically store it in GeoJSON format to achieve unified management and index mapping of the macro boundary layer and the detailed grid layer, so as to support the standardized and modular management of multi-granularity three-dimensional watershed areas in cities and improve the retrieval efficiency and scalability of the system. BRIEF DESCRIPTION OF THE DRAWINGS

[0083] The drawings described herein are used to provide a further understanding of the embodiments of the present invention, constitute a part of this application, and do not constitute a limitation of the embodiments of the present invention. In the drawings:

[0084] Figure 1 It is a framework diagram of the method of the present invention;

[0085] Figure 2 This is a schematic diagram of the results of extracting road feature boundary lines in an example area of ​​the present invention;

[0086] Figure 3 This is a schematic diagram of the result of constructing the detailed grid layer of a single-layer roof building in the example area of ​​the present invention;

[0087] Figure 4 This is a schematic diagram of the result of constructing the detailed mesh layer of a double-roof building in the example area of ​​the present invention;

[0088] Figure 5 Schematic diagram of the construction result of the road detail grid layer in the example area of ​​the present invention;

[0089] Figure 6 This is a schematic diagram of the construction results of the detailed grid layer of the overpass in the example area of ​​the present invention;

[0090] Figure 7 This is a schematic diagram of the result of constructing the common ground detail grid layer in the example area of ​​the present invention;

[0091] Figure 8 This is a diagram of the physical storage architecture of data in GeoJSON format for the example area of ​​the present invention. DETAILED DESCRIPTION

[0092] In order to make the objectives, technical solutions and advantages of the present invention more clearly understood, the present invention is further described in detail below in conjunction with examples and drawings. The exemplary embodiments of the present invention and their descriptions are only used to explain the present invention and are not intended to limit the present invention.

[0093] Now take a typical waterlogging-prone area as an example, the target area has complete various hydrological geographical elements, high-precision three-dimensional data and rainfall observation data of rainfall stations and the like. The total area of the area is about 4470.31m 2 , the terrain is relatively small, and it is a relatively closed and independent catchment area.

[0094] As Figure 1 shown, the embodiment of the present application discloses a city multi-granularity catchment area division method based on high-precision three-dimensional data, mainly including the following steps:

[0095] S1: obtaining city high-precision three-dimensional data by using laser point cloud technology;

[0096] S2: classifying and processing the point cloud data by using the progressive triangulation network filtering method, and extracting the boundary lines of the building, road and overpass elements with hydrological regulation in vector form;

[0097] S3: performing hydrological analysis based on high-precision DEM to extract the catchment area, and performing spatial merging with the vector boundary line of S2 to construct a macro boundary layer;

[0098] S4: adopting a four-tree grid subdivision and merging mechanism to dynamically divide the irregular grid to construct a detail grid layer;

[0099] S5: based on the city three-dimensional vertical layering feature, further dividing the detail grid layer constructed by S4 into runoff, confluence and runoff-confluence grid layers, and performing runoff-confluence calculation in combination with rainfall input driven by a four-dimensional matrix of rainfall intensity;

[0100] S6: integrating the macro boundary layer of S3 and the detail grid layer of S5 in combination with the geometric shape features and hydrological processes of the building, road, overpass and ordinary ground elements, assigning corresponding hydrological parameters, and forming a city multi-granularity three-dimensional catchment area;

[0101] S7: realizing efficient storage and rapid retrieval of the city multi-granularity three-dimensional catchment area by using an object-oriented spatial data organization and GeoJSON storage method;

[0102] Further, the S2 includes the following steps:

[0103] S21: dividing the point cloud data into a plurality of grid blocks, selecting the lowest point in each grid block as an initial ground seed point to establish an initial ground triangulation network, realizing classification of the ground points, and the specific formula for dividing the grid block is:

[0104]

[0105] In the formula, row is the number of rows of the grid block, col is the number of columns of the grid block, l is the length of the point cloud, D is the diagonal size of the largest ground object, and p is the point cloud density;

[0106] S22: Determine whether the non-seed point P in the grid block is a ground point. If the two determination parameters of the included angle θ of the point P with the triangle and the distance d of the point P to the triangular face are less than the set iteration angle and distance threshold, the point P is determined to be a ground point, and the ground triangular network is reconstructed. The specific calculation formula is:

[0107]

[0108] In the formula, θ is the included angle of the point P with the triangular face, V i and V j are two vertices of the triangular face, d is the distance of the point P to the triangular face, V k is the vertex of the triangular face, P is the point cloud, and T is the triangular face.

[0109] S23: Continuously repeat the determination condition of S22, filter the ground points and non-ground points through iterative processing, and classify the building point cloud, road point cloud and overpass point cloud in the point cloud data until all the point cloud data in the grid block are processed.

[0110] S24: Convert the point cloud data classified in S23 into a raster surface, and extract the boundary lines of the building, road and overpass hydrographic geographical elements through vector conversion, as shown in Figure 2 .

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

[0112] S31: Based on the laser point cloud data, obtain the ordinary ground point cloud after removing the buildings, roads and overpasses, and generate a high-precision DEM by using IDW inverse distance weighted interpolation;

[0113] S32: Adjust the elevation value of the grid in the depression area to be the same as the lowest elevation grid in the surrounding area by identifying the depression area in the DEM where the surrounding pixel elevations are higher than the center pixel, and filling the depression area;

[0114] S33: Use the D8 algorithm to calculate the flow direction of the DEM after filling the depression area according to the elevation difference between adjacent grids;

[0115] S34: Start from each grid cell of the DEM, calculate the number of cells flowing into the grid along the flow direction, mark the grid with a flow accumulation amount reaching a set threshold as part of the river, and identify the river network from the DEM to generate a river network;

[0116] S35: Based on the extracted river network, the outlet of each catchment is determined, and according to the water flow direction data of each grid cell, the upstream grid cells are tracked, each grid cell is assigned to a unique catchment, and the exact range of each catchment is determined;

[0117] S36: The catchment area extracted in S35 is spatially merged with the vector boundary line of the hydrological geographical element extracted in S2 to construct a macro boundary layer.

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

[0119] S41: The point cloud of the secondary layer of the double-layer roof building is interpolated to fill the part blocked by the top roof; the point cloud of the ordinary ground is interpolated to generate the terrain data under the overpass. After obtaining the complete building, road, overpass and ordinary ground point cloud data, inverse distance weighted interpolation method is used to convert to regular grid to generate DEM;

[0120] S42: Based on the DEM, the terrain slope and curvature information are extracted, and the grid subdivision index is calculated, and the specific calculation formula is:

[0121]

[0122] In the above formula, C t is the grid subdivision index, S s is the slope of the current grid, C c is the curvature, S max and C max are the maximum values of the slope and curvature in the region respectively;

[0123] S43: According to the DEM, an initial quadtree structure is constructed, each node corresponds to a rectangular grid, the initial resolution of the grid is set to I0, and the elevation value is obtained by bilinear interpolation on the DEM, thereby forming a regular initial grid structure, and the initial quadtree depth H0 is calculated according to the following formula:

[0124]

[0125] In the above formula, max(grid x ,grid y ) is the maximum pixel number of the DEM data in the horizontal and vertical directions, H0 is the initial quadtree depth, and I0 is the initial resolution;

[0126] S44: According to the grid subdivision index C tThe initial resolution S0 is used to determine whether the initial grid meets the terrain accuracy requirements. If not, the grid is subdivided. The average grid subdivision index of the area covered by each quadtree node is calculated. If it is greater than the grid subdivision threshold corresponding to the initial resolution, the current node is subdivided into 4 subgrids and the depth is increased by 1. Otherwise, the grid remains independent and the quadtree grid is continuously subdivided until the subdivision termination condition is met or the maximum depth is reached. The elevation of the subdivided subgrids is also assigned using bilinear interpolation.

[0127] S45: Based on the S44 subdivision, the grid merging index is used to determine whether the grid meets the low complexity area. If it meets the requirements, the grid is merged to optimize the computational efficiency of the simple terrain area while ensuring accuracy. The grid merging index of each quadtree node is calculated. If it is less than the merging threshold, the current node is merged with its parent node or adjacent low merging index nodes until all mergable grid cells are completed. The specific calculation formula is:

[0128]

[0129] E t =|S slope -S mean |

[0130] In the above formula, S slope is the slope gradient within the current quadtree grid cell, S mean is the average slope gradient of the current grid and its adjacent grids, w is the elevation, a and b are the horizontal and vertical coordinates on the two-dimensional plane, E t is the grid merging index.

[0131] S46: Extract all leaf nodes from the subdivided and merged quadtree structure. These leaf nodes together constitute a non-overlapping non-regularized grid set, thereby constructing a detail grid layer.

[0132] Furthermore, the S5 includes the following steps:

[0133] S51: The irregular grid cells corresponding to the area that directly receives precipitation and participates in runoff calculation are divided into runoff grid layers. The runoff calculation adopts the Horton model. According to the rainfall duration and the set Horton model parameters, the infiltration rate of each runoff grid cell is calculated, and the runoff is calculated by combining the rainfall intensity four-dimensional matrix.

[0134] Furthermore, the S51 includes the following steps:

[0135] S511: For the phenomena of rainfall stratification, runoff asynchronous superposition and vertical dynamic coupling caused by the three-dimensional structure of building multi-layer roof, overpass up and down space, a four-dimensional matrix of rainfall intensity including four coordinated dimensions of plane dimension, vertical dimension, time dimension and characteristic dimension is introduced as the driving input, wherein the plane dimension represents the position of the rainfall grid in the horizontal direction, the vertical dimension represents the different rain layers in the three-dimensional space, the time dimension represents the time evolution of the whole process of the rainfall event, and the characteristic dimension represents the physical characteristics of the rainfall intensity of the rainfall itself;

[0136] S512: According to the time span of the rainfall data from July 1, 2016 to July 2, 2016, the rainfall data is segmented according to the set time step of 1h, and the whole time interval is evenly divided into 24 blocks, each block representing the rainfall field at a certain time. If there is no rainfall observation data at a certain time, the corresponding rainfall intensity is set to 0;

[0137] S513: The time sequence rainfall data of the Dinghuimen rainfall station site is directly assigned as the rainfall of the whole study area, and converted into the rainfall grid of the study area at the corresponding time. For small-scale urban areas, due to the small area range, it can be considered that the rainfall received is the average, and the rainfall grid values of the whole study area are directly assigned to the rainfall data of a single rainfall station. For large-scale urban areas, the area range is large, and the rainfall is non-uniformly distributed. A spline interpolation spatial interpolation algorithm is used to interpolate the rainfall intensity to obtain the spatial distribution of the rainfall intensity in the whole area;

[0138] S514: For the multi-level rain structure of different hydrological geographical elements, the rainfall grid values of different rain layers at each time are dynamically adjusted. For ordinary ground and roads, the rainfall is directly received; for buildings, the roof of single-story building and the top roof of double-story building directly receive the rainfall, and part of the rainfall of the secondary roof of double-story building is blocked by the top roof and does not directly receive the rainfall, so the rainfall grid value of the blocked part needs to be adjusted to 0; for overpass, the bridge surface directly receives the rainfall, and the rainfall of the space under the bridge is blocked by the bridge surface and does not directly receive the rainfall, so the rainfall grid value under the bridge needs to be adjusted to 0. After the adjustment, the rainfall grids of different rain layers of various hydrological geographical elements at 12 o'clock on July 1, 2016 are spliced together as the 12th block of the four-dimensional matrix of rainfall intensity;

[0139] S515: S512 and S513 are repeatedly executed to obtain the rainfall grids at 24 times, and all the rainfall grids at different times are superimposed to obtain the four-dimensional matrix of rainfall intensity of the study area. Finally, the four-dimensional matrix of rainfall intensity is stored in the NetCDF format.

[0140] S52: the irregular grid cells corresponding to the areas that do not directly receive precipitation but receive water input from the upper runoff-producing grid layer or other surrounding land surface are divided into the confluence grid layer, and the interlayer confluence in the surface layer is calculated based on water transfer on the vertical stereoscopic structure;

[0141] Further, the S52 includes the following steps:

[0142] S521: the intralayer grid confluence calculation in the surface confluence is based on the irregular grid structure in the detailed grid layer, a two-dimensional shallow water equation is used to calculate the exchange of water flow between the grids according to the water level difference between adjacent grids, and the specific calculation formula is:

[0143]

[0144] In the above formula, t is time, q is flow vector, h is water depth, u and v are flow velocities in the horizontal and vertical directions respectively, z is the elevation of the riverbed bottom surface, u and v are flow velocities in the horizontal and vertical directions respectively, n is the Manning coefficient, f and g are water flux vectors in the horizontal and vertical directions respectively, S is the source term, R is the source term vector of runoff, S b is the source term vector of the bottom slope, S f is the source term vector of friction, and n is the Manning roughness coefficient.

[0145] S522: the interlayer grid confluence calculation in the surface confluence is performed through the common edge of the grid edge overlap in the vertical direction between the runoff-producing grid layer and the confluence grid layer for interlayer water exchange, and the specific calculation formula is:

[0146] h n = h x,y + Δt (F x-1 / 2,y -F x+1 / 2,y +F x,y-1 / 2 -F x,y+1 / 2 ) / A x,y

[0147] In the above formula, F x-1 / 2,y is the water flux at the left boundary of the grid horizontal axis, F x,y+1 / 2 is the water flux at the upper boundary of the grid vertical axis, Δt is the time step, Δi and Δj are the sizes of the grid unit in the horizontal and vertical directions respectively, (uh) and (uh) are the flows in the horizontal and vertical directions respectively, (x, y) is the row and column index of the grid unit, A x,y is the area of the (x, y) grid unit, h n is the water depth after water exchange of the current grid;

[0148] S53: The irregular grid unit corresponding to the area with the function of runoff generation and convergence is divided into a runoff generation and convergence grid layer. Some areas directly receive precipitation and participate in runoff calculation, and also receive water input from the upper runoff generation grid layer or the surrounding grid in the same layer. Through the calculation of coupling runoff generation and convergence, the total water quantity and flow characteristics of the runoff generation and convergence grid layer are obtained.

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

[0150] S61: The building element is accurately defined by the macro boundary layer to outline the building profile, and the single-layer roof is modified by the runoff generation grid layer to modify the elevation of the grid, as shown in Figure 3 The top roof of the double-layer roof uses the runoff generation grid layer to finely depict the building roof divide line to simulate the flow, distribution, water accumulation and drainage path of rainwater on the top roof. The secondary roof blocked by the top roof receives water from the top roof or the surrounding grid in the same layer through the convergence grid layer. The unblocked part of the secondary roof considers the combined action of direct rainfall input and top roof convergence through the runoff and convergence grid layer, as shown in Figure 4 (a) Top roof runoff generation grid layer (b) Secondary roof runoff and convergence grid layer.

[0151] S62: The road element is described by the macro boundary layer to describe the road edge profile and is divided from the surrounding ordinary ground. The elevation fluctuation inside the road is finely depicted by the runoff generation grid layer to simulate the real water flow evolution path, as shown in Figure 5

[0152] S63: The overpass element is described by the macro boundary layer to describe the bridge edge profile. The bridge surface uses the runoff generation grid layer to accurately depict the slope change details. The bridge underpass receives the water vertically transmitted from the upper bridge surface runoff generation grid layer and the lateral flow from the surrounding ground through the convergence grid layer, as shown in Figure 6 (a) Overpass bridge surface runoff generation grid layer; (b) Overpass bridge underpass convergence grid layer.

[0153] S64: The ordinary ground element is described by the macro boundary layer to describe the ground edge profile. The terrain change is captured by the runoff generation grid layer, as shown in Figure 7

[0154] S65: Based on the geometric shape characteristics of each element, the macro boundary layer is given boundary control coefficient, impermeable percentage and depression storage parameters. The boundary control coefficient is used to define the blocking and guiding effect of the hydrological geographical element boundary on water flow. The impermeable percentage reflects the water permeability of the area. The depression storage parameter represents the rainwater interception capacity of the area. The hydrological parameter value range of the macro boundary layer and the detail grid layer is shown in the following table.

[0155]

[0156] ​​S66: Based on the micro-hydrological process of each element, the infiltration coefficient, surface friction coefficient and evaporation coefficient of the detailed grid layer are given. The infiltration coefficient describes the water permeability of the ground surface, the surface friction coefficient affects the flow rate of surface runoff, and the evaporation coefficient reflects the evaporation loss of ground water. The hydrological parameters of the macro boundary layer and the detailed grid layer of each hydrological geographical element are as shown in the following table:

[0157]

[0158] S67: Integrate the macro boundary layer and the detailed grid layer data of all hydrological geographical elements to form a multi-granularity three-dimensional catchment division. The entire region is divided into 435 multi-granularity three-dimensional catchment areas, each of which contains a macro boundary layer and a detailed grid layer. The detailed grid layer is divided into 84796 irregular grid cells of different resolutions.

[0159] Further, the S7 comprises the following steps:

[0160] S71: Construct an object-oriented spatial data organization method to realize the index mapping of the macro boundary layer and the detailed grid layer in space and attributes. Each hydrological geographical element corresponds to a multi-granularity three-dimensional catchment area, and each multi-granularity three-dimensional catchment area has a unique identifier ID. The macro boundary layer and the detailed grid layer included therein also have unique identifiers ID. The boundary point pointer records the vector boundary coordinate sequence. The surface area pointer records the grid information in the detailed grid layer, including the grid elevation and other hydrological properties of the grid. The hydrological element relationship marks which hydrological geographical element the three-dimensional catchment area is divided based on. The flow direction relationship marks the outlet flow direction of the catchment area, which may point to a pipe point or another three-dimensional catchment area.

[0161] S72: Based on the GeoJSON format, the data is physically stored and the corresponding mapping relationship is stored. The storage architecture presents an inverted tree structure, realizing standardized interoperability. The entry file is index.json, and the type attribute of the root object of the file is FeatureCollection. The features attribute stores all the GeoJSON objects corresponding to the multi-granularity three-dimensional catchment area. Each GeoJSON object has a type of Feature, and the geometry attribute stores the corresponding GeoJSON object of the macro boundary layer. Each GeoJSON object has a type of Polygon, and the properties attribute stores the grid data ASCII file pointer corresponding to the detailed grid layer, as shown in Figure 8 .

[0162] 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. In conclusion, the content of the present description should not be understood as a limitation of the present application.

Claims

1. A method for dividing urban watersheds into multiple granularities based on high-precision three-dimensional data, characterized in that: The following steps are involved: S1: Use laser point cloud technology to obtain high-precision three-dimensional data of the city; S2: The point cloud data is classified using the progressive triangulation filtering method to extract the boundary lines of buildings, roads, and overpasses with hydrological regulation functions in vector form; S3: Perform hydrological analysis based on high-precision DEM to extract the catchment area, spatially merge it with the vector boundary lines of S2, and construct a macro boundary layer; S4: Uses quadtree mesh subdivision and merging mechanism to dynamically divide irregular meshes to construct detail mesh layers; S5: Based on the three-dimensional vertical stratification characteristics of the city, the detailed grid layer constructed in S4 is further divided into runoff generation, runoff confluence, and runoff generation and confluence grid layers. The runoff generation and confluence calculation is performed in combination with the rainfall input driven by the four-dimensional rainfall intensity matrix. S6: Combining the geometric features and hydrological processes of buildings, roads, overpasses, and common ground elements, the macro boundary layer of S3 is integrated with the detailed grid layer of S5, and corresponding hydrological parameters are assigned to form a multi-granularity three-dimensional urban watershed. S7: Through object-oriented spatial data organization and GeoJSON storage methods, efficient storage and fast retrieval of urban multi-granularity three-dimensional watershed areas are achieved.

2. The method for dividing urban watersheds into multiple granularities based on high-precision three-dimensional data according to claim 1, characterized in that: Boundary line extraction in S2 includes the following steps: S21: Divide the point cloud data into several grid blocks, select the lowest elevation point in each grid block as the initial ground seed point to establish the initial ground triangulation. The specific formula for dividing the grid blocks is: In the above formula, row is the number of rows in the grid block, col is the number of columns in the grid block, l is the length of the point cloud, D is the diagonal size of the largest feature, and ρ is the point cloud density; S22: Determine whether the non-seed point P in the grid block is a ground point. If the angle θ between the point P and the triangle and the distance d between the point P and the triangle are less than the set iteration angle and distance thresholds, then the point P is determined to be a ground point and the ground triangulation network is reconstructed. The specific calculation formula is: In the above formula, θ is the angle between point P and the triangle, V i and V j are the two vertices of the triangle, d is the distance from point P to the triangle, V k is the vertex of the triangle, P is the point cloud, and T is the triangle; S23: filtering ground points and non-ground points by iteratively processing S22 until all point cloud data in the grid block are processed, and finally classifying the point cloud data into building point clouds, road point clouds, and overpass point clouds; S24: The point cloud data obtained by classification in S23 is converted into a raster surface, and the boundary lines of the hydrological and geographical elements of buildings, roads and overpasses are extracted through vector conversion.

3. The method for dividing urban watersheds into multiple granularities based on high-precision three-dimensional data according to claim 1, characterized in that: The macro boundary layer construction in S3 includes the following steps: S31: Based on the laser point cloud data, obtain the ordinary ground point cloud after removing buildings, roads and overpasses, and use IDW inverse distance weighted interpolation to generate a high-precision DEM; S32: By identifying the depression area in the DEM where the surrounding pixel elevation is higher than the central pixel, the elevation value of the grid in the depression is adjusted to the same as the surrounding lowest elevation grid to fill the depression; S33: Using the D8 algorithm, the water flow direction is calculated for the DEM after filling based on the elevation difference between adjacent grids; S34: Starting from each grid cell of the DEM, the number of cells flowing into the grid from upstream is accumulated along the flow direction, and the grid whose cumulative flow reaches the set threshold is marked as part of the river. The river network is identified from the DEM, thereby generating a river network; S35: Based on the extracted river network, the outlet of each catchment area is determined, and according to the water flow direction data of each grid cell, the upstream grid cells are tracked, and each grid cell is assigned to a unique catchment area, thereby defining the exact scope of each catchment area; S36: Spatially merge the catchment area extracted by S35 with the vector boundary lines of the hydrological and geographical elements extracted by S2 to construct a macro boundary layer.

4. The method for dividing urban watersheds into multiple granularities based on high-precision three-dimensional data according to claim 1, characterized in that: The non-regular grid division in S4 includes the following steps: S41: Interpolate the secondary point cloud of the double-roofed building to fill the part blocked by the top roof; interpolate the ordinary ground point cloud to generate the terrain data under the overpass. After obtaining the complete building, road, overpass and ordinary ground point cloud data, use the inverse distance weighted interpolation method to convert it into a regular grid to generate a DEM; S42: Extract terrain slope and curvature information based on DEM and calculate grid subdivision index. The specific calculation formula is: In the above formula, C t is the grid subdivision index, S s is the slope of the current grid, C c is the curvature, S max and C max are the maximum values ​​of slope and curvature in the region, respectively; S43: Construct an initial quadtree structure based on the DEM. Each node corresponds to a rectangular grid. The initial resolution of the grid is set to I0. Its elevation value is obtained by bilinear interpolation of the DEM to form a regularized initial grid structure. The specific calculation formula for the initial quadtree depth H0 is: In the above formula, max(grid x ,grid y ) is the maximum number of pixels in the horizontal and vertical directions of the DEM data, H0 is the initial quadtree depth, and I0 is the initial resolution; S44: Calculate the average grid subdivision index of the area covered by each quadtree node. If it is greater than the grid subdivision threshold corresponding to the initial resolution, subdivide the current node into 4 subgrids and increase the depth by 1. Otherwise, keep the grid independent and continue to subdivide the quadtree grid until the subdivision termination condition is met or the maximum depth is reached. The elevation of the subdivided subgrids is also assigned by bilinear interpolation. S45: Based on the S44 subdivision, the grid merging index of each quadtree node is calculated. If it is less than the merging threshold, the current node is merged with its parent node or adjacent nodes with low merging indexes until all merging grid cells are completed. The specific calculation formula is: E t =|S slope -S mean | In the above formula, S slope is the slope gradient within the current quadtree grid cell, S mean is the average slope gradient of the current grid and its adjacent grids, w is the elevation, a and b are the horizontal and vertical coordinates on the two-dimensional plane, respectively, E t is the grid merging index; S46: Extract all leaf nodes from the subdivided and merged quadtree structure. These leaf nodes together constitute a non-overlapping non-regularized grid set, thereby constructing a detail grid layer.

5. The method for dividing urban watersheds into multiple granularities based on high-precision three-dimensional data according to claim 1, characterized in that: The hierarchical structure division in S5 includes the following steps: S51: The irregular grid cells corresponding to the area that directly receives precipitation and participates in runoff calculation are divided into runoff grid layers. The runoff calculation adopts the Horton model. According to the rainfall duration and the set Horton model parameters, the infiltration rate of each runoff grid cell is calculated, and the runoff is calculated by combining the rainfall intensity four-dimensional matrix. S52: The irregular grid cells corresponding to the area that does not directly receive precipitation but receives water input from the upper runoff grid layer or other surrounding surfaces are divided into confluence grid layers. Based on the water transfer in the vertical three-dimensional structure, the confluence within and between surface layers is calculated; S53: The irregular grid cells corresponding to areas with both runoff generation and runoff confluence are divided into runoff generation and confluence grid layers. Some areas directly receive precipitation to participate in runoff generation calculations, and also receive water input from the upper runoff generation grid layer or the surrounding grids in the same layer. By coupling the calculation of runoff generation and confluence, the total water volume and water flow characteristics of the runoff generation and confluence grid layer are obtained.

6. The method for dividing urban watersheds into multiple granularities based on high-precision three-dimensional data according to claim 1, characterized in that: The multi-granularity three-dimensional catchment area division in S6 includes the following steps: S61: Building elements are precisely defined using a macro-boundary layer. A runoff raster layer is used to modify the grid elevation of single-story roofs. A runoff raster layer is used to finely characterize the roof watershed of a double-story rooftop to simulate the flow, distribution, accumulation, and drainage path of rainwater on the top rooftop. Sub-roofs shielded by the top roof receive water from the top rooftop or surrounding grids on the same layer through a confluence raster layer. The unshielded portion of the sub-roof is treated through a combined consideration of direct rainfall input and top roof runoff through a confluence raster layer. S62: The road element uses a macro boundary layer to describe the road edge contour and the boundary between it and the surrounding ordinary ground. The runoff raster layer is used to finely depict the elevation fluctuations inside the road to simulate the actual water flow evolution path. S63: The overpass element uses a macro boundary layer to describe the edge contour of the bridge deck. The bridge deck uses a runoff grid layer to accurately depict the details of slope changes. The bridge receives water vertically transferred from the runoff grid layer above the bridge deck through a confluence grid layer, as well as lateral overflow from the surrounding surface. S64: Common ground features describe the ground edge contours through a macro boundary layer and capture terrain changes through a runoff raster layer; S65: Based on the geometric characteristics of each element, the macro boundary layer boundary control coefficient, impervious percentage and depression water storage parameters are assigned; S66: Based on the micro-hydrological process of each element, the infiltration coefficient, surface friction coefficient and evaporation coefficient parameters are assigned to the detail grid layer; S67: Integrate the macro boundary layer and detailed grid layer data of all hydrological and geographical elements to form a multi-granularity three-dimensional watershed division.

7. The method for dividing urban watersheds into multiple granularities based on high-precision three-dimensional data according to claim 5, characterized in that: The four-dimensional matrix of rainfall intensity in S51 includes the following steps: S511: To address the phenomena of rainfall stratification, asynchronous superposition of runoff generation, and vertical dynamic coupling caused by the three-dimensional structures of multi-story roofs of buildings and upper and lower spaces of overpasses, a four-dimensional rainfall intensity matrix is ​​introduced, which includes four coordinated dimensions: plane dimension, vertical dimension, time dimension, and characteristic dimension. S512: Based on the time span of the rainfall data, the rainfall data is evenly segmented according to the set time step. Each segment corresponds to the rainfall field at a specific time. If there is no rainfall observation data at a certain time, the corresponding rainfall intensity is set to 0; S513: Convert the rainfall data at each moment into the rainfall grid of the study area at the corresponding moment, assign a single rainfall value to the small-scale area, and use spline interpolation for the large-scale area; S514: Dynamically adjust the rainfall grid values ​​for each rainfall-receiving layer at each moment, targeting the multi-layered rainfall structure of different hydrological and geographical elements. Except for the secondary rooftops and the space under the bridge that are blocked by the top roof, the corresponding grid rainfall intensity is adjusted to 0. No adjustment is required for other elements. Finally, the rainfall grids for different rainfall-receiving layers are spliced ​​into a complete three-dimensional rainfall distribution. S515: Repeat S512 and S513 to obtain rainfall grids at n moments, and superimpose the rainfall grids at all moments into a four-dimensional matrix of rainfall intensity in the study area. Finally, the four-dimensional matrix of rainfall intensity is stored in NetCDF format.

8. The method for dividing urban watersheds into multiple granularities based on high-precision three-dimensional data according to claim 5, characterized in that: The confluence calculation in S52 includes the following steps: S521: Surface flow calculation within a grid layer is based on the irregular grid structure in the detail grid layer. It uses the two-dimensional shallow water equation to calculate the water flow exchange between grids according to the water level difference between adjacent grids. The specific calculation formula is: In the above formula: t is time, q is flow vector, h is water depth, u and v are flow velocities in the horizontal and vertical directions respectively, z is riverbed elevation, u and v are flow velocities in the horizontal and vertical directions respectively, n is Manning coefficient, f and g are water flux vectors in the horizontal and vertical directions respectively, S is source term, including R is source term vector of runoff, S b is the source vector of the bottom slope, S f is the source term vector of friction, n is the Manning roughness coefficient; S522: Calculation of inter-layer grid confluence in surface runoff. Water exchange between layers is performed through the common edges of the runoff-producing grid layer and the confluence grid layer in the vertical direction. The specific calculation formula is: h n =h x,y +Δt(F x-1 / 2,y -F x+1 / 2,y +F x,y-1 / 2 -F x,y+1 / 2 ) / A x,y In the above formula, F x-1 / 2,y is the water flux at the left boundary of the grid in the horizontal direction, F x,y+1 / 2 The water flux at the boundary in the vertical direction of the grid, Δt is the time step, Δi and Δj are the dimensions of the grid unit in the horizontal and vertical directions, (uh) and (vh) are the flow rates in the horizontal and vertical directions, (x, y) are the row and column indices of the grid unit, A x,y is the (x,y) grid cell area, h n It is the water depth after the water volume of the current grid is exchanged.

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, it implements a method for dividing urban multi-granularity watershed areas based on high-precision three-dimensional data as described in any one of claims 1 to 8.

10. A computer-readable storage medium having computer instructions stored thereon, characterized in that: When the computer instructions are executed by a processor, a method for dividing urban multi-granularity watershed areas based on high-precision three-dimensional data as described in any one of claims 1 to 8 is implemented.

Images