A ground point automatic filtering method for lidar point cloud

Abstract

The application provides a ground point automatic filtering method for a laser radar point cloud, relates to the technical field of point cloud data processing and terrain digital modeling, and comprises the following steps: carrying out voxelization processing on the point cloud, constructing a base centroid point layer, generating a near-ground digital elevation model, extracting and optimizing a non-ground map patch, constructing a semantic centroid point layer, and establishing a KD-tree index to extract a ground point. Through the combination of maximum value filtering and map patch elevation feature judgment, accurate identification and error correction of the non-ground map patch are realized, a centroid point layer is further constructed, a bidirectional spatial index structure is used to mark the ground point, 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 structure analysis, and has wide engineering application value.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of point cloud data processing and digital terrain modeling, in particular to a ground point automatic filtering method for laser radar point cloud. BACKGROUND

[0002] With the development of laser scanning technology, multi-temporal laser pulse signals reflected from the ground or target object surface are accurately captured and decoded into a large number of discrete point clouds, which can efficiently depict the three-dimensional structural features and surface morphological features of the scanned target. Point clouds have unique advantages in describing ground information, especially in complex terrain and multi-layer vegetation covered areas. Point clouds can not only describe the micro-terrain features under the forest in detail, but also accurately reflect the vertical structure features of the forest canopy, which is of great significance to high-precision terrain mapping, terrain feature extraction and forest canopy biomass assessment in vegetation covered areas. The spatial distribution of ground point clouds clearly depicts the terrain information such as slope and aspect, which can provide a scientific basis for agricultural machinery operation path planning; it can also be used to find suitable building areas, evaluate land development potential, especially in mountainous and hilly terrain complex areas; and provide support for orchard planting site selection. However, due to the mixed ground points and non-ground points (such as buildings, trees, power lines, etc.) in point clouds under complex scenes, it is a technical problem to distinguish the two types of points. Accurate identification and extraction of ground points are crucial for high-precision terrain modeling, and extraction of non-ground point structural features is also important for three-dimensional reconstruction of objects.

[0003] Ground filters (GFs) can be used to distinguish ground points and non-ground points from point clouds. GFs usually distinguish the two types of point clouds based on the fact that commonly used bandwidth laser pulses cannot penetrate soil layers or opaque building surfaces. Therefore, the ground point group extracted by GFs is usually located in the lower region of the local longitudinal section of the point cloud, covering the ground, underlying vegetation or building surface. Many different GF algorithms have been proposed, which can be divided into two categories according to their working principles: rule-based and machine learning-based algorithms, each of which has its own advantages and disadvantages. While meeting the accuracy requirements of GF results, users tend to use rule-based GFs because they do not require a large amount of training cost. The design of rule-based GFs is divided into slope threshold-based, morphological filter-based, spatial interpolation function-based and "cloth" filter-based schemes. However, each type of existing rule-based GF has its own difficult-to-apply scenarios, and the robustness of the algorithm needs to be improved, and automatic operation cannot be realized.

[0004] Due to the limited penetration ability of laser to high-density vegetation layer, the ground points of point cloud in dense forest and building-dense area are prone to be missed. During the ALS data acquisition process, the small scanning angle and sparse flight path setting are prone to cause incomplete description of the ground near the side of the building and the eaves, which increases the GF error. 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 flight height. Therefore, the GF must have sufficient robustness and accuracy to cope with the point cloud data obtained by scanning in different scenes. When optimizing the rule-based GF, reducing the number of parameters and enhancing its robustness in different scenes are the keys. It is a challenging task to design a GF that is accurate, easy to use and robust in complex scenes. SUMMARY

[0005] In order to overcome the shortcomings of the prior art, the purpose of the present application is to provide a ground point automatic filtering method for lidar point cloud, which can accurately and robustly automatically extract ground points from diversified scanning scenes. At the same time, the intermediate results in the analysis process of the GF can also assist in realizing instance segmentation of building point cloud to support building height measurement and morphological feature extraction.

[0006] To achieve the above purpose, the present application provides the following scheme:

[0007] A ground point automatic filtering method for lidar point cloud, comprising:

[0008] Voxelizing the input lidar point cloud, determining the voxel size according to the sampling statistical result, and after removing the empty voxels in the lidar point cloud according to the voxel size, extracting the centroid point set of non-empty voxels;

[0009] Based on the centroid point set of non-empty voxels, using the method of iteratively adjusting the elevation of the centroid point to construct the base centroid point layer, and taking the centroid point set of non-empty voxels as a reference, obtaining the near-ground centroid point set composed of nearest neighbors in the centroid point set of non-empty voxels, and constructing the digital elevation model of near-ground centroid points according to the near-ground centroid 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 based on the elevation features of the map patches 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 patch according to the comparison result of the maximum elevation of the patch edge and the average elevation of the filling area;

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

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

[0013] Respectively constructing KD-tree index structures for the laser radar point cloud and the centroid point layer, taking the centroid points in the centroid point layer as references in the three-dimensional coordinate space, using the KD-tree index to obtain the nearest neighbor centroid point of each point in the input laser radar point cloud, and if the nearest neighbor centroid point belongs to a ground centroid point, marking the point in the input laser radar point cloud corresponding to the nearest neighbor centroid point as a ground point.

[0014] Preferably, the voxel size is obtained by calculating the average Euclidean distance between randomly selected 1% points in the input laser radar 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 step of adjusting the elevation value of each centroid point is: selecting adjacent centroid points with lower elevations than the centroid point within the columnar neighborhood range of the centroid point, calculating the median value of the adjacent centroid points, and updating the elevation of the current centroid point with the median value until the nearest neighbor spatial index is stable.

[0016] Preferably, the digital elevation model is a two-dimensional elevation grid constructed according to the spatial resolution based on the near-ground centroid point set.

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

[0018] Performing maximum value filtering and threshold binarization processing on the digital elevation model to generate a binary map marking the elevation mutation area; the processing formula of the binary map is:

[0019]

[0020] wherein R HC is the binary map, NGNEM max is the maximum value filtered digital elevation model, NGNEM is the digital elevation model, and W is the voxel size.

[0021] Performing grid filling and elevation judgment based on the image patches in the binary image to generate an initial non-ground covered area binary image.

[0022] Preferably, based on the elevation characteristics of the image patches in the initial non-ground covered area binary image, the maximum elevation of the image patch edge and the average elevation of the filling area are extracted, and the non-ground image patch is extracted according to the comparison result of the maximum elevation of the image patch edge and the average elevation of the filling area, including:

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

[0024] Performing a grid filling operation on the marked image patch to obtain a complete internal filling area;

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

[0026] If the maximum elevation of the image patch edge is greater than the average elevation of the filling area, the image patch is determined as a non-ground image patch in the initial non-ground covered area binary image.

[0027] Preferably, the step of correcting the mis-recognized image patch includes:

[0028] If the image patch is marked as a non-ground image patch before optimization, but all the corresponding regions of the image patch in the elevation model after marking are not completely contained in the non-ground image patch mask, and the elevation of the image patch has a continuous change trend between its adjacent ground image patches in space, the image patch is re-marked as a ground image patch.

[0029] Preferably, 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 a completed elevation value to construct a continuous centroid point layer.

[0030] Preferably, the KD-tree index structure is constructed based on the laser radar 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 a 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 point layer, the elevation values of the centroid points in the centroid point layer are uniformly reduced by 2 times the size of the body element.

[0032] According to the specific embodiments provided by the present application, the following technical effects are disclosed:

[0033] This invention provides an automated ground point filtering method for lidar point clouds, comprising: performing voxelization on the input lidar point cloud; determining the voxel size based on sampling statistics; removing empty voxels from the lidar point cloud according to the voxel size; extracting a set of centroids of non-empty voxels; constructing a base centroid layer based on the set of centroids of non-empty voxels using an iterative adjustment method for centroid elevation; obtaining a set of near-geological centroids composed of nearest neighbors from the set of centroids of non-empty voxels using the set of centroids of non-empty voxels as a reference; constructing a digital elevation model of near-geological centroids based on the set of near-geological centroids; performing maximum value filtering and threshold binarization on the digital elevation model to generate an initial non-ground coverage area binarized map; extracting the maximum elevation of the patch edges and the average elevation of the filled area based on the patch elevation features in the initial non-ground coverage area binarized map; and extracting non-ground areas based on the comparison results of the maximum elevation of the patch edges and the average elevation of the filled area. The non-ground patches are subjected to morphological closing operations and edge pixel completion processing. Misidentified patches and patches with missing edges are corrected through pixel consistency and elevation continuity detection within the patches, forming an optimized non-ground patch mask. The non-ground patch mask is converted into a centroid set containing both ground and non-ground semantics. Centroids in this set are marked as ground centroids and non-ground centroids according to their pixel values. Missing points in the ground centroids are identified and filled in, generating a centroid layer to support ground centroid extraction. A KD-tree index structure is constructed for both the lidar point cloud and the centroid layer. In the three-dimensional coordinate space, using the centroids in the centroid layer as references, the nearest neighbor centroid for each point in the input lidar point cloud is obtained using the KD-tree index. If the nearest neighbor centroid is a ground centroid, the corresponding point in the input lidar point cloud is marked as a ground point. This invention can accurately and robustly extract ground points automatically from diverse scanning scenes. Furthermore, the intermediate results of this ground truth (GF) during the analysis process can assist in instance segmentation of building point clouds, supporting building height measurement and morphological feature extraction. Attached Figure Description

[0034] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the drawings used in the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0035] Figure 1 A flowchart of the method provided in an embodiment of the present invention;

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

[0037] Figure 3 This is a schematic diagram illustrating the construction of a digital elevation model of a point cloud base layer and near-geological core points provided in an embodiment of the present invention;

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

[0039] Figure 5 This is a schematic diagram of the ground point extraction process provided in an embodiment of the present invention. Detailed Implementation

[0040] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort 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, which can accurately and robustly extract ground points automatically from diverse scanning scenes. Furthermore, the intermediate results of this automated ground point filtering method during the analysis process can also assist in instance segmentation of building point clouds, supporting building height measurement and morphological feature extraction.

[0042] To make the above-mentioned objects, features and advantages of the present invention more apparent and understandable, the present invention will be further described in detail below with reference to the accompanying drawings and specific embodiments.

[0043] Figure 1 The method flowchart provided in the embodiments of the present invention is as follows: Figure 1 As shown, this invention provides an automated ground point filtering method for lidar point clouds, comprising:

[0044] The input lidar point cloud is processed into voxels. The voxel size is determined based on the sampling statistics. After removing empty voxels in the lidar point cloud according to the voxel size, the set of centroids of the non-empty voxels is extracted.

[0045] Based on the centroid set of non-empty volume elements, an iterative method for adjusting the centroid elevation is used to construct a base centroid layer. Using the centroid set of non-empty volume elements as a reference, a near-geological centroid set composed of the nearest neighbors is obtained from the centroid set of non-empty volume elements, and a digital elevation model of the near-geological centroid is constructed based on the near-geological centroid set.

[0046] The digital elevation model is subjected to maximum value filtering and threshold binarization to generate an initial non-ground coverage area binarized map. Based on the elevation features of the patches in the initial non-ground coverage area binarized map, the maximum elevation of the patch edge and the average elevation of the filled area are extracted. The non-ground patches are extracted based on the comparison results between the maximum elevation of the patch edge and the average elevation of the filled area.

[0047] Morphological closing operations and edge pixel completion are performed on non-ground patches. Misidentified patches and patches with missing edges are corrected by detecting pixel consistency and elevation continuity within the patches, thus forming an optimized non-ground patch mask.

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

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

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

[0051] (1) Point cloud volume elementation and near-geocentric point set extraction

[0052] LWF first voxels the input point cloud, which effectively reduces the computational cost of point-by-point analysis. The voxel size is represented by W. W is determined by randomly selecting a subset of points representing 1% of the total number of points in the input point cloud. Then, the average Euclidean distance between each selected point and its nearest neighbor is calculated. To reduce the number of empty volume elements while preserving as much structural detail as possible in the point cloud volumetric representation, W was determined through parametric ablation experiments.

[0053]

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

[0055] Because multiple echo laser pulses can penetrate plant canopies but not the ground, they can be obtained from the volumetric centroid set (ND) of the input point cloud. CE The centroid subsets of near-ground volume elements are extracted 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 subsets in building and vegetation-covered areas are always missing in the ND... CE In addition, some ND coverings on roof and canopy surfaces CE Easily mistaken for near-surface centroids. Some NDs with anomalous heights CE This is due to the presence of noisy points in the input point cloud. Removing point cloud noise and misclassified near-ground centroids is crucial for the accurate operation of this invention.

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

[0057] Extracting ND from the input point cloud CE Then, gradually reduce ND. CE The elevation of the input point cloud is used to construct the basal centroid layer, which supports the extraction of ground points. First, for ND... CE Construct a two-dimensional KD-tree spatial index with a search radius set to To obtain the local centroid subset C of cylindrical distribution along the z-axis of centroid G. p Comparison with C p Compare the elevations of other centroids with the elevation of G. If C p If there are centroids whose elevation is lower than G, then the median elevation of these centroids is calculated and assigned to the elevation of G. ND is then adjusted point by point using this method. CE The elevation is used to form the base layer of the input point cloud. If in C p In the analysis, the spatial index of the nearest neighbor centroid of the target centroid G differs between the N+1th and Nth iterations. Therefore, according to C... p The elevation of G is adjusted using the median elevation of the lower subset of points until the two indices are the same. Finally, the basal centroid layer ND of the input point cloud is obtained. BL Besides elevation, ND BL The (x,y) coordinates of each point in the ND CE same.

[0058] For ND BL The target centroid Q in the ND is obtained by using a three-dimensional KD-tree search engine in ND. CE Find its nearest centroid NQ To obtain ND CE Near-Geocentric Point Set (ND) GD ), ND GD It is ND CE A subset of ND. GD Some centroids from the rooftops and forest canopy surfaces are still retained. To remove these error points, based on ND... CE The elevation information was used to construct a near-geological center digital elevation model (NGNEM) at Wm resolution. The elevation changes described by the NGNEM were used to drive the detection of ND. GD Non-ground centroids in the data. This is based on ND. CE To construct NGNEM is to prevent ND GD Retaining only partially misaligned centroids covering rooftops and tree canopies reduces the completeness of identifying non-ground-covered areas from NGNEM.

[0059] (3) Ground coverage area identification

[0060] In medium-to-high density lidar point clouds, the height difference between ground points and rooftop points typically exceeds W. Local point cloud elevation variations are continuous near cliffs, ridges, and steep slopes. This continuous topographic variation provides a basis for delineating ground and non-ground regions (e.g., rooftops and trees) in the NGNEM. This invention applies a maximum value filter with an analysis window size of 5×5 to analyze the NGNEM, obtaining the NGNEM... max Therefore, using equation (2) to analyze NGNEM max Binarization is performed to obtain R. HC .

[0061]

[0062] The edges of non-ground objects in R HC The outline is drawn in the middle. Therefore, the Union-Find algorithm is used to refine R. HC The polygon in the image is uniquely labeled with a value, namely R. HC Each individual patch shown is assigned a unique number. Then, for a tagged patch P... L Perform raster fill analysis. If P L If it can be filled, then calculate the maximum elevation (H) of the edge of the patch. edge ) and the average elevation (H) of the internal filling area fill If H fill ≥H edge Then it can be determined that P L This area is not covered by ground objects. Otherwise, P L Located in a steep slope or plateau area with continuous elevation changes.

[0063] Because the elevation information described by NGNEM may contain some misclassifications from rooftops or tree canopies, ND GD Centroids are difficult to detect in NGNEM patches using only hole-filling algorithms. To address this issue, morphological closing operations are applied to analyze labeled patches in NGNEM to obtain NGNEM data. BI The window size w used for morphological closing analysis. close for:

[0064]

[0065] Where, N close This represents the number of pixels filled within each infillable patch. Closed patches of non-ground objects in NGNEM... BI Cells in the middle are marked with a value of "1", while the rest are marked with a value of "0". This is to avoid some non-terrestrial patches near the NGNEM edge being included in the NGNEM. BI NGNEM lost in China 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, thus obtaining the NGNEM. edge Then, Equation (4) is used to calculate ED. P .

[0066]

[0067] If ED P If a pixel in a patch has two non-zero values, it is considered a non-terrestrial patch located at the edge of the NGNEM.

[0068] These non-terrestrial patches located at the edge of NGNEM are merged into NGNEM. BI middle.

[0069] In NGNEM BI In NGNEM, patches located on slopes may be misidentified as non-ground types. To address this issue, this invention assigns unique values ​​to patches in the NGNEM based on elevation changes. In the annotation results, areas at the same elevation in the NGNEM will form multiple patches. If, in the labeled NGNEM, all pixels within the coverage area of ​​a target patch are within the NGNEM... BI If all pixels in the NGNEM have a value of "1", then the pixel is classified as a non-terrestrial patch in the NGNEM. Otherwise, the patch is reassigned to "0" in the labeled NGNEM, and is considered a terrestrial patch. Binarizing the NGNEM with uniquely labeled patches yields the NGNEM. NG .

[0070] Next, morphological dilatation analysis was applied to NGNEM. NG The image features in NGNEM. NG Morphological dilatation results and NGNEM NG Subtraction yields the outer boundary cell set for each patch. This boundary pixel set is used to extract elevation changes near the edge of a single patch from the NGNEM. For a patch, if the maximum elevation of its edge cell coverage area exceeds the average elevation within its coverage area, it is classified as a ground patch on a slope. After identifying ground patches on slopes, this invention completes the NGNEM... NG Optimization.

[0071] (4) Ground point extraction preparation

[0072] According to NGNEM NG The row and column indices, elevation information, and spatial resolution W of each pixel are used to transform the center point of each pixel back into a spatially discrete centroid. In the transformation result, the type of the centroid is determined according to the NGNEM. NG The pixel value is determined as follows: the centroid TP with a value of "1". N Centroid TP, located in a non-terrestrial patch, has a value of "0". G Located within a ground patch. Each TP N and TP G The elevation of the point is determined by ND. CE Centroids with the same (x,y) index are identified. To improve the accuracy of GF, ND is also utilized. CE (x,y) index recognition of the centroid point ND GD The missing (x,y) index is used to determine the centroids that need to be added. For each new centroid, its elevation is assigned to its 24 nearest neighbors TP. G The average elevation of the points was used as the basis for constructing a centroid layer located below the input point cloud coverage area.

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

[0074] A KD-tree spatial index (OKDT) was constructed for the input point cloud, and another KD-tree spatial index (TKDT) was constructed for the near-ground centroid layer to form a bidirectional KD-tree spatial index. Prior to this step, the elevation of the near-ground centroid set was reduced by 2W to avoid extracting point clouds from the surfaces of small non-ground objects (e.g., vehicles parked on the road). For each point in the OKDT index, its nearest centroid was searched. Then, the number of points in the input point cloud sharing the same nearest neighbor centroid was counted, and their average (M, M≥1) was calculated. For each TP in the input point cloud... GThe node searches for M neighboring points. In the input point cloud, if a point's nearest centroid belongs to TP... N Then it is located in an area not covered by ground. For areas near TP G If a node has M points, and its nearest centroid belongs to TP, then... G Record its index in OKDT to obtain a list G index According to G index It can extract ground points from the input point cloud.

[0075] As an optional implementation, this embodiment utilizes a UAV or airborne LiDAR scanner to scan and acquire discrete point clouds within the flight path coverage area according to a 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 in this invention. The point cloud is voxelized, and empty voxels are removed. The size of each voxel is set to W. The centroids of each non-empty voxel are extracted to reduce the computational cost of subsequent analysis. After extracting the set of centroid points of the voxels, the elevation of the centroid points is gradually reduced to construct the base layer of the input point cloud. (ND) CE Construct a two-dimensional KD-tree spatial index with a search radius set to To obtain the cylindrical local centroid subset C along the z-axis of each centroid. p Comparison with C p The elevations of other centroids and the target centroid. If C p If the elevation of any centroid point is lower than that of the target centroid point, then the median elevation of these centroid points is calculated and assigned to the elevation of G. ND is then adjusted point by point using this method. CE The elevation is used to form the base layer of the input point cloud. If in C p In the analysis, the spatial index of the nearest neighbor centroid of the target centroid differs between the N+1th and Nth iterations. Therefore, according to C... p The elevation of G is adjusted using the median elevation of the lower subset of points until the two indices are the same. This ultimately yields the basal centroid layer ND of the input point cloud. BL .

[0076] For ND BL The target centroid Q in the model is obtained using a 3D KD-tree search engine in ND. CE Find its nearest centroid N Q To obtain ND CE Near-Geocentric Point Set (ND) GD To remove ND GD Some centroids from rooftops and forest canopy surfaces are still preserved, based on ND. CEThe elevation information was used to construct a near-geological center digital elevation model (NGNEM) at a resolution of Wm. The NGNEM was analyzed using a maximum value filter with a 5×5 window size, yielding the NGNEM. max For NGNEM max Binarization is performed to obtain R. HC The edges of non-ground objects are in R. HC It is outlined in the middle.

[0077] Drawing inspiration from the Union-Find algorithm for R HC The polygon in the image is uniquely labeled with a value, namely R. HC Each individual patch shown is assigned a unique number. Then, for a labeled patch P... L Perform raster fill analysis. Based on P L The maximum elevation (H) of the edge of the map patch edge ) and the average elevation (H) of the internal filling area fill ) size, determine P L Does it belong to a non-ground object or a ground-covered area? When H fill <H edge P L Located on a steep slope or plateau area with continuous elevation changes. Otherwise, P L Located in an area not covered by ground objects.

[0078] Morphological closing operations are applied to analyze labeled patches in NGNEM to obtain NGNEM data. BI Closed patches of non-ground objects in NGNEM BI Cells in the middle are marked with a value of "1", while the rest are marked with a value of "0". This is to avoid some non-terrestrial patches near the NGNEM edge being included in the NGNEM. BI NGNEM lost in China 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, thus obtaining the NGNEM. edge ED is calculated accordingly P If ED P If a pixel within a patch has fewer than two non-zero values, it is considered a non-terrestrial patch located at the edge of the NGNEM. These non-terrestrial patches located at the edge of the NGNEM are then merged into the NGNEM. BI middle.

[0079] To reduce NGNEM BIThe proportion of patches located on slopes being misidentified as non-ground types is addressed in this invention by uniquely labeling patches in the NGNEM based on elevation changes, thus creating multiple patches within the same elevation area. If, in the labeled NGNEM, all pixels within the coverage area of ​​a target patch are within the NGNEM... BI If all pixels in the NGNEM have a value of "1", then the pixel is classified as a non-terrestrial patch in the NGNEM. Otherwise, the patch is reassigned to "0" in the labeled NGNEM, i.e., it is considered a terrestrial patch. The labeled NGNEM is then binarized to obtain the NGNEM. NG Applying morphological dilatation analysis to NGNEM NG The image patches in the image, and compared with NGNEM NG Subtraction yields the set of outer boundary pixels for each patch. This set of boundary pixels for each patch is used to extract elevation variations near the edges of individual patches from the NGNEM. For a patch, if the maximum elevation of its edge pixel coverage area exceeds the average elevation within its coverage area, it is classified as a ground patch on a slope. After identifying ground patches on slopes, the NGNEM is completed. NG optimization.

[0080] According to NGNEM NG The row and column indices, elevation information, and spatial resolution W of each pixel are used to transform the center point of each pixel back into a spatially discrete centroid. To improve the accuracy of the GF (Geometric Grid), the ND (Digital Neutral Density) is also utilized. CE (x,y) index recognition of the centroid point ND GD The missing (x, y) index is used to determine the centroid that needs to be added. For each new node, its elevation is assigned to its 24 nearest neighbors TP. G The average elevation of the points was determined. Based on this, a centroid layer was constructed below the area covered by the input point cloud. A KD-tree spatial index (OKDT) was constructed for the input point cloud, and another KD-tree spatial index (TKDT) was constructed for the centroid layer to form a bidirectional KD-tree spatial index. Prior to this step, the elevation of the centroid points in the centroid layer was reduced by 2W to avoid extracting point clouds from the surface 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). For each TP in the input point cloud... G The node searches for M neighboring points. In the input point cloud, if a point's nearest centroid belongs to TP... N Then it is located in an area not covered by ground. For areas near TP GIf a node has M points, and its nearest centroid belongs to TP, then... G Record its index in OKDT to obtain a list G index According to G index Finally, ground points are extracted from the input point cloud.

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

[0083] This document uses specific examples to illustrate the principles and implementation methods of the present invention. The descriptions of the above embodiments are only for the purpose of helping to understand the method and core ideas of the present invention. Furthermore, those skilled in the art will recognize that, based on the ideas of the present invention, there will be changes in the specific implementation methods and application scope. Therefore, the content of this specification should not be construed as a limitation of the present invention.

Claims

1. An automated ground point filtering method for lidar point clouds, characterized in that, include: The input lidar point cloud is processed into voxels. The voxel size is determined based on the sampling statistics. After removing empty voxels in the lidar point cloud according to the voxel size, the set of centroids of the non-empty voxels is extracted. Based on the set of centroids of the non-empty volume element, an iterative method for adjusting the elevation of the centroids is used to construct a base centroid layer. With the set of centroids of the non-empty volume element as a reference, a set of near-geological centroids is obtained from the set of centroids of the non-empty volume element, which is the nearest neighbor of the target centroid. A digital elevation model of the near-geological centroids is constructed based on the set of near-geological centroids. The digital elevation model is subjected to maximum value filtering and threshold binarization to generate an initial non-ground coverage area binarized map. Based on the elevation features of the patches in the initial non-ground coverage area binarized map, the maximum elevation of the patch edge and the average elevation of the filled area are extracted. Non-ground patches are extracted based on the comparison results of the maximum elevation of the patch edge and the average elevation of the filled area. Morphological closing operations and edge pixel completion processing are performed on the non-ground patches. Misidentified patches and patches with missing edges are corrected by detecting pixel consistency and elevation continuity within the patches, thus forming an optimized non-ground patch mask. The non-ground patch mask is converted into a centroid set containing ground and non-ground semantics. The centroids in the centroid set containing ground and non-ground semantics are marked as ground centroids and non-ground centroids according to their pixel values. The missing points in the ground centroids are identified and filled in to generate a centroid layer to support ground centroid extraction. A KD-tree index structure is constructed for the lidar point cloud and the centroid layer respectively. In the three-dimensional coordinate space, with the centroid in the centroid layer as a reference, the nearest neighbor centroid of each point in the input lidar point cloud is obtained using the KD-tree index. If the nearest neighbor centroid belongs to the ground centroid, the point in the input lidar point cloud corresponding to the nearest neighbor centroid is marked as a ground point. When constructing the base centroid layer, the elevation value of each centroid is adjusted iteratively. The steps for adjusting the elevation value of each centroid are as follows: select neighboring centroids with lower elevations within the columnar neighborhood of the centroid, calculate the median value of the neighboring centroids, and update the elevation of the current centroid with the median value until the nearest neighbor spatial index is stable.

2. The automated ground point filtering method for lidar point clouds according to claim 1, characterized in that, The volume element size is calculated by taking the average Euclidean distance between 1% of the randomly selected point pairs in the input lidar point cloud.

3. The automated ground point filtering method for lidar point clouds according to claim 1, characterized in that, The digital elevation model is a two-dimensional elevation grid constructed based on the set of near-geological centers according to spatial resolution.

4. The automated ground point filtering method for lidar point clouds according to claim 1, characterized in that, The digital elevation model is subjected to maximum value filtering and threshold binarization to generate an initial binarized map of the non-ground coverage area, including: The digital elevation model is subjected to maximum value filtering and threshold binarization to generate a binary map marking regions of abrupt elevation changes; the processing formula for the binary map is as follows: in, For the binary image, This is the digital elevation model after maximum value filtering. For digital elevation models, For volume element dimensions; Based on the patches in the binary map, raster filling and elevation determination are performed to generate an initial binary map of the non-ground coverage area.

5. The automated ground point filtering method for lidar point clouds according to claim 4, characterized in that, Based on the elevation features of the patches in the initial binary map of the non-ground covered area, the maximum elevation of the patch edges and the average elevation of the filled area are extracted. Then, based on the comparison between the maximum elevation of the patch edges and the average elevation of the filled area, non-ground patches are extracted, including: Perform a unique value labeling on each patch in the binary graph; Perform a raster fill operation on the marked patches to obtain the complete internal fill area; Calculate the maximum elevation of the edge of the map patch and the average elevation of the filled area; If the maximum elevation of the edge of the patch is greater than the average elevation of the filled area, then the patch is determined to be a non-ground patch in the binarized map of the initial non-ground covered area.

6. The automated ground point filtering method for lidar point clouds according to claim 1, characterized in that, The steps to correct misidentified patches include: If a patch is marked as a non-ground patch before optimization, but the area corresponding to all pixels in the marked elevation model is not completely included in the non-ground patch mask, and the elevation of the patch has a continuous variation trend with its adjacent ground patches in space, then the patch is remarked as a ground patch.

7. The automated ground point filtering method for lidar point clouds according to claim 1, characterized in that, When generating the centroid layer, for missing ground centroids, the average elevation of the 24 nearest neighbor labeled ground centroids is used as the fill-in elevation value to construct a continuous centroid layer.

8. The automated ground point filtering method for lidar 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 layer, respectively, and is used to obtain the nearest centroid of each input point in the three-dimensional coordinate space; if the nearest centroid belongs to the ground centroid, the corresponding input point is marked as a ground point.

9. The automated ground point filtering method for lidar point clouds according to claim 1, characterized in that, Before constructing the KD-tree index structure of the centroid layer, the elevation values ​​of each centroid in the centroid layer are uniformly reduced by twice the volume element size.

Images