A Multi-Scale PTD Filtering Method Based on Terrain Compensation

Abstract

This invention discloses a multi-scale PTD filtering method based on terrain compensation. The method comprises the following steps: extracting the last echo point cloud from a lidar radar and initially removing gross errors using statistical outlier filtering; constructing a multi-scale grid structure according to a proportional attenuation coefficient, adjusting the original seed point selection principle, and progressively selecting seed points layer by layer using a TIN filtering method, employing an adaptive boundary expansion strategy to suppress non-ground point interference and edge information loss; constructing a TIN by combining seed points and extended points, and performing TIN filtering on potential seed points in the next smaller grid size; marking the non-ground points after TIN filtering as unclassified points, performing secondary filtering based on terrain specificity, and iteratively approximating the real ground; this method solves the problems of excessive smoothing and terrain feature loss in traditional PTD algorithms at discontinuous terrain edges, steep slopes, and mountain tops.

Description

Technical Field

[0001] This invention belongs to the field of topographic mapping technology, specifically relating to a multi-scale PTD filtering method based on topographic compensation. Background Technology

[0002] As a core data acquisition tool for active remote sensing, lidar technology has established an irreplaceable position in fields such as topographic mapping, disaster assessment, and 3D modeling due to its all-weather operation capabilities, multi-level vegetation penetration characteristics, and millimeter-level 3D coordinate resolution. However, the mixed noise interference and multi-source ground cover (such as buildings and vegetation) in the raw point cloud data make the separation of ground points from non-ground points (i.e., ground filtering) a crucial and challenging aspect of lidar data applications. Existing filtering methods can be systematically categorized into five types: slope analysis-based, surface interpolation, iterative encryption, digital morphology, and clustering segmentation methods. Among them, the progressive encryption irregular triangular mesh algorithm, as a typical representative of iterative encryption algorithms, constructs the terrain surface by iteratively encryption of initial ground seed points. It is widely adopted due to its strong adaptability to terrain continuity and its ability to utilize more contextual information, and has been integrated into the commercial software TerraSolid. The main steps are as follows: assign points to a regular grid, select the lowest value of the grid as the initial ground point, construct a TIN with these points, calculate the distance s between the unclassified point and the TIN surface and the angle θ between the point and the vertex of the TIN surface. If the threshold condition is met, update it as a ground point and add it to the TIN. The TIN is continuously iterated and encrypted in this process until all points are classified.

[0003] Although PTD achieved the best performance among the eight filtering algorithms in the Sithole and Vosselman comparison, it still has certain limitations. For example, the iterative process of the PTD algorithm leads to low computational efficiency; it is difficult to retain fine terrain details in complex terrain; and it is often unable to effectively remove low points attached to the ground, such as points under vegetation cover, leading to misclassification problems.

[0004] In conclusion, the progressive irregular triangular mesh encryption filtering algorithm needs further improvement. Summary of the Invention

[0005] The purpose of this invention is to provide a multi-scale PTD filtering method based on terrain compensation, which can obtain a more accurate seed point set, preserve the terrain details of complex terrain areas to the maximum extent, and improve computational efficiency and the accuracy of ground point extraction.

[0006] The technical solution adopted in this invention is a multi-scale PTD filtering method based on terrain compensation, comprising the following steps:

[0007] 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;

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

[0009] 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;

[0010] Step S4: Use the ground points after TIN filtering as seed points to test the next layer. Perform secondary filtering on non-ground points based on terrain specificity. Update the points that pass the compensation judgment as seed points, and mark the rest as non-ground points.

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

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

[0013] Step S11: Retain the point cloud of the last echo and remove non-surface signals from the intermediate echoes;

[0014] Step S12: 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(μ,σ), with μ+d std ×σ is the threshold value. In drone / vehicle-mounted LiDAR scenarios, the noise ratio is 0.5%~2%. d is set... std =2.0~2.5 Balanced noise removal and terrain preservation.

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

[0016] 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.

[0017] 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 W of the topmost mesh is determined by the largest non-ground object in the scene. minRelated to the point density, the mesh size W of the intermediate layer k Determined by the proportional attenuation coefficient s, the increase factor of the seed points between adjacent layers is s = 2 / 3.

[0018] 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, and the calculation formula is as follows:

[0019] ;

[0020] 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;

[0021] Step S24, Extension Point Distribution: Extend a new ring of seed points around the existing seed points, ensuring that all points are encompassed within the outer rectangle; the extended points are evenly distributed at equal intervals according to the grid size. Elevation calculations are performed based on potential non-ground point interference from the potential seed points, using the following method:

[0022] 1) Target center coordinates and neighborhood filtering: Set the target center coordinates as... The corresponding point cloud set is P, and 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 analyze variance using a sliding window:

[0026] , ;

[0027] In the formula, and For high program columns middle Elevation value Calculate the mean and variance of elevations; retain those that satisfy the following conditions. contiguous window set Finally, it was calculated that in a certain Interpolated elevation at point :

[0028] ;

[0029] like Then take the global minimum value of the elevation of that layer. ,in Let be the set of elevation values ​​for the k-th layer.

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

[0031] 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. ;

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

[0033] ;

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

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

[0036] ;

[0037] 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 vertex The coordinate values.

[0038] Furthermore, the specific steps for terrain-specific compensation in step S4 are as follows:

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

[0040] Step S42: Use the least squares method to fit the local surface to estimate ground points, and classify unclassified points using a dynamic elevation difference threshold.

[0041] 1) Using the unclassified point as the center, retrieve at least 6 ground points within a 5×5 area and perform surface fitting;

[0042] 2) The equation of the two-dimensional polynomial surface is:

[0043] ;

[0044] Represented as a matrix equation:

[0045] ;

[0046] in:

[0047] , , ;

[0048] When m≥6 and When it is invertible, the optimal parameters need to be solved using the least squares method:

[0049] ;

[0050] 3) Calculate the fitted elevation difference d between each ground point and its corresponding fitted surface equation. i The mean d is obtained. mean Based on the standard error σ, the height difference threshold d is determined for each point. mean +3σ;

[0051] Step 43: Compensation algorithm for slopes and faults:

[0052] 1) Calculate the slope of each point relative to P by statistically analyzing the ground points in the eight neighboring areas around the unclassified point P;

[0053] 2) Set the slope threshold θ r = 10°, used as the standard to distinguish between slopes and flat areas;

[0054] 3) Count the number of ground points n that satisfy θ < θᵣ flat And the number n of ground points satisfying θ ≥ θᵣ slope The formulas for fault and slope compensation are as follows:

[0055] .

[0056] The improvements of this invention mainly focus on three aspects:

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

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

[0059] 3) Accurate compensation of terrain features is achieved through secondary filtering, iteratively approximating the real ground.

[0060] This invention proposes a multi-scale PTD filtering algorithm based on terrain compensation, integrating three main aspects: multi-scale seed point selection, adaptive boundary expansion, and terrain-specific compensation. Validation based on the ISPRS standard dataset shows that this method significantly outperforms the classic PTD and other improved PTD algorithms in terms of total error (2.96%) and Kappa coefficient (89.36%). This invention optimizes the spatial distribution of seed points through the proportional attenuation coefficient of the layered 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 densification, reducing the Type I error (TI) to 2.06%, achieving a 70.5% performance improvement over the classic PTD. The use of extended grid interpolation and neighborhood variance constraints effectively solves the problem of boundary omission in fault and steep slope areas (such as S11 and S41), avoiding the misjudgment of narrow triangles caused by corner point insertion in traditional methods. Adaptive boundary expansion enhances the protection of edge information. Finally, the secondary filtering (surface fitting and slope filtering) recovers the missed ground points, significantly improving the problems of over-smoothing and feature loss in complex terrains caused by the classic PTD.

[0061] This invention combines PTD with surface fitting and slope filtering. Through innovative methods such as multi-scale grids, terrain compensation, adaptive boundary expansion, and dynamic threshold control, it obtains a more accurate seed point set and preserves terrain details in complex areas to the maximum extent. This significantly improves the accuracy and robustness of ground point extraction. 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

[0062] Figure 1 This is a flowchart of the steps of the multi-scale PTD filtering algorithm based on terrain compensation described in this invention;

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

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

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

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

[0067] 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.

[0068] Example 1

[0069] A multi-scale PTD filtering method based on terrain compensation, such as Figure 1 As shown, it includes:

[0070] 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;

[0071] Specifically, such as Figure 2 As shown, the specific steps for data preprocessing include:

[0072] Step S11: Retain the point cloud of the last echo: The lidar echo sequence contains multiple reflection signals, among which the last echo becomes the key data source for surface elevation modeling because it reaches the ground surface; therefore, firstly, only the point cloud data of the last echo is retained, and non-surface signals in the intermediate echoes are removed.

[0073] Step S12: Efficiently retrieve the k nearest neighbors of each point using a spatial index based on a KD tree; for each point, calculate the average distance di within its neighborhood, and assume that the distance between neighboring points follows a Gaussian distribution N(μ,σ). Using μ+dstd×σ as the threshold, in UAV / vehicle-mounted LiDAR scenarios, the noise ratio is usually 0.5%~2%. Experiments show that setting dstd =2.0~2.5 can balance noise removal and terrain preservation.

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

[0075] Specifically, the steps for selecting ground seed points in step S2 are as follows:

[0076] 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 it with the third lowest point until the condition is met.

[0077] 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 The ratio of the attenuation coefficient s is determined by the ratio of the attenuation coefficient s; the number of seed points between adjacent layers increases by a factor of 1.5 (i.e., s = 2 / 3).

[0078] 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:

[0079] ;

[0080] 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;

[0081] 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.

[0082] Calculation of extension point elevation:

[0083] 1) Target center coordinates and neighborhood filtering: Set the target center coordinates as...

[0084] The corresponding point cloud set is P, and the filtering criteria are as follows:

[0085] ;

[0086] 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;

[0087] 2) Dense region detection and interpolation: For the filtered high-order sequences Sort in ascending order and analyze variance using a sliding window:

[0088]

[0089] ;

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

[0091] Retain satisfaction contiguous window set Finally, it was calculated that in a certain Interpolated elevation at point :

[0092] ;

[0093] like Then take the global minimum value. ,in Let be the set of elevation values ​​for the k-th layer.

[0094] 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;

[0095] Specifically, such as Figure 4 As shown, the specific steps of TIN filtering in step S3 are as follows:

[0096] 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. ;

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

[0098] ;

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

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

[0101] ;

[0102] 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 vertex The coordinate values.

[0103] Step S4: Use the ground points after TIN filtering as seed points to test the next layer. Perform secondary filtering on non-ground points based on terrain specificity. Update the points that pass the compensation judgment as seed points, and mark the rest as non-ground points.

[0104] Specifically, such as Figure 5 As shown, the specific steps of terrain-specific compensation in step S4 are as follows:

[0105] Step S41: Divide the non-ground points obtained by each layer of TIN filtering into unclassified points; perform secondary filtering using fitting and slope methods to select effective surface points for supplementation; the compensation process is divided into depression point and convex hump point compensation, slope point and fault point compensation according to the terrain characteristics.

[0106] Step S42: In mountainous and hilly areas, some unclassified points located in convex areas of the terrain may be misclassified as non-ground points because their angle with the TIN triangle exceeds the maximum angle threshold. Similarly, points located in depressions at the foot of mountains may also be lost due to angles exceeding the threshold. The least squares method is used to fit the local surface to estimate ground points. Unclassified points are then classified using a dynamic elevation difference threshold.

[0107] 1) To ensure the accuracy of surface fitting, expand the domain scope; retrieve at least 6 ground points within a 5×5 area centered on the unclassified point, and perform surface fitting;

[0108] 2) The equation of the two-dimensional polynomial surface is:

[0109] ;

[0110] Represented as a matrix equation:

[0111] ;

[0112] in:

[0113] , , ;

[0114] When m≥6 and When it is invertible, the optimal parameters need to be solved using the least squares method:

[0115] ;

[0116] 3) Calculate the fitted elevation difference d between each ground point and its corresponding fitted surface equation. i The mean d is obtained. mean The standard error σ is determined based on 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σ. Therefore, the elevation difference threshold for each point is determined as d. mean +3σ;

[0117] Step 43: Compensation algorithm for slopes and faults:

[0118] 1) Calculate the slope of each point relative to P by statistically analyzing the ground points in the eight neighboring areas around the unclassified point P;

[0119] 2) Set a slope threshold θᵣ = 10° as the standard for distinguishing between slopes and flat areas;

[0120] 3) Count the number n of ground points that satisfy θ < θᵣ. flat And the number n of ground points satisfying θ ≥ θᵣ slope The formulas for fault and slope compensation are as follows:

[0121] ;

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

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

[0124] 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.

[0125] Table 1

[0126] This invention integrates four key parameters: maximum mesh size W max Angle threshold Δthr based on TIN filtering, and maximum distance constraint S at each level. max and minimum distance threshold S min W max The size is usually determined based on the dimensions of the largest building, specifically its diagonal length, typically set at 30-40 meters. For mountainous areas with no buildings and significant terrain undulations, W... max The value should be appropriately reduced. The `thr` parameter is set based on the maximum slope of the terrain in the point cloud data; for mountainous areas, it is typically set to 30°-88°, while for flat areas, it is set to less than 30°. max This refers to the distance threshold between the point to be classified and the corresponding triangle in the initial triangulation network, typically set between 1.2 and 1.5 meters. min This refers to the distance threshold between the point to be classified and the corresponding triangle in the top-level triangular mesh, which should be less than S. max Generally, the minimum distance between a non-ground point and a ground point is used as a reference, and is usually set to 0.4-0.7 meters.

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

[0128] 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. This invention overcomes the inherent limitations of the traditional PTD algorithm 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 to compensate for elevation discontinuities in the terrain.

[0129] Table 2

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

[0131] Table 3

[0132]

[0133] The improved PTD algorithm proposed in this invention integrates three main aspects: multi-scale seed point selection, adaptive boundary expansion, and terrain-specific compensation. Validation based on the ISPRS standard dataset shows that this method significantly outperforms the classic PTD and other improved PTD algorithms in terms of total error (2.96%) and Kappa coefficient (89.36%). This invention optimizes the spatial distribution of seed points through the proportional attenuation coefficient of the layered 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 densification, reducing the Type I error (TI) to 2.06%, achieving a 70.5% performance improvement over the classic PTD. The use of extended grid interpolation and neighborhood variance constraints effectively solves the problem of boundary omission in fault and steep slope areas (such as S11 and S41), avoiding the misjudgment of narrow triangles caused by corner point insertion in traditional methods. Adaptive boundary expansion enhances the protection of edge information. Finally, the secondary filtering (surface fitting and slope filtering) recovers the missed ground points, significantly improving the problems of over-smoothing and feature loss in complex terrains caused by the classic PTD.

[0134] 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 multi-scale PTD filtering method based on terrain compensation, 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, and preprocess the data; 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 of mesh size; Step S4: Use the ground points after TIN filtering as seed points to test the next layer. Perform secondary filtering on non-ground points based on terrain specificity. Update the points that pass the compensation judgment 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 of terrain-specific compensation 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: Use the least squares method to fit the local surface to estimate ground points, and classify unclassified points using a dynamic elevation difference threshold. 1) Using the unclassified point as the center, retrieve at least 6 ground points within a 5×5 area and perform surface fitting; 2) The equation of the two-dimensional polynomial surface is: ； Represented as a matrix equation: ； in: , , ； When m≥6 and When it is invertible, the optimal parameters need to be solved using the least squares method: ； 3) Calculate the fitted elevation difference d between each ground point and its corresponding fitted surface equation. i The mean d is obtained. mean Based on the standard error σ, the height difference threshold d is determined for each point. mean +3σ; Step 43: Compensation algorithm for slopes and faults: 1) Calculate the slope of each point relative to P by statistically analyzing the ground points in the eight neighboring areas around the unclassified point P; 2) Set the slope threshold θ r = 10°, used as the standard to distinguish between slopes and flat areas; 3) Count the number of ground points n that satisfy θ < θᵣ flat And the number n of ground points satisfying θ ≥ θᵣ slope The formulas for fault and slope compensation are as follows: 。 2. The multi-scale PTD filtering method based on terrain compensation according to claim 1, characterized in that, The specific steps of data preprocessing in step S1 are as follows: Step S11: Retain the point cloud of the last echo and remove non-surface signals from the intermediate echoes; Step S12: 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(μ,σ), with μ+d std ×σ is the threshold value. In drone / vehicle-mounted LiDAR scenarios, the noise ratio is 0.5%~2%. d is set... std =2.0~2.5 Balanced noise removal and terrain preservation.

3. The multi-scale PTD filtering method based on terrain compensation according to claim 1, characterized in that, 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, the lowest layer L max Corresponding maximum grid size W max Only high-confidence ground points are selected as initial seeds, with a maximum grid size W. max The size W of the topmost mesh is determined by the largest non-ground object in the scene. min Related to the point density, the mesh size W of the intermediate layer k Determined by the proportional attenuation coefficient s, the increase factor of the grid seed points 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, and 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, Extension Point Distribution: Extend a new ring of seed points around the existing seed points, ensuring that all points are encompassed within the outer rectangle; the extended points are evenly distributed at equal intervals according to the grid size. Elevation calculations are performed based on the non-ground point interference included in the potential seed points, using the following method: 1) Target center coordinates and neighborhood filtering: Set the target center coordinates as... The corresponding point cloud set is P, and 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 analyze variance using a sliding window: ， ； In the formula, and For high program columns middle Elevation value Calculate the mean and variance of elevations; retain those that satisfy the following conditions. contiguous window set Finally, it was calculated that in a certain Interpolated elevation at point : ； like Then take the global minimum value of the elevation of that layer. ,in Let be the set of elevation values ​​for the k-th layer.

4. The multi-scale PTD filtering method based on terrain compensation according to 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 point to the vertex of the TIN face. If both the angle and the distance are less than a 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 vertex The coordinate values.

Images