A LiDAR Point Cloud Clustering and Simplification Method Considering Terrain Features

Abstract

The present invention discloses a LiDAR point cloud clustering and simplification method considering terrain features. Aiming at the problems existing in the current ground point cloud simplification methods, such as poor applicability in complex environments and easy loss of terrain detail features, the following solutions are proposed: First, the K-means algorithm is used to segment the point cloud into initial point cloud clusters, and then the point cloud clusters are further subdivided according to the terrain complexity information of each cluster; Next, the terrain feature points are identified by means of the point cloud normal vector information and the elevation difference of the edge points between adjacent clusters; Finally, by retaining the boundary feature points of the point cloud region, the boundary contraction of the original point cloud is prevented. In addition, the present invention is also experimentally compared with traditional methods. Under the same point cloud simplification ratio, the accuracy of the digital elevation model generated by the method of the present invention and the accuracy of its derivatives (including average slope and terrain roughness) are significantly better than those of traditional methods, and the terrain feature information is better retained.

Description

Technical Field

[0001] The present invention belongs to the technical field of point cloud simplification, and relates to a LiDAR point cloud clustering and simplification method considering terrain features. Background Art

[0002] Accurate and effective three-dimensional spatio-temporal information is an essential important support for major national needs such as the construction of new infrastructure, the construction of real-scene three-dimensional China, and natural resource management and monitoring. In recent years, the innovative development of various earth observation technologies has improved the perception ability of the whole space and the whole time domain. In particular, the three-dimensional point cloud data acquisition method represented by airborne lidar technology (LiDAR) provides a new technology for intelligent surveying and mapping. High-quality and refined point cloud data can accurately express various spatial elements of the natural ground surface. However, when the airborne LiDAR system acquires ground points with spatial information and attribute information, the data acquisition specification is generally designed according to the point density requirements in complex terrain areas, resulting in excessive redundancy of the acquired ground point cloud in other flat terrain areas, which seriously restricts the storage, transmission, and analysis efficiency of terrain information.

[0003] Point cloud simplification is a prerequisite for the efficient transmission and multi-scale application of massive airborne LiDAR ground point clouds. Therefore, how to achieve the automation and intelligence of massive ground point clouds with high complexity and polymorphism, and meet the requirements of high precision and high efficiency for geoscience analysis has become an urgent problem to be solved. At present, scholars have conducted in-depth research on the simplification of airborne LiDAR ground point clouds, and classified them into 5 categories according to the working principles of the current mainstream methods: random downsampling method, voxel grid downsampling method, curvature sampling method, point reduction method, and point addition method. The random downsampling method randomly selects a certain number of sampling points; the voxel grid downsampling method uses the original point cloud to construct a three-dimensional voxel grid, and takes the centroid of each voxel point set as the sampling point of the voxel.

[0004] The above two algorithms are simple and efficient, but the simplification results only reduce the point cloud density and are difficult to accurately express the terrain feature structure. The curvature sampling method first calculates the curvature value of each point cloud, and then deletes the points with curvature less than a certain threshold; however, this method is prone to data holes in flat areas, resulting in the loss of terrain detail information.

[0005] Based on this, the point reduction method and the point addition method are widely adopted. For example, prior art documents have proposed to perform point deletion operations on all point clouds by constructing a TIN iteratively. In each iteration, all point clouds are evaluated in sequence, and points with a value less than a specified tolerance are removed. This method has good performance, but its computational efficiency is low and it is not applicable to airborne LiDAR point cloud data with a huge amount of data. In the research on the point addition method, the classic maximum Z tolerance method selects the point with the largest deviation from the current TIN surface from the point cloud to be evaluated, and reconstructs the TIN again using the updated key point set. Although this method takes into account the global terrain features, it is prone to losing the shape and topological relationships at river network features, and the features at micro terrains are always ignored.

[0006] Subsequently, prior art documents have proposed a greedy multiquadric (MQ) method. This method first selects initial key points in the terrain for interpolation, and then iteratively compares the elevation differences between the interpolation surface generated by MQ and the corresponding point set to evaluate the importance of all candidate points. This method uses a non-linear interpolation method to improve the interpolation accuracy of the reference surface, but reduces the computational efficiency of selecting key points. On this basis, prior literature has also proposed to use the thin plate spline (TPS) interpolation method to select terrain key points. Compared with the greedy MQ method, this method not only improves the simplification efficiency, but also preserves the river network features in the terrain; however, this method has a high dependence on the reliability of the terrain feature extraction algorithm. In addition, prior art documents have also proposed a coarse-to-fine point cloud reduction method, which mainly generates a refined point cloud with data density varying with the terrain surface complexity by distinguishing local terrain complexity; however, this method ignores the feature points on the terrain skeleton line and is prone to terrain distortion problems.

[0007] In summary, the current point reduction method and point addition method used in airborne LiDAR ground point cloud simplification mainly achieve point cloud reduction by iteratively deleting redundant points in the original point cloud or retaining key points. However, this process highly depends on the accurate calculation of the interpolation method, and uses the elevation difference between points and the interpolation surface to identify key points, which is prone to losing the detailed features of local terrain. Summary of the Invention

[0008] The object of the present invention is to propose a LiDAR point cloud clustering and simplification method considering terrain features to solve the problems existing in the current ground point cloud simplification methods, such as inaccurate identification of terrain feature points and easy loss of terrain detailed features.

[0009] To achieve the above object, the present invention adopts the following technical solutions:

[0010] A LiDAR point cloud clustering and simplification method considering terrain features, comprising the following steps:

[0011] Step 1. Using the LiDAR ground point cloud as the original point cloud data, the original point cloud data is divided into multiple initial point cloud clusters by using the K-means algorithm; each initial point cloud cluster is further subdivided into multiple point cloud sub-clusters according to the terrain complexity of the point cloud cluster;

[0012] Step 2. First, according to the position information characteristics of the terrain feature lines in each point cloud sub-cluster, the terrain feature points in each point cloud sub-cluster are searched, and the found terrain feature points are used as the representative points of the point cloud sub-cluster;

[0013] If there are no terrain feature points in the point cloud sub-cluster, the centroid point of the point cloud sub-cluster is used as the representative point of the point cloud sub-cluster;

[0014] Step 3. Considering the problem of the boundary contraction of the original point cloud caused by the clustering simplification method, the boundary feature points in the original point cloud data are identified by using the two-dimensional Douglas-Peucker method;

[0015] The final ground simplified point cloud includes the representative points of each point cloud sub-cluster and the boundary feature points.

[0016] The present invention has the following advantages:

[0017] As shown above, the present invention describes a LiDAR point cloud clustering and simplification method considering terrain features. Aiming at the problems of inaccurate identification of terrain feature points and easy loss of terrain detail features in the current airborne LiDAR ground point cloud simplification method, an adaptive point cloud clustering considering terrain complexity is first adopted to generate point cloud sub-clusters; then the point cloud sub-clusters are described by identifying the normal vector mutation feature points, elevation mutation feature points and centroid points in the point cloud sub-clusters; finally, the boundary key points of the point cloud are captured to prevent the boundary contraction of the original point cloud. Compared with the traditional method, at the same point cloud simplification ratio, the accuracy of the digital elevation model (DEM) generated by the present invention and the accuracy of its derivatives (including average slope and terrain roughness) are significantly better than those of the traditional method, and the terrain feature information is better retained, providing an effective technical support for the reduction of remote sensing point cloud big data. Description of the Drawings

[0018] Figure 1 It is a flow chart of the LiDAR point cloud clustering and simplification method considering terrain features in the embodiment of the present invention.

[0019] Figure 2 It is a flow chart of point cloud clustering based on terrain complexity in the embodiment of the present invention.

[0020] Figure 3 It is a comparison chart of terrain feature lines and TPI in the embodiment of the present invention.

[0021] Figure 4 It is a schematic diagram of initial point cloud clusters assigned with TPI in the embodiment of the present invention.

[0022] Figure 5 Schematic diagram of point cloud clustering considering terrain complexity in the embodiment of the present invention.

[0023] Figure 6 Flow chart of terrain feature point recognition in the embodiment of the present invention.

[0024] Figure 7 Schematic diagram of terrain feature point recognition in the normal vector mutation area in the embodiment of the present invention.

[0025] Figure 8 Schematic diagram of terrain feature point recognition in the elevation mutation area in the embodiment of the present invention.

[0026] Figure 9 Schematic diagram of terrain boundary missing caused by clustering simplification in the embodiment of the present invention.

[0027] Figure 10 Schematic diagram of two-dimensionalization of boundary point cloud in the embodiment of the present invention.

[0028] Figure 11 Schematic diagram of simplified point cloud distribution in different terrain feature areas in the embodiment of the present invention.

[0029] Figure 12 Schematic diagram of reference DEMs of six groups of data.

[0030] Figure 13 Schematic diagram of comparison of RMSE and MAE of four methods in each region under different point cloud simplification ratios;

[0031] Among them, the left axis corresponds to the bar chart, and the right axis corresponds to the line chart.

[0032] Figure 14 DEM schematic diagram of four methods for Data2 at a simplification rate of 1%.

[0033] Figure 15 DEM schematic diagram of four methods for Data6 at a simplification rate of 1%.

[0034] Figure 16 Schematic diagram of terrain roughness and average slope of six methods in various research regions;

[0035] Among them, the left axis corresponds to the bar chart, and the right axis corresponds to the line chart. Detailed implementation manners

[0036] The present invention relates to a method for clustering and simplifying LiDAR point clouds considering terrain features. The method for clustering and simplifying LiDAR point clouds considering terrain features mainly includes the following steps: (1) adaptively clustering the point clouds using terrain complexity to reduce the loss of detailed features caused by insufficient point density at terrain features; (2) identifying key points located at terrain features within the point cloud clusters according to the different characteristics of terrain features such as ridges, valleys, and cliffs to reduce the loss of terrain features caused by inaccurate identification of key points; (3) avoiding the problem of boundary reduction caused by point cloud simplification by retaining the boundary feature points of the point cloud.

[0037] The present invention will be further described in detail below with reference to the accompanying drawings and specific embodiments:

[0038] As Figure 1 shown, a method for clustering and simplifying LiDAR point clouds considering terrain features includes the following steps:

[0039] Step 1. Using the LiDAR ground point cloud as the original point cloud data, dividing the original point cloud data into multiple initial point cloud clusters by using the K-means algorithm; further dividing each initial point cloud cluster into multiple point cloud sub-clusters according to the terrain complexity of the point cloud cluster.

[0040] Traditional point cloud simplification methods based on three-dimensional voxel grids divide the point cloud data into equally sized voxels and select a point within the voxel to replace all the point clouds within the voxel. Although this method has high computational efficiency, it is prone to accuracy loss in complex terrain feature areas. The present invention clusters the point clouds to be simplified using the similarity of adjacent point sets to generate point cloud clusters, and uses the cluster as the processing unit. Considering the large amount of airborne LiDAR ground point clouds, the present invention is based on the K-means method with high computational efficiency for its rapid division. However, this method only utilizes the spatial attributes of the point cloud and does not consider the terrain complexity. Therefore, the present invention proposes a method for clustering point clouds considering terrain complexity, and the specific process is as Figure 2 shown.

[0041] First, use the K-means method to divide the point cloud into multiple point cloud clusters.

[0042] Then, further divide each cluster into more sub-clusters based on the terrain complexity of the point cloud cluster. Among them, the higher the complexity of the cluster, the more the corresponding number of sub-clusters, and thus the more points retained after replacing the sub-clusters with key points.

[0043] Specifically, step 1 is as follows:

[0044] Step 1.1. Initial point cloud cluster division.

[0045] Divide the original point cloud data using the K-means algorithm, and determine the number k of initial point cloud clusters according to formula (1) first .

[0046] k first = [m * t scale (1)

[0047] In the formula, [·] is the ceiling function, m is the number of simplified point clouds, and t scale is the point cloud cluster subdivision ratio, and t scale ∈(0, 1); the smaller the value of t scale , the more sub - clusters there are in the terrain complex area.

[0048] Step 1.2. Terrain complexity calculation.

[0049] Compared with the massive and disordered scattered point clouds, it is faster to extract terrain features using a regular grid. Therefore, the present invention uses the terrain position index TPI to describe the terrain complexity, and TPI is shown in formula (2).

[0050]

[0051] In the formula, |·| is the absolute value symbol, Z0 represents the elevation value of the current grid, and Z i represents the elevation value of the i - th grid adjacent to the current grid, where i ∈ 1, 2, ···, 8.

[0052] As Figure 3 shown, the terrain position index TPI can better describe terrain feature information, such as ridge lines, valley lines, etc. Among them, Figure 3 (a) represents the contour line and terrain feature line, Figure 3 (b) represents the corresponding TPI value.

[0053] Step 1.3. Assigning terrain complexity to point cloud clusters.

[0054] In view of the accuracy loss caused by rasterizing the point cloud to calculate TPI, to prevent under - subdivision of point cloud clusters, the present invention assigns the maximum value in the TPI grid points where each initial point cloud cluster is located to the initial point cloud cluster, as Figure 4 shown.

[0055] Among them, the larger the TPI value, the more complex the terrain where the point cloud cluster is located.

[0056] Step 1.4. Point cloud cluster subdivision.

[0057] According to the proportion of the terrain complexity of each point cloud cluster, that is, the TPI value in all point cloud clusters, determine the further subdivision number of the initial point cloud cluster, and complete the division of sub - clusters with the help of K - means.

[0058] The further subdivision number C k of the k - th initial point cloud cluster is calculated by the following formula:

[0059]

[0060] In the formula, [·] is the ceiling symbol, and TPI k is the TPI index of the k-th cluster. The larger the value of this TPI index, the more sub-clusters the cluster is divided into, and the more corresponding cluster centroid points are in the complex terrain area, such as Figure 5 shown.

[0061] Figure 5 (a) shows the distribution of the initial point cloud clusters after clustering the original point cloud data using K-means.

[0062] Figure 5 (b) shows the distribution of the point cloud sub-clusters after the initial point cloud clusters are subdivided by adding terrain complexity.

[0063] Figure 5 (c) shows the centroid points of the point cloud sub-clusters.

[0064] Step 2. First, according to the characteristics of the position information of the terrain feature lines in each point cloud sub-cluster, find the terrain feature points in each point cloud sub-cluster, and use the found terrain feature points as the representative points of the point cloud sub-clusters.

[0065] If there are no terrain feature points in the point cloud sub-cluster, use the centroid point of the point cloud sub-cluster as the representative point of the point cloud sub-cluster.

[0066] Although directly using the cluster centroid as the representative point of the point cloud cluster is simple, when the centroid point does not fall on the terrain feature, it will cause the problem that the simplified points are difficult to accurately describe the terrain feature information.

[0067] The present invention adopts the strategy of using the terrain feature points in the point cloud sub-cluster to replace the centroid point to represent the point cloud sub-cluster. According to the actual terrain feature line types, the present invention proposes two different terrain feature point recognition schemes, such as Figure 6 shown.

[0068] Step 2.1. For terrain features such as ridges and valleys, the normal mutation terrain feature points are often located at the normal vector mutation points, such as Figure 7 shown. Therefore, by comparing the normal vector differences within each sub-cluster, it is judged whether there is a normal vector mutation feature within the sub-cluster.

[0069] The point cloud normal vector is calculated by the principal component analysis (PCA) method based on an adaptive neighborhood.

[0070] The specific steps for identifying the feature points in the normal vector mutation area are as follows:

[0071] First, calculate the normal information of all points in the point cloud sub-cluster, and then use the K-means method to divide the point cloud sub-cluster into two clusters according to the normal vector ( Figure 7The normal angle θ between the surfaces represented by the two clusters is calculated for the hollow triangles and hollow squares in the middle. C .

[0072] If θ C is less than or equal to the preset normal angle threshold θ th , the normal vectors within the point cloud sub-cluster are similar, and there is no sudden change in the normal vectors.

[0073] If θ C is greater than θ th , the adjacent points of the two clusters ( Figure 7 the solid squares in the middle) are marked as candidate feature points, and the point closest to the centroid of the point cloud sub-cluster among all candidate feature points is used as the terrain feature point with sudden change in normal vector for the point cloud sub-cluster.

[0074] Step 2.2. At terrains with drastic elevation changes (such as steep cliffs), the elevation between adjacent point clouds is discontinuous.

[0075] The K-means method divides the point cloud using the three-dimensional coordinates of the point cloud, and adjacent points with obvious elevation differences will be divided into different clusters. Therefore, the terrain feature points in the elevation mutation area should be located in the edge area of the point cloud sub-cluster.

[0076] The fracture terrain feature points located in the elevation mutation area can be identified by using the elevation differences between the adjacent edge points of adjacent point cloud sub-clusters.

[0077] The specific steps for identifying terrain feature points in the elevation mutation area are as follows:

[0078] I. Using the centroid of the point cloud sub-cluster as a node, the adjacency relationship of all point cloud sub-clusters is described using a graph structure, as shown in Figure 8 (a).

[0079] II. Using the α-shape method to find the edge point set of each point cloud sub-cluster, as shown in Figure 8 (b).

[0080] III. According to the proximity relationship between point cloud sub-clusters, search for the closest point between the edge points of the point cloud sub-cluster and the edge point set of its neighboring sub-cluster point by point. Finally, each edge point P of the point cloud sub-cluster cb ( Figure 8 the hexagonal points in (c)) is matched with an edge point P of a neighboring sub-cluster nb ( Figure 8 the square points in (c)).

[0081] IV. Calculate the elevation difference d nb between the corresponding two points of P cb and P i .

[0082] As shown in Figure 8 (d), if d iare all less than the preset elevation difference threshold f b , then there is no elevation mutation between the point cloud sub-cluster and its neighboring sub-cluster. Otherwise, take the P i corresponding to the maximum value of d cb point as the fracture terrain feature point of this sub-cluster.

[0083] Step 2.3. If there are no normal mutation terrain feature points in Step 2.1 and fracture terrain feature points in Step 2.2 for the point cloud sub-cluster, it indicates that the area where the point cloud sub-cluster is located is relatively flat. The closer its representative point is to the centroid, the stronger its representation ability for the local terrain surface. Then use the centroid point of the current point cloud sub-cluster as the representative point of the point cloud sub-cluster.

[0084] Step 3. The feature points selected based on the clustering method are likely to cause the loss of boundary point clouds, resulting in the contraction of the original point cloud boundary, that is, the simplified point coverage area is smaller than the area of the original point cloud, as Figure 9 shown. Therefore, the present invention thins the boundary point clouds based on the two-dimensional Douglas-Peucker method and retains the key points therein to prevent the contraction of the original point cloud boundary. The specific steps are as follows:

[0085] Step 3.1. Extract the boundary point set {P j |j = 0, 1, …, n} of the original point cloud data.

[0086] Among them, P j represents the j-th boundary point, and n + 1 represents the total number of boundary point clouds.

[0087] Select any one of them as the origin point, sort all the boundary point sets in the counterclockwise direction according to their adjacent relationships, and convert the three-dimensional boundary point set into two-dimensional coordinates according to formula (4), as Figure 10 shown.

[0088]

[0089] Among them, (x j , y j , z j ) represents the three-dimensional coordinates of the j-th point in the three-dimensional boundary point set; (x j+1 , y j+1 , z j+1 ) represents the three-dimensional coordinates of the (j + 1)-th point in the three-dimensional boundary point set sorted according to its adjacent relationship.

[0090] X j is the position of the j-th boundary point on the abscissa in the two-dimensional coordinates, and (X j+1 , Y j+1 ) represents the position of the (j + 1)-th boundary point in the two-dimensional coordinates; (X1, Y1) represents the position of the first boundary point in the two-dimensional coordinates, and (X1, Y1) = (0, 0).

[0091] Step 3.2. Mark the boundary point set of the two-dimensional coordinates as {P j '|j = 0, 1, …, n}.

[0092] Take the boundary points P0' and P end ' in the boundary point set of the two-dimensional coordinates as the baseline points, and the line connecting the baseline points P0' and P end ' as the baseline, calculate the perpendicular distance from all boundary points P i ' to the baseline, and find the maximum value d max .

[0093] Among them, the two-dimensional coordinates of P0' are (X1, Y1); the two-dimensional coordinates of P end ' are (X j+1 , Y j+1 ).

[0094] If the maximum value d max is less than the preset height difference d th , then delete all the intermediate boundary point clouds and only retain the two endpoints of the baseline; if the maximum value d max is greater than or equal to the preset height difference d th , mark the boundary point corresponding to d max as the new baseline point.

[0095] Step 3.3. Use the new baseline points and the baseline points P0' and P end ' to divide the boundary point set of the two-dimensional coordinates into two groups of boundary point sets, and repeat Step 3.2 for the newly divided boundary point sets to determine whether there are new baseline points in the boundary point sets;

[0096] If there are new baseline points, continue to divide the boundary point sets, otherwise retain the baseline points of the boundary point sets;

[0097] Finally, until the maximum value d max in all boundary point sets is less than the preset height difference d th , that is, all boundary point sets do not need to be further divided, then retain the three-dimensional boundary point clouds corresponding to all baseline points;

[0098] Take the retained three-dimensional boundary point clouds as the boundary feature points of the original point cloud data.

[0099] Finally, the ground simplified point cloud includes the representative points of each point cloud sub-cluster and the boundary feature points.

[0100] Figure 11 Among them, ridge lines, valley lines ( Figure 11 (a)), elevation mutation points ( Figure 11 (b)) all have feature points, and the feature point clouds of the study area boundary are also well retained, indicating that the simplified points have good terrain feature retention ability.

[0101] In addition, in order to verify the effectiveness of the LiDAR point cloud clustering and simplification method considering terrain features of the present invention, experiments were also carried out. In the experiments, six groups of high-density airborne LiDAR point cloud instance data were used as the research object, including different landform types such as plains, mountains and hills, and there are terrain features such as multiple ridge lines, valley lines and steep cliffs, as Figure 12 shown.

[0102] Among them, the overall terrain undulation of Data1 and Data 2 is relatively small. The overall area of Data1 is flat, and there is an obvious hillock in the upper left position. The area of Data2 has obvious river network features; the overall terrain undulation of Data3 and Data 4 is relatively large, and the ridge line and valley line features are obvious; the terrain in the areas of Data5 and Data 6 is discontinuous, and there is a large elevation fracture terrain (such as terrain features such as steep cliffs). Since the acquired original point cloud contains ground points and non-ground points, first, semi-automatic filtering and visual discrimination correction are carried out on it through commercial software TerraSolid, and then the ground points are accurately extracted.

[0103] Among them, the statistical information of the six groups of ground point data is shown in Table 1.

[0104] Table 1 Terrain statistical information of six groups of data

[0105]

[0106] To verify the effectiveness and applicability of the method of the present invention, the method of the present invention is compared and analyzed with the random method (Random), voxel grid method (VG), curvature sampling method (Curvature) in the open-source PCL (point cloud library) and the maximum Z tolerance method (Max-Z) in ArcGIS. The above five methods are all tested by retaining the same number of simplified points, and the selected key point numbers are set to 0.1%, 0.2%, 0.4%, 0.6%, 0.8% and 1% of the total number of ground points respectively. In order to quantitatively evaluate the performance of the simplification method, the root mean square error (RMSE) and mean absolute error (MAE) are used to compare the differences between the simplified point cloud and the original point cloud to generate DEM. At the same time, the average slope and terrain roughness (K) are used to evaluate the influence of the point cloud reduction method on the loss of terrain surface morphology. The calculation formulas of each accuracy index are as follows:

[0107]

[0108]

[0109]

[0110]

[0111] wherein, H i and E i represent the i-th DEM grid point generated by the original LiDAR points and the simplified point cloud, N is the number of DEM grids, and S i is the slope of each grid, A i is the surface area of each grid, and A p is the projected area.

[0112] All DEM resolutions generated by the present invention are 1 m. In order to test the best performance of the algorithm, the parameters of the present invention are set according to the following criteria: the grid size l is based on being able to extract the main terrain features, and the fracture height difference threshold f b should be greater than the elevation difference of the minimum terrain fracture, and the normal deviation threshold θ th is set to be greater than the normal mutation size of the regional terrain features, the boundary height difference threshold d th is set according to the undulation degree of the regional boundary, and the point cloud cluster subdivision ratio t scale is 0.5.

[0113] Based on this, the parameter set adopted by the method of the present invention for processing 6 groups of data is:

[0114] l = 4 m, f b = 10 m, θ th = 30°, d th = 4 m, t scale = 0.5.

[0115] Such as Figure 13 (a) to Figure 13As shown in Fig. (f), the present invention compares the RMSE and MAE of five point cloud reduction methods in each region under different point cloud simplification ratios. The results show that regardless of the method, the more simplified points are retained, the higher the accuracy of the generated DEM, that is, the RMSE and MAE decrease with the increase of the simplification rate. Among the five methods, Random, Curvature, and Max-Z have lower accuracy. Due to the uniform distribution of the simplified point cloud, VG performs relatively better. In contrast, the method of the present invention has the best accuracy and is the most robust, especially in regions Data1-4, where its RMSE and MAE are significantly lower than those of other methods. This is mainly due to the fact that this method adopts point cloud clustering considering terrain complexity and the identification of feature points in regions with normal mutations, which can retain more simplified points in complex terrain areas while identifying terrain feature points within the regions. In the cliff regions (Data5, 6), since the method of the present invention can identify feature points at fractured terrains, its RMSE and MAE are significantly lower than those of the Random, Curvature, and Max-Z methods. Table 2 lists the average errors of each method under different simplification ratios. The results show that the average RMSE and MAE of the method of the present invention in each region are significantly lower than those of the traditional methods. Among them, the average RMSE of the method of the present invention for the six groups of data is reduced by 37.8%, 40.3%, 12.1%, and 51.8% respectively compared with Random, Curvature, VG, and Max-Z, and the average MAE is reduced by 28.5%, 39.8%, 9.6%, and 52.2% respectively.

[0116] Table 2 Average RMSE (m) and MAE (m) of the five methods for each experimental area under different simplification ratios

[0117]

[0118]

[0119] Compared with the original terrain surface, the simplified ground surface will be distorted with the reduction of the simplified point cloud. Therefore, under the same simplification ratio, the method that can highlight and retain terrain features has better simplification performance. Therefore, taking Data2 and Data6 as examples, the shaded relief maps of the DEMs generated by the five methods at a 1% point cloud reduction rate are compared, as shown in Figure 14 (a) to Figure 14 (f), 15(a) to Figure 15As shown in (f). Among the five methods, Curvature has excessive smoothing in flat areas caused by data holes. This is because refined points tend to gather in areas with larger curvatures, and there are discontinuities in the terrain feature lines at the river valleys of Data2, and serious distortion of features at the fractured terrain of Data6. The Random method shows serious blurring of terrain features and local smoothing in all regions. The VG method performs slightly better than the former overall, but due to its inability to capture terrain features, some tributaries on the left side and mountain roads on the right side in Data2 are lost, and there are obvious sawteeth at the steep cliffs in the middle of Data6. The Max-Z method is prone to pothole filling problems and has serious smoothing in non-feature regions, which does not match the original terrain features. Overall, the method of the present invention has the best performance, and it better retains the detailed features of the river network in Data2 and the terrain feature information at the steep cliffs in Data6.

[0120] Figure 16 (a) to Figure 16 (a) to (f) show the average slope and terrain roughness of the five methods at different simplification ratios. The results show that there are significant differences between the average slope and terrain roughness of the Random, Curvature, and VG methods and the reference values. The average slope and terrain roughness of the method of the present invention and the Max-Z method are higher than those of other methods. However, the former takes into account the detailed features in the terrain while retaining the terrain skeleton features, so it is superior to the latter. At the same simplification rate, the method of the present invention can always maintain a relatively high average slope and terrain roughness. Overall, the method of the present invention is more consistent with the original terrain and is significantly superior to other methods. As shown in Table 3, the average slope of the method of the present invention is increased by 2.3%, 1.7%, and 2.6% respectively compared with Radom, VG, and Curvature, and the average terrain roughness is increased by 1.2%, 0.9%, and 1.0% respectively, indicating that it has a strong ability to maintain terrain features. Since the Max-Z method retains more redundant points in complex terrains, increasing the undulation degree of the local terrain, the average slope and average terrain roughness in the Data1, 3, and 6 regions are closer to the true values.

[0121] Table 3 Terrain roughness and average slope of the five methods for six regions at different simplification ratios

[0122]

[0123]

[0124] In summary, six groups of high-density airborne LiDAR ground point clouds were selected as the research objects, and the method of the present invention was compared with the random method, voxel grid method, Curvature method, and maximum Z tolerance method. According to the analysis of the data of the six groups of examples, it shows that the calculation accuracy of the method of the present invention is significantly better than that of the Random, VG, Curvature, and Max-Z methods. The average RMSE of the DEM constructed by it is at least 12.1% lower than that of other methods, and the average MAE is at least 9.6% lower. Compared with the traditional point cloud simplification method, the average terrain roughness and average slope of the method of the present invention are more consistent with the reference values. The former is increased by 0.9%-1.2%, and the latter is increased by 1.7%-2.6%, indicating that the method of the present invention can better preserve the terrain feature information.

[0125] Of course, the above description is only the preferred embodiment of the present invention. The present invention is not limited to listing the above embodiments. It should be noted that all equivalent substitutions and obvious deformation forms made by any person skilled in the art under the teaching of this specification fall within the substantial scope of this specification and should be protected by the present invention.

Claims

1. A LiDAR point cloud clustering and simplification method considering terrain features, characterized in that It includes the following steps: Step 1. Using the LiDAR ground point cloud as the original point cloud data, divide the original point cloud data into multiple initial point cloud clusters by using the K-means algorithm; further divide each initial point cloud cluster into multiple point cloud sub-clusters according to the terrain complexity of the point cloud cluster; Step 2. First, according to the position information characteristics of the terrain feature lines in each point cloud sub-cluster, find the terrain feature points in each point cloud sub-cluster, and use the found terrain feature points as the representative points of the point cloud sub-clusters; If there are no terrain feature points in the point cloud sub-cluster, use the centroid point of the point cloud sub-cluster as the representative point of the point cloud sub-cluster; Step 3. Considering the problem of the boundary contraction of the original point cloud caused by the clustering simplification method, use the two-dimensional Douglas-Peucker method to identify the boundary feature points in the original point cloud data; The final ground simplified point cloud includes the representative points of each point cloud sub-cluster and the boundary feature points; The specific content of the said Step 1 is as follows: Step 1.

1. Initial point cloud cluster division; Divide the original point cloud data using the K-means algorithm and determine the number of initial point cloud clusters k according to formula (1). first ; k first = [m * t scale (1) where [·] is the ceiling symbol, m is the number of simplified point clouds, and t scale is the point cloud cluster subdivision ratio, t scale ∈(0, 1); Step 1.

2. Terrain complexity calculation; Use the terrain position index TPI to describe the terrain complexity, and TPI is shown in formula (2); where |·| is the absolute value symbol, Z0 represents the elevation value of the current grid cell, and Z i represents the elevation value of the i-th grid cell adjacent to the current grid cell, where i ∈ 1, 2, ···, 8; Step 1.

3. Assign the terrain complexity to the point cloud cluster; Assign the maximum value in the TPI grid points where each initial point cloud cluster is located to the initial point cloud cluster; Step 1.

4. Point cloud cluster subdivision; According to the proportion of the terrain complexity (i.e., the TPI value) of each point cloud cluster in all point cloud clusters, determine the number of further subdivisions of the initial point cloud cluster, and use the K-means to complete the division of the point cloud cluster; The further subdivision number C of the k-th initial point cloud cluster k The calculation formula is as follows: where, TPI k is the TPI index of the k-th cluster, and [·] is the ceiling symbol.

2. The LiDAR point cloud clustering and simplification method considering terrain features according to claim 1, characterized in that The specific content of the said Step 2 is as follows: Step 2.

1. For the terrain features of ridges and valleys, by comparing the differences in normal vectors within each point cloud sub-cluster, judge whether there are normal vector mutation features within the point cloud sub-cluster, and identify the normal vector mutation terrain feature points located at ridges and valleys; Step 2.

2. In the terrain where the elevation changes violently, the elevation between adjacent point clouds is discontinuous. Use the elevation difference between the adjacent edge points of adjacent point cloud sub-clusters to identify the fracture terrain feature points located in the elevation mutation area; Step 2.

3. If the point cloud sub-cluster does not have the normal vector mutation terrain feature points in Step 2.1 and the fracture terrain feature points in Step 2.2, use the centroid point of the current point cloud sub-cluster as the representative point of the point cloud sub-cluster.

3. The LiDAR point cloud clustering and simplification method considering terrain features according to claim 2, characterized in that In the said Step 2.1, the specific steps for identifying the terrain feature points in the normal vector mutation area are as follows: First, calculate the normal information of all points within the point cloud sub-cluster. Then, use the K-means method to divide the point cloud sub-cluster into two clusters based on the normal vectors, and calculate the normal angle θ of the surfaces represented by the two clusters. C ; If θ C is less than or equal to the preset normal angle threshold θ th , the normal vectors within the point cloud sub-cluster are similar and there is no sudden change in the normal vector; If θ C is greater than the preset normal angle threshold θ th , then mark the adjacent points of the two clusters as candidate feature points, and take the point closest to the centroid of the point cloud sub-cluster among all candidate feature points as the normal mutation terrain feature point of the point cloud sub-cluster.

4. The LiDAR point cloud clustering and simplification method considering terrain features according to claim 2, characterized in that In the said Step 2.2, the specific steps for identifying the terrain feature points in the elevation mutation area are as follows: I. Using the centroid of the point cloud sub-cluster as a node, use the graph structure to describe the adjacency relationship of all point cloud sub-clusters; II. Use the α-shape method to find the edge point set of each point cloud sub-cluster; III. Search for the nearest point between the edge points of the point cloud sub-cluster and the set of edge points of its neighboring point cloud sub-clusters point by point according to the proximity relationship between the point cloud sub-clusters. Finally, each edge point P of the point cloud sub-cluster cb is matched with an edge point P of a neighboring point cloud sub-cluster nb ; IV. Calculate P nb and P cb The elevation difference d between two points i ; If the elevation difference d i is less than the preset elevation difference threshold f b , then there is no elevation mutation between the point cloud sub-cluster and its neighboring point cloud sub-clusters; otherwise, the edge point P corresponding to the maximum value of d i is used as the fracture terrain feature point of this point cloud sub-cluster. cb ​ 5. The LiDAR point cloud clustering and simplification method considering terrain features according to claim 1, characterized in that The specific content of the said Step 3 is as follows: Step 3.

1. Extract the boundary point set {P j | j = 1, …, n}; where P j represents the j-th boundary point, and n represents the total number of boundary point clouds; Select any one of the boundary points from the boundary point set as the origin, sort all the boundary point sets in a counterclockwise direction according to their adjacent relationships, and convert the three-dimensional boundary point set into two-dimensional coordinates according to formula (4); where (x j , y j , z j ) represents the three-dimensional coordinates of the j-th boundary point in the three-dimensional boundary point set; (x j+1 , y j+1 , z j+1 ) represents the three-dimensional coordinates of the (j + 1)-th boundary point sorted by the adjacent relationship in the three-dimensional boundary point set; X j is the abscissa position of the j-th boundary point in the two-dimensional coordinate. (X j+1 , Y j+1 ) represents the position of the (j + 1)-th boundary point in the two-dimensional coordinate; (X1, Y1) represents the position of the 1st boundary point in the two-dimensional coordinate, and (X1, Y1) = (0, 0); Step 3.

2. Mark the boundary point set of the two-dimensional coordinates as {P j '|j = 1, …, n}; Take the boundary points P0' and P end ' in the two-dimensional coordinate boundary point set as the baseline points. The line connecting the baseline points P0' and P end ' is used as the baseline. Calculate the perpendicular distance from all boundary points P i ' to the baseline, and find the maximum value d among the perpendicular distances max ; Among them, the two-dimensional coordinates of P0' are (X1, Y1); P end 's two-dimensional coordinates are (X j+1 , Y j+1 ); If the maximum value d max is less than the preset height difference d th , then all the intermediate boundary point clouds are deleted, and only the two endpoints of the baseline are retained; if the maximum value d max is greater than or equal to the preset height difference d th , then d max the corresponding boundary points are marked as new baseline points. Step 3.

3. Using the new baseline points and the baseline points P0' and P end ', divide the boundary point set of the two-dimensional coordinates into two groups of boundary point sets, and repeat Step 3.2 for the newly divided boundary point sets to determine whether there are new baseline points in the boundary point sets; If there are new baseline points, continue to divide the boundary point set; otherwise, retain the baseline points of the boundary point set; Finally, until the maximum value d within all boundary point sets max is less than the preset height difference d th , that is, all boundary point sets do not need to be further divided, then retain the three-dimensional boundary point cloud corresponding to all baseline points; Take the retained three-dimensional boundary point cloud as the boundary feature points of the original point cloud data.

Images