Point cloud hierarchical filtering method for complex terrain areas supplemented by multi-feature clustering

Abstract

The present invention discloses a method for hierarchical filtering of point clouds in complex terrain areas assisted by multi-feature clustering. The method comprises the following steps: Step 1. Gridding the original point cloud data, and sequentially performing multi-feature clustering on the point clouds in each independent grid to improve the effective division of adjacent objects / ground in complex environments, and mapping the obtained clustering results to three-dimensional space to obtain the three-dimensional coordinate distribution of the point cloud clusters; Step 2. Using the criterion of taking fractured terrain into consideration, the ground clusters are identified from the three-dimensional coordinate distribution of the point cloud clusters, thereby avoiding the loss of terrain details and improving the recognition rate of the ground point clusters; Step 3. Using the captured ground points in the ground clusters as ground seed points, an initial ground reference surface is constructed, and the ground points in the original point cloud are further refined in combination with multi-scale hierarchical filtering. The present invention has high filtering accuracy in complex terrain areas such as terrain fractures, and effectively distinguishes between object points and ground points under complex terrain.

Description

Technical Field

[0001] The invention relates to a point cloud hierarchical filtering method for complex terrain areas assisted by multi-feature clustering. Background Art

[0002] As an active remote sensing earth observation method, airborne laser radar (Light Detection and Ranging, LiDAR) technology plays a supporting role in the acquisition of large-scale real-life three-dimensional geographic information.

[0003] With the rapid development of airborne LiDAR systems, LiDAR point clouds have been widely used in many fields such as real-life 3D construction, urban intelligent management, and ecological and environmental protection. However, the original LiDAR point cloud is mixed with ground information and non-ground information such as vegetation, buildings, and vehicles. Therefore, point cloud filtering is the premise and key step for applications in the above fields.

[0004] In the past two decades, scholars at home and abroad have proposed a variety of filtering algorithms. In summary, they can be classified into five categories: morphology-based filtering, interpolation-based filtering, slope-based filtering, block-based filtering, and machine learning-based filtering. Among them, compared with other filtering algorithms, interpolation-based filtering has better overall applicability and can effectively filter out non-ground points.

[0005] However, due to the high complexity, interlacing, and polymorphism of terrain and objects in real scenes, as well as the similarity in elevation between some buildings located on hills and terrain, traditional interpolation filtering algorithms are still unable to effectively distinguish between ground points and object points in complex terrain areas, which greatly affects the filtering quality. In addition, traditional interpolation filtering algorithms only use the three-dimensional spatial coordinate information of point clouds, and lack full utilization of physical attribute information such as echo intensity, color, and echo number. Therefore, traditional interpolation filtering algorithms have low filtering accuracy in complex environments such as terrain fractures, and it is difficult to effectively distinguish between object points and ground points in complex terrain. Summary of the invention

[0006] The purpose of the present invention is to propose a point cloud hierarchical filtering method for complex terrain areas assisted by multi-feature clustering, so as to solve the technical problem of low filtering accuracy of existing filtering methods in complex terrain areas such as terrain breaks.

[0007] In order to achieve the above object, the present invention adopts the following technical scheme:

[0008] The point cloud hierarchical filtering method for complex terrain areas assisted by multi-feature clustering includes the following steps:

[0009] Step 1. Grid the original point cloud data, perform multi-feature clustering on the point cloud in each independent grid in turn, obtain the multi-feature clustering result, map the clustering result to three-dimensional space, and obtain the three-dimensional coordinate distribution of the point cloud cluster;

[0010] Step 2. Identify ground clusters from the three-dimensional coordinate distribution of point cloud clusters using a criterion that takes into account fracture terrain;

[0011] Step 3. Use the ground points in the ground cluster captured in step 2 as ground seed points to construct the initial ground reference surface, and further refine the ground points in the original point cloud by combining the multi-scale hierarchical filtering method.

[0012] The present invention has the following advantages:

[0013] As described above, the present invention describes a method for hierarchical filtering of point clouds in complex terrain areas assisted by multi-feature clustering. First, multi-feature clustering of point cloud geometry and physical information is combined to improve the effective division of adjacent objects / ground in complex environments; then, ground point clusters are identified using criteria that take into account fractured terrain to avoid loss of terrain details and improve the recognition rate of ground point clusters; finally, the height difference from the point to the corresponding interpolation point is used instead of the height difference from the point to the ground reference surface grid, overcoming the influence of spatial position errors of point clouds on steep slopes. The present invention has high filtering accuracy in complex terrain areas such as terrain fractures, and effectively distinguishes between object points and ground points under complex terrain. BRIEF DESCRIPTION OF THE DRAWINGS

[0014] Figure 1 It is a flow chart of a point cloud hierarchical filtering method for complex terrain areas assisted by multi-feature clustering in an embodiment of the present invention;

[0015] Figure 2 It is a schematic diagram of the interconnected area between the ground and the objects;

[0016] Figure 3 It is the relationship diagram between the sampling point and its corresponding 3×3 neighboring grid;

[0017] Figure 4 This is the result diagram of multi-feature clustering;

[0018] Figure 5 Schematic diagram of the three-dimensional coordinate distribution of the point cloud cluster;

[0019] Figure 6 Schematic diagram of ground cluster identification taking into account fractured terrain;

[0020] Figure 7 It is a schematic diagram of the traditional MHF;

[0021] Figure 8 Schematic diagram of MHF taking into account spatial position error in an embodiment of the present invention;

[0022] Fig. 9 This is an experimental data image of the gully terrain;

[0023] Fig.10 This is the experimental data image of the fracture terrain;

[0024] Fig.11 This is an image of experimental data of steep slope terrain;

[0025] Fig.12 This is the experimental data image of the fault steep slope terrain;

[0026] Fig.13 Schematic diagram of digital elevation model (DEM) generated based on reference ground point data of gully terrain;

[0027] Fig.14 Schematic diagram of digital elevation model (DEM) generated based on reference ground point data of fractured terrain;

[0028] Fig.15 Schematic diagram of a digital elevation model (DEM) generated based on reference ground point data of steep slope terrain;

[0029] Fig.16 Schematic diagram of a digital elevation model (DEM) generated based on reference ground point data of the fault scarp terrain;

[0030] Fig.17 Schematic diagram of the comparison of the total error of the method of the present invention and six representative algorithms;

[0031] Fig.18 A schematic diagram showing the comparison of the Kappa coefficients of the method of the present invention and six representative algorithms;

[0032] Fig.19 Schematic diagram of the comparison of mean error of DEM constructed by the method of the present invention and six representative algorithms;

[0033] Fig. 20 A schematic diagram showing the comparison of the average errors of DEMs constructed by the method of the present invention and six representative algorithms;

[0034] Fig.21 This is a schematic diagram of the reference DEM;

[0035] Fig. 22 Schematic diagram of DEMs after Plot4 filtering by the Proposed algorithm;

[0036] Fig.23 Schematic diagram of DEMs after Plot4 filtering by MHF algorithm;

[0037] Fig.24 Schematic diagram of DEMs after Plot4 filtering by SMRF algorithm;

[0038] Fig.25 Schematic diagram of DEMs after Plot4 filtering by CSF algorithm;

[0039] Fig.26 Schematic diagram of DEMs after Plot4 filtering by MSF algorithm;

[0040] Fig. 27 Schematic diagram of DEMs after Plot4 filtering by SBF algorithm;

[0041] Fig.28 Schematic diagram of DEMs after Plot4 filtering using the PTD algorithm. DETAILED DESCRIPTION

[0042] In view of the high complexity, overlap and polymorphism of objects in real scenes, in order to make full use of the geometric and physical information of point clouds, the present invention constructs a hierarchical point cloud filtering method assisted by multi-feature clustering.

[0043] Among them, point cloud clustering can better preserve the boundaries of objects / ground, thereby improving the ability to cope with complex environments.

[0044] like Figure 1 As shown in FIG, the hierarchical point cloud filtering method assisted by multi-feature clustering mainly includes four steps: data preprocessing, multi-feature point cloud clustering, ground cluster identification considering fracture terrain, and multi-scale hierarchical filtering (MHF).

[0045] The implementation process of the hierarchical point cloud filtering method assisted by multi-feature clustering is roughly as follows:

[0046] First, non-last echo points are eliminated based on echo information; then the point cloud is gridded, and multi-feature clustering is performed on the point cloud in each grid; then a criterion that takes into account fractured terrain is used to identify ground clusters, and finally MHF is used to further identify ground points.

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

[0048] Step 1. Data preprocessing and multi-feature point cloud clustering.

[0049] When a laser pulse propagates from the air to the ground, if it encounters multiple reflecting surfaces, multiple echoes will be generated.

[0050] Usually, flat surfaces such as the ground and building roofs only produce a single echo, while vertical structures such as vegetation and building edges will produce at least two echoes. Generally, except for the last echo, the probability that the other echoes in multiple echoes are ground points is low.

[0051] Before point cloud filtering, this embodiment first removes point clouds of non-last echoes such as building edges from the original point cloud according to echo information to reduce the amount of model calculation, thereby improving calculation efficiency.

[0052] Traditional interpolation filtering algorithms perform well in relatively simple scenarios, but they suffer from problems such as low accuracy and loss of terrain details in complex environments, which pose severe challenges to high-precision point cloud filtering.

[0053] Point cloud clustering can cluster connected terrain / ground object points into the same cluster according to specific discrimination rules, which can effectively make up for the loss of terrain details in interpolation filtering algorithms in complex scenes such as terrain fractures.

[0054] Among traditional clustering methods, DBSCAN clustering is a density-based spatial clustering method that uses the density connectivity of clusters to cluster point clouds into clusters of arbitrary shapes. It has the advantages of not requiring the number of clusters to be specified in advance and being able to find outliers while clustering. Based on this, the present invention uses the DBSCAN clustering method to cluster point clouds.

[0055] First, the original point cloud is gridded.

[0056] Theoretically, when processing high-density, large-scale point cloud data, DBSCAN has problems such as low computational efficiency, poor robustness, and poor clustering effect, and cannot meet the aggregation needs of homogeneous points in a large range of point clouds.

[0057] Therefore, in this embodiment, the original point cloud data is first gridded (the grid size is set to Ws), and then the point cloud in each independent grid is clustered with multiple features in turn. The advantages of this method are:

[0058] This can avoid the clustering confusion problem caused by similar features of points of different categories in a large-scale scene, and thus achieve efficient aggregation of point clouds with similar features; moreover, since only neighborhood points are operated, the clustering efficiency can be greatly improved.

[0059] Multi-feature point cloud clustering.

[0060] Usually, the point density of airborne LiDAR point clouds does not change much, and buildings and the ground are connected by tree point clouds, such as Figure 2 As shown, due to the influence of terrain undulation, two points with the same elevation may correspond to the ground and feature points respectively, and the elevations of the same type of feature points at different locations are also different.

[0061] Therefore, clustering point clouds in complex environments using only three-dimensional coordinate information can easily lead to erroneous clustering results.

[0062] Based on this, the present invention proposes to perform multi-feature point cloud clustering by combining point cloud geometry and physical information.

[0063] First, the ground reference surface is fitted using the least squares method, and the fitting residual of each point is calculated as the auxiliary feature of the point cloud.

[0064] In order to accurately describe the terrain surface, the lowest point in each grid is taken as the initial ground point, and the initial ground points of 9 adjacent grids are used to fit the tangent plane to calculate the fitting residual of each point in the grid.

[0065] like Figure 3 As shown, a sampling point is taken from the DTM grid. The point P i Located in the (I, J) grid of the DTM surface.

[0066] Among them, (I, J) represents a grid with the horizontal coordinate I and the vertical coordinate J.

[0067] Then its 9 nearest neighbor grids are (I-1, J-1), (I-1, J), (I-1, J+1), (I, J-1), (I, J), (I, J+1), (I+1, J-1), (I+1, J) and (I+1, J+1).

[0068] Set point P i The corresponding tangent plane equation is z = Ax + By + C, then P i The fitted residual of can be estimated as:

[0069] r i =z i -(Ax i +By i +C) (1)

[0070] Among them, (x i ,y i , z i ) is point P i The three-dimensional coordinates of , A, B, and C are the coefficients of the plane equation.

[0071] Since dense low vegetation has similar elevation and fitting residual information as the ground, low vegetation is easily mixed with ground points into a cluster during point cloud clustering. Given that the same ground object often has similar or continuous RGB color information, and the color variation between different ground objects is large, the present invention also uses color information as an auxiliary feature for point cloud clustering.

[0072] In order to simplify the calculation and highlight the difference between the object and the ground, the RGB color information is converted into a single-value visible-band difference vegetation index (VDVI).

[0073] The calculation formula of VDVI is:

[0074] Where R, G, and B represent the red, green, and blue color values ​​of the point cloud, respectively.

[0075] This embodiment uses the matrix [elevation, fitting residual, VDVI] as the point cloud virtual coordinates (z, h, VDVI) to perform DBSCAN clustering based on Euclidean distance, where the elevation is the z coordinate in the spatial geometric information.

[0076] The specific steps of DBSCAN clustering in this embodiment are as follows:

[0077] Step 1.3.1. Mark all points in the original point cloud data as unclassified points, and calculate the Euclidean distance between any two points in the unclassified points: The Euclidean distance calculation formula between any two points in the clustering process is as follows:

[0078]

[0079] Among them, d ij Represents the Euclidean distance between point i and point j.

[0080] z i 、z j Represents the elevation of point i and point j, r i 、r j Represent the fitting residual values ​​of point i and point j respectively, VDVI i 、VDVI j They represent the VDVI values ​​at point i and point j respectively. α, β, and γ are the proportional factor coefficients. α, β, and γ are all greater than 0.

[0081] Step 1.3.2. Randomly select a point P from the unclassified points and mark it as a classified point. Then search the ε neighborhood N of the point and calculate the number of points in N. If the number of points in the set N is greater than the minimum number of points P, ts , then mark P as a core point, create a new cluster C, and add point P to cluster C.

[0082] Among them, ε represents the radius of the search neighborhood point; N is the set of points in the neighborhood; P ts Indicates the minimum number of points.

[0083] Step 1.3.3. Traverse the ε neighborhood M of all unclassified points q in N. If the number of points in M ​​is greater than P ts , then add the point in M ​​to N, and if point q does not belong to any cluster, then add q to cluster C.

[0084] Step 1.3.4. Repeat step 1.3.3 until no new points are added to cluster C, then output cluster C and return to step 1.3.2; if no new cluster is output, the remaining points are classified as noise points until there are no points marked as unclassified, then go to step 1.3.5.

[0085] Step 1.3.5. Output result: The clusters are divided into C′={C1,C2,C3,C4,……,C k}.

[0086] Among them, C′ represents the set of point cloud clusters, C1, C2, C3, C4, ..., C k They represent the divided point cloud clusters respectively.

[0087] The results after multi-feature clustering are as follows Figure 4 As shown by Figure 4 The clustering results show that the multi-feature clustering fully integrates multiple features such as elevation, fitting residual, VDVI, etc., and can effectively aggregate point cloud clusters with similar features.

[0088] However, it is impossible to determine the category of clusters (ground or non-ground clusters) by relying solely on multi-feature clustering results. Therefore, it is necessary to map the clustering results into a three-dimensional space such as Figure 5 As shown, the three-dimensional coordinate distribution of the point cloud cluster is obtained.

[0089] Then, based on the three-dimensional coordinate distribution of the point cloud cluster and its geometric characteristics, the category of the point cloud cluster is determined.

[0090] Step 2. Ground cluster identification taking into account fractured terrain.

[0091] Traditional methods generally judge the category of point cloud clusters based on the curvature characteristics and whether they contain the lowest point of the grid, that is, the cluster to which the lowest point of the grid belongs or the cluster whose curvature σ is less than the point cloud cluster threshold σ th The cluster is determined to be a ground cluster.

[0092] Among them, the calculation process of the point cloud cluster curvature σ is as follows:

[0093]

[0094] Where λ1, λ2, λ3 (λ1<λ2<λ3) are the point cloud cluster eigenvalues ​​calculated by the PCA method.

[0095] However, the above methods are prone to missing ground clusters at terrain features or misclassifying non-ground clusters as ground clusters.

[0096] like Figure 6 As shown in the figure, the traditional method classifies the ground cluster B into a non-ground cluster due to the large height difference between the ground cluster B and the lowest point g2 of the grid where it is located, while the non-ground cluster C is misclassified as a ground cluster due to its small curvature.

[0097] To this end, in addition to directly designating the cluster containing the lowest point of the grid as the ground cluster, this embodiment identifies the ground cluster based on the height difference between the point cloud cluster and the lowest points of 8 adjacent grids, the horizontal distance to the edge of the adjacent grid, and the curvature feature criterion.

[0098] That is, the ground cluster GP must satisfy:

[0099] In the formula, σ i is the curvature of the i-th point cloud cluster, σ th is the point cloud cluster curvature threshold.

[0100] CP i represents the i-th point cloud cluster, is the average elevation of the i-th point cloud cluster, h gj represents the lowest point elevation of the jth neighboring grid, where i = 1, 2, ..., j = 1, ..., 8, h th Indicates the preset height difference threshold.

[0101] l i is the distance from the midpoint of the i-th point cloud cluster to the edge of the adjacent grid, l th Indicates the preset distance threshold.

[0102] More specifically, the ground cluster identification process considering fractured terrain is as follows:

[0103] First, the point cloud cluster curvature is less than σ th The point cloud cluster is judged as a candidate ground cluster.

[0104] Then the height difference between the average elevation of the point cloud cluster and the lowest points of the eight adjacent grids is calculated.

[0105] If any of the eight height differences calculated is less than the preset height difference threshold h th , then further determine the minimum value l of the distance from the point in the point cloud cluster to the edge of the adjacent grid min (For example Figure 6 l b and l c ).

[0106] If l min Less than the distance threshold l th , then the point cloud cluster is judged as a ground cluster, otherwise, the point cloud cluster is a non-ground cluster.

[0107] like Figure 6 In the example, due to the minimum distance l from the point in point cloud cluster B to the edge of the adjacent grid b Less than the distance threshold l th , the point cloud cluster B is identified as the ground cluster, and l c Greater than the distance threshold l th , therefore, the point cloud cluster C is judged as a non-ground cluster.

[0108] The above criteria can accurately identify the ground clusters in the terrain fracture area ( Figure 6 The point cloud cluster B in the figure can also avoid the non-ground cluster with smaller curvature ( Figure 6The point cloud cluster C) in the figure is misjudged as a ground cluster, which greatly improves the recognition rate of the ground cluster.

[0109] Step 3. Multi-scale hierarchical filtering.

[0110] Although the ground clusters that take terrain breaks into account can accurately identify a large number of ground points, due to the complexity of the point cloud scene, some terrain detail features are ignored because they cannot be clustered into point cloud blocks.

[0111] Therefore, in order to further identify detailed ground points, this embodiment uses the points in the captured ground cluster as ground seed points to construct an initial ground reference surface, and performs multi-scale hierarchical filtering (MHF) on this basis.

[0112] Traditional MHF uses the vertical distance from the point to be classified to the ground reference grid as the criterion for identifying ground points. However, when the grid resolution is coarse and the terrain slope at the grid is large, this criterion is significantly affected by the spatial position error. Figure 7 shown.

[0113] To this end, this embodiment proposes an MHF that takes into account the spatial position error, and the specific process is as follows:

[0114] Step 3.1. First, use the ground points in the ground cluster as ground seed points to construct the initial ground reference surface.

[0115] Step 3.2. Use thin plate spline interpolation to construct a ground reference surface with a resolution of R at each layer, and use the four adjacent ground reference surface grid center points to bilinearly interpolate the points in the non-ground cluster to calculate the interpolated elevation of the point.

[0116] Subtract the interpolated elevation from the elevation of the point. If the difference is less than the preset elevation difference threshold, it is judged as a ground point.

[0117] Step 3.3. Use the filtered ground points to update the ground reference plane in step 3.2 until there are no ground points added to this layer. The updated ground reference plane is used for filtering the next layer, and the above step 3.2 is repeated until the filtering is completed.

[0118] In this embodiment, three-layer filtering is preferably used, but it is certainly not limited to three-layer filtering. It should be noted that as the number of filtering levels increases in this embodiment, the resolution of the ground reference surface and the preset height difference threshold are gradually improved.

[0119] The effectiveness of the method of the present invention is further described in detail below through experiments.

[0120] Currently, the evaluation data used by most algorithms is the benchmark dataset released by the International Society of Photogrammetry and Remote Sensing (ISPRS), but this dataset only contains point cloud geometric three-dimensional coordinates and classification information.

[0121] With the development of point cloud acquisition equipment, the density of point cloud data has increased dramatically. At the same time, the high-precision registration of the 3D point cloud data scanned by LiDAR and the 2D image data can obtain a color point cloud with physical information. However, the point density of the ISPRS benchmark dataset is much lower than the currently produced point cloud data, and it lacks other auxiliary physical information such as color and echo number.

[0122] Therefore, the embodiment of the present invention selects four sets of high-density airborne LiDAR point cloud data with complex terrain, mixed buildings and vegetation as the research objects, and compares the filtering results of the new method with six representative filtering algorithms, including:

[0123] Multi-scale hierarchical filtering (MHF), improved simple morphology (SMRF), cloth filter (CSF), slope filter (SBF), progressively refined triangulation filter (PTD) and moving surface fitting filter (MSF).

[0124] The area of ​​the four experimental groups is 500×500m 2 , including a variety of terrain features, such as gully terrain, fracture terrain, steep slope terrain and fracture steep slope terrain, respectively Figures 9 to 12 shown.

[0125] In order to ensure the quality of the reference point cloud data, Terrascan software was first used to automatically obtain the initial filtering results, and then the misclassified points were removed through manual verification.

[0126] Digital elevation models (DEMs) generated from reference ground point data, such as Figures 13 to 16 Table 1 describes in detail the statistical information of the terrain characteristics, ground points, non-ground points and average height of objects of the four groups of experimental data.

[0127] Table 1 Statistics of the example dataset

[0128]

[0129] In order to quantitatively evaluate the filtering results, the present invention uses the confusion matrix (type I error TI, type II error T.II and total error TE) and the Kappa coefficient k as accuracy indicators. The calculation formulas of the four evaluation indicators are:

[0130] TI=b / (a+b)×100% (6)

[0131] T.II=c / (c+d)×100% (7)

[0132] TE=(b+c) / e×100% (8)

[0133] k=(p0-p c ) / (1-p c)×100% (9)

[0134] Among them, a represents the number of correctly classified ground points; b represents the number of ground points mistakenly classified as non-ground points; c represents the number of non-ground points mistakenly classified as ground points; d represents the number of correctly classified non-ground points.

[0135] e=a+b+c+d; p0=(a+b) / e;p c =((a+b)×(a+c)+(c+d)×(b+d)) / e 2 .

[0136] In order to evaluate the algorithm's ability to restore terrain, the DEM errors of each filtering method were further calculated.

[0137] First, the DEM is generated by interpolating the filtered ground points through the terrain-to-raster method, and then the DEM of each method is calculated relative to the reference DEM (such as Figures 13 to 16 )’s RMSE and MAE:

[0138]

[0139]

[0140] In the formula, e i is the elevation difference between the reference DEM and the interpolated DEM at the i-th grid, and n is the number of DEM grids.

[0141] The new method in this embodiment involves five parameters: initial resolution R of the ground reference surface, maximum window resolution W s , initial height difference threshold d th , height difference threshold h for identifying ground point clusters th and horizontal distance threshold l th .

[0142] In order to test the best performance of the new method, the optimal parameters corresponding to the four groups of experimental data and the filtering result accuracy are obtained by the joint debugging method in this embodiment as shown in Table 2.

[0143] Table 2 Optimal parameters and filtering accuracy of four sets of example data

[0144]

[0145] The results show that the new method achieves the smallest total error (3.59%) and the highest Kappa coefficient (92.58%) in Plot1, because Plot1 has a gentle terrain, sparse vegetation, and low filtering difficulty; while in Plot4, the filtering accuracy is slightly lower than other areas due to the large terrain slope and low vegetation attached (total error is 4.36%, Kappa coefficient is 87.82%). On average, the total error and Kappa coefficient of the new method are 4.03% and 90.39%, respectively, with standard deviations of 0.41 and 2.26, which proves that the new method is robust in dealing with various types of terrain. Compared with Plot1 and Plot2, the new method has a larger type I error than type II error in Plot3 and Plot 4, because the number of ground points in this area is far less than the number of non-ground points (Table 1), and a small number of misclassified ground points can easily lead to larger type I errors. However, the average Type II error of the new method differs from the Type I error by only 0.36%, indicating that the new method effectively suppresses the Type II error while reducing the Type I error.

[0146] The comparison of the total error and Kappa coefficient of the new method with six representative algorithms shows that Fig.17 and Fig.18 shown.

[0147] From the comparison results, the filtering accuracy and stability of the method of the present invention are the best. From the perspective of total error, the new method has the lowest total filtering error in different scenarios, and the smallest error floating range, followed by CSF and SBF, indicating that the new method has the best total error stability. From the Kappa coefficient, it can be seen that the new method has the highest Kappa value, followed by SMRF and SBF, which is the worst. On average, the total error of the new method is 4.03%, which is 66.7% higher than the MSF algorithm with the worst accuracy and 20.3% higher than the SMRF algorithm with the second highest accuracy; the Kappa coefficient is 90.39%, which is at least 3.0% higher than the other six methods.

[0148] The accuracy of DEM constructed by seven methods under different terrains is shown as follows: Fig.19 and Fig. 20 As shown. By comparison, it is found that the new method generates the highest accuracy of DEM in all scenarios. Among them, the accuracy of the new method in Plot2 is significantly higher than that of other methods. Compared with other areas, the mean error and average error of each method for constructing DEM in Plot4 are larger. This is because the ground points in this area are sparse, and a small number of missed ground points will lead to a large terrain loss. On average, the mean error of DEM generated by the new method is 71.5% lower than that of PTD with the worst accuracy, and 31.5% lower than that of MHF with the best accuracy; the average error is 39.3% lower than that of MHF with the best accuracy, and 78.2% lower than that of SBF with the worst accuracy.

[0149] The DEM generated by Plot4 filtering using seven filtering methods is shown in the figure below. Figures 21 to 28 The terrain in this area is complex, with a lot of low vegetation and buildings, and there is an obvious terrain break in the southwest corner, making filtering difficult.

[0150] The results show that the DEM generated by the proposed method is the closest to the reference DEM, accurately filters out buildings in complex terrain, has a high recognition rate for ground points in fractured terrain and local convex areas, and only has a few abnormal convexities on steep slopes (such as Fig. 22 This may be caused by the fact that a small amount of low vegetation on the steep slope is similar to the ground points and is mistakenly classified as ground points. In contrast, the DEM generated by MHF is rough on the steep slope terrain and has poor ability to maintain the fractured terrain (such as Fig.23 The dashed ellipse). SMRF retains too many non-ground points, resulting in a rough DEM ( Fig.24 Solid ellipse), but some terrain details are lost on the fractured terrain ( Fig.24 The CSF has a large type I error on slopes, which means that it will filter out too many terrain details, resulting in smooth terrain ( Fig.25 The dotted ellipse). MSF, SBF and PTD ( Figure 26-28 ) have large Type II errors, resulting in a rough DEM surface.

[0151] In order to further verify the influence of the innovative points of the new method on the filtering results, the present invention takes plot2, plot3 and plot4 as examples to compare the filtering accuracy of the new method with and without adding the corresponding innovative points (ground point cluster recognition taking into account terrain breaks, multi-feature point cloud clustering, and point cloud filtering taking into account spatial position errors), as shown in Table 3.

[0152] Table 3 Comparison of filtering accuracy of the new method with and without adding innovative points (%)

[0153]

[0154]

[0155] Taking Plot2 as an example, the present invention verifies the effect of ground point cluster recognition criteria with and without considering fracture terrain on the preservation of fracture terrain details. The results show (Table 3) that compared with the method without considering fracture terrain, the I, II and total errors of ground point clusters selected by the discrimination criteria taking fracture terrain into account are reduced by 74.3%, 25.6% and 52.5% respectively, and the Kappa coefficient is increased by 11.0%. The comparison of the filtering effects of using and not using multi-feature point cloud clustering with Plot3 shows that (Table 3) the former reduces the I, II and total errors by 17.9%, 8.8% and 11.6% respectively compared with the latter, and the Kappa coefficient is increased by 1.7%. The comparison of the filtering effects of the new method considering and not considering the spatial position error with Plot4 shows that (Table 3) compared with the method without considering the spatial position error, the I type error of the improved method increases slightly, but the II type and total errors are reduced by 23.7% and 11.0% respectively, and the Kappa coefficient is increased by 1.4%. Comparison of the new method with the results of six traditional filtering methods using four sets of experimental data showed that the method of the present invention had the best filtering performance in all experimental areas, with an average Kappa coefficient of 90.39% and an average total error of 4.03%, which was at least 20.3% higher than the other six representative filtering methods; at the same time, the average Type I and Type II errors of the new method differed by only 0.36%, reflecting the robustness of the method in processing point clouds in complex environments.

[0156] In summary, the new method has the smallest average total error, the best filtering performance and the highest stability.

[0157] Of course, the above description is only a preferred embodiment of the present invention, and the present invention is not limited to the above embodiments. It should be noted that all equivalent substitutions and obvious deformation forms made by any technician familiar with the field under the guidance of this specification fall within the essential scope of this specification and should be protected by the present invention.

Claims

1. A point cloud hierarchical filtering method for complex terrain areas assisted by multi-feature clustering, characterized in that: The steps include: Step 1. Grid the original point cloud data, perform multi-feature clustering on the point cloud in each independent grid in turn, obtain the multi-feature clustering result, map the clustering result to three-dimensional space, and obtain the three-dimensional coordinate distribution of the point cloud cluster; Step 1.

1. First, use the least squares method to fit the ground reference surface and calculate the fitting residual of each point; Step 1.

2. Convert the RGB color information of the point cloud into a single-value visible light band difference vegetation index VDVI; Step 1.

3. Use the matrix [elevation, fitting residual, VDVI] as the point cloud virtual coordinates (z, h, VDVI) and perform DBSCAN clustering based on Euclidean distance; Step 2. Identify ground clusters from the three-dimensional coordinate distribution of point cloud clusters using a criterion that takes into account fracture terrain; The ground cluster is identified by the height difference between the point cloud cluster and the lowest point of the eight adjacent grids, the horizontal distance to the edge of the adjacent grid, and the curvature characteristic criterion, that is, the ground cluster GP must satisfy the following formula; In the formula, CP i represents the i-th point cloud cluster; σ i is the curvature of the i-th point cloud cluster, σ th is the point cloud cluster curvature threshold; is the average elevation of the i-th point cloud cluster, h gj represents the lowest point elevation of the jth neighboring grid, i=1,2,…,j=1,…,8, h th Indicates the preset height difference threshold; l i is the distance from the midpoint of the i-th point cloud cluster to the edge of the adjacent grid, l th Indicates the preset distance threshold; Step 3. Use the ground points in the ground cluster captured in step 2 as ground seed points to construct an initial ground reference surface, and further refine the ground points in the original point cloud by combining a multi-scale hierarchical filtering method; Step 3.

1. First, use the ground points in the ground cluster as ground seed points to construct the initial ground reference surface; Step 3.

2. Use thin plate spline interpolation to construct a ground reference surface with a resolution of R at each layer, and use the four adjacent ground reference surface grid center points to perform bilinear interpolation on the points in the non-ground cluster to calculate the interpolated elevation of the point; Subtract the interpolated elevation from the elevation of the point. If the difference is less than the preset elevation difference threshold, it is determined to be a ground point. Step 3.

3. Use the filtered ground points to update the ground reference plane in step 3.2 until there are no ground points added to this layer. The updated ground reference plane is used for filtering the next layer, and the above step 3.2 is repeated until the filtering is completed.

2. The complex terrain area point cloud hierarchical filtering method according to claim 1 is characterized in that: In step 1, before the original point cloud data is gridded, Firstly, the point cloud data of non-last echo among multiple echoes are removed from the original point cloud data according to the echo information.

3. The complex terrain area point cloud hierarchical filtering method according to claim 1 is characterized in that: In step 1.3, the specific process of DBSCAN clustering is as follows: Step 1.3.

1. Mark all points in the original point cloud data as unclassified points, and calculate the Euclidean distance between any two points in the unclassified points: The Euclidean distance calculation formula between any two points in the clustering process is as follows: Among them, d ij Represents the Euclidean distance between point i and point j; z i 、z j Represents the elevation of point i and point j, r i 、r j Represent the fitting residual values ​​of point i and point j respectively, VDVI i 、VDVI j Represent the VDVI values ​​of point i and point j respectively, α, β, γ are the proportional factor coefficients, and α, β, γ are all greater than 0; Step 1.3.

2. Randomly select a point P from the unclassified points and mark it as a classified point. Then search the ε neighborhood N of the point P and calculate the number of points in the neighborhood N. If the number of points in the neighborhood N is greater than the minimum number of points P ts , then mark point P as a core point, create a new cluster C, and add point P to cluster C; Among them, ε represents the radius of the search neighborhood point, N is the set of points in the ε neighborhood, P ts Indicates the minimum number of points; Step 1.3.

3. Traverse the ε neighborhood M of all unclassified points q in N. If the number of points in the neighborhood M is greater than the minimum number of points P ts , then add the point in M ​​to N, and if point q does not belong to any cluster, then add point q to cluster C; Step 1.3.

4. Repeat step 1.3.3 until no new points are added to cluster C, then output cluster C and return to step 1.3.2; if no new cluster is output, the remaining points are classified as noise points until there are no points marked as unclassified, then go to step 1.3.5; Step 1.3.

5. Output result: The clusters are divided into C′={C1,C2,C3,C4,……,C k }; Among them, C′ represents the set of point cloud clusters, C1, C2, C3, C4, ..., C k They represent the divided point cloud clusters respectively.

4. The complex terrain area point cloud hierarchical filtering method according to claim 3 is characterized in that: In step 1.1, the calculation process of the fitting residual of each point is: taking the lowest point of each grid as the initial ground point, and using the initial ground points of 9 adjacent grids to fit the tangent plane to calculate the fitting residual of each point in the grid.

5. The complex terrain area point cloud hierarchical filtering method according to claim 4 is characterized in that: Set a point P i If it is located in the (I, J) grid on the DTM surface, its 9 nearest neighboring grids are (I-1, J-1), (I-1, J), (I-1, J+1), (I, J-1), (I, J), (I, J+1), (I+1, J-1), (I+1, J), and (I+1, J+1); Set point P i The corresponding tangent plane equation is z = Ax + By + C, then P i The fitted residual r i for: r i =z i -(Ax i +By i +C); Where (x i ,y i , z i ) is point P i The three-dimensional coordinates of , A, B, and C are the coefficients of the plane equation.

6. The complex terrain area point cloud hierarchical filtering method according to claim 3 is characterized in that: In step 1.2, the calculation formula of VDVI is: Where R, G, and B represent the red, green, and blue color values ​​of the point cloud, respectively.

7. The complex terrain area point cloud hierarchical filtering method according to claim 1 is characterized in that: In step 2, the specific process of identifying the category of the point cloud cluster from the three-dimensional coordinate distribution of the point cloud cluster is as follows: Step 2.

1. First, cluster the point cloud with a curvature less than σ th The point cloud cluster is judged as a candidate ground cluster; Step 2.

2. Then calculate the height difference between the average elevation of the point cloud cluster and the lowest points of the eight adjacent grids; If any of the eight height differences calculated is less than the preset height difference threshold h th , then further determine the minimum value l of the distance from the point in the point cloud cluster to the edge of the adjacent grid min ; If l min Less than the distance threshold l th , then the point cloud cluster is judged as a ground cluster; otherwise, the point cloud cluster is a non-ground cluster.

8. The complex terrain area point cloud hierarchical filtering method according to claim 1 is characterized in that: The calculation process of the point cloud cluster curvature σ is as follows: Where λ1, λ2, λ3 represent the eigenvalues ​​of the point cloud cluster calculated by the PCA method, where λ1<λ2<λ3.

Images