Automatic ground point filtering method for laser radar point cloud

Abstract

The invention provides a ground point automatic filtering method for laser radar point cloud, and relates to the technical field of point cloud data processing and terrain digital modeling. The method comprises the steps of carrying out voxelization processing on point cloud, constructing a base centroid point layer, generating a near-earth digital elevation model, extracting and optimizing non-ground pattern spots, constructing a semantic centroid point layer and establishing a KD-tree index to extract ground points. Through combination of maximum value filtering and pattern spot elevation feature judgment, accurate identification and error correction of non-ground pattern spots are realized, a centroid point layer is further constructed, ground points are marked by adopting a bidirectional spatial index structure, and the filtering precision and robustness are improved. The method is suitable for automatic ground point extraction under complex ground surface conditions, can support subsequent building instance segmentation, height measurement and structural analysis, and has wide engineering application value.

Description

Technical Field

[0001] The present invention relates to the technical field of point cloud data processing and terrain digital modeling, and in particular to an automatic ground point filtering method for laser radar point clouds. Background Art

[0002] With the advancement of laser scanning technology, multi-temporal laser pulse signals reflected from the ground or target surface are accurately captured and decoded into a large number of discrete point clouds, effectively depicting the three-dimensional structural features and surface morphology of the scanned target. Point clouds offer unique advantages in depicting surface information, especially in areas with complex terrain and multi-layered vegetation. Point clouds not only provide detailed characterization of the microtopography beneath the forest but also accurately reflect the vertical structure of the forest canopy. This is of great significance for high-precision topographic mapping of vegetated areas, topographic feature extraction, and forest canopy biomass assessment. The spatial distribution of ground point clouds clearly depicts topographic information such as slope and aspect, providing a scientific basis for path planning for agricultural machinery operations. They can also be used to identify suitable construction areas, assess land development potential, especially in complex mountainous and hilly terrain, and support the selection of orchard planting sites. However, because point clouds in complex scenes often contain a mixture of ground points and non-ground points (such as buildings, trees, and power lines), distinguishing these two types of points is a technical challenge. Accurately identifying and extracting ground points is crucial for high-precision terrain modeling, while extracting structural features from non-ground points is also crucial for 3D object reconstruction.

[0003] Ground filters (GFs) can be used to distinguish ground points from non-ground points in a point cloud. GFs typically distinguish between the two types of point clouds based on the fact that laser pulses with common bandwidths cannot penetrate soil layers or opaque building surfaces. Therefore, the ground point clusters extracted using GFs are typically located in the lower region of a local longitudinal profile of the point cloud, encompassing the ground, underlying vegetation, or building surfaces. Numerous different GF algorithms have been proposed, which can be categorized into two types based on their operating principles: rule-based and machine learning-based. Each approach has its own advantages and disadvantages. While meeting the accuracy requirements of GF results, rule-based GFs are more popular because they do not require significant training costs. Rule-based GF designs can be further categorized into slope thresholding, morphological filtering, spatial interpolation functions, and "cloth" filtering. However, existing rule-based GFs have their own limitations, their robustness needs to be improved, and automated operation is not feasible.

[0004] Due to the limited penetration of lasers into dense vegetation layers, point clouds in dense forests and densely built-up areas are prone to missing ground points. During ALS data acquisition, smaller scanning angles and sparse flight path settings can easily lead to incomplete depictions of the ground near the sides and eaves of buildings in the point cloud, increasing GF errors. In addition, sensor errors and environmental factors (such as wind) can cause low-value noise points to appear below the ground, or high-value noise points to appear near the sensor's flight altitude. Therefore, GF must be robust and accurate enough to cope with point cloud data acquired by scanning different scenes. When optimizing rule-based GF, reducing the number of parameters and enhancing its robustness in different scenarios are key. Designing a GF that is accurate, easy to use, and robust in complex scenes is a valuable but challenging task. Summary of the Invention

[0005] To overcome the shortcomings of the existing technology, this invention aims to provide an automated ground point filtering method for LiDAR point clouds, enabling accurate and robust automated extraction of ground points from diverse scanning scenarios. Furthermore, the intermediate results of this GF analysis process can also assist in instance segmentation of building point clouds, supporting building height measurement and morphological feature extraction.

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

[0007] An automated ground point filtering method for a lidar point cloud, comprising:

[0008] Perform voxel processing on the input LiDAR point cloud, determine the voxel size according to the sampling statistics, remove empty voxels in the LiDAR point cloud according to the voxel size, and extract the centroid point set of non-empty voxels;

[0009] Based on the centroid point set of the non-empty voxels, a method of iteratively adjusting the centroid point elevation is adopted to construct a base centroid point layer, and using the centroid point set of the non-empty voxels as a reference, a near-geocentroid point set consisting of nearest neighbors is obtained in the centroid point set of the non-empty voxels, and a digital elevation model of the near-geocentroid point is constructed based on the near-geocentroid point set;

[0010] Performing maximum filtering and threshold binarization processing on the digital elevation model to generate an initial non-ground coverage area binarization map, and extracting the maximum elevation of the edge of the patch and the average elevation of the filled area based on the elevation features of the patch in the initial non-ground coverage area binarization map, and extracting the non-ground patch based on the comparison result of the maximum elevation of the edge of the patch and the average elevation of the filled area;

[0011] Performing morphological closing operations and edge pixel completion processing on the non-ground patches, and correcting misidentified patches and edge missing patches through intra-patch pixel consistency and elevation continuity detection to form an optimized non-ground patch mask;

[0012] Converting the non-ground patch mask into a set of centroid points containing ground and non-ground semantics, marking the centroid points in the set of centroid points containing ground and non-ground semantics as ground centroid points and non-ground centroid points according to pixel values, identifying and completing missing points in the ground centroid points, and generating a centroid point layer for supporting ground centroid point extraction;

[0013] A KD-tree index structure is constructed for the lidar point cloud and the centroid point layer respectively. In the three-dimensional coordinate space, the centroid point in the centroid point layer is used as a reference, and the nearest neighbor centroid point of each point in the input lidar point cloud is obtained using the KD-tree index. If the nearest neighbor centroid point belongs to the ground centroid point, the point in the input lidar point cloud corresponding to the nearest neighbor centroid point is marked as a ground point.

[0014] Preferably, the voxel size is obtained by calculating the average Euclidean distance between 1% of randomly selected point pairs in the input lidar point cloud.

[0015] Preferably, when constructing the base centroid point layer, the elevation value of each centroid point is adjusted in an iterative manner; the steps for adjusting the elevation value of each centroid point are: selecting neighboring centroid points with elevations lower than the centroid point within the columnar neighborhood of the centroid point, calculating the median value of the neighboring centroid points, and using the median value to update the elevation of the current centroid point until the nearest neighbor spatial index is stable.

[0016] Preferably, the digital elevation model is a two-dimensional elevation grid constructed based on the near-geocentric point set and in accordance with spatial resolution.

[0017] Preferably, performing maximum filtering and threshold binarization processing on the digital elevation model to generate an initial non-ground coverage area binarization map includes:

[0018] The digital elevation model is subjected to maximum filtering and threshold binarization processing to generate a binary map marking elevation mutation areas; the processing formula of the binary map is:

[0019]

[0020] Among them, R HC For the binary image, NGNEM max is the maximum filtered digital elevation model, NGNEM is the digital elevation model, and W is the voxel size;

[0021] Grid filling and elevation judgment are performed based on the spots in the binary map to generate an initial non-ground coverage area binary map.

[0022] Preferably, based on the patch elevation features in the initial non-ground coverage area binary map, extracting the maximum elevation of the patch edge and the average elevation of the filling area, and extracting the non-ground patch according to the comparison result of the maximum elevation of the patch edge and the average elevation of the filling area, comprises:

[0023] Perform unique value marking on each patch in the binary image;

[0024] Perform grid filling operation on the marked spots to obtain the complete internal filling area;

[0025] Calculate the maximum elevation of the patch edge and the average elevation of the filling area;

[0026] If the maximum elevation of the patch edge is greater than the average elevation of the filling area, the patch is determined to be a non-ground patch in the initial non-ground coverage area binarization map.

[0027] Preferably, the step of correcting misidentified spots includes:

[0028] If a patch is marked as a non-ground patch before optimization, but the corresponding areas of all pixels in the elevation model after marking are not completely included in the non-ground patch mask, and the elevation of the patch has a continuous change trend in space with its adjacent ground patches, then the patch will be re-marked as a ground patch.

[0029] Preferably, when generating a centroid point layer, for missing ground centroid points, the average elevation value of the 24 nearest marked ground centroid points is used as the complementary elevation value to construct a continuous centroid point layer.

[0030] Preferably, the KD-tree index structure is constructed based on the lidar point cloud and the centroid point layer respectively, and is used to obtain the nearest neighbor centroid point of each input point in the three-dimensional coordinate space; if the nearest neighbor centroid point belongs to the ground centroid point, the corresponding input point is marked as a ground point.

[0031] Preferably, before constructing the KD-tree index structure of the centroid layer, the elevation value of each centroid in the centroid layer is uniformly adjusted down by 2 times the voxel size.

[0032] According to the specific embodiments provided by the present invention, the present invention discloses the following technical effects:

[0033] The present invention provides a ground point automatic filtering method for a laser radar point cloud, comprising: voxelizing the input laser radar point cloud, determining the voxel size according to sampling statistics results, and removing empty voxels in the laser radar point cloud according to the voxel size, and then extracting a centroid point set of non-empty voxels; based on the centroid point set of the non-empty voxels, adopting an iterative adjustment method of the centroid point elevation to construct a base centroid point layer, and using the centroid point set of the non-empty voxels as a reference, obtaining a near-geocentroid point set consisting of the nearest neighbors in the centroid point set of the non-empty voxels, and constructing a digital elevation model of the near-geocentroid point based on the near-geocentroid point set; performing maximum filtering and threshold binarization on the digital elevation model to generate an initial non-ground coverage area binarization map, and based on the patch elevation features in the initial non-ground coverage area binarization map, extracting the maximum elevation of the patch edge and the average elevation of the filling area, and extracting the non-ground coverage area based on the comparison result of the maximum elevation of the patch edge and the average elevation of the filling area. The invention relates to a method for extracting a non-ground patch; performing morphological closing operation and edge pixel completion processing on the non-ground patch, and correcting misidentified patches and edge missing patches through pixel consistency and elevation continuity detection within the patch to form an optimized non-ground patch mask; converting the non-ground patch mask into a centroid point set containing ground and non-ground semantics, marking the centroid points in the centroid point set containing ground and non-ground semantics as ground centroid points and non-ground centroid points according to pixel values, identifying missing points in the ground centroid points and completing them to generate a centroid point layer for supporting ground centroid point extraction; constructing a KD-tree index structure for the lidar point cloud and the centroid point layer respectively, and in the three-dimensional coordinate space, using the centroid point in the centroid point layer as a reference, using the KD-tree index to obtain the nearest neighbor centroid point of each point in the input lidar point cloud; if the nearest neighbor centroid point belongs to the ground centroid point, the point in the input lidar point cloud corresponding to the nearest neighbor centroid point is marked as a ground point. The present invention can accurately and robustly extract ground points from a variety of scanned scenes. At the same time, the intermediate results of the GF analysis process can also assist in the instance segmentation of building point clouds to support building height measurement and morphological feature extraction. BRIEF DESCRIPTION OF THE DRAWINGS

[0034] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the following briefly introduces the drawings required for use in the embodiments. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying any creative work.

[0035] Figure 1 A flow chart of a method provided by an embodiment of the present invention;

[0036] Figure 2 A technical roadmap of the present invention provided for embodiments of the present invention;

[0037] Figure 3 A schematic diagram of constructing a digital elevation model of a point cloud base layer and a near-geocenter point provided by an embodiment of the present invention;

[0038] Figure 4 A schematic diagram of a ground coverage area identification process provided by an embodiment of the present invention;

[0039] Figure 5 A schematic diagram of the ground point extraction process provided by an embodiment of the present invention. DETAILED DESCRIPTION

[0040] The following will clearly and completely describe the technical solutions in the embodiments of the present invention in conjunction with the accompanying drawings. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without making creative efforts are within the scope of protection of the present invention.

[0041] The purpose of this invention is to provide an automated ground point filtering method for LiDAR point clouds, capable of accurately and robustly extracting ground points from diverse scanning scenarios. Furthermore, the intermediate results of this ground point filtering method during analysis can also assist in instance segmentation of building point clouds, supporting building height measurement and morphological feature extraction.

[0042] In order to make the above-mentioned objects, features and advantages of the present invention more obvious and easy to understand, the present invention is further described in detail below with reference to the accompanying drawings and specific embodiments.

[0043] Figure 1 A flow chart of the method provided in the embodiment of the present invention is shown in FIG. Figure 1 As shown, the present invention provides a method for automatically filtering ground points for a lidar point cloud, comprising:

[0044] The input LiDAR point cloud is voxelized, the voxel size is determined based on the sampling statistics, and the empty voxels in the LiDAR point cloud are removed based on the voxel size, and the centroid point set of non-empty voxels is extracted;

[0045] Based on the centroid point set of non-empty voxels, the base centroid point layer is constructed by iteratively adjusting the centroid point elevation. The centroid point set of non-empty voxels is used as a reference to obtain the near-geocentroid point set consisting of the nearest neighbors in the centroid point set of non-empty voxels, and the digital elevation model of the near-geocentroid point is constructed based on the near-geocentroid point set.

[0046] Perform maximum filtering and threshold binarization processing on the digital elevation model to generate an initial non-ground coverage area binarization map. Based on the patch elevation characteristics in the initial non-ground coverage area binarization map, extract the maximum elevation of the patch edge and the average elevation of the filling area. Extract the non-ground patch based on the comparison result of the maximum elevation of the patch edge and the average elevation of the filling area.

[0047] Perform morphological closing operations and edge pixel completion on non-ground patches, and correct misidentified patches and edge missing patches through intra-patch pixel consistency and elevation continuity detection to form an optimized non-ground patch mask.

[0048] The non-ground patch mask is converted into a set of centroid points containing ground and non-ground semantics. The centroid points in the set of centroid points containing ground and non-ground semantics are marked as ground centroid points and non-ground centroid points according to pixel values. The missing points in the ground centroid points are identified and completed to generate a centroid point layer to support ground centroid point extraction.

[0049] KD-tree index structures are constructed for the lidar point cloud and centroid point layer respectively. In the three-dimensional coordinate space, the centroid point in the centroid point layer is used as a reference to obtain the nearest neighbor centroid point of each point in the input lidar point cloud using the KD-tree index. If the nearest neighbor centroid point belongs to the ground centroid point, the point in the input lidar point cloud corresponding to the nearest neighbor centroid point is marked as a ground point.

[0050] like Figures 2 to 5 As shown in FIG, after voxelizing the input point cloud, the present invention extracts near-ground voxels and their centroid point sets to depict ground fluctuation trends. Then, the present invention applies a composite morphological filter to mark areas covered by non-ground objects in the centroid point digital elevation model. Finally, the present invention extracts ground points from the input point cloud using a bidirectional k-dimensional tree (KD-tree) spatial engine constructed between the centroid point set and the input point cloud. The specific invention content is as follows:

[0051] (1) Point cloud voxelization and near-geocenter point set extraction

[0052] LWF first voxelizes the input point cloud, which can effectively reduce the computational cost of point-by-point analysis. The voxel size is denoted as W. A subset of points that occupies 1% of the number of points in the input point cloud is randomly selected to determine W. Then, the average Euclidean distance between each selected point and its nearest neighbor is calculated. In order to reduce the number of empty voxels and keep the structural details of the point cloud voxelization results as much as possible, the parameter ablation experiment was performed to determine W as

[0053]

[0054] After removing empty voxels, the centroid of each non-empty voxel is extracted to reduce the computational cost of subsequent analysis while preserving the rough structure and distribution pattern of the point cloud.

[0055] Since the multi-echo laser pulse can penetrate the plant canopy but cannot penetrate the ground, the voxel centroid point set (ND CE ) to extract the centroid subset of near-ground voxels one by one. Considering that laser pulses in the green, near-infrared, and short-wave infrared bands cannot penetrate buildings and have limited penetration into dense forest canopies, the near-ground centroid point set in areas covered by buildings and plants is always missing in the ND CE In addition, some ND covering the roof and canopy surface CE It is easy to be mistaken for a point near the ground centroid. Some NDs with abnormal heights CE This is due to the presence of noise points in the input point cloud. Removing point cloud noise and misclassified near-ground centroid points is crucial for the accurate operation of the present invention.

[0056] (2) Construction of the base layer of the point cloud coverage area

[0057] Extract ND from input point cloud CE Then gradually reduce ND CE The elevation of ND is used to construct the base centroid point layer of the input point cloud to support the extraction of ground points. CE Construct a two-dimensional KD-tree spatial index with a search radius set to To obtain the local centroid subset C along the z-axis cylindrical distribution of the centroid G p . Compare C p The elevations of other centroid points in are compared with the elevation of G. If C p If there are centroids whose elevation is lower than G, calculate the median elevation of these centroids and assign it to the elevation of G. CE to form the base layer of the input point cloud. p In the example, the spatial index of the nearest neighbor centroid of the target centroid G is different between the N+1 and N iterations, then according to C p The median elevation of the lower point subset in G is adjusted until the two indices are the same. Finally, the base centroid point layer ND of the input point cloud can be obtained. BL In addition to elevation, ND BL The (x,y) coordinates of each point in ND CE same.

[0058] For ND BL The target centroid point Q in the present invention uses a three-dimensional space KD-tree search engine in ND CE Find its nearest neighbor centroid point NQ , to obtain ND CE The near-geocenter point set (ND GD ), ND GD It is ND CE Subset of ND GD There are still some centroid points from the roof and forest canopy surface. In order to remove these error points, based on ND CE The elevation information of the ground is used to build a near-geocenter digital elevation model (NGNEM) with a resolution of Wm, and the elevation changes described by the NGNEM are used to drive the detection of ND GD The non-ground centroid point in . Here is based on ND CE To build NGNEM is to prevent ND GD Only the misclassified centroid points that are partially covered on the roof and tree canopy surfaces are retained, which will reduce the completeness of the subsequent identification of non-ground cover areas from NGNEM.

[0059] (3) Ground coverage area identification

[0060] In medium and high density LiDAR point clouds, the height difference between ground points and roof points usually exceeds W. The local point cloud elevation change trend near cliffs, ridges, and steep slopes is continuous. This continuous terrain change property provides the basis for dividing the ground and non-ground areas (e.g., roofs and trees) in NGNEM. The present invention uses a maximum filter with an analysis window size of 5×5 to analyze NGNEM and obtains NGNEM max Based on this, we use formula (2) to calculate NGNEM max Binarize and get R HC .

[0061]

[0062] The edge of non-ground objects is in R HC Therefore, we use the Union-Find algorithm to find the HC The spots in the image are marked with unique values, which is R HC Each independent patch shown in is assigned a unique number. Then, for a labeled patch P L Perform grid filling analysis. If P L Can be filled, then calculate the maximum elevation of the edge of the patch (H edge ) and the average elevation of the internal filling area (H fill ). If H fill ≥H edge , then we can judge P L belongs to the area not covered by ground objects. Otherwise, P L Located on steep slopes or terraced areas with continuous elevation changes.

[0063] Since the elevation information described by NGNEM may contain some misclassified ND from roofs or tree crowns, GD It is difficult to detect the spots where these error centroid points are located by only using the hole filling algorithm. To solve this problem, the morphological closing operation is applied to analyze the marked spots in NGNEM to obtain the NGNEM BI The window size w used for morphological closing analysis close for:

[0064]

[0065] Among them, N close The number of pixels filled in each fillable patch. BI In order to avoid some non-ground patches near the edge of NGNEM from being reflected in the NGNEM, the pixels are marked as “1” and the rest are marked as “0”. BI Lost in NGNEM BI The pixel values ​​of the first row, first column, last row, and last column are temporarily reset to zero. The pixel values ​​of these edge rows and columns are set to "1", "2", "3", and "4" respectively, so that NGNEM edge Then, calculate formula (4) to get ED P .

[0066]

[0067] If ED P If a patch has less than two non-zero values ​​inside the pixel, it is considered as a non-ground patch located at the edge of the NGNEM.

[0068] Merge these non-ground patches located at the edge of NGNEM into NGNEM BI middle.

[0069] In NGNEM BI In the NGNEM, patches on slopes may be mistakenly identified as non-ground types. To solve this problem, the present invention labels patches in NGNEM with unique values ​​based on elevation changes. In the labeling results, areas at the same elevation in NGNEM will form multiple patches. If all pixels within the coverage area of ​​a target patch in the labeled NGNEM are in the same elevation as the target patch, the target patch will be marked as non-ground type. BI If all pixels in the image are "1", they are classified as non-ground patches in NGNEM. Otherwise, the patch is reassigned to "0" in the marked NGNEM, that is, it is considered a ground patch. After binarization of the NGNEM with the patch unique value, the NGNEM can be obtained. NG .

[0070] Next, morphological dilation analysis was applied to NGNEM. NG The NGNEM NG The morphological expansion results are consistent with those of NGNEM. NG By subtracting the pixels from each other, we can obtain the outer boundary pixel set of each patch. Each pixel boundary set will be used to extract the elevation change near the edge of a single patch from the NGNEM. For a patch, if the maximum elevation of its edge pixel coverage area exceeds the average elevation of its coverage area, it will be classified as a ground patch on a slope. After identifying the ground patch on the slope, the present invention completes the NGNEM. NG Optimization.

[0071] (4) Preparation for ground point extraction

[0072] According to NGNEM NG The row and column index, elevation information and spatial resolution W of the pixel in the image are used to convert the center point of each pixel back to the spatially discrete centroid point. In the conversion result, the type of the centroid point is based on the NGNEM NG The pixel value in is determined as follows: the centroid point TP with a value of "1" N The centroid point TP with a value of "0" in the non-ground patch G Located in the ground patch. Each TP N and TP G The elevation of the point is given by ND CE In order to improve the accuracy of GF, ND is also used. CE The (x,y) index of the centroid identifies the ND GD The missing (x,y) index is used to determine the centroid points that need to be added. For each newly added centroid point, its elevation is assigned to its 24 nearest neighbor TPs. G Based on this, a centroid point layer is constructed below the area covered by the input point cloud.

[0073] (5) Ground point extraction based on bidirectional spatial KD-tree index

[0074] A KD-tree spatial index (OKDT) is built for the input point cloud, and another KD-tree spatial index (TKDT) is built for the near-ground centroid point layer to form a bidirectional KD-tree spatial index. Before this step, the elevation of the near-ground centroid point set is reduced by 2W to avoid extracting point clouds from the surfaces of small non-ground objects (such as vehicles parked on the road). For each point in the OKDT index, its nearest centroid point is searched. Then, the number of points in the input point cloud that share the same nearest neighbor centroid point is counted, and their average value (M, M≥1) is calculated. For each TP in the input point cloud GThe node searches for M neighboring points. If the nearest neighbor centroid of a point in the input point cloud belongs to TP N , then it is located in a non-ground coverage area. G M points of the node, if the nearest centroid of a point belongs to TP G , record its index in OKDT to obtain a list G index According to G index , ground points can be extracted from the input point cloud.

[0075] As an optional implementation, this embodiment uses a drone or an aerial laser radar scanner to scan and obtain discrete point clouds within the route coverage area according to the set flight path and scanning angle. The point cloud of the target area is cropped, and the cropped result is input into the ground point filter proposed by the present invention. The point cloud is voxelized, and the empty voxels after the point cloud voxelization are removed. The size of the voxel is set to W, and the centroid of each non-empty voxel is extracted to reduce the computational cost of the subsequent analysis process. After extracting the centroid point set of the voxel, the elevation of the centroid point is gradually lowered to construct the base layer of the input point cloud. It is ND CE Construct a two-dimensional KD-tree spatial index with a search radius set to To obtain the cylindrical distribution of local centroid subset C along the z axis of each centroid point p . Compare C p The elevation of other centroid points in the target centroid. If C p If there are centroids whose elevation is lower than the target centroid, calculate the median elevation of these centroids and assign it to the elevation of G. CE to form the base layer of the input point cloud. p In the analysis, the spatial index of the nearest neighbor centroid of the target centroid is different between the N+1 and N iterations. p The median elevation of the lower point subset in G is adjusted until the two indices are the same. Finally, the base centroid point layer ND of the input point cloud is obtained. BL .

[0076] For ND BL The target centroid point Q in ND is obtained by using the three-dimensional KD-tree search engine. CE Find its nearest neighbor centroid point N Q , to obtain ND CE The near-geocenter point set (ND GD ). To remove ND GD Some of the remaining centroid points from the roof and forest canopy surfaces are based on ND CEThe elevation information of the geocentric point is used to construct a near-geocentric digital elevation model (NGNEM) with a resolution of Wm. The NGNEM is analyzed using a maximum filter with a window size of 5×5 to obtain the NGNEM. max About NGNEM max Binarize and get R HC The edges of non-ground objects are in R HC was outlined.

[0077] Using Union-Find algorithm to analyze R HC The spots in the image are marked with unique values, which is R HC Each individual patch shown in is assigned a unique number. Then, for a labeled patch P L Perform grid filling analysis. L The maximum elevation of the edge of the patch (H edge ) and the average elevation of the internal filling area (H fill ) size, judge P L Belongs to non-ground objects or ground coverage areas. fill <H edge , P L Located on a steep slope or terraced area with continuous elevation changes. Otherwise, P L Located in an area not covered by ground objects.

[0078] Apply morphological closing operation to analyze the marked patches in NGNEM to obtain NGNEM BI Closed patches of non-ground objects in NGNEM BI In order to avoid some non-ground patches near the edge of NGNEM from being reflected in the NGNEM, the pixels are marked as “1” and the rest are marked as “0”. BI Lost in NGNEM BI The pixel values ​​of the first row, first column, last row, and last column are temporarily reset to zero. The pixel values ​​of these edge rows and columns are set to "1", "2", "3", and "4" respectively, so that NGNEM edge , based on which ED is calculated P If ED P If a pixel in a map has less than two non-zero values, it is considered as a non-ground map located at the edge of NGNEM. These non-ground maps located at the edge of NGNEM are merged into NGNEM. BI middle.

[0079] To reduce NGNEM BIThe present invention labels the patches in NGNEM with unique values ​​based on elevation changes, so that the areas at the same elevation in NGNEM form multiple patches. If all pixels within the coverage area of ​​a target patch in the marked NGNEM are in the same elevation in NGNEM, BI If all pixels in the NGNEM are “1”, they are classified as non-ground patches in the NGNEM. Otherwise, the patches are reassigned to “0” in the marked NGNEM, that is, they are regarded as ground patches. After binarization of the marked NGNEM, the NGNEM is obtained. NG Applying morphological expansion analysis to NGNEM NG and compared with NGNEM NG By subtracting the pixels from each patch, we can obtain the outer boundary pixel set of each patch. Each pixel boundary set will be used to extract the elevation change near the edge of a single patch from the NGNEM. For a patch, if the maximum elevation of its edge pixel coverage area exceeds the average elevation of its coverage area, it will be classified as a ground patch on a slope. After identifying the ground patches on the slope, the NGNEM is completed. NG optimization.

[0080] According to NGNEM NG The row and column index, elevation information and spatial resolution W of the pixel in the image are used to convert the center point of each pixel back to the spatially discrete centroid point. In order to improve the accuracy of GF, ND is also used. CE The (x,y) index of the centroid identifies the ND GD The missing (x,y) index in the ,is used to determine the centroid point that needs to be added. For each newly added node, its elevation is assigned to its 24 nearest neighbor TPs G The average elevation of the points. Based on this, a centroid point layer is constructed below the area covered by the input point cloud. A KD-tree spatial index (OKDT) is constructed for the input point cloud, and another KD-tree spatial index (TKDT) is constructed for the centroid point layer to form a bidirectional KD-tree spatial index. Before this step, the centroid point elevation of the centroid point layer is reduced by 2W to avoid extracting point clouds from the surfaces of small non-ground objects (such as vehicles parked on the road).

[0081] For each point in the OKDT index, search for its nearest centroid. Then, count the number of points in the input point cloud that share the same nearest neighbor centroid and calculate their average (M, M≥1). G The node searches for M neighboring points. If the nearest neighbor centroid of a point in the input point cloud belongs to TP N , then it is located in a non-ground coverage area. GM points of the node, if the nearest centroid of a point belongs to TP G , record its index in OKDT to obtain a list G index According to G index , and finally extract ground points from the input point cloud.

[0082] The various embodiments in this specification are described in a progressive manner, and each embodiment focuses on the differences from other embodiments. The same or similar parts between the various embodiments can be referenced to each other.

[0083] This document uses specific examples to illustrate the principles and implementation methods of the present invention. The above examples are only intended to help understand the method and core concept of the present invention. At the same time, those skilled in the art will find that the specific implementation methods and application scopes may vary based on the concept of the present invention. In summary, the contents of this specification should not be construed as limiting the present invention.

Claims

1. A method for automatic filtering of ground points in lidar point clouds, characterized in that: include: Perform voxel processing on the input LiDAR point cloud, determine the voxel size according to the sampling statistics, remove empty voxels in the LiDAR point cloud according to the voxel size, and extract the centroid point set of non-empty voxels; Based on the centroid point set of the non-empty voxels, a method of iteratively adjusting the centroid point elevation is adopted to construct a base centroid point layer, and using the centroid point set of the non-empty voxels as a reference, a near-geocentroid point set consisting of nearest neighbors is obtained in the centroid point set of the non-empty voxels, and a digital elevation model of the near-geocentroid point is constructed based on the near-geocentroid point set; Performing maximum filtering and threshold binarization processing on the digital elevation model to generate an initial non-ground coverage area binarization map, and extracting the maximum elevation of the edge of the patch and the average elevation of the filled area based on the elevation features of the patch in the initial non-ground coverage area binarization map, and extracting the non-ground patch based on the comparison result of the maximum elevation of the edge of the patch and the average elevation of the filled area; Performing morphological closing operations and edge pixel completion processing on the non-ground patches, and correcting misidentified patches and edge missing patches through intra-patch pixel consistency and elevation continuity detection to form an optimized non-ground patch mask; Converting the non-ground patch mask into a set of centroid points containing ground and non-ground semantics, marking the centroid points in the set of centroid points containing ground and non-ground semantics as ground centroid points and non-ground centroid points according to pixel values, identifying and completing missing points in the ground centroid points, and generating a centroid point layer for supporting ground centroid point extraction; A KD-tree index structure is constructed for the lidar point cloud and the centroid point layer respectively. In the three-dimensional coordinate space, the centroid point in the centroid point layer is used as a reference, and the nearest neighbor centroid point of each point in the input lidar point cloud is obtained using the KD-tree index. If the nearest neighbor centroid point belongs to the ground centroid point, the point in the input lidar point cloud corresponding to the nearest neighbor centroid point is marked as a ground point.

2. The method for automatic filtering of ground points for laser radar point clouds according to claim 1, characterized in that: The voxel size is obtained by calculating the average Euclidean distance between 1% of randomly selected point pairs in the input lidar point cloud.

3. The automatic filtering method for ground points of laser radar point cloud according to claim 1, characterized in that: When constructing the base centroid point layer, the elevation value of each centroid point is adjusted in an iterative manner; the steps for adjusting the elevation value of each centroid point are: selecting neighboring centroid points with lower elevations than the centroid point within the columnar neighborhood of the centroid point, calculating the median value of the neighboring centroid points, and using the median value to update the elevation of the current centroid point until the nearest neighbor spatial index is stable.

4. The method for automatic filtering of ground points for laser radar point clouds according to claim 1, characterized in that: The digital elevation model is a two-dimensional elevation grid constructed based on the near-geocentric point set and in accordance with spatial resolution.

5. The automatic filtering method for ground points of laser radar point cloud according to claim 1, characterized in that: Performing maximum filtering and threshold binarization processing on the digital elevation model to generate an initial non-ground cover area binarization map, including: The digital elevation model is subjected to maximum filtering and threshold binarization processing to generate a binary map marking elevation mutation areas; the processing formula of the binary map is: Among them, R HC For the binary image, NGNEM max is the maximum filtered digital elevation model, NGNEM is the digital elevation model, and W is the voxel size; Grid filling and elevation judgment are performed based on the spots in the binary map to generate an initial non-ground coverage area binary map.

6. The method for automatic filtering of ground points for laser radar point clouds according to claim 5, characterized in that: Based on the patch elevation features in the initial non-ground coverage area binary map, extracting the maximum elevation of the patch edge and the average elevation of the filling area, and extracting the non-ground patch according to the comparison result of the maximum elevation of the patch edge and the average elevation of the filling area, including: Perform unique value marking on each patch in the binary image; Perform grid filling operation on the marked spots to obtain the complete internal filling area; Calculate the maximum elevation of the patch edge and the average elevation of the filling area; If the maximum elevation of the patch edge is greater than the average elevation of the filling area, the patch is determined to be a non-ground patch in the initial non-ground coverage area binarization map.

7. The method for automatic filtering of ground points for laser radar point clouds according to claim 1, characterized in that: The steps to correct misidentified spots include: If a patch is marked as a non-ground patch before optimization, but the corresponding areas of all pixels in the elevation model after marking are not completely included in the non-ground patch mask, and the elevation of the patch has a continuous change trend in space with its adjacent ground patches, then the patch will be re-marked as a ground patch.

8. The method for automatic filtering of ground points for laser radar point clouds according to claim 1, characterized in that: When generating the centroid point layer, for the missing ground centroid points, the average elevation of the 24 nearest neighbor marked ground centroid points is used as the complementary elevation value to construct a continuous centroid point layer.

9. The method for automatic filtering of ground points for laser radar point clouds according to claim 1, characterized in that: The KD-tree index structure is constructed based on the lidar point cloud and the centroid point layer respectively, and is used to obtain the nearest neighbor centroid point of each input point in the three-dimensional coordinate space; if the nearest neighbor centroid point belongs to the ground centroid point, the corresponding input point is marked as a ground point.

10. The method for automatic filtering of ground points for laser radar point clouds according to claim 1, characterized in that: Before constructing the KD-tree index structure of the centroid layer, the elevation value of each centroid in the centroid layer is uniformly reduced by 2 times the voxel size.

Images