Intelligent management and control method and system for low-altitude economy based on territorial space information model

Abstract

The application discloses a low-altitude economy intelligent management and control method and system based on a territorial space information model, and comprises the following steps: S1, acquiring city three-dimensional point clouds and territorial space planning data, preprocessing to generate a to-be-recognized point cloud set and extracting land use states and building height initial values; S2, using a dynamic interest region sampler and PointNet++ to extract features from the to-be-recognized point cloud set and output a three-dimensional semantic point set; S3, generating a three-dimensional obstacle voxel grid according to voxelization of the three-dimensional semantic point set; S4, calling CA-Markov and a continuous height function to generate a height evolution prediction sequence; S5, mapping the height evolution prediction sequence to generate a future three-dimensional obstacle voxel grid; S6, constructing a time series three-dimensional risk grid based on the current and future voxels and calculating a feasible flight path; S7, performing time availability calculation on the feasible flight path and outputting a low-altitude flight path planning result. The application realizes accurate identification of multi-period three-dimensional obstacles and continuous availability planning of a flight path, and improves the safety and timeliness of low-altitude operation.

Description

Technical Field

[0001] This invention relates to the field of intelligent planning of low-altitude air routes, and in particular to a method and system for intelligent economic management and control of low-altitude air routes based on a land spatial information model. Background Technology

[0002] With the development of the low-altitude economy, the operational demand for various drones and low-altitude vehicles in urban airspace is constantly increasing, placing higher demands on refined low-altitude route planning and dynamic control. Existing low-altitude control technologies typically rely on two-dimensional land use data or simple three-dimensional building data, which can only represent static spatial structures and cannot accurately identify the true three-dimensional morphology of various types of obstacles in the urban environment. Three-dimensional point cloud recognition technology has been applied to urban spatial analysis, but it suffers from insufficient sampling efficiency and inadequate local feature extraction in large-scale scenarios, making it difficult to meet the needs of refined analysis of high-density obstacle areas. Existing three-dimensional scene recognition results mostly focus on the current moment and lack the ability to model future spatial evolution trends, making low-altitude route planning only possible under static spatial conditions and unable to adapt to the changes in airspace brought about by urban construction and development.

[0003] Land use change prediction typically employs the traditional CA-Markov model. While this model can simulate the evolution of land use categories, height information is discretized into a limited number of levels, failing to represent the continuous variation patterns of actual building heights. Land use change results are difficult to effectively integrate with three-dimensional obstacle space, making it impossible to construct a three-dimensional risk representation with a time dimension. In route planning, existing methods mostly use static three-dimensional grids to generate feasible paths, lacking analysis of time-series risks. This leads to planning results potentially conflicting with newly formed obstacles in the future, resulting in insufficient reliability.

[0004] Therefore, how to provide a method and system for intelligent economic management and control of low-altitude areas based on a land spatial information model is a problem that urgently needs to be solved by those skilled in the art. Summary of the Invention

[0005] One objective of this invention is to propose a low-altitude economic intelligent management and control method and system based on a land spatial information model. This invention utilizes PointNet++, CA-Markov evolution prediction, and continuous altitude function modeling to achieve the collaborative processing of accurate current obstacle identification and future obstacle evolution prediction, thereby completing low-altitude route planning. It has the advantages of high identification accuracy, strong prediction capability, and high route stability.

[0006] The low-altitude economic intelligent control method based on a land spatial information model according to an embodiment of the present invention includes the following steps: S1. Acquire urban 3D point cloud data and land spatial planning data. Perform preprocessing on the urban 3D point cloud data to generate a set of points to be identified. Extract the initial values ​​of land use status and building height from the land spatial planning data. S2. Construct a PointNet++ model containing a dynamic region of interest (ROI) sampler, perform feature extraction and ROI determination on the point cloud to be identified, and generate an enhanced point cloud through the dynamic ROI sampler. S3. Perform point-level classification processing on the interest enhancement point cloud, output a three-dimensional semantic point set, and generate a three-dimensional obstacle object pixel grid based on the three-dimensional semantic point set. S4. Construct a CA-Markov model using a continuous height function based on land use status and initial building height values, set the height evolution step size, and generate a height evolution prediction sequence based on the input of land use status and initial building height values ​​according to the CA-Markov model using a continuous height function. S5. Map the height evolution prediction sequence to the corresponding spatial location of the 3D obstacle object pixel grid to generate the future 3D obstacle object pixel grid. S6. Construct a time-series three-dimensional risk grid based on the three-dimensional obstacle object pixel grid and the future three-dimensional obstacle object pixel grid, and perform low-altitude route feasibility calculation based on the time-series three-dimensional risk grid to generate a candidate route set. S7. Perform time-dimensional stability calculations on the candidate route set and output the route with the highest stability as the low-altitude route planning result.

[0007] Optionally, the dynamic region of interest sampler in S2 specifically includes: S21. Calculate the local neighborhood point set for each point in the point cloud to be identified, record the number of points as the local point density, record the absolute value of the height difference of the points as the local height difference, and record the spatial curvature as the local spatial curvature. S22. Calculate the global average point density of all points in the point cloud to be identified, and initialize the interest weight values ​​of all points to zero. S23. Set local point density threshold, height difference threshold, spatial curvature threshold, and local point density maximum value determination condition for all points in the point cloud to be identified, and set interest weight increment to form an interest weight matrix. S24. Perform sampling based on the interest weight matrix, mark points with interest weight values ​​greater than the preset sampling threshold as interest retention points, mark points with interest weight values ​​less than or equal to the preset sampling threshold as interest rejection points, and form a primary interest point set from all interest retention points. S25. Perform neighborhood uniformity calculation on the primary interest point set. Calculate the uniformity adjustment coefficient for each point in the primary interest point set based on the spatial distribution differences of the local neighborhood point set. Multiply the uniformity adjustment coefficient by the interest weight value to obtain the updated interest weight value. S26. Perform final sampling based on the updated interest weight values. Mark the points whose updated interest weight values ​​are greater than the preset final sampling threshold as the final retained points, and mark the points whose updated interest weight values ​​are less than or equal to the preset final sampling threshold as the final discard points. Combine all the final retained points into an interest enhancement point cloud.

[0008] Optionally, the PointNet++ model specifically includes: Input the point cloud to be identified and the point cloud of interest enhancement into the input layer of the PointNet++ model. Perform set difference operation on the point cloud to be identified, and form the basic point cloud with the points that do not belong to the point cloud of interest enhancement. The first hierarchical sampling and the first grouping operation are called on the interest enhancement point cloud. The first hierarchical sampling sets the sampling radius and the number of sampling points. The first grouping operation divides the interest enhancement point cloud into several first local point sets. Several scale feature extraction structures are called on each first local point set, and several first sampling scales are set to generate a first local feature set. The second hierarchical sampling and second grouping operations are called on the basic point cloud. The second hierarchical sampling sets a sampling radius greater than the first sampling radius and sets a number of sampling points less than the number of the first sampling points. The second grouping operation divides the basic point cloud into several second local point sets. Several scale feature extraction structures are called on each second local point set, and several second sampling scales are set to generate a second local feature set. The first local feature set and the second local feature set are input into the feature upsampling structure. The feature propagation method of spatial neighborhood interpolation is called to interpolate the two types of local features to all points in the point cloud set to be identified according to the spatial location of the point cloud to form a point-by-point feature set. Call the point-by-point classification operation on the point-by-point feature set to output a three-dimensional semantic point set.

[0009] Optionally, S3 specifically includes: Spatial voxelization is performed on the 3D semantic point set, the voxel edge length is set, and the 3D space is divided into regular voxel grids according to the voxel edge length; each point in the 3D semantic point set is mapped to the corresponding voxel position according to spatial coordinates, and points belonging to the building, structure, tower, and power line categories are marked as obstacle points; voxels containing obstacle points are marked as obstacle voxels, and voxels not containing obstacle points are marked as non-obstacle voxels; all obstacle voxels are arranged in spatial index order to form a 3D obstacle voxel grid.

[0010] Optionally, S4 specifically includes: The initial building height value is written into the height array point by point according to the spatial location index. The height values ​​of adjacent time and adjacent spatial location in the height array are called to generate a continuous sequence of initial height values. The CA-Markov model is invoked to perform evolutionary calculations on land use status, and the land use status transition probabilities in the evolution results are written into a probability array according to spatial location. For each spatial location, the continuous height function is called, taking the initial value of continuous height, the land use state transition probability of the corresponding location in the probability array, the continuous height variable within the neighborhood obtained by neighborhood search, and the set height evolution step size as the function input to calculate the height change. The change in altitude is added to the value of the continuous altitude variable to form the continuous altitude variable for the next period; The continuous altitude variables for all periods are arranged in spatial index order to form an altitude evolution prediction sequence.

[0011] Optionally, the CA-Markov model specifically includes: The land use status is written into the land use array point by point according to the spatial location index, and the land use status of different periods is used as the input of several periods of land use. Perform state transition probability calculation on the land use array, and generate a land use state transition probability matrix according to the transition frequency between different land use categories; Invoke the neighborhood search operation to retrieve the land use status within a preset neighborhood for each spatial location, forming a neighborhood land use set. The corresponding transition probabilities in the state transition probability matrix and the land use states in the neighboring land use set are numerically calculated according to the set spatial weights to generate the land use state for the next period.

[0012] Optionally, S5 specifically includes: S51. Call the height comparison operation on each voxel in the 3D obstacle object voxel grid to compare the continuous height variables of the corresponding spatial position in the height prediction array with the voxel space height range. S52. When the continuous height variable is within the height range of the voxel space, mark the voxel as a future obstacle voxel; when the continuous height variable is not within the height range of the voxel space, keep the voxel as a non-obstacle voxel. S53. Arrange all future obstacle voxels in spatial index order to form a future three-dimensional obstacle voxel grid.

[0013] Optionally, S7 specifically includes: S71. Mark the obstacle voxel risk value in the three-dimensional obstacle voxel grid and the future three-dimensional obstacle voxel grid as 1, and the non-obstacle voxel risk value as 0, and merge them into a time series three-dimensional risk grid. S72. Write the voxels with a risk value of 0 into the set of passable voxels; call the spatial neighborhood expansion operation in the set of passable voxels to expand the spatial connectivity path voxel by voxel from the starting voxel to the ending voxel, and record the voxel sequence that satisfies the continuity of time series and spatial connectivity as a feasible route. S73. Arrange all feasible routes in the order of generation to form a candidate route set.

[0014] Optionally, the stability calculation in S7 specifically includes: For each voxel in the route, the time series risk statistics operation is called to read all time slice risk values ​​of the current voxel in the time series three-dimensional risk grid. The sum of the time slice risk values ​​is recorded as the stability score, and the one with the lowest stability score is the best. If the stability is the same, the shortest path is output.

[0015] The low-altitude economic intelligent control system based on a land spatial information model according to an embodiment of the present invention includes the following modules: The data acquisition module is used to acquire urban 3D point cloud data and land spatial planning data, and to perform preprocessing on the urban 3D point cloud data to generate a set of points to be identified. The point cloud recognition module is used to perform dynamic region of interest sampling and PointNet++ recognition operations on the point cloud set to be recognized, output a 3D semantic point set and generate a 3D obstacle object pixel grid; The height evolution module is used to construct a continuous initial height sequence based on land use status and initial building height values, and to generate a height evolution prediction sequence by calling the CA-Markov model and the continuous height function. The future obstacle generation module is used to map the height evolution prediction sequence to the corresponding spatial location of the 3D obstacle object pixel grid to generate the future 3D obstacle object pixel grid. The risk construction module is used to construct a time-series 3D risk raster based on a 3D obstacle object pixel raster and a future 3D obstacle object pixel raster. The route generation and selection module is used to perform low-altitude route feasibility calculations based on a time-series three-dimensional risk grid, generate a set of candidate routes, and output low-altitude route planning results based on stability.

[0016] The beneficial effects of this invention are: (1) By introducing a dynamic region of interest sampler and the PointNet++ model, this invention can perform structured recognition of large-scale urban 3D point clouds, effectively improving the extraction accuracy of high-density areas and complex obstacles, and realizing a refined 3D representation of various types of obstacles such as buildings, towers, and power lines, thus solving the problem of insufficient obstacle recognition accuracy in the prior art.

[0017] (2) This invention combines the CA-Markov model with a continuous height function to construct a height evolution prediction mechanism that can change over time, thereby realizing a continuous expression of future building height changes and urban spatial morphology changes. This solves the problem that traditional discrete height models cannot reflect real height changes, and allows future obstacles to be accurately predicted.

[0018] (3) This invention constructs a time-series three-dimensional risk grid based on current and future three-dimensional obstacle object pixel grids, and performs time-series feasibility calculation and stability assessment on candidate routes, thereby realizing the improvement from static route planning to dynamic route control. The generated routes have high stability and low conflict risk, improving the safety and reliability of low-altitude economic operation. Attached Figure Description

[0019] The accompanying drawings are provided to further illustrate the invention and form part of the specification. They are used in conjunction with embodiments of the invention to explain the invention and do not constitute a limitation thereof. In the drawings: Figure 1 This is a flowchart of the low-altitude economic intelligent management and control method based on the national land spatial information model proposed in this invention; Figure 2 This is a structural diagram of the low-altitude economic intelligent control system based on the national land spatial information model proposed in this invention; Figure 3 This is a flowchart illustrating the future obstacle generation process of the low-altitude economic intelligent management and control method based on the national land spatial information model proposed in this invention. Detailed Implementation

[0020] The present invention will now be described in further detail with reference to the accompanying drawings. These drawings are simplified schematic diagrams, illustrating only the basic structure of the invention, and therefore only show the components relevant to the invention.

[0021] refer to Figures 1-3 The low-altitude economic intelligent management and control method based on the national land spatial information model includes the following steps: S1. Acquire urban 3D point cloud data and land spatial planning data. Perform preprocessing on the urban 3D point cloud data to generate a set of points to be identified. Extract the initial values ​​of land use status and building height from the land spatial planning data. S2. Construct a PointNet++ model containing a dynamic region of interest (ROI) sampler, perform feature extraction and ROI determination on the point cloud to be identified, and generate an enhanced point cloud through the dynamic ROI sampler. S3. Perform point-level classification processing on the interest enhancement point cloud, output a three-dimensional semantic point set, and generate a three-dimensional obstacle object pixel grid based on the three-dimensional semantic point set. S4. Construct a CA-Markov model using a continuous height function based on land use status and initial building height values, set the height evolution step size, and generate a height evolution prediction sequence based on the input of land use status and initial building height values ​​according to the CA-Markov model using a continuous height function. S5. Map the height evolution prediction sequence to the corresponding spatial location of the 3D obstacle object pixel grid to generate the future 3D obstacle object pixel grid. S6. Construct a time-series three-dimensional risk grid based on the three-dimensional obstacle object pixel grid and the future three-dimensional obstacle object pixel grid, and perform low-altitude route feasibility calculation based on the time-series three-dimensional risk grid to generate a candidate route set. S7. Perform time-dimensional stability calculations on the candidate route set and output the route with the best stability as the low-altitude route planning result.

[0022] In this embodiment, S1 specifically includes: S11. Acquire urban 3D point cloud data and land spatial planning data; urban 3D point cloud data is obtained from airborne lidar scanning, oblique photogrammetry, and ground 3D scanning; land spatial planning data is obtained from the land spatial information model. When preprocessing urban 3D point cloud data, the point cloud data is denoised to eliminate isolated and outlier points; coordinate system transformation is performed on the point cloud coordinates to unify them to the urban spatial reference coordinate system; sampling processing is performed on the point cloud density; normalization processing is performed on the point cloud height; after preprocessing, a set of points to be identified is formed. When extracting land use status from territorial spatial planning data, land use category information is read according to the zoning code of the planning map patch, and the land use category associated with each spatial grid location is recorded as the land use status; when extracting the initial value of building height, the building height control value corresponding to each spatial location is recorded as the initial value of building height according to the planned building height control line or the planned building strength index.

[0023] In this embodiment, the dynamic region of interest sampler in S2 specifically includes: S21. Calculate the local neighborhood point set for each point in the point cloud to be identified, record the number of points as the local point density, record the absolute value of the height difference of the points as the local height difference, and record the spatial curvature as the local spatial curvature. S22. Calculate the global average point density of all points in the point cloud to be identified, and initialize the interest weight values ​​of all points to zero. S23. Set local point density threshold, height difference threshold, spatial curvature threshold, and local point density maximum value determination condition for all points in the point cloud to be identified, and set interest weight increment to form an interest weight matrix. S24. Sampling is performed based on the interest weight matrix. Points with interest weight values ​​greater than a preset sampling threshold are marked as interest retention points, and points with interest weight values ​​less than or equal to the preset sampling threshold are marked as interest rejection points. All interest retention points are combined into a primary interest point set. In this embodiment, the sampling threshold is selected as 3. S25. Perform neighborhood uniformity calculation on the primary interest point set. Calculate the uniformity adjustment coefficient for each point in the primary interest point set based on the spatial distribution differences of its local neighborhood point set. Multiply the uniformity adjustment coefficient by the interest weight value to obtain the updated interest weight value. In this embodiment, the uniformity adjustment coefficient is calculated as follows: For each point in the primary interest point set, call the radius neighborhood retrieval method to retrieve neighborhood points within a fixed spatial radius. Calculate the spatial coordinate difference of the neighborhood points as the neighborhood point distance. Sum all the neighborhood point distances and divide by the number of neighborhood points to obtain the average neighborhood distance. Set the uniformity benchmark distance to 1, perform a ratio operation between the average neighborhood distance and the uniformity benchmark distance, and record the ratio as the uniformity adjustment coefficient. S26. Perform final sampling based on the updated interest weight values. Mark points whose updated interest weight values ​​are greater than a preset final sampling threshold as final retained points, and mark points whose updated interest weight values ​​are less than or equal to the preset final sampling threshold as final discard points. Combine all the final retained points into an interest enhancement point agglomeration. In this embodiment, the final sampling threshold is set to 4.

[0024] In this embodiment, S21 specifically includes: S211. Take each point as the center point, call the fixed radius neighborhood retrieval method to retrieve all points within the spatial range of the center point, and record the retrieved point set as a local neighborhood point set. S212. Record the number of points in the local neighborhood point set as the local point density; S213. Calculate the difference between the height value of each point in the local neighborhood point set and the height value of the center point, and perform an absolute value operation on the difference. Record the maximum value among all the absolute value operation results as the local height difference. S214. Perform surface fitting processing based on the local neighborhood point set, call the least squares plane fitting method to obtain the local fitting plane, calculate the angle between the normal vector of the local fitting plane and the normal vector of the neighborhood points, and record the average value of the angle calculation as the local spatial curvature.

[0025] In this embodiment, S23 specifically includes: S231. Set the local point density threshold to the global average point density, and increase the interest weight value of points whose local point density is greater than the global average point density by 2. S232. Set the height difference threshold to the median of all local height differences, and increase the interest weight value of points whose absolute local height difference is greater than the height difference threshold by 1. S233. Set the spatial curvature threshold to the average value of all local spatial curvatures, and increase the interest weight value of points whose local spatial curvature values ​​are greater than the spatial curvature threshold by 1. S234. Set the condition for determining the maximum local point density as having the maximum local point density in the local neighborhood point set, and increase the interest weight value of the points that satisfy the condition by 1.

[0026] In this embodiment, the PointNet++ model specifically includes: Input the point cloud to be identified and the point cloud of interest enhancement into the input layer of the PointNet++ model. Perform set difference operation on the point cloud to be identified, and form the basic point cloud with the points that do not belong to the point cloud of interest enhancement. The first hierarchical sampling and first grouping operations are invoked on the interest enhancement point cloud. The first hierarchical sampling sets the sampling radius and the number of sampling points. The first grouping operation divides the interest enhancement point cloud into several first local point sets. Several scale feature extraction structures are invoked on each first local point set, and several first sampling scales are set to generate a first local feature set. In this embodiment, the first sampling radius is set to 10m, the number of first sampling points is 4096, and the first sampling scale is three scales: 10m, 30m, and 60m. The second hierarchical sampling and second grouping operations are invoked on the base point cloud. The second hierarchical sampling sets a sampling radius greater than the first sampling radius and a number of sampling points less than the number of the first sampling points. The second grouping operation divides the base point cloud into several second local point sets. Several scale feature extraction structures are invoked on each second local point set, and several second sampling scales are set to generate a second local feature set. In this embodiment, the second sampling radius is 30m, the number of second sampling points is 2048, and the second sampling scales are 30m and 90m. The first local feature set and the second local feature set are input into the feature upsampling structure. The feature propagation method of spatial neighborhood interpolation is called to interpolate the two types of local features to all points in the point cloud set to be identified according to the spatial position of the point cloud to form a point-by-point feature set. In this embodiment, in the spatial neighborhood interpolation parameters, K=3, that is, the number of neighborhoods is 3. The interpolation weight adopts the inverse distance weight based on the Euclidean distance of the neighborhood points. When calculating the weight value, the distance is added to the bias amount of 0.01. A point-by-point classification operation is called on the point-by-point feature set to output a three-dimensional semantic point set. In this embodiment, the point-by-point classification operation has a structure consisting of three fully connected layers, with output dimensions of 128, 64, and the number of categories, respectively; the activation function used is ReLU.

[0027] In this embodiment, S3 specifically includes: Spatial voxelization is performed on the 3D semantic point set, and the voxel edge length is set. The 3D space is divided into a regular voxel grid according to the voxel edge length. Each point in the 3D semantic point set is mapped to the corresponding voxel position according to spatial coordinates. Points belonging to the building, structure, pole, and power line categories are marked as obstacle points. Voxels containing obstacle points are marked as obstacle voxels, and voxels not containing obstacle points are marked as non-obstacle voxels. All obstacle voxels are arranged in spatial index order to form a 3D obstacle voxel grid. In this embodiment, the voxel edge length is set to 5m.

[0028] In this embodiment, S4 specifically includes: The initial building height value is written into the height array point by point according to the spatial location index. The interpolation calculation method is called to generate a continuous sequence of initial height values ​​for adjacent time and adjacent spatial location in the height array. The CA-Markov model is invoked to perform evolutionary calculations on land use status, and the land use status transition probabilities in the evolution results are written into a probability array according to spatial location. For each spatial location, the continuous height function is called, taking the initial value of continuous height, the land use state transition probability of the corresponding location in the probability array, the continuous height variable within the neighborhood obtained by neighborhood search, and the set height evolution step size as the function input to calculate the height change. The change in altitude is added to the value of the continuous altitude variable to form the continuous altitude variable for the next period; The continuous altitude variables for all periods are arranged in spatial index order to form an altitude evolution prediction sequence.

[0029] In this embodiment, the continuous height function specifically includes: When calculating the change in altitude, the following method is used: the current altitude is used as the initial continuous altitude value; the average altitude of the neighborhood is obtained by averaging the continuous altitude values ​​of the neighborhood; the land use state transition probability is used as the land use probability; the altitude difference coefficient is set to 0.3, the land use coefficient to 5, and the altitude step size to 1; and the calculation is performed according to the following formula: Altitude change = 0.3 * (neighborhood average altitude) (Current location height) + 5·Land use probability· 1.

[0030] In this embodiment, the CA-Markov model specifically includes: The land use status is written into the land use array point by point according to the spatial location index, and the land use status of different periods is used as the input of several periods of land use. Perform state transition probability calculation on the land use array, and generate a land use state transition probability matrix according to the transition frequency between different land use categories; Invoke the neighborhood search operation to retrieve the land use status within a preset neighborhood for each spatial location, forming a neighborhood land use set. The corresponding transition probabilities in the state transition probability matrix and the land use states in the neighboring land use set are numerically calculated according to the set spatial weights to generate the land use state for the next period.

[0031] In this embodiment, the CA-Markov model specifically includes: The study area is divided into regular grids according to a preset spatial resolution; the land use category of each grid unit is written into the land use array row by row according to the spatial location index, and the land use data collected in different years are written into the land use array in sequence to form several periods of land use input data; the land use category is recorded using an integer encoding method. After reading several periods of land use arrays, the transfer statistics of land use status between any two adjacent periods are performed; the ratio of the number of transfers to the total number of samples is used as the transfer probability. For each spatial location, a neighborhood search operation is invoked, using a Moore 3×3 neighborhood structure. Taking the current grid cell as the center, the land use categories of the eight surrounding grid cells are retrieved, and the nine retrieved land use categories are written into the neighborhood land use set. If the center cell is located in the boundary area, the land use category of the nearest valid grid cell is added to the missing location. When updating land use status, a spatial weight is assigned to each location in the neighboring land use set; the spatial weight of the central location is set to 0.4, and the spatial weights of the other eight neighboring locations are all set to 0.075. The CA-Markov update rule performs numerical calculations on the transition probability of the current land use category in the state transition probability matrix and the category value in the neighboring land use set according to spatial weights to form a land use score array; the land use category with the largest value in the land use score array is taken as the land use state of the next period, and the land use states of all spatial locations are written into the land use state evolution sequence in spatial index order; After obtaining the land use status for the next period, the transition probability corresponding to the current land use category in the state transition probability matrix is ​​written into the probability array; the value of the probability array is used as the land use probability input for the continuous height function to calculate the height change; the continuous height function uses the land use probability provided by the probability array, the continuous height variable in the neighborhood range, and the height step size as input data to perform the height change calculation.

[0032] In this embodiment, S5 specifically includes: S51. Call the height comparison operation on each voxel in the 3D obstacle object voxel grid to compare the continuous height variables of the corresponding spatial position in the height prediction array with the voxel space height range. S52. When the continuous height variable is within the height range of the voxel space, mark the voxel as a future obstacle voxel; when the continuous height variable is not within the height range of the voxel space, keep the voxel as a non-obstacle voxel. S53. Arrange all future obstacle voxels in spatial index order to form a future three-dimensional obstacle voxel grid.

[0033] In this embodiment, S7 specifically includes: S71. Mark the obstacle voxel risk value in the three-dimensional obstacle voxel grid and the future three-dimensional obstacle voxel grid as 1, and the non-obstacle voxel risk value as 0, and merge them into a time series three-dimensional risk grid. S72. Write the voxels with a risk value of 0 into the set of passable voxels; call the spatial neighborhood expansion operation in the set of passable voxels to expand the spatial connectivity path voxel by voxel from the starting voxel to the ending voxel, and record the voxel sequence that satisfies the continuity of time series and spatial connectivity as a feasible route. S73. Arrange all feasible routes in the order of generation to form a candidate route set.

[0034] In this embodiment, the obstacle voxel risk values ​​of the three-dimensional obstacle voxel grid and the future three-dimensional obstacle voxel grid are set to 1, the non-obstacle voxel risk values ​​are set to 0, and the current risk value and the predicted risk value are written into two time slices in spatial index order to form a time series three-dimensional risk grid. Voxels with a risk value of 0 are written into the set of passable voxels. Voxels in the set of passable voxels undergo spatial neighborhood expansion according to the 26-neighborhood structure, expanding outward from the starting voxel. During the expansion process, the risk value of each voxel at its corresponding spatial position in the next time slice is checked. When the risk value is 0, the voxel is written into the expansion queue. When the risk value is 1, the expansion in the current direction is terminated. When the voxels in the extended queue match the endpoint voxel, the voxel index sequence from the starting voxel to the endpoint voxel is recorded as a feasible route; all voxel sequences that satisfy spatial connectivity and temporal continuity are arranged in the recording order to form a candidate route set.

[0035] In this embodiment, the stability calculation in S7 specifically includes: For each voxel in the route, the time series risk statistics operation is called to read all time slice risk values ​​of the current voxel in the time series three-dimensional risk grid. The sum of the time slice risk values ​​is recorded as the stability score, and the one with the lowest stability score is the best. If the stability is the same, the shortest path is output.

[0036] The low-altitude economic intelligent control system based on a land spatial information model according to an embodiment of the present invention includes the following modules: The data acquisition module is used to acquire urban 3D point cloud data and land spatial planning data, and to perform preprocessing on the urban 3D point cloud data to generate a set of points to be identified. The point cloud recognition module is used to perform dynamic region of interest sampling and PointNet++ recognition operations on the point cloud set to be recognized, output a 3D semantic point set and generate a 3D obstacle object pixel grid; The height evolution module is used to construct a continuous initial height sequence based on land use status and initial building height values, and to generate a height evolution prediction sequence by calling the CA-Markov model and the continuous height function. The future obstacle generation module is used to map the height evolution prediction sequence to the corresponding spatial location of the 3D obstacle object pixel grid to generate the future 3D obstacle object pixel grid. The risk construction module is used to construct a time-series 3D risk raster based on a 3D obstacle object pixel raster and a future 3D obstacle object pixel raster. The route generation and selection module is used to perform low-altitude route feasibility calculations based on a time-series three-dimensional risk grid, generate a set of candidate routes, and output low-altitude route planning results based on stability.

[0037] Figure 2This is a schematic diagram of the low-altitude economic intelligent control system based on a territorial spatial information model proposed in this invention. Arrows pointing from the data acquisition module to the point cloud recognition module indicate that preprocessed 3D point cloud data and territorial spatial planning data are input to the point cloud recognition module. Arrows pointing from the point cloud recognition module to the altitude evolution module and the future obstacle generation module indicate that 3D obstacle pixel grids and related spatial information serve as the basic input for subsequent spatiotemporal prediction calculations. Arrows pointing from the altitude evolution module to the future obstacle generation module indicate that the altitude evolution prediction sequence is used to generate future 3D obstacle pixel grids. Arrows pointing from the future obstacle generation module to the risk construction module indicate that current and future obstacle information jointly participate in the construction of a time-series 3D risk grid. Arrows pointing from the risk construction module to the route generation and selection module indicate that risk grid data is used for route feasibility calculations and stability assessments, thereby completing low-altitude route planning.

[0038] Figure 3 This is a schematic diagram of the future obstacle generation process proposed in this invention. Arrows pointing from the height evolution prediction sequence generation step to the height comparison step indicate that the predicted height data is input into the voxel height range determination process. Arrows pointing from the height comparison step to the future obstacle voxel labeling processing and non-obstacle preservation processing indicate that different processing paths are executed based on the height determination results. Arrows pointing from each determination path to the future 3D obstacle voxel raster generation step indicate that the determination results for each spatial location are finally aggregated to form the future 3D obstacle voxel raster. The arrow connections clearly represent the data flow direction and logical processing relationships in the future obstacle generation process.

[0039] Example 1: To verify the feasibility of this invention in practice, it was applied to a typical urban area as a verification region. A verification set was constructed by collecting urban 3D point cloud data and land use planning data from consecutive years to test the continuous availability of the low-altitude flight path planning method of this invention over time. This urban area contains densely built-up areas, as well as power corridors, tree cover, and low-rise building clusters, with significant variations in altitude and land use between different years. Traditional static 3D obstacle modeling methods struggle to maintain flight path effectiveness in such dynamic environments, frequently exhibiting problems such as path breaks, reduced passable voxels, and flight path failures in certain years. This invention enhances the sampling capability of key point cloud components through a dynamic region of interest sampler, improves 3D semantic recognition accuracy through PointNet++, and constructs altitude evolution trends using a CA-Markov model and a continuous altitude function, generating more usable low-altitude flight paths under historical multi-period changes.

[0040] In the validation set, ten typical regions were selected as samples. Preprocessing and interest-enhanced sampling were performed on the point cloud for each period, so that high-density, high-curvature, and high-height-difference regions received higher weight sampling, making the local features extracted by PointNet++ more stable. The identified categories such as buildings, structures, and power lines were converted into voxelized obstacle representations to form a 3D obstacle voxel grid for each period.

[0041] To construct a continuous trend of altitude variation over time, this invention inputs multi-period altitude data into a continuous altitude function and combines it with CA-Markov land use state transition probabilities for altitude evolution fitting, making the altitude trend more consistent with the continuity of annual changes. For example, in the validation data of region A04, the actual average altitude variation over the years is 6.0 meters, and the average altitude prediction error of this invention is 1.4 meters, while the error of the traditional method increases to 3.8 meters due to the use of discrete altitude levels. In region A06, the altitude variation is even greater, with the error of the traditional method reaching as high as 4.2 meters, while the error of this invention is only 1.5 meters, indicating that this invention has a stronger ability to fit the altitude evolution trend.

[0042] In the time-series 3D risk construction step, obstacle voxel grids from multiple periods are aligned to determine whether a route remains within a passable voxel in each period. The "time availability score" is further defined as the number of time slices in which a path remains passable across all time slices.

[0043] Comparative analysis across ten regions shows that the time availability of this invention is superior to that of the traditional model in all regions. For example, in region A05, the traditional method scores 18 in time availability, while this invention achieves 31, indicating that the routes generated by this invention remain navigable in the vast majority of years. In regions with more drastic dynamic changes, such as A06, the traditional method only has 9 available time slots, while this invention still guarantees 23 available time slots, demonstrating better robustness and cross-time period adaptability. In regions with moderate changes, such as A03 and A08, this invention also shows significant improvement, achieving time availability scores of 28 and 29 respectively, while the traditional method only achieves 14 and 15.

[0044] Overall, in areas with significant altitude variations, this invention reduces obstacle misjudgments through more refined altitude fitting, maintaining stable path continuity over time. In areas with less or moderate altitude variations, the paths generated by this invention remain passable for more years, demonstrating higher temporal availability. Combining results on altitude prediction error, obstacle fitting accuracy, and the number of consecutive passable time slices, this invention outperforms traditional methods in historical multi-period data validation. Experiments show that the temporal availability score of this invention is generally more than twice that of traditional methods, validating its multi-period robustness and low-altitude route continuity capabilities, demonstrating its application value in dynamically evolving urban environments.

[0045] Table 1: Comparison of time availability between the present invention and traditional methods under multiple validation sets

[0046] As shown in Table 1, the present invention achieves higher time availability scores in all regions, indicating that the generated routes remain continuously passable across more historical time slices. For example, in region A05, the present invention achieves a time availability score of 31, while the traditional method only achieves 18. This demonstrates that the traditional method experiences path failures due to obstacle prediction biases in multiple years, while the present invention can stably maintain path passability. In regions with varying altitudes, such as A06, the traditional method has an altitude prediction error of 4.2 meters, resulting in a time availability score of only 9, while the present invention has an error of only 1.5 meters, improving the time availability score to 23, showcasing the advantages of continuous altitude modeling and temporal risk analysis. Overall, the present invention improves the average time availability by more than 50% across ten regions, validating the stability and reliability of the routes in dynamically changing urban environments.

[0047] The above description is only a preferred embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any equivalent substitutions or modifications made by those skilled in the art within the scope of the technology disclosed in the present invention, based on the technical solution and inventive concept of the present invention, should be covered within the scope of protection of the present invention.

Claims

1. A method for intelligent economic management and control of low-altitude areas based on a land spatial information model, characterized in that, The steps include the following: S1. Acquire urban 3D point cloud data and land spatial planning data. Perform preprocessing on the urban 3D point cloud data to generate a set of points to be identified. Extract the initial values ​​of land use status and building height from the land spatial planning data. S2. Construct a PointNet++ model containing a dynamic region of interest (ROI) sampler, perform feature extraction and ROI determination on the point cloud to be identified, and generate an enhanced point cloud through the dynamic ROI sampler. S3. Perform point-level classification processing on the interest enhancement point cloud, output a three-dimensional semantic point set, and generate a three-dimensional obstacle object pixel grid based on the three-dimensional semantic point set. S4. Construct a CA-Markov model using a continuous height function based on land use status and initial building height values, set the height evolution step size, and generate a height evolution prediction sequence based on the input of land use status and initial building height values ​​according to the CA-Markov model using a continuous height function. S5. Map the height evolution prediction sequence to the corresponding spatial location of the 3D obstacle object pixel grid to generate the future 3D obstacle object pixel grid. S6. Construct a time-series three-dimensional risk grid based on the three-dimensional obstacle object pixel grid and the future three-dimensional obstacle object pixel grid, and perform low-altitude route feasibility calculation based on the time-series three-dimensional risk grid to generate a candidate route set. S7. Perform time-dimensional stability calculations on the candidate route set and output the route with the best stability as the low-altitude route planning result.

2. The low-altitude economic intelligent control method based on a land spatial information model according to claim 1, characterized in that, The dynamic region of interest sampler in S2 specifically includes: S21. Calculate the local neighborhood point set for each point in the point cloud to be identified, record the number of points as the local point density, record the absolute value of the height difference of the points as the local height difference, and record the spatial curvature as the local spatial curvature. S22. Calculate the global average point density of all points in the point cloud to be identified, and initialize the interest weight values ​​of all points to zero. S23. Set local point density threshold, height difference threshold, spatial curvature threshold, and local point density maximum value determination condition for all points in the point cloud to be identified, and set interest weight increment to form an interest weight matrix. S24. Perform sampling based on the interest weight matrix, mark points with interest weight values ​​greater than the preset sampling threshold as interest retention points, mark points with interest weight values ​​less than or equal to the preset sampling threshold as interest rejection points, and form a primary interest point set from all interest retention points. S25. Perform neighborhood uniformity calculation on the primary interest point set. Calculate the uniformity adjustment coefficient for each point in the primary interest point set based on the spatial distribution differences of the local neighborhood point set. Multiply the uniformity adjustment coefficient by the interest weight value to obtain the updated interest weight value. S26. Perform final sampling based on the updated interest weight values. Mark the points whose updated interest weight values ​​are greater than the preset final sampling threshold as the final retained points, and mark the points whose updated interest weight values ​​are less than or equal to the preset final sampling threshold as the final discard points. Combine all the final retained points into an interest enhancement point cloud.

3. The low-altitude economic intelligent control method based on a land spatial information model according to claim 2, characterized in that, The PointNet++ model specifically includes: Input the point cloud to be identified and the point cloud of interest enhancement into the input layer of the PointNet++ model. Perform set difference operation on the point cloud to be identified, and form the basic point cloud with the points that do not belong to the point cloud of interest enhancement. The first hierarchical sampling and the first grouping operation are called on the interest enhancement point cloud. The first hierarchical sampling sets the sampling radius and the number of sampling points. The first grouping operation divides the interest enhancement point cloud into several first local point sets. Several scale feature extraction structures are called on each first local point set, and several first sampling scales are set to generate a first local feature set. The second hierarchical sampling and second grouping operations are called on the basic point cloud. The second hierarchical sampling sets a sampling radius greater than the first sampling radius and sets a number of sampling points less than the number of the first sampling points. The second grouping operation divides the basic point cloud into several second local point sets. Several scale feature extraction structures are called on each second local point set, and several second sampling scales are set to generate a second local feature set. The first local feature set and the second local feature set are input into the feature upsampling structure. The feature propagation method of spatial neighborhood interpolation is called to interpolate the two types of local features to all points in the point cloud set to be identified according to the spatial location of the point cloud to form a point-by-point feature set. Call the point-by-point classification operation on the point-by-point feature set to output a three-dimensional semantic point set.

4. The low-altitude economic intelligent control method based on a territorial spatial information model according to claim 3, characterized in that, S3 specifically includes: Spatial voxelization is performed on the 3D semantic point set, the voxel edge length is set, and the 3D space is divided into regular voxel grids according to the voxel edge length; each point in the 3D semantic point set is mapped to the corresponding voxel position according to spatial coordinates, and points belonging to the building, structure, tower, and power line categories are marked as obstacle points; voxels containing obstacle points are marked as obstacle voxels, and voxels not containing obstacle points are marked as non-obstacle voxels; all obstacle voxels are arranged in spatial index order to form a 3D obstacle voxel grid.

5. The low-altitude economic intelligent control method based on the land spatial information model according to claim 4, characterized in that, S4 specifically includes: The initial building height value is written into the height array point by point according to the spatial location index, and a continuous sequence of initial height values ​​is generated for the height values ​​of adjacent time and adjacent spatial location in the height array. The CA-Markov model is invoked to perform evolutionary calculations on land use status, and the land use status transition probabilities in the evolution results are written into a probability array according to spatial location. For each spatial location, the continuous height function is called, taking the initial value of continuous height, the land use state transition probability of the corresponding location in the probability array, the continuous height variable within the neighborhood obtained by neighborhood search, and the set height evolution step size as the function input to calculate the height change. The change in altitude is added to the value of the continuous altitude variable to form the continuous altitude variable for the next period; The continuous altitude variables for all periods are arranged in spatial index order to form an altitude evolution prediction sequence.

6. The low-altitude economic intelligent control method based on a land spatial information model according to claim 5, characterized in that, The CA-Markov model specifically includes: The land use status is written into the land use array point by point according to the spatial location index, and the land use status of different periods is used as the input of several periods of land use. Perform state transition probability calculation on the land use array, and generate a land use state transition probability matrix according to the transition frequency between different land use categories; Invoke the neighborhood search operation to retrieve the land use status within a preset neighborhood for each spatial location, forming a neighborhood land use set. The corresponding transition probabilities in the state transition probability matrix and the land use states in the neighboring land use set are numerically calculated according to the set spatial weights to generate the land use state for the next period.

7. The low-altitude economic intelligent control method based on a land spatial information model according to claim 6, characterized in that, S5 specifically includes: S51. Call the height comparison operation on each voxel in the 3D obstacle object voxel grid to compare the continuous height variables of the corresponding spatial position in the height prediction array with the voxel space height range. S52. When the continuous height variable is within the height range of the voxel space, mark the voxel as a future obstacle voxel; when the continuous height variable is not within the height range of the voxel space, keep the voxel as a non-obstacle voxel. S53. Arrange all future obstacle voxels in spatial index order to form a future three-dimensional obstacle voxel grid.

8. The low-altitude economic intelligent control method based on a territorial spatial information model according to claim 7, characterized in that, Specifically, S7 includes: S71. Mark the obstacle voxel risk value in the three-dimensional obstacle voxel grid and the future three-dimensional obstacle voxel grid as 1, and the non-obstacle voxel risk value as 0, and merge them into a time series three-dimensional risk grid. S72. Write the voxels with a risk value of 0 into the set of passable voxels; call the spatial neighborhood expansion operation in the set of passable voxels to expand the spatial connectivity path voxel by voxel from the starting voxel to the ending voxel, and record the voxel sequence that satisfies the continuity of time series and spatial connectivity as a feasible route. S73. Arrange all feasible routes in the order of generation to form a candidate route set.

9. The low-altitude economic intelligent control method based on a land spatial information model according to claim 8, characterized in that, The stability calculation in S7 specifically includes: For each voxel in the route, the time series risk statistics operation is called to read all time slice risk values ​​of the current voxel in the time series three-dimensional risk grid. The sum of the time slice risk values ​​is recorded as the stability score, and the one with the lowest stability score is the best. If the stability is the same, the shortest path is output.

10. A low-altitude economic intelligent control system based on a land spatial information model, applied to the low-altitude economic intelligent control method based on a land spatial information model as described in any one of claims 1 to 9, characterized in that, Includes the following modules: The data acquisition module is used to acquire urban 3D point cloud data and land spatial planning data, and to perform preprocessing on the urban 3D point cloud data to generate a set of points to be identified. The point cloud recognition module is used to perform dynamic region of interest sampling and PointNet++ recognition operations on the point cloud set to be recognized, output a 3D semantic point set and generate a 3D obstacle object pixel grid; The height evolution module is used to construct a continuous initial height sequence based on land use status and initial building height values, and to generate a height evolution prediction sequence by calling the CA-Markov model and the continuous height function. The future obstacle generation module is used to map the height evolution prediction sequence to the corresponding spatial location of the 3D obstacle object pixel grid to generate the future 3D obstacle object pixel grid. The risk construction module is used to construct a time-series 3D risk raster based on a 3D obstacle object pixel raster and a future 3D obstacle object pixel raster. The route generation and selection module is used to perform low-altitude route feasibility calculations based on a time-series three-dimensional risk grid, generate a set of candidate routes, and output low-altitude route planning results based on stability.

Images