Multi-scale hierarchical efficient filtering method considering topographic features

Abstract

The invention discloses a multi-scale hierarchical efficient filtering method considering topographic features, which can efficiently extract a large number of uniformly distributed accurate ground seed points and ensure the accuracy of an initial topographic surface by introducing a method of combining a sliding grid technology and an M estimation sample consistency algorithm. Then, constructing a ground reference surface grid according to the ground points subjected to iterative recognition, and quickly generating a ground reference surface by adopting a grid-based global weighted finite difference TPS method; and finally, establishing a self-adaptive threshold value considering topographic features to capture ground points, thereby effectively improving the point cloud filtering precision in a complex scene. The filtering method is stable in performance and good in filtering effect in a large-scale scene with complex terrain and mixed ground features, and the filtering method is high in filtering precision, high in calculation efficiency and high in terrain adaptive capacity. Meanwhile, the filtering method is high in operation speed and has the performance advantage of processing large complex scenes.

Description

Technical Field

[0001] The invention belongs to the technical field of point cloud filtering and relates to a multi-scale hierarchical efficient filtering method taking into account terrain features. Background Art

[0002] With the rapid development of advanced multi-source mapping technologies in the air, sky and ground, laser scanning (Light Detection And Ranging, LiDAR) technology provides the most effective technical means for obtaining three-dimensional geographic information with its advantages of high speed, high efficiency and high precision. However, since the original LiDAR point cloud is mixed with ground information and ground object information such as vegetation, buildings, and vehicles, high-precision point cloud filtering is the premise and key to three-dimensional geographic information modeling.

[0003] At present, scholars at home and abroad have proposed a variety of point cloud filtering algorithms. According to the working principles of the existing algorithms, they can be roughly divided into five categories: morphology-based filtering, surface-based filtering, slope-based filtering, block-based filtering and machine learning-based filtering. Among them, morphology-based filtering implements an open operation by eroding and dilating the terrain to filter out non-ground points in the local window. However, this type of filtering has poor accuracy in complex terrain areas with mixed terrain objects. Surface-based filtering generates a reference surface according to the interpolation model, and then iteratively classifies the unclassified point set, but this type of filtering is inefficient for large-scale data calculations and has poor terrain details in steep slope areas. Slope-based filtering extracts ground points by comparing the slope difference between two points, and classifies points with a slope greater than a specific threshold as non-ground points. However, this type of filtering is sensitive to parameter settings and is difficult to apply to sudden terrain. Block-based filtering uses a clustering segmentation method to generate a point cloud block set, and then extracts ground blocks with the help of geometric and topological information. It can effectively retain ground points with terrain discontinuities, but its effectiveness is highly dependent on segmentation performance and has poor accuracy in areas with significant terrain changes. Machine learning-based filtering formulates the filtering task as a binary classification problem and uses sufficient labeled samples and advanced machine learning models to extract ground points from the original point cloud. However, the performance of this type of filtering method is limited by the quality and quantity of training samples, and the training time is lengthy.

[0004] In summary, existing research has failed to solve the problems of low point cloud filtering accuracy, slow calculation speed and poor terrain adaptation under large-scale data conditions in complex scene areas. Therefore, it is necessary to propose a multi-scale hierarchical efficient filtering method that takes into account terrain characteristics to solve the above technical problems existing in the existing technology. Summary of the invention

[0005] The purpose of the present invention is to provide a multi-scale hierarchical efficient filtering method that takes into account terrain characteristics, so as to solve the problems of low point cloud filtering accuracy, slow calculation speed and poor terrain adaptability under large-scale data conditions in complex scene areas.

[0006] To achieve the above object, the present invention provides the following technical solutions:

[0007] A multi-scale hierarchical efficient filtering method taking into account terrain characteristics comprises the following steps:

[0008] I. Use radius filtering and elevation histogram method to remove local and global outliers in the original point cloud;

[0009] II. Based on the point cloud data after removing outliers in step I, determine the rectangular study area; further use the sliding grid technology and M-estimation sample consistency algorithm to obtain selected ground seed points from the point cloud data of the rectangular study area;

[0010] III. Construct a ground reference surface based on the selected ground seed points obtained in step II. The specific steps are as follows:

[0011] III.1. Construct a grid with resolution h based on selected ground seed points to cover the entire rectangular study area;

[0012] III.2. Use the lowest selected ground seed point elevation value in each grid as the grid value;

[0013] III.3. Use global weighted finite difference TPS to interpolate the mesh surface, while reducing the impact of misclassified non-ground points and obtaining a high-quality ground reference surface;

[0014] IV. using cubic spline interpolation to calculate the elevation difference between the point to be classified and the corresponding position in the ground reference surface constructed in step III, and determining the point whose elevation difference is less than a preset filtering threshold as a ground point;

[0015] V. Repeat steps III.2 to IV on the ground reference plane with a resolution of h until no ground points are identified;

[0016] VI. Set the resolution h = h / 2, the initial elevation threshold t0 = t0 + 0.1, and repeat steps III.1 to V until the predetermined ground reference surface resolution is reached.

[0017] Preferably, the step I specifically comprises:

[0018] First, a kd-tree is constructed to organize the point cloud. Then, the number of neighboring points of each point within a specified range is searched. If the number of neighboring points of a point is less than a given threshold, the point is marked as a local outlier and removed.

[0019] Then, an elevation histogram is built for the remaining points, and global outlier points below the specified elevation threshold are removed.

[0020] Preferably, the step II is specifically:

[0021] First, based on the maximum value of the x-coordinate and the maximum value of the y-coordinate in the point cloud data after the outliers are removed in step I, the length and width of the rectangular research area are determined;

[0022] Then, a square grid with a length and width of W is established in the rectangular area, and the square grid is slid from left to right and from bottom to top with a step length of S, and the lowest point in each square grid is selected as the initial ground seed point;

[0023] Finally, the initial ground seed points are selected using the M-estimation sample consistency algorithm to obtain selected ground seed points.

[0024] Preferably, in step II, the length or width W of the square grid is greater than the minimum side length of the largest non-ground object in the rectangular study area; the step length S is the minimum side length of the largest non-ground object in the rectangular study area. times.

[0025] Preferably, in step II, in the process of selecting the lowest point in each sliding grid as the initial ground seed point, when the lowest point of the current grid overlaps with the lowest point of the previous grid, the next lowest point is selected as the initial ground seed point of the current grid.

[0026] Preferably, in step II, the specific steps of selecting the initial ground seed points using the M estimation sample consistency algorithm are:

[0027] II.1. For any point P among the initial ground seed points, search for the K neighboring points closest to point P;

[0028] II.2. Randomly select three non-collinear neighbor points from the K neighbor points as the minimum sample set, and use the least squares method to calculate the fitting plane parameters corresponding to the minimum sample set;

[0029] II.3. Preset the distance threshold, calculate the distance between the remaining K-3 neighbor points and the fitting plane, and mark the neighbor points whose distance to the fitting plane is less than the preset distance threshold as in-plane points;

[0030] Use the minimum cost function to iteratively obtain more in-plane points and the best fitting plane parameters;

[0031] II.4. Neighbor points whose distance to the fitting plane is greater than the preset distance threshold are marked as non-ground points and removed.

[0032] Preferably, in step III.1, when setting the first level of filtering, the resolution h is taken as half of the average point spacing of the selected ground seed points, and when filtering the next level, the resolution h continues to be halved until the predetermined ground resolution is reached.

[0033] Preferably, in step III.3, the specific process of using global weighted finite difference TPS to interpolate the grid surface and obtain a high-quality ground reference surface is:

[0034] The TPS objective function is expressed as follows using the second-order finite difference operator:

[0035]

[0036] In the formula, f=[f 11 ···f 1m ;f 21 ···f 2m ;···;f n1 ···f n-m ] T , f ij is the simulated value at the grid (i, j), n and m are the number of rows and columns of the grid in the study area respectively; W is a diagonal matrix diag(w i ), which contains the data Z i The corresponding weight w i ∈[0,1], where data Z i is the filling value for each grid; λ is the smoothing parameter determined by generalized cross validation; D = [D1, 2D2, D3] T , where D1, D2 and D3 are the second-order finite difference f xx 、f xy and f yy Matrix expression of ;

[0037] Based on the DCT method and the finite difference strategy, the following formula can be obtained by solving formula (1):

[0038] f (k+1) =IDCT(ΓDCT(W(zf (k) )+f (k) )) (2)

[0039] Where k is the number of iterations, DCT represents discrete cosine transform, IDCT represents inverse discrete cosine transform; Γ is a diagonal matrix, in which each non-zero element is a function of the smoothing parameter λ and the eigenvalue D;

[0040] In the first iteration, the grid without sample values ​​is set to w = 0, and the grid with sample values ​​is set to w = 1; in subsequent iterations, w i The calculation formula is:

[0041]

[0042] In the formula, r i is the standardized residual, expressed as

[0043] e i is the simulated residual, h i is the diagonal element of the hat matrix H, denoted as h i =h=Tr(H) / mn; σ is the robust estimate of the residual standard deviation, expressed as σ=1.843MAD, where MAD is the median of the absolute deviation.

[0044] Preferably, step IV comprises the following steps:

[0045] IV.1. Search for the 8 nearest ground points around the point to be classified and use cubic spline interpolation to calculate the interpolation elevation of the point;

[0046] IV.2. Calculate the elevation difference between the interpolated elevation of the point to be classified and the known elevation. If the elevation difference is less than the set filtering threshold, the point to be classified is marked as a ground point; if the elevation difference is greater than the set filtering threshold, the point to be classified is marked as a non-ground point; the calculation formula of the filtering threshold is: th = t0 + c × (t c +t s );

[0047] In the formula, c is the coefficient, t0, t c ,t s They are the initial elevation threshold, ground curvature value, and slope value respectively;

[0048] t c The value of is obtained by calculating the curvature value of each grid cell in the ground reference surface and normalizing it;

[0049] t s The value of is obtained by calculating the slope value of each grid in the ground reference surface and normalizing it.

[0050] Compared with the prior art, the present invention has the following beneficial effects:

[0051] As described above, the present invention discloses a multi-scale hierarchical efficient filtering method that takes terrain features into account. This method firstly uses a sliding grid technique combined with an M-estimation sample consistency algorithm to efficiently extract a large number of evenly distributed precise ground seed points, thereby ensuring the accuracy of the initial terrain surface; secondly, a ground reference surface grid is constructed based on iteratively identified ground points, and a grid-based global weighted finite difference TPS method is used to quickly generate a ground reference surface; finally, an adaptive threshold that takes terrain features into account is established to capture ground points, thereby effectively improving the accuracy of point cloud filtering in complex scenes. The filtering method of the present invention is robust in large-scale scenes with complex terrain and mixed objects, and has a good filtering effect. The filtering method has high filtering accuracy, fast computational efficiency, and strong terrain adaptability; at the same time, the filtering method runs quickly and has the performance advantage of processing large and complex scenes. BRIEF DESCRIPTION OF THE DRAWINGS

[0052] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings required for use in the embodiments or the description of the prior art are briefly introduced below.

[0053] Figure 1 It is a flow chart of a multi-scale hierarchical efficient filtering method taking into account terrain characteristics in the present invention;

[0054] Figure 2 It is a schematic diagram of selecting ground seed points based on the combination of sliding grid technology and M estimation sample consistency algorithm in an embodiment of the present invention;

[0055] Figure 3 A schematic diagram of constructing a ground reference plane in an embodiment of the present invention;

[0056] Figure 4 A schematic diagram showing a comparison of identifying ground points based on fixed threshold filtering and based on adaptive threshold filtering in an embodiment of the present invention;

[0057] Figure 5 are three sample areas Test I, Test II and Test III of the OpenGF data set in an embodiment of the present invention;

[0058] Figure 6 A schematic diagram showing the accuracy comparison between the method of the present invention and other filtering methods;

[0059] Figure 7 A schematic diagram showing the error comparison between the method of the present invention and other filtering methods;

[0060] Figure 8 Schematic diagram of the comparison of DEMs generated after Test III filtering by various filtering algorithms in the embodiments of the present invention. DETAILED DESCRIPTION

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

[0062] like Figure 1 As shown, this embodiment describes a multi-scale hierarchical efficient filtering method taking into account terrain features, comprising the following steps:

[0063] I. Data preprocessing.

[0064] As the laser pulse propagates from the air to the ground, a large number of non-ground objects (such as vegetation and building edges) usually generate two or more echoes, and the probability of non-last echoes being ground points is low; therefore, this embodiment eliminates non-last echo point clouds to improve the algorithm calculation efficiency.

[0065] Since the algorithm in this embodiment uses the lowest point in the local range as the ground seed point, it is particularly sensitive to low outliers; therefore, this embodiment first uses radius filtering to remove local outliers of each point in the original point cloud, and then uses the elevation histogram to remove global outliers of the original point cloud. The specific process is as follows:

[0066] First, a kd-tree is constructed to organize the point cloud, and the number of neighboring points of each point within a specified range is searched. If the number of neighboring points of a point is less than a given threshold, the point is marked as a local outlier and removed.

[0067] Then, an elevation histogram is built for the remaining points, and global outlier points below the specified elevation threshold are removed.

[0068] II. Selection of ground seed points.

[0069] The traditional virtual grid method can only obtain a small number of initial ground seed points, and there is a lack of necessary structural information between adjacent grids. Therefore, this embodiment constructs a sliding grid method for ground point search, such as Figure 2 (a); First, based on the point cloud data after removing outliers in step I, determine the rectangular study area; then use the sliding grid technology and M estimation sample consistency algorithm to obtain selected ground seed points from the point cloud data of the rectangular study area. The specific process is as follows:

[0070] First, based on the maximum value of the x-coordinate and the maximum value of the y-coordinate in the point cloud data after the outliers are removed in step I, the length and width of the rectangular research area are determined;

[0071] Then, a square grid with a length and width of W is established in the rectangular area, and the square grid is slid from left to right and from bottom to top with a step length of S, and the lowest point in each square grid is selected as the initial ground seed point; wherein W is greater than the minimum side length of the largest non-ground object in the study area, and in this embodiment, W is slightly greater than the minimum side length of the largest non-ground object in the study area; S is the minimum side length of the largest non-ground object in the rectangular study area. times.

[0072] This method inevitably leads to repeated lowest points. To ensure that the lowest point of each grid is not repeated, when sliding the grid to obtain the lowest point, if the lowest point index of the current grid is repeated with that of the previous grid, the next lowest point is selected as the initial ground seed point of the current grid.

[0073] Finally, the M-estimation sample consistency algorithm is used to select the initial ground seed points to obtain the selected ground seed points, such as Figure 2 (b), the specific steps are:

[0074] II.1. For any point P among the initial ground seed points, search for the K neighboring points closest to point P;

[0075] II.2. Randomly select three non-collinear neighbor points from the K neighbor points as the minimum sample set, and use the least squares method to calculate the fitting plane parameters corresponding to the minimum sample set;

[0076] II.3. Calculate the distance between the remaining K-3 neighbor points and the fitting plane and set the distance threshold. The neighbor points whose distance is less than the distance threshold are marked as in-plane points. Use the minimum cost function to iterate to obtain more in-plane points and the best fitting plane parameters.

[0077] II.4. Neighboring points whose distance is greater than the distance threshold are marked as non-ground points and removed.

[0078] III. Construction of ground reference surface.

[0079] The core of point cloud filtering based on ground reference surface is to construct ground reference surface by interpolating ground points. Therefore, the choice of interpolation method directly determines the accuracy and computational efficiency of filtering results. Among them, TPS is widely used in terrain modeling due to its advantages such as good stability of calculation results and strong fitting ability. However, the computational efficiency of TPS based on point-by-point local analysis is very slow in large-scale scenarios. Therefore, this embodiment constructs a grid-based global weighted finite difference TPS method for ground reference surface construction, such as Figure 3 The specific steps are as follows:

[0080] III.1. Ground reference surface grid construction.

[0081] A grid with a resolution of h is constructed based on the selected ground seed points to cover the entire rectangular study area, such as Figure 3 As shown in (a), when setting the first level of filtering, the resolution h is half of the average point spacing of the selected ground seed points, and when filtering the next level, the resolution h continues to be halved until the predetermined ground reference surface resolution is reached.

[0082] III.2. Grid value filling.

[0083] Since the present embodiment uses the selected ground seed points as the data source to construct the reference surface, a single grid may contain multiple selected ground seed points. Therefore, the present embodiment uses the elevation value of the lowest point in the grid to fill each grid value update, thereby effectively avoiding the error accumulation caused by the misclassified non-ground points in the iteration process. Figure 3 (b) as shown.

[0084] III.3. Ground reference surface simulation.

[0085] Under the condition of mesh surface, the TPS objective function is expressed as follows using the second-order finite difference operator:

[0086]

[0087] In the formula, f=[f 11 ···f 1m ;f 21 ···f 2m ;···;f n1 ···f n-m ] T , f ij is the simulated value at the grid (i, j), n and m are the number of rows and columns of the grid in the study area respectively; λ is the smoothing parameter determined by generalized cross validation; D = [D1, 2D2, D3] T , where D1, D2 and D3 are the second-order finite difference f xx 、f xy and f yy The matrix expression of .

[0088] By finding the extreme value of f with the above formula and setting it to zero, we can get: f-z+λD T Df=0, therefore, f=H -1 z, where H = I + λD T D.

[0089] Solving the above formula based on the DCT method and finite difference strategy yields:

[0090] f = IDCT(ΓDCT(z))

[0091] Where Γ is a diagonal matrix, in which each non-zero element is a function of the smoothing parameter λ and the eigenvalue D; DCT represents discrete cosine transform, and IDCT represents inverse discrete cosine transform.

[0092] like Figure 3 As shown in (c), in order to reduce the impact of misclassified non-ground points on the construction of the reference surface, this embodiment uses global weighted finite difference TPS to interpolate the grid surface, which is specifically expressed as:

[0093]

[0094] Where W is a diagonal matrix diag(w i ), which contains the data Z i The corresponding weight w i ∈[0,1], where Z i is the filling value for each grid; f=[f 11 ···f 1m ;f 21 ···f 2m ;···;f n1 ···f n-m ] T , f ij is the simulated value at the grid (i, j), n and m are the number of rows and columns of the grid in the study area respectively; λ is the smoothing parameter determined by generalized cross validation; D = [D1, 2D2, D3] T , where D1, D2 and D3 are the second-order finite difference f xx 、f xy and f yy The matrix expression of .

[0095] Similarly, by finding the extreme value of f in formula (1) and setting it to zero, we can obtain: (W+λD T D) f = Wz;

[0096] Therefore, f (k+1) =H -1 (W(zf (k) )+f (k) ), where k is the number of iterations, H = I + λDTD;

[0097] Similarly, based on the DCT method and the finite difference strategy, we can solve formula (1) to obtain:

[0098] f (k+1) =IDCT(ΓDCT(W(zf (k) )+f (k) )) (2)

[0099] Where k is the number of iterations; Γ is a diagonal matrix, in which each non-zero element is a function of the smoothing parameter λ and the eigenvalue D; DCT represents discrete cosine transform; and IDCT represents inverse discrete cosine transform.

[0100] In the first iteration, the grid without sample values ​​is set to w = 0, and the grid with sample values ​​is set to w = 1; in subsequent iterations, w i The calculation formula is:

[0101]

[0102] In the formula, r i is the standardized residual, expressed as

[0103] e i is the simulated residual, h i is the diagonal element of the hat matrix H, denoted as h i =h=Tr(H) / mn; σ is the robust estimate of the residual standard deviation, expressed as σ=1.843MAD, where MAD is the median of the absolute deviation.

[0104] Theoretically, the computational complexity of solving formula (2) is about O((m·n)log(m·n)), while the TPS computational complexity of the analytical formula is O((m·n) 3 ), so the computational complexity of the former is about two orders of magnitude lower than that of the latter.

[0105] IV. Ground point recognition.

[0106] In the process of identifying ground points, the point to be classified P is usually used i The elevation difference (elevation difference) of the grid corresponding to the ground reference surface is used as the criterion to identify the ground point. When the grid resolution is large and the terrain changes significantly at the grid, the point to be classified P i The elevation difference to the center of the grid is significantly affected by the spatial position error. Therefore, this embodiment uses cubic spline interpolation to calculate the elevation difference between the point to be classified and the corresponding position in the ground reference surface constructed in step III, and determines the point with an elevation difference less than a preset filtering threshold as a ground point. The specific steps are as follows:

[0107] IV.1. Searching for the point P to be classified i The interpolation elevation of the eight nearest ground points is calculated using cubic spline interpolation;

[0108] IV.2. Calculate the points to be classified P i The elevation difference between the interpolated elevation and the known elevation. If the elevation difference is less than the set filtering threshold, the point to be classified P i Marked as a ground point; if the elevation difference is greater than the set filtering threshold, the point to be classified P iMark as non-ground point.

[0109] Affected by the undulation of the terrain, using only a single height difference threshold in steep slope areas can easily lead to ground point recognition errors. Therefore, the present invention designs a filtering threshold that takes into account the complexity of the terrain, and the formula is:

[0110] th=t0+c×(t c +t s )

[0111] In the formula, c is the coefficient, t0, t c ,t s are the initial elevation threshold, the ground curvature value, and the slope value, respectively. The smoothness characteristic of TPS itself will cause the generated ground reference surface to be smooth. Therefore, this embodiment calculates the curvature value of each grid cell in the ground reference surface and normalizes it to obtain t c Since even a small horizontal or vertical displacement may produce significant errors in steep terrain areas, this embodiment calculates the slope value of each grid on the ground reference surface and normalizes it as t s In view of the diverse terrain features, the accuracy of ground point recognition in complex scenes can be significantly enhanced by adaptively adjusting the filtering threshold strategy.

[0112] like Figure 4 As shown in (a), ground points are identified based on fixed threshold filtering. Affected by the spatial position error and terrain complexity, point P1 among the points to be classified is located outside the fixed threshold as a ground point and is mistakenly classified as a non-ground point, while point P2 is located within the fixed threshold as a non-ground point and is mistakenly classified as a ground point.

[0113] like Figure 4 As shown in (b), ground points are identified based on adaptive threshold filtering; by eliminating spatial position errors and adopting a filtering threshold that takes into account the complexity of the terrain, P1 can be correctly classified as a ground point and P2 can be correctly classified as a non-ground point.

[0114] V. Inner iteration.

[0115] Repeat steps III.2 to IV.2 on the ground reference plane with a resolution of h until there are no points to be classified, that is, no ground points are identified.

[0116] VI. External iteration.

[0117] In step V, when there are no points to be classified, the next level of ground reference surface construction begins.

[0118] In which, the resolution of the next level of ground reference surface is set to h=h / 2, and the initial elevation threshold t0=t0+0.1; steps III.1 to V are repeated until a predetermined ground reference surface resolution or a predetermined level is reached.

[0119] The effectiveness of the multi-scale hierarchical efficient filtering method proposed in the present invention that takes into account terrain characteristics is verified below.

[0120] First, the ISPRS publicly released benchmark dataset is used as the experimental data. This dataset uses the Optech ALTM scanner to collect samples and provides 15 samples from 7 regions, including various objects such as buildings, vegetation, bridges, and terrain features such as rugged terrain, steep slopes, and terrain discontinuities, as shown in Table 1.

[0121] Since its release, the ISPRS benchmark dataset has been widely used in point cloud filtering algorithm research. In this regard, this embodiment compares the calculation results with the common filtering methods for processing the ISPRS benchmark data.

[0122] In addition, with the rapid development of hardware equipment performance, the density of point cloud data has increased dramatically. Therefore, in order to reflect the performance of the algorithm in the present invention in processing large-scale data, the ultra-large-scale point cloud dataset OpenGF is selected as the research object, and the filtering results of the algorithm in the present invention are compared with three representative filtering algorithms, including improved simple morphology (SMRF), cloth simulation filter (CSF) and slope filter (SBF). Figure 5 ,OpenGF includes 3 sample areas (Test I, Test II and TestIII), covering different terrain types; Figure 5 (a) is Test I, Figure 5 (b) Test II, Figure 5 (c) is Test III.

[0123] Table 1 Statistical information of the experimental data set in the present invention

[0124]

[0125] In order to quantitatively evaluate the filtering results, this embodiment uses four accuracy indicators, including type I error (TI), type II error (T.II), total error (TE) and Kappa coefficient (K). The specific calculation formula is:

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

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

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

[0129] K=(P0-P c ) / (1-P c )×100%

[0130] Where a represents the number of correctly classified ground points; b represents the number of ground points misclassified as non-ground points; c represents the number of non-ground points misclassified as ground points; d represents the number of correctly classified non-ground points; e = a + b + c + d; p0 = (a + b) / e; p c =[(a+b)×(a+c)+(c+d)×(b+d)] / e 2 .

[0131] The method of the present invention includes four key parameters: the window size (W) in the initial ground seed point selection, the predetermined grid resolution h, the initial elevation threshold (t0) and the coefficient (c) in the terrain adaptive threshold. The optimal parameters of the method of the present invention and the corresponding filtering accuracy are shown in Table 2. The results show that the average total error and Kappa coefficient of the method of the present invention are 2.31% and 92.25% respectively, and the medians are 1.90% and 93.75% respectively, which indicates that the method of the present invention has good performance on the ISPRS benchmark dataset. Among the 15 samples, samp61 has the smallest total error and samp31 has the highest Kappa coefficient, because the samp61 and samp31 samples have flat terrain and simple scenes.

[0132] Table 2 Parameter settings and calculation accuracy of the proposed method for ISPRS data set

[0133]

[0134] Comparison of the results of the method of the present invention with those of common filtering methods shows that, among the 15 samples, the method of the present invention has the highest accuracy in 8 samples (samp11, 12, 21, 22, 51, 52, 53, 61), and the other 7 samples are also close to the highest accuracy of the comparison algorithm, as shown in Table 3. The kappa coefficient of the method of the present invention is 92.3%, which is at least 1.65% higher than the accuracy of the previous method.

[0135] Table 3 Kappa coefficient (%) comparison between the method of the present invention and other algorithms (bold represents the best result)

[0136]

[0137] Method 1: PINGEL T J, CLARKE K C, MCBRIDE W A. An improved simple morphological filter for the terrain classification of airborne LIDAR data[J]. ISPRS Journal of Photogrammetry and Remote Sensing, 2013, 77:21-30.

[0138] Method 2: HU Han, DING Yulin, ZHU Qing, et al. An adaptive surface filter for airborne laser scanning point clouds by means of regularization and bending energy[J]. ISPRS Journal of Photogrammetry and Remote Sensing, 2014, 92:98-111.

[0139] Method 3: ZHANG Wuming, QI Jianbo, WAN Peng, et al. An Easy-to-Use Airborne LiDAR Data Filtering Method Based on Cloth Simulation[J]. Remote Sensing, 2016, 8(6):501.

[0140] Method 4: KUMAR B, YADAV M, LOHANI B, et al. A two-stage algorithm for ground filtering of airborne laser scanning data[J]. International Journal of Remote Sensing, 2018, 39(20):6757-6783.

[0141] Method 5: S, CORDERO M, MIRANDA D. Hybrid Overlap Filter for LiDAR Point Clouds Using Free Software[J]. Remote Sensing, 2020, 12(7):1051.

[0142] Method 6: CHEN Chuanfa, CHANG Bingtao, LI Yanyan, et al.Filtering airborneLiDAR point clouds based on a scale-irrelevant and terrain-adaptive approach[J].Measurement, 2021,171:108756.

[0143] Table 4 shows the optimal parameters and filtering accuracy of the method of the present invention in OpenGF. On average, the total error and standard deviation of the method of the present invention are low, and the Kappa coefficient is high, which proves that the method of the present invention still has excellent performance when processing large scenes. Among the three samples, the method of the present invention achieves the highest Kappa coefficient in the Test I sample area and obtains the lowest total error in the Test III sample.

[0144] Table 4 Optimal parameters and filtering accuracy of the method of the present invention for OpenGF data

[0145]

[0146] The comparison with three representative methods shows that the method of the present invention has the lowest total error and the highest Kappa coefficient among the three samples. From the perspective of total error, the total filtering error of the method of the present invention is the lowest in different scenarios, and the error floating range is the smallest; SMRF is second, and SBF is the worst. From the Kappa coefficient, it can be seen that the Kappa value of the method of the present invention is the highest, followed by SMRF, and SBF is the worst. On average, the total error of the method of the present invention is 2.74%, which is 74.8% higher than the SBF method with the worst accuracy, and 26.3% higher than the SMRF method with the second highest accuracy; the Kappa coefficient is 93.74%, which is at least 2.2% higher than the other three methods. Figure 6 (a) and 6(b).

[0147] like Figure 7As shown (blue is type I error, red is type II error), it can be seen from the error distribution comparison diagram that the method of the present invention has the best filtering effect in Test I, Test II and Test III samples. In the Test I sample, the method of the present invention and the SMRF method have the highest accuracy, while CSF has a large type I error in mountainous areas, and SBF has a large type II error in building areas. Although the method of the present invention has a large type II error in the Test II sample, its error remains the lowest compared to the other three methods, and the three compared methods have large type II errors in large buildings. The Test III sample contains dense vegetation on steep slopes, but the method of the present invention can accurately capture ground points on steep slopes with sharp changes in elevation, with almost no errors in this area, while traditional methods have errors of varying degrees. In general, the method of the present invention is robust in large-scale scenes with complex terrain and mixed objects, and has a good filtering effect.

[0148] Figure 8 The DEM results of four different filtering methods in Test III are shown. The natural scene of this sample is extremely complex, including dense vegetation coverage and rapid changes in elevation, which leads to a significant sparsity in the distribution of ground points, posing a huge challenge to existing filtering methods. The experimental results show that the DEM of the proposed method is highly consistent with the reference DEM, and can accurately capture ground points in environments with dense vegetation and steep terrain. Its ground point recognition ability is particularly outstanding in areas with significant terrain undulations and steep slopes (such as Figure 8 (b)). In contrast, the DEMs generated using SMRF and CSF are relatively rough in steep terrain areas and have large Type I errors, resulting in the loss of some key terrain details (e.g. Figure 8 (c) and Figure 8 (d)). However, the SBF filter has the worst effect, retaining too many non-ground points, resulting in a large number of convexities in the DEM (such as Figure 8 (e)).

[0149] In addition, computational efficiency is an important indicator to measure the ease of use of an algorithm. In this regard, Table 5 gives the computational time of each filtering method for processing OpenGF data. The results show that, on average, SBF is the fastest, because its principle is simple and filtering is mainly based on the neighborhood search of the points to be classified. The computational speed of the method of the present invention is second, because it constructs a ground reference surface through finite differences. SMRF has a low computational efficiency, because although its processing object is a rasterized point cloud, the complex window changes increase the computational time. CSF has the slowest computational efficiency, because this method uses a cloth model to iteratively simulate the ground.

[0150] Table 5 Computation time of each filtering method for processing OpenGF data (s)

[0151]

[0152] In view of the challenges of low filtering accuracy, slow computational efficiency and poor terrain adaptation in complex scenes in existing filtering methods, this paper proposes a multi-scale hierarchical efficient filtering method that takes into account terrain features. This method aims to efficiently extract a large number of uniformly distributed accurate ground seed points by introducing a sliding grid technology combined with an M-estimation sample consistency algorithm to ensure the accuracy of the initial terrain surface; then, a ground reference surface grid is constructed based on the iteratively identified ground points, and a grid-based global weighted finite difference TPS method is used to quickly generate a ground reference surface; finally, an adaptive threshold that takes into account terrain features is established to capture ground points, thereby effectively improving the point cloud filtering accuracy in complex scenes. The method is tested using the ISPRS benchmark dataset and the OpenGF large-scale dataset. The experimental results show that the average total error and kappa coefficient of the method of the present invention on the ISPRS dataset are 2.31% and 92.25% respectively, and the overall performance is better than 23 filtering methods published in the past 15 years. In addition, the processing accuracy of the OpenGF dataset by the method of the present invention is better than that of three representative filtering methods, and the accuracy is improved by at least 26.3%; at the same time, the running speed of the method of the present invention is relatively fast, reflecting the performance advantage of processing large and complex scenes. In future research, the method of the present invention can be combined with clustering technology to extract more ground seed points, thereby improving the performance of the algorithm.

[0153] The embodiments of the present invention are only used to illustrate the technical solutions of the present invention rather than to limit the present invention. It can be understood by those skilled in the art that various changes, modifications, substitutions and variations may be made to the embodiments without departing from the principles and spirit of the present invention. The scope of the present invention is defined by the attached claims and their equivalents.

Claims

1. A multi-scale hierarchical efficient filtering method taking into account terrain characteristics, characterized in that: The steps include: I. Use radius filtering and elevation histogram method to remove local and global outliers in the original point cloud; II. Based on the point cloud data after removing outliers in step I, determine the rectangular study area; further use the sliding grid technology and M-estimation sample consistency algorithm to obtain selected ground seed points from the point cloud data of the rectangular study area; III. Construct a ground reference surface based on the selected ground seed points obtained in step II. The specific steps are as follows: III.

1. Construct a grid with resolution h based on selected ground seed points to cover the entire rectangular study area; III.

2. Use the lowest selected ground seed point elevation value in each grid as the grid value; III.

3. Use global weighted finite difference TPS to interpolate the mesh surface, while reducing the impact of misclassified non-ground points and obtaining a high-quality ground reference surface; IV. using cubic spline interpolation to calculate the elevation difference between the point to be classified and the corresponding position in the ground reference surface constructed in step III, and determining the point whose elevation difference is less than a preset filtering threshold as a ground point; V. Repeat steps III.2 to IV on the ground reference plane with a resolution of h until no ground points are identified; VI. Set the resolution h = h / 2, the initial elevation threshold t0 = t0 + 0.1, and repeat steps III.1 to V until the predetermined ground reference surface resolution is reached.

2. The multi-scale hierarchical efficient filtering method taking into account terrain characteristics according to claim 1, characterized in that: The step I is specifically as follows: First, a kd-tree is constructed to organize the point cloud. Then, the number of neighboring points of each point within a specified range is searched. If the number of neighboring points of a point is less than a given threshold, the point is marked as a local outlier and removed. Then, an elevation histogram is built for the remaining points, and global outlier points below the specified elevation threshold are removed.

3. The multi-scale hierarchical efficient filtering method taking into account terrain characteristics according to claim 1 is characterized in that: The step II is specifically: First, based on the maximum value of the x-coordinate and the maximum value of the y-coordinate in the point cloud data after the outliers are removed in step I, the length and width of the rectangular research area are determined; Then, a square grid with a length and width of W is established in the rectangular area, and the square grid is slid from left to right and from bottom to top with a step length of S, and the lowest point in each square grid is selected as the initial ground seed point; Finally, the initial ground seed points are selected using the M-estimation sample consistency algorithm to obtain selected ground seed points.

4. The multi-scale hierarchical efficient filtering method taking into account terrain characteristics according to claim 3 is characterized in that: In step II, the length or width W of the square grid is greater than the minimum side length of the largest non-ground object in the rectangular study area; the step length S is the minimum side length of the largest non-ground object in the rectangular study area. times.

5. The multi-scale hierarchical efficient filtering method taking into account terrain characteristics according to claim 3 is characterized in that: In step II, in the process of selecting the lowest point in each sliding grid as the initial ground seed point, when the lowest point of the current grid overlaps with the lowest point of the previous grid, the next lowest point is selected as the initial ground seed point of the current grid.

6. The multi-scale hierarchical efficient filtering method taking into account terrain characteristics according to claim 3 is characterized in that: In step II, the specific steps of selecting the initial ground seed points using the M estimation sample consistency algorithm are as follows: II.

1. For any point P among the initial ground seed points, search for the K neighboring points closest to point P; II.

2. Randomly select three non-collinear neighbor points from the K neighbor points as the minimum sample set, and use the least squares method to calculate the fitting plane parameters corresponding to the minimum sample set; II.

3. Preset the distance threshold, calculate the distance between the remaining K-3 neighbor points and the fitting plane, and mark the neighbor points whose distance to the fitting plane is less than the preset distance threshold as in-plane points; Use the minimum cost function to iteratively obtain more in-plane points and the best fitting plane parameters; II.

4. Neighbor points whose distance to the fitting plane is greater than the preset distance threshold are marked as non-ground points and removed.

7. The multi-scale hierarchical efficient filtering method taking into account terrain characteristics according to claim 1, characterized in that: In step III.1, when setting the first level of filtering, the resolution h is taken as half of the average point spacing of the selected ground seed points, and when filtering the next level, the resolution h continues to be halved until the predetermined ground reference surface resolution is reached.

8. The multi-scale hierarchical efficient filtering method taking into account terrain characteristics according to claim 1, characterized in that: In step III.3, the specific process of using global weighted finite difference TPS to interpolate the grid surface and obtain a high-quality ground reference surface is as follows: The TPS objective function is expressed as follows using the second-order finite difference operator: In the formula, f=[f 11 ···f 1m ;f 21 ···f 2m ;···;f n1 ···f n-m ] T , f ij is the simulated value at the grid (i, j), n and m are the number of rows and columns of the grid in the study area respectively; W is a diagonal matrix diag(w i ), which contains the data Z i The corresponding weight w i ∈[0,1], where data Z i is the filling value for each grid; λ is the smoothing parameter determined by generalized cross validation; D = [D1, 2D2, D3] T , where D1, D2 and D3 are the second-order finite difference f xx 、f xy and f yy Matrix expression of ; Based on the DCT method and the finite difference strategy, the following formula can be obtained by solving formula (1): f (k+1) =IDCT(ΓDCT(W(z-f (k) )+f (k) )) (2) Where k is the number of iterations, DCT represents discrete cosine transform, IDCT represents inverse discrete cosine transform; Γ is a diagonal matrix, in which each non-zero element is a function of the smoothing parameter λ and the eigenvalue D; In the first iteration, the grid without sample values ​​is set to w = 0, and the grid with sample values ​​is set to w = 1; in subsequent iterations, w i The calculation formula is: In the formula, r i is the standardized residual, expressed as e i is the simulated residual, h i is the diagonal element of the hat matrix H, denoted as h i =h=Tr(H) / mn; σ is the robust estimate of the residual standard deviation, expressed as σ=1.843MAD, where MAD is the median of the absolute deviation.

9. The multi-scale hierarchical efficient filtering method taking into account terrain characteristics according to claim 1, characterized in that: The step IV comprises the following steps: IV.

1. Search for the 8 nearest ground points around the point to be classified and use cubic spline interpolation to calculate the interpolation elevation of the point; IV.

2. Calculate the elevation difference between the interpolated elevation of the point to be classified and the known elevation. If the elevation difference is less than the set filtering threshold, the point to be classified is marked as a ground point; if the elevation difference is greater than the set filtering threshold, the point to be classified is marked as a non-ground point; the calculation formula of the filtering threshold is: th = t0 + c × (t c +t s ); In the formula, c is the coefficient, t0, t c ,t s They are the initial elevation threshold, ground curvature value, and slope value respectively; t c The value of is obtained by calculating the curvature value of each grid cell in the ground reference surface and normalizing it; t s The value of is obtained by calculating the slope value of each grid in the ground reference surface and normalizing it.

Images