A point cloud ground filtering algorithm combining multi-scale PTD and TPS interpolation

Abstract

The application discloses a point cloud ground filtering algorithm combining multi-scale PTD and TPS interpolation, and the method is realized by the following steps: rough error points are preliminarily removed by using statistical outlier filtering; a multi-scale grid structure is constructed according to a proportional attenuation coefficient, the screening principle of original seed points is adjusted, seed points are selected layer by layer in a progressive manner in combination with a TIN filtering method, a self-adaptive boundary expansion strategy is adopted to inhibit non-ground point interference and edge information loss; a TIN is constructed in combination with seed points and extended points, and TIN filtering is performed on potential seed points in a smaller grid size of the next layer; non-ground points after TIN filtering are marked as unclassified points, a local surface is fitted by using TPS interpolation, unclassified points are classified by using the residual error between points and the constructed surface, ground points are recovered, and the real ground is iteratively approximated. The method solves the problems of excessive smoothing and terrain feature loss of traditional PTD algorithms in the edge of discontinuous terrain, steep slopes and mountain top regions.

Description

Technical Field

[0001] This invention belongs to the technical field of point cloud data ground filtering methods, specifically relating to a point cloud ground filtering algorithm that combines multi-scale PTD and TPS interpolation. Background Technology

[0002] LiDAR, as an advanced remote sensing technology, has been widely used in various fields such as topographic mapping, 3D modeling, and disaster assessment due to its high precision, high efficiency, and all-weather operation capabilities. LiDAR scanning can acquire a large amount of 3D point cloud data, reflecting detailed information about the Earth's surface. However, raw LiDAR point cloud data often contains a large amount of noise and non-ground points (such as buildings and vegetation). Accurately distinguishing between ground points and non-ground points and effectively filtering them is one of the key technical challenges in LiDAR data processing.

[0003] Existing point cloud ground filtering methods can be mainly divided into five categories: slope analysis-based methods, surface interpolation methods, digital morphology methods, iterative densification methods, and clustering-based segmentation methods. Iterative densification algorithms, such as the Progressive Triangulated Irregular Network (PTD) method, have proven effective in processing large-scale point cloud data. They iteratively construct ground seed points and gradually refine the ground model using triangular meshes, effectively adapting to complex terrain and irregular surfaces. However, the PTD algorithm still faces the following problems: low computational efficiency, insufficient detail preservation in complex terrain, poor performance on low-lying points (such as ground points under vegetation cover), and a tendency to misclassify. Summary of the Invention

[0004] The purpose of this invention is to provide a point cloud ground filtering algorithm that combines multi-scale PTD and TPS interpolation, which can improve the extraction accuracy of ground points and retain more terrain details when dealing with large-scale complex terrain, thereby effectively overcoming the limitations of existing algorithms.

[0005] The technical solution adopted in this invention is a point cloud ground filtering algorithm that combines multi-scale PTD and TPS interpolation, comprising the following steps:

[0006] Step S1: Extract the point cloud of the last echo of the lidar and use statistical outlier filtering to initially remove gross errors;

[0007] Step S2: Construct a multi-scale grid structure according to the proportional attenuation coefficient, screen out potential ground seed points and expand the boundary;

[0008] Step S3: Construct a TIN by combining seed points and extension points, and perform TIN filtering on potential seed points in the next layer with smaller grid size;

[0009] Step S4: Use the ground points after TIN filtering as seed points to test the next layer. Fit the local surface for non-ground points based on TPS interpolation. Update the points smaller than the threshold as seed points and mark the rest as non-ground points.

[0010] Step S5: If the top layer is not reached, return to step S3; otherwise, output the final ground point set and the algorithm terminates.

[0011] Furthermore, the specific steps of data preprocessing in step S1 above are as follows:

[0012] Step S11: Retrieve the k nearest neighbors of each point using a KD-tree-based spatial index. For each point, calculate the average distance d within its neighborhood. i And assume that the distance between neighboring points follows a Gaussian distribution N(μ,σ);

[0013] Step S12, with μ+d std ×σ is the threshold value. In UAV / vehicle-mounted LiDAR scenarios, the noise ratio is 0.5%~2%. Experiments show that setting d... std =2.0~2.5 Balanced noise removal and terrain preservation.

[0014] Furthermore, the specific steps for selecting ground seed points in step S2 above are as follows:

[0015] Step S21, Low-position noise point interference suppression: Based on local minimum constraints, the lowest point in the grid is selected as the initial ground point. The seed point verification method is adopted. By comparing the elevation difference between the lowest point and the second lowest point, if the difference is less than the preset threshold Δh, the lowest point is considered to be a reliable seed point. Otherwise, the second lowest point is a point to be confirmed. The same method is used to compare with the third lowest point until the condition is met.

[0016] Step S22, Hierarchical Structure Construction: In the lowest level L max (corresponding to the maximum grid size W) max In this process, only high-confidence ground points are selected as initial seeds, and the maximum grid size W is used. max The size of the topmost mesh is determined by the largest non-ground object in the scene to avoid misselection of ground object edges; min Related to point density; the mesh size W of the intermediate layer k Determined by the proportional attenuation coefficient s; using a bottom-up approach, each layer corresponds to an initial grid of seed points with different resolutions, and the seed point increment between adjacent layers is s = 2 / 3;

[0017] Step S23, Distance Threshold Adjustment: Introduce the maximum distance S max and minimum distance S min Two parameters control the thresholds for the top and bottom layers, respectively, and the distance threshold S for the intermediate layers. kAdjustments are made according to a linear decreasing rule, and the calculation formula is as follows:

[0018] ;

[0019] In the formula, S k S is the distance threshold corresponding to the k-th layer. max and S min These are the specified maximum and minimum distance thresholds, L. max and L min These are the maximum and minimum grid sizes, L. k This represents the grid size corresponding to the k-th layer;

[0020] Step S24, Distribution of extended points: Extend a new ring of seed points around the existing seed points to ensure that all points are included within the outer rectangle. The extended points are evenly distributed at equal intervals according to the grid size, and the elevation is calculated based on the non-ground point interference that may be included in the potential seed points.

[0021] Calculation of extension point elevation:

[0022] 1) Let the coordinates of the target center be... Search for valid grids containing seed points within the eight neighborhoods; the corresponding seed point set is: The filtering criteria are as follows:

[0023] ;

[0024] In the formula, For the filtered point cloud set, The target center point The set of neighborhood points, The average distance between neighboring points and the target center;

[0025] 2) Dense region detection and interpolation: For the filtered high-order sequences Sort in ascending order and perform ANOVA based on a sliding window:

[0026] ;

[0027] ;

[0028] In the formula, and For high program columns middle Elevation value The calculated mean and variance of elevation;

[0029] Retain satisfaction contiguous window set Finally, it was calculated that in a certain Interpolated elevation at point for ,in Given an elevation variance threshold;

[0030] like Then take the global minimum value of the elevation of that layer. ,in This is the set of elevation values ​​for this layer.

[0031] Furthermore, the specific steps of TIN filtering in step S3 above are as follows:

[0032] S31. Construct an initial irregular triangular mesh using selected seed points. This mesh consists of triangles, with each vertex of a triangle serving as a seed point. Calculate the triangle plane formula based on its coordinates. ;

[0033] S32. For each unclassified point p, calculate its distance to the TIN surface:

[0034] ;

[0035] In the formula, A, B, C, and D are the coefficients obtained by solving the TIN surface equation;

[0036] S33. Calculate the angle between the point and the line connecting the vertex of the TIN face. If both the angle and the distance are less than the preset threshold, the point is classified as a ground point.

[0037] ;

[0038] In the formula, For the first Point and TIN face vertex The angle between the line connecting it and its perpendicular projection line on the TIN plane. For point The vertical distance to the corresponding TIN plane. , and For the first The coordinates of the point , and As vertices The coordinate values.

[0039] Furthermore, the specific steps for complex terrain restoration in step S4 above are as follows:

[0040] Step S41: Divide the non-ground points obtained by each layer of TIN filtering into unclassified points;

[0041] Step S42: In mountainous and hilly areas, based on known ground points, local TPS interpolation is used to reconstruct the terrain surface, and the residuals between the points and the constructed surface are used to classify the unclassified points.

[0042] 1) Search for ground points within a 5×5 area centered on non-ground points. Perform surface fitting when at least 8 ground points are found.

[0043] 2) The TPS interpolation method obtains a smooth interpolation function by minimizing the energy function. This method combines polynomial terms and radial basis functions to fit given discrete data points; assuming that... Known data points ,in The interpolated surface can be represented as:

[0044] ;

[0045] In the formula, It is a polynomial trend function with degree k, used to capture the overall trend of terrain change. basis functions The weighting coefficients, It is a radial basis function. It is the Euclidean distance between the point to be interpolated and the known point;

[0046] 3) To ensure the interpolation function A square-integrable second derivative must satisfy the following constraints:

[0047] , ;

[0048] A regularization term is introduced, which weights the diagonal of the kernel matrix to balance the influence of local errors;

[0049] Construct a system of linear equations to solve for the weighting coefficients and ,

[0050] ;

[0051] Where K is the kernel matrix calculated from the distances between known data points. is the regularization coefficient, P is a matrix generated by a polynomial function, and z is the elevation vector of a known point.

[0052] 4) Calculate the difference di between the true value of each ground point and the fitted elevation of its corresponding surface, and obtain the mean value d. mean Based on the standard error σ and the characteristic that a large number of discrete ground points follow a normal distribution under natural conditions, the elevation difference threshold for each point is determined as d. mean +3σ.

[0053] Compared with the prior art, the improvements of this invention are mainly concentrated in three aspects:

[0054] 1) Adjust the original seed point selection criteria and combine the TIN filtering method to select seed points layer by layer;

[0055] 2) An adaptive boundary expansion strategy is adopted to suppress non-ground point interference and edge information loss;

[0056] 3) For points lost in complex terrain, TPS interpolation is used to fit the local surface, recover the ground points, and iteratively approximate the real ground.

[0057] The improved PTD algorithm proposed in this invention combines three main aspects: multi-scale seed point selection, adaptive boundary expansion, and TPS differential completion. Validation on the ISPRS standard dataset shows that the proposed method significantly outperforms the classic PTD and other improved PTD algorithms in terms of total error and Kappa coefficient. The spatial distribution of seed points is optimized through the proportional attenuation coefficient of the hierarchical grid and the maximum window constraint. Compared to the traditional quadtree structure, this strategy suppresses the intrusion of low-elevation noise at the initial level and gradually refines the terrain during iterative encryption, reducing the Type I error to 2.28%, a performance improvement of 53.7% compared to the classic PTD. The use of extended grid interpolation and neighborhood variance constraints avoids the misjudgment of narrow triangles caused by corner point insertion in traditional methods. Adaptive boundary expansion enhances the protection of edge information. Finally, secondary filtering recovers missed ground points through TPS interpolation fitting, significantly improving the over-smoothing and feature loss problems of the classic PTD in complex terrain.

[0058] This invention significantly improves the accuracy and robustness of ground point extraction through innovative methods such as multi-scale grids, TPS interpolation, adaptive boundary expansion, and dynamic threshold control. In particular, when dealing with complex terrain, it can more accurately recover ground features, reduce misclassification, and improve the overall performance of the algorithm. Attached Figure Description

[0059] Figure 1 This is a flowchart of the point cloud ground filtering algorithm combining multi-scale PTD and TPS interpolation as described in this invention;

[0060] Figure 2 This is a flowchart of the specific steps in step S1 of the method described in this invention;

[0061] Figure 3 This is a flowchart of the specific steps in step S2 of the method described in this invention;

[0062] Figure 4 This is a flowchart of the specific steps in step S3 of the method described in this invention;

[0063] Figure 5 This is a flowchart of the specific steps in step S4 of the method described in this invention. Detailed Implementation

[0064] The present invention will be further explained and described below with reference to the accompanying drawings to enable those skilled in the art to better understand it.

[0065] Example 1

[0066] A point cloud ground filtering algorithm that combines multi-scale PTD and TPS interpolation, such as Figure 1 As shown, it includes:

[0067] Step S1: Extract the point cloud of the last echo of the lidar and use statistical outlier filtering to initially remove gross errors, and preprocess the data;

[0068] Specifically, such as Figure 2 As shown, the specific method for data preprocessing is as follows:

[0069] Step S11: Efficiently retrieve the k nearest neighbors of each point using a KD-tree-based spatial index. For each point, calculate the average distance d within its neighborhood. i And assume that the distance between neighboring points follows a Gaussian distribution N(μ,σ);

[0070] Step S12, with μ+d std ×σ is the threshold value. In drone / vehicle-mounted LiDAR scenarios, the noise ratio is typically 0.5%~2%. Setting d... std =2.0~2.5 can balance noise removal and terrain preservation.

[0071] Step S2: Construct a multi-scale grid structure according to the proportional attenuation coefficient, screen out potential ground seed points and expand the boundary;

[0072] Specifically, such as Figure 3 As shown, the specific steps for selecting ground seed points are as follows:

[0073] Step S21, Low-position noise point interference suppression: Based on local minimum constraints, the lowest point in the grid is selected as the initial ground point. If the preprocessing fails to completely eliminate low-position noise points, these points are significantly more harmful than high-position noise points and are easily misjudged as terrain depression features. To suppress low-position noise point interference, a seed point verification method is adopted. By comparing the elevation difference between the lowest point and the second lowest point, if the difference is less than the preset threshold Δh, the lowest point is considered a reliable seed point; otherwise, the second lowest point is a point to be confirmed. The same method is used to compare with the third lowest point until the conditions are met.

[0074] Step S22, Hierarchical Structure Construction: In the lowest level L max (corresponding to the maximum grid size W) maxIn this process, only high-confidence ground points are selected as initial seeds, and the maximum grid size W is used. max The size of the topmost mesh is determined by the largest non-ground object in the scene to avoid misselection of ground object edges; min Related to point density; the mesh size W of the intermediate layer k The ratio attenuation coefficient s determines the initial grid of seed points for each layer, which is generated from the bottom up. In the traditional quadtree structure, each layer is generated by sampling down from the previous layer by a factor of 2, while this method controls the increase of the seed points between adjacent layers to 1.5 times, i.e., s=2 / 3.

[0075] Step S23, Distance Threshold Adjustment: Introduce the maximum distance S max and minimum distance S min Two parameters control the thresholds for the top and bottom layers respectively; the distance threshold S for intermediate layers... K Adjustments are made according to a linear decreasing rule, and the calculation formula is as follows:

[0076] ;

[0077] In the formula, S k S is the distance threshold corresponding to the k-th layer. max and S min These are the specified maximum and minimum distance thresholds, L. max and L min These are the maximum and minimum grid sizes, L. k This represents the grid size corresponding to the k-th layer;

[0078] Step S24, Distribution of Extended Points: To avoid misjudgment and error accumulation caused by the narrow triangles generated by the boundary, a new ring of seed points is extended outside the existing seed points to ensure that all points are included within the outer rectangle; the extended points are evenly distributed at equal intervals according to the grid size, and the elevation is calculated based on the non-ground point interference that may be included in the potential seed points.

[0079] Calculation of extension point elevation:

[0080] 1) Let the coordinates of the target center be... Search for valid grids containing seed points within the eight neighborhoods; the corresponding seed point set is: The filtering criteria are as follows:

[0081] ;

[0082] In the formula, For the filtered point cloud set, The target center point The set of neighborhood points, The average distance between neighboring points and the target center;

[0083] 2) Dense region detection and interpolation: For the filtered high-order sequences Sort in ascending order and perform ANOVA based on a sliding window:

[0084] ;

[0085] ;

[0086] In the formula, and For high program columns middle Elevation value The calculated mean and variance of elevation;

[0087] Retain satisfaction contiguous window set Finally, it was calculated that in a certain Interpolated elevation at point for ,in Given an elevation variance threshold;

[0088] like Then take the global minimum value of the elevation of that layer. ,in This is the set of elevation values ​​for this layer.

[0089] Step S3: Construct a TIN by combining seed points and extension points, and perform TIN filtering on potential seed points in the next layer with smaller grid size;

[0090] Specifically, such as Figure 4 As shown, the specific steps of TIN filtering are as follows:

[0091] Step S31: Construct an initial irregular triangular mesh using the selected seed points. This mesh consists of triangles, with each vertex of a triangle serving as a seed point. Calculate the triangle plane formula based on its coordinates. ;

[0092] Step S32: For each unclassified point p, calculate its distance to the TIN surface:

[0093] ;

[0094] In the formula, A, B, C, and D are the coefficients obtained by solving the TIN surface equation;

[0095] Step S33: Calculate the angle between the point and the line connecting the vertex of the TIN face. If both the angle and the distance are less than the preset threshold, the point is classified as a ground point.

[0096] ;

[0097] In the formula, For the first Point and TIN face vertex The angle between the line connecting it and its perpendicular projection line on the TIN plane. For point The vertical distance to the corresponding TIN plane. , and For the first The coordinates of the point , and As vertices The coordinate values.

[0098] Step S4: Use the ground points after TIN filtering as seed points to test the next layer. Fit the local surface for non-ground points based on TPS interpolation. Update the points smaller than the threshold as seed points and mark the rest as non-ground points.

[0099] Specifically, such as Figure 5 As shown, the specific steps for complex terrain restoration are as follows:

[0100] Step S41: Divide the non-ground points obtained by each layer of TIN filtering into unclassified points;

[0101] Step S42: In mountainous and hilly areas, some unclassified points located in convex areas of the terrain may be misclassified as non-ground points due to the inherent defects of PTD filtering because their angle with the TIN triangle exceeds the maximum angle threshold. Similarly, points located in depressions at the foot of mountains may also lose effective points because their angle exceeds the threshold. To address this, a ground point recovery method is proposed: based on known ground points, local TPS interpolation is used to reconstruct the terrain surface, and the residual between the points and the constructed surface is used to classify the unclassified points.

[0102] 1) The computational complexity of TPS is directly related to the number of ground points involved in the calculation, so local interpolation is used to improve efficiency; ground points within a 5×5 range are retrieved with non-ground points as the center, and surface fitting is performed when no less than 8 ground points are retrieved.

[0103] 2) The TPS interpolation method obtains a smooth interpolation function by minimizing the energy function; this method combines polynomial terms and radial basis functions to fit given discrete data points, assuming that... Known data points ,in The interpolated surface can be represented as:

[0104] ;

[0105] In the formula, It is a polynomial trend function with degree k, used to capture the overall trend of terrain change. basis functions The weighting coefficients, It is a radial basis function. It is the Euclidean distance between the point to be interpolated and the known point;

[0106] 3) To ensure the interpolation function A square-integrable second derivative must satisfy the following constraints:

[0107] , ;

[0108] In TPS interpolation, a regularization term is usually introduced to improve stability and prevent overfitting. This regularization term weights the diagonal of the kernel matrix, thereby balancing the influence of local errors. In the process of ground point classification and restoration, the regularization term controls the smoothness of the interpolation surface, avoiding unnatural fluctuations in details.

[0109] Construct a system of linear equations to solve for the weighting coefficients and ,

[0110] ;

[0111] In the formula, K is the kernel matrix calculated from the distances between known data points. is the regularization coefficient, P is a matrix generated by a polynomial function, and z is the elevation vector of a known point.

[0112] 4) Calculate the difference di between the true value of each ground point and the fitted elevation of its corresponding surface, and obtain the mean value d. mean Based on the standard error σ, and considering the characteristic that a large number of discrete ground points follow a normal distribution under natural conditions, while non-ground points are error points with a standard deviation greater than 3σ, the elevation difference threshold for each point is determined as d. mean +3σ.

[0113] Step S5: If the top layer is not reached, return to step 3; otherwise, output the final ground point set and the algorithm terminates.

[0114] To illustrate the filtering effect of the point cloud filtering algorithm provided in this embodiment of the invention, a comparative experiment was conducted on 15 sets of reference data published by the International Society for Photogrammetry and Remote Sensing, comparing the method of this invention with five existing PTD correlation filtering algorithms.

[0115] Table 1 shows the calculation process of the quantitative indicators for evaluating the filtering effect of this invention. Where TI represents Type I error, T.II represents Type II error, TE represents total error, and κ represents the Kappa coefficient. The smaller the values ​​of Type I error, Type II error, and total error, the smaller the error caused by filtering. The larger the Kappa coefficient value, the higher the classification consistency.

[0116] Table 1

[0117]

[0118] After systematic evaluation on 15 benchmark datasets, as shown in Table 2, the optimized parameter configuration achieved an average total error of 3.06% and an average Kappa coefficient of 88.9%. Notably, samples S21, S31, and S42 (all located in flat terrain) performed particularly well, with Kappa coefficients exceeding 95% and total errors below 2%. This aligns with existing research findings that the filtering algorithm performs excellently in flat areas, but its reliability decreases in steep or complex terrain.

[0119] In contrast, samples S11, S52, and S53 have lower Kappa coefficients (<85%). These regions are characterized by significant slopes, abrupt elevation changes, and low vegetation on the slopes, making it difficult to remove many low outliers. Furthermore, S52 and S53 suffer from severe class imbalance, with non-ground points accounting for only 10.51% and 4.04% of the dataset, respectively (as shown in Table 3). Due to the low fault tolerance, even a small number of misclassified non-ground points significantly increase Class II errors. The Kappa coefficients of the remaining samples remain stable (85-95%), demonstrating the adaptability of this invention in different scenarios. Our invention overcomes the inherent limitations of traditional PTD algorithms through a two-stage filtering strategy. The first stage divides the original point cloud into a multi-layered structure, with point cloud filtering at different resolutions associated with different distance parameters to identify reliable ground seed points, effectively reducing TI and enhancing the algorithm's robustness in heterogeneous environments. The second stage adaptively reclassifies non-ground points, thereby restoring the terrain with elevation discontinuities.

[0120] Table 2 Parameter settings and error distribution

[0121]

[0122] This invention was compared with previous filters using six PTD-related algorithms developed between 2000 and 20223 and the commercial software TerraSolid. Performance was evaluated based on TE, as shown in Table 3. In 15 test samples, this method had the lowest mean total error and achieved the best accuracy in 9 samples. Compared to classic PTD, this method achieved better results in almost all scenarios. It provided approximately 36.5% performance improvement in TE and a 4.66% increase in the Kappa coefficient.

[0123] Table 3 Performance evaluation of PTD-related filtering methods

[0124]

[0125] The improved PTD algorithm proposed in this invention combines three main aspects: multi-scale seed point selection, adaptive boundary expansion, and TPS differential completion. Validation based on the ISPRS standard dataset shows that the proposed method significantly outperforms the classic PTD and other improved PTD algorithms in terms of total error and Kappa coefficient. This paper optimizes the spatial distribution of seed points through the proportional attenuation coefficient of the hierarchical grid and the maximum window constraint. Compared to the traditional quadtree structure, this strategy suppresses the intrusion of low-elevation noise at the initial level and gradually refines the terrain during iterative encryption, reducing the Type I error to 2.28%, achieving a performance improvement of 53.7% compared to the classic PTD. The use of extended grid interpolation and neighborhood variance constraints avoids the misjudgment of narrow triangles caused by corner point insertion in traditional methods. Adaptive boundary expansion enhances the protection of edge information. Finally, secondary filtering recovers missed ground points through TPS interpolation fitting, significantly improving the over-smoothing and feature loss problems of the classic PTD in complex terrain.

[0126] The embodiments described above are merely preferred embodiments of the present invention and are not intended to limit the scope of the present invention. Various modifications and improvements made by those skilled in the art to the technical solutions of the present invention without departing from the spirit and principles thereof should fall within the protection scope defined by the claims of the present invention.

Claims

1. A point cloud ground filtering algorithm combining multi-scale PTD and TPS interpolation, characterized in that, Includes the following steps: Step S1: Extract the point cloud of the last echo of the lidar and use statistical outlier filtering to initially remove gross errors; Step S2: Construct a multi-scale grid structure according to the proportional attenuation coefficient, screen out potential ground seed points and expand the boundary; Step S3: Construct a TIN by combining seed points and extension points, and perform TIN filtering on potential seed points in the next layer with smaller grid size; Step S4: Use the ground points after TIN filtering as seed points to test the next layer. Fit the local surface for non-ground points based on TPS interpolation. Update the points smaller than the threshold as seed points and mark the rest as non-ground points. Step S5: If the top layer has not been reached, return to step S3; otherwise, output the final ground point set and the algorithm terminates. The specific steps for selecting ground seed points in step S2 are as follows: Step S21, Low-position noise point interference suppression: Based on local minimum constraints, the lowest point in the grid is selected as the initial ground point. The seed point verification method is adopted. By comparing the elevation difference between the lowest point and the second lowest point, if the difference is less than the preset threshold Δh, the lowest point is considered to be a reliable seed point. Otherwise, the second lowest point is a point to be confirmed. The same method is used to compare with the third lowest point until the condition is met. Step S22, Hierarchical Structure Construction: In the lowest level L max In the middle, only high-confidence ground points are selected as initial seeds, and the lowest layer L... max Corresponding maximum grid size W max Maximum grid size W max The size of the topmost mesh is determined by the largest non-ground object in the scene to avoid misselection of ground object edges; min Related to point density; the mesh size W of the intermediate layer k Determined by the proportional attenuation coefficient s; using a bottom-up approach, each layer corresponds to an initial grid of seed points with different resolutions, and the seed point increment between adjacent layers is s = 2 / 3; Step S23, Distance Threshold Adjustment: Introduce the maximum distance S max and minimum distance S min Two parameters control the thresholds for the top and bottom layers respectively. The distance thresholds for intermediate layers are adjusted according to a linear decreasing rule. The calculation formula is as follows: ； In the formula, S k S is the distance threshold corresponding to the k-th layer. max and S min These are the specified maximum and minimum distance thresholds, L. max and L min These are the maximum and minimum grid sizes, L. k This represents the grid size corresponding to the k-th layer; Step S24, Distribution of extended points: Extend a new ring of seed points around the existing seed points to ensure that all points are included within the outer rectangle. The extended points are evenly distributed at equal intervals according to the grid size, and the elevation is calculated based on the non-ground point interference that may be included in the potential seed points. Calculation of extension point elevation: 1) Let the coordinates of the target center be... Search for valid grids containing seed points within the eight neighborhoods; the corresponding seed point set is: The filtering criteria are as follows: ； In the formula, For the filtered point cloud set, The target center point The set of neighborhood points, The average distance between neighboring points and the target center; 2) Dense region detection and interpolation: For the filtered high-order sequences Sort in ascending order and perform ANOVA based on a sliding window: ； ； In the formula, and For high program columns middle Elevation value The calculated mean and variance of elevation; Retain satisfaction contiguous window set Finally, it was calculated that in a certain Interpolated elevation at point for ,in Given an elevation variance threshold; like Then take the global minimum value of the elevation of that layer. ,in This is the set of elevation values ​​for this layer.

2. The point cloud ground filtering algorithm combining multi-scale PTD and TPS interpolation as described in claim 1, characterized in that, The specific steps of data preprocessing in step S1 are as follows: Step S11: Retrieve the k nearest neighbors of each point using a KD-tree-based spatial index. For each point, calculate the average distance d within its neighborhood. i And assume that the distance between neighboring points follows a Gaussian distribution N(μ,σ); Step S12, with μ+d std ×σ is the threshold value. In UAV / vehicle-mounted LiDAR scenarios, the noise ratio is 0.5%~2%. Experiments show that setting d... std =2.0~2.5 Balanced noise removal and terrain preservation.

3. The point cloud ground filtering algorithm combining multi-scale PTD and TPS interpolation as described in claim 1, characterized in that, The specific steps of TIN filtering in step S3 are as follows: S31. Construct an initial irregular triangular mesh using selected seed points. This mesh consists of triangles, with each vertex of a triangle serving as a seed point. Calculate the triangle plane formula based on its coordinates. ; S32. For each unclassified point p, calculate its distance to the TIN surface: ； In the formula, A, B, C, and D are the coefficients obtained by solving the TIN surface equation; S33. Calculate the angle between the point and the line connecting the vertex of the TIN face. If both the angle and the distance are less than the preset threshold, the point is classified as a ground point. ； In the formula, For the first Point and TIN face vertex The angle between the line connecting it and its perpendicular projection line on the TIN plane. For point The vertical distance to the corresponding TIN plane. , and For the first The coordinates of the point , and As vertices The coordinate values.

4. The point cloud ground filtering algorithm combining multi-scale PTD and TPS interpolation as described in claim 1, characterized in that, The specific steps for complex terrain restoration in step S4 are as follows: Step S41: Divide the non-ground points obtained by each layer of TIN filtering into unclassified points; Step S42: In mountainous and hilly areas, based on known ground points, local TPS interpolation is used to reconstruct the terrain surface, and the residuals between the points and the constructed surface are used to classify the unclassified points. 1) Search for ground points within a 5×5 area centered on non-ground points. Perform surface fitting when at least 8 ground points are found. 2) The TPS interpolation method obtains a smooth interpolation function by minimizing the energy function. This method combines polynomial terms and radial basis functions to fit given discrete data points; assuming that... Known data points ,in The interpolated surface can be represented as: ； In the formula, It is a polynomial trend function with degree k, used to capture the overall trend of terrain change. basis functions The weighting coefficients, It is a radial basis function. It is the Euclidean distance between the point to be interpolated and the known point; 3) To ensure the interpolation function A square-integrable second derivative must satisfy the following constraints: , ； A regularization term is introduced, which weights the diagonal of the kernel matrix to balance the influence of local errors; Construct a system of linear equations to solve for the weighting coefficients and , ； Where K is the kernel matrix calculated from the distances between known data points. is the regularization coefficient, P is a matrix generated by a polynomial function, and z is the elevation vector of a known point. 4) Calculate the difference di between the true value of each ground point and the fitted elevation of its corresponding surface, and obtain the mean value d. mean Based on the standard error σ and the characteristic that a large number of discrete ground points follow a normal distribution under natural conditions, the elevation difference threshold for each point is determined as d. mean +3σ.

Images