CN113344808A - Terrain-adaptive specified-density airborne LiDAR point cloud simplification method - Google Patents
Terrain-adaptive specified-density airborne LiDAR point cloud simplification method Download PDFInfo
- Publication number
- CN113344808A CN113344808A CN202110582542.4A CN202110582542A CN113344808A CN 113344808 A CN113344808 A CN 113344808A CN 202110582542 A CN202110582542 A CN 202110582542A CN 113344808 A CN113344808 A CN 113344808A
- Authority
- CN
- China
- Prior art keywords
- points
- point cloud
- grid
- elevation
- point
- 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.)
- Granted
Links
- 238000000034 method Methods 0.000 title claims abstract description 56
- 230000003044 adaptive effect Effects 0.000 claims abstract description 10
- 238000012545 processing Methods 0.000 claims abstract description 10
- 239000011159 matrix material Substances 0.000 claims description 21
- 230000009467 reduction Effects 0.000 claims description 8
- 238000005070 sampling Methods 0.000 claims description 8
- 230000003746 surface roughness Effects 0.000 claims description 8
- 238000004364 calculation method Methods 0.000 claims description 6
- 239000013256 coordination polymer Substances 0.000 claims description 6
- 150000001875 compounds Chemical class 0.000 claims description 3
- 230000008520 organization Effects 0.000 claims description 3
- 239000000126 substance Substances 0.000 claims description 3
- 238000004519 manufacturing process Methods 0.000 abstract description 6
- 238000002474 experimental method Methods 0.000 abstract description 2
- 238000010586 diagram Methods 0.000 description 5
- 230000009897 systematic effect Effects 0.000 description 3
- 238000000354 decomposition reaction Methods 0.000 description 2
- 230000007547 defect Effects 0.000 description 2
- 230000000694 effects Effects 0.000 description 2
- 238000007726 management method Methods 0.000 description 2
- 238000004321 preservation Methods 0.000 description 2
- 238000012887 quadratic function Methods 0.000 description 2
- 230000037303 wrinkles Effects 0.000 description 2
- ATJFFYVFTNAWJD-UHFFFAOYSA-N Tin Chemical compound [Sn] ATJFFYVFTNAWJD-UHFFFAOYSA-N 0.000 description 1
- 238000007792 addition Methods 0.000 description 1
- 238000005452 bending Methods 0.000 description 1
- 230000000052 comparative effect Effects 0.000 description 1
- 238000013500 data storage Methods 0.000 description 1
- 239000006185 dispersion Substances 0.000 description 1
- 238000009826 distribution Methods 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 230000014759 maintenance of location Effects 0.000 description 1
- 238000013507 mapping Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 230000008569 process Effects 0.000 description 1
- 238000011160 research Methods 0.000 description 1
Images
Classifications
-
- G06T5/90—
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/60—Analysis of geometric attributes
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
- Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
- Y02A90/10—Information and communication technologies [ICT] supporting adaptation to climate change, e.g. for weather forecasting or climate simulation
Abstract
The invention relates to a terrain adaptive specified density airborne LiDAR point cloud simplification method. The method comprises the steps of dividing an original point cloud into a virtual grid, dividing the grid into three different surface types including a flat area, a rough area and a non-rough area according to the undulation degree and roughness of the interior of the grid, and then carrying out different thinning strategies on the different surface types to obtain the point cloud with the specified density. The method adopts different thinning strategies aiming at different types of regions with different tables, can fully reserve the key points of topographic relief, prevent the generation of too low local density and cavities, and avoid the distortion of the ground model caused by the missing of the key points. The method has the advantages of simple parameter setting, less parameter quantity and certain terrain self-adaptive capacity, and can effectively solve the problem of difficult threshold setting caused by the mixed point cloud data of different surface morphologies in the same area in actual production. Experiments prove that the method has high speed and high precision for processing the point cloud data with large data volume and can meet the actual production requirement.
Description
Technical Field
The invention belongs to the field of LiDAR data processing and application, and particularly relates to a terrain-adaptive specified-density airborne LiDAR point cloud simplification method.
Background
With the continued advancement of hardware technology, the density of acquired LiDAR point cloud data has reached tens of points per square meter. Although the high-density point cloud is helpful for more detailed description of topographic features, a large amount of data redundancy is brought to a relatively flat area, and the redundant data directly influences the efficiency of data storage, processing and display speed, so that the existing point cloud processing algorithm and software are challenged, and the requirement on hardware is also increased. The point cloud rarefying is to regularly screen the original point cloud on the basis of ensuring to meet the production precision or the actual requirement, and to furthest reduce the number of the point cloud. According to the rules adopted by different algorithms, point cloud thinning can be summarized into a random sampling algorithm and an algorithm considering terrain features. The random sampling algorithm mainly comprises a rarefaction algorithm based on a regular grid and a rarefaction algorithm based on a system. The former divides the entire data range into a large number of small grids, and randomly selects one data point or an average of several data points in the grid to replace all data in the grid. The latter is that when loading the original point cloud, setting the sampling interval N, and taking the first of N point cloud data or randomly extracting a data point from the N point clouds. The two algorithms are simple and quick as well, but the topographic features cannot be well reserved, so that the precision of point cloud data is reduced. In the algorithm considering the terrain features, parameters such as elevation, altitude difference, curvature, gradient and the like are mostly adopted. There are the ddr (data Density reduction) rarefaction algorithm proposed by Pamelas, the rarefaction algorithm based on point cloud dispersion proposed by slow scene et al, the rarefaction algorithm based on gradient proposed by zhiyu et al, the constrained TIN node rarefaction method proposed by yangming et al, and the rarefaction algorithm adopting poisson disk sampling in geodesic space proposed by hough et al. In the production process of mapping products based on point cloud data, how to keep key points in the point cloud data while reducing the amount and redundancy of the point cloud data and enable the precision of the point cloud data to reach the standard becomes the key point of research.
Disclosure of Invention
Aiming at the defects of the prior art, the invention provides a terrain self-adaptive specified-density airborne LiDAR point cloud simplification method. Dividing a virtual grid into original point clouds, calculating the topographic relief degree in the grid, dividing the grid into a flat area and a non-flat area according to a topographic relief degree threshold value, then calculating the topographic roughness of the grid in the non-flat area, further dividing the grid into a rough area and a non-rough area according to a set threshold value, and adopting a random thinning method with specified density aiming at the flat area; aiming at the non-flat area, calculating the cloud curvature of points in the grid, reserving high curvature points and elevation extreme points, and obtaining point cloud with the specified density by adopting different thinning strategies for different areas of the non-flat area according to the size relation between the reserved high curvature points and the number of the elevation extreme points and the number of the points required by the single grid under the specified density.
In order to achieve the aim, the technical scheme provided by the invention is a terrain adaptive specified density airborne LiDAR point cloud simplification method, which comprises the following steps:
step 1, establishing a two-dimensional grid point cloud index;
step 2, calculating the topographic relief degree in the grid;
step 3, dividing the grid into a flat area and a non-flat area according to the topographic relief degree threshold value alpha;
step 4, calculating the roughness of the terrain in the grid aiming at the non-flat area;
step 5, further dividing the non-flat area into a rough area and a non-rough area according to the roughness threshold beta;
step 6, aiming at the surface grid of the flat area, calculating the number of points required by the grid according to the specified density, and obtaining point cloud data after rarefaction by adopting random downsampling;
step 7, calculating the curvature of the point cloud in the grid network aiming at the point cloud in the non-flat area, and reserving an elevation extreme point and a high curvature point;
and 8, adopting different thinning strategies for different areas of the non-flat area according to the relationship between the number of points reserved in the step 7 and the number of points required by the single grid under the specified density.
And in the step 1, the point cloud is projected to a two-dimensional plane according to XY coordinates of the point cloud data, a regular grid is divided, and the organization and management of the point cloud data are realized by using the row and column numbers of the grid.
In step 2, the topographic relief degree is a difference between a maximum elevation value and a minimum elevation value in a local range, and the topographic relief degree in the local range is reflected, and the topographic relief degree is calculated by the following method:
dH＝Hmax-Hmin (1)
in the formula (d)HIs the local extent of relief, Hmax、HminRespectively, the maximum value and the minimum value of the elevation of the surface point in the local range.
In step 3, the region with the topographic relief degree greater than α in the grid is marked as a flat region, the region with the topographic relief degree less than α is marked as a non-flat region, and α is an empirical value.
In step 4, the surface roughness represents the roughness or the wrinkle degree of the local terrain, and the surface roughness is calculated based on the elevation fitting error of each ground point in the local range, wherein the calculation method comprises the following steps:
wherein Rough is surface roughness, ZiIs the measured elevation of the ith point, n is the number of points in the local range, Z (x)i,yi) For the fitting elevation of the ith point on the quadric surface, the quadric surface can be obtained by performing quadratic fitting on local point cloud data in the grid, and the calculation mode is as follows:
Z＝a0+a1x+a2y+a3x2+a4y2+a5xy (3)
wherein Z represents the elevation of a point in the local area, x and y are the plane coordinates of the point in the local area, a0、a1、a2、a3、a4And a5Are coefficients of a quadratic function.
For the grid area with redundant data, the number of points is generally more than 6, redundant observation exists, so the least square method is adopted for solving, namely the sum of squares of the error between the fitted elevation and the actually measured elevation is the minimum, and the formula is as follows:
in the formula, Z (x)i,yi) As the fitted elevation of the ith point, ZiAnd n is the number of points in the local range.
Establishing an error equation, and solving a coefficient matrix of a quadric function, wherein the formula is as follows:
V＝BX-L (5)
wherein the content of the first and second substances,
X＝(BTB)-1BTL (6)
in the formula, V represents an error matrix between the fitting elevation and the actually measured elevation, B represents a coordinate combination term corresponding to each coefficient in formula (3), X is a coefficient matrix of a quadric function to be solved in formula (3), and L is an actually measured elevation matrix.
In step 5, the area with the terrain roughness greater than β in the grid is marked as a rough area, the area with the terrain roughness less than β is marked as a non-rough area, and β is an empirical value.
Moreover, the local elevation extreme points in step 7 include two elevation maximum points and minimum points in the grid, which can be directly calculated according to the elevation of discrete points, where a high curvature point is a point with a curvature value greater than a threshold γ, γ is an empirical value, and a target point P is assumed to bei＝(xi,yi,zi)THas a covariance matrix of CPiAnd then:
in the formula (I), the compound is shown in the specification,
is the central position of the target point and its neighborhood points, N is the number of neighborhood points of the target point, CPiIs a symmetric semi-positive definite matrix.
Performing eigenvalue decomposition on the covariance matrix to obtain three eigenvalues lambda1、λ2、λ3(λ1≥λ2≥λ3) And its corresponding feature vector e1、e2、e3。e3For the estimated normal vector, λ3Describing the variation of the curved surface along the normal vector direction, the curvature curve is calculated as follows:
although the accuracy of curvature estimation based on eigenvalues is slightly lower, the processing speed is faster, and the method has better applicability to large-range redundant airborne LiDAR ground point cloud data, so that the curvature characteristics of the laser foot points are estimated by adopting the formula (9), and the neighborhood type is selected as the K neighborhood.
And in step 8, if the number of the preserved elevation extreme points and the high curvature points is greater than the number of the points required by the unit grid under the specified density, preserving the high curvature points and the elevation extreme points to the specified number in the order from the large curvature to the small curvature for the non-rough area, and preserving all the curvature points and the elevation extreme points for the rough area. And if the number of the elevation extreme points and the high curvature points is less than the number of the points required by the unit grid under the specified density, calculating the point number difference between the elevation extreme points and the high curvature points and filling the point number difference and the required point number in the rest points in a random sampling mode. And obtaining the point cloud after rarefaction after processing all the point cloud data.
Compared with the prior art, the invention has the following advantages:
(1) the accuracy is as follows: the method adopts different thinning strategies aiming at different terrain areas, fully retains the key points of terrain fluctuation, prevents the generation of too low local density and cavities, and simultaneously avoids the distortion of a ground model caused by the loss of the key points;
(2) ease of use: the method has the advantages that the parameter setting is simple, the number of parameters is small, the method has certain terrain self-adaption capability, and the problem of difficult threshold setting caused by the fact that point cloud data of different surface forms in the same area are mixed in actual production can be effectively solved;
(3) the method is high in efficiency, and experiments prove that the method is high in processing speed and high in precision aiming at the point cloud data with large data volume, and can meet the actual production requirements.
Drawings
FIG. 1 is a technical flow diagram of an embodiment of the present invention.
Fig. 2 is a comparison diagram of experimental results of different methods for mountain region data according to an embodiment of the present invention, where fig. 2(a) is original point cloud data, fig. 2(b) is point cloud data after being thinned out by a system thinning-out method, fig. 2(c) is point cloud data after being thinned out by a curvature feature global method, and fig. 2(d) is point cloud data after being thinned out by the method of the present invention.
Fig. 3 is a comparison diagram of local detail of experimental results of mountain region data according to an embodiment of the present invention, where fig. 3(a) is a local detail diagram of point cloud data after being thinned out by a system thinning-out method, and fig. 3(b) is a local detail diagram of point cloud data after being thinned out by a method according to the present invention.
Fig. 4 is a comparison graph of local details of the experimental results of the urban area data 1 according to the embodiment of the present invention, where fig. 4(a) is the original point cloud data, and fig. 4(b) is the point cloud data after thinning by the method of the present invention.
Fig. 5 is a comparison graph of local details of the experimental results of the urban area data 2 according to the embodiment of the present invention, in which fig. 5(a) is the original point cloud data, and fig. 5(b) is the point cloud data after thinning by the method of the present invention.
Fig. 6 is a comparison graph of the local detail of the experimental result of the hilly area data 1 according to the embodiment of the present invention, wherein fig. 6(a) is the original point cloud data, and fig. 6(b) is the point cloud data after thinning by the method of the present invention.
Fig. 7 is a comparison graph of the local detail of the experimental result of the hilly area data 1 according to the embodiment of the present invention, wherein fig. 7(a) is the original point cloud data, and fig. 7(b) is the point cloud data after thinning by the method of the present invention.
Fig. 8 is a comparison graph of the local detail of the experimental result of the hilly area data 2 according to the embodiment of the present invention, in which fig. 8(a) is the original point cloud data, and fig. 8(b) is the point cloud data after thinning by the method of the present invention.
Fig. 9 is a comparison graph of the local detail of the experimental result of the hilly area data 2 according to the embodiment of the present invention, in which fig. 9(a) is the original point cloud data, and fig. 9(b) is the point cloud data after being thinned out by the method of the present invention.
Fig. 10 is a comparison graph of the local detail of the experimental result of the hilly area data 3 according to the embodiment of the present invention, in which fig. 10(a) is the original point cloud data, and fig. 10(b) is the point cloud data after being thinned out by the method of the present invention.
Detailed Description
Aiming at the defects of the prior art, the invention provides a terrain self-adaptive specified-density airborne LiDAR point cloud simplification method. Dividing a virtual grid into original point clouds, calculating the topographic relief degree in the grid, dividing the grid into a flat area and a non-flat area according to a topographic relief degree threshold value, then calculating the topographic roughness of the grid in the non-flat area, further dividing the grid into a rough area and a non-rough area according to a set threshold value, and adopting a random thinning method with specified density aiming at the flat area; aiming at the non-flat area, calculating the cloud curvature of points in the grid, reserving high curvature points and elevation extreme points, and obtaining point cloud with the specified density by adopting different thinning strategies for different areas of the non-flat area according to the size relation between the reserved high curvature points and the number of the elevation extreme points and the number of the points required by the single grid under the specified density.
The technical solution of the present invention is further explained with reference to the drawings and the embodiments.
As shown in fig. 1, the process of the embodiment of the present invention includes the following steps:
step 1, establishing a two-dimensional grid point cloud index. And projecting the point cloud to a two-dimensional plane according to XY coordinates of the point cloud data, dividing a regular grid, and realizing the organization and management of the point cloud data by using the row and column numbers of the grid.
And 2, calculating the undulation degree of the terrain in the grid. The topographic relief degree is the difference between the maximum elevation value and the minimum elevation value of a local range and reflects the topographic relief of the local range, and the topographic relief degree is calculated by the following method:
dH＝Hmax-Hmin (1)
in the formula (d)HIs the local extent of relief, Hmax、HminRespectively, the maximum value and the minimum value of the elevation of the surface point in the local range.
And 3, dividing the grid into a flat area and a non-flat area according to the topographic relief degree threshold alpha. If the topographic relief degree in the grid is larger than alpha, the grid is marked as a flat area, otherwise, the grid is marked as a non-flat area, and alpha is an empirical value.
And 4, calculating the terrain roughness of the grid aiming at the non-flat area. The surface roughness expresses the roughness or the wrinkle degree of a local terrain, and the surface roughness is calculated on the basis of elevation fitting errors of each surface point in a local range, wherein the calculation method comprises the following steps:
wherein Rough is surface roughness, ZiIs the measured elevation of the ith point, n is the number of points in the local range, Z (x)i,yi) For the fitting elevation of the ith point on the quadric surface, the quadric surface can be obtained by performing quadratic fitting on local point cloud data in the grid, and the calculation mode is as follows:
Z＝a0+a1x+a2y+a3x2+a4y2+a5xy (3)
wherein Z represents the elevation of a point in the local area, x and y are the plane coordinates of the point in the local area, a0、a1、a2、a3、a4And a5Are coefficients of a quadratic function.
For the grid area with redundant data, the number of points is generally more than 6, redundant observation exists, so the least square method is adopted for solving, namely the sum of squares of the error between the fitted elevation and the actually measured elevation is the minimum, and the formula is as follows:
in the formula, Z (x)i,yi) As the fitted elevation of the ith point, ZiAnd n is the number of points in the local range.
Establishing an error equation, and solving a coefficient matrix of a quadric function, wherein the formula is as follows:
V＝BX-L (5)
wherein the content of the first and second substances,
X＝(BTB)-1BTL (6)
in the formula, V represents an error matrix between the fitting elevation and the actually measured elevation, B represents a coordinate combination term corresponding to each coefficient in formula (3), X is a coefficient matrix of a quadric function to be solved in formula (3), and L is an actually measured elevation matrix.
And 5, further dividing the non-flat area into a rough area and a non-rough area according to the roughness threshold beta. If the roughness of the terrain in the grid is larger than beta, marking the terrain as a rough area, otherwise marking the terrain as a non-rough area, wherein the beta is an empirical value.
And 6, aiming at the surface grid of the flat area, calculating the number of points required by the grid according to the specified density, and obtaining the point cloud data after rarefaction by adopting random downsampling.
And 7, calculating the curvature of the point cloud in the grid aiming at the point cloud in the non-flat area, and reserving the elevation extreme point and the high curvature point.
The local elevation extreme points comprise elevation maximum points and minimum points in the grid, reflect the variation range of the elevation of the local area, and can be directly calculated according to the elevation of discrete points.
The curvature of the point reflects the bending degree of the terrain curved surface at the point, so that the high curvature point is an important characteristic point for expressing that local terrain is obviously changed, and the set of the high curvature point can depict a key terrain skeleton of the whole area. Taking the point with curvature value larger than threshold gamma as high curvature point, and assuming target point Pi＝(xi,yi,zi)THas a covariance matrix of CPiThen, then：
In the formula (I), the compound is shown in the specification,
is the central position of the target point and its neighborhood points, N is the number of neighborhood points of the target point, CPiIs a symmetric semi-positive definite matrix.
Performing eigenvalue decomposition on the covariance matrix to obtain three eigenvalues lambda1、λ2、λ3(λ1≥λ2≥λ3) And its corresponding feature vector is e1、e2、e3。e3For the estimated normal vector, λ3Describing the variation of the curved surface along the normal vector direction, and calculating the curvature curve as follows:
although the accuracy of curvature estimation based on eigenvalues is slightly lower, the processing speed is faster, with better applicability to large-scale redundant airborne LiDAR ground point cloud data. Therefore, the curvature characteristics of the laser foot points are estimated by adopting a formula (9), and K neighborhoods are selected according to the neighborhood types.
And 8, adopting different thinning strategies for different areas of the non-flat area according to the relationship between the number of points reserved in the step 7 and the number of points required by the single grid under the specified density.
If the number of the preserved elevation extreme points and the high curvature points is more than the number of the points required by the single grid under the specified density, the high curvature points and the elevation extreme points are preserved to the specified number of the high curvature points and the elevation extreme points in the non-rough area from the large curvature to the small curvature, and all the curvature points and the elevation extreme points are preserved to the rough area. And if the number of the elevation extreme points and the high curvature points is less than the number of the points required by the unit grid under the specified density, calculating the point number difference between the elevation extreme points and the high curvature points and filling the point number difference and the required point number in the rest points in a random sampling mode. And obtaining the point cloud after rarefaction after processing all the point cloud data.
FIG. 2 is a graph of comparative results of different thinning scenarios, wherein the results of the systematic thinning method, while substantially preserving the overall topographical features, do not characterize a portion of the valley regions; although a plurality of topographic feature points are reserved in the method for carrying out global thinning based on curvature features, a large-range data hole appears; compared with the prior art, the thinning method provided by the invention not only well keeps the topographic characteristics, especially the characteristics of partial valley areas are more obvious after thinning, but also avoids the occurrence of data holes, and the obtained thinning effect is best.
Fig. 3 is a comparison graph of local details of point cloud data after the systematic thinning method and the method of the present invention are thinned, and it can be found by observing the local details that the thinning method of the present invention identifies points in the valley region as important topographic feature points and performs targeted preservation, while the systematic thinning method is interval sampling and the preservation of the local topographic feature points is random and blind. Under the condition of the same point cloud retention rate, the thinning method can give consideration to the terrain features, the capability of retaining key terrain feature points is stronger, and the point location distribution is more reasonable.
Fig. 4-10 show the experimental results of various point cloud data, which shows that the invention has good experimental effects, no matter the building contour points in urban areas or the key points in the terrain feature areas such as ridges, valleys and hills in mountainous areas are well preserved, and meanwhile, the number of point clouds is greatly reduced, the redundancy of data is greatly reduced, and the invention is proved to have wide applicability and practicability.
The specific embodiments described herein are merely illustrative of the spirit of the invention. Various modifications or additions may be made to the described embodiments or alternatives may be employed by those skilled in the art without departing from the spirit or ambit of the invention as defined in the appended claims.
Claims (8)
1. A terrain adaptive specified density airborne LiDAR point cloud simplification method is characterized by comprising the following steps:
step 1, establishing a two-dimensional grid point cloud index;
step 2, calculating the topographic relief degree in the grid;
step 3, dividing the grid into a flat area and a non-flat area according to the topographic relief degree threshold value alpha;
step 4, calculating the roughness of the terrain in the grid aiming at the non-flat area;
step 5, further dividing the non-flat area into a rough area and a non-rough area according to the roughness threshold beta;
step 6, aiming at the surface grid of the flat area, calculating the number of points required by the grid according to the specified density, and obtaining point cloud data after rarefaction by adopting random downsampling;
step 7, calculating the curvature of the point cloud in the grid network aiming at the point cloud in the non-flat area, and reserving an elevation extreme point and a high curvature point;
and 8, adopting different thinning strategies for different areas of the non-flat area according to the relationship between the number of points reserved in the step 7 and the number of points required by the single grid under the specified density.
2. The terrain adaptive, density-specified airborne LiDAR point cloud reduction method of claim 1, wherein: and step 1, projecting the point cloud to a two-dimensional plane according to XY coordinates of the point cloud data, dividing a regular grid, and realizing organization and management of the point cloud data by using row and column numbers of the grid.
3. The terrain adaptive, density-specified airborne LiDAR point cloud reduction method of claim 1, wherein: the topographic relief degree in the step 2 is the difference between the maximum elevation value and the minimum elevation value of the local range, and reflects the topographic relief of the local range, and the topographic relief degree calculation method comprises the following steps:
dH＝Hmax-Hmin (1)
in the formula (d)HIs the local extent of relief, Hmax、HminRespectively, the maximum value and the minimum value of the elevation of the surface point in the local range.
4. The terrain adaptive, density-specified airborne LiDAR point cloud reduction method of claim 1, wherein: in the step 3, the area with the topographic relief degree larger than alpha in the grid is marked as a flat area, the area with the topographic relief degree smaller than alpha is marked as a non-flat area, and alpha is an empirical value.
5. The terrain adaptive, density-specified airborne LiDAR point cloud reduction method of claim 1, wherein: the method for calculating the surface roughness in the step 4 comprises the following steps:
wherein Rough is surface roughness, ZiIs the measured elevation of the ith point, n is the number of points in the local range, Z (x)i,yi) For the fitting elevation of the ith point on the quadric surface, the quadric surface can be obtained by performing quadratic fitting on local point cloud data in the grid, and the calculation mode is as follows:
Z＝a0+a1x+a2y+a3x2+a4y2+a5xy (3)
wherein Z represents the elevation of a point in the local area, x and y are the plane coordinates of the point in the local area, a0、a1、a2、a3、a4And a5Is the coefficient of the quadric function;
for the grid area with data redundancy, a least square method is adopted for solving, namely the sum of squares of the error between the fitted elevation and the actually measured elevation is the minimum, and the formula is as follows:
in the formula, Z (x)i,yi) As the fitted elevation of the ith point, ZiThe measured elevation of the ith point is obtained, and n is the number of points in a local range;
establishing an error equation, and solving a coefficient matrix of a quadric function, wherein the formula is as follows:
V＝BX-L (5)
wherein the content of the first and second substances,
X＝(BTB)-1BTL (6)
in the formula, V represents an error matrix between the fitting elevation and the actually measured elevation, B represents a coordinate combination term corresponding to each coefficient in formula (3), X is a coefficient matrix of a quadric function to be solved in formula (3), and L is an actually measured elevation matrix.
6. The terrain adaptive, density-specified airborne LiDAR point cloud reduction method of claim 1, wherein: in the step 5, the region with the terrain roughness larger than beta in the grid is marked as a rough region, the region with the terrain roughness smaller than beta is marked as a non-rough region, and beta is an empirical value.
7. The terrain adaptive, density-specified airborne LiDAR point cloud reduction method of claim 1, wherein: the local elevation extreme points in the step 7 comprise elevation maximum points and minimum points in the grid, and can be directly calculated according to the elevations of discrete points; the high curvature point is a point where the curvature value is greater than a threshold value γ, γ being an empirical value, assuming a target point Pi＝(xi,yi,zi)THas a covariance matrix of CPiAnd then:
in the formula (I), the compound is shown in the specification,
is the central position of the target point and its neighborhood points, N is the number of neighborhood points of the target point, CPiIs a symmetric semi-positive definite matrix;
the curvature curve is calculated as follows:
in the formula, λ1、λ2、λ3Is the three eigenvalues of the covariance matrix, and1≥λ2≥λ3。
8. the terrain adaptive, density-specified airborne LiDAR point cloud reduction method of claim 1, wherein: in the step 8, if the number of the preserved elevation extreme points and the high curvature points is more than the number of the points required by the single grid under the specified density, the high curvature points and the elevation extreme points are preserved to the specified number in the order of the curvature from large to small in the non-rough area, and all the curvature points and the elevation extreme points are preserved in the rough area; and if the number of the elevation extreme points and the high curvature points is less than that of the points required by the unit grid under the specified density, calculating the point number difference between the elevation extreme points and the high curvature points and the required points, filling up the rest points in a random sampling mode, and obtaining the point cloud after rarefaction after processing all the point cloud data.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110582542.4A CN113344808B (en) | 2021-05-27 | 2021-05-27 | Terrain-adaptive specified-density airborne LiDAR point cloud simplification method |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110582542.4A CN113344808B (en) | 2021-05-27 | 2021-05-27 | Terrain-adaptive specified-density airborne LiDAR point cloud simplification method |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN113344808A true CN113344808A (en) | 2021-09-03 |
| CN113344808B CN113344808B (en) | 2022-06-07 |
Family
ID=77471669
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202110582542.4A Active CN113344808B (en) | 2021-05-27 | 2021-05-27 | Terrain-adaptive specified-density airborne LiDAR point cloud simplification method |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN113344808B (en) |
Cited By (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN114299235A (en) * | 2021-12-31 | 2022-04-08 | 中铁二院工程集团有限责任公司 | DOM (document object model) manufacturing method based on color point cloud |
| CN114332631A (en) * | 2022-01-12 | 2022-04-12 | 中铁二院工程集团有限责任公司 | LiDAR point cloud data extraction method suitable for mountainous area dangerous rockfall |
| CN114648621A (en) * | 2022-04-06 | 2022-06-21 | 重庆市勘测院(重庆市地图编制中心) | Ground point cloud rapid filtering method, device, equipment and storage medium |
Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20130144565A1 (en) * | 2011-12-01 | 2013-06-06 | Harris Corporation | Accuracy-based significant point derivation from dense 3d point clouds for terrain modeling |
| CN110853139A (en) * | 2019-07-02 | 2020-02-28 | 中国人民解放军战略支援部队信息工程大学 | Multi-beam sounding data simplification method and device |
| CN110992479A (en) * | 2020-03-05 | 2020-04-10 | 浙江交工集团股份有限公司 | High-roughness three-dimensional curved surface fitting method suitable for scattered point clouds |
| CN111369436A (en) * | 2020-02-27 | 2020-07-03 | 山东科技大学 | Airborne LiDAR point cloud rarefying method considering multi-terrain features |
-
2021
- 2021-05-27 CN CN202110582542.4A patent/CN113344808B/en active Active
Patent Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20130144565A1 (en) * | 2011-12-01 | 2013-06-06 | Harris Corporation | Accuracy-based significant point derivation from dense 3d point clouds for terrain modeling |
| CN110853139A (en) * | 2019-07-02 | 2020-02-28 | 中国人民解放军战略支援部队信息工程大学 | Multi-beam sounding data simplification method and device |
| CN111369436A (en) * | 2020-02-27 | 2020-07-03 | 山东科技大学 | Airborne LiDAR point cloud rarefying method considering multi-terrain features |
| CN110992479A (en) * | 2020-03-05 | 2020-04-10 | 浙江交工集团股份有限公司 | High-roughness three-dimensional curved surface fitting method suitable for scattered point clouds |
Non-Patent Citations (3)
| Title |
|---|
| 杜浩等: "顾及地形特征的LiDAR点云数据抽稀算法", 《测绘科学》 * |
| 钱金菊等: "机载LiDAR点云数据抽稀算法研究述评", 《测绘通报》 * |
| 陈佩奇等: "一种城区LiDAR点云数据的抽稀算法", 《遥感技术与应用》 * |
Cited By (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN114299235A (en) * | 2021-12-31 | 2022-04-08 | 中铁二院工程集团有限责任公司 | DOM (document object model) manufacturing method based on color point cloud |
| CN114332631A (en) * | 2022-01-12 | 2022-04-12 | 中铁二院工程集团有限责任公司 | LiDAR point cloud data extraction method suitable for mountainous area dangerous rockfall |
| CN114332631B (en) * | 2022-01-12 | 2023-04-18 | 中铁二院工程集团有限责任公司 | LiDAR point cloud data extraction method suitable for mountainous area dangerous rockfall |
| CN114648621A (en) * | 2022-04-06 | 2022-06-21 | 重庆市勘测院(重庆市地图编制中心) | Ground point cloud rapid filtering method, device, equipment and storage medium |
Also Published As
| Publication number | Publication date |
|---|---|
| CN113344808B (en) | 2022-06-07 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN113344808B (en) | Terrain-adaptive specified-density airborne LiDAR point cloud simplification method | |
| CN102663237B (en) | Point cloud data automatic filtering method based on grid segmentation and moving least square | |
| CN108830931B (en) | Laser point cloud simplification method based on dynamic grid k neighborhood search | |
| CN113781667B (en) | Three-dimensional structure simplified reconstruction method and device, computer equipment and storage medium | |
| CN109685080B (en) | Multi-scale plane extraction method based on Hough transformation and region growth | |
| CN103701466A (en) | Scattered point cloud compression algorithm based on feature reservation | |
| CN105118090A (en) | Adaptive point-cloud filtering method for complex terrain structure | |
| CN111340723B (en) | Terrain-adaptive airborne LiDAR point cloud regularization thin plate spline interpolation filtering method | |
| CN109410342A (en) | A kind of point cloud compressing method retaining boundary point | |
| CN113158892B (en) | Face recognition method irrelevant to textures and expressions | |
| CN115222625A (en) | Laser radar point cloud denoising method based on multi-scale noise | |
| CN110910462B (en) | Point cloud light weight method based on feature calculation and storage medium | |
| CN109671155A (en) | Surface mesh method for reconstructing, system and relevant device based on point cloud data | |
| CN116310849B (en) | Tree point cloud monomerization extraction method based on three-dimensional morphological characteristics | |
| CN112669462B (en) | Model processing method and system suitable for scanning point cloud | |
| CN112614216B (en) | Variable-curvature self-adaptive point cloud data down-sampling method | |
| CN112687002B (en) | Three-dimensional geological model grid optimization method | |
| CN116721228B (en) | Building elevation extraction method and system based on low-density point cloud | |
| CN112767429A (en) | Ground-snow surface point cloud rapid segmentation method | |
| CN113744389B (en) | Point cloud simplifying method for complex part curved surface feature preservation | |
| CN114228154B (en) | Gradient void structure modeling slicing method and system based on three-dimensional section characteristics | |
| CN113470177A (en) | Three-dimensional model geometric self-adaptive simplification method in GIS system | |
| CN115841538B (en) | Multivariable control DEM terrain feature line visual blanking method | |
| CN113192172B (en) | Airborne LiDAR ground point cloud simplification method | |
| CN117078754A (en) | Light-weight pose estimation method and system based on directivity feature selection |
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 |