CN109785261B - Airborne LIDAR three-dimensional filtering method based on gray voxel model - Google Patents
Airborne LIDAR three-dimensional filtering method based on gray voxel model Download PDFInfo
- Publication number
- CN109785261B CN109785261B CN201910027077.0A CN201910027077A CN109785261B CN 109785261 B CN109785261 B CN 109785261B CN 201910027077 A CN201910027077 A CN 201910027077A CN 109785261 B CN109785261 B CN 109785261B
- Authority
- CN
- China
- Prior art keywords
- voxel
- gray
- ground
- elevation
- airborne lidar
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
Images
Landscapes
- Image Processing (AREA)
Abstract
The invention provides an airborne LIDAR three-dimensional filtering method based on a gray level voxel model, and relates to the technical field of remote sensing data processing. The method comprises the steps of firstly, reading original airborne LIDAR point cloud data to form an original airborne LIDAR point cloud data set; regularizing airborne LIDAR point cloud data into a gray voxel model, wherein the gray is discretization representation of the average intensity of laser points in voxels; and selecting a ground seed voxel based on the lowest local elevation characteristic of the ground point, and marking a non-0 value voxel which is three-dimensionally communicated with the ground seed voxel and has a voxel value close to the terrain slope as a ground voxel. According to the airborne LIDAR three-dimensional filtering method based on the gray voxel model, the theory of three-dimensional connected region construction is taken as a basis, the intensity and terrain gradient information is used as auxiliary information for point cloud data filtering, more effective information can be provided for distinguishing ground targets from non-ground targets, the filtering precision is improved, and the three-dimensional filtering method based on the gray voxel model is expanded and is suitable for more complex scenes.
Description
Technical Field
The invention relates to the technical field of remote sensing data processing, in particular to an airborne LIDAR three-dimensional filtering method based on a gray voxel model.
Background
The data filtering of an airborne laser radar (LIDAR) is a key technology in the field of current point cloud data processing, and the filtering precision directly affects the quality of products such as a Digital Elevation Model (DEM) and the like and the subsequent precision of object classification such as buildings and vegetation. Therefore, the method has important significance for the research of the LIDAR point cloud data filtering technology.
Scholars at home and abroad carry out deep research around point cloud data filtering, and a great number of filtering algorithms are provided, such as filtering methods based on interpolation, gradient, mathematical morphology, three-dimensional connected set construction and cluster segmentation, and the like. Data structures used by these filtering algorithms include grid meshes, irregular triangulated networks (TINs), voxels, and discrete point clouds, among others. Grid meshes and TINs are 2.5-dimensional data structures with which to express airborne LIDAR point cloud data will result in information loss and further impact the integrity of filtering results based on such data structures. The spatial structure and topological information in the discrete point cloud are difficult to utilize, so that the filtering algorithm based on the data structure is difficult to design and low in efficiency; the sizes of all nodes of the octree are different, and the adjacency relation among all nodes is difficult to establish, so that the difficulty of the design of the filtering algorithm based on the data structure is increased. The filtering algorithm based on the voxel data structure can well avoid the defects of the method, because: (1) the voxel structure is a true three-dimensional data structure; (2) The inter-voxel inside the system contains adjacent and topological information. Therefore, the filtering algorithm based on the voxel structure is simpler to design. However, currently existing filtering algorithms based on voxel data structures all use binary voxel structures (i.e. the voxel assignment is based on whether the voxel contains laser point assignment 1 and 0, wherein 1 value voxel corresponds to the target, and 0 value voxel corresponds to the background, respectively) (coloring Wang, yan Xu, and Yu Li. The filtering algorithm considers that the ground target forms a three-dimensional connected set, but when other targets (such as low vegetation and building facades) are connected with the ground target to form the three-dimensional connected set, the algorithm cannot effectively distinguish the two targets. In addition, the algorithm does not take into account terrain slope information, thereby resulting in a lower filtering accuracy of the algorithm.
Disclosure of Invention
The technical problem to be solved by the invention is to provide an airborne LIDAR three-dimensional filtering method based on a gray voxel model aiming at the defects of the prior art, and solve the problem of effectively distinguishing connected ground and non-ground targets.
In order to solve the technical problems, the technical scheme adopted by the invention is as follows: an airborne LIDAR three-dimensional filtering method based on a gray voxel model comprises the following steps:
step 1: reading original airborne LIDAR point cloud data to form an original airborne LIDAR point cloud data set;
step 2: regularizing an original airborne LIDAR point cloud dataset into a gray voxel model, wherein the gray is discretized representation of the average intensity of laser points in voxels, and the method specifically comprises the following steps:
step 2.1: removing abnormal data from an original airborne LIDAR point cloud data set to obtain a removed abnormal data set;
step 2.1.1: counting the frequency of each laser point elevation value in the original airborne LIDAR point cloud data set, and displaying the counting result in a histogram mode in a visualized manner;
step 2.1.2: determining a maximum elevation threshold T corresponding to real terrain and ground features he And a minimum elevation threshold T le ；
Step 2.1.3: aiming at each laser point in the original airborne LIDAR point cloud data set, if the elevation value of each laser point is higher than the highest elevation threshold value T he Or below a minimum elevation threshold T le If the laser point is the elevation abnormal data, removing the laser point, otherwise, keeping the laser point, and obtaining a removed elevation abnormal data set;
step 2.1.4: counting and eliminating the frequency of the intensity value of each laser point in the elevation abnormal data set, and displaying the counting result in a histogram mode in a visualized manner;
step 2.1.5: determining a maximum intensity threshold T corresponding to real terrain and terrain hi And a minimum intensity threshold T li ；
Step 2.1.6: aiming at each laser point in the removed elevation abnormal data set, if the intensity value is higher than the highest intensity threshold value T hi Or below a minimum intensity threshold T hi If the laser point is abnormal intensity data, rejecting the laser point, otherwise, keeping the laser point, and finally obtaining a rejected elevation and abnormal intensity data set;
step 2.2: regularizing the removed abnormal data set into a gray voxel model;
step 2.2.1: representing the spatial extent of the dataset by an axial parallel bounding box which rejects the abnormal dataset;
step 2.2.2: calculating the resolution of the volume element, namely the size of the volume element, wherein the resolution in the directions of x, y and z is determined according to the average point distance of the laser points in the abnormal data set removed;
step 2.2.3: dividing the axial parallel bounding boxes from which the abnormal data sets are removed according to the resolution ratios in the x direction, the y direction and the z direction to obtain three-dimensional grids, and calling each three-dimensional grid unit as a volume element;
step 2.2.4: mapping each laser point in the abnormal data set to a three-dimensional grid, assigning values to each voxel according to the intensity attribute of the laser points contained in the voxel, and discretizing the values of the voxels to {0, \ 8230;, 255} to obtain a gray voxel model;
and step 3: based on the theory of three-dimensional connected region construction, the ground voxel in the gray voxel model is detected, and the specific method comprises the following steps:
step 3.1: based on the local elevation minimum characteristic of the ground point, taking the voxel with the lowest non-0 value in the height from the gray voxel model as a ground seed voxel set;
step 3.1.1: in the horizontal direction, the gray voxel model is partitioned by using the set grid size, and the voxel with the lowest non-0 value in the elevation in each block is taken as a seed voxel;
step 3.1.2: searching in each block for a height difference with the ground seed smaller than a height difference threshold value T e The non-0 value voxel is a ground seed voxel, and an encrypted ground seed voxel set V is obtained s Wherein, s =1, 2.;
step 3.2: marking a ground seed voxel and a three-dimensional connected region formed by voxels which are three-dimensionally connected with the ground seed voxel and have similar gradient and gray level as a ground body metadata set, and finishing airborne LIDAR three-dimensional filtering based on a gray level voxel model;
step 3.2.1: calculating a frequency histogram of gray values of voxels with non-0 values in the gray voxel model, fitting the gray frequency histogram by using a Gaussian mixture model, and determining a gray value distribution range corresponding to a ground target;
step 3.2.2: for any ground seed voxel V s Traversing the voxel V of the gray voxel model and the current ground seed voxel s Three dimensional connectivityThe gray value is positioned in the gray range corresponding to the ground target, and the local terrain gradient is smaller than the gradient threshold value T s All unlabeled voxels of (2), and labeling as L g Until all the ground seed voxels V are marked s A three-dimensional connected region of (a), a ground voxel set; the gradient threshold value is determined in a self-adaptive mode, and the specific method comprises the following steps: for the current seed voxel V s Detecting whether the space neighborhood contains marked ground voxels, if so, determining the maximum gradient value among the marked ground voxels as a terrain gradient threshold value T in the neighborhood s (ii) a Otherwise, the terrain grade threshold is set to 90 °.
Adopt the produced beneficial effect of above-mentioned technical scheme to lie in: according to the airborne LIDAR three-dimensional filtering method based on the gray voxel model, the theory of three-dimensional connected region construction is taken as a basis, the strength and terrain gradient information is used as auxiliary information for point cloud data filtering, more effective information can be provided for distinguishing ground targets from non-ground targets, the filtering precision is improved, and the three-dimensional filtering method based on the gray voxel model is expanded and is suitable for more complex scenes.
Drawings
Fig. 1 is a flowchart of an airborne LIDAR three-dimensional filtering method based on a grayscale voxel model according to an embodiment of the present invention;
fig. 2 is a schematic diagram of original airborne LIDAR point cloud data provided by an embodiment of the present invention, where (a) is Csite1 point cloud data, (b) is Csite2 point cloud data, (c) is Csite3 point cloud data, (d) is Csite4 point cloud data, (e) is Fsite5 point cloud data, (f) is Fsite6 point cloud data, and (g) is Fsite7 point cloud data;
FIG. 3 is a detailed flow chart for regularizing an original airborne LIDAR point cloud dataset into a grayscale voxel model according to an embodiment of the present invention;
FIG. 4 is a schematic diagram of calculating a convex hull and an area by laser point projection according to an embodiment of the present invention;
FIG. 5 is a flowchart of ground voxel detection on grayscale voxel dataset according to an embodiment of the present invention;
FIG. 6 is a grayscale frequency histogram of voxels with non-0 values in a three-dimensional voxel matrix V corresponding to a sample samp41 according to an embodiment of the present invention;
FIG. 7 is a diagram illustrating an arbitrary seed voxel V in a depth-first traversal gray voxel model according to an embodiment of the present invention s Is shown in (a).
Detailed Description
The following detailed description of embodiments of the present invention is provided in connection with the accompanying drawings and examples. The following examples are intended to illustrate the invention but are not intended to limit the scope of the invention.
In this embodiment, a certain group of remote sensing data is taken as an example, and the airborne LIDAR three-dimensional filtering method based on the gray voxel model is used to perform filtering analysis on the group of data.
An airborne LIDAR three-dimensional filtering method based on a gray voxel model, as shown in fig. 1, includes the following steps:
step 1: and reading the original airborne LIDAR point cloud data to form an original airborne LIDAR point cloud data set.
In this embodiment, 15 sample data provided by the third working group of the International Society for Photogrammetry and Remote Sensing (isps) and specially used for a filter algorithm test are used as experimental data to check the validity and feasibility of the method. As shown in fig. 2, the sample data is cut from 7 test data, which includes scenes with different complexity and situations where filtering may encounter difficulties, such as the influence of rough difference points, complex features, features connected to the ground, vegetation on slopes or low, discontinuities in terrain, etc., and table 1 lists the characteristics and basic parameters of each sample data in this embodiment. The sample records the first and last echoes and intensity information of the point cloud (the sample data does not contain the intensity information of the point cloud data, and the intensity information is extracted from the corresponding 7 test data).
TABLE 1 characteristics and basic parameters of sample data
In the embodiment, reference data (15 sample data are accurately classified into ground points and non-ground points) provided by the ISPRS is adopted as standard data to evaluate the algorithm precision.
In this embodiment, defining original airborne LIDAR point cloud data P = { P = { (P) } i (x i ，y i ，z i ) I =1, \ 8230;, n }, where i is an index of the original airborne LIDAR point cloud data, n is the number of the original airborne LIDAR point cloud data, p i Is the ith original airborne LIDAR point cloud data with coordinates of (x) i ，y i ，z i )。
And 2, step: regularizing an original airborne LIDAR point cloud data set into a gray voxel model, as shown in fig. 3, the specific method is as follows:
step 2.1: and rejecting abnormal data from the original airborne LIDAR point cloud data set to obtain a rejected abnormal data set.
Step 2.1.1: counting the frequency of each laser point elevation value in the original airborne LIDAR point cloud data set, and displaying the counting result in a histogram mode in a visualized manner;
step 2.1.2: determining the highest elevation threshold T corresponding to the real terrain and ground features he And a minimum elevation threshold T le ；
Step 2.1.3: aiming at each laser point in the original airborne LIDAR point cloud data set, if the elevation value is higher than the highest elevation threshold value T he Or below a minimum elevation threshold T le If the laser point is the elevation abnormal data, removing, otherwise, keeping the laser point to obtain a removed elevation abnormal data set;
step 2.1.4: counting and eliminating the frequency of the intensity value of each laser point in the elevation abnormal data set, and displaying the counting result in a histogram mode in a visualized manner;
step 2.1.5: determining a maximum intensity threshold T corresponding to real terrain and terrain hi And a minimum intensity threshold T li ；
Step 2.1.6: aiming at each laser point in the removed elevation abnormal data set, if the intensity value is higher than the highest intensity threshold value T hi Or below the lowestIntensity threshold T hi And if not, the laser point is reserved, and finally a data set with the abnormal elimination elevation and intensity is obtained.
In this embodiment, the maximum elevation threshold T he Minimum elevation threshold T le Maximum intensity threshold T hi And a minimum intensity threshold T li The values of all the points are constants, and the values are determined according to the spatial distribution condition of the original airborne LIDAR point cloud data.
In this embodiment, the removed abnormal data set is denoted as Q = { Q = i′ (x i′ ，y i′ ，z i′ ) I '=1, \ 8230;, t }, where i' is the index of the laser spot in the rejected abnormal data set, t is the number of laser spots in the rejected abnormal data set, q i′ The ith 'laser point in the abnormal data set is removed, and the coordinate of the ith' laser point is (x) i′ ，y i′ ，z i′ )。
Step 2.2: and regularizing the abnormal data eliminating set into a gray voxel model.
Step 2.2.1: the spatial extent of the dataset is represented by the axially parallel bounding box that culls the outlier dataset.
In this embodiment, the spatial range of the abnormal data set Q to be removed may be determined by an axial parallel Bounding Box (Axis-Aligned Bounding Box, AABB), AABB = { (x, y, z) | x min ≤x≤x max ，y min ≤y≤y max ，z min ≤z≤z max In the formula (x) max ，y max ，z max ) And (x) min ，y min ，z min ) Representing the maximum and minimum values of the x, y and z coordinates in the culled outlier dataset, respectively.
Step 2.2.2: and calculating the resolution of the volume element, namely the size of the volume element, wherein the resolution in the directions of x, y and z is determined according to the average point distance of the laser points in the abnormal data set.
In the present embodiment, the formula for calculating the resolutions Δ x, Δ y, and Δ z of voxels in the x, y, and z directions is as shown in formula (1):
wherein S is xy ＝{(x i′ ，y i′ )，i′＝1，…，t}Sxz＝{(x i′ ，z i′ )，i′＝1，…，t}、S yz ＝{(y i′ ，z i′ ) I' =1, \ 8230;, t } are two-dimensional point sets obtained by eliminating the projection of the abnormal data set Q on the XOY, XOZ and YOZ planes respectively, C (·) is the convex hull of the corresponding point set, and a (·) is the area of the corresponding convex hull, as shown in fig. 4.
Step 2.2.3: and dividing the axial parallel bounding boxes according to the resolution ratios in the x direction, the y direction and the z direction to obtain three-dimensional grids, wherein each three-dimensional grid unit is called a voxel.
In this embodiment, the axially parallel bounding box may be divided into three-dimensional grids based on the voxel resolutions (Δ x, Δ y, Δ z), which may be represented by a three-dimensional voxel matrix. Let V be the voxel set in the three-dimensional voxel matrix, as shown in equation (2):
V＝{v j (r j ，c j ，l j )，j＝1，…，m} (2)
where j is the voxel index; m is the voxel number; v. of j Is the voxel value of the jth voxel; (r) j ，c j ，l j ) Is the coordinate of the jth voxel in the voxel array, r j ，c j ，l j Respectively the row, column and layer number of the coordinates.
Step 2.2.4: mapping each laser point in the abnormal data set to a three-dimensional grid, assigning values to each voxel according to the intensity attribute of the laser points contained in the voxel, and discretizing the values of each voxel to {0, \ 8230;, 255}, thereby obtaining a gray voxel model;
in this embodiment, each laser point in the abnormal data set Q is mapped to a 3D grid, and then each voxel is assigned a value according to the intensity attribute of the laser point contained in the voxel. Wherein, the voxel containing the laser points is assigned as the laser point intensity value (if a plurality of laser points are contained, the intensity value of the laser point with the lowest elevation is assigned), the voxel not containing the laser points is assigned as 0, and the assignment of each voxel is further discretized to {0, \8230;, 255}, so as to obtain the value of each voxel. Therefore, a gray voxel model is obtained, and regularization of abnormal data set elimination is completed.
And step 3: and (3) detecting the ground voxel in the gray voxel model based on a three-dimensional connected region construction theory, wherein the specific flow is shown in fig. 5.
Step 3.1: based on the local elevation minimum characteristic of the ground point, taking the voxel with the lowest non-0 value in the gray voxel model as a seed voxel set V s Wherein s =1,2.
Step 3.1.1: and in the horizontal direction, the gray voxel model is partitioned by using the set grid size, and the voxel with the lowest elevation in each block and a non-0 value is taken as a seed voxel.
In this embodiment, the size of the mesh needs to be determined according to the complexity of the terrain, and if the terrain is more complex, the value of the mesh size is smaller, otherwise, the mesh size can be appropriately increased. However, the minimum mesh size cannot be smaller than the size of the largest target structure (e.g., building) in the original airborne LIDAR point cloud dataset, otherwise points inside the largest structure would be misinterpreted as ground seed voxels.
Step 3.1.2: searching in each block for height difference smaller than height difference threshold T e The non-0 value voxel is a ground seed, and an encrypted ground seed set V is obtained s Wherein s =1,2.;
in this embodiment, the height difference threshold T e The value of the constant is determined according to the spatial distribution condition of the original airborne LIDAR point cloud data.
Step 3.2: marking a ground seed voxel and a three-dimensional connected region formed by voxels which are three-dimensionally connected with the ground seed voxel and have similar gradient and gray level as a ground body metadata set, and finishing airborne LIDAR three-dimensional filtering based on a gray level voxel model;
step 3.2.1: calculating a frequency histogram of gray values of voxels of non-0 values in the gray voxel Model, fitting the frequency histogram with a Gaussian Mixture Model (GMM), and determining a gray value distribution range corresponding to the ground target.
In the present embodiment of the present invention,taking sample samp41 as an example, the frequency of the gray value of the non-0 value voxel in its corresponding V is counted and displayed in the form of a histogram, as shown in fig. 6. As can be seen from fig. 6, the gray scale distribution therein exhibits multimodality. If the multimodal distribution is assumed to be normal multimodal, the histogram in fig. 1 can be fitted with a gaussian mixture model to obtain distribution characteristic parameters (mean, standard deviation) of each gaussian distribution. The object of the invention is to filter the ground voxels, so that the distribution characteristic parameters mu, sigma and the distribution range [0, 20 ] of the Gaussian distribution corresponding to the ground]Is used for determining the gray scale range of the ground target, and the detailed scheme is as follows: diameter of order of μm l ×σ＝0，μ+m r X σ =20, m can be determined l And m r From this, the gray scale range of the ground target, i.e., [ mu ] -m σ, mu + m σ, can be determined]Wherein the multiplier m = max { m } l ，m r Get m from m l And m r Is determined based on the symmetry of the gaussian distribution. In addition, it should be noted that if there is a negative number in the grayscale range, the minimum value is set to 0 but 0 is not included.
Step 3.2.2: for any ground seed volume element V s Traversing the voxel V of the gray voxel model and the current ground seed voxel s Three-dimensional communication, the gray value is located in the gray range corresponding to the ground target, and the local terrain gradient is smaller than the gradient threshold value T s All unlabeled voxels of (1), and labeled L g Until all ground seed voxels V have been marked s A three-dimensional connected region of (a), a ground voxel set; the gradient threshold value is determined in a self-adaptive mode, and the specific method comprises the following steps: for the current seed voxel V s Detecting whether the space neighborhood contains marked ground voxels, if so, determining the maximum gradient value among the marked ground voxels as a terrain gradient threshold value T in the neighborhood s (ii) a Otherwise, the terrain slope threshold is set to 90 °.
In this embodiment, depth-first traversal of any seed voxel V in the grayscale voxel model s The specific method of three-dimensional connected regions of (2) is shown in fig. 7.
In this embodiment, based on V s And marked ground in its spatial vicinityThe formula of the voxel calculation gradient E is shown in formula (3):
wherein (r) p ，c p ，l p ) Is a V s Coordinates; (r) e ，c e ，l e ) And t is the marked ground voxel coordinate in the spatial neighborhood.
In this embodiment, different filtering results may be obtained by applying different neighborhood sizes in the above marking process. The best neighborhood scale will be determined experimentally.
The ground data obtained by the filtering method proposed by the invention is expressed by voxels, while the reference data is in the form of discrete LIDAR laser points. In order to compare with reference data and further evaluate the accuracy of the method provided by the invention, in this embodiment, firstly, laser points in ground voxels detected by the method are counted, then, the laser points are compared with the reference data, and then, the accuracy of the filtering algorithm is quantitatively evaluated by indexes such as type I errors (dividing ground point errors into non-ground point proportions), type II errors (dividing non-ground point errors into ground point proportions), total errors (dividing laser point proportions by errors) and Kappa coefficients.
In this embodiment, when the neighborhood dimensions are 6, 18, 26, 51, 56, and 64, respectively, the method of the present invention is applied to filter 15 test data, and Kappa coefficients of corresponding filtering results are shown in table 2. The data in this table is intended to examine the effect of different domain sizes on the accuracy of the filtering results and thus determine the optimal neighborhood size.
TABLE 2 Kappa coefficient of filtering methods for different neighborhood sizes
As can be seen from table 2, the average Kappa coefficients in the neighborhood of 6, 18, 26, 51, 56 and 64 are 71.64%, 78.73%, 80.81%, 84.34%, 83.29% and 83.28%, respectively. This indicates that: (1) The 51 neighborhood corresponds to the largest Kappa coefficient, and thus the 51 neighborhood is the best neighborhood size from the Kappa coefficient index. (2) An increase in the size of the neighborhood does not imply a necessary improvement in the filtering accuracy. This is because: information of the ground object may be propagated through connectivity, gray scale, and terrain slope strength similarities defined in the gray scale three-dimensional voxel data set. Taking the 6 neighborhoods as an example, the information of the ground object can only be propagated to the 6 voxels above, below, left, right, front and back thereof, thereby causing the seed voxel to be propagated only to the flat terrain around the seed voxel. If the neighborhood size can be increased, such as 26 neighborhoods, the increased propagation direction compared to 6 neighborhoods may include more ground voxels, thereby improving the accuracy of the filtering. This may explain why the Kappa coefficient for the filtering algorithm for the 26 neighborhood versus the 6 neighborhood is higher. However, if the neighborhood scale is too large, some non-ground voxels may be misclassified as ground voxels, resulting in a reduction in filtering accuracy. This may explain why the accuracy of the 56 neighborhood versus the 51 neighborhood is instead degraded.
In this embodiment, the method of the present invention is applied to perform quantitative evaluation on the filtering precision of 15 sample data in the neighborhood size of 51 using the reference data as a standard, and the precision of the filtering result is shown in table 3.
TABLE 3 precision of the filtering results
| Sample(s) | Type I error (%) | Type II error (%) | Total error (a)％) | Kappa coefficient (%) |
| samp11 | 13.49 | 19.80 | 16.18 | 66.86 |
| samp12 | 3.72 | 3.54 | 3.63 | 92.74 |
| samp21 | 1.12 | 3.47 | 1.64 | 95.25 |
| samp22 | 3.32 | 5.53 | 4.01 | 90.71 |
| samp23 | 2.78 | 5.37 | 4.00 | 91.96 |
| samp24 | 2.32 | 9.38 | 4.26 | 89.21 |
| samp31 | 0.94 | 1.94 | 1.40 | 97.18 |
| samp41 | 2.11 | 3.99 | 3.05 | 93.89 |
| samp42 | 1.24 | 5.99 | 0.79 | 98.10 |
| samp51 | 3.51 | 15.85 | 6.20 | 81.62 |
| samp52 | 1.03 | 31.12 | 4.19 | 75.28 |
| samp53 | 1.58 | 48.23 | 3.46 | 52.93 |
| samp54 | 1.91 | 6.88 | 4.58 | 90.83 |
| samp61 | 0.33 | 31.62 | 1.40 | 76.30 |
| samp71 | 2.36 | 20.28 | 4.39 | 77.95 |
As can be seen from table 3: the error distribution ranges of the class I error, the class II error and the total error are 0-14%, 1-49% and 0-17%; the average Kappa coefficient of the filtering can reach 84.34%. Thereby verifying the effectiveness of the method provided by the invention.
In this embodiment, the total error comparison results of the three-dimensional filtering method based on the binary voxel model, which are proposed by the method of the present invention and by Wang et al (bearing Wang, yan Xu, and Yu li. Intrinsic LIDAR point closed volume with its 3d group filtering application [ j ]. Photomeric Engineering and Remote Sensing,2017, 83 (2): 95-107), are shown in table 4.
TABLE 4 Total error comparison of the method of the present invention with a three-dimensional filtering method based on a binary voxel model
| Sample(s) | Total error of Wang et al method (%) | Total error of the method of the invention (%) |
| samp11 | 19.49 | 16.18 |
| samp12 | 4.02 | 3.63 |
| samp21 | 2.05 | 1.64 |
| samp22 | 4.97 | 4.01 |
| samp23 | 5.91 | 4.00 |
| samp24 | 6.34 | 4.26 |
| samp31 | 1.58 | 1.40 |
| samp41 | 2.17 | 3.05 |
| samp42 | 1.07 | 0.79 |
| samp51 | 8.09 | 6.20 |
| samp52 | 4.90 | 4.19 |
| samp53 | 3.46 | 3.46 |
| samp54 | 5.62 | 4.58 |
| samp61 | 1.93 | 1.40 |
| samp71 | 5.42 | 4.39 |
| Average out | 5.13 | 4.21 |
From table 4 it can be seen that: the total error of the method is obviously lower than that of the existing three-dimensional filtering algorithm based on the binary voxel model, namely the method has higher precision.
Finally, it should be noted that: the above examples are only intended to illustrate the technical solution of the present invention, but not to limit it; although the present invention has been described in detail with reference to the foregoing embodiments, it will be understood by those of ordinary skill in the art that: the technical solutions described in the foregoing embodiments may still be modified, or some or all of the technical features may be equivalently replaced; such modifications and substitutions do not depart from the spirit of the corresponding technical solutions and scope of the present invention as defined in the appended claims.
Claims (3)
1. An airborne LIDAR three-dimensional filtering method based on a gray voxel model is characterized in that: the method comprises the following steps:
step 1: reading original airborne LIDAR point cloud data to form an original airborne LIDAR point cloud data set;
step 2: regularizing an original airborne LIDAR point cloud dataset into a gray voxel model, wherein the gray is discretized representation of the average intensity of laser points in voxels, and the method specifically comprises the following steps:
step 2.1: rejecting abnormal data from an original airborne LIDAR point cloud data set to obtain a rejected abnormal data set;
step 2.2: regularizing the removed abnormal data set into a gray voxel model;
and 3, step 3: based on the theory of three-dimensional connected region construction, the ground voxel of the gray voxel model is detected, and the specific method comprises the following steps:
step 3.1: based on the local elevation minimum characteristic of the ground point, taking the voxel with the lowest non-0 value in the height from the gray voxel model as a ground seed voxel set;
step 3.1.1: in the horizontal direction, the gray voxel model is partitioned by using the set grid size, and the voxel with the lowest elevation in each block and a non-0 value is taken as a seed voxel;
step 3.1.2: searching for the height difference between each block and the ground seed less than the height difference threshold T e The non-0 value voxel is a ground seed voxel, and an encrypted ground seed voxel set V is obtained s Wherein s =1,2, \8230;
step 3.2: marking a ground seed voxel and a three-dimensional connected region formed by voxels which are three-dimensionally connected with the ground seed voxel and have similar gradient and gray level as a ground body metadata set, and finishing airborne LIDAR three-dimensional filtering based on a gray level voxel model;
step 3.2.1: calculating a frequency histogram of gray values of voxels with non-0 values in the gray voxel model, fitting the gray frequency histogram by using a Gaussian mixture model, and determining a gray value distribution range corresponding to a ground target;
step 3.2.2: for any ground seed volume element V s Traversing the voxel V of the gray voxel model and the current ground seed voxel s Three-dimensional communication, the gray value is located in the gray range corresponding to the ground target, and the local terrain gradient is smaller than the gradient threshold value T s All unlabeled voxels of (1), and labeled L g Until all the ground seed voxels V are marked s A three-dimensional connected region of (a), a ground voxel set;
the determination of the gradient threshold value in step 3.2.2 is adaptive, and the specific method is as follows: for the current seed voxel V s Detecting whether the space neighborhood contains marked ground voxels, if so, determining the maximum gradient value among the marked ground voxels as a terrain gradient threshold value T in the neighborhood s (ii) a Otherwise, the terrain grade threshold is set to 90 °.
2. The method of claim 1, wherein the airborne LIDAR three-dimensional filtering method based on a grayscale voxel model comprises: the specific method of the step 2.1 comprises the following steps:
step 2.1.1: counting the frequency of each laser point elevation value in the original airborne LIDAR point cloud data set, and displaying the counting result in a histogram visualization manner;
step 2.1.2: determining a maximum elevation threshold T corresponding to real terrain and ground features he And a minimum elevation threshold T le ；
Step 2.1.3: aiming at each laser point in the original airborne LIDAR point cloud data set, if the elevation value of each laser point is higher than the highest elevation threshold value T he Or below a minimum elevation threshold T le If the laser point is the elevation abnormal data, removing the laser point, otherwise, keeping the laser point, and obtaining a removed elevation abnormal data set;
step 2.1.4: counting and eliminating the frequency of the intensity value of each laser point in the elevation abnormal data set, and displaying the counting result in a histogram mode in a visualized manner;
step 2.1.5: determining a maximum intensity threshold T corresponding to the real terrain and the ground features hi And a minimum intensity threshold T li ；
Step 2.1.6: aiming at each laser point in the removed elevation abnormal data set, if the intensity value is higher than the highest intensity threshold value T hi Or below a minimum intensity threshold T hi And if not, retaining the laser point, and finally obtaining a removed elevation and intensity abnormal data set.
3. The method of claim 1, wherein the method comprises: the specific method of the step 2.2 comprises the following steps:
step 2.2.1: representing the spatial extent of the dataset by an axial parallel bounding box which rejects the abnormal dataset;
step 2.2.2: calculating the resolution of the volume element, namely the size of the volume element, wherein the resolution in the directions of x, y and z is determined according to the average point distance of the laser points in the abnormal data set removed;
step 2.2.3: dividing the axial parallel bounding boxes with the abnormal data sets removed according to the resolution ratios in the x direction, the y direction and the z direction to obtain three-dimensional grids, and calling each three-dimensional grid unit as a voxel;
step 2.2.4: and mapping each laser point in the abnormal data set to a three-dimensional grid, assigning values to each voxel according to the intensity attribute of the laser points contained in the voxel, and discretizing the values of each voxel to {0, \ 8230;, 255} to obtain a gray voxel model.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910027077.0A CN109785261B (en) | 2019-01-11 | 2019-01-11 | Airborne LIDAR three-dimensional filtering method based on gray voxel model |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910027077.0A CN109785261B (en) | 2019-01-11 | 2019-01-11 | Airborne LIDAR three-dimensional filtering method based on gray voxel model |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN109785261A CN109785261A (en) | 2019-05-21 |
| CN109785261B true CN109785261B (en) | 2023-03-17 |
Family
ID=66500206
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201910027077.0A Active CN109785261B (en) | 2019-01-11 | 2019-01-11 | Airborne LIDAR three-dimensional filtering method based on gray voxel model |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN109785261B (en) |
Families Citing this family (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110517220A (en) * | 2019-06-10 | 2019-11-29 | 长安大学 | A kind of surface of aggregate quantity detection method based on laser three-D data |
| CN112099046B (en) * | 2020-09-16 | 2023-05-16 | 辽宁工程技术大学 | Airborne LIDAR three-dimensional plane detection method based on multi-value voxel model |
Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN105844602A (en) * | 2016-04-01 | 2016-08-10 | 辽宁工程技术大学 | Airborne LIDAR point cloud 3D filtering method based on volume elements |
| CN106022259A (en) * | 2016-05-20 | 2016-10-12 | 江苏得得空间信息科技有限公司 | Laser-point cloud based method for extracting mountainous road by use of three-dimensional characteristic description model |
| CN106529469A (en) * | 2016-11-08 | 2017-03-22 | 华北水利水电大学 | Unmanned aerial vehicle airborne LiDAR point cloud filtering method based on adaptive gradient |
| CN108109139A (en) * | 2017-12-18 | 2018-06-01 | 辽宁工程技术大学 | Airborne LIDAR three-dimensional building object detecting method based on gray scale volume element model |
| US10002407B1 (en) * | 2017-08-11 | 2018-06-19 | Intermap Technologies Inc. | Method and apparatus for enhancing 3D model resolution |
Family Cites Families (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US8553942B2 (en) * | 2011-10-21 | 2013-10-08 | Navteq B.V. | Reimaging based on depthmap information |
| CN103679655A (en) * | 2013-12-02 | 2014-03-26 | 河海大学 | LiDAR point cloud filter method based on gradient and area growth |
| CN106157309B (en) * | 2016-07-04 | 2019-03-22 | 南京大学 | A kind of airborne LiDAR ground point cloud filtering method based on virtual seed point |
| CN108805154A (en) * | 2017-08-23 | 2018-11-13 | 辽宁工程技术大学 | A kind of geological fault recognition methods based on space clustering |
-
2019
- 2019-01-11 CN CN201910027077.0A patent/CN109785261B/en active Active
Patent Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN105844602A (en) * | 2016-04-01 | 2016-08-10 | 辽宁工程技术大学 | Airborne LIDAR point cloud 3D filtering method based on volume elements |
| CN106022259A (en) * | 2016-05-20 | 2016-10-12 | 江苏得得空间信息科技有限公司 | Laser-point cloud based method for extracting mountainous road by use of three-dimensional characteristic description model |
| CN106529469A (en) * | 2016-11-08 | 2017-03-22 | 华北水利水电大学 | Unmanned aerial vehicle airborne LiDAR point cloud filtering method based on adaptive gradient |
| US10002407B1 (en) * | 2017-08-11 | 2018-06-19 | Intermap Technologies Inc. | Method and apparatus for enhancing 3D model resolution |
| CN108109139A (en) * | 2017-12-18 | 2018-06-01 | 辽宁工程技术大学 | Airborne LIDAR three-dimensional building object detecting method based on gray scale volume element model |
Non-Patent Citations (3)
| Title |
|---|
| Sagi Filin 等.Segmentation of airborne laser scanning data using a slope adaptive neighborhood.2006,第60卷(第2期),71-80. * |
| 史建青.机载LiDAR在省级基础测绘中若干关键技术研究.2015,(第06期),A008-8. * |
| 王丽英 等.强度体元基元下的机载LiDAR 3D滤波.2019,第21卷(第12期),1945-1954. * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN109785261A (en) | 2019-05-21 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN108109139B (en) | Airborne LIDAR three-dimensional building detection method based on gray voxel model | |
| CN104240251B (en) | Multi-scale point cloud noise detection method based on density analysis | |
| CN105844602B (en) | A kind of airborne LIDAR point cloud three-dimensional filtering method based on volume elements | |
| CN111145174B (en) | 3D target detection method for point cloud screening based on image semantic features | |
| CN104573705B (en) | A kind of clustering method of building Point Cloud of Laser Scanner | |
| CN108074232B (en) | Voxel segmentation-based airborne LIDAR building detection method | |
| CN111340875B (en) | Space moving target detection method based on three-dimensional laser radar | |
| CN102944174A (en) | Point cloud data processing method and system | |
| CN104049245A (en) | Urban building change detection method based on LiDAR point cloud spatial difference analysis | |
| CN105260737A (en) | Automatic laser scanning data physical plane extraction method with multi-scale characteristics fused | |
| CN112099046B (en) | Airborne LIDAR three-dimensional plane detection method based on multi-value voxel model | |
| CN106157309A (en) | A kind of airborne LiDAR ground point cloud filtering method based on virtual Seed Points | |
| CN110047036B (en) | Polar grid-based ground laser scanning data building facade extraction method | |
| CN111008989B (en) | Airborne multispectral LIDAR three-dimensional segmentation method based on multivalued voxels | |
| CN109785261B (en) | Airborne LIDAR three-dimensional filtering method based on gray voxel model | |
| CN106778749A (en) | Based on the touring operating area boundary extraction method that concentration class and Delaunay triangles are reconstructed | |
| CN106952242A (en) | A kind of progressive TIN point cloud filtering method based on voxel | |
| CN112986964A (en) | Photon counting laser point cloud self-adaptive denoising method based on noise neighborhood density | |
| CN108765446B (en) | Power line point cloud segmentation method and system based on random field and random forest | |
| CN115861566A (en) | LiDAR point cloud building integrity enhancement method based on intensity voxel model | |
| CN115063698A (en) | Automatic identification and information extraction method and system for slope surface deformation crack | |
| Zhang et al. | Automatic extraction of DTM from low resolution DSM by twosteps semi-global filtering | |
| CN113836484A (en) | Self-adaptive point cloud rarefying method based on path point adjacent domain and ground filtering | |
| CN117292181A (en) | Sheet metal part hole group classification and full-size measurement method based on 3D point cloud processing | |
| CN117036971A (en) | Method for extracting airborne LiDAR data building under self-adaptive local spatial spectrum consistency |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |