Terrain model construction-oriented airborne LiDAR point cloud rapid thinning method

Abstract

The invention discloses a terrain model construction-oriented airborne LiDAR point cloud rapid extraction method, and relates to the technical field of airborne LiDAR point cloud data processing.The extraction method comprises the specific steps of S100, multi-scale grid point cloud preliminary screening: setting a small grid size and a large grid size. According to the method, the original ground point cloud is subjected to double gridding thinning, the data volume of the point cloud is effectively reduced, rich topographic feature information is provided for subsequent global topographic model construction, and compared with a traditional single gridding thinning method, the point cloud data of a region with severe topographic relief can be more accurately reserved by the method, and the method has the advantages of being high in robustness, high in robustness and the like. Meanwhile, point cloud boundary points and severe topographic relief areas are extracted through large-grid-size point cloud extraction and serve as seed points for subsequent construction of the global topographic model, and the construction efficiency and precision of the topographic model are further improved.

Description

Technical Field

[0001] The present invention relates to the technical field of airborne LiDAR point cloud data processing, and in particular to a fast thinning method of airborne LiDAR point cloud for terrain model construction. Background Art

[0002] With the continuous development of remote sensing technology, airborne LiDAR (Light Detection and Ranging) technology has become an important technical means in the fields of terrain surveying and mapping, urban planning, and forestry resource investigation. LiDAR technology can accurately measure the distance, position, and three-dimensional shape of the target by emitting lasers to the target and receiving the reflected signals, thereby generating high-precision point cloud data. These point cloud data contain rich terrain information and are essential for building an accurate terrain model. In the process of terrain model construction, the processing of point cloud data is a key link, which directly affects the accuracy and efficiency of the final terrain model. At present, airborne LiDAR point cloud processing technology for terrain model construction has made significant progress, including point cloud filtering, classification, registration, and thinning. Among them, point cloud thinning is an important means to reduce the amount of data and improve processing efficiency. On the premise of maintaining terrain characteristics, it reduces the number of point clouds and reduces the computational complexity of subsequent terrain modeling.

[0003] Although the existing airborne LiDAR point cloud thinning technology has improved the data processing efficiency to a certain extent, there are still some shortcomings. Traditional point cloud thinning methods often adopt simple grid thinning or random thinning strategies. Although these methods can reduce the number of point clouds, they are prone to lose terrain feature information, resulting in a decrease in the accuracy of the constructed terrain model. In addition, some complex thinning methods, such as thinning methods based on clustering analysis, although they can better retain terrain features, have high computational complexity and are difficult to meet the real-time processing requirements of large-scale point cloud data. Therefore, how to achieve rapid thinning of point cloud data while maintaining terrain features has become a technical problem that needs to be urgently solved in the current field of terrain model construction. The shortcomings of traditional technologies are mainly reflected in the contradiction between terrain feature retention and data processing efficiency. How to find a balance between the two is the focus of current research. Summary of the invention

[0004] The purpose of the present invention is to make up for the shortcomings of the prior art and provide a method for rapid thinning of airborne LiDAR point clouds for terrain model construction. This method achieves rapid thinning of point cloud data and effective retention of terrain features through the steps of preliminary thinning of point clouds based on multi-scale grids, global terrain model construction, local encryption and redundant point removal. Compared with traditional methods, the present invention not only significantly improves the data processing efficiency, but also ensures the accuracy and detail expression of the constructed terrain model.

[0005] To solve the above technical problems, the present invention provides the following technical solutions: An airborne LiDAR point cloud rapid thinning method for terrain model construction, and the specific steps of the thinning method are as follows:

[0006] S100, multi-scale grid point cloud preliminary screening: Set the small grid size and the large grid size, and perform grid thinning on the original ground point cloud respectively to obtain the thinned point cloud with the small grid size and the thinned point cloud with the large grid size;

[0007] S200, global terrain model construction: Extract the point cloud boundary points and the areas with severe terrain undulations based on the thinned point cloud with the large grid size. According to the mapping relationship between the large grid size and the small grid size, extract the small grid points corresponding to the severe areas from the index matrix as terrain feature points. Use the extracted point cloud boundary points and terrain feature points as seed points, and construct a Delaunay triangulation through the point-by-point insertion method to form a global terrain model;

[0008] S300, recursive encryption and iterative optimization: Search for deviation points based on recursive triangle division for each triangle in the global terrain model, insert the searched deviation points and their corresponding fixed points into the global terrain model for encryption. For the new triangles where the deviation points are located in the encrypted global terrain model, re-search for deviation points based on recursive triangle division and encrypt them, and iterate until no new deviation points are generated;

[0009] S400, normal vector redundant point removal: Traverse each vertex in the triangulation, calculate the normal vector of the triangle where it is located, set the maximum normal vector angle threshold between triangles, and delete the redundant points smaller than the maximum normal vector angle threshold from the triangulation. Use the remaining triangulation vertices to form the thinned ground point cloud.

[0010] Furthermore, the specific steps of the grid thinning in the S100, multi-scale grid point cloud preliminary screening are as follows:

[0011] S101, set the small grid size as , and the large grid size as , obtain the minimum and maximum values of the point cloud data in the , directions, denoted as , , , ;

[0012] S102, calculate the number of grid rows and the number of columns corresponding to the small grid and the large grid respectively. The calculation formula is: , where is the currently used grid size, Denotes the ceiling operation;

[0013] S103, Initialize the grid storage structure to create a two-dimensional matrix , the size of which is , is the number of grid rows, is the number of grid columns. Set the value of all elements in the array to 1, indicating that the grid does not currently retain any points;

[0014] S104, Traverse the original point cloud and perform the following processing on each point in the original ground point cloud sequentially;

[0015] Calculate the grid row and column numbers where the current point is located , the formula is: , where denotes the floor operation, and check the value of grid ; grid is a two-dimensional matrix;

[0016] If grid =1, store the point number of the current point into grid ;

[0017] If grid ≠1, obtain the elevation of the currently retained point in this grid , and compare it with the elevation of the current point ; if , then replace the original value in grid with the point number of the current point ; if

[0018] S105, Traverse the grid storage structure grid, extract the points corresponding to the non-1 point numbers stored in it from the original point cloud, and form a new point cloud data, which is the downsampled point cloud under the current grid size;

[0019] S106, Replace the grid sizes in steps S102 - S105 with the small grid size and the large grid size respectively, and repeat steps S102 - S105 to obtain the downsampled point cloud with the small grid size and the downsampled point cloud with the large grid size. ​​​​​​

[0020] Furthermore, the specific steps for extracting boundary points using the Alphashapes algorithm in the global terrain model construction of S200 are as follows:

[0021] S201, thinning the point cloud with a large grid size From three-dimensional space Projected onto the horizontal plane to obtain a two-dimensional point set , and set the radius parameter of Alphashapes according to the large grid size , with a value of ; ;

[0022] S202, perform Delaunay triangulation on the two-dimensional point set to generate a Delaunay triangulation network DT, and perform the following screening on the Alpha edges and Alpha triangles;

[0023] For each edge in the Delaunay triangulation network , calculate its corresponding circumradius , if , then this edge is an Alpha edge and is retained; otherwise it is excluded;

[0024] A triangle composed of three Alpha edges is an Alpha triangle and is retained; otherwise it is excluded;

[0025] S203, among the retained Alpha triangles, count all the edges that belong to only a single Alpha triangle, which are the boundary edges, and use the endpoints of all the boundary edges to form a boundary point set is the endpoint of a certain boundary edge ;

[0026] S204, calculate the minimum bounding rectangle of the extracted boundary point set B, and verify whether it completely contains the target survey area. If the survey area is not covered, increase the value of α, and repeat steps S202 - S203 until the inclusiveness requirement is met.

[0027] Furthermore, the specific steps for determining the area with severe terrain undulation in the global terrain model construction of S200 are as follows:

[0028] (1) For each point in the point cloud thinned with a large grid size , select its 25 nearest points to form a neighborhood , estimate the normal vector of the point based on principal component analysis, and perform regularization processing on the calculated normal vector ;

[0029] (2) For each point and its neighborhood within , , calculate the included angle between its normal vector and , and the formula is: , since the normal vector has been normalized, , so , use the inverse distance weighted method to calculate the weighted average of the included angles of the normal vectors within the neighborhood of point , and set the plane distance between point and point to be , and the calculation formula is: ; ;

[0030] (3) According to the actual terrain conditions, set the threshold = of the weighted average of the normal vector included angles, and traverse each point in the downsampled point cloud with a large grid size . If the weighted average of its normal vector included angles is greater than the threshold , it is considered that the grid area where this point is located is a region with severe terrain undulation.

[0031] Further, in the S200, when constructing the global terrain model, the normal vector of point is estimated based on principal component analysis. Principal component analysis calculates the covariance matrix of the point set within the neighborhood, performs eigenvalue decomposition on the covariance matrix to obtain three eigenvalues and , as well as the corresponding eigenvectors , where, and , the eigenvector corresponding to the smallest eigenvalue is the normal vector of point .

[0032] Further, in the S200, when determining the three-dimensional points in the downsampled point cloud with a small grid size within the region with severe terrain undulation, let the row number of the grid in the region with severe terrain undulation in the large grid be and the column number be , the large grid size be , the small grid size be , and the corresponding small grid row number range and column number range The calculation formula is: , where represents rounding down, represents rounding up.

[0033] Furthermore, the specific steps of the recursive triangular division deviation point search in the recursive encryption and iterative optimization of the S300 are as follows:

[0034] S311. For the constructed Delaunay triangulation network, obtain the vertex coordinates of each triangle. Let the three vertices of the triangle be , , . Set the recursive termination conditions, the triangle area threshold and the maximum number of recursions ;

[0035] S312. Use the vector cross product formula to calculate the triangle area. The formula is: , where , ; Fit the triangle plane using the least squares method. The formula is: , where , , are coefficients;

[0036] S313. Let the centroid coordinates of the triangle be . The calculation formula is: , , . Substitute the centroid coordinates into the plane equation to obtain the model elevation . Obtain the actual point cloud elevation at the centroid from the point cloud data, and calculate the elevation deviation . The formula is: ;

[0037] S314. If the triangle area or the number of recursions reaches , stop the recursive division of this triangle; if , divide the triangle into three sub - triangles evenly from the triangle centroid to obtain new vertices and triangles. For the newly generated sub - triangles, repeat steps S312 - S314 until the recursive conditions are no longer met;

[0038] S315. Record all the points with elevation deviation values greater than , which are the deviation points.

[0039] Furthermore, the encryption of the global terrain model in the recursive encryption and iterative optimization of the S300:

[0040] S321. Let , create a new empty set of terrain encryption points;

[0041] S322, for the th triangle in the triangle set , respectively obtain the elevation deviation point set and elevation fixed point set of the triangle according to a specific deviation point search method based on recursive triangle division;

[0042] S323, perform deduplication on the elevation fixed point set of the triangle , and remove the duplicate elevation fixed points;

[0043] S324, add all the point elements in the elevation deviation point set of the triangle and the deduplicated elevation fixed point set to the created terrain encryption point set;

[0044] S325, let , repeat steps S322 - S325, and perform deviation point search, fixed point deduplication, and encryption point insertion operations based on recursive triangle division on all triangles in the terrain model in sequence until all triangles are traversed;

[0045] S326, sequentially determine the triangle positions of each point in the terrain encryption point set in the global terrain model, and encrypt these points into the terrain model triangular network according to the Delaunay triangle principle.

[0046] Furthermore, the specific steps for deleting redundant points with an included angle less than the maximum normal vector included angle threshold in the normal vector redundant point deletion in S400 are as follows:

[0047] S401, let , as the vertex index of the currently traversed triangle, start processing from the first vertex;

[0048] S402, for the th triangle vertex in the terrain model , query all triangles with as the vertex through establishing a vertex triangle index table , where is the number of associated triangles;

[0049] S403, for each associated triangle , calculate its normal vector , assuming the triangle vertices are , , , the calculation formula is: , where, , , and normalize the normal vector;

[0050] S404. For the normal vectors of all associated triangles , calculate the included angle. The formula is: , where and , record the maximum value among all the included angles ;

[0051] S405. Compare the maximum normal vector included angle value with the preset threshold . If , it is determined that the vertex is located in a flat terrain area and its contribution to terrain undulation is small, then delete this vertex from the triangular mesh. If , retain the vertex ;

[0052] S406. Let , and repeat steps S402 - S405 until all triangle vertices in the terrain model are traversed.

[0053] Compared with the prior art, the proposed airborne LiDAR point cloud rapid thinning method for terrain model construction has the following beneficial effects:

[0054] First, by simultaneously setting a small grid size and a large grid size, the present invention performs double - grid thinning on the original ground point cloud, which not only effectively reduces the amount of point cloud data, but also provides rich terrain feature information for subsequent global terrain model construction. Compared with the traditional single - grid thinning method, the method of the present invention can more accurately retain the point cloud data in areas with severe terrain undulation, avoiding the loss of terrain features. At the same time, by thinning the point cloud with the large grid size to extract the point cloud boundary points and areas with severe terrain undulation as the seed points for subsequent global terrain model construction, the construction efficiency and accuracy of the terrain model are further improved. When dealing with large - scale airborne LiDAR point cloud data, the present invention can significantly improve the data processing speed and terrain model construction quality, providing a more efficient and accurate technical means for the fields of terrain mapping and urban planning.

[0055] Second, through recursive triangular partitioning and deviation point search, the present invention locally densifies the terrain model, ensuring a fine expression of terrain features. At the same time, by calculating the maximum normal vector angle threshold between triangles, redundant points are intelligently removed, further optimizing the point cloud data structure and reducing unnecessary computational effort. This not only improves the accuracy and detail expressiveness of the terrain model but also reduces the computational burden of subsequent terrain analysis and visualization applications. Compared with traditional point cloud thinning methods, the method of the present invention achieves a more compact and efficient expression of point cloud data while maintaining terrain features, laying a solid foundation for the further development of the terrain model construction field.

[0056] Other advantages, objectives, and features of the present invention will be described to some extent in the subsequent specification, and to some extent, will be obvious to those skilled in the art based on the study of the following text, or can be taught from the practice of the present invention. BRIEF DESCRIPTION OF THE DRAWINGS

[0057] To more clearly illustrate the technical solutions in the embodiments of the present invention or in the prior art, the following will briefly introduce the drawings required for use in the description of the embodiments or the prior art. Obviously, the following-described drawings are only some embodiments of the present invention, and those of ordinary skill in the art can obtain other drawings based on these drawings without creative efforts.

[0058] Figure 1 It is a flowchart of a fast point cloud thinning method for airborne LiDAR for terrain model construction;

[0059] Figure 2 It is a schematic diagram showing the relationship between normal vector change and terrain undulation;

[0060] Figure 3 It is a schematic diagram for detecting elevation deviation points in four-layer recursive triangular partitioning;

[0061] Figure 4 It is a schematic diagram of triangular partitioning;

[0062] Figure 5 It is a schematic diagram for encrypting height difference deviation points and their corresponding fixed points. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0063] To further elaborate on the technical means and effects adopted by the present invention to achieve the predetermined invention objective, the following, in combination with the accompanying drawings and preferred embodiments, details the specific implementation manners, structures, features, and their effects of the present invention as follows.

[0064] Embodiment 1:

[0065] Processing of airborne LiDAR point clouds of urban building complex terrains.

[0066] The central area of ​​a city contains complex features such as high-rise buildings, roads, squares and green spaces. It is necessary to use airborne LiDAR point clouds to construct a 3D terrain model. It is required to retain the key features of building outlines and road boundaries, and at the same time remove redundant point clouds to optimize model storage and computing efficiency.

[0067] Mesh parameter setting: Set small mesh size 0.05m, large grid size is 1m, where and The value of is determined according to the actual terrain accuracy requirements: 0.05m-0.2m, used to preserve the detailed features of the bottom of the building and the edge of the road; for 10-20 times of the original value, used to extract the main contour of the ground object. In specific implementation, the final value can be determined by point cloud density analysis or experimental comparison optimization. , The minimum and maximum values ​​in the direction are denoted as , , , , calculate the number of grid rows corresponding to the small grid and the large grid respectively and number of columns , the formula is: ,in Take separately and , initialize the two-dimensional matrix grid (size is ), the element value defaults to (indicates no point), traverse each point of the original point cloud , calculate the row and column number of the grid: , , if grid[ ][ ]=1, the current point The point number is stored in grid[ ][ ]; if grid[ ][ ]≠1, get the current reserved point of the grid Elevation and with the current point Elevation For comparison; if , then use the current point Replace grid[ ][ ]; if , no processing is done and the original reserved points are retained. and Repeat the above steps to obtain the small grid sparse point cloud (including the bottom of the building and the details of the road edge) and the large grid sparse point cloud (including the main outline of the urban objects). Figure 1 .

[0068] Boundary point extraction (based on Alphashapes algorithm) to thin out the point cloud of large grid Project to plane, get a two-dimensional point set Set Alphashapes radius parameter ,right Perform Delaunay triangulation to generate a triangulated network , for each edge of the triangulated network , calculate the radius of its circumscribed circle ;like , keep it as Alpha edge; otherwise, remove it; the triangle composed of three Alpha edges is Alpha triangle, keep it; otherwise, remove it, count all the edges (boundary edges) that belong only to a single Alpha triangle, and their endpoints constitute the boundary point set The endpoint of a boundary edge , represents the outer contour of the urban area (such as road red lines, building complex boundaries), such as Figure 2 .

[0069] Identification of areas with dramatic terrain fluctuations (based on normal vector angle) for large grid point clouds Each point , select the nearest 25 neighboring points , calculate the normal vector by principal component analysis : Calculate the covariance matrix of the neighborhood point set ,right Perform eigenvalue decomposition and take the minimum eigenvalue The corresponding eigenvector As the normal vector ,calculate With neighboring points The normal vector angle , and the weighted average is calculated by the inverse distance weighted method, the formula is: ,in for and If the plane distance Threshold ,determination The area is a drastic undulation (such as building facades, steps), such as Figure 2 .

[0070] Feature point extraction and triangular mesh construction: According to the mapping relationship between the large grid and the small grid, determine the row and column range of the small grid corresponding to the severely undulating area: , where are the row and column numbers of the severe area in the large grid. Extract the small grid points within this range as terrain feature points (such as building corners and balcony edges). Use the boundary point set and the terrain feature points as seed points, and construct a Delaunay triangular mesh through the point-by-point insertion method to form a global terrain model containing urban structure features.

[0071] Deviation point search (based on recursive triangle partitioning): For each triangle in the triangular mesh , calculate the area of the triangle using the vector cross product formula. The formula is: , where , ; Fit the triangle plane using the least squares method. The formula is: , where , , are the coefficients. The formula is: , , . Substitute the barycentric coordinates into the plane equation to obtain the model elevation , and obtain the actual point cloud elevation at the centroid from the point cloud data , calculate the elevation deviation . The formula is: . If the triangle area or the recursion count reaches , stop the recursive partitioning of the triangle; if Figure 3 , and if , divide the triangle into three sub-triangles from the triangle centroid. For example, Figure 4 , obtain new vertices and triangles. For the newly generated sub-triangles, repeat the above steps until there are no new deviation points. Record all deviation points and their corresponding fixed points (such as triangle vertices), insert them into the triangular mesh for encryption. For example, Figure 5 , update the terrain model, and repeat the deviation point search for the newly generated triangles until no new deviation points are generated.

[0072] Normal vector redundant point removal: Traverse the vertices of the triangular mesh , query all triangles with as the vertex through the vertex triangle index table . For each triangle , calculate the normal vector: where , and normalize . Calculate the maximum value of the included angle of all associated triangle normal vectors: If the threshold is 5°, it is determined that if it is located in a flat area (such as a road surface or a square), the vertex is deleted; otherwise, it is retained (such as the vertex of a building facade).

[0073] In summary, for the point cloud processing of urban building groups, the present invention retains the details of building corners and road edges through small grids, constructs the urban framework through large grids, combines the Alphashapes algorithm to extract the boundaries of building groups, uses the angle between normal vectors to identify the characteristic areas of building facades, refines the complex areas of building facades as needed during recursive encryption, reduces the calculation for flat road areas, and finally eliminates redundant points on the road surface through the normal vector threshold. This solution significantly reduces the data volume while ensuring the geometric accuracy of the urban three-dimensional model, is applicable to the fields of urban three-dimensional modeling and intelligent transportation planning, and effectively balances the needs of model detail retention and data reduction.

[0074] Embodiment 2:

[0075] Downsampling and modeling of airborne LiDAR point clouds for mountainous terrain.

[0076] The terrain of a certain mountainous area is complex, with steep slopes and areas with severe undulations such as ravines. It is necessary to use airborne LiDAR point clouds to construct a high-precision terrain model while reducing data redundancy to improve calculation efficiency.

[0077] Set the small grid size (used to capture details of building corners and road cracks), and the large grid size (used to quickly construct the overall urban framework). The steps for determining the small grid size and the large grid size are as follows: statistically calculate the average point spacing of the original point cloud. According to the detail retention requirements, set to 1.5 times the average point spacing, and set the large grid size based on a 20-fold ratio of the small grid size. Considering both efficiency and the integrity of the main outline, verify through point cloud density analysis or experimental comparison to determine the final value. And calculate the extreme values of the point cloud data in , directional extreme values , , , , and calculate the number of grid rows and columns: , where is the current grid size (small grid or large grid). Initialize the two-dimensional matrix grid, traverse the original point cloud, and calculate the grid row and column numbers for each point : , if the grid is empty ( ), then deposit the point; otherwise, compare the elevations and retain the point with the lower elevation. Respectively use and Repeat the above steps to obtain the thinned point clouds of the small grid and the large grid.

[0078] Project the large grid point cloud onto a plane, set the radius parameter using the Alphashapes algorithm , screen the Alpha edges and Alpha triangles of the Delaunay triangulation, count the edges that belong to only a single triangle as boundary edges, and the endpoints are the boundary points , terrain undulation identification: for each point in the large grid point cloud , select the 25 nearest neighborhood points , calculate the normal vector through the principal component analysis (the eigenvector corresponding to the minimum eigenvalue), calculate the weighted average of the normal vector angles within the neighborhood: , if threshold , it is determined as a severely undulating area. According to the mapping relationship between the large and small grids, extract the small grid points corresponding to the severe area as terrain feature points, and together with the boundary points, use them as seed points to construct a Delaunay triangulation by the point-by-point insertion method.

[0079] For each triangle △ABC in the triangulation, the formula is: , let the barycentric coordinates of the triangle be , and the calculation formula is: , , , substitute the barycentric coordinates into the plane equation to obtain the model elevation , obtain the actual point cloud elevation at the barycenter from the point cloud data , calculate the elevation deviation , and the formula is: ; if and the recursive termination condition (area threshold or maximum recursion times) is not reached, divide the triangle into three sub-triangles, and repeat the above steps until there are no new deviation points.

[0080] Traverse the vertices of the triangulation , query all the triangles with as vertices, calculate the normal vectors of each triangle , calculate the normal vector angle, and the formula is: , compare the maximum normal vector angle value with the preset threshold , if , delete the vertex and retain the points with significant terrain features.

[0081] In summary, for the airborne LiDAR point cloud processing of mountainous terrain, the present invention realizes the stratification of point cloud density through multi-scale grid preliminary screening, accurately extracts the feature points of the boundary and the severely undulating area by using the Alphashapes algorithm and normal vector analysis, constructs a global terrain model by combining the Delaunay triangulation. During the recursive densification process, through the dual threshold control of triangle area and elevation deviation, it ensures that the details of steep slopes and gullies are effectively retained, while the elimination of redundant points with normal vectors eliminates the redundant points in flat areas, avoiding the damage to terrain features caused by traditional thinning, and improving the modeling efficiency of complex terrain through hierarchical processing, which is applicable to scenarios with high requirements for terrain details in mountain disaster monitoring and forest resource surveys.

[0082] The above are only the preferred embodiments of the present invention, and do not impose any form of limitation on the present invention. Although the present invention has been disclosed above with preferred embodiments, it is not intended to limit the present invention. Any person skilled in the art can make some changes or modifications to the equivalent embodiments with equivalent changes by using the technical content disclosed above within the scope of the technical solution of the present invention. However, as long as it does not depart from the content of the technical solution of the present invention, any brief modification, equivalent change and modification made to the above embodiments according to the technical essence of the present invention still fall within the scope of the technical solution of the present invention.

Claims

1. An airborne LiDAR point cloud rapid thinning method for terrain model construction, characterized in that The specific steps of the thinning method are as follows: S100, multi-scale grid point cloud preliminary screening: Set the small grid size and the large grid size, and perform grid thinning on the original ground point cloud respectively to obtain the thinned point cloud with the small grid size and the thinned point cloud with the large grid size; S200, global terrain model construction: Extract the point cloud boundary points and the areas with drastic terrain undulations based on the thinned point cloud with the large grid size. According to the mapping relationship between the large grid size and the small grid size, extract the small grid points corresponding to the drastic areas from the index matrix as terrain feature points. Use the extracted point cloud boundary points and terrain feature points as seed points, and construct a Delaunay triangulation by the point-by-point insertion method to form a global terrain model; S300, recursive encryption and iterative optimization: Search for deviation points based on recursive triangle division for each triangle in the global terrain model, insert the searched deviation points and their corresponding fixed points into the global terrain model for encryption. For the new triangles where the deviation points are located in the encrypted global terrain model, re-search for and encrypt the deviation points based on recursive triangle division, and iterate until no new deviation points are generated; S400, removal of normal vector redundant points: Traverse each vertex in the triangulation, calculate the normal vector of the triangle where it is located, set the maximum normal vector angle threshold between triangles, and delete the redundant points smaller than the maximum normal vector angle threshold from the triangulation. Use the remaining triangulation vertices to form the thinned ground point cloud.

2. A fast thinning method for airborne LiDAR point clouds for terrain model construction according to claim 1, characterized in that, In the S100, the specific steps of grid thinning in the multi-scale grid point cloud preliminary screening are as follows: S101, set the small grid size as , and the large grid size as . Obtain the minimum and maximum values of the point cloud data in the and directions, denoted as , , , ; S102, calculate the number of grid rows and columns corresponding to the small grid and the large grid respectively and columns , and the calculation formula is: , where is the currently used grid size, represents the ceiling operation; S103, Initialize the grid storage structure to create a two-dimensional matrix , whose size is , is the number of grid rows, is the number of grid columns. Set the value of all elements in the array to 1, indicating that the grid does not currently hold any points; S104, traverse the original point cloud and perform the following processing on each point in the original ground point cloud in turn Perform the following processing; Calculate the current point The row and column numbers of the grid where it is located , the formula is: , where Indicates the floor operation, check grid value, grid is a two-dimensional matrix;​ If grid = 1, store the point number of the current point into grid ;​​ If grid ≠ 1, obtain the elevation of the current reserved point and compare it with the elevation of the current point . If , replace the original value in grid with the point number of the current point . If , do nothing and keep the original reserved point; ; if , do not process and keep the original reserved point;​​ S105, traverse the grid storage structure grid, extract the points corresponding to the non-1 point numbers stored in it from the original point cloud, and form a new point cloud data, which is the thinned point cloud under the current grid size; S106. Replace the grid sizes in steps S102 - S105 with a small grid size and a large grid size respectively, and repeat steps S102 - S105 to obtain a downsampled point cloud with a small grid size and a downsampled point cloud with a large grid size.

3. The airborne LiDAR point cloud rapid thinning method for terrain model construction according to claim 2, characterized in that, The S100, the small grid size in the preliminary screening of the multi-scale grid point cloud is 0.05 m, and the large grid size is 1 m, and the decimated point cloud points of the large grid size are 1 / 400 of the small grid size.

4. A fast thinning method for airborne LiDAR point clouds for terrain model construction according to claim 1, characterized in that, In the S200, the specific steps of using the Alphashapes algorithm to extract boundary points in the global terrain model construction are as follows: S201, thin out the point cloud with a large grid size From three-dimensional space Project onto the horizontal plane to obtain a two-dimensional point set , according to the large grid size Set the radius parameter of Alphashapes , with a value of ; S202, perform Delaunay triangulation on the two-dimensional point set to generate a Delaunay triangulation network DT, and perform the following screening on the Alpha edges and Alpha triangles; For each edge in the Delaunay triangulation , calculate its corresponding circumradius . If , then this edge is an Alpha edge and is retained; otherwise, it is removed. The triangle composed of three Alpha edges is an Alpha triangle and is retained; otherwise, it is excluded; S203. In the remaining Alpha triangles, count all the edges that belong to only a single Alpha triangle, which are the boundary edges, and use the endpoints of all the boundary edges to form a boundary point set is the endpoint of a certain boundary edge ; S204, calculate the minimum bounding rectangle of the extracted boundary point set B, and verify whether it completely covers the target survey area. If the survey area is not covered, increase the α value, and repeat the steps of S202 - S203 until the inclusiveness requirement is met.

5. A rapid thinning method for airborne LiDAR point clouds for terrain model construction according to claim 1, characterized in that In the S200, the specific steps of determining the areas with drastic terrain undulations in the global terrain model construction are as follows: (1) For each point in the downsampled point cloud with a large grid size , select its 25 nearest points to form a neighborhood , and estimate the normal vector of the point based on principal component analysis , and regularize the calculated normal vector ; ​ (2)For the point and its neighborhood For each point in , calculate the angle between its normal vector and . The formula is: . Use the inverse distance weighted method to calculate the weighted average of the normal vector angles in the neighborhood of point . Let the plane distance between point and point be . The calculation formula is: ; (3) Set the threshold of the weighted average of the normal vector angles according to the actual terrain conditions = , traverse the downsampled point cloud with a large grid size for each point , the weighted average of its normal vector angles is greater than the threshold , and the grid area where the point is located is considered a region with severe terrain undulations.

6. A rapid thinning method for airborne LiDAR point clouds for terrain model construction according to claim 5, characterized in that The S200, the point estimated based on principal component analysis in global terrain model construction Normal vector , principal component analysis calculates the covariance matrix of the point set within the neighborhood , performs eigenvalue decomposition on the covariance matrix to obtain three eigenvalues and , as well as the corresponding eigenvectors , where, and , the eigenvector corresponding to the minimum eigenvalue is the normal vector of the point .​ 7. A rapid thinning method for airborne LiDAR point clouds for terrain model construction according to claim 1, characterized in that, For the S200, the determination of the three-dimensional points in the terrains with drastic undulations in the small grid size thinning point cloud during the global terrain model construction. Let the row number of the grid in the terrain with drastic undulations in the large grid be , the column number be , the large grid size be , the small grid size be , and the calculation formulas for the corresponding small grid row number range and column number range are as follows: , where represents rounding down, and represents rounding up.

8. A fast thinning method for airborne LiDAR point clouds for terrain model construction according to claim 1, characterized in that In the S300, the specific steps of searching for deviation points based on recursive triangle division of triangles in the recursive encryption and iterative optimization are as follows: S311. For the constructed Delaunay triangulation, obtain the vertex coordinates of each triangle. Let the three vertices of the triangle be , , . Set the recursive termination conditions, the triangle area threshold and the maximum number of recursions . S312, calculate the area of the triangle using the vector cross product formula, and the formula is: , where , ; Fitting a triangular plane using the least squares method, the formula is: , where , , are coefficients; S313. Set the centroid coordinates of the triangle as , and the calculation formula is: , , . Substitute the centroid coordinates into the plane equation to obtain the model elevation . Obtain the actual point cloud elevation at the centroid from the point cloud data . Calculate the elevation deviation . The formula is: ; S314, if the area of the triangle or the number of recursions reaches , stop the recursive partitioning of the triangle; if , divide the triangle into three sub - triangles evenly from the centroid of the triangle, obtain new vertices and triangles, and for the newly generated sub - triangles, repeat steps S312 - S314 until the recursive condition is no longer satisfied; S315, record all points with elevation deviation values greater than as deviation points.

9. A fast thinning method for airborne LiDAR point clouds for terrain model construction according to claim 1, characterized in that In the S300, the encryption of the global terrain model in the recursive encryption and iterative optimization: S321, let , create a new empty set of terrain encryption points; S322, for the th triangle in the triangle set, respectively obtain the elevation deviation point set and the elevation fixed point set of the triangle according to a specific deviation point search method based on recursive triangle partitioning; S323. Dedup the set of elevation fixed points of the triangle and remove the duplicate elevation fixed points therefrom; S324, add all the point elements in the elevation deviation point set of the triangle and the de-duplicated elevation fixed point set to the created terrain encryption point set; S325, let , repeat steps S322 - S325 to perform operations of deviation point search based on recursive triangle division, fixed point duplicate removal, and encrypted point insertion on all triangles in the terrain model in sequence until all triangles are traversed; S326, sequentially determine the positions of the points in the terrain encryption point set in the triangles in the global terrain model, and encrypt these points into the terrain model triangulation according to the Delaunay triangle principle.

10. A fast thinning method for airborne LiDAR point clouds for terrain model construction according to claim 1, characterized in that, In the S400, the specific steps of deleting the redundant points smaller than the maximum normal vector angle threshold from the triangulation in the removal of normal vector redundant points are as follows: S401, let , as the index of the triangle vertex being traversed currently, start processing from the first vertex; S402. For the th triangular vertex in the terrain model, by establishing a vertex triangle index table, query all triangles with as vertices, where is the number of associated triangles; Among them, is the number of associated triangles; S403. For each associated triangle , calculate its normal vector . Let the vertices of the triangle be , , . The calculation formula is: , where , . And normalize the normal vector. S404, for the normal vectors of all associated triangles , calculate the included angle, and the formula is: , where and , record the maximum value among all the included angles ; S405, take the maximum normal vector angle value and compare it with the preset threshold . If , it is determined that the vertex is located in a flat terrain area and its contribution to terrain undulation is small, so delete this vertex from the triangular mesh. If , retain the vertex ; S406, let , and repeat steps S402 - S405 until all triangular vertices in the terrain model are traversed.

Images