CN111209828B - Method and system for extracting building roof vertex from airborne laser radar point cloud - Google Patents

Abstract

The invention discloses a method and a system for extracting building roof vertexes from airborne laser radar point clouds, wherein the method for extracting the building roof vertexes divides the airborne laser radar point clouds in a survey area into a plurality of clusters through three-dimensional Euclidean clustering, calculates the Laplacian of DSM for the cluster with larger area, and divides the points in the clusters into the building roof vertexes if the pixel occupation ratio of the Laplacian smaller than a flatness threshold value is higher; for clusters with smaller areas, points on the roof of the building are extracted by adopting an improved RANSAC algorithm. The method can solve the problems of high calculation complexity and low execution efficiency when the plane roof point of the building is extracted, can also avoid classifying non-building tree points into the roof top point of the building, and improves the efficiency and the accuracy of extracting the roof point of the building.

Description

Method and system for extracting building roof vertex from airborne laser radar point cloud

Technical Field

The invention belongs to the technical field of airborne laser radar point cloud data processing, and particularly relates to a method and a system for extracting points falling on a building roof from airborne laser radar point cloud.

Background

Airborne LiDAR (laser radar) is one of the most efficient And fastest-developed large-area surveying And mapping means in the current surveying And mapping field. The method is characterized in that the distance measurement is carried out by using a laser to emit and receive high-energy laser pulses, a GNSS (Global Navigation Satellite System) receiver gives a real-time position of the laser, an INS (Inertial Navigation System) gives a real-time three-dimensional attitude of the laser, and a three-dimensional coordinate of a scattering surface can be calculated by using a vector formula and coordinate conversion. By means of mechanical scanning devices, airborne LiDAR can achieve large area measurements. Currently, the transmission frequency of airborne LiDAR is commonly up to hundreds of KHZ and beyond, and even 2000 KHZ. This means that a lidar can acquire hundreds of thousands or even millions of points a second, without taking into account multi-echoes and no echoes. Thus, airborne LiDAR can obtain a large number of points, called point clouds, in a short amount of time.

In the preliminarily processed airborne LiDAR point cloud, point attributes provided include echo intensity, echo for the next time, data acquisition time, and the like, in addition to point coordinates. However, the system lacks semantic information, cannot give the physical properties of a scattering surface on which laser pulses are directed, and does not know whether the measurement is ground, buildings, vegetation, birds, or the like. While in constructing a DEM (Digital elevation model) points on the ground must be used, in constructing a three-dimensional building model points on the building must be used, and in forest mapping points on trees must be used. Therefore, classifying point clouds to determine on which basic ground object type the point is located is one of the most important fundamental jobs in airborne LiDAR data processing.

On-board LiDAR point clouds, particularly those covering residential areas, often contain points that fall on buildings. The points falling on the building are extracted, and the method can be used for three-dimensional modeling of the building, DSM (Digital Surface Model) generation, city planning and the like. Building point extraction is a hot problem for point cloud data processing.

Researchers have proposed many methods to extract building points from airborne LiDAR point clouds, both point cloud-based and image-based. Image-based methods (see xu literature, yangbishen, welfare, et al. LiDAR point cloud data building and crown extraction [ J ]. survey and drawing, 2013, 42 (1): 51-58.) in a multi-marker process interpolate the point cloud into an image, then typically extract the pixels containing the building by some image analysis method, thereby enabling the extraction of building points and the tracking of boundaries. Point cloud-based methods generally extract points located on a local plane from discrete point clouds as roof vertices using RANSAC (random sampling consistency, see patent CN105139379B) algorithm or Hough transform (see yuan chenxin, Hough transform-based building point cloud data feature extraction and three-dimensional modeling, 2019, the university of east china physics, university) or seed point-based region growing algorithm (see patents CN104036544B and CN 107944383 a).

However, both image-based and point cloud-based methods have significant problems. The image-based method usually loses information in the original point cloud, and particularly, the vertical structure information of ground objects reflected by the vertical distribution of the point cloud has a poor effect on extracting buildings covered by trees. In addition, the flat-topped vegetation which is manually trimmed is easily divided into roof peaks by mistake. The algorithm based on region growing needs to give a normal vector included angle threshold and a curvature threshold, and values of the two thresholds are influenced by data quality and roof flatness and are difficult to determine simply. More importantly, the roof points extracted by the method are not complete, and a large number of roof points which are missed to be detected often exist. The Hough transformation extraction plane needs three parameters, the Hough transformation of the three parameters is very complex, the three parameters need to be partitioned in a three-dimensional Hough space,used as a partition of an accumulator to vote. The granularity of spatial partitioning directly affects the effect of subsequent plane detection, and errors often exist in the detection result. The computational complexity is high in point cloud size, requiring a large amount of storage space. The RANSAC-based method has the basic idea that 3 points are randomly extracted from a point cloud to construct a plane model, and whether other points can be well fitted by using the plane is tested. The above process is repeatedly performed many times. Assuming that the point cloud has N points, randomly taking 3 points to have

A non-repeating combination. As N increases, the number of combinations increases rapidly. If N is up to ten thousand,

will be an astronomical number. The number of random samples is limited, which results in a significantly reduced probability of finding a robust model. Even if the random sampling times are increased to a limited extent, the problem is not solved significantly, but the calculation time is prolonged significantly. From this point of view, the RANSAC method is not suitable for the case of a large number of points. Furthermore, RANSAC theoretically may derive some models even for randomly distributed points. For example, in forest point clouds, RANSAC can also be used to detect some planes, and the points on the planes may be enough. In this case, it is difficult for the conventional RANSAC to ensure that the plane comes from a building rather than a non-building.

In order to overcome the defects of RANSAC, the RANSAC method is urgently needed to be improved, so that the RANSAC method is more suitable for carrying out plane extraction on massive airborne LiDAR point clouds, distinguishing forests and building roofs, and improving the extraction efficiency and accuracy.

Disclosure of Invention

The purpose of the invention is as follows: the invention aims to solve the problems of high calculation complexity and low execution efficiency when plane roof points of a building are extracted by using a RANSAC method in airborne LiDAR point cloud, and avoid classifying non-building points into building roof points so as to improve the efficiency and accuracy of roof point extraction of the building.

The technical scheme is as follows: the invention discloses a method for extracting building roof points from airborne laser radar point cloud, which comprises the following steps:

(1) acquiring airborne laser radar point cloud of a survey area; performing three-dimensional Euclidean clustering on point classes possibly comprising a building roof; calculating the average ground clearance of the point clouds in each cluster, processing each cluster with the average ground clearance larger than the height threshold according to the steps 2-4, and extracting the top point of the building house;

(2) generating a current cluster C_pCorresponding DSM, obtain C_pNumber of active pixels N in DSM_validCalculating C_pProjected area A on XY plane_C；

(3) If C is present_pProjected area A of_CGreater than or equal to the area threshold, calculating Laplacian L in DSM_pLess than a flatness threshold L_TOf the DSM active pixel total number N_validRatio P of_LaplacianIf said ratio is greater than a flat-top proportion threshold P_{Laplacian_T}Then cluster C is currently_pAll points in (a) are points on the roof of the building;

(4) if C is present_pProjected area A of_CIf the area is smaller than the area threshold, judging according to the following steps:

(4.1) translating the current cluster C_pTo obtain C_tmp，C_tmpThe center of gravity of (2) is at the origin of coordinates; initializing a current number of samples N_S0; initialization C_tmpIf all the points are not extracted, emptying the effective model sequence and the optimal model sequence;

(4.2) if the current sampling time N_SUp to the maximum number of samples N_{S_T}Skipping to step 4.5;

from C_tmpThe point p in which the random extraction state is not extracted₀ A 1 is to p₀Is modified to extracted, find C_tmpInner distance p₀Three nearest neighbourhood points p₁、p₂、p₃(ii) a If p is₁、p₂、p₃Distance p₀All the distances of (a) are sufficiently small, and four points can be fitted by using planes, and fitting is calculatedPlane P_fit(ii) a Otherwise, the current sampling times are N_SAdding one, and re-executing the step 4.2;

(4.3) calculation of C_tmpFrom each point to plane P_fitThe distance of (d); setting the distance less than a flatness threshold L_TThe points of (2) constitute an interior point set C_inlier(ii) a If set C_inlierThe number of the medium points is less than the threshold value of the minimum clustering points, and the current sampling times N_SAdding one, and re-executing the step 4.2; else for set C_inlierPerforming three-dimensional Euclidean clustering again;

(4.4) if only one cluster is obtained in step 4.3, the plane P is assigned_fitAdding an effective model sequence;

current number of samples N_SAdding one, and re-executing the step 4.2;

(4.5) traversing the effective model sequence, searching the effective model with the largest number of interior points as an optimal model, and adding the optimal model sequence; from C_tmpRemoving points in the internal point set corresponding to the optimal model, and initializing the current sampling times N if the remaining points are greater than the threshold of the minimum clustering points_s0; initialization C_tmpIf all the points are not extracted, emptying the effective model sequence and the optimal model sequence, and skipping to the step 4.2;

(4.6) calculation of C_pMinimum distance d between each point in the interior and the plane in the optimal model sequence_minD is mixing_minLess than a flatness threshold L_TThe points of (2) constitute a set C_{inlier_all}(ii) a If C is present_{inlier_all}Number of points C_pThe proportion of the number of the medium points is greater than the threshold value P_{inlier_all_T}Then C will be_{inlier_all}The point in (1) is divided into building roof points.

On the other hand, the invention also discloses a system for realizing the method for extracting the roof points of the building, which comprises the following steps:

the airborne laser radar point cloud clustering module is used for acquiring airborne laser radar point cloud of the survey area; clustering the acquired point cloud by adopting three-dimensional Euclidean clustering; calculating the average ground clearance of the point clouds in each cluster, and acquiring clusters with the average ground clearance larger than a height threshold; DSMA generating module for generating corresponding DSM for each cluster acquired by the airborne laser radar point cloud clustering module, acquiring the effective pixel number of the cluster in the DSM, and calculating the projection area A of the cluster in the XY plane_C(ii) a And according to A_CDividing the clusters into a large-area cluster and a facet cluster;

the large-area clustering classification module is used for judging whether points in the large-area accumulation classes are points on the roof of the building: the judging method comprises the following steps: calculating DSM for large area clustering, and counting Laplace operator L_pLess than a flatness threshold L_TTotal number of pixel occupied effective pixels N_validRatio P of_LaplacianIf said ratio is greater than a flat-top proportion threshold P_{Laplacian_T}Then cluster the current C_pAll the points in the building are divided into building roof points;

the small-area clustering classification module is used for extracting points on the roof of the building in the small-area clustering classification, and the extraction method comprises the following steps:

(4.1) translating the current small-area cluster C_pTo obtain C_tmp，C_tmpThe center of gravity of (2) is at the origin of coordinates; initializing a current number of samples N_S0; initialization C_tmpIf all the points are not extracted, emptying the effective model sequence and the optimal model sequence;

(4.2) if the current sampling time N_sUp to the maximum number of samples N_{s_T}Skipping to step 4.5;

from C_tmpThe point p in which the random extraction state is not extracted₀ A 1 is to p₀Is modified to extracted, find C_tmpInner distance p₀Three nearest neighbourhood points p₁、p₂、p₃(ii) a If p is₁、p₂、p₃Distance p₀Are sufficiently small, and four points can be fitted with a plane to calculate a fitted plane P_fit(ii) a Otherwise, the current sampling times are N_SAdding one, and re-executing the step 4.2;

(4.3) calculation of C_tmpFrom each point to plane P_fitThe distance of (d); setting the distance less than a flatness threshold L_TThe points of (2) constitute an interior point set C_inlier(ii) a If set C_inlierThe number of the medium points is less than the threshold value of the minimum clustering points, and the current sampling times N_SAdding one, and re-executing the step 4.2; else for set C_inlierPerforming three-dimensional Euclidean clustering again;

(4.4) if only one cluster is obtained in step 4.3, the plane P is assigned_fitAdding an effective model sequence;

current number of samples N_SAdding one, and re-executing the step 4.2;

(4.5) traversing the effective model sequence, searching the effective model with the largest number of interior points as an optimal model, and adding the optimal model sequence; from C_tmpRemoving points in the internal point set corresponding to the optimal model, and initializing the current sampling times N if the remaining points are greater than the threshold of the minimum clustering points_S0; initialization C_tmpIf all the points are not extracted, emptying the effective model sequence and the optimal model sequence, and skipping to the step 4.2;

(4.6) calculation of C_pMinimum distance d between each point in the interior and the plane in the optimal model sequence_minD is mixing_minLess than a flatness threshold L_TThe points of (2) constitute a set C_{inlier_all}(ii) a If C is present_{inlier_all}Number of points C_pThe proportion of the number of the medium points is greater than the threshold value P_{inlier_all_T}Then C is_{inlier_all}A point in (1) is a point on the roof of a building.

Has the advantages that: compared with the prior art, the method for extracting the building roof vertex from the airborne laser radar point cloud disclosed by the invention has the following advantages: 1. the invention adopts three-dimensional Euclidean clustering, can effectively divide the roof of a building and adjacent trees into different clusters under most conditions, and can eliminate most of non-roof tree points by carrying out Euclidean clustering on the points in the plane model after the points in the plane model are extracted, thereby improving the accuracy rate of roof point extraction; (ii) a 2. The method adopts a strategy of dividing and treating the large roof and the small roof, avoids high calculation amount and potential missing detection possibility caused by plane detection of the large roof by using RANSAC, and maintains high accuracy; 3. for the facet accumulation class, using improvementsThe RANSAC algorithm of (1); given the number of points in a cluster as N, the number of potential planar models is set from

Down to N. If N is less than the specified operation times, the random point taking actually only needs to be operated for N times; 4. the method is simple and feasible in theory and high in calculation efficiency. For 300 ten thousand points, this can typically be done in ten seconds. Multiple runs may be performed to more fully extract the roof vertices.

Drawings

FIG. 1 is a flow chart of a method for extracting a roof point of a building from an airborne laser radar point cloud according to the present disclosure;

FIG. 2 is a schematic diagram of Euclidean clustering;

FIG. 3 is a schematic diagram of construction of a DSM;

FIG. 4 is a diagram illustrating a judgment of a randomly selected point p₀Whether or not to contact 3 surrounding neighborhood points p₁、p₂、p₃A schematic diagram that can be approximately fitted with a plane;

FIG. 5 is a schematic diagram showing the components of a system for extracting the building roof vertex from the airborne lidar point cloud according to the present invention;

FIG. 6 is a diagram illustrating the effect of the present invention.

Detailed Description

The invention is further elucidated with reference to the drawings and the detailed description.

As shown in fig. 1, the invention discloses a method for extracting building roof vertex from airborne laser radar point cloud, comprising the following steps:

step 1, acquiring airborne laser radar point cloud of a survey area; clustering the acquired point cloud by adopting three-dimensional Euclidean clustering; calculating the average ground clearance of the point clouds in each cluster, processing each cluster with the average ground clearance larger than the height threshold according to the steps 2-4, and extracting the top point of the building house;

as shown in fig. 2, the result of euclidean clustering is to make the minimum distance d between different clusters equal to or greater than the search radius R. Parameters of the euclidean clustering algorithm include: search radius R, clustering minimumNumber of points threshold N_minClustering maximum point threshold N_maxThese parameter values affect the outcome of the clustering. In the present invention, the above parameters are set as follows:

the search radius R ranges from (R)_min,R_max) R minimum value R_minInversely related to the laser pulse foot point density, in this embodiment,

d is the density of laser pulse foot points, which can be obtained from the aviation design. If R is less than R_minThis may result in too few clusters, and a roof may be divided into multiple clusters. R maximum value R_maxThe minimum distance between adjacent buildings may be taken, such as 1.5m, 2.0m, etc. If R exceeds R_maxA plurality of roofs, trees, and the like may be gathered together, and accurate division cannot be achieved.

N_min、N_maxPositively correlated to the minimum and maximum area of the building roof. Setting the minimum area of the roof of a building as A_minMaximum area of A_maxThen N is_min、N_maxThe calculation formula of (A) is as follows:

N_min＝A_minD-Q

N_max＝A_maxD+Q

where Q is an error adjustment constant, in this embodiment Q is 50. Q is used to accommodate A_minD、A_maxD and the number of real minimum and maximum points of the roof.

Calculating the average ground clearance of the point clouds in the cluster, comprising:

(1.1) acquiring the ground height of each point in the cluster;

the ground height of one point can be obtained by constructing a two-dimensional index on the existing known ground points and carrying out nearest neighbor query, and can also be obtained by querying DEM data of other sources. If the known ground point calculation is used, a two-dimensional index T needs to be constructed for the known ground point by utilizing the XY coordinates of the known ground point_grd。T_grdEither KD or octree indexing. Setting clusters to be computedThe coordinate of the nth point is (X)_n,Y_n,H_n) By using T_grdQuery for distance (X) in XY plane_n,Y_n) The number of the nearest ground points G, and then obtaining the three-dimensional coordinates of G, and taking the elevation of G as the ground height H corresponding to the nth point in the cluster_Tn；

(1.2) calculating the height from the ground of each point in the cluster: △ H_n＝H_n-H_Tn；

Average ground clearance of point clouds within cluster:

n is the total number of points within the cluster.

Since a building roof is generally at least a certain height above the ground, e.g. 2.5m, if Δ H_meanLess than a given height threshold Δ H_{mean_T}Then the cluster is deemed not to represent a building if Δ H_meanIs greater than or equal to Δ H_{mean_T}Processing according to the step 2-4 to determine points belonging to the building roof;

in this example,. DELTA.H_{mean_T}The value is 2.5 m.

Step 2, generating a current cluster C_pCorresponding DSM, obtain C_pNumber of active pixels N in DSM_validCalculating C_pArea A in XY plane_c(ii) a The method comprises the following specific steps:

(2.1) calculation of C_pMinimum X coordinate of midpoint X_minMaximum X coordinate X_maxMinimum Y coordinate Y_minMaximum Y coordinate Y_max(ii) a Will be expressed by (X)_min,Y_min)、(X_max,Y_max) Rectangular region partitioning into N for diagonal vertices_rowLine, N_colA grid of column, grid size R, wherein:

N_row＝int(Y_max-Y_min)/R+1

N_col＝int(X_max-X_min)/R+1

r is the search radius of the three-dimensional Euclidean cluster, and int (·) is an integer function;

(2.2) adding C_pThe point in (C) is assigned to the pixel in the grid according to its X, Y coordinate, resulting in the current cluster C_pA corresponding DSM;

the nth point is mapped into pixel (i, j),

wherein (X)_n,Y_n,H_n) Is the coordinate of the nth point;

(2.3) traversing the pixels in the grid; if a pixel has a point, the pixel is identified as an effective pixel, the height index of the point is taken as the value of DSM, and the DSM value of the pixel without the point is set as H_invalid(ii) a As shown in fig. 3, a white pixel is a pixel where no point falls, and a gray pixel is a pixel where a point falls; the height indexes are as follows: the minimum elevation value of the falling point, or the maximum elevation value of the falling point, or the average elevation value of the falling point;

(2.4) traversing DSM pixel by pixel, and counting the number N of effective pixels_valid；C_pProjected area A on XY plane_CComprises the following steps: a. the_c＝N_validR²。

If A is_CGreater than or equal to a specified area threshold A_{C_T}Then, consider C_pClustering in large area; otherwise C_pClustering for small area.

Step 3, calculating Laplacian L in DSM for large-area clustering_pLess than a flatness threshold L_TTotal number of pixel occupied effective pixels N_validRatio P of_LaplacianIf said ratio is greater than a flat-top proportion threshold P_{Laplacian_T}Then cluster C is currently_pAll points in (a) are points on the roof of the building;

laplacian L of DSM at pixel (i, j)_p(i, j) is:

flatness threshold L_TOf values and point cloudsThe relative error is positively correlated, and in this example, 0.10-0.20m is taken. Let Laplace operator L_pValue less than L_THas N pixels_planeAnd then:

if P is_LaplacianGreater than or equal to a given flat top ratio threshold P_{Laplacian_T}Then, consider C_pIs a building consisting of flat tops, and is characterized in that_pThe interior points are all divided into building points. Roof with larger area and simpler structure, P_{Laplacian_T}The larger the value can be, the value is 0.6 in this embodiment.

Step 4, if clustering the area A of the projection area_CJudging C by adopting an improved RANSAC algorithm when the area is smaller than the area threshold value_pThe midpoint classification, the steps are as follows:

(4.1) translating the current cluster C_pTo obtain C_tmpLet C be_tmpThe center of gravity of (2) is at the origin of coordinates;

compute cluster C_pCenter of gravity O of_p(X_Op,Y_Op,H_Op)：

C_tmpThe coordinates of the nth point in (1) are: (X)_n-X_Op,Y_n-Y_Op,H_n-H_Op) Then C is_tmpThe center of gravity of (2) is at the origin of coordinates;

initializing a current number of samples N_S0; initialization C_tmpIf all the points are not extracted, emptying the effective model sequence and the optimal model sequence;

(4.2) if the current sampling time N_SUp to the maximum number of samples N_{S_T}Skipping to step 4.5; maximum number of samples N_{s_T}The value is a preset value, the value of the preset value does not exceed the number N of the points in the current cluster at most, and the value in the embodiment is 100;

from C_tmpThe point p in which the random extraction state is not extracted₀A 1 is to p₀Is modified to extracted, find C_tmpInner distance p₀Three nearest neighbourhood points p₁、p₂、p₃(ii) a If p is₁、p₂、p₃Distance p₀Are sufficiently small, and four points can be fitted with a plane to calculate a fitted plane P_fit(ii) a Otherwise, the current sampling times are N_SAdding one, and re-executing the step 4.2;

in this embodiment, p is judged₁、p₂、p₃And p₀Whether the Euclidean distances are all smaller than R or not is judged, if not, the distances are not small enough, the neighborhood relationship is not stable enough, and the four points cannot be used for stably estimating the local plane; if yes, judging whether the four points can be approximately fitted by using a plane, wherein the judgment method can adopt least square fitting to form the plane and then judge p₁、p₂、p₃Whether the maximum distance to the plane is small enough or not can be judged based on the plane included angle. In view of the simplicity and high efficiency of the method based on the planar angle, the method will be described with reference to fig. 4.

p₁、p₂、p₃And p₀Forming three triangles, i.e. Δ p₀p₁p₂、Δp₀p₂p₃、Δp₀p₃p₁. The normal vectors, Δ p, of these triangles are calculated₀p₁p₂Normal vector N of₀₁₂Comprises the following steps: n is a radical of₀₁₂＝(p₁-p₀)×(p₂-p₀) (ii) a To N₀₁₂Modular normalization is performed. Δ p₀p₂p₃Normal vector N of₀₂₃、Δp₀p₃p₁Normal vector N of₀₃₁Calculations and normalization were performed using a similar method.

Calculating N₀₁₂And N₀₂₃Angle of (theta)₀：θ₀＝cos^-1(N₀₁₂·N₀₂₃)；

N₀₁₂And N₀₃₁Angle of (theta)₁、N₀₂₃And N₀₃₁Angle of (theta)₂Calculations were performed using a similar method. cos (chemical oxygen demand)^-1Is an inverse cosine function.

If theta is greater than theta₀、θ₁、θ₂Are both close to 0 or 180, consider p₀、p₁、p₂And p₃Can use a certain plane P_fitFitting is carried out; otherwise, considering that the four points can not be fitted by using a plane;

calculating P_fitIs given by P_fitThe equation of (a) is: ax + by + ch + d is 0;

a. b and a are components of the plane normal vector on three coordinate axes respectively, and d is ═ ax₀+by₀+ch₀)，(x₀,y₀,h₀) Is p₀The coordinates of (a).

(4.3) calculation of C_tmpFrom each point to plane P_fitWherein the nth point is to P_fitDistance d of_fitnComprises the following steps:

(x_tmpn,y_tmpn,h_tmpn) Is C_tmpCoordinates of the nth point;

setting the distance less than a flatness threshold L_TThe points of (2) constitute an interior point set C_inlier(ii) a If set C_inlierThe number of medium points is less than the threshold value N of the minimum clustering points_minCurrent number of samples N_sAdding one, and re-executing the step 4.2; else for set C_inlierThree-dimensional euclidean clustering is performed again, where the parameters of euclidean clustering are set as follows:

the search radius is R, the threshold value of the minimum clustering points is 1, and the threshold value of the maximum clustering points is a set C_inlierTotal number of points N in_inlier；

(4.4) if only one cluster is obtained in step 4.3, the plane P_fitAdding the effective model into an effective model sequence as an effective model;

current number of samples N_sAdding one, and re-executing the step 4.2;

(4.5) traversing the effective model sequence, searching the effective model with the largest number of interior points as an optimal model, and adding the optimal model sequence; if the model with the largest number of interior points has a plurality of interior points, respectively calculating a set C_inlierTaking the model with the minimum standard deviation as the optimal model, and adding the optimal model sequence;

from C_tmpRemoving points in the internal point set corresponding to the optimal model, if the residual points are more than the clustering minimum point threshold value N_minThen initialize the current sampling number N_s0; initialization C_tmpIf all the points are not extracted, emptying the effective model sequence and the optimal model sequence, and skipping to the step 4.2;

(4.6) calculation of C_pMinimum distance d between each point in the interior and the plane in the optimal model sequence_minD is mixing_minLess than a flatness threshold L_TThe points of (2) constitute a set C_{inlier_all}(ii) a If C is present_{inlier_all}Number of points C_pThe proportion of the number of the medium points is greater than the threshold value P_{inlier_all_T}Then, consider C_{inlier_all}The building roof is divided into points in the building roof into building roof points. In this example, P_{inlier_all_T}The value of (A) is 0.50.

And (4) skipping to the step (2) to judge and process the next cluster until all the clusters divided in the step (1) are processed.

The embodiment also discloses a system for implementing the method for extracting the roof point of the building from the point cloud of the airborne laser radar, as shown in fig. 5, the method includes:

the airborne laser radar point cloud clustering module is used for acquiring airborne laser radar point cloud of the survey area according to the step 1; clustering the acquired point cloud by adopting three-dimensional Euclidean clustering; calculating the average ground clearance of the point clouds in each cluster, and acquiring clusters with the average ground clearance larger than a height threshold;

the DSM generation module is used for generating a corresponding DSM for each cluster acquired by the airborne laser radar point cloud clustering module according to the step 2, acquiring the number of effective pixels clustered in the DSM and calculating the projection area of the cluster in an XY plane;

the large-area clustering classification module is used for judging whether the points in the large-area accumulation classes are the top points of the building houses or not according to the step 3;

and the small-area cluster classification module is used for judging whether the points in the small-area cluster are the top points of the building house or not according to the step 4.

FIG. 6 is a diagram illustrating the effect of the present invention, wherein FIG. 6- (a) is a high resolution optical image of a place, FIG. 6- (b) is input point cloud data, wherein the point cloud has been preliminarily classified, and FIG. 6- (c) is the classified point cloud data. Comparing fig. 6- (a) and 6- (c), it can be seen that the method of the present invention has good extraction effect for large and small buildings, flat-topped and non-flat-topped buildings, and building roof top with or without tree shelter.

The invention provides a method for extracting building roof vertexes from airborne LiDAR point clouds based on a divide-and-conquer strategy, and particularly relates to a method for quickly extracting planes from small-area clustering point clouds by using an improved RANSAC method and judging whether a model is effective or not by using clustering numbers for model interior points. There are numerous ways of implementing the technical solution and the above are only preferred embodiments of the present invention. It should be noted that numerous modifications and adaptations can be made by those skilled in the art without departing from the principles of the present invention. Such modifications and refinements are also to be considered within the scope of the present invention. All the components not specified in the present embodiment can be realized by the prior art.

Claims (10)

1. The method for extracting the building roof from the airborne laser radar point cloud is characterized by comprising the following steps of:

(1) acquiring airborne laser radar point cloud of a survey area; clustering the acquired point cloud by adopting three-dimensional Euclidean clustering; calculating the average ground clearance of the point clouds in each cluster, processing each cluster with the average ground clearance larger than the height threshold according to the steps (2) to (4), and extracting the top point of the building house;

(2) generating a current cluster C_pCorresponding DSM, obtain C_pNumber of active pixels N in DSM_validCalculating C_pProjected area A on XY plane_C；

(3) If cluster C_pProjected area A_CGreater than or equal to the area threshold, calculating Laplacian L in DSM_pLess than a flatness threshold L_TTotal number of pixel occupied effective pixels N_validRatio P of_LaplacianIf said ratio is greater than a flat-top proportion threshold P_{Laplacian_T}Then cluster the current C_pAll the points in the building are divided into building roof points;

(4) if clustering projected area A_CIf the area is smaller than the area threshold, the judgment is carried out according to the following steps:

(4.1) translating the current cluster C_pTo obtain C_tmp，C_tmpThe center of gravity of (2) is at the origin of coordinates; initializing a current number of samples N_S0; initialization C_tmpIf all the points are not extracted, emptying the effective model sequence and the optimal model sequence;

(4.2) if the current sampling time N_SUp to the maximum number of samples N_{S_T}Jumping to the step (4.5);

from C_tmpThe point p in which the random extraction state is not extracted₀A 1 is to p₀Is modified to extracted, find C_tmpInner distance p₀Three nearest neighbourhood points p₁、p₂、p₃(ii) a If p is₁、p₂、p₃Distance p₀Is sufficiently small, and four points can be approximated by plane fitting to calculate a fitting plane P_fit(ii) a Otherwise, the current sampling times are N_SAdding one, and re-executing the step (4.2);

(4.3) calculation of C_tmpFrom each point to plane P_fitThe distance of (d); setting the distance less than a flatness threshold L_TThe points of (2) constitute an interior point set C_inlier(ii) a If set C_inlierThe number of the medium points is less than the threshold value of the minimum clustering points, and the current sampling times N_SAdding one, and re-executing the step (4.2); else for set C_inlierPerforming three-dimensional Euclidean clustering again;

(4.4) if only one cluster is obtained in step (4.3), the plane P is divided into two_fitAdding an effective model sequence;

current number of samples N_SAdding one, and re-executing the step (4.2);

(4.5) traversing the effective model sequence, searching the effective model with the largest number of interior points as an optimal model, and adding the optimal model sequence; from C_tmpRemoving points in the internal point set corresponding to the optimal model, and initializing the current sampling times N if the remaining points are greater than the threshold of the minimum clustering points_S0; initialization C_tmpIf all the points are not extracted, emptying the effective model sequence and the optimal model sequence, and skipping to the step (4.2);

(4.6) calculation of C_pMinimum distance d between each point in the interior and the plane in the optimal model sequence_minD is mixing_minLess than a flatness threshold L_TThe points of (2) constitute a set C_{inlier_all}(ii) a If C is present_{inlier_all}Number of points C_pThe proportion of the number of the medium points is greater than the threshold value P_{inlier_all_T}Then C will be_{inlier_all}The point in (1) is divided into building roof points.

2. The method for extracting the apex of the building according to claim 1, wherein the parameters of the three-dimensional euclidean cluster in the step (1) are set as follows:

the search radius R ranges from (R)_min，R_max) Wherein

D is the density of the laser pulse foot spots, R_maxIs the minimum distance between adjacent buildings;

clustering minimum point threshold N_min＝A_minD-Q

Threshold value N of maximum points of cluster_max＝A_maxD+Q

Wherein Q is an error regulation constant, A_minFor the minimum area of the roof of a building, A_maxThe largest area of the roof of the building.

3. The method for extracting the apex of the building according to claim 1, wherein the step (1) of calculating the average height from the ground of the point cloud in the cluster comprises:

obtaining the ground height of each point in the cluster: the nth point has the coordinate of (X)_n，Y_n，H_n) Corresponding to ground height H_TnIs the distance (X) on the XY plane_n，Y_n) Elevation of the nearest ground point;

calculate the terrain clearance for each point in the cluster: Δ H_n＝H_n-HT_n；

Average ground clearance of point clouds within cluster:

n is the total number of points within the cluster.

4. The method for extracting the apex of the building house according to claim 1, wherein the step (2) specifically comprises:

(2.1) calculation of C_pMinimum X coordinate of midpoint X_minMaximum X coordinate X_maxMinimum Y coordinate Y_minMaximum Y coordinate Y_max(ii) a Will be expressed by (X)_min，Y_min)、(X_max，Y_max) Rectangular region partitioning into N for diagonal vertices_rowLine, N_colA grid of column, grid size R, wherein:

N_row＝int(Y_max-Y_min)/R+1

N_col＝int(X_max-X_min)/R+1

r is the search radius of the three-dimensional Euclidean cluster, and int (·) is an integer function;

(2.2) adding C_pThe points in (C) are assigned to pixels in the grid to obtain the current cluster C_pA corresponding DSM;

the nth point is mapped into pixel (i, j),

wherein (X)_n，Y_n，H_n) Is the coordinate of the nth point;

(2.3) traversing the pixels in the grid; if a pixel has a point, the height index of the point is taken as the value of DSM, the pixel is an effective pixel, and the DSM value of the pixel without the point is set as H_invalid(ii) a The height indexes are as follows: the minimum elevation value of the falling point, or the maximum elevation value of the falling point, or the average elevation value of the falling point;

(2.4) traversing DSM pixel by pixel, and counting the number N of effective pixels_valid；C_pProjected area A on XY plane_CComprises the following steps: a. the_c＝N_validR²。

5. The method of claim 4, wherein the step (3) DSM extracts the Laplacian L at pixel (i, j)_p(i, j) is:

6. system for follow extraction building roof summit in the airborne laser radar point cloud, its characterized in that includes:

the airborne laser radar point cloud clustering module is used for acquiring airborne laser radar point cloud of the survey area; clustering the acquired point cloud by adopting three-dimensional Euclidean clustering; calculating the average ground clearance of the point clouds in each cluster, and acquiring clusters with the average ground clearance larger than a height threshold;

a DSM generation module for generating corresponding DSM for each cluster acquired by the airborne laser radar point cloud clustering module and acquiring the DSMThe effective pixel number of the cluster in DSM, calculating the projection area A of the cluster in XY plane_C(ii) a And according to A_CDividing the clusters into a large-area cluster and a facet cluster;

the large-area clustering classification module is used for judging whether points in the large-area accumulation classes are points on the roof of the building: the judging method comprises the following steps: for large area clustered DSM, statistical Laplacian L_pLess than a flatness threshold L_TOf the DSM active pixel total number N_validRatio P of_LaplacianIf said ratio is greater than a flat-top proportion threshold P_{Laplacian_T}Then cluster the current C_pAll the points in the building are divided into building roof points;

the small-area clustering classification module is used for extracting points on the roof of the building in the small-area clustering classification, and the extraction method comprises the following steps:

(4.1) translating the current small-area cluster C_pTo obtain C_tmpLet C be_tmpThe center of gravity of (2) is at the origin of coordinates; initializing a current number of samples N_S0; initialization C_tmpIf all the points are not extracted, emptying the effective model sequence and the optimal model sequence;

(4.2) if the current sampling time N_SUp to the maximum number of samples N_{S_T}Jumping to the step (4.5);

from C_tmpThe point p in which the random extraction state is not extracted₀A 1 is to p₀Is modified to extracted, find C_tmpInner distance p₀Three nearest neighbourhood points p₁、p₂、p₃(ii) a If p is₁、p₂、p₃Distance p₀Is sufficiently small, and four points can be approximated by plane fitting to calculate a fitting plane P_fit(ii) a Otherwise, the current sampling times are N_SAdding one, and re-executing the step (4.2);

(4.3) calculation of C_tmpFrom each point to plane P_fitThe distance of (d); setting the distance less than a flatness threshold L_TThe points of (2) constitute an interior point set C_inlier(ii) a If set C_inlierThe number of the medium points is less than the threshold value of the minimum clustering points, and the current sampling times N_SAdding one, and re-executing the step (4.2); else for set C_inlierPerforming three-dimensional Euclidean clustering again;

(4.4) if only one cluster is obtained in step (4.3), the plane P is divided into two_fitAdding an effective model sequence;

current number of samples N_SAdding one, and re-executing the step (4.2);

(4.5) traversing the effective model sequence, searching the effective model with the largest number of interior points as an optimal model, and adding the optimal model sequence; from C_tmpRemoving points in the internal point set corresponding to the optimal model, and initializing the current sampling times N if the remaining points are greater than the threshold of the minimum clustering points_S0; initialization C_tmpIf all the points are not extracted, emptying the effective model sequence and the optimal model sequence, and skipping to the step (4.2);

(4.6) calculation of C_pMinimum distance d between each point in the interior and the plane in the optimal model sequence_minD is mixing_minLess than a flatness threshold L_TThe points of (2) constitute a set C_{inlier_all}(ii) a If C is present_{inlier_all}Number of points C_pThe proportion of the number of the medium points is greater than the threshold value P_{inlier_all_T}Then C is_{inlier_all}A point in (1) is a point on the roof of a building.

7. The system for extracting the vertices of the building house from the airborne lidar point cloud of claim 6, wherein the parameters of the three-dimensional euclidean clusters in the airborne lidar point cloud clustering module and the facet accumulation class classification module are set as:

the search radius R ranges from (R)_min，R_max) Wherein

D is the density of the laser pulse foot spots, R_maxIs the minimum distance between adjacent buildings;

clustering minimum point threshold N_min＝A_minD-Q

Threshold value N of maximum points of cluster_max＝A_maxD+Q

Wherein Q is an error regulation constant, A_minFor the minimum area of the roof of a building, A_maxThe largest area of the roof of the building.

8. The system of claim 6, wherein the calculating the average ground clearance of the point clouds within the clusters in the airborne lidar point cloud clustering module comprises:

obtaining the ground height of each point in the cluster: the nth point has the coordinate of (X)_n，Y_n，H_n) Height of ground H_TnIs the distance (X) on the XY plane_n，Y_n) Elevation of the nearest ground point;

calculate the terrain clearance for each point in the cluster: Δ H_n＝H_n-H_Tn；

Average ground clearance of point clouds within cluster:

n is the total number of points within the cluster.

9. The system for extracting building roof points from an onboard lidar point cloud of claim 6, wherein the DSM generation module generates cluster C_pCorresponding DSM, get C_pEffective pixel number in DSM and calculate cluster C_pThe steps of projecting the area on the XY plane and dividing the clusters into large area clusters and small area clusters include:

(2.1) calculation of C_pMinimum X coordinate of midpoint X_minMaximum X coordinate X_maxMinimum Y coordinate Y_minMaximum Y coordinate Y_max(ii) a Will be expressed by (X)_min，Y_min)、(X_max，Y_max) Rectangular region partitioning into N for diagonal vertices_rowLine, N_colA grid of column, grid size R, wherein:

N_row＝int(Y_max-Y_min)/R+1

N_col＝int(X_max-X_min)/R+1

r is the search radius of the three-dimensional Euclidean cluster, and int (·) is an integer function;

(2.2) adding C_pThe points in (C) are assigned to pixels in the grid to obtain the current cluster C_pA corresponding DSM;

the nth point is mapped into pixel (i, j),

wherein (X)_n，Y_n，H_n) Is the coordinate of the nth point;

(2.3) traversing the pixels in the grid; if a pixel has a point, the pixel is identified as an effective pixel, the height index of the point is taken as the value of DSM, and the DSM value of the pixel without the point is set as H_invalid(ii) a The height indexes are as follows: the minimum elevation value of the falling point, or the maximum elevation value of the falling point, or the average elevation value of the falling point;

(2.4) traversing DSM pixel by pixel, and counting the number N of effective pixels_valid；C_pProjected area A on XY plane_CComprises the following steps: a. the_c＝N_validR²；

If clustering projected area A_CIf the area is larger than or equal to the area threshold value, clustering C_pClustering in large area;

if clustering projected area A_CLess than the area threshold, cluster C_pClustering for small area.

10. The system for extracting building roof points from an onboard lidar point cloud of claim 6, wherein the large area cluster classification module calculates the Laplacian L of DSM at pixel (i, j)_p(i, j) is:

Images