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

Abstract

This invention relates to a method for urban multi-granularity catchment area delineation based on high-precision 3D data, comprising the following steps: acquiring high-precision 3D urban data using laser point cloud technology; classifying the point cloud data using a progressive triangulation filtering method to extract vector boundary lines of building, road, and overpass elements; performing hydrological analysis based on a high-precision DEM to extract catchment areas and merge them with boundary lines to construct a macroscopic boundary layer; dynamically dividing irregular grids using a quadtree grid subdivision and merging mechanism to construct a detailed grid layer; dividing the detailed grid layer into runoff generation, runoff confluence, and runoff-confluence raster layers; integrating the macroscopic boundary layer and the detailed grid layer, assigning hydrological parameters to form a multi-granularity 3D urban catchment area delineation; and achieving efficient storage and rapid retrieval through object-oriented spatial data organization and GeoJSON storage. This invention improves the precision of catchment area delineation and the accuracy of dynamic simulation in complex urban environments.

Description

Technical Field

[0001] This invention belongs to the fields of geographic information system technology and urban hydrology, specifically relating to a method for dividing urban multi-granularity catchment areas based on high-precision three-dimensional data. Background Technology

[0002] With the acceleration of urbanization, the underlying surface environment of urban three-dimensional space is becoming increasingly complex. Buildings, roads, overpasses and ordinary ground together form a three-dimensional water catchment network, which significantly changes the original hydrological cycle process of the city and profoundly affects surface infiltration, water flow path and water accumulation distribution.

[0003] However, traditional two-dimensional catchment area delineation methods ignore the complexity of urban three-dimensional space, cannot accurately reflect the vertical distribution of hydrogeographic elements such as buildings and overpasses in the city, lack the expression of the structural complexity of geographical space, and cannot accurately depict the subtle undulations of urban topography and the internal complex structure of hydrogeographic elements. Traditional methods are based on relatively fixed catchment area delineation methods and simplified hydrological parameter settings, making it difficult to reflect the significant differences in elevation and permeability of different areas of the city. Traditional methods lack the ability to capture the vertical stratification of rainfall transfer mechanisms between different spatial levels, making it difficult to accurately simulate the effects of hydrogeographic elements on water flow obstruction and guidance, aggregation and dispersion.

[0004] Therefore, this invention proposes a method for dividing urban multi-granularity catchment areas based on high-precision three-dimensional data. This method integrates the detailed characterization of complex urban underlying surface structures, quadtree grid subdivision and merging mechanisms, runoff generation and runoff calculation driven by a four-dimensional matrix of rainfall intensity, and water volume interaction mechanisms between irregular grids. Using hydrogeographic elements as hydrological calculation units, it divides multi-granularity three-dimensional catchment areas that take into account both macroscopic boundaries and microscopic details. This effectively improves the precision of catchment area division and the accuracy of dynamic simulation in complex urban environments, and has strong adaptability and scalability. Summary of the Invention:

[0005] To address the problems existing in the above-mentioned technologies, this invention provides a method for urban multi-granularity catchment area delineation based on high-precision three-dimensional data. The purpose is to make full use of high-precision three-dimensional data to delineate three-dimensional catchment areas that can express the three-dimensional spatial characteristics of the city, multi-scale hydrological response mechanisms, and heterogeneous distribution elements, providing a new method to improve the precision of catchment area delineation and the accuracy of dynamic simulation in complex urban environments.

[0006] To achieve the above-mentioned objectives, the technical solution adopted by this invention is as follows:

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

[0008] S1: Acquire high-precision 3D urban data using laser point cloud technology;

[0009] S2: The point cloud data is classified using the progressive triangulation filtering method, and the boundary lines of buildings, roads and overpasses with hydrological regulation functions are extracted in vector form;

[0010] S3: Hydrological analysis based on high-precision DEM to extract catchment areas, spatially merge with the vector boundary line of S2, and construct a macroscopic boundary layer;

[0011] S4: Employs a quadtree mesh subdivision and merging mechanism to dynamically divide irregular meshes in order to construct detailed mesh layers;

[0012] S5: Based on the three-dimensional vertical layering features of the city, the detailed mesh layer constructed in S4 is further divided into runoff generation, runoff confluence, and runoff generation and confluence grid layers. Runoff generation and confluence calculation is performed by combining the rainfall input driven by the four-dimensional matrix of rainfall intensity.

[0013] S6: Combining the geometric features of buildings, roads, overpasses and ordinary ground elements with hydrological processes, the macro boundary layer of S3 and the detailed grid layer of S5 are integrated and assigned corresponding hydrological parameters to form a multi-granular three-dimensional water catchment area in the city.

[0014] S7: Achieve efficient storage and fast retrieval of urban multi-granular three-dimensional catchment areas through object-oriented spatial data organization and GeoJSON storage methods;

[0015] Furthermore, S2 includes the following steps:

[0016] S21: Divide the point cloud data into several grid blocks, and 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 as follows:

[0017]

[0018] 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.

[0019] S22: Determine whether a non-seed point P in a grid block is a ground point. If the angle θ between point P and the triangle and the distance d from point P to the triangle face are both less than the set thresholds for iteration angle and distance, then point P is determined to be a ground point, and the ground triangulation is reconstructed. The specific calculation formula is as follows:

[0020]

[0021] In the above formula, θ is the angle between point P and the triangular face, and V iand V j Let P be the two vertices of the triangle, d be the distance from point P to the triangle, and V be the distance from point P to the triangle. k P is the vertex of the triangle, P is the point cloud, and T is the triangle face;

[0022] S23: By iteratively processing S22, ground points and non-ground points are filtered until all point cloud data within the grid block has been processed, and finally, building point clouds, road point clouds and overpass point clouds are classified from the point cloud data.

[0023] S24: Convert the point cloud data obtained from the classification in S23 into a raster surface, and extract the boundary lines of hydrogeographic elements such as buildings, roads and overpasses through vector transformation.

[0024] Furthermore, S3 includes the following steps:

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

[0026] S32: By identifying depression areas in the DEM where the elevation of surrounding pixels is higher than that of the center pixel, the elevation value of the raster within the depression is adjusted to be the same as that of the lowest elevation raster around it, and the depression is filled.

[0027] S33: Using the D8 algorithm, the water flow direction is calculated for the DEM after filling depressions based on the elevation difference between adjacent grids;

[0028] S34: Starting from each grid cell of the DEM, the number of cells flowing into the grid from upstream is accumulated along the direction of water flow. Grids whose cumulative flow reaches a set threshold are marked as part of the river. The river network is identified from the DEM, thereby generating the river network.

[0029] S35: Based on the extracted river network, identify the outlet of each catchment area, and track the upstream grid cells according to the water flow direction data of each grid cell, assign each grid cell to a unique catchment area, thereby defining the exact range of each catchment area.

[0030] S36: Spatially merge the vector boundary lines of the catchment area extracted in S35 and the hydrogeographic elements extracted in S2 to construct a macro boundary layer.

[0031] Furthermore, S4 includes the following steps:

[0032] S41: By interpolating the secondary point cloud of double-roofed buildings, the portion obscured by the top roof is filled; by interpolating ordinary ground point cloud, terrain data under the overpass is generated. After acquiring complete point cloud data of buildings, roads, overpasses, and ordinary ground, inverse distance weighted interpolation is used to convert them into regular grids to generate a DEM.

[0033] S42: Extract terrain slope and curvature information based on DEM, and calculate grid subdivision index. The specific calculation formula is as follows:

[0034]

[0035] In the above formula, C t It is a grid subdivision index, S s C is the slope of the current grid. c It is curvature, S max and C max These are the maximum values ​​of slope and curvature within the region, respectively.

[0036] S43: Construct an initial quadtree structure based on the DEM. Each node corresponds to a rectangular grid. The initial grid resolution is set to I0, and its elevation values ​​are obtained by bilinear interpolation of the DEM, thus forming a regularized initial grid structure. The specific formula for calculating the initial quadtree depth H0 is as follows:

[0037]

[0038] In the above formula, max(grid x ,grid y H0 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.

[0039] 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, the current node is subdivided into 4 subgrids and the depth is increased by 1. Otherwise, the independent state of the grid is maintained and the quadtree grid is continuously subdivided until the subdivision termination condition is met or the maximum depth is reached. The elevation of the subgrid after subdivision is also assigned by bilinear interpolation.

[0040] S45: Based on the subdivision in S44, calculate the mesh merging index for each quadtree node. If it is less than the merging threshold, merge the current node with its parent node or adjacent nodes with lower merging indices until all mergeable mesh cells have been completed. The specific calculation formula is as follows:

[0041]

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

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

[0044] S46: Extract all leaf nodes from the subdivided and merged quadtree structure. These leaf nodes together form a set of non-overlapping irregular meshes, thus constructing the detailed mesh layer.

[0045] Furthermore, S5 includes the following steps:

[0046] S51: The irregular grid cells corresponding to the areas that directly receive precipitation and participate in runoff generation calculation are divided into runoff generation grid layers. The runoff generation calculation adopts the Horton model. Based on the rainfall duration and the set Horton model parameters, the infiltration rate of each runoff generation grid cell is calculated, and the runoff generation is calculated in combination with the four-dimensional matrix of rainfall intensity.

[0047] Furthermore, S51 includes the following steps:

[0048] S511: In response to the phenomena of stratified rainfall reception, asynchronous superposition of runoff, and vertical dynamic coupling caused by the three-dimensional structure of multi-story building roofs and overpasses, a four-dimensional rainfall intensity matrix with four collaborative dimensions including planar dimension, vertical dimension, time dimension, and feature dimension is introduced.

[0049] S512: Based on the time span of the rainfall data, the rainfall data is evenly divided into segments 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.

[0050] S513: Convert the rainfall data at each time point into a rainfall raster of the study area at the corresponding time point. Assign a single rainfall value to small-scale areas and use spline interpolation for large-scale areas.

[0051] S514: For the multi-layered rain-receiving structure of different hydrogeographic elements, dynamically adjust the rainfall grid values ​​of different rain-receiving layers at each moment. Except for the secondary roof part and the space under the bridge that are blocked by the top roof, the corresponding grid rainfall intensity is adjusted to 0. Other elements do not need to be adjusted. Finally, the rainfall grids of different rain-receiving layers are spliced ​​into a complete three-dimensional rainfall distribution.

[0052] S515: Repeatedly execute S512 and S513 to obtain rainfall rasters at n time points, and superimpose the rainfall rasters at all time points into a four-dimensional matrix of rainfall intensity in the study area. Finally, store the four-dimensional matrix of rainfall intensity in NetCDF format.

[0053] S52: Irregular grid cells corresponding to regions that do not directly receive precipitation but receive water input from the upper runoff grid layer or other surrounding surfaces are divided into runoff grid layers. Based on water transfer on the vertical three-dimensional structure, the runoff within and between the surface layers is calculated.

[0054] Furthermore, S52 includes the following steps:

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

[0056]

[0057]

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

[0059] S522: Interlayer grid runoff calculation in surface runoff, using the common edge of the vertical overlap between the runoff-generating grid layer and the runoff-collecting grid layer to perform interlayer water exchange. The specific calculation formula is as follows:

[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 It is the water flux at the left boundary of the grid along the horizontal axis, F x,y+1 / 2The water flux at the boundary along the vertical axis of the grid is given by Δt, where Δt is the time step, Δi and Δj are the dimensions of the grid cell in the horizontal and vertical directions, respectively, (uh) and (vh) are the flow rates in the horizontal and vertical directions, respectively, and (x, y) are the row and column indices of the grid cell. x,y h is the area of ​​the (x,y) grid cell. n It is the water depth after the current grid water exchange;

[0062] S53: Irregular grid cells corresponding to regions that have both runoff generation and runoff collection functions are divided into runoff generation and runoff collection grid layers. Some regions 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 of the same layer. By coupling the calculations of runoff generation and runoff collection, the total water volume and flow characteristics of the runoff generation and runoff collection grid layer are obtained.

[0063] Furthermore, S6 includes the following steps:

[0064] S61: Building elements are precisely defined by a macro boundary layer. For single-layer roofs, the elevation of the grid is modified by a flow generation grid layer. For double-layer roofs, the top layer uses a flow generation grid layer to finely depict the roof water distribution line to simulate the flow, distribution, accumulation, and drainage path of rainwater on the top layer. The second-layer roof, which is blocked by the top layer, receives water from the top layer or the surrounding grids of the same layer through a flow collection grid layer. The unblocked part of the second-layer roof takes into account the combined effect of direct rainfall input and flow collection from the top layer through a flow generation and collection grid layer.

[0065] S62: Road elements are described by a macroscopic boundary layer to define the road edge contour and demarcate it from the surrounding ordinary ground. A runoff-generating grid layer is used to finely depict the elevation undulations within the road to simulate the actual water flow evolution path.

[0066] S63: The overpass elements describe the edge contour of the bridge deck through a macro boundary layer, and the bridge deck uses a runoff generation grid layer to accurately depict the details of slope changes; the water volume vertically transferred from the upper bridge deck runoff generation grid layer and the lateral runoff from the surrounding ground surface are received under the bridge through a runoff collection grid layer.

[0067] S64: Ordinary ground features describe the ground edge contour through a macro boundary layer and capture terrain changes through a runoff raster layer;

[0068] S65: Based on the geometric morphological characteristics of each element, assign the macro boundary layer boundary control coefficient, impermeability percentage, and depression water storage parameters;

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

[0070] S67: Integrate macroscopic boundary layer and detailed grid layer data of all hydrogeographic elements to form a multi-granular three-dimensional catchment area division.

[0071] Furthermore, S7 includes the following steps:

[0072] S71: Construct an object-oriented spatial data organization method, treating each multi-granularity 3D catchment area as an object with a unique identifier, and realize the index mapping and synchronous expression of macro boundary layer and detailed mesh layer in space and attributes;

[0073] S72: Uses GeoJSON format for physical data storage, storing multi-granularity catchment area data in an inverted tree structure.

[0074] An electronic device includes a memory, a processor, and a computer program stored in the memory and executable on the processor. When the processor executes the program, it implements the aforementioned method for urban multi-granularity catchment area delineation based on high-precision three-dimensional data.

[0075] A computer-readable storage medium having computer instructions stored thereon, which, when executed by a processor, implement the aforementioned method for urban multi-granularity catchment area delineation based on high-precision three-dimensional data.

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

[0077] 1. A multi-granularity three-dimensional catchment area delineation method integrating macro boundary layer and detailed grid layer is proposed for the first time. It takes into account both the overall control capability and the local fine expression capability of complex urban space. By assigning differentiated parameters to each hydrological unit, the accurate expression of hydrological simulation is achieved.

[0078] 2. Breaking through the limitations of traditional two-dimensional representation methods, this method utilizes laser point cloud data and progressive triangulation filtering to achieve full-space three-dimensional representation of hydrogeographic elements such as buildings, roads, overpasses, and ordinary ground, accurately depicting the vertical distribution characteristics of urban space.

[0079] 3. To address the vertical stratification characteristics of complex urban underlying surfaces, three functional units are constructed: runoff generation grid layer, runoff collection grid layer, and runoff generation and collection grid layer. By dividing irregular grids, the details of hydrological processes in different areas are reflected. The influence of the characteristics of three-dimensional hydrogeographic elements such as buildings, roads, overpasses, and ordinary ground on hydrological processes is combined to achieve multi-granular spatial modeling of hydrological processes and support refined flow path simulation.

[0080] 4. Establish a water exchange mechanism between irregular grids, and propose to construct a vertical water exchange path by identifying the common edges between different levels of grid cells in the runoff generation grid layer and the runoff collection grid layer;

[0081] 5. To overcome the rainfall input error caused by traditional rainfall homogenization processing, a four-dimensional rainfall intensity matrix is ​​proposed, which comprehensively considers four dimensions: plane, time, vertical and feature, and effectively simulates the process of rainfall reception, diversion and vertical transmission by the three-dimensional spatial structure of multi-story roofs and overpasses.

[0082] 6. An object-oriented spatial data modeling approach is proposed, and physical storage is performed in GeoJSON format to achieve unified management and index mapping of macro boundary layer and detailed mesh layer, so as to support the standardized and modular management of multi-granularity three-dimensional watersheds in cities and improve the retrieval efficiency and scalability of the system. Attached Figure Description

[0083] The accompanying drawings, which are included to provide a further understanding of embodiments of the invention and form part of this application, do not constitute a limitation thereof. In the drawings:

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

[0085] Figure 2 This is a schematic diagram of the road feature boundary line extraction results in the example area of ​​this invention;

[0086] Figure 3 This is a schematic diagram of the detailed mesh layer construction result of a single-story roof building in the example area of ​​this invention;

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

[0088] Figure 5 A schematic diagram of the road detail mesh layer construction result in an example area of ​​this invention;

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

[0090] Figure 7 This is a schematic diagram of the construction result of a common ground detail mesh layer in the example area of ​​this invention;

[0091] Figure 8 This is a diagram of the physical storage architecture for GeoJSON formatted data in the example area of ​​this invention. Detailed Implementation

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

[0093] Taking a typical flood-prone area as an example, the target area possesses complete hydrological and geographical features, high-precision three-dimensional data, and rainfall observation data from rain gauge stations. The total area of ​​this region is approximately 4470.31 m². 2 The terrain is relatively gentle, making it a relatively enclosed and independent catchment area.

[0094] like Figure 1 As shown in the figure, the present invention discloses a method for urban multi-granularity catchment area delineation based on high-precision three-dimensional data, which mainly includes the following steps:

[0095] S1: Acquire high-precision 3D urban data using laser point cloud technology;

[0096] S2: The point cloud data is classified using the progressive triangulation filtering method, and the boundary lines of buildings, roads and overpasses with hydrological regulation functions are extracted in vector form;

[0097] S3: Hydrological analysis based on high-precision DEM to extract catchment areas, spatially merge with the vector boundary line of S2, and construct a macroscopic boundary layer;

[0098] S4: Employs a quadtree mesh subdivision and merging mechanism to dynamically divide irregular meshes in order to construct detailed mesh layers;

[0099] S5: Based on the three-dimensional vertical layering features of the city, the detailed mesh layer constructed in S4 is further divided into runoff generation, runoff confluence, and runoff generation and confluence grid layers. Runoff generation and confluence calculation is performed by combining the rainfall input driven by the four-dimensional matrix of rainfall intensity.

[0100] S6: Combining the geometric features of buildings, roads, overpasses and ordinary ground elements with hydrological processes, the macro boundary layer of S3 and the detailed grid layer of S5 are integrated and assigned corresponding hydrological parameters to form a multi-granular three-dimensional water catchment area in the city.

[0101] S7: Achieve efficient storage and fast retrieval of urban multi-granular three-dimensional catchment areas through object-oriented spatial data organization and GeoJSON storage methods;

[0102] Further, S2 includes the following steps:

[0103] S21: Divide the point cloud data into several grid blocks, and select the lowest elevation point in each grid block as the initial ground seed point to establish an initial ground triangulation network, thereby classifying the ground points. The specific formula for dividing the grid blocks is as follows:

[0104]

[0105] 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.

[0106] S22: Determine whether a non-seed point P in a grid block is a ground point. If the angle θ between point P and the triangle and the distance d from point P to the triangle face are both less than the set thresholds for iteration angle and distance, then point P is determined to be a ground point, and the ground triangulation is reconstructed. The specific calculation formula is as follows:

[0107]

[0108] In the above formula, θ is the angle between point P and the triangular face, and V i and V j Let P be the two vertices of the triangle, d be the distance from point P to the triangle, and V be the distance from point P to the triangle. k P is the vertex of the triangle, P is the point cloud, and T is the triangle face;

[0109] S23: Repeat the S22 judgment condition continuously, and encrypt layer by layer through iterative processing to filter ground points and non-ground points until all point cloud data in the grid block has been processed. Finally, the point cloud data is classified into building point cloud, road point cloud and overpass point cloud.

[0110] S24: Convert the point cloud data obtained from the classification in S23 into a raster surface, and extract the boundary lines of hydrogeographic features such as buildings, roads, and overpasses through vector transformation, such as... Figure 2 As shown.

[0111] Further, S3 includes the following steps:

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

[0113] S32: By identifying depression areas in the DEM where the elevation of surrounding pixels is higher than that of the center pixel, the elevation value of the raster within the depression is adjusted to be the same as that of the lowest elevation raster around it, and the depression is filled.

[0114] S33: Using the D8 algorithm, the water flow direction is calculated for the DEM after filling depressions based on the elevation difference between adjacent grids;

[0115] S34: Starting from each grid cell of the DEM, the number of cells flowing into the grid upstream is accumulated along the direction of water flow. Grids whose accumulated flow reaches a set threshold are marked as part of the river. The river network is identified from the DEM, thereby generating the river network.

[0116] S35: Based on the extracted river network, identify the outlet of each catchment area, and track the upstream grid cells according to the water flow direction data of each grid cell, assign each grid cell to a unique catchment area, thereby defining the exact range of each catchment area.

[0117] S36: Spatially merge the vector boundary lines of the catchment area extracted in S35 and the hydrogeographic elements extracted in S2 to construct a macro boundary layer.

[0118] Further, S4 includes the following steps:

[0119] S41: By interpolating the secondary point cloud of double-roofed buildings, the portion obscured by the top roof is filled; by interpolating ordinary ground point cloud, terrain data under the overpass is generated. After acquiring complete point cloud data of buildings, roads, overpasses, and ordinary ground, inverse distance weighted interpolation is used to convert them into regular grids to generate a DEM.

[0120] S42: Extract terrain slope and curvature information based on DEM, and calculate grid subdivision index. The specific calculation formula is as follows:

[0121]

[0122] In the above formula, C t It is a grid subdivision index, S s C is the slope of the current grid. c It is curvature, S max and C max These are the maximum values ​​of slope and curvature within the region, respectively.

[0123] S43: Construct an initial quadtree structure based on the DEM. Each node corresponds to a rectangular grid. The initial grid resolution is set to I0, and its elevation values ​​are obtained by bilinear interpolation of the DEM, thus forming a regularized initial grid structure. The specific formula for calculating the initial quadtree depth H0 is as follows:

[0124]

[0125] In the above formula, max(grid x ,grid y H0 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.

[0126] S44: Based on grid subdivision index C tGiven an initial resolution S0, determine if the initial grid meets the terrain accuracy requirements. If not, subdivide the grid. Calculate the average grid subdivision index for 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, maintain the independent state of the grid and continuously 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 using bilinear interpolation.

[0127] S45: Building upon the subdivision in S44, a mesh merging index is used to determine whether a mesh meets the criteria for a low-complexity region. If it does, the meshes are merged, thus optimizing computational efficiency for simple geometric regions while maintaining accuracy. The mesh merging index for 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 indices until all mergeable mesh cells have been completed. The specific calculation formula is as follows:

[0128]

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

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

[0131] S46: Extract all leaf nodes from the subdivided and merged quadtree structure. These leaf nodes together form a set of non-overlapping irregular meshes, thus constructing the detailed mesh layer.

[0132] Further, S5 includes the following steps:

[0133] S51: The irregular grid cells corresponding to the areas that directly receive precipitation and participate in runoff generation calculation are divided into runoff generation grid layers. The runoff generation calculation adopts the Horton model. Based on the rainfall duration and the set Horton model parameters, the infiltration rate of each runoff generation grid cell is calculated, and the runoff generation is calculated in combination with the four-dimensional matrix of rainfall intensity.

[0134] Further, S51 includes the following steps:

[0135] S511: In response to phenomena such as stratified rainfall reception, asynchronous superposition of runoff, and vertical dynamic coupling caused by three-dimensional structures such as multi-layered roofs of buildings and the vertical space of overpasses, a four-dimensional rainfall intensity matrix with four collaborative dimensions, including planar dimension, vertical dimension, time dimension, and feature dimension, is introduced as the driving input. The planar dimension represents the position of the rainfall grid in the horizontal direction, the vertical dimension represents the different rain-receiving layers in three-dimensional space, the time dimension represents the temporal evolution of the entire rainfall event, and the feature dimension represents the physical characteristics of the rainfall intensity itself.

[0136] S512: Based on the rainfall data spanning from July 1st to July 2nd, 2016, the rainfall data is segmented according to a set time step of 1 hour. The entire time interval is evenly divided into 24 blocks, each representing the rainfall field at a specific moment. If there is no rainfall observation data at a certain moment, the corresponding rainfall intensity is set to 0.

[0137] S513: The time-series rainfall data from the Dinghuaimen rain gauge station are directly assigned to the rainfall data of the entire study area, converting it into a rainfall raster for the corresponding time in the study area. For small-scale urban areas, due to the small area, the rainfall can be considered to be average, and all rainfall raster values ​​in the study area are directly assigned to the rainfall data of a single rain gauge station. For large-scale urban areas, the area is large and the rainfall is non-uniformly distributed. Spline interpolation spatial interpolation algorithm is used to interpolate the rainfall intensity to obtain the spatial distribution of rainfall intensity throughout the entire area.

[0138] S514: For the multi-layered rainfall-receiving structure of different hydrogeographic elements, dynamically adjust the rainfall grid values ​​for different rainfall-receiving layers at each moment. For ordinary ground and roads, rainfall is received directly. For buildings, the roofs of single-story buildings and the top roofs of double-story buildings directly receive rainfall. For double-story buildings, some rainfall on the second-floor roof is blocked by the top roof and does not directly receive rainfall; the rainfall grid value for the blocked portion needs to be adjusted to 0. For overpasses, the bridge deck directly receives rainfall, while rainfall in the space beneath the bridge is blocked by the bridge deck and does not directly receive rainfall; the rainfall grid value for the bridge surface needs to be adjusted to 0. After adjustment, the rainfall grids of different rainfall-receiving layers for various hydrogeographic elements at 12:00 on July 1, 2016, are stitched together as the 12th block of the four-dimensional rainfall intensity matrix.

[0139] S515: Repeatedly execute S512 and S513 to obtain rainfall rasters at 24 time points, and superimpose the rainfall rasters at all time points into a four-dimensional matrix of rainfall intensity in the study area. Finally, store the four-dimensional matrix of rainfall intensity in NetCDF format.

[0140] S52: Irregular grid cells corresponding to regions that do not directly receive precipitation but receive water input from the upper runoff grid layer or other surrounding surfaces are divided into runoff grid layers. Based on water transfer on the vertical three-dimensional structure, the runoff within and between the surface layers is calculated.

[0141] Further, S52 includes the following steps:

[0142] S521: Intra-layer grid runoff calculation in surface runoff, based on the irregular grid structure in the detailed grid layer, using the two-dimensional shallow water equation, calculates the exchange of water flow between grids according to the water level difference between adjacent grids. The specific calculation formula is as follows:

[0143]

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

[0145] S522: Interlayer grid runoff calculation in surface runoff, using the common edge of the vertical overlap between the runoff-generating grid layer and the runoff-collecting grid layer to perform interlayer water exchange. The specific calculation formula is as follows:

[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 It is the water flux at the left boundary of the grid along the horizontal axis, F x,y+1 / 2 The water flux at the boundary along the vertical axis of the grid is given by Δt, where Δt is the time step, Δi and Δj are the dimensions of the grid cell in the horizontal and vertical directions, respectively, (uh) and (uh') are the flow rates in the horizontal and vertical directions, respectively, and (x, y) are the row and column indices of the grid cell. x,y h is the area of ​​the (x,y) grid cell. n It is the water depth after the current grid water exchange;

[0148] S53: Irregular grid cells corresponding to regions that have both runoff generation and runoff collection functions are divided into runoff generation and runoff collection grid layers. Some regions 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 of the same layer. By coupling the calculations of runoff generation and runoff collection, the total water volume and flow characteristics of the runoff generation and runoff collection grid layer are obtained.

[0149] Further, S6 includes the following steps:

[0150] S61: Building elements are precisely defined by a macroscopic boundary layer, and single-story roofs have their grid elevation modified by a flow generation grid layer, such as... Figure 3 As shown. The top layer of the double-roof uses a flow-generating grid layer to meticulously depict the building's roof water distribution line, simulating the flow, distribution, accumulation, and drainage paths of rainwater on the top layer roof. The second-layer roof, shielded by the top layer roof, receives water from the top layer roof or surrounding grids on the same layer through a flow-generating grid layer. The unshielded portion of the second-layer roof, through the flow-generating and flow-collecting grid layer, comprehensively considers the combined effects of direct rainfall input and flow collection from the top layer roof, such as... Figure 4 As shown; (a) Top roof flow generation grid layer (b) Secondary roof flow generation grid layer.

[0151] S62: Road elements are described by a macroscopic boundary layer, defining the road edge contour and separating it from the surrounding ordinary ground. A runoff-generating grid layer finely depicts the elevation variations within the road to simulate the actual water flow path, such as... Figure 5 As shown;

[0152] S63: The interchange elements describe the bridge deck edge contour through a macroscopic boundary layer, while the bridge deck uses a runoff-generating grid layer to precisely depict the details of slope changes; beneath the bridge, a runoff-collecting grid layer receives the water vertically transferred from the upper bridge deck runoff-generating grid layer, as well as the lateral runoff from the surrounding surface, such as... Figure 6 As shown; (a) the runoff generation grid layer on the overpass deck; (b) the runoff collection grid layer under the overpass.

[0153] S64: Ordinary ground features describe the ground edge contour through a macroscopic boundary layer and capture terrain changes through a runoff raster layer, such as... Figure 7 As shown;

[0154] S65: Based on the geometric morphological characteristics of each element, a boundary control coefficient, impermeability percentage, and depression water storage parameter are assigned to the macro boundary layer. The boundary control coefficient is used to define the boundary of hydrogeographic elements and their role in blocking and guiding water flow. The impermeability percentage reflects the permeability of the area. The depression water storage parameter characterizes the area's ability to retain rainwater. The range of values ​​for the hydrological parameters of the macro boundary layer and the detailed grid layer is shown in the table below.

[0155]

[0156] S66: Based on the micro-hydrological processes of each element, the detailed grid layer is assigned parameters such as infiltration coefficient, surface friction coefficient, and evaporation coefficient. The infiltration coefficient describes the permeability of the surface, the surface friction coefficient affects the flow rate of surface runoff, and the evaporation coefficient reflects the evaporation loss of surface water. The values ​​of the hydrological parameters of the macro-boundary layer and detailed grid layer of each hydrogeographic element are shown in the table below.

[0157]

[0158] S67: Integrating the macro-boundary layer and detailed grid layer data of all hydrogeographic elements to form a multi-granularity three-dimensional catchment area 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 84,796 irregular grid units of different resolutions.

[0159] Further, S7 includes the following steps:

[0160] S71: Construct an object-oriented spatial data organization method to realize the spatial and attribute index mapping between the macro boundary layer and the detailed mesh layer. Each hydrogeographic element corresponds to a multi-granularity three-dimensional catchment area. Each multi-granularity three-dimensional catchment area has a unique identifier ID, and its included macro boundary layer and detailed mesh layer also have unique identifier IDs. Boundary point pointers record the vector boundary coordinate sequence; area pointers record the information of each grid in the detailed mesh layer, including the grid elevation and other hydrological attributes; hydrogeographic element relationship marks which hydrogeographic element the three-dimensional catchment area is based on; flow direction relationship marks the flow direction of the catchment area outlet, which may point to a pipe point or another three-dimensional catchment area.

[0161] S72: Physically stores data and corresponding mapping relationships based on GeoJSON format. The storage architecture presents an inverted tree structure, achieving standardized interoperability. Its entry file is `index.json`, the root object's `type` property is `FeatureCollection`, and the `features` property stores the GeoJSON objects corresponding to all multi-granularity 3D catchment areas. Each GeoJSON object has a `type` of `Feature`, and its `geometry` property stores the GeoJSON objects corresponding to the macroscopic boundary layers it contains. Each GeoJSON object also has a `type` of `Polygon`, and its `properties` property stores the ASCII file pointers to the raster data corresponding to the detailed mesh layers it contains, such as... Figure 8 As shown.

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

Claims

1. A method for urban multi-granularity catchment area delineation based on high-precision three-dimensional data, characterized in that, Includes the following steps: S1: Acquire high-precision 3D urban data using laser point cloud technology; S2: The point cloud data is classified using the progressive triangulation filtering method, and the boundary lines of buildings, roads and overpasses with hydrological regulation functions are extracted in vector form; S3: Hydrological analysis based on high-precision DEM is performed to extract the catchment area, which is then spatially merged with the vector boundary line of S2 to construct a macroscopic boundary layer; S4: Employs a quadtree mesh subdivision and merging mechanism to dynamically divide irregular meshes in order to construct detailed mesh layers; S5: Based on the three-dimensional vertical layering features of the city, the detailed mesh layer constructed in S4 is further divided into runoff generation, runoff confluence, and runoff generation and confluence raster layers. Runoff generation and confluence calculation is performed by combining the rainfall input driven by the four-dimensional matrix of rainfall intensity. S6: Combining the geometric features of buildings, roads, overpasses and ordinary ground elements with hydrological processes, the macro boundary layer of S3 and the detailed grid layer of S5 are integrated and assigned corresponding hydrological parameters to form a multi-granular three-dimensional water catchment area in the city. S7: Achieve efficient storage and fast retrieval of urban multi-granular three-dimensional catchment areas through object-oriented spatial data organization and GeoJSON storage methods; The irregular mesh generation in S4 includes the following steps: S41: By interpolating the secondary point cloud of the double-roofed building, the part obscured by the top roof is filled; by interpolating the ordinary ground point cloud, the terrain data under the overpass is generated. After obtaining complete building, road, overpass and ordinary ground point cloud data, the inverse distance weighted interpolation method is used 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 as follows: In the above formula, C t It is a grid subdivision index, S s C is the slope of the current grid. c It is curvature, S max and C max These are the maximum values ​​of slope and curvature within the region, respectively. S43: Construct an initial quadtree structure based on the DEM. Each node corresponds to a rectangular grid. The initial grid resolution is set to I0, and its elevation values ​​are obtained by bilinear interpolation of the DEM, thus forming a regularized initial grid structure. The specific formula for calculating the initial quadtree depth H0 is as follows: In the above formula, This refers to the maximum number of pixels in the horizontal and vertical directions of the DEM data. It is the initial quadtree depth. This 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, the current node is subdivided into 4 subgrids and the depth is increased by 1. Otherwise, the independent state of the grid is maintained and the quadtree grid is continuously subdivided until the subdivision termination condition is met or the maximum depth is reached. The elevation of the subgrid after subdivision is also assigned by bilinear interpolation. S45: Based on the subdivision in S44, calculate the mesh merging index for each quadtree node. If it is less than the merging threshold, merge the current node with its parent node or adjacent nodes with lower merging indices until all mergeable mesh cells have been completed. The specific calculation formula is as follows: E t =|S slope -S mean | In the above formula, S is the slope gradient within the current quadtree grid cell. mean The slope gradient is the average gradient of the current grid and its neighboring grids, where w is the elevation. and These are the horizontal and vertical coordinates on a two-dimensional plane. It is a grid merging indicator; S46: Extract all leaf nodes from the subdivided and merged quadtree structure. These leaf nodes together form a set of non-overlapping irregular meshes, thereby constructing the detailed mesh layer.

2. The method for urban multi-granularity catchment area delineation 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, and 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 as follows: In the above formula, `col` is the row number of the grid block, and `col` is the column number of the grid block. It is the length of the point cloud. It is the diagonal size of the largest feature. It is point cloud density; S22: Determine whether a non-seed point P in a grid block is a ground point. If the angle between point P and the triangle is... If the two decision parameters, the distance d from point P to the triangular face, are less than the set thresholds for the iteration angle and distance, then point P is determined to be a ground point, and the ground triangulation is reconstructed. The specific calculation formula is as follows: In the above formula, It is the angle between point P and the triangle face. and Let P be the two vertices of the triangle, d be the distance from point P to the triangle, and V be the distance from point P to the triangle. k It is the vertex of the triangle. It's a point cloud. It is a triangular face; S23: By iteratively processing S22, ground points and non-ground points are filtered until all point cloud data within the grid block has been processed, and finally, building point clouds, road point clouds and overpass point clouds are classified from the point cloud data. S24: Convert the point cloud data obtained from the classification in S23 into a raster surface, and extract the boundary lines of hydrogeographic elements such as buildings, roads and overpasses through vector transformation.

3. The method for urban multi-granularity catchment area delineation based on high-precision three-dimensional data according to claim 1, characterized in that, The construction of the macroscopic boundary layer in S3 includes the following steps: S31: Based on laser point cloud data, obtain ordinary ground point cloud after removing buildings, roads and overpasses, and generate high-precision DEM using IDW inverse distance weighted interpolation; S32: By identifying depression areas in the DEM where the elevation of surrounding pixels is higher than that of the center pixel, the elevation value of the raster within the depression is adjusted to be the same as that of the lowest elevation raster around it, and the depression is filled. S33: The D8 algorithm is used to calculate the water flow direction of the DEM after filling the depressions 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 upstream is accumulated along the direction of water flow. Grids whose cumulative flow reaches a set threshold are marked as part of the river. The river network is identified from the DEM, thereby generating the river network. S35: Based on the extracted river network, identify the outlet of each catchment area, and track the upstream grid cells according to the water flow direction data of each grid cell, assign each grid cell to a unique catchment area, thereby defining the exact range of each catchment area. S36: Spatially merge the vector boundary lines of the catchment areas extracted in S35 and the hydrogeographic elements extracted in S2 to construct a macro boundary layer.

4. The method for urban multi-granularity catchment area delineation based on high-precision three-dimensional data according to claim 1, characterized in that, The hierarchical structure partitioning in S5 includes the following steps: S51: The irregular grid cells corresponding to the areas that directly receive precipitation and participate in runoff generation calculation are divided into runoff generation grid layers. The runoff generation calculation adopts the Horton model. Based on the rainfall duration and the set Horton model parameters, the infiltration rate of each runoff generation grid cell is calculated, and the runoff generation is calculated in combination with the four-dimensional matrix of rainfall intensity. S52: Irregular grid cells corresponding to regions that do not directly receive precipitation but receive water input from the upper runoff grid layer or other surrounding surfaces are divided into runoff grid layers. Based on water transfer on the vertical three-dimensional structure, the runoff within and between the surface layers is calculated. S53: Irregular grid cells corresponding to regions that have both runoff generation and runoff collection functions are divided into runoff generation and runoff collection grid layers. Some regions 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 of the same layer. By coupling the calculations of runoff generation and runoff collection, the total water volume and flow characteristics of the runoff generation and runoff collection grid layer are obtained.

5. The method for urban multi-granularity catchment area delineation based on high-precision three-dimensional data according to claim 1, characterized in that, The multi-granularity three-dimensional catchment area delineation in S6 includes the following steps: S61: Building elements are precisely defined by a macro boundary layer. For single-layer roofs, the elevation of the grid is modified by a flow generation grid layer. For double-layer roofs, the top roof uses a flow generation grid layer to finely depict the roof water distribution line to simulate the flow, distribution, accumulation, and drainage path of rainwater on the top roof. The second-layer roof, which is blocked by the top roof, receives water from the top roof or the surrounding grids of the same layer through a flow collection grid layer. The unblocked part of the second-layer roof takes into account the combined effect of direct rainfall input and flow collection from the top roof through a flow generation and collection grid layer. S62: Road elements describe the road edge contour through a macro boundary layer to delineate the boundary with the surrounding ordinary ground, and finely depict the elevation undulations inside the road through a runoff-generating grid layer to simulate the evolution path of real water flow. S63: The overpass elements describe the edge contour of the bridge deck through a macro boundary layer. The bridge deck uses a runoff generation grid layer to accurately depict the details of slope changes. Under the bridge, the runoff generation grid layer receives the water volume vertically transmitted from the upper bridge deck runoff generation grid layer as well as the lateral runoff from the surrounding ground surface. S64: Ordinary ground features describe the ground edge contour through a macro boundary layer and capture terrain changes through a runoff raster layer; S65: Based on the geometric morphological characteristics of each element, assign the macro boundary layer boundary control coefficient, impermeability percentage, and depression water storage parameters; S66: Based on the micro-hydrological processes of each element, parameters such as infiltration coefficient, surface friction coefficient, and evaporation coefficient are assigned to the detailed grid layer; S67: Integrate macroscopic boundary layer and detailed grid layer data of all hydrogeographic elements to form a multi-granular three-dimensional catchment area division.

6. The method for urban multi-granularity catchment area delineation 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: In response to the phenomena of stratified rainfall reception, asynchronous superposition of runoff, and vertical dynamic coupling caused by the three-dimensional structure of multi-story building roofs and overpasses, a four-dimensional rainfall intensity matrix with four collaborative dimensions including planar dimension, vertical dimension, time dimension, and feature dimension is introduced. S512: Based on the time span of the rainfall data, the rainfall data is evenly divided into segments 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 time point into a rainfall raster of the study area at the corresponding time point. Assign a single rainfall value to small-scale areas and use spline interpolation for large-scale areas. S514: For the multi-layered rain-receiving structure of different hydrogeographic elements, dynamically adjust the rainfall grid values ​​of different rain-receiving layers at each moment. Except for the secondary roof part and the space under the bridge that are blocked by the top roof, the corresponding grid rainfall intensity is adjusted to 0. Other elements do not need to be adjusted. Finally, the rainfall grids of different rain-receiving layers are spliced ​​into a complete three-dimensional rainfall distribution. S515: Repeatedly execute S512 and S513 to obtain rainfall rasters at n time points, and superimpose the rainfall rasters at all time points into a four-dimensional matrix of rainfall intensity in the study area. Finally, store the four-dimensional matrix of rainfall intensity in NetCDF format.

7. The method for urban multi-granularity catchment area delineation based on high-precision three-dimensional data according to claim 5, characterized in that, The flow calculation in S52 includes the following steps: S521: Intra-layer grid runoff calculation in surface runoff, based on the irregular grid structure in the detailed grid layer, using the two-dimensional shallow water equation, calculates the exchange of water flow between grids according to the water level difference between adjacent grids. The specific calculation formula is as follows: In the above formula: t represents time. Here, h is the flow vector, and h is the water depth. and These represent the flow velocities in the horizontal and vertical directions, respectively. Let be the elevation of the riverbed bottom, u and v be the flow velocities in the transverse and longitudinal directions, respectively, n be the Manning coefficient, and f and g be the water flux vectors in the transverse and longitudinal directions, respectively. quantity, For source terms, including S is the source term vector of runoff. b Let be the source term vector of the bottom slope. Let n be the source term vector of friction, and n be the Manning roughness coefficient. S522: Interlayer grid runoff calculation in surface runoff uses the common edge of the vertical overlap between the runoff-generating grid layer and the runoff-collecting grid layer to perform interlayer water exchange. The specific calculation formula is as follows: In the above formula, It is the water flux at the left boundary of the grid along the horizontal axis. Water flux at the boundary along the vertical axis of the grid. For time step, and These are the dimensions of the grid cells in the horizontal and vertical directions, respectively. and These are the flow rates in the horizontal and vertical directions, respectively. These are the row and column indices of the grid cells. yes Grid cell area, It is the water depth after the current grid water exchange.

8. An electronic device comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, characterized in that, When the processor executes the program, it implements a method for urban multi-granularity catchment area delineation based on high-precision three-dimensional data as described in any one of claims 1 to 7.

9. A computer-readable storage medium storing computer instructions thereon, characterized in that, When executed by the processor, the computer instructions implement a method for urban multi-granularity catchment area delineation based on high-precision three-dimensional data as described in any one of claims 1-7.

Images