A ground point cloud filtering method based on curvature dynamic segmentation and a storage medium

Abstract

The application discloses a ground point cloud filtering method based on curvature dynamic segmentation and a storage medium, belongs to the technical field of point cloud data processing, and carries out point cloud preprocessing and blocking; curvature dynamic segmentation and TIN cluster semantic analysis based denoising are carried out; for steep terrain areas, curvature threshold recursive segmentation sub-blocks are carried out; for flat areas, seed points with terrain representation are obtained by screening based on global elevation minimum point constraint and linear fitting. Then, progressive triangulation network iteration encryption optimization is carried out. Finally, a ground DEM based on grid filling is constructed, and a bilinear interpolation filtering based on the ground DEM is carried out, the difference between the actual elevation of the point cloud and the DEM elevation is calculated, if the difference is less than the final distance threshold, the ground point is marked, and the ground classification of the whole point cloud is completed. The application specifically solves the core problems of large filtering data volume, low processing efficiency and poor adaptability to complex terrain of the traditional triangulation network, and has high filtering precision and engineering practicability.

Description

Technical Field

[0001] This invention belongs to the technical field of point cloud data processing, specifically relating to a ground point cloud filtering method and storage medium based on curvature dynamic segmentation. Background Technology

[0002] LiDAR technology, with its advantages of high precision and high efficiency in acquiring three-dimensional spatial information, has become an indispensable core technology support for fields such as topographic mapping, smart city construction, digital twin modeling, and geological disaster monitoring. Among them, ground filtering, as a key core link in LiDAR point cloud data preprocessing, has the core objective of accurately separating clean ground points from point cloud data that mixes ground points (such as soil, roads, and flat land) and non-ground points (such as buildings, trees, and obstacles). This provides high-quality basic data for subsequent tasks such as digital elevation model (DEM) generation, terrain analysis, and 3D scene reconstruction. The processing effect directly determines the accuracy and reliability of downstream applications.

[0003] Currently, as LiDAR technology expands into large-scale, high-density, and complex scenarios, higher demands are placed on the accuracy, robustness, terrain adaptability, and processing efficiency of ground filtering. However, existing ground filtering methods can no longer meet the high-precision filtering requirements of complex terrains (such as urban building complexes and mountainous ravines) and large-scale point cloud data. Specifically, in complex scenarios, existing technologies have the following shortcomings:

[0004] (1) Weak adaptability to complex terrain; In scenarios such as densely built-up urban areas, mountainous and hilly areas, and vegetation-covered areas, non-ground points and ground points often have spatial adhesion and high overlap, and the curvature of the terrain changes drastically. The fixed parameters (such as slope threshold, seed point spacing, and structural element size) relied upon by traditional methods are difficult to adapt to diverse terrain features, and problems such as "filtering wave" (falsely deleting effective surface points) or "missed filtering" (not removing non-ground points) are likely to occur.

[0005] (2) The selection of seed points lacks a scientific optimization mechanism; the seed points of existing mainstream methods such as progressive triangulation (TIN) mostly adopt fixed spacing or random sampling strategies, which cannot dynamically adjust the selection logic according to the curvature of the local terrain, resulting in limited ground fitting accuracy.

[0006] (3) Insufficient robustness; When faced with noisy points, low-density areas or dense distribution of non-ground points in LiDAR point cloud data, the existing methods have weak anti-interference ability and the filtering accuracy is significantly reduced.

[0007] (4) Efficiency and accuracy are difficult to balance; high-precision methods (such as deep learning semantic segmentation and fine geometric constraint methods) have high computational complexity and long inference time, making them difficult to adapt to the needs of large-scale real-world data processing, while efficient methods (such as simple morphology and fixed threshold methods) cannot guarantee segmentation accuracy in complex scenarios, resulting in a core contradiction between accuracy and efficiency. Summary of the Invention

[0008] The purpose of this invention is to provide a ground point cloud filtering method and storage medium based on curvature dynamic segmentation, aiming to solve any of the above-mentioned problems.

[0009] This invention is mainly achieved through the following technical solutions:

[0010] A ground point cloud filtering method based on curvature dynamic segmentation includes the following steps:

[0011] Step S1: Point cloud preprocessing and segmentation to obtain point cloud blocks;

[0012] Step S2: Initial screening of potential ground points and dynamic segmentation based on curvature: For steep terrain areas, recursively segment into sub-blocks based on curvature thresholds, and extract the lowest elevation point in each sub-block, summarizing them into an initial seed point set; for flat areas, based on global lowest elevation point constraints and linear fitting, select seed points that are representative of the terrain.

[0013] Then, based on TIN cluster semantic analysis and denoising, the optimized seed point set is obtained;

[0014] Step S3: Iteratively encrypt and optimize the triangulation, and output the encrypted triangulation and the optimized seed point set;

[0015] Step S4: Construct a ground DEM based on grid filling, and perform bilinear interpolation filtering based on the ground DEM to calculate the difference between the actual elevation of the point cloud and the elevation of the DEM. If the difference is less than the final distance threshold, it is marked as a ground point, thus completing the full point cloud ground classification.

[0016] To better implement the present invention, step S1 further includes the following steps:

[0017] Step S11: Adaptive downsampling based on point cloud density; if the density of the original point cloud is greater than the density threshold, random downsampling is performed; otherwise, the original point cloud is directly retained.

[0018] Step S12: Divide the point cloud into L×L blocks evenly to generate several point cloud blocks.

[0019] To better realize the present invention, step S2 further includes the following steps:

[0020] Step S21: Preliminary screening of potential ground points; First, the point cloud blocks are divided into grids, the lowest elevation point in each grid and the global lowest elevation point of the point cloud block are extracted, and the angle α between the line connecting the lowest elevation point and the global lowest elevation point and the horizontal plane is determined. If α is less than the set angle threshold, it is selected as a potential ground point.

[0021] Step S22: Project the potential ground points in the X / Y directions and calculate their curvature. Recursively divide the potential ground points into sub-blocks according to the curvature threshold, extract the lowest elevation point in the final sub-block, and form an initial seed point set.

[0022] Step S23: Denoise the initial seed point set based on TIN cluster semantic analysis; construct an initial triangular network (TIN) based on the initial seed point set, traverse all triangular units in the initial TIN, filter out suspected triangular units by the slope or maximum elevation difference of the triangles, and cluster them to obtain suspected triangular clusters; then, combine semantic features to remove high-position noise points and building-type pseudo-ground points.

[0023] Step S24: Based on the denoised seed point set, perform two-dimensional projection and linear fitting to remove residual noise points that deviate from the ground trend, and obtain the optimized seed point set.

[0024] To better realize the present invention, step S22 further includes the following steps:

[0025] Step A1: Using the center coordinates (X_center, Y_center) of the point cloud block of the potential ground points as the reference, extract the point cloud subsets [X_center-L / 2, X_center+L / 2] and [Y_center-L / 2, Y_center+L / 2] in the X-axis and Y-axis directions respectively;

[0026] Step A2: Project the point cloud subset in the X direction onto the XZ plane, project the point cloud subset in the Y direction onto the YZ plane, and sort them in ascending order of X / Y coordinates to obtain a continuous sequence of Z values ​​as X / Y change.

[0027] Step A3: Calculate the terrain curvature κ in the X-axis and Y-axis directions for each point using the numerical differentiation method. x and κᵧ;

[0028] If κ x If the curvature threshold κ0 is greater than the curvature threshold κ0, then the X coordinate of the point is marked as the X-direction dividing line; if κᵧ > curvature threshold κ0, then the Y coordinate of the point is marked as the Y-direction dividing line; the current point cloud is divided into several sub-blocks according to the dividing lines;

[0029] Step A4: Repeat steps A2 and A3 to dynamically recursively divide each sub-block until a stopping condition is met. The stopping condition is: κ of all points within the sub-block.x ≤κ0 and κᵧ≤κ0; or the size of the sub-block in the X / Y direction is less than or equal to the preset minimum size S0;

[0030] Step A5: Extract the lowest elevation point within each sub-block and summarize them into an initial seed point set.

[0031] To better implement the present invention, further, in step S23, if all triangles of the suspected triangle cluster share a common vertex and the triangles are connected end to end, then the suspected triangle cluster is determined to be a high-noise cluster, its common vertex is a high-noise point, and the high-noise point is removed from the initial seed point set.

[0032] If the architectural semantic features of a suspected triangle cluster satisfy the following conditions: ① the number of triangles in the cluster is greater than or equal to the minimum number of triangles N1 in the architectural triangle cluster; ② the average angle between the normal vector of all triangles in the cluster and the vertical direction is greater than 60°; ③ the aspect ratio of the minimum bounding box of the cluster is greater than the minimum aspect ratio λ of the minimum bounding box of the architectural triangle cluster and the area of ​​the bounding box is greater than or equal to the minimum area S1 of the minimum bounding box of the architectural triangle cluster; ④ the ratio of the total area of ​​triangles in the cluster to the area of ​​the cluster's convex hull is less than the threshold μ; then the suspected triangle cluster is determined to be an architectural pseudo-ground cluster, and all seed points in the cluster are removed from the initial seed point set.

[0033] To better implement the present invention, step S24 further includes the following steps:

[0034] Step B1: Project the semantically denoised seed points onto the XY plane, divide them into segments with an interval width I along the X-axis, extract the lowest elevation point in each segment, and form a core feature point set.

[0035] Step B2: Calculate the slope K between two adjacent points in the core feature point set, mark the locations where the absolute value of the slope is greater than the slope threshold K1 and the difference between adjacent slopes is greater than the slope difference threshold K2 as terrain change points, and remove the corresponding high-position noise points;

[0036] Step B3: Fit a straight line to adjacent point pairs in the core feature point set, calculate the vertical distance d from all points between the two points to the straight line, mark the points whose d > distance threshold D3 and remove them;

[0037] Step B4: Divide the denoised 2D point cloud into blocks twice according to the preset step size S. For sub-segments with a number of points greater than or equal to the minimum number of points N4, fit a straight line using the Tukey robust least squares method. If the absolute value of the slope of the fitted line is greater than the terrain steepness threshold K3, mark outlier noise points according to the 3σ criterion. If the absolute value of the slope is less than the terrain steepness threshold K3, retain all points in the sub-segment. Map the points marked as retained back to the 3D point cloud through indexing to obtain the optimized seed point set.

[0038] To better realize the present invention, step S3 further includes the following steps:

[0039] Step S31: Traverse all point cloud blocks and mark point cloud blocks with a seed point count greater than 0 as valid blocks;

[0040] Step S32: Based on the effective blocks, extract the X / Y coordinates of the seed points and construct a KD-Tree. Search the points in the neighborhood of each point, retain the point with the lowest elevation in the neighborhood as the seed point, merge the seed points of the current effective block and the eight neighboring effective blocks to generate an extended seed point set, and construct the Delaunay triangulation of the extended seed point set.

[0041] Step S33: Calculate the XYZ bounding box of the potential ground point cloud within the current valid block, divide the X / Y plane grid, traverse all triangles of the triangulation, map the triangle index to the corresponding grid and remove duplicates, generate a global mapping table of grid-triangle, and match potential ground points with triangulation cells.

[0042] Step S34: Determine the ground attributes of potential ground points based on the slope of the triangular network; traverse each potential ground point that matches a triangular network unit, calculate the normal vector and slope of the triangular plane, and mark the potential ground points corresponding to the triangular planes that meet the threshold range as ground points;

[0043] Step S35: Identify pseudo-ground triangles, remove the vertex with the highest elevation, and add the potential ground point farthest from the triangle below it as a real ground seed point, and update the seed point set;

[0044] Step S36: Repeat steps S32 to S35 until the iteration termination condition is met, and output the final encrypted triangular network and the optimized seed point set.

[0045] To better implement the present invention, further, in step S34, if the approximate planar slope threshold θ3 < the slope of the triangular plane ≤ the steep terrain slope threshold θ2, and if the directed distance d1 from the potential ground point to the triangular plane < the dynamic distance threshold d... 10 If the angle between the potential ground point and the triangular plane is less than the dynamic angle threshold θ4, then the potential ground point is marked as a ground point.

[0046] If the slope of the triangular plane is ≤ θ3, then the potential ground points below the triangular plane are directly marked as ground points, and the points located above the triangular plane and inside the triangle with the required angle are marked as ground points.

[0047] If the slope of the triangular plane is greater than the steep terrain slope threshold θ2, then the highest vertex of the triangle is identified, a mirror point of the potential ground point is generated, and the distance d2 from the mirror point to the triangular plane and the included angle α2 are calculated. If d2 < dynamic distance threshold d 10 If α2 < dynamic angle threshold θ4, then the original potential ground point is marked as a ground point.

[0048] To better realize the present invention, step S4 further includes the following steps:

[0049] Step S41: Calculate the full point cloud bounding box, initialize the DEM grid according to the preset resolution, traverse the ground seed points, assign each seed point to the corresponding grid according to the X / Y coordinates, extract the lowest elevation point in each grid as the initial DEM grid point; mark the grid with valid points as filled, and mark the grid without data as to be filled.

[0050] Step S42: Collect the coordinates of the blank grid, search for the eight-direction valid grid of the blank grid in stages, fill it with linear interpolation and verify the validity of the Z value; summarize the filled grid points, filter invalid Z values ​​and generate an irregular DEM grid.

[0051] Step S43: Convert the irregular DEM grid into a regular DEM grid; traverse the valid nodes of the irregular DEM grid, and calculate the Z value of the regular DEM grid node based on the concavity and convexity of the quadrilateral formed by the four neighbors of the DEM grid; for convex quadrilaterals, calculate the Z value of the nearest neighbor point using bilinear interpolation; for non-convex quadrilaterals, calculate the Z value of the nearest neighbor point, and finally generate a complete regular DEM grid.

[0052] Step S44: Perform bilinear interpolation based on the regular DEM grid to calculate the difference between the actual elevation of the point cloud and the DEM elevation. If the difference is less than the final distance threshold, it is marked as a ground point, and the ground classification is completed.

[0053] A computer-readable storage medium having a computer program stored thereon that, when executed by a processor, implements the above-described ground point cloud filtering method based on curvature dynamic segmentation.

[0054] The beneficial effects of this invention are as follows:

[0055] (1) This invention specifically addresses the core problems of traditional triangulation filtering, such as large data volume, low processing efficiency, and poor adaptability to complex terrain. In the point cloud preprocessing stage, this invention significantly reduces the data volume through downsampling and block parallel strategies; in the seed point extraction stage, it ensures the quality of seed points by relying on curvature dynamic segmentation and a two-level denoising mechanism; in the triangulation iteration stage, it improves computational efficiency through rasterization for rapid matching; and in the DEM construction and interpolation stage, it replaces the traditional full iteration with gridded processing, achieving the dual goals of "efficiency improvement + accuracy assurance".

[0056] (2) This invention significantly improves the ground extraction accuracy of LiDAR point clouds in complex terrains through a collaborative innovative design that optimizes seed points by dynamic curvature segmentation and semantic denoising based on TIN clusters and bilinear interpolation verification. It maintains good robustness to scenes with noisy points, low-density areas, and densely distributed non-ground points. Furthermore, through block-based parallel computation, it achieves high-efficiency point cloud processing, possessing high filtering accuracy and engineering practicality. This invention provides reliable technical support for the in-depth application of LiDAR data; it achieves accurate and efficient separation of ground points and non-ground points in large-scale LiDAR point cloud data, providing high-quality basic data support for downstream applications such as digital elevation model generation, terrain analysis, and 3D scene reconstruction.

[0057] (3) This invention first adapts to terrain abrupt changes through curvature recursive segmentation, and then combines TIN cluster semantic analysis for denoising, two-dimensional projection and linear fitting for refined denoising, successively eliminating high-position noise points, building-type pseudo-ground points and residual noise points, solving the problems of low quality and large noise interference of traditional seed points, and obtaining a high-quality initial seed point set, laying the foundation for rapid construction of triangulation networks. Specifically, in steep terrain areas, this invention achieves multiple densification effects of seed points through dynamic recursive segmentation, ensuring the accuracy of terrain representation of steep terrain; in flat areas, this invention retains seed points with terrain representativeness through global minimum elevation point constraints and linear fitting screening, avoiding redundant points from interfering with the efficiency of subsequent triangulation network construction.

[0058] (4) Based on the semantic analysis of TIN clusters and the slope mutation-robust fitting denoising mechanism, this invention effectively removes building-type pseudo-ground noise, low vegetation noise and isolated high-position noise points in the original point cloud. The final optimized seed point set not only accurately matches the real terrain undulation features, but also minimizes the interference of non-ground noise, laying a high-quality data foundation for subsequent triangular network iteration encryption and DEM construction.

[0059] (5) This invention generates high-precision regular DEM quickly through a phased grid filling and concavity / convexity fitting strategy. The bilinear interpolation classification based on DEM avoids the triangular mesh matching of the entire point cloud one by one, greatly improving processing efficiency and outputting high-precision ground point cloud. Attached Figure Description

[0060] Figure 1 This is a flowchart of the ground point cloud filtering method based on curvature dynamic segmentation according to the present invention;

[0061] Figure 2 This is a schematic diagram of obtaining seed points through dynamic region segmentation of the point cloud in Example 2;

[0062] Figure 3 This is a schematic diagram illustrating the extraction of suspected triangular clusters from the triangular mesh in Example 2;

[0063] Figure 4 This is a schematic diagram of high-bit noise extraction in Example 2;

[0064] Figure 5 This is a schematic diagram of building noise extraction in Example 2;

[0065] Figure 6 This is a schematic diagram of the original point cloud data in Example 2;

[0066] Figure 7 This is a diagram showing the effect of seed points extracted after dynamic segmentation and denoising in Example 2;

[0067] Figure 8 This is a schematic diagram of point cloud extraction after ground filtering in Example 2. Detailed Implementation

[0068] Example 1:

[0069] A ground point cloud filtering method based on curvature dynamic segmentation, such as Figure 1 As shown, it includes the following steps:

[0070] Step 1: Preprocess the point cloud and divide it into blocks to obtain point cloud blocks.

[0071] Step 1.1: Determine the relationship between the original LiDAR point cloud density and the preset density threshold ρ0. If the point cloud density exceeds the threshold, perform random downsampling; otherwise, retain the point cloud directly.

[0072] Step 1.2: Divide the processed point cloud into uniform blocks of a fixed size of L×L to obtain point cloud blocks.

[0073] The core optimization objective of Step 1 is to significantly reduce the amount of data and accelerate subsequent triangulation processing; reduce redundant data through density adaptive downsampling; and decompose large-scale point clouds into sub-blocks suitable for processing by combining uniform block division and merging mechanisms to avoid efficiency bottlenecks caused by overall processing, while balancing processing speed and terrain feature integrity.

[0074] Step 2: Initially screen potential ground points and perform dynamic segmentation based on curvature. Then, denoise the points based on TIN cluster semantic analysis to obtain an optimized seed point set.

[0075] Step 2.1: Divide the point cloud into blocks according to the grid, extract the lowest elevation point of the grid and the lowest global elevation point of the block, and filter potential ground points by the angle between the line connecting the two points and the horizontal line.

[0076] Step 2.2: Project the potential ground point cloud in the X / Y direction and calculate its curvature. Recursively divide the potential ground point cloud into sub-blocks according to the curvature threshold and extract the lowest elevation point in each sub-block to form an initial seed point set.

[0077] Step 2.3: Construct an initial triangulation network (TIN) based on the initial seed points, filter suspected triangle clusters by the slope and elevation difference of the triangles, and remove high-position noise points and building-type pseudo-ground points by combining semantic features;

[0078] Step 2.4: Perform two-dimensional projection and linear fitting on the denoised seed points to remove residual noise points that deviate from the ground trend, and obtain the optimized seed point set.

[0079] The optimization point of step 2 is to construct a seed point purification system of "curvature dynamic segmentation + two-level denoising". First, the curvature recursive segmentation is used to adapt to terrain changes. Then, TIN cluster semantic analysis is used for denoising, and two-dimensional projection and linear fitting are used for refined denoising. High-position noise points, building pseudo-ground points and residual noise points are removed in turn to solve the problems of low quality and large noise interference of traditional seed points. A high-quality initial seed point set is obtained, which lays the foundation for rapid construction of triangulation.

[0080] Step 3: Iterative encryption and optimization of the progressive triangulation network;

[0081] Step 3.1: Mark the point cloud blocks with a seed point count greater than 0 as valid blocks, and initialize the iteration parameters and termination conditions;

[0082] Step 3.2: Search for points in the neighborhood of each point, retain the point with the lowest elevation in the neighborhood as the seed point, merge the seed points of the effective block and the eight-neighbor block to generate an extended seed point set, and construct the Delaunay triangulation of the extended seed point set.

[0083] Step 3.3: Calculate the XYZ bounding box of the potential ground point cloud within the current valid block, divide the X / Y plane grid, and determine the grid range covered by the point cloud; traverse all triangles of the triangulation, calculate their XY plane bounding boxes, determine the grid range covered, map the triangle index to the corresponding grid and remove duplicates, generate a grid-triangle global mapping table, and quickly match potential ground points with triangles.

[0084] Step 3.4: Determine the ground attributes of potential ground points based on the planar characteristics of the triangle; for steep terrain, use mirror points for assistance in determination.

[0085] Step 3.5: Identify pseudo-ground triangles, remove unreasonable vertices and supplement with real ground seed points, and update the seed point set;

[0086] Step 3.6: Calculate the rate of change of the number of seed points. If the iteration termination condition is met, stop and output the final encrypted triangular network and the optimized seed point set; otherwise, proceed to step 3.2.

[0087] Step 4: Construct a grid-filled ground DEM and perform bilinear interpolation filtering to complete the full point cloud ground classification.

[0088] Step 4.1: Calculate the full point cloud bounding box, initialize the DEM grid at the preset resolution, assign seed points and mark the filled grid;

[0089] Step 4.2: Search for the eight-directional valid grid of the blank grid in stages, fill it with linear interpolation and verify the validity of the Z value;

[0090] Step 4.3: Determine the concavity or convexity of the quadrilateral formed by the four neighboring regions of the DEM grid, and generate a regular DEM grid by bilinear interpolation or the nearest neighbor method;

[0091] Step 4.4: Perform bilinear interpolation based on the regular DEM grid to calculate the difference between the actual elevation of the point cloud and the DEM elevation, and complete the ground classification according to the threshold.

[0092] The core optimization objective of step 4 is to replace the traditional full-scale iteration of the triangulation network, further accelerate the classification process, and quickly generate a high-precision regular DEM through a phased grid filling and concavity / convexity fitting strategy. The bilinear interpolation classification based on the DEM avoids the triangulation network from matching the entire point cloud one by one, greatly improving processing efficiency and outputting a high-precision ground point cloud.

[0093] Example 2:

[0094] A ground point cloud filtering method based on curvature dynamic segmentation.

[0095] like Figure 1 As shown, it is mainly divided into four parts: point cloud preprocessing, seed point optimization and extraction, triangular network iterative encryption, and DEM construction and interpolation filtering. Each step is closely connected to achieve accurate extraction of ground points.

[0096] Step 1, point cloud preprocessing;

[0097] Step 1.1, Determining and Downsampling the Original LiDAR Point Cloud Density; The original LiDAR point cloud data is collected by airborne or ground-based lidar equipment, covering mixed features such as terrain, buildings, and vegetation. The point cloud density is calculated (point cloud density = total number of points / actual area of ​​the region), with a preset density threshold ρ0. If the original point cloud density > ρ0, a random downsampling algorithm is used to retain the target density ρ, avoiding data redundancy that could reduce processing efficiency; if the density ≤ ρ0, the original point cloud is directly retained without performing downsampling.

[0098] Step 1.2, Uniformly divide and merge point cloud: Divide the downsampled (or undownsampled) point cloud into uniform blocks of a fixed size L×L to generate several independent point cloud blocks.

[0099] Step 2, Seed point optimization extraction based on curvature dynamic segmentation;

[0100] Step 2.1, preliminary screening of potential ground points; set the grid resolution r, and divide each valid block into grids. Extract the lowest elevation point Pij within each grid, and simultaneously extract the global lowest elevation point Pi of that block; calculate the angle α between the line connecting Pij and Pi and the horizontal plane, and set an angle threshold θ. If α < θ, it indicates that Pij has good terrain continuity with the global lowest elevation point and is included in the potential ground point set; if α ≥ θ, it is determined to be a high-altitude non-ground point and is not included for the time being.

[0101] Step 2.2, dynamic region segmentation based on curvature;

[0102] (1) such as Figure 2 As shown, using the center coordinates (X_center, Y_center) of each potential ground point block as a reference, subsets of the point cloud are extracted along the X-axis ([X_center-L / 2, X_center+L / 2]) and along the Y-axis ([Y_center-L / 2, Y_center+L / 2]). The X-axis subset is projected onto the XZ plane, and the Y-axis subset is projected onto the YZ plane. These subsets are then sorted in ascending order of X / Y coordinates to obtain a continuous sequence of Z values ​​as a function of X / Y.

[0103] (2) The topographic curvature κ (topographic curvature κ in the X direction) of each point is calculated using the numerical differential method. x The topographic curvature κᵧ in the Y direction is calculated using the following formula: κᵧ x =|d²Z / dX²| / [(1+(dZ / dX)²)^(3 / 2)] (The principle of terrain curvature calculation is the same in the X and Y directions), and the curvature threshold κ0 is set. If a certain point κ x If κᵧ > κ0, mark the X coordinate of the point as the X-direction dividing line; if κᵧ > κ0, mark the Y coordinate of the point as the Y-direction dividing line.

[0104] (3) Divide the current block into multiple sub-blocks according to the dividing line;

[0105] (4) Repeat the above (1)~(3) process of "projection-calculation of curvature-marking of dividing lines-segmentation" for each sub-block to achieve dynamic recursive segmentation. Continue until any of the following stopping conditions are met: ① κ of all points in the sub-block x ≤κ0 and κᵧ≤κ0 (the terrain within the sub-block is continuous without abrupt changes); ② The X / Y dimensions of the sub-block are both ≤ the preset minimum dimension S0.

[0106] (5) Finally, extract the lowest elevation point in each final sub-block and summarize it into an initial seed point set.

[0107] Step 2.3, Denoising of seed points based on TIN cluster semantic analysis;

[0108] (1) Based on the initial seed point set, the Delaunay triangulation algorithm is used to construct the initial triangulation network (TIN). The empty circle property and the minimum angle maximization property of the triangulation network TIN are used to ensure the reliability of the terrain geometric representation. Set the slope threshold θ1 and the elevation difference threshold h1, and traverse all triangular elements in the triangulation network TIN: if the slope of the triangle is greater than θ1, or the maximum elevation difference between the three vertices of the triangle is greater than h1, it is marked as a "suspected triangulation element".

[0109] (2) such as Figure 3 As shown, based on the adjacency relationship of triangles, suspected triangular units are clustered to form several "suspect triangular clusters". Semantic determination and denoising are then performed on each suspected cluster.

[0110] like Figure 4 As shown, if all triangles in a suspected cluster share a common vertex and the triangles are connected end to end (adjacent triangles share a side), it is determined to be a "high-order noise cluster", and its common vertex is a high-order noise point, which is removed from the initial seed point set;

[0111] like Figure 5 As shown, if a suspected cluster satisfies the following architectural semantic features: ① The number of triangles in the cluster is greater than or equal to the minimum number of triangles N1 in the architectural triangle cluster; ② The average angle between the normal vector of all triangles in the cluster and the vertical direction is greater than 60°; ③ The aspect ratio of the minimum bounding box of the cluster is greater than the minimum aspect ratio λ of the minimum bounding box of the architectural triangle cluster and the area of ​​the bounding box is greater than or equal to the minimum area S1 of the minimum bounding box of the architectural triangle cluster; ④ The ratio of the total area of ​​triangles in the cluster to the area of ​​the cluster's convex hull is less than the threshold μ, then it is determined to be an "architectural pseudo-ground cluster", and all seed points in the cluster are removed from the initial seed point set.

[0112] Step 2.4, refine and denoise the seed points based on two-dimensional projection-slope change-linear fitting;

[0113] If the number of points in the current seed point sub-block is less than the minimum number of point clouds N2 in the sub-block, return directly; otherwise, project the semantically denoised seed points onto the XY plane, divide them into segments with an interval width I along the X-axis, extract the lowest elevation point in each segment, and form a core feature point set.

[0114] Calculate the slope K = ΔZ / ΔX between adjacent points in the core feature point set. Set slope threshold K1 and slope difference threshold K2. Mark locations where the absolute slope value > K1 and the adjacent slope difference > K2 as terrain abrupt change points and remove the corresponding high-position noise points. Fit a straight line to adjacent point pairs in the core feature point set, calculate the perpendicular distance d from all points between the two points to the line, set a distance threshold D3, and mark points where d > D3 as "removed".

[0115] The denoised 2D point cloud is divided into blocks a second time according to a preset step size S. For blocks with a number of points greater than or equal to the minimum number of points N4, Tukey's robust least squares method is used to fit a straight line. If the absolute value of the slope of the fitted line is greater than the terrain steepness threshold K3, outlier noise points are marked according to the 3σ criterion (distance ≥ μ + 3σ). If the absolute value of the slope is less than the terrain steepness threshold K3, all points in the block are retained. The points marked as "retained" are mapped back to the 3D point cloud through indexing to obtain the optimized seed point set.

[0116] like Figure 6 and Figure 7 As shown, based on the above processing, this invention achieves multi-densification of seed points in steep terrain areas (such as slopes, valleys, and other places with abrupt curvature changes) through dynamic recursive segmentation. Each tiny continuous terrain unit can extract a core seed point, ensuring the accuracy of terrain representation for steep terrain. In flat areas (such as plains, gentle slopes, and other places with gentle curvature), the invention retains representative seed points through global elevation minimum point constraints and linear fitting, avoiding redundant points from interfering with the efficiency of subsequent triangulation network construction. Simultaneously, relying on TIN cluster semantic analysis and a slope abrupt change-robust fitting denoising mechanism, this invention effectively removes building-related pseudo-ground noise, low-lying vegetation noise, and isolated high-position noise points from the original point cloud. The resulting optimized seed point set accurately matches the real terrain undulation features while minimizing the interference of non-ground noise, laying a high-quality data foundation for subsequent triangulation network iteration and DEM construction.

[0117] Step 3: Iterative encryption and optimization of the progressive triangular network;

[0118] Step 3.1: Validity marking of point cloud blocks, and initialization of iteration parameters and termination conditions;

[0119] Traverse all point cloud blocks. If the number of optimized seed points in a block is greater than 0, mark it as a "valid block"; if there are no seed points in a block, mark it as an "invalid block".

[0120] Initialize the number of iterations k=1, and set the iteration termination conditions: ① the rate of change of the number of seed points of all valid blocks ≤ Δk; ② the number of iterations k≥k_max; the iteration terminates when either condition is met.

[0121] Step 3.2: Segmented seed point grid selection and neighborhood TIN construction;

[0122] For each valid block, a minimum distance threshold for the grid is set, the X / Y coordinates of the seed point are extracted to construct a KD-Tree, and the points in the neighborhood of each point are searched. The point with the lowest elevation in the neighborhood is retained as the seed point, and the remaining points are assigned to the potential ground point set.

[0123] Calculate the row and column indices of the current block, define the offsets of the eight neighborhoods (top, bottom, left, right, and four diagonals), merge the seed points of the current valid block with the seed points of the neighboring valid blocks, generate an extended seed point set, and construct a Delaunay triangulation based on this set.

[0124] Step 3.3: Fast matching of potential ground points with triangular rasterization;

[0125] Calculate the XYZ bounding box of the potential ground point cloud within the current block, set the grid size s, divide the X / Y plane grid, and determine the grid range covered by the point cloud. Traverse all triangles of the updated triangulation, calculate their XY plane bounding boxes, determine the grid range covered, map the triangle indices to the corresponding grids and remove duplicates, and generate a grid-triangle global mapping table.

[0126] Step 3.4: Determine the ground attributes of potential ground points based on TIN plane features;

[0127] Iterate through each potential ground point that matches the triangle, and calculate the triangle's plane normal and slope:

[0128] If the slope of the triangle is less than or equal to the steep terrain slope threshold θ2 (the value range is less than 90°, and 88° is taken in this embodiment): ① If the slope is greater than the approximate plane slope threshold θ3 (which can be considered as the minimum slope of the plane, and 2° is taken in this embodiment), calculate the directed distance d1 from the potential ground point to the triangular plane, and set the dynamic distance threshold d. 10 (Value range 0.3-1m, 0.6m in this embodiment), if d1 < dynamic distance threshold d 10 Furthermore, if the angle between a potential ground point and the triangular plane is less than θ4 (ranging from 5° to 15°, and 10° in this embodiment), it is marked as a ground point; ② If the slope is less than or equal to θ3, potential ground points below the triangular plane are directly marked as ground points, and points above the plane that are located inside the triangle and whose angles meet the requirements are also marked as ground points;

[0129] If the triangle slope > θ2: Identify the highest vertex of the triangle, generate a mirror image of the potential ground point (mirror image X / Y coordinates = 2 × highest vertex X / Y coordinates - potential ground point X / Y coordinates, Z coordinate is the same as the potential ground point), calculate the distance d2 from the mirror image point to the triangle plane and the included angle α2. If d2 < dynamic distance threshold d 10 And α2 < dynamic angle threshold θ4, mark the original potential ground point as a ground point.

[0130] Step 3.5: Pseudo-ground triangle identification and seed point optimization and update;

[0131] Determine whether the number of potential ground points n1 below each triangular unit and the average directed distance d_avg are both greater than the threshold. If so, it is determined to be a pseudo-ground triangle.

[0132] For a pseudo-ground triangle, remove its unreasonable vertices (the vertices with the highest elevation) and add the potential ground point (already identified as a ground point) that is farthest away from the triangle as a new seed point.

[0133] Step 3.6, Iteration termination determination and block state update;

[0134] Add all potential ground points marked as ground points to the seed point set, retain the remaining potential ground points, and update the potential ground point set. If the rate of change of the number of seed points in the current block Δ < Δk (the ratio of the number of newly added seed points in this iteration to the total number of seed points after this iteration), then mark the block as invalid and no longer participate in subsequent iterations. If the termination condition is not met, return to step 3.2 to continue iterating; if the termination condition is met, terminate the iteration and output the final ground seed points.

[0135] Step 4: Construct a grid-filled ground DEM and perform bilinear interpolation filtering to complete the full point cloud ground classification;

[0136] Step 4.1, DEM grid initialization and seed point grid allocation;

[0137] Calculate the XYZ bounding box of the entire point cloud and expand the grid boundary outward by a preset expansion size r_dem to avoid classification omissions caused by insufficient grid range in the edge point cloud. Traverse the final ground seed points and assign each seed point to the corresponding grid according to its X / Y coordinates. Extract the lowest elevation point in each grid as the initial DEM grid point, and mark the grid with valid points as "filled" and the grid without data as "to be filled".

[0138] Step 4.2: Collect the coordinates of blank grid cells. In the first stage, the maximum search distance is set to a preset proximity, and valid grid cells are searched in eight directions. For blank grid cells with paired valid grid cells, the Z-value is calculated using linear interpolation (verified to be within a reasonable range of the original point cloud elevation) and then filled. For blank grid cells not filled in the first stage, the search distance is increased in the second stage, and the above logic is repeated. Finally, the filled grid points are summarized, invalid Z-values ​​are filtered out, and an irregular DEM grid is generated.

[0139] Step 4.3: Convert irregular grid to regular grid and fit Z-value;

[0140] Traverse the valid nodes of the irregular DEM grid and determine the concavity / convexity of the quadrilateral using the four neighboring grid points as units; for convex quadrilaterals, use bilinear interpolation to calculate the Z value of the regular grid nodes; for non-convex quadrilaterals or abnormal interpolation parameters, take the Z value of the nearest neighbor point, and finally generate a complete regular DEM grid.

[0141] Step 4.4, full point cloud ground classification based on bilinear interpolation;

[0142] Initialize the point cloud ground marker array (default is non-ground), determine the regular grid index to which each point belongs and verify the validity of the neighborhood Z value; calculate the DEM elevation value corresponding to the point through bilinear interpolation, and solve for the difference ΔZ between the actual point cloud elevation and the DEM elevation. If ΔZ ≤ the final distance threshold ΔZ0, then mark it as a "ground point"; otherwise, mark it as a "non-ground point".

[0143] This invention specifically addresses the core problems of traditional triangulation filtering, such as large data volume, low processing efficiency, and poor adaptability to complex terrain. In the point cloud preprocessing stage, this invention significantly reduces the data volume through downsampling and block-based parallel strategies; in the seed point extraction stage, it ensures seed point quality by relying on dynamic curvature segmentation and a two-stage denoising mechanism; in the triangulation iteration stage, it improves computational efficiency through rasterization for rapid matching; and in the DEM construction and interpolation stage, it replaces traditional full-scale iteration with gridded processing, achieving the dual goals of "efficiency improvement + accuracy assurance."

[0144] like Figure 6 and Figure 8 As shown, based on point cloud data of complex scenes containing steep terrain, flat areas, building clusters, and low vegetation, this invention significantly optimizes the classification effect of point clouds. In steep terrain areas, due to the dense layout of seed points, the ground points are extracted completely without obvious breaks, accurately restoring the topographic undulation features of slopes and valleys. In flat areas, representative ground points are retained without redundant data interference. Building noise and low vegetation noise have been removed, and there are no false ground point residues. This fully verifies the practicality and superiority of this invention in complex scenes, and can provide high-quality basic data support for downstream applications such as digital elevation model generation, terrain analysis, and 3D scene reconstruction.

[0145] The above description is merely a preferred embodiment of the present invention and is not intended to limit the present invention in any way. Any simple modifications or equivalent changes made to the above embodiments based on the technical essence of the present invention shall fall within the protection scope of the present invention.

Claims

1. A ground point cloud filtering method based on curvature dynamic segmentation, characterized in that, Includes the following steps: Step S1: Point cloud preprocessing and segmentation to obtain point cloud blocks; Step S2: Initial screening of potential ground points and dynamic segmentation based on curvature: For steep terrain areas, recursively segment into sub-blocks based on curvature thresholds, and extract the lowest elevation point in each sub-block, summarizing them into an initial seed point set; For flat areas, seed points with terrain representativeness are selected based on the constraint of the lowest global elevation point and linear fitting. Then, based on TIN cluster semantic analysis and denoising, the optimized seed point set is obtained; Step S3: Iteratively encrypt and optimize the triangulation, and output the encrypted triangulation and the optimized seed point set; Step S4: Construct a ground DEM based on grid filling, and perform bilinear interpolation filtering based on the ground DEM to calculate the difference between the actual elevation of the point cloud and the elevation of the DEM. If the difference is less than the final distance threshold, it is marked as a ground point, thus completing the full point cloud ground classification.

2. The ground point cloud filtering method based on curvature dynamic segmentation according to claim 1, characterized in that, Step S1 includes the following steps: Step S11: Adaptive downsampling based on point cloud density; if the density of the original point cloud is greater than the density threshold, random downsampling is performed; otherwise, the original point cloud is directly retained. Step S12: Divide the point cloud into L×L blocks evenly to generate several point cloud blocks.

3. The ground point cloud filtering method based on curvature dynamic segmentation according to claim 1, characterized in that, Step S2 includes the following steps: Step S21: Preliminary screening of potential ground points; First, the point cloud blocks are divided into grids, the lowest elevation point in each grid and the global lowest elevation point of the point cloud block are extracted, and the angle α between the line connecting the lowest elevation point and the global lowest elevation point and the horizontal plane is determined. If α is less than the set angle threshold, it is selected as a potential ground point. Step S22: Project the potential ground points in the X / Y directions and calculate their curvature. Recursively divide the potential ground points into sub-blocks according to the curvature threshold, extract the lowest elevation point in the final sub-block, and form an initial seed point set. Step S23: Denoise the initial seed point set based on TIN cluster semantic analysis; construct an initial triangular network (TIN) based on the initial seed point set, traverse all triangular units in the initial TIN, filter out suspected triangular units by the slope or maximum elevation difference of the triangles, and cluster them to obtain suspected triangular clusters; then, combine semantic features to remove high-position noise points and building-type pseudo-ground points. Step S24: Based on the denoised seed point set, perform two-dimensional projection and linear fitting to remove residual noise points that deviate from the ground trend, and obtain the optimized seed point set.

4. The ground point cloud filtering method based on curvature dynamic segmentation according to claim 3, characterized in that, Step S22 includes the following steps: Step A1: Using the center coordinates (X_center, Y_center) of the point cloud block of the potential ground points as the reference, extract the point cloud subsets [X_center-L / 2, X_center+L / 2] and [Y_center-L / 2, Y_center+L / 2] in the X-axis and Y-axis directions respectively; Step A2: Project the point cloud subset in the X direction onto the XZ plane, project the point cloud subset in the Y direction onto the YZ plane, and sort them in ascending order of X / Y coordinates to obtain a continuous sequence of Z values ​​as X / Y change. Step A3: Calculate the terrain curvature κ in the X-axis and Y-axis directions for each point using the numerical differentiation method. x and κᵧ; If κ x If the curvature threshold κ0 is greater than the curvature threshold κ0, then the X coordinate of the point is marked as the X-direction dividing line; if κᵧ > curvature threshold κ0, then the Y coordinate of the point is marked as the Y-direction dividing line; the current point cloud is divided into several sub-blocks according to the dividing lines; Step A4: Repeat steps A2 and A3 to dynamically recursively divide each sub-block until a stopping condition is met. The stopping condition is: κ of all points within the sub-block. x ≤κ0 and κᵧ≤κ0; or the size of the sub-block in the X / Y direction is less than or equal to the preset minimum size S0; Step A5: Extract the lowest elevation point within each sub-block and summarize them into an initial seed point set.

5. A ground point cloud filtering method based on curvature dynamic segmentation according to claim 3, characterized in that, In step S23, if all triangles in a suspected triangle cluster share a common vertex and the triangles are connected end to end, then the suspected triangle cluster is determined to be a high-noise cluster, its common vertex is a high-noise point, and the high-noise point is removed from the initial seed point set. If the architectural semantic features of a suspected triangle cluster satisfy the following conditions: ① the number of triangles in the cluster is greater than or equal to the minimum number of triangles N1 in the architectural triangle cluster; ② the average angle between the normal vector of all triangles in the cluster and the vertical direction is greater than 60°; ③ the aspect ratio of the minimum bounding box of the cluster is greater than the minimum aspect ratio λ of the minimum bounding box of the architectural triangle cluster and the area of ​​the bounding box is greater than or equal to the minimum area S1 of the minimum bounding box of the architectural triangle cluster; ④ the ratio of the total area of ​​triangles in the cluster to the area of ​​the cluster's convex hull is less than the threshold μ; then the suspected triangle cluster is determined to be an architectural pseudo-ground cluster, and all seed points in the cluster are removed from the initial seed point set.

6. The ground point cloud filtering method based on curvature dynamic segmentation according to claim 3, characterized in that, Step S24 includes the following steps: Step B1: Project the semantically denoised seed points onto the XY plane, divide them into segments with an interval width I along the X-axis, extract the lowest elevation point in each segment, and form a core feature point set. Step B2: Calculate the slope K between two adjacent points in the core feature point set, mark the locations where the absolute value of the slope is greater than the slope threshold K1 and the difference between adjacent slopes is greater than the slope difference threshold K2 as terrain change points, and remove the corresponding high-position noise points; Step B3: Fit a straight line to adjacent point pairs in the core feature point set, calculate the vertical distance d from all points between the two points to the straight line, mark the points whose d > distance threshold D3 and remove them; Step B4: Divide the denoised 2D point cloud into blocks twice according to the preset step size S. For sub-segments with a number of points greater than or equal to the minimum number of points N4, fit a straight line using the Tukey robust least squares method. If the absolute value of the slope of the fitted line is greater than the terrain steepness threshold K3, mark outlier noise points according to the 3σ criterion. If the absolute value of the slope is less than the terrain steepness threshold K3, retain all points in the sub-segment. Map the points marked as retained back to the 3D point cloud through indexing to obtain the optimized seed point set.

7. A ground point cloud filtering method based on curvature dynamic segmentation according to any one of claims 3-6, characterized in that, Step S3 includes the following steps: Step S31: Traverse all point cloud blocks and mark point cloud blocks with a seed point count greater than 0 as valid blocks; Step S32: Based on the effective blocks, extract the X / Y coordinates of the seed points and construct a KD-Tree. Search the points in the neighborhood of each point, retain the point with the lowest elevation in the neighborhood as the seed point, merge the seed points of the current effective block and the eight neighboring effective blocks to generate an extended seed point set, and construct the Delaunay triangulation of the extended seed point set. Step S33: Calculate the XYZ bounding box of the potential ground point cloud within the current valid block, divide the X / Y plane grid, traverse all triangles of the triangulation, map the triangle index to the corresponding grid and remove duplicates, generate a global mapping table of grid-triangle, and match potential ground points with triangulation cells. Step S34: Determine the ground attributes of potential ground points based on the slope of the triangular network; traverse each potential ground point that matches a triangular network unit, calculate the normal vector and slope of the triangular plane, and mark the potential ground points corresponding to the triangular planes that meet the threshold range as ground points; Step S35: Identify pseudo-ground triangles, remove the vertex with the highest elevation, and add the potential ground point farthest from the triangle below it as a real ground seed point, and update the seed point set; Step S36: Repeat steps S32 to S35 until the iteration termination condition is met, and output the final encrypted triangular network and the optimized seed point set.

8. A ground point cloud filtering method based on curvature dynamic segmentation according to claim 7, characterized in that, In step S34, if the approximate planar slope threshold θ3 < the slope of the triangular plane ≤ the steep terrain slope threshold θ2, and if the directed distance d1 from the potential ground point to the triangular plane < the dynamic distance threshold d... 10 If the angle between the potential ground point and the triangular plane is less than the dynamic angle threshold θ4, then the potential ground point is marked as a ground point. If the slope of the triangular plane is ≤ θ3, then the potential ground points below the triangular plane are directly marked as ground points, and the points located above the triangular plane and inside the triangle with the required angle are marked as ground points. If the slope of the triangular plane is greater than the steep terrain slope threshold θ2, then the highest vertex of the triangle is identified, a mirror point of the potential ground point is generated, and the distance d2 from the mirror point to the triangular plane and the included angle α2 are calculated. If d2 < dynamic distance threshold d 10 If α2 < dynamic angle threshold θ4, then the original potential ground point is marked as a ground point.

9. A ground point cloud filtering method based on curvature dynamic segmentation according to claim 1, characterized in that, Step S4 includes the following steps: Step S41: Calculate the full point cloud bounding box, initialize the DEM grid according to the preset resolution, traverse the ground seed points, assign each seed point to the corresponding grid according to the X / Y coordinates, extract the lowest elevation point in each grid as the initial DEM grid point; mark the grid with valid points as filled, and mark the grid without data as to be filled. Step S42: Collect the coordinates of the blank grid, search for the eight-direction valid grid of the blank grid in stages, fill it with linear interpolation and verify the validity of the Z value; summarize the filled grid points, filter invalid Z values ​​and generate an irregular DEM grid. Step S43: Convert the irregular DEM grid into a regular DEM grid; traverse the valid nodes of the irregular DEM grid, and calculate the Z value of the regular DEM grid node based on the concavity and convexity of the quadrilateral formed by the four neighbors of the DEM grid; for convex quadrilaterals, calculate the Z value of the nearest neighbor point using bilinear interpolation; for non-convex quadrilaterals, calculate the Z value of the nearest neighbor point, and finally generate a complete regular DEM grid. Step S44: Perform bilinear interpolation based on the regular DEM grid to calculate the difference between the actual elevation of the point cloud and the DEM elevation. If the difference is less than the final distance threshold, it is marked as a ground point, and the ground classification is completed.

10. A computer-readable storage medium having a computer program stored thereon, characterized in that, When executed by a processor, the program implements a ground point cloud filtering method based on curvature dynamic segmentation as described in any one of claims 1-9.

Images