Land utilization planning auxiliary method and system based on 3D terrain model

Abstract

The invention relates to the technical field of three-dimensional geographic modeling, in particular to a land utilization planning auxiliary method and system based on a 3D terrain model, and the method comprises the following steps: obtaining point cloud and elevation information, converting the information into a terrain curved surface, discretizing the information into a voxel grid, analyzing a planning index, converting the planning index into a geometric constraint boundary, and traversing to construct a planning voxel field. Calculating earthwork digging and filling requirements and sunshine hours to generate an evaluation matrix, correcting the building voxel distribution and finely adjusting the terrain elevation. According to the method, terrain refined expression is realized by using an octree segmentation algorithm, indexes such as a plot ratio are converted into three-dimensional space geometric constraints, site suitability is evaluated through connectivity analysis, and building layout and earthwork balance are automatically optimized by adopting morphological operation, so that the problems of lack of space intuition and error accumulation in traditional two-dimensional planning are solved; the planning calculation precision and decision scientificity under the complex terrain are remarkably improved, and dynamic adjustment and visual expression of the scheme are achieved.

Description

Technical Field

[0001] This invention relates to the field of three-dimensional geographic modeling technology, and in particular to a land use planning auxiliary method and system based on a 3D terrain model. Background Technology

[0002] The field of 3D geographic modeling technology involves using computer graphics to construct, represent, and process 3D digital models of real-world geographic environments. Through the combination of digital elevation models and texture mapping technology, it achieves virtual reconstruction and visualization of landforms, feature distribution, and spatial structure. Traditional land use planning support methods rely primarily on a combination of 2D planar map data and manual field surveys. Planners determine land attributes by overlaying and analyzing 2D vector data from different layers and then use manually measured geographic parameters such as slope and aspect to classify land parcels and calculate indicators.

[0003] Traditional land use planning mainly relies on a combination of two-dimensional planar map data and manual on-site surveys. Planners determine land attributes by overlaying and analyzing two-dimensional vector data from different layers, and divide plots into functional areas and calculate indicators based on manually measured slope and aspect geographical parameters. Due to the lack of intuitive support from three-dimensional spatial information, it is difficult to accurately control the impact of terrain elevation changes on building layout in planning and design in areas with complex terrain. Furthermore, manual measurement data has the potential for error accumulation, making it impossible to adjust planning schemes and terrain conditions in real time. This reduces the scientific rigor and practicality of land use planning schemes, making it difficult to meet the current growing demand for refined spatial control. Summary of the Invention

[0004] The purpose of this invention is to address the shortcomings of existing technologies by proposing a land use planning auxiliary method and system based on a 3D terrain model.

[0005] To achieve the above objectives, the present invention adopts the following technical solution: a land use planning auxiliary method based on a 3D terrain model, comprising the following steps:

[0006] S1: Obtain raw point cloud data and digital elevation information, use interpolation logic to transform the raw point cloud data and digital elevation information into a three-dimensional terrain surface, use an octree segmentation algorithm to discretize the three-dimensional terrain surface into a three-dimensional voxel grid, calculate coordinates, elevation and slope, and generate a terrain voxel dataset.

[0007] S2: Analyze the planning indicators to extract the upper limit of floor area ratio, building density and building height limit parameters, convert the upper limit of floor area ratio, building density and building height limit parameters into geometric constraint boundaries, traverse the terrain voxel dataset, detect whether the voxel nodes are located within the geometric constraint boundaries and mark the status, and construct the planning voxel field;

[0008] S3: Perform connectivity analysis on the set of constructable voxels for the planned voxel field, calculate the neighborhood elevation difference and calculate the earthwork excavation and filling demand in combination with the site elevation, calculate the cumulative sunshine hours using the sunshine analysis logic, and generate the earthwork and sunshine evaluation matrix.

[0009] S4: Based on the earthwork and sunshine evaluation matrix, identify areas with insufficient sunshine and excessive earthwork. Use morphological erosion and expansion calculations to correct the building voxel distribution and fine-tune the terrain voxel elevation to optimize the earthwork balance. Perform visualization conversion on the optimized voxel model to generate building layout and site elevation decision data.

[0010] As a further aspect of the present invention, step S1 specifically comprises:

[0011] S11: Acquire the original point cloud data including site spatial coordinates and reflection intensity, and simultaneously acquire the digital elevation information covering the same area. Use the inverse distance weighted interpolation algorithm to fuse discrete point data and construct the three-dimensional terrain surface.

[0012] S12: Set a minimum voxel resolution threshold, use an octree recursive segmentation strategy to hierarchically divide the bounding box occupied by the three-dimensional terrain surface, stop segmentation when the spatial size of the child node reaches the resolution threshold, and establish the spatial unit in the leaf node as a three-dimensional voxel mesh.

[0013] S13: Extract the geometric center of the surface fragments included in each voxel grid as voxel coordinates, calculate the average Z value of the fragment vertices as voxel elevation, and determine the slope by analyzing the surface tilt angle based on the fragment normal vector. Store the voxel objects including coordinates, elevation and slope attributes in a structured manner to generate the terrain voxel dataset.

[0014] As a further aspect of the present invention, step S2 specifically includes:

[0015] S21: Analyze land use planning control documents, extract key numerical indicators, and identify the plot ratio limit for land development intensity, building density for building footprint coverage, and building height limit parameters for vertical space restrictions;

[0016] S22: Construct a three-dimensional spatial envelope based on the land parcel boundary and the building height limit parameters, and modify the bottom boundary of the envelope by combining the urban road setback distance to generate the geometric constraint boundary;

[0017] S23: Traverse the voxel units in the terrain voxel dataset one by one, use the ray crossing method to detect the spatial positional relationship between the voxel center point and the geometric constraint boundary, mark the voxels located inside the boundary as buildable, and mark the voxels located outside the boundary or intersecting with it as non-buildable, thereby constructing the planned voxel field.

[0018] As a further aspect of the present invention, step S3 specifically comprises:

[0019] S31: Use a breadth-first search algorithm to traverse the voxels marked as buildable in the planned voxel field, identify voxel clusters with adjacent relationships to form connected components, remove isolated connected components with a volume smaller than the preset minimum building unit threshold, and establish an effective set of buildable voxels.

[0020] S32: Obtain the preset site design elevation, calculate the vertical difference between the bottom elevation of each voxel in the effective buildable voxel set and the site design elevation, calculate the individual excavation and filling volume based on the voxel bottom area integral, and sum them up to obtain the earthwork excavation and filling requirements.

[0021] S33: Calculate the solar altitude angle and azimuth angle based on the latitude and longitude of the project location and the sampling time step. Perform shadow ray projection detection on each voxel. Calculate the sum of time segments that are not blocked by the surrounding terrain or neighboring voxels to obtain the cumulative sunshine hours. Combine the data dimension of the earthwork excavation and filling demand to generate the earthwork and sunshine evaluation matrix.

[0022] As a further aspect of the present invention, step S4 specifically comprises:

[0023] S41: Traverse the earthwork and sunshine evaluation matrix, filter the shadow voxels whose cumulative sunshine hours are lower than the national standard threshold, and the terrain voxels whose absolute value of earthwork excavation and filling demand exceeds the engineering cost limit, and identify and mark them as areas with insufficient sunshine and excessive earthwork that need to be optimized.

[0024] S42: Perform morphological erosion operation on the building voxel group marked as insufficient sunlight area to reduce the volume of the shading source, perform smooth expansion operation on the terrain voxel in the earthwork over-limit area to reshape the local terrain slope, and iteratively adjust the voxel state until the preset evaluation convergence condition is met to complete the building voxel distribution correction and terrain voxel elevation fine adjustment.

[0025] S43: Read the optimized voxel model state data, convert the voxel mesh into a polygon mesh entity, use the rendering engine to map material and lighting properties, and output the building layout and site elevation decision data.

[0026] As a further aspect of the present invention, the slope calculation process includes:

[0027] Obtain the coordinates of the three vertices of the triangular facet inside the voxel mesh, construct the normal vector of the triangular facet using vector product operation, and calculate the cosine value of the angle between the normal vector and the vertical axis of the world coordinate system.

[0028] The tilt angle of the land surface relative to the horizontal plane is obtained by analyzing the inverse cosine function, and the angle value is normalized and mapped to a preset slope grading interval. The grading result is written into the terrain voxel dataset as the slope attribute value of the voxel.

[0029] As a further aspect of the present invention, the process of generating the earthwork and sunshine evaluation matrix includes:

[0030] Obtain the numerical matrix of the earthwork excavation and filling demand and the cumulative sunshine hours, introduce the engineering cost weighting factor and the living comfort weighting factor, calculate the development suitability score of each spatial unit according to the weighted comprehensive evaluation model, and construct the comprehensive evaluation matrix calculation logic as shown below: ;

[0031] in, Representing coordinates A comprehensive score assessing the suitability of phytoestrogens for development. The measured cumulative sunshine duration representing the location of a voxel. This represents the minimum threshold number of hours of sunshine specified by the sunshine standard. The actual elevation value representing the voxel. This represents the benchmark design elevation set in the plan. This represents the maximum allowable elevation difference of the site. This represents the weighting factor for residential comfort. This represents the weighting factor for the project cost.

[0032] As a further aspect of the present invention, the calculation process for accumulated sunshine hours includes:

[0033] Using the winter solstice as the base date, a solar position vector sequence is constructed according to a preset time sampling interval, and a virtual ray is emitted from the geometric center of each voxel to be detected in the direction indicated by the solar position vector sequence;

[0034] The system detects whether there are any blocking voxels marked as non-transparent on the virtual ray path. If there are no obstructions on the path, the system accumulates the sunshine duration corresponding to the time sampling interval. After traversing all sampling times throughout the day, the system summarizes the total effective sunshine duration of the voxels as the accumulated sunshine duration.

[0035] As a further aspect of the present invention, the execution process of morphological erosion and dilation operations includes:

[0036] Define a structural element that includes a central voxel and its 26 neighboring voxels. During the erosion stage, if all voxels within the coverage area of ​​the structural element do not meet the retention conditions, the state of the central voxel is flipped from buildable to non-buildable, thereby removing edge noise from the building voxel group.

[0037] During the expansion phase, if there is a voxel in the neighborhood of the central voxel that satisfies the target elevation continuity condition, the state of the central voxel is assimilated into the state of the neighboring voxel, thereby filling the gaps inside the building layout and smoothing the abrupt areas of the terrain surface.

[0038] A land use planning assistance system based on a 3D terrain model, the system being used to implement the aforementioned land use planning assistance method based on a 3D terrain model, the system comprising:

[0039] The terrain discretization construction module is used to acquire raw point cloud data and digital elevation information, generate a three-dimensional terrain surface using interpolation logic, and calculate coordinates, elevation and slope using an octree segmentation algorithm to output a terrain voxel dataset.

[0040] The planning constraint mapping module is used to extract the upper limit of floor area ratio, building density and building height limit parameters, convert them into geometric constraint boundaries, and perform traversal detection and state labeling on the terrain voxel dataset to generate a planning voxel field.

[0041] The site suitability analysis module is used to perform connectivity analysis on the planned voxel field, calculate the earthwork excavation and filling demand and cumulative sunshine hours, and generate an earthwork and sunshine evaluation matrix by integrating each indicator.

[0042] The spatial morphology optimization module is used to identify areas with insufficient sunlight and excessive earthwork, and uses morphological erosion and dilation calculations to correct voxel state and elevation information, visualize the model, and output building layout and site elevation decision data.

[0043] Compared with the prior art, the advantages and positive effects of the present invention are as follows:

[0044] In this invention, a three-dimensional terrain surface is constructed by acquiring original point cloud and elevation information and discretized into a voxel mesh. The planning indicators are analyzed and transformed into geometric constraint boundaries, and the voxel node status is traversed and detected. A planning voxel field containing coordinate and slope information is constructed. Connectivity analysis is performed on the buildable voxels, and earthwork excavation and filling requirements are calculated in combination with the site elevation. An evaluation matrix is ​​generated using solar radiation analysis logic. Based on the evaluation results, morphological erosion and dilation operations are used to correct the building voxel distribution. The elevation of the terrain voxels is fine-tuned to optimize the earthwork balance, thereby achieving accurate planning and automatic optimization in three-dimensional space and effectively improving the scientific nature and decision-making efficiency of land use planning schemes. Attached Figure Description

[0045] Figure 1 This is the main flowchart of the land use planning auxiliary method of the present invention;

[0046] Figure 2 This is a flowchart illustrating the generation process of the terrain voxel dataset for this invention.

[0047] Figure 3 This is a flowchart illustrating the construction process of the voxel field in this invention.

[0048] Figure 4 This is a flowchart of the earthwork and sunshine evaluation matrix generation process of the present invention;

[0049] Figure 5 This is a flowchart of the data output for building layout and site elevation decision-making in this invention. Detailed Implementation

[0050] To make the objectives, technical solutions, and advantages of this invention clearer, the software-based technical solution is described in detail below with reference to system architecture diagrams and embodiments. It should be understood that the specific embodiments described herein are only for explaining the technical solutions of this invention and do not constitute a limitation on the scope of protection.

[0051] In the description of this invention, the system architecture relationships or data processing flows indicated by terms such as "layer," "module," "interface," "data flow," "client," and "server" are all defined based on the architecture diagram or flowchart corresponding to the embodiments. This way of describing is only used to clearly illustrate the logical relationships between the elements in the technical solution, and not to limit the physical deployment form. The term "multiple" includes two or more technical units, including but not limited to multiple data nodes, processing threads, service instances, or functional components and other scalable elements. The specific number is determined according to the actual business scenario and needs to be specifically specified.

[0052] Please see Figure 1 and Figure 2 This invention provides a technical solution: a land use planning auxiliary method based on a 3D terrain model, comprising the following steps:

[0053] S1: Obtain raw point cloud data and digital elevation information, use interpolation logic to transform the raw point cloud data and digital elevation information into a three-dimensional terrain surface, use an octree segmentation algorithm to discretize the three-dimensional terrain surface into a three-dimensional voxel mesh, calculate coordinates, elevation and slope, and generate a terrain voxel dataset.

[0054] S11: Acquire raw point cloud data including site spatial coordinates and reflection intensity, and simultaneously acquire digital elevation information covering the same area. Use inverse distance weighted interpolation algorithm to fuse discrete point data and construct a three-dimensional terrain surface.

[0055] S12: Set the minimum voxel resolution threshold, use an octree recursive segmentation strategy to hierarchically divide the bounding box occupied by the three-dimensional terrain surface, stop segmentation when the spatial size of the child node reaches the resolution threshold, and establish the spatial unit in the leaf node as a three-dimensional voxel mesh.

[0056] S13: Extract the geometric center of the surface fragments included in each voxel grid as voxel coordinates, calculate the average Z value of the fragment vertices as voxel elevation, and determine the slope by analyzing the surface tilt angle based on the fragment normal vector. Store the voxel objects including coordinates, elevation and slope attributes in a structured manner to generate a terrain voxel dataset.

[0057] The process of calculating slope includes:

[0058] Obtain the coordinates of the three vertices of the triangular facet inside the voxel mesh, construct the normal vector of the triangular facet using vector product operation, and calculate the cosine value of the angle between the normal vector and the vertical axis of the world coordinate system.

[0059] The tilt angle of the land surface relative to the horizontal plane is obtained by analyzing the inverse cosine function, and the angle value is normalized and mapped to a preset slope grading interval. The grading result is written into the terrain voxel dataset as the slope attribute value of the voxel.

[0060] First, establish a high-precision terrain data acquisition environment. Select the project site as the target area and deploy an airborne LiDAR scanning system, such as the Riegl VZ-400i 3D laser scanner, setting the laser pulse repetition frequency to [value missing]. The scanning angular resolution is This is to obtain raw point cloud data of the Earth's surface. Raw point cloud data not only contains spatial three-dimensional coordinates... It also includes reflection intensity information, and the data collection density needs to be controlled at no less than [a certain value]. To ensure the capture of micro-features of the land surface, digital elevation information (DEM) covering the same area is obtained by simultaneously measuring control points using real-time dynamic carrier phase differential technology (RTK-GPS), with both horizontal and vertical accuracy controlled within a specified range. Within.

[0061] In the data preprocessing stage, statistical filtering is performed on the raw point cloud data to remove floating noise. The number of neighboring points to search is set. Standard deviation multiple threshold Calculate the distance from each point to its The average distance between the neighboring points is used to eliminate outliers whose average distance is greater than a multiple of the global average distance plus the standard deviation. Then, the fusion logic described in S11 is executed, and a 3D terrain surface is constructed using an inverse distance-weighted interpolation (IDW) algorithm. Specifically, the interpolation search radius is set. Distance weighted exponent For the target grid point Search its radius Effective discrete points within the range ,calculate Points and each Euclidean distance of points ,Will and Sum of products divided by The sum of these factors allows for the analytical extraction of a continuous three-dimensional terrain surface. .

[0062] The aforementioned statistical filtering refers to a point cloud denoising algorithm. It calculates the average distance distribution between the query point and its neighboring point set, assuming that the distance distribution conforms to a Gaussian distribution. Based on the mean and standard deviation, it defines points at the tail of the distribution as outliers and removes them.

[0063] Next, the discretization process described in S12 is executed. First, the three-dimensional terrain surface is calculated. The axis-aligned bounding box is used to determine its maximum and minimum boundary points in the world coordinate system. and Set the minimum voxel resolution threshold. This serves as the termination condition for octree partitioning. The root node of the octree is initialized as the entire bounding box. A recursive partitioning strategy is used: check if the spatial region bounded by the current node contains terrain surface fragments; if it does, and the node's spatial size (side length) is greater than... Then the node along Each of the three axes is divided into 8 child nodes. This process is repeated until the child node size is equal to or smaller than the resolution threshold, at which point the segmentation stops. At this point, all spatial cells located at leaf nodes and containing terrain information are established as a 3D voxel mesh. .

[0064] The aforementioned axis-aligned bounding box refers to a cuboid whose faces are parallel to the coordinate axes, used to enclose a three-dimensional geometric object. Its boundary range is determined by the minimum and maximum coordinate values ​​of the object in the X, Y, and Z axes.

[0065] Then, the attribute calculation and structured storage described in S13 are performed. For each established voxel mesh... Extract the set of terrain surface fragments enclosed within it. Calculate the geometric centroid of all fragment vertices in this set and use it as the voxel coordinate system. Simultaneously, iterate through all sampling points within the fragment and calculate its... The arithmetic mean of the axis coordinates is assigned as the voxel elevation. .

[0066] For the calculation of slope properties, the following specific vector analysis process is performed: Obtain the three non-collinear vertices of the triangular facets inside the voxel mesh. , , Constructing vectors with vector The normal vector of the triangular facet is calculated using the vector product operation. Then, the normal vector is calculated. Perpendicular axis vector of world coordinate system cosine value of the angle between The specific calculation is as follows: and The dot product is divided by the product of the magnitudes of the two. Based on the inverse cosine function... The angle of inclination of the Earth's surface relative to the horizontal plane is obtained through analysis.

[0067] To achieve data structuring, angle values Normalization is mapped to a preset slope gradation interval. The preset interval is divided into: flat slope gentle slope Middle slope ,steep slope and steep slopes The slopes are coded into levels 1 through 5. The calculated grading results are used as the slope attribute values. Write to the voxel object. Ultimately, this will contain a unique identifier (ID), voxel coordinates, and elevation. and slope The voxel objects are serialized and stored as OpenVDB format files to generate a terrain voxel dataset.

[0068] Table 1 shows some voxel attribute examples generated after processing of collected point cloud data, verifying the accuracy of data conversion.

[0069] Table 1. Attribute Sampling Table for Terrain Voxel Dataset:

[0070]

[0071] Referring to Table 1, the data in the table shows that this embodiment can accurately capture subtle changes in terrain undulations and transform complex geometric surfaces into discrete voxels with semantic attributes, providing a standardized data foundation for subsequent planning and analysis.

[0072] Please see Figure 1 and Figure 3 S2: Extract the floor area ratio limit, building density and building height limit parameters from the planning indicators, transform the floor area ratio limit, building density and building height limit parameters into geometric constraint boundaries, traverse the terrain voxel dataset, detect whether the voxel nodes are within the geometric constraint boundaries and mark the status, and construct the planning voxel field.

[0073] S21: Analyze land use planning control documents, extract key numerical indicators, and identify parameters related to plot development intensity, building density, building footprint coverage, and building height restrictions.

[0074] S22: Construct a three-dimensional spatial envelope based on the land parcel boundary and building height limit parameters, and modify the bottom boundary of the envelope by combining the urban road setback distance to generate a geometric constraint boundary.

[0075] S23: Iterate through the voxel cells in the terrain voxel dataset one by one, use the ray crossing method to detect the spatial positional relationship between the voxel center point and the geometric constraint boundary, mark the voxels located inside the boundary as buildable, and mark the voxels located outside the boundary or intersecting with it as non-buildable, thereby constructing the planning voxel field.

[0076] First, the indicator analysis described in S21 is executed. The control detailed planning documents for the project area are accessed, and key values ​​regarding the development intensity of the land parcel are extracted through the attribute query interface. For example, in this embodiment of a mixed-use commercial and residential land project, the maximum floor area ratio is extracted. Building density limit and building height limit parameters Simultaneously, the land parcel boundary line (land parcel boundary coordinate sequence) and surrounding road classification information are read.

[0077] Next, execute S22 to construct the geometric constraint boundary. First, construct the polygon base based on the plot boundary coordinate sequence. According to urban planning management technical regulations, adjust the envelope base boundary in conjunction with urban road setback distances. If the plot's north side is adjacent to a main road, set the setback distance as follows: The eastern side adjacent to the secondary road should yield. A polygon shrinkage algorithm is used to shrink the red-line polygon inward by a corresponding distance to generate the actual buildable land surface. Subsequently, based on the building height restriction parameters... ,by As the bottom surface, along The axis is stretched vertically upwards, and the construction height is [missing information]. The three-dimensional spatial envelope is then generated. Furthermore, if there are solar radiation cone restrictions or airport airspace restrictions, a corresponding oblique surface needs to be introduced to perform a Boolean subtraction operation on this envelope to remove the excess portion, ultimately generating a precise geometric constraint boundary. .

[0078] Based on this, execute S23 to construct the planned voxel field. Load the terrain voxel dataset generated in S1, and iterate through the voxel units within it one by one. Extract the center coordinates of each voxel. The ray-crossing method is used to detect its relationship with the geometric constraint boundary. Spatial location relationship. The specific operation is: from... Emit a ray in any direction and count the ray's intersection with the target direction. The number of intersections of the surface mesh. If the number of intersections is odd, the voxel is determined to be inside the boundary; if the number is even, it is determined to be outside.

[0079] The aforementioned ray intersection method refers to a computational geometry algorithm that determines whether a point is located inside a polygon by emitting rays from a test point and counting the number of intersections between rays and the polygon boundary, based on the Odd-Even Rule.

[0080] Voxel status is marked based on the determination result:

[0081] 1. Constructable state: When the voxel center is located Internally, and the slope attribute value of the voxel is less than the preset construction slope threshold (e.g., When the slope is not steep, mark the status bit. ;

[0082] 2. Non-buildable state: When the voxel is located in External, or intersecting with the boundary but with most of the volume on the outside, or although inside but with a steep slope (e.g.) When construction is deemed unsuitable, a status flag is set. .

[0083] To handle "suspended" voxels (i.e., air voxels located above the terrain but within the height limit), this embodiment expands the planned voxel field into a full-space voxel field. Using the terrain voxels as a base, empty voxels are filled upwards until the height limit is reached. For these empty voxels located above the terrain and within the envelope, ray detection is also performed. If they are inside the envelope, they are marked as potential building voxels, thus forming a complete planning voxel field that includes both the terrain layer and the potential building layer. This process ensures the continuity and compliance of the planning space in three dimensions.

[0084] Please see Figure 1 and Figure 4 S3: Perform connectivity analysis on the set of constructable voxels for the planned voxel field, calculate the elevation difference of the neighborhood and calculate the earthwork excavation and filling demand by combining the site elevation, calculate the cumulative sunshine hours using the sunshine analysis logic, and generate the earthwork and sunshine evaluation matrix.

[0085] S31: Use a breadth-first search algorithm to traverse the voxels marked as buildable in the planning voxel field, identify voxel clusters with adjacent relationships to form connected components, remove isolated connected components with a volume smaller than the preset minimum building unit threshold, and establish an effective set of buildable voxels.

[0086] S32: Obtain the preset site design elevation, calculate the vertical difference between the bottom elevation of each voxel in the effective buildable voxel set and the site design elevation, calculate the individual excavation and filling volume based on the voxel bottom area integral, and sum them up to obtain the earthwork excavation and filling requirements.

[0087] S33: Calculate the solar altitude angle and azimuth angle based on the latitude and longitude of the project location and the sampling time step. Perform shadow ray projection detection on each voxel and sum the time segments that are not blocked by the surrounding terrain or neighboring voxels to obtain the cumulative sunshine hours. Combine the data dimension of earthwork excavation and filling demand to generate an earthwork and sunshine evaluation matrix.

[0088] The calculation process for cumulative sunshine hours includes:

[0089] Using the winter solstice as the base date, a solar position vector sequence is constructed according to a preset time sampling interval, and virtual rays are emitted from the geometric center of each voxel to be detected in the direction indicated by the solar position vector sequence;

[0090] The system detects whether there are any blocking voxels marked as non-transparent on the virtual ray path. If there are no obstructions on the path, it accumulates the sunshine duration corresponding to the time sampling interval. After traversing all sampling times throughout the day, it summarizes the total effective sunshine duration of the voxels as the cumulative sunshine duration.

[0091] The process of generating the earthwork and sunshine evaluation matrix includes:

[0092] A numerical matrix of earthwork excavation and filling demand and cumulative sunshine hours is obtained. Engineering cost weighting factors and residential comfort weighting factors are introduced. Based on a weighted comprehensive evaluation model, the development suitability score of each spatial unit is calculated, and the comprehensive evaluation matrix calculation logic is constructed as shown below:

[0093] ;

[0094] in, Representing coordinates A comprehensive score assessing the suitability of phytoestrogens for development. The measured cumulative sunshine duration representing the location of a voxel. This represents the minimum threshold number of hours of sunshine specified by the sunshine standard. The actual elevation value representing the voxel. This represents the benchmark design elevation set in the plan. This represents the maximum allowable elevation difference of the site. Represents the weighting factor of living comfort. This represents the weighting factor for project cost.

[0095] First, perform the connectivity analysis described in S31. Use a breadth-first search (BFS) algorithm to traverse the planned voxel field. Initialize all constructible voxels ( The access flag for any unvisited buildable voxel is set to False. Starting with the first unvisited buildable voxel, a queue is created, and voxels are pushed into the queue. In the loop, the voxel at the head of the queue is popped, and all voxels within its 26-neighborhood (including face, edge, and corner adjacencies) are searched. If a neighboring voxel is also buildable and unvisited, it is added to the queue and marked as visited, belonging to the same connected component ID. When the queue is empty, a connected component search is complete. This process is repeated until all nodes are traversed. Subsequently, the number of voxels contained in each connected component is counted, combined with the voxel volume (e.g., ...). Calculate the physical volume. Set the preset minimum building unit threshold to [value]. Remove isolated connected components with a volume smaller than the threshold to establish a valid set of constructible voxels. .

[0096] The breadth-first search algorithm mentioned above is an algorithm used to traverse or search a tree or graph. It starts from the root node and traverses the nodes of the tree along the width of the tree. The algorithm terminates when all nodes have been visited.

[0097] Next, execute S32 to calculate the required earthwork excavation and filling volume. Obtain the baseline design elevation set in the plan. In this embodiment, it is set to Traversal Voxels belonging to the surface layer were used to extract their bottom elevation. (Right now Subtract half of the diagonal vertical component). Calculate the vertical difference. .like This indicates that excavation is required, and the excavation volume is specified. ;like This indicates that fill is required, and the fill volume is specified. .in voxel base area ( The total earthwork excavation and filling requirements are obtained by summing the calculation results of all surface voxels.

[0098] Simultaneously execute S33 and sunshine calculation processes. Set latitude and longitude as (North Latitude) East longitude Time step Using the winter solstice as the base date, the time range is from 09:00 to 15:00. The Solar Position Algorithm (PSA) is used to calculate each sampling time. solar altitude angle With azimuth Construct a sequence of solar position vectors .for Each voxel in the array, from its geometric center along A virtual ray is emitted in the opposite direction. Ray tracing is used to detect the ray path: if the path intersects with other non-transparent voxels (terrain or identified surrounding building voxels), it is considered "occlusion"; otherwise, it is considered "lit". The time step of all "lit" moments is summed to obtain the cumulative sunshine hours. .

[0099] Finally, the earthwork and sunlight evaluation matrices are generated. A weighted comprehensive evaluation model is introduced, and weighting factors for residential comfort are set. Engineering cost weighting factor Set a minimum threshold for sunshine hours. Maximum allowable elevation difference of the site Regarding coordinates For the voxel at that location, substitute the parameters obtained above into the formula:

[0100] ;

[0101] in, Representing coordinates Comprehensive score for the suitability of phytoestrogens for development; Measured cumulative sunshine hours representing the location of a voxel; This represents the minimum threshold number of hours of sunshine specified by the sunshine standard. The actual elevation value representing the voxel; Represents the benchmark design elevation set in the plan; This represents the maximum permissible elevation difference for the site used in normalization. This represents the weighting factor for residential comfort, used to adjust the importance of sunlight factors in the evaluation; This represents the weighting factor of the project cost, used to adjust the importance of earthwork cost factors in the evaluation.

[0102] A practical example is provided: Select a voxel. Its coordinates are Actual cumulative sunshine hours (i.e., 90 minutes), actual elevation Known .

[0103] Substitute the parameters into the formula:

[0104] First item (sunlight gain):

[0105] ;

[0106] Second item (penalty for earthwork):

[0107] ;

[0108] Final score: .

[0109] The results indicate that although this voxel has some sunshine, its overall suitability score is low due to insufficient sunshine and increased earthwork costs caused by a certain elevation difference. The system traversed all voxels to complete the calculation and generated earthwork and sunshine evaluation matrices. Table 2 shows a comparison of the calculation results for different voxels.

[0110] Table 2 Sampling Table of Earthwork and Sunlight Evaluation Matrix Data

[0111] See Table 2. The data in this table shows the differences in the development value of voxels under different sunlight and topographic conditions. High-scoring areas (such as A) are the preferred areas for subsequent building layout.

[0112] Please see Figure 1 and Figure 5 S4: Based on the earthwork and sunshine evaluation matrix, identify areas with insufficient sunshine and excessive earthwork. Use morphological erosion and expansion calculations to correct the building voxel distribution and fine-tune the terrain voxel elevation to optimize the earthwork balance. Perform visualization transformation on the optimized voxel model to generate building layout and site elevation decision data.

[0113] S41: Traverse the earthwork and sunshine evaluation matrix, filter shadow voxels whose cumulative sunshine hours are lower than the national standard threshold, and topographic voxels whose absolute value of earthwork excavation and filling demand exceeds the engineering cost limit, identify and mark areas with insufficient sunshine and excessive earthwork that need to be optimized.

[0114] S42: Perform morphological erosion operation on the building voxel group marked as insufficient sunlight area to reduce the volume of the shading source, perform smooth expansion operation on the terrain voxel in the earthwork over-limit area to reshape the local terrain slope, and iteratively adjust the voxel state until the preset evaluation convergence condition is met to complete the building voxel distribution correction and terrain voxel elevation fine adjustment.

[0115] The execution process of morphological erosion and dilation operations includes:

[0116] Define a structural element that includes a central voxel and its 26 neighboring voxels. During the erosion stage, if all voxels within the coverage area of ​​the structural element do not meet the retention conditions, the state of the central voxel is flipped from buildable to non-buildable, thereby removing edge noise from the building voxel group.

[0117] During the expansion phase, if there is a voxel in the neighborhood of the central voxel that satisfies the target elevation continuity condition, the state of the central voxel is assimilated into the state of the neighboring voxel, thereby filling the gaps inside the building layout and smoothing the abrupt areas of the terrain surface.

[0118] S43: Read the optimized voxel model state data, convert the voxel mesh into a polygon mesh entity, use the rendering engine to map material and lighting properties, and output building layout and site elevation decision data.

[0119] First, perform the region identification described in S41. Iterate through the evaluation matrix generated in S3 and set the filtering criteria:

[0120] 1. Areas with insufficient sunlight: Filter (Right now And it is a set of voxels of the proposed building;

[0121] 2. Areas exceeding earthwork limits: Screening A set of terrain voxels (within the engineering cost limit). Voxels that meet the above conditions are labeled with specific tags for subsequent processing.

[0122] Next, morphological optimization of S42 is performed.

[0123] For areas with insufficient sunlight, morphological erosion is performed on the building voxel group that blocks sunlight. The structural element is defined as... The cube core. During the erosion phase, for the central voxels marked as buildable. The system checks 26 neighboring voxels within the coverage area of ​​its structural element. If the number of voxels in the neighborhood that are in a "non-buildable" state exceeds a preset threshold (e.g., 5), it indicates that the voxel is located at the edge of the building and may cause unnecessary occlusion or is too fragmented. In this case, [the system will proceed as follows]. The state is flipped to unbuildable. This operation removes edge noise from the building's volumetric elements, effectively reducing sharp edges and increasing light penetration.

[0124] The aforementioned morphological erosion operation refers to a mathematical morphological operation based on set theory. By using structuring elements to probe images or voxel models, it removes pixels or voxels at the boundaries of objects, thereby shrinking the object volume and eliminating small protrusions.

[0125] For areas with excessive earthwork, a smoothing expansion operation is performed to reshape the terrain. For the central voxels of elevation abrupt changes... If there exists a neighborhood that satisfies the "target elevation continuity condition" (i.e., the neighborhood average elevation is closer), If the voxels are modified, then modify them. The elevation value is such that it converges to the neighborhood average.

[0126] Specific example: Elevation of a certain voxel Its 8-neighborhood average elevation is The coefficient of thermal expansion is The corrected elevation is This process is repeated iteratively until the elevation change rate of all voxels is less than [a certain value]. This achieves terrain smoothing and earthwork balance.

[0127] Finally, the S43 visualization transformation is executed. The optimized voxel model state data is read, and the moving cubes algorithm is used to convert the voxel mesh into polygonal mesh entities. The algorithm traverses the voxel field, calculates isosurfaces based on the voxel vertex state values, and generates triangular patches. The generated mesh is loaded using a rendering engine (such as Unity3D or Unreal Engine), and materials are mapped according to the voxel slope attributes (grass material for flat ground, rock material for steep slopes) and lightmaps are mapped according to the sunlight attributes.

[0128] The aforementioned moving cube algorithm is a computer graphics algorithm used to extract isosurfaces from a three-dimensional scalar field (such as voxel data), determine the topological structure of isosurfaces within voxels using a lookup table method, and generate corresponding triangular patches, thereby reconstructing a three-dimensional surface model.

[0129] The output decision data includes:

[0130] 1. Building layout model: The final building volume OBJ file after sunlight optimization;

[0131] 2. Site Elevation Map: Topographic elevation CSV data optimized for earthwork balance, containing the suggested design elevation for each grid point. Through the above steps, an optimal land use planning scheme is ultimately formed that satisfies planning regulations while achieving maximum sunlight and minimum earthwork at the physical level.

[0132] A land use planning assistance system based on a 3D terrain model is provided. This system is used to execute the aforementioned land use planning assistance method based on a 3D terrain model. The system includes:

[0133] The terrain discretization construction module is used to acquire raw point cloud data and digital elevation information, generate a three-dimensional terrain surface using interpolation logic, and calculate coordinates, elevation and slope using an octree segmentation algorithm to output a terrain voxel dataset.

[0134] The planning constraint mapping module is used to extract the upper limit of floor area ratio, building density and building height limit parameters, convert them into geometric constraint boundaries, and perform traversal detection and state labeling on the terrain voxel dataset to generate a planning voxel field.

[0135] The site suitability analysis module is used to perform connectivity analysis on the planning voxel field, calculate the earthwork excavation and filling demand and cumulative sunshine hours, and generate an earthwork and sunshine evaluation matrix by integrating each indicator.

[0136] The spatial morphology optimization module is used to identify areas with insufficient sunlight and excessive earthwork. It uses morphological erosion and dilation calculations to correct voxel state and elevation information, visualizes the model, and outputs building layout and site elevation decision data.

[0137] The above embodiments illustrate preferred embodiments of the present invention. Any equivalent adjustments to the technical solution based on software engineering methods are within the scope of protection, including but not limited to: implementing algorithm logic using different programming languages, refactoring functional modules into services, adjusting data interaction protocols, and optimizing resource scheduling strategies. Any implementation scheme derived from reasonable modifications to the data processing flow, service call chain, or system architecture layer without departing from the core technology of the present invention should be considered within the protection scope defined by the technical solution of the present invention.

Claims

1. A 3D terrain model-based land use planning assistance method, characterized by, The method comprises the following steps: S1: obtaining original point cloud data and digital elevation information, converting the original point cloud data and the digital elevation information into a three-dimensional terrain surface by using interpolation logic, discretizing the three-dimensional terrain surface into a three-dimensional voxel grid by using an octree segmentation algorithm, calculating coordinates, elevations and slopes, and generating a terrain voxel dataset; S2: analyzing planning indicators to extract a floor area ratio upper limit, a building density and a building height limit parameter, converting the floor area ratio upper limit, the building density and the building height limit parameter into a geometric constraint boundary, traversing the terrain voxel dataset, detecting whether a voxel node is located within the geometric constraint boundary and marking a state, and constructing a planning voxel field; S3: performing connectivity analysis on a buildable voxel set of the planning voxel field, calculating a neighborhood elevation difference and combining site elevation to calculate earthwork excavation and filling demand, calculating cumulative sunshine hours by using sunshine analysis logic, and generating an earthwork and sunshine evaluation matrix; S4: identifying sunshine deficiency and earthwork over-limit areas according to the earthwork and sunshine evaluation matrix, modifying building voxel distribution and fine-tuning terrain voxel elevation to optimize earthwork balance by using morphological erosion and dilation operations, and generating building layout and site elevation decision data by visualizing the optimized voxel model.

2. The 3D terrain model-based land use planning assistance method according to claim 1, characterized by, The step S1 specifically comprises: S11: obtaining the original point cloud data including site spatial coordinates and reflection intensity, and synchronously obtaining the digital elevation information covering the same area, fusing discrete point data by using an inverse distance weighted interpolation algorithm, and constructing the three-dimensional terrain surface; S12: setting a minimum voxel resolution threshold, adopting an octree recursive segmentation strategy to hierarchically divide a space bounding box occupied by the three-dimensional terrain surface, stopping segmentation when the spatial size of a child node reaches the resolution threshold, and establishing a space unit at a leaf node as a three-dimensional voxel grid; S13: extracting the geometric center of a surface patch element included in each voxel grid as a voxel coordinate, calculating the average Z value of a patch vertex as a voxel elevation, determining a slope based on a patch normal vector, and structuring storage of a voxel object including coordinate, elevation and slope attributes, thereby generating the terrain voxel dataset.

3. The 3D terrain model-based land use planning assistance method according to claim 1, characterized by, The step S2 specifically comprises: S21: analyzing a land use planning control planning file, extracting key numerical indicators, identifying the floor area ratio upper limit related to plot development intensity, the building density related to building footprint coverage, and the building height limit parameter related to vertical space restriction; S22: constructing a three-dimensional space envelope body according to a plot red line range and the building height limit parameter, and combining a city road retreat distance to correct the envelope body bottom boundary, thereby generating the geometric constraint boundary; S23: traversing voxel units in the terrain voxel dataset one by one, detecting the spatial position relationship between a voxel center point and the geometric constraint boundary by using a ray intersection method, marking voxels located inside the boundary as buildable, marking voxels located outside the boundary or intersecting the boundary as non-buildable, and thereby constructing the planning voxel field. 4.The 3D terrain model based land use planning assistance method of claim 1, wherein, The step S3 specifically comprises: S31: A breadth-first search algorithm is used to traverse the voxels marked as buildable in the planning voxel field, identify voxel clusters with adjacent relationships to form connected domains, eliminate isolated connected domains with a volume smaller than a preset minimum building unit threshold, and establish an effective buildable voxel set; S32: A preset site design elevation is obtained, the vertical difference between the bottom elevation of each voxel in the effective buildable voxel set and the site design elevation is calculated, the single excavation and filling volume is calculated based on voxel bottom area integration, and the total earthwork excavation and filling demand is obtained by accumulation; S33: Based on the latitude and longitude of the project location and the sampling time step, the solar elevation angle and azimuth angle are calculated, shadow ray projection detection is performed on each voxel, the total sum of time segments that are not blocked by surrounding terrain or adjacent voxels is counted to obtain the cumulative sunshine hours, and the earthwork and sunshine evaluation matrix is generated in combination with the data dimension of the total earthwork excavation and filling demand. 5.The 3D terrain model based land use planning assistance method according to claim 1, wherein, The steps of S4 are specifically: S41: The earthwork and sunshine evaluation matrix is traversed, shadow voxels with a cumulative sunshine hours lower than the national standard threshold and terrain voxels with an absolute value of the total earthwork excavation and filling demand exceeding the engineering cost limit value are screened and marked as sunshine insufficient and earthwork over-limit areas that need to be optimized; S42: Morphological erosion operation is performed on the building voxel group marked as a sunshine insufficient area to reduce the blocking source volume, and smoothing dilation operation is performed on the terrain voxels in the earthwork over-limit area to reshape the local terrain slope, and the voxel state is adjusted through iteration until the preset evaluation convergence condition is met, thereby completing the distribution correction of the building voxels and the fine adjustment of the terrain voxel elevation; S43: The optimized voxel model state data is read, the voxel grid is converted into a polygon grid entity, the rendering engine is used to map the material and lighting properties, and the building layout and site elevation decision data are output. 6.The 3D terrain model based land use planning assistance method according to claim 2, wherein, The calculation process of the slope includes: The three vertex coordinates of the triangular facet inside the voxel grid are obtained, the normal vector of the triangular facet is constructed by using vector product operation, and the cosine value of the included angle between the normal vector and the vertical axis of the world coordinate system is calculated; The inclination angle value of the ground surface relative to the horizontal plane is obtained by inverse cosine function analysis, and the angle value is normalized and mapped to the preset slope classification interval, and the classification result is written into the terrain voxel dataset as the slope attribute value of the voxel. 7.The 3D terrain model based land use planning assistance method according to claim 4, wherein, The generation process of the earthwork and sunshine evaluation matrix includes: The numerical matrix of the total earthwork excavation and filling demand and the cumulative sunshine hours is obtained, the engineering cost weight factor and the residential comfort weight factor are introduced, the development suitability score of each spatial unit is calculated based on the weighted comprehensive evaluation model, and the following comprehensive evaluation matrix calculation logic is constructed: ； wherein, representative coordinates a development suitability composite score for the voxel, a measured accumulated sunshine hours representative of the voxel location, a minimum threshold of sunshine hours representative of the sunshine standard, a real elevation value representative of the voxel, a reference design elevation representative of the planning setting, a maximum elevation difference value allowed by the site, a residential comfort weight factor representative of, a construction cost weight factor representative of. 8.The 3D terrain model based land use planning assistance method of claim 4, wherein, The calculation process of the cumulative sunshine hours includes: Taking the winter solstice day as the reference date, a solar position vector sequence is constructed according to a preset time sampling interval, and a virtual ray is emitted from the geometric center of each voxel to be detected in the direction indicated by the solar position vector sequence; If there is no blocking voxel marked as non-transparent attribute on the virtual ray path, the time sampling interval corresponding to the sunshine duration is accumulated, and after all the sampling time points in a day are traversed, the effective illumination total duration of the voxel is obtained as the cumulative sunshine duration. 9.The 3D terrain model based land use planning assistance method of claim 5, wherein, The execution process of the morphological erosion and dilation operation includes: A structure element including a center voxel and 26-neighborhood voxels is defined, in the erosion stage, if all voxels in the coverage range of the structure element do not satisfy the retention condition, the state of the center voxel is flipped from buildable to unbuildable, thereby stripping the edge noise points of the building voxel group; In the dilation stage, if there is a voxel satisfying the target elevation continuity condition in the neighborhood of the center voxel, the state of the center voxel is assimilated to the state of the neighborhood voxel, thereby filling the gap inside the building layout and smoothing the mutation area of the terrain surface.

10. A 3D terrain model-based land use planning assistance system, characterized by, The system is used to implement the 3D terrain model based land use planning auxiliary method of any one of claims 1-9, and the system includes: A terrain discretization construction module is configured to obtain original point cloud data and digital elevation information, generate a three-dimensional terrain surface by using interpolation logic, calculate coordinates, elevations and slopes by using an octree segmentation algorithm, and output a terrain voxel dataset; A planning constraint mapping module is configured to extract upper limit of volume rate, building density and building height limit parameters, convert them into geometric constraint boundaries, and perform traversal detection and state marking on the terrain voxel dataset to generate a planning voxel field; A site suitability analysis module is configured to perform connectivity analysis on the planning voxel field, calculate earthwork excavation and filling demand and cumulative sunshine duration, and generate earthwork and sunshine evaluation matrices by comprehensively considering each index; A spatial form optimization module is configured to identify the sunshine insufficient and earthwork over-limit areas, correct voxel state and elevation information by using morphological erosion and dilation operation, perform visual processing on the model, and output building layout and site elevation decision data.

Images