A fast thinning method for airborne LiDAR point clouds for terrain model construction

Abstract

The present invention discloses a rapid thinning method for airborne LiDAR point clouds for terrain model construction, and relates to the technical field of airborne LiDAR point cloud data processing. The specific steps of the thinning method are: S100, multi-scale grid point cloud initial screening: setting a small grid size and a large grid size. The present invention performs double gridding and thinning on the original ground point cloud by simultaneously setting the small grid size and the large grid size, 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 gridding and thinning method, the method of the present invention can more accurately retain point cloud data in areas with severe terrain fluctuations, avoiding the loss of terrain features. At the same time, by thinning the point cloud with a large grid size, point cloud boundary points and areas with severe terrain fluctuations are extracted as seed points for subsequent construction of the global terrain model, further improving the construction efficiency and accuracy of the terrain model.

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 clouds 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 topographic surveying and mapping, urban planning, and forestry resource surveys. 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 accurate terrain models. 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, alignment, and thinning. Among them, point cloud thinning is an important means to reduce data volume and improve processing efficiency. While maintaining terrain characteristics, it reduces the computational complexity of subsequent terrain modeling by reducing the number of point clouds.

[0003] Although the existing airborne LiDAR point cloud thinning technology has improved data processing efficiency to a certain extent, it still has 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 cluster analysis, although they can better preserve 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 solved urgently 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 existing technology 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 preservation 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 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: a fast thinning method for airborne LiDAR point clouds for terrain model construction, wherein the specific steps of the thinning method are as follows:

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

[0007] S200, global terrain model construction: Based on the large grid size sparse point cloud, point cloud boundary points and areas with severe terrain fluctuations are extracted. Based on the mapping relationship between the large grid size and the small grid size, the small grid points corresponding to the severe areas are extracted from the index matrix as terrain feature points. The extracted point cloud boundary points and terrain feature points are used as seed points, and a Delaunay triangulation is constructed through point-by-point interpolation to form a global terrain model.

[0008] S300, recursive encryption and iterative optimization: performing a deviation point search based on recursive triangulation on each triangle in the global terrain model, inserting the searched deviation points and their corresponding fixed points into the global terrain model for encryption, and re-performing a deviation point search and encryption based on recursive triangulation on the new triangles in the encrypted global terrain model where the deviation points are located, and iterating until no new deviation points are generated;

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

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

[0011] S101, set the small grid size to , the large grid size is , get point cloud data in 、 The minimum and maximum values ​​in the direction are recorded as 、 、 、 ;

[0012] S102, respectively calculate the number of grid rows corresponding to the small grid and the large grid and number of columns , the calculation formula is: ,in, is the currently used grid size, Indicates rounding up operation;

[0013] S103, initializing 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. Setting the value of all elements in the array to 1 indicates that the grid currently does not retain any points.

[0014] S104, traverse the original point cloud and perform a search on each point in the original ground point cloud Perform the following processing;

[0015] Calculate the current point The row and column number of the grid , the formula is: ,in, Indicates rounding down operation, check grid[ ][ ] value, grid is a two-dimensional matrix;

[0016] If grid[ ][ ]=1, the current point The point number is stored in grid[ ][ ];

[0017] If grid[ ][ ]≠1, get the current reserved point of the grid Elevation and with the current point Elevation To compare; if , then use the current point Replace the dot grid[ ][ ] in the original value; if , no processing is done and the original reserved point is retained;

[0018] S105, traversing the grid storage structure grid, extracting the points corresponding to the non-1 point numbers stored therein from the original point cloud, and forming new point cloud data, i.e., the sparse point cloud at the current grid size;

[0019] S106, change the grid size of steps S102-S105 Replace with smaller grid sizes and large grid sizes , repeat steps S102-S105 to obtain a small grid size sparse point cloud and a large grid size sparse point cloud.

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

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

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

[0023] For each edge in the Delaunay triangulation , calculate the corresponding circumscribed circle radius ,like , then the edge is an Alpha edge and is retained; otherwise, it is removed;

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

[0025] S203: In the retained Alpha triangle, count all edges that belong only to a single Alpha triangle, i.e., the boundary edges, and use the endpoints of all boundary edges to form a boundary point set. The endpoint of a boundary edge ;

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

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

[0028] (1) Thinning point clouds for large grid sizes Every point in , select the 25 nearest points to form a neighborhood , based on principal component analysis to estimate the points Normal vector , and the calculated normal vector Perform regularization processing;

[0029] (2) For point and its neighborhood Every point within , , calculate its normal vector and The angle between , the formula is: , since the normal vector has been regularized, ,so , using the inverse distance weighted method to calculate the point Weighted average of the normal vector angles within the neighborhood , set up a point with dot The plane distance between , the calculation formula is: ;

[0030] (3) According to the actual terrain conditions, set the threshold of the weighted average of the normal vector angle = , traverse the large grid size to extract the thin point cloud Every point in , the weighted average of the normal vector angles Greater than threshold , it is considered that the grid area where the point is located is an area with severe terrain fluctuations.

[0031] Furthermore, in said S200, the global terrain model is constructed based on the principal component analysis estimation point Normal vector , principal component analysis calculates the covariance matrix of the point set in the neighborhood , perform eigenvalue decomposition on the covariance matrix and obtain three eigenvalues and , and the corresponding eigenvector , among which, and , the minimum eigenvalue The corresponding eigenvector That is the point Normal vector .

[0032] Furthermore, in the above S200, the three-dimensional points in the area with severe terrain fluctuations in the small grid size sparse point cloud in the global terrain model construction are determined, and the row number of the grid in the area with severe terrain fluctuations in the large grid is set as , column number is , the large grid size is , the small grid size is , corresponding to the small grid row number range , column number range The calculation formula is: ,in, Indicates rounding down. Indicates rounding up.

[0033] Furthermore, the specific steps of searching for deviation points of recursive triangulation of triangles in the recursive encryption and iterative optimization in S300 are as follows:

[0034] S311, for the constructed Delaunay triangulation, obtain the vertex coordinates of each triangle, and set the three vertices of the triangle to be 、 、 , set the recursive termination condition, triangle area threshold and the maximum number of recursions ;

[0035] S312, use the vector cross product formula to calculate the area of ​​the triangle, the formula is: ,in , ; Use the least squares method to fit the triangular plane. The formula is: ,in 、 、 is the coefficient;

[0036] S313, let the coordinates of the triangle barycenter be , the calculation formula is: , , , substitute the coordinates of the center of gravity into the plane equation to obtain the model elevation , get the actual point cloud elevation at the center of gravity from the point cloud data , calculate the elevation deviation , the formula is: ;

[0037] S314, if the area of ​​the triangle Or the number of recursions reaches , stop the recursive division of the triangle; if , divide the triangle into three sub-triangles from the triangle centroid, obtain new vertices and triangles, and repeat steps S312,-S314 for the newly generated sub-triangles until the recursive condition is no longer met;

[0038] S315, all elevation deviation values ​​greater than The point is recorded as the deviation point.

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

[0040] S321, order , create a new empty terrain encryption point set;

[0041] S322, for the triangle set triangles , according to a specific deviation point search method based on recursive triangulation, the elevation deviation point set and elevation fixed point set of the triangle are obtained respectively;

[0042] S323, triangle The set of elevation fixed points is deduplicated to remove duplicate elevation fixed points;

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

[0044] S325, order , repeating steps S322-S325, performing deviation point search, fixed point deduplication and encryption point insertion operations based on recursive triangulation on all triangles in the terrain model, until all triangles are traversed;

[0045] S326 , sequentially determining the triangular position of each point in the terrain densification point set in the global terrain model, and densifying these points into the terrain model triangulation network according to the Delaunay triangle principle.

[0046] Furthermore, the specific steps of removing redundant points smaller than the maximum normal vector angle threshold from the triangulated network in S400 are as follows:

[0047] S401, order , as the vertex index of the currently traversed triangle, starting from the first vertex;

[0048] S402, the terrain model triangle vertices By establishing a vertex triangle index table, query all A triangle with vertices ,in is the number of associated triangles;

[0049] S403, for each associated triangle , calculate its normal vector , let the triangle vertices be 、 、 , the calculation formula is: ,in, , , and normalize the normal vector;

[0050] S404, normal vectors of all associated triangles , calculate the angle, the formula is: ,in and , record the maximum value of all angles ;

[0051] S405, the maximum normal vector angle value With preset threshold Compare, if , determine the vertex Located in a flat terrain area, its contribution to the terrain undulation is small, so the vertex is deleted from the triangulation network. , retain the vertex ;

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

[0053] Compared with the existing technology, this method for rapid thinning of airborne LiDAR point clouds for terrain model construction has the following beneficial effects:

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

[0055] 2. The present invention locally encrypts the terrain model through recursive triangulation and deviation point search, ensuring the 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 calculations. This not only improves the accuracy and detail expression 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 further development in the field of terrain model construction.

[0056] Other advantages, objects and features of the present invention will be described in part in the following description and, in part, will be apparent to those skilled in the art based on an examination of the following or may be learned from the practice of the invention. BRIEF DESCRIPTION OF THE DRAWINGS

[0057] To more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the following briefly introduces the drawings required for use in the embodiments or the description of the prior art. Obviously, the drawings described below are only some embodiments of the present invention. Those skilled in the art can also derive other drawings based on these drawings without inventive effort.

[0058] Figure 1 A flowchart of a fast thinning method for airborne LiDAR point clouds for terrain model construction;

[0059] Figure 2 Schematic diagram of the relationship between normal vector change and terrain undulation;

[0060] Figure 3 Schematic diagram of elevation deviation point detection using four-layer recursive triangulation;

[0061] Figure 4 This is a schematic diagram of triangulation;

[0062] Figure 5 This is an encrypted schematic diagram of the height difference deviation points and their corresponding fixed points. DETAILED DESCRIPTION

[0063] In order to further illustrate the technical means and effects adopted by the present invention to achieve the predetermined purpose of the invention, the specific implementation methods, structures, features and effects of the present invention are described in detail below in conjunction with the accompanying drawings and preferred embodiments.

[0064] Example 1:

[0065] Airborne LiDAR point cloud processing of urban building terrain.

[0066] The central area of ​​a certain 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. The key features of building outlines and road boundaries must be retained, while redundant point clouds must be eliminated to optimize model storage and computing efficiency.

[0067] Grid parameter setting: set small grid 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 is 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. And obtain point cloud data in 、 The minimum and maximum values ​​in the direction are recorded 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 To compare; if , then use the current point Replace the dot grid[ ][ ] in the original value; 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, and obtain 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 , retain as Alpha edge; otherwise, remove; the triangle composed of three Alpha edges is an Alpha triangle and is retained; otherwise, remove, count all 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 line, building group boundary), such as Figure 2 .

[0069] Identification of areas with severe 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 calculate the weighted average using the inverse distance weighted method, the formula is: ,in for and If the plane distance Threshold ,determination The area is a drastically undulating area (such as building facades, steps), such as Figure 2 .

[0070] Feature point extraction and triangulation network construction determine the range of small grid rows and columns corresponding to the violently undulating area based on the mapping relationship between large grids and small grids: ,in The row and column numbers of the intense area in the large grid are extracted, and the small grid points in the range are used as terrain feature points (such as building corners and balcony edges). The terrain feature points are used as seed points, and the Delaunay triangulation is constructed by point-by-point interpolation method to form a global terrain model that includes urban structural features.

[0071] Deviation point search (based on recursive triangulation) for each triangle in the triangle network , use the vector cross product formula to calculate the area of ​​the triangle, the formula is: ,in , ; Use the least squares method to fit the triangular plane. The formula is: ,in 、 、 is the coefficient, the formula is: , , , substitute the coordinates of the center of gravity into the plane equation to obtain the model elevation , get the actual point cloud elevation at the center of gravity from the point cloud data , calculate the elevation deviation , the formula is: , if the area of ​​the triangle Or the number of recursions reaches , stop the recursive division of the triangle; if Figure 3 ,like , divide the triangle into three sub-triangles from the triangle center, such as Figure 4 , get new vertices and triangles, repeat the above steps for the newly generated sub-triangles until there are no new deviation points, and record all The deviation points and their corresponding fixed points (such as triangle vertices) are inserted into the triangulated network for encryption, such as 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 elimination: traversing the triangulated network vertices , query all the vertex triangle index tables A triangle with vertices , for each triangle , calculate the normal vector: in , and Normalize and calculate the maximum angle of all associated triangle normal vectors: like Threshold 5°, judgment If the vertex is located in a flat area (such as a road surface or a square), delete it; otherwise, retain it (such as a vertex on a building facade).

[0073] In summary, for point cloud processing of urban building complexes, the present invention uses small grids to retain building corners and road edge details, uses large grids to construct the city framework, combines the Alphashapes algorithm to extract the boundaries of building complexes, uses normal vector angles to identify feature areas of building facades, refines complex areas of building facades as needed during recursive encryption, reduces calculations for flat road areas, and finally eliminates redundant road surface points through normal vector thresholds. This solution significantly reduces the amount of data while ensuring the geometric accuracy of urban three-dimensional models. It is suitable for urban three-dimensional modeling and smart transportation planning, and effectively balances the needs of retaining model details and streamlining data.

[0074] Example 2:

[0075] Airborne LiDAR point cloud thinning and modeling of mountainous terrain.

[0076] A certain mountainous area has complex terrain, with steep slopes, gullies and violent undulations. It is necessary to use airborne LiDAR point clouds to build a high-precision terrain model while reducing data redundancy to improve computing efficiency.

[0077] Set the small grid size (used to capture details of building corners and road cracks), large grid size (used to quickly build the overall framework of the city), where the steps for determining the small and large grid sizes are: counting the average point spacing of the original point cloud, and according to the detail retention requirements, Set to 1.5 times the average point spacing, large grid size Based on the 20-fold scale setting of the small grid size, taking into account both efficiency and the integrity of the main body outline, the final value is determined through point cloud density analysis or experimental comparison. 、 Extreme value of direction 、 、 、 , calculate the number of grid rows and columns: ,in Initialize the two-dimensional matrix grid for the current grid size (small grid or large grid), traverse the original point cloud, and for each point Calculate its grid row and column numbers: , if the grid is empty ( ), then store the point; otherwise, compare the elevations and keep the lower elevation point, respectively. and Repeat the above steps to obtain the sparse point clouds of small grids and large grids.

[0078] Project the large mesh point cloud onto Plane, using Alphashapes algorithm to set radius parameters , filter the Alpha edges and Alpha triangles of the Delaunay triangulation, count the edges belonging to only a single triangle as boundary edges, and the endpoints are the boundary points , identification of areas with dramatic terrain fluctuations: for each point in the large grid point cloud , select the nearest 25 neighboring points , calculate the normal vector by principal component analysis (The eigenvector corresponding to the minimum eigenvalue), calculate the weighted average of the normal vector angles in the neighborhood: ,like Threshold , it is determined to be a violently undulating area. According to the mapping relationship between large and small grids, the small grid points corresponding to the violent area are extracted as terrain feature points, and together with the boundary points, they are used as seed points. The Delaunay triangulation is constructed by the point-by-point insertion method.

[0079] For each triangle △ABC in the triangulation network, the formula is: , let the coordinates of the triangle barycenter be , the calculation formula is: , , , substitute the coordinates of the center of gravity into the plane equation to obtain the model elevation , get the actual point cloud elevation at the center of gravity from the point cloud data , calculate the elevation deviation , the formula is: ;like If the recursive termination condition (area threshold or maximum number of recursions) is not reached, the triangle is divided into three sub-triangles and the above steps are repeated until there are no new deviation points.

[0080] Traversing the triangulation vertices , query all For the triangle with vertices, calculate the normal vector of each triangle , calculate the normal vector angle, the formula is: , the maximum normal vector angle value With preset threshold Compare, 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 point cloud density stratification through multi-scale grid initial screening, uses the Alphashapes algorithm and normal vector analysis to accurately extract feature points of boundaries and sharply undulating areas, and combines the Delaunay triangulation to construct a global terrain model. During the recursive encryption process, dual threshold control of triangle area and elevation deviation is used to ensure that the details of steep slopes and gullies are effectively retained, and the removal of redundant normal points eliminates redundant points in flat areas, avoiding the destruction of terrain features by traditional thinning. The efficiency of complex terrain modeling is improved through hierarchical processing, which is suitable for scenarios with high requirements for terrain details such as mountain disaster monitoring and forestry resource surveys.

[0082] The above description is merely a preferred embodiment of the present invention and does not constitute any form of limitation to the present invention. Although the present invention has been disclosed as above in terms of a preferred embodiment, it is not intended to limit the present invention. Any person skilled in the art can, without departing from the scope of the technical solution of the present invention, make some changes or modifications to equivalent embodiments using the technical contents disclosed above. However, any brief modifications, equivalent changes and modifications made to the above embodiments based on the technical essence of the present invention without departing from the content of the technical solution of the present invention are still within the scope of the technical solution of the present invention.

Claims

1. A fast thinning method for airborne LiDAR point clouds for terrain model construction, characterized by: The specific steps of the thinning method are: S100, multi-scale grid point cloud initial screening: Set the small grid size and large grid size, and perform grid thinning on the original ground point cloud to obtain the small grid size thinned point cloud and the large grid size thinned point cloud respectively; S200, global terrain model construction: Based on the large grid size sparse point cloud, point cloud boundary points and areas with severe terrain fluctuations are extracted. Based on the mapping relationship between the large grid size and the small grid size, the small grid points corresponding to the severe areas are extracted from the index matrix as terrain feature points. The extracted point cloud boundary points and terrain feature points are used as seed points, and a Delaunay triangulation is constructed through point-by-point interpolation to form a global terrain model. S300, recursive encryption and iterative optimization: performing a deviation point search based on recursive triangulation on each triangle in the global terrain model, inserting the searched deviation points and their corresponding fixed points into the global terrain model for encryption, and re-performing a deviation point search and encryption based on recursive triangulation on the new triangles in the encrypted global terrain model where the deviation points are located, and iterating until no new deviation points are generated; The specific steps of searching for deviation points of the recursive triangulation of the triangle are: S311, for the constructed Delaunay triangulation, obtain the vertex coordinates of each triangle, and set the three vertices of the triangle as A(x a ,y a ,z a )、B(x b ,y b ,z b )、C(x c ,y c ,z c ), set the recursive termination condition, the triangle area threshold S threshold and the maximum number of recursions N max ; S312, use the vector cross product formula to calculate the area of ​​the triangle, the formula is: in The least squares method is used to fit the triangular plane. The formula is: z = ax + by + c, where a, b, and c are coefficients. S313, let the coordinates of the triangle center be G(x g ,y g ,z g ), the calculation formula is: Substitute the coordinates of the center of gravity into the plane equation to obtain the model elevation z model , obtain the actual point cloud elevation z at the center of gravity from the point cloud data point , calculate the elevation deviation Δh, the formula is: Δh=|z model -z point |; S314, if the area of ​​the triangle S≤S threshold Or the number of recursions reaches N max , stop the recursive division of the triangle; if Δh>Δh threshold , divide the triangle into three sub-triangles from the triangle centroid, obtain new vertices and triangles, and repeat steps S312,-S314 for the newly generated sub-triangles until the recursive condition is no longer met; S315: All elevation deviations greater than Δh tresold The point is recorded as the deviation point; Encryption of the global terrain model: S321, let i = 0, create a new empty terrain encryption point set; S322, for the i-th triangle T in the triangle set i , according to a specific deviation point search method based on recursive triangulation, the elevation deviation point set and elevation fixed point set of the triangle are obtained respectively; S323, triangle T i The set of elevation fixed points is deduplicated to remove duplicate elevation fixed points; S324, triangle T i Add all the point elements in the elevation deviation point set and the deduplicated elevation fixed point set to the created terrain encryption point set; S325, let i = i + 1, repeat steps S322-S325, and perform deviation point search, fixed point deduplication, and encryption point insertion operations based on recursive triangulation on all triangles in the terrain model until all triangles are traversed; S326, sequentially determining the triangular position of each point in the terrain densification point set in the global terrain model, and densifying these points into the terrain model triangulation network according to the Delaunay triangulation principle; S400, normal vector redundant point elimination: traverse each vertex in the triangulation network, calculate the normal vector of the triangle in which it is located, set the maximum normal vector angle threshold between triangles, delete the redundant points smaller than the maximum normal vector angle threshold from the triangulation network, and use the remaining triangulation network vertices to form a thinned ground point cloud.

2. The method for rapid thinning of airborne LiDAR point clouds for terrain model construction according to claim 1, characterized in that: The specific steps of grid thinning in the multi-scale grid point cloud initial screening in S100 are: S101, set the small grid size to g1, the large grid size to g2, obtain the minimum and maximum values ​​of the point cloud data in the x and y directions, recorded as x min 、x max 、y min 、y max ; S102, respectively calculate the number of grid rows L and the number of columns N corresponding to the small grid and the large grid, and the calculation formula is: Where g is the current grid size, Indicates rounding up operation; S103, initializing the grid storage structure to create a two-dimensional matrix grid, whose size is L×N, where L is the number of grid rows and N is the number of grid columns. The values ​​of all elements in the array are set to 1, indicating that the grid currently does not retain any points; S104, traverse the original point cloud and perform a search on each point P in the original ground point cloud. i (x i ,y i ,z i ) perform the following processing; Calculate the current point P i The row and column number of the grid (l i ,n i ), the formula is: in, Indicates rounding down operation, check grid[l i ][n i ] value, grid is a two-dimensional matrix; If grid[l i ][n i ]=1, the current point P i The point number is stored in grid[l i ][n i ]; If grid[l i ][n i ]≠1, get the current reserved point P of the grid j The elevation z j and with the current point P i The elevation z i Compare; if z i <z j , then use the current point P i Replace the dot grid[l i ][n i ] in the original value; if z i ≥z j , no processing is done and the original reserved point is retained; S105, traversing the grid storage structure grid, extracting the points corresponding to the non-1 point numbers stored therein from the original point cloud, and forming new point cloud data, i.e., the sparse point cloud at the current grid size; S106, replacing the grid size g in steps S102-S105 with a small grid size g1 and a large grid size g2 respectively, and repeating steps S102-S105 to obtain a small grid size sparse point cloud and a large grid size sparse point cloud.

3. The method for rapid thinning of airborne LiDAR point clouds for terrain model construction according to claim 2, characterized in that: In the S100 , in the initial screening of the multi-scale grid point cloud, the small grid size g1 is 0.05 m, the large grid size g2 is 1 m, and the number of thinned point cloud points of the large grid size is 1 / 400 of the small grid size.

4. The method for rapid thinning of airborne LiDAR point clouds for terrain model construction according to claim 1, characterized in that: The specific steps of extracting boundary points using the Alphashapes algorithm in the above-mentioned S200 in constructing the global terrain model are as follows: S201, project the large grid size sparse point cloud P2 from the three-dimensional space (X, Y, Z) to the horizontal XY plane to obtain a two-dimensional point set Set the radius parameter α of Alphashapes according to the large grid size g2, and the value is (1.0~1.5); S202, for a two-dimensional point set Perform Delaunay triangulation to generate Delaunay triangulation network DT, and perform the following screening on Alpha edges and Alpha triangles; For each edge e in the Delaunay triangulation, calculate the corresponding circumscribed circle radius r e , if r e ≤α, then the edge is an Alpha edge and is retained; otherwise, it is removed; A triangle composed of three Alpha edges is an Alpha triangle and is retained; otherwise, it is discarded; S203: In the retained Alpha triangle, count all edges that belong only to a single Alpha triangle, i.e., boundary edges, and use the endpoints of all boundary edges to form a boundary point set B = {p∈P2|p is an endpoint of a boundary edge}; S204, calculate the minimum enclosing rectangle of the extracted boundary point set B, and verify whether it completely contains the target measurement area. If the measurement area is not covered, increase the α value and repeat steps S202-S203 until the inclusion requirement is met.

5. The method for rapid thinning of airborne LiDAR point clouds for terrain model construction according to claim 1, characterized in that: The specific steps of determining the area with severe terrain fluctuations in the construction of the global terrain model in S200 are as follows: (1) For each point P in the large grid size sparse point cloud P2 i , select the 25 nearest points to form the neighborhood N i , based on principal component analysis to estimate point P i Normal vector N i , and the calculated normal vector N i Perform regularization processing; (2) For point P i and its neighborhood N i Every point P in j , j∈N i , calculate its normal vector N i and N j The angle θ between ij , the formula is: Calculate point P using the inverse distance weighted method i Weighted average of the normal vector angles within the neighborhood Set point P i and point P j The plane distance between them is d ij , the calculation formula is: (3) According to the actual terrain conditions, the threshold value of the weighted average of the normal vector angle is set to T = 5°, and each point P in the thinned point cloud P2 with a large grid size is traversed. i , the weighted average of the normal vector angles If the value is greater than the threshold T, the grid area where the point is located is considered to be an area with severe terrain fluctuations.

6. The method for rapid thinning of airborne LiDAR point clouds for terrain model construction according to claim 5, characterized in that: In step S200, the global terrain model is constructed based on the principal component analysis to estimate the point P i Normal vector N i , principal component analysis calculates the covariance matrix C of the point set in the neighborhood, performs eigenvalue decomposition on the covariance matrix, and obtains three eigenvalues ​​λ1, λ2 and λ3, and the corresponding eigenvectors v1, v2, v3, where λ1≥λ2≥λ3, and the eigenvector v3 corresponding to the minimum eigenvalue λ3 is the point P i Normal vector N i .

7. The method for rapid thinning of airborne LiDAR point clouds for terrain model construction according to claim 1, characterized in that: In the above S200, the three-dimensional points in the area with severe terrain fluctuations in the sparse point cloud of small grid size in the construction of the global terrain model are determined. The row number of the grid in the area with severe terrain fluctuations in the large grid is set as 1. i,2 , column number is n i,2 , the large grid size is g2, the small grid size is g1, and the corresponding small grid row number range is [l min1 ,l max1 ]、Column number range[n min1 ,n max1 The calculation formula of ] is: in, Indicates rounding down. Indicates rounding up.

8. The method for rapid thinning of airborne LiDAR point clouds for terrain model construction according to claim 1, characterized in that: The specific steps of removing redundant points smaller than the maximum normal vector angle threshold from the triangulated network in S400 are as follows: S401, let i = 0, as the vertex index of the triangle currently traversed, and start processing from the first vertex; S402, for the i-th triangle vertex P in the terrain model, i (x i ,y i ,z i ), by establishing a vertex triangle index table, query all i The triangle with vertices {△ i1 ,△ i2 ,…,△ ik }, where k is the number of associated triangles; S403, for each associated triangle △ ij , calculate its normal vector N ij , let the triangle vertex be A(x a ,y a ,z a )、B(x b ,y b ,z b )、C(x c ,y c ,z c ), the calculation formula is: in, And normalize the normal vector; S404, for all associated triangle normal vectors {N i1 ,N i2 ,…,N ik }, calculate the angle, the formula is: θ mn =arccos(N im ·N in ), where m, n∈[1,k] and m≠n, record the maximum value θ among all angles max ; S405, the maximum normal vector angle value θ max Compared with the preset threshold of 5°, if θ max <5°, determine vertex P i Located in a flat terrain area, its contribution to the terrain undulation is small, and the vertex is deleted from the triangulation network. max ≥5°, retain vertex P i ; S406 , let i=i+1, and repeat steps S402 - S405 until all triangle vertices in the terrain model are traversed.

Images