Method for automatically identifying building roof structure based on airborne point cloud

Abstract

A method for automatically identifying a building roof structure based on airborne point cloud relates to the field of surveying and mapping, and is characterized by comprising the following steps: laser point cloud preprocessing, building house extraction, marshalling building roof patch laser point cloud segmentation, roof point cloud fitting plane equation set, and roof structure type generation by a plurality of fitting planes. And verifying the house structure type. According to the method, the airborne radar is adopted, the high-precision three-dimensional point cloud of the urban building roof patch can be rapidly obtained in a large area, the method has the advantages of being high in precision and high in penetrability, meanwhile, the building roof patch can be segmented and a roof plane equation can be fitted by adopting an improved algorithm, and the method is high in practicability. The method is a high-precision and high-efficiency data acquisition means for quickly acquiring the urban three-dimensional model structure type.

Description

Technical Field

[0001] The present invention relates to the field of surveying and mapping. Specifically, an airborne radar can be used to quickly and widely acquire high-precision three-dimensional point clouds of rooftop patches of urban buildings. The method has the advantages of high precision and high penetration. An improved algorithm can be used to segment building rooftop patches and fit roof plane equations. The method is a high-precision and high-efficiency data acquisition method for quickly acquiring structural types of urban three-dimensional models. The method is a method for automatically identifying building rooftop structures based on airborne point clouds. Background Art

[0002] In today's era, technology is developing rapidly and urban development needs are becoming increasingly diverse. Against this dual background, the real-life 3D city model (LOD1.3 level) came into being.

[0003] In terms of the background of technological development, surveying and mapping technology has made significant progress. Advanced laser scanning technology can accurately obtain three-dimensional information of various objects in the city, and photogrammetry technology can obtain rich image data from different angles, all of which provide a solid data foundation for building high-precision three-dimensional urban models. The rapid improvement of computer hardware performance has made large-scale data processing more efficient and faster. In the process of building a three-dimensional urban model, massive point cloud data and image data can be quickly processed and analyzed, creating conditions for the rapid construction and optimization of the model. At the same time, the accuracy of positioning technologies such as GNSS continues to improve, providing reliable guarantees for the accurate spatial positioning of various elements in the three-dimensional urban model, ensuring the accuracy and reliability of the model in terms of geographical location.

[0004] In terms of the background of urban development needs, with the rapid expansion and continuous evolution of cities, urban planning and management urgently need more intuitive and accurate three-dimensional models. The urban three-dimensional model (LOD1.3 level) can clearly display the city's building distribution, spatial pattern and other information, helping government departments to better grasp the context of urban development, formulate scientific and reasonable planning schemes, and effectively improve the utilization efficiency of urban space resources.

[0005] Based on technological development and urban development, the country is accelerating the construction and improvement of a new basic surveying and mapping system, and establishing a new business structure with modern geodetic datum, real-life 3D, and spatiotemporal big data platforms as the main service systems. The main products of new basic surveying and mapping are transforming from basic scale topographic map data to real-life 3D data. On February 24, 2022, the General Office of the Ministry of Natural Resources issued the "Notice on Comprehensively Promoting the Construction of Real-life 3D China" (hereinafter referred to as the "Notice"), which clarified the goals, tasks and division of labor for the construction of real-life 3D China. The "Notice" put forward two major construction goals for the construction of real-life 3D China. In the same year, the "Work Plan for Rapid Construction of Urban 3D Models (Draft for Comments)" compiled by the Department of Land Surveying and Mapping of the Ministry of Natural Resources proposed: By the end of 2024, the latest technical methods and data sources will be used to achieve the coverage of urban development boundaries of cities at or above the prefecture level across the country by urban 3D models (LOD1.3 level), and provide a unified 3D spatial positioning framework and analysis basis for economic and social development and natural resource management. Efficient and fast model construction, accurate model attribute information and exquisite texture feature presentation are key steps in the development of 3D city scene construction technology. Among them, the building structure type is a required item of attribute information, which is related to the completeness of the city 3D model (LOD1.3 level) results.

[0006] In summary, how to quickly and accurately obtain the building structure information of the urban three-dimensional model has become an urgent problem to be solved. At present, the mainstream idea is to directly obtain the building structure type based on large-scale topographic map data. This method requires basic data support and has the problem of not updating in time. There is no method based on airborne laser radar point cloud data to automatically obtain the building roof structure type information through intelligent algorithms to improve the automation and intelligence level of house structure type extraction. Summary of the invention

[0007] In order to overcome the drawbacks of existing technical solutions, airborne radar can be used to quickly obtain high-precision three-dimensional point clouds of urban building roof surfaces over a large area. It has the advantages of high precision and high penetration. At the same time, the improved algorithm can segment the building roof surfaces and fit the roof plane equation. It is a high-precision and efficient data acquisition method for quickly obtaining the structural type of urban three-dimensional models.

[0008] A method for automatically identifying a building roof structure based on an airborne point cloud, the method comprising the following steps: laser point cloud preprocessing, building house extraction, segmentation and grouping of building roof surface laser point clouds, roof point cloud fitting plane equation set, multiple fitting planes generating roof structure type, house structure type verification, wherein:

[0009] Laser house structure type verification: classify the building roof structure type, record it in the corresponding building outline attribute, obtain the ridge line based on the segmented and grouped roof patch point cloud and building outline, predict the structure type based on geometric features and ridge line, judge the predicted structure type and compare it with the structure type in the corresponding building outline attribute, and modify the inconsistent structure type;

[0010] Point cloud preprocessing: The original laser point cloud data is obtained through the mutual coordination and combined operation between the laser transmitter, inertial navigation and GNSS equipment. The laser point cloud is matched with the point cloud of a single flight and a single flight strip to obtain the registration point cloud of each flight. Then, the rigid registration point cloud of each flight is regionally optimized to obtain the optimized point cloud of the entire survey area. The optimized point cloud of the survey area is then converted from geodetic coordinates to the standard national 2000 coordinate system and the local city coordinate system to obtain the target point cloud of the survey area. The target point cloud of the survey area is converted into point cloud data in an internationally accepted format.

[0011] Building extraction: remove noise points from the preprocessed point cloud, extract the ground point cloud, calculate the point cloud distance and normal vector, filter the low-value point cloud according to the distance and normal vector factors, delete the clustered tree point cloud around the building, and obtain the building roof point cloud data;

[0012] Segmentation and grouping of building roof laser point cloud: Randomly select initial seed points and basic points from the building roof point cloud to construct a plane by normal vector and distance determination. For the points that cannot construct a plane, use them as other points to determine the distance and normal vector between the plane and determine the inner points and outer points. After grouping the inner points, perform a cycle of inner point selection on different planes and group them in order.

[0013] Roof point cloud fitting plane equation set: Based on the outline of a single building, extract the laser point cloud of a single roof surface after segmentation and grouping. According to the grouping number of the laser point cloud of a single roof surface, select the laser point cloud of the roof surface in the same group as the plane fitting object, perform feature centering on the fitting object, establish the covariance matrix to obtain the eigenvalue and its corresponding eigenvector, and establish the plane equation set and normal vector set for the single fitting object after traversal processing;

[0014] Generate roof structure type from multiple fitting planes: Take the building outline as the range, obtain the polygonal area of ​​the building roof patch, remove the plane equations and corresponding normal vectors with an area less than the threshold, obtain a new set of plane equations and normal vectors, and determine the number of normal vectors, perform normal vector normalization, area weight analysis, and normal vector comparison analysis to generate the house structure type;

[0015] Verification of building structure type: classify the building roof structure type and record it in the corresponding building outline attribute. Use the segmented and grouped roof patch point cloud and building outline as the data basis to obtain the ridge line. Predict the structure type based on the geometric features and ridge line. Determine the predicted structure type and compare it with the structure type in the corresponding building outline attribute. Modify the inconsistent structure type.

[0016] A method for automatically identifying a building roof structure based on airborne point clouds, wherein the specific steps of extracting the building roof structure are:

[0017] Eliminate noise points: Remove extremely high and extremely low discrete points and noise points from the general format point cloud data to obtain a three-dimensional laser point cloud including topography, buildings and surrounding small pieces, trees, power lines, etc. within the target measurement area;

[0018] Extract ground point cloud: Use opening and closing operations on the three-dimensional laser points within the target measurement area to preliminarily determine the ground points, remove isolated high points and fill small depressions to obtain a relatively smooth surface model, compare the point cloud with the constructed surface model, and reclassify the ground points and non-ground points; points that are close to the surface model are determined as ground points; points that are far away are determined as non-ground points. Update the surface model based on the new classification results, and repeat the above process until the classification results are stable or a certain number of iterations are reached, and obtain a three-dimensional laser point cloud containing ground points and non-ground points;

[0019] Calculate the point cloud distance and normal vector: Calculate the distance from each point to the ground point directly below and assign it to the point cloud attribute; the calculation formula is:

[0020] Ax+By+Cz+D=0

[0021]

[0022] (A, B, C) represents the normal vector of the plane formed by nearby ground points,

[0023] D represents the distance from the origin of the coordinate axis (0,0,0) to the plane.

[0024] (x, y, z) represents the coordinates of nearby ground points in the three coordinate axis directions.

[0025] (x0, y0, z0) represents the coordinates of the point cloud in the three coordinate axis directions;

[0026] d represents the distance from the point to the ground point directly below;

[0027] Then, using the nearby ground points as reference and the neighborhood points as support, we select 1-5 times the average distance between the nearby ground points as the neighborhood radius value, calculate the normal vector from each point in the point cloud to the nearby ground points, and assign it to the point attributes; the calculation formula is:

[0028] Ax+By+Cz+D=0

[0029]

[0030] Ax+By+Cz+D=0 represents the plane equation of the ground point,

[0031] (A, B, C) represents the normal vector of the plane formed by nearby ground points,

[0032] D represents the distance from the origin of the coordinate axis (0,0,0) to the plane.

[0033] (x1, y1, z1) represents the coordinates of any point in the three coordinate axes, namely point Q.

[0034] (x, y, z) represents the coordinates of any point on the ground in the three coordinate axis directions, namely point P.

[0035] represent exist projection on

[0036] Filter low vegetation point cloud: set the distance threshold by the minimum height difference between the building roof and the nearby ground, and obtain the distance value to the nearby ground point by traversing each point in the laser point cloud. Use the distance threshold as a reference to remove the laser point cloud of low vegetation, small objects close to the ground, and wall facades to obtain the processed laser point cloud; the reference range of the distance value is between 0.5 meters and 2 meters;

[0037] Eliminate irrelevant point clouds around buildings: set the density radius to calculate the number of points around each point in the point cloud, obtain the density of the point, determine the relationship between the density of the point and the set core point threshold, and determine it as a core point when the number of points is greater than or equal to the set core point threshold. When the density of the point is less than the set core point threshold, the point is within a certain core point density radius and is determined as a boundary point to form a cluster. The cluster category is determined according to the shape and density of the trees, and the clusters belonging to the tree point cloud are eliminated to obtain the building roof point cloud data. The reference range of the density radius is between 0.5 meters and 3 meters.

[0038] A method for automatically identifying a building roof structure based on an airborne point cloud, wherein the specific steps of segmenting and grouping the building roof surface laser point cloud are:

[0039] Select the initial seed point and base point: randomly select a point from the building roof point cloud as the initial seed point, randomly select another two points from the remaining point cloud as the base points of the plane, and add them to the current plane model;

[0040] Judge the normal vector and distance and construct the plane: judge the distance and normal vector consistency between the base point and the seed point. When the distance is less than or equal to the set basic distance threshold and the normal vector consistency, construct the plane. When the distance is greater than the set basic distance threshold or the normal vector is inconsistent, discard the base point and reselect the base point and calculate the distance and normal vector consistency between the new base point and the seed point. Loop until the base point and the seed point are all constructed with the plane, and cache the points that cannot be constructed as other points; the reference range of the basic distance threshold is between 10 meters and 100 meters, and the plane equation calculation formula is:

[0041] Ax+By+Cz+D=0

[0042] D=-(Ax4+By4+Cz4)

[0043] Represents the normal vector of the plane formed by the top surface of the building roof, that is,

[0044] (x, y, z) represents the coordinates of nearby ground points in the three coordinate axis directions.

[0045] (x4, y4, z4) represents the coordinates of the seed point in the three coordinate axis directions.

[0046] (x2, y2, z2) and (x3, y3, z3) represent the coordinates of two base points in the three coordinate axis directions.

[0047] D represents the distance from the origin of the coordinate axis (0,0,0) to the plane;

[0048] Get the inner points of the roof patch: Analyze and obtain the point cloud of the target object surface. The surface point cloud has up and down fluctuations. Set the fluctuation threshold, traverse the other points in the laser point cloud except the seed point and the basic point, obtain the distance between each other point and the constructed plane model, and determine the consistency between the normal vector of the other point and the normal vector of the constructed plane. When the distance is less than the fluctuation threshold and the normal vector of the other point is determined to be consistent with the normal vector of the constructed plane, the other point is classified as an inner point. When the distance is greater than or equal to the fluctuation threshold or the normal vector of the other point is determined to be inconsistent with the normal vector of the constructed plane, the other point is classified as an outer point. The reference range of the fluctuation threshold is between 0.10 meters and 0.30 meters.

[0049] Count the number of inliers that meet the conditions: Set a threshold for the number of inliers based on the minimum number of point clouds required for the plane model, and count the relationship between the number of inliers and the threshold. When the number of inliers is less than the threshold, all inliers are classified as non-planar discrete points. When the number of inliers is greater than or equal to the threshold, all inliers are classified as the same plane points to form the same plane point cloud. The same plane point clouds are grouped and assigned group numbers according to the order of judgment and screening.

[0050] Loop to select points in different planes: remove the remaining point cloud of the same plane points that have been assigned numbers, start from the step of selecting initial seed points and basic points to the step of counting the number of qualified internal points, select points in different planes and assign group numbers in the order of selection, loop this step until the remaining points cannot construct a plane, and end.

[0051] A method for automatically identifying a building roof structure based on an airborne point cloud, wherein the specific steps of fitting the roof point cloud to a set of plane equations are as follows:

[0052] Extracting laser point cloud of single roof patch: Based on the outline of a single building, extract the segmented point cloud of the building roof, and select the laser point cloud of the roof patch in the same group as the plane fitting object set according to the group number of the laser point cloud of the single roof patch. When the number in the attribute of the fitting object set is equal to 1, directly process the fitting object; when the number in the attribute of the fitting object set is greater than 1, process them as fitting objects in the order of the numbers.

[0053] Feature centering: After receiving a single fitting object, subtract the mean value M = (x5, y5, z5) of the point cloud in each dimension from the coordinate value of each dimension of the fitting object to obtain the feature-centered fitting object P. The mean value of each dimension after the transformation is zero. The calculation formula is:

[0054]

[0055] x i ,y i ,z i Represents the coordinates of the centralized point cloud;

[0056] (x5, y5, z5) represents the mean of the point cloud in this dimension, i.e., M;

[0057] p represents the feature-centered fitting object;

[0058] Establish the covariance matrix: calculate the feature-centered fitting object, calculate the covariance matrix, and obtain its eigenvalues ​​λ1, λ2, λ3 and eigenvectors V1, V2, V3. The calculation formula is:

[0059]

[0060] n represents the number of points in the fitting object,

[0061] C represents the covariance matrix of the three coordinates,

[0062] cov(x,y) represents the covariance of the x-coordinate and the y-coordinate.

[0063] cov(x,y) represents the variance of the x coordinate,

[0064] x i represents the x-coordinate of the i-th point in the fitting object,

[0065] y i Represents the y coordinate of the i-th point in the fitting object,

[0066] z i Represents the z coordinate of the i-th point in the fitting object;

[0067] Establish a set of plane equations and a set of normal vectors: select the eigenvector λ3 corresponding to the minimum eigenvalue in the covariance matrix to determine the plane equation; the eigenvector V3 corresponding to the minimum eigenvalue λ3 is the direction perpendicular to the fitting plane, V3 = (a, b, c) is the plane normal vector, the plane normal vector (a, b, c,) and the average position M of the point cloud determine the D1 value, D1 = -a*x5-b*y5-c*z5, and the final plane equation of the roof patch is: ax+by+cz+D1=0; the plane equations and normal vectors of the single fitting object after the traversal cycle are integrated to form a one-to-one corresponding set of plane equations and normal vectors.

[0068] A method for automatically identifying a building roof structure based on an airborne point cloud, wherein the specific steps of generating the roof structure type from the multiple fitting planes are as follows:

[0069] Calculate the area of ​​the building roof patch: Take the building outline as the range, combine the plane equation set of the building roof patch and the corresponding numbered point cloud, obtain the coordinates of the corner points of each plane, and calculate the polygonal area of ​​the building roof patch. The formula is:

[0070]

[0071] (x i ,y i ,z i ) represents the coordinates of the i-th point on the plane;

[0072] Represents the difference product of the vectors of the three corner points of a single plane.

[0073] (n1, n2, n3) represents the coordinates of the difference product of the vectors of the three corner points of a single plane in the directions of the three coordinate axes;

[0074] n represents the number of corner points of a single plane,

[0075] Represents the plane area composed of the first point, the i-th point, and the i+1-th point of a single plane.

[0076] S represents the area of ​​a single polygon;

[0077] Construct a valid building roof normal vector set: According to the polygonal area of ​​a single plane corresponding to the normal vector set, obtain the area mean and median, select the minimum of the area mean and median as the effective area threshold, remove the plane equations and corresponding normal vectors whose areas are less than the effective area threshold, and obtain a new plane equation set and normal vector set;

[0078] Determine the number of normal vectors: Count the number of normal vectors in the new normal vector set and perform determination processing. When it is determined that the number of normal vectors is equal to 1 and the absolute value of the normal vector Z value is less than the set normal vector Z threshold, a flat roof house structure type is directly generated. When the absolute value of the normal vector Z value is greater than or equal to the set normal vector Z threshold, a single slope roof house structure type is directly generated. When it is determined that the number of normal vectors is equal to 2, an area weight analysis process is performed. When it is determined that the number of normal vectors is greater than or equal to 3, a normal vector normalization process is performed. The reference range of the normal vector Z threshold is between 0 meters and 0.1 meters.

[0079] Normalization processing of normal vectors: screen the new normal vector set in the same direction, remove the normal vectors with small polygonal areas, and simultaneously delete the single plane equation corresponding to the normal vector in the plane equation set to obtain a new normal vector set and plane equation set; when the number of new normal vectors is equal to 1, the flat roof house structure type is directly generated when the absolute value of the normal vector Z value is less than the set normal vector Z threshold, and the single slope roof house structure type is directly generated when the absolute value of the normal vector Z value is greater than or equal to the set threshold; when the number of new normal vectors is equal to 2 and the number of new normal vectors is greater than or equal to 3, process the normal vector results respectively;

[0080] Area weight analysis and processing generates house structure type: when the number of normal vectors is equal to 2, determine the relationship between the plane areas corresponding to the two normal vectors; when the difference between the plane areas corresponding to the two normal vectors is greater than or equal to the set area threshold, obtain the normal vector corresponding to the larger polygon area, and directly determine the normal vector Z value. When the absolute value of the normal vector Z value is less than the set normal vector Z threshold, directly generate the flat roof house structure type; when the absolute value of the normal vector Z value is greater than or equal to the set normal vector Z threshold, directly generate the single slope roof house structure type; when the difference between the plane areas corresponding to the two normal vectors is less than the set area threshold, calculate the difference product of the two new normal vectors. When the absolute value of the difference product Z value is less than the set normal vector Z threshold, directly generate the mansard roof house structure type; when the absolute value of the difference product Z value is greater than or equal to the set normal vector Z threshold, directly generate the single slope roof house structure type; the reference range value of the area threshold is 3km 2 -5km 2 ;

[0081] Normal vector comparison analysis generates house structure type: when the number of normal vectors after normalization is greater than or equal to 3, they are sorted in descending order according to the polygon area value, and the normal vector with the largest polygon area is selected as the basic vector; when the difference between the current two areas is greater than or equal to 2 times the area threshold, the basic vector Z value is directly determined. When the absolute value of the basic vector Z value is less than the set normal vector Z threshold, the flat roof house structure type is directly generated; when the absolute value of the basic vector Z value is greater than or equal to the set normal vector Z threshold, the single slope roof house structure type is directly generated; when the difference between the current two areas is less than 2 times the area threshold, when the absolute value of the basic vector Z value is greater than or equal to the set normal vector Z threshold, the difference product of other normal vectors and the basic vector is traversed, when the absolute value of the difference product Z value is less than the set normal vector Z threshold, the human slope roof house structure type is directly generated, when the absolute value of the difference product Z value is greater than or equal to the set normal vector Z threshold, the single slope roof house structure type is directly generated; when the absolute value of the basic vector Z value is less than the set normal vector Z threshold, the flat roof house structure type is directly generated.

[0082] A method for automatically identifying a building roof structure based on an airborne point cloud, wherein the specific steps of verifying the building structure type are:

[0083] Classify the roof structure: by fitting the roof structure type of the raw building, assign 0 to the flat roof, 1 to the single-slope roof, and 2 to the gable roof, and record it in the corresponding building outline attribute;

[0084] Roof line extraction: Based on the segmented and grouped roof patch point cloud and building outline, obtain the seed point, traverse other points of the building roof patch, obtain the regional points and boundary points, and merge and fit to obtain the ridge line equation;

[0085] Judgment based on geometric features: when the ridge line equation cannot be fitted, the building structure type is a flat roof; when the ridge line equation is fitted, the point clouds on the left and right sides of the ridge line equation are judged. When there is a point on one side of the ridge line equation and the ridge line coincides with the edge of the roof, the building structure type is a single-slope roof; when there are points on both sides of the ridge line equation and the ridge line coincides with the edge of the roof, the building structure type is a human-slope roof;

[0086] Structural type modification: Compare and analyze the building structure type obtained with the corresponding building outline attribute structure type; when the two structure types are the same, directly determine that the structure type is correct; when the two structure types are different, add a building outline attribute field, fill in the structure type obtained by the verification method, and modify the structure type.

[0087] A method for automatically identifying a building roof structure based on airborne point clouds, wherein the specific steps of extracting the roof ridge line are:

[0088] Seed point selection: Based on the segmented roof patch point cloud and building outline, the height range of the roof point cloud [h min ,h max ], select a height close to the highest point h of the roof patch max The point is taken as the seed point, where h min Represents the lowest point of the roof patch;

[0089] Get regional points: Set the initial growth region R, including the seed point. Traverse the remaining roof patch points and calculate the distance and normal vector angle between each point and the seed point; when the distance is less than or equal to the set distance threshold dth and the normal vector angle is greater than or equal to the set normal vector angle threshold cosθth, d(P i ,P j )≤dth and ≥cosθth, place the point in the growth region R. When the distance is greater than the set distance threshold or the normal vector angle is greater than the set normal vector angle threshold, discard the point.

[0090] Boundary point determination: When the growth area R cannot continue to increase, check the boundary point of R. Select the boundary point P b , so that P b Belongs to R, there is a point P nb does not belong to R, so that P b To P nb Distance L(P b ,P nb ) is less than or equal to d neighbor , its neighbor is a distance slightly larger than dth.

[0091] Merge fitting: Get multiple growth areas (R1, R2, R3…R n ), calculate the adjacent growth area R i and R j The distance davg(R i ,R j ), determine davg(R i ,R j ) and the set merging threshold; when davg(R i ,R j ) is less than or equal to the merge threshold, then merge R i and R j , and fit the combined points to obtain the straight line equation of the ridge line. Its calculation formula is:

[0092]

[0093] n bi and n bj Represents R i and R j The number of boundary points,

[0094] and Represents a boundary point.

[0095] It can be seen from this that:

[0096] The method in the embodiment of the present invention uses airborne radar to quickly obtain high-precision three-dimensional point clouds of urban building roof surfaces over a large area, and has the advantages of high precision and high penetration. At the same time, an improved algorithm can be used to segment building roof surfaces and fit roof plane equations. It is a high-precision and efficient data acquisition method for quickly obtaining urban three-dimensional model structure types. BRIEF DESCRIPTION OF THE DRAWINGS

[0097] Figure 1 An overall flow chart of a method for automatically identifying a building roof structure based on an airborne point cloud is provided in an embodiment of the present invention;

[0098] Figure 2 A schematic flow chart of a building house extraction step in a method for automatically identifying a building roof structure based on an airborne point cloud provided in an embodiment of the present invention;

[0099] Figure 3 A schematic flow chart of the steps of segmenting and grouping the laser point cloud of the roof surface of a building in a method for automatically identifying the roof structure of a building based on an airborne point cloud provided in an embodiment of the present invention;

[0100] Figure 4A schematic flow chart of a set of plane equation fitting steps of a roof point cloud in a method for automatically identifying a building roof structure based on an airborne point cloud provided in an embodiment of the present invention;

[0101] Figure 5 A schematic flow chart of steps for generating roof structure types from multiple fitting planes in a method for automatically identifying building roof structures based on airborne point clouds provided in an embodiment of the present invention;

[0102] Figure 6 A schematic flow chart of a building structure type verification step in a method for automatically identifying a building roof structure based on an airborne point cloud provided in an embodiment of the present invention;

[0103] Figure 7 A schematic flow chart of a ridge line extraction step in a method for automatically identifying a building roof structure based on an airborne point cloud provided in an embodiment of the present invention;

[0104] Figure 8 A diagram showing ground points around a building and low-valued point clouds in a method for automatically identifying a building roof structure based on an airborne point cloud provided in an embodiment of the present invention;

[0105] Fig. 9 A point cloud presentation diagram of a building roof surface in a method for automatically identifying a building roof structure based on an airborne point cloud provided in an embodiment of the present invention;

[0106] Fig.10 A diagram showing the segmentation and grouping of laser point clouds of building roof panels in a method for automatically identifying building roof structures based on airborne point clouds provided in an embodiment of the present invention;

[0107] Fig.11 A diagram showing the fitting plane and intersection lines and points of a building roof patch in a method for automatically identifying a building roof structure based on an airborne point cloud provided in an embodiment of the present invention;

[0108] Fig.12 A diagram presenting building outlines and structure types in a method for automatically identifying building roof structures based on airborne point clouds provided in an embodiment of the present invention. DETAILED DESCRIPTION

[0109] In order to enable those skilled in the art to better understand the present invention, the present invention will be described in detail below with reference to the accompanying drawings and specific embodiments. The exemplary embodiments and descriptions of the present invention are used to explain the present invention but are not intended to limit the present invention.

[0110] Embodiment 1:

[0111] Figure 1A flow chart of a method for automatically identifying a building roof structure based on an airborne point cloud provided in this embodiment is shown in FIG. Figure 1 As shown, a method for automatically identifying building roof structures based on airborne point clouds includes the following steps: laser point cloud preprocessing, building house extraction, segmentation and grouping of building roof surface laser point clouds, roof point cloud fitting plane equation set, multiple fitting planes to generate roof structure types, house structure type verification, wherein:

[0112] Point cloud preprocessing: The original laser point cloud data is obtained through the mutual coordination and combined operation between the laser transmitter, inertial navigation and GNSS equipment. The laser point cloud is matched with the point cloud of a single flight and a single flight strip to obtain the registration point cloud of each flight. Then, the rigid registration point cloud of each flight is regionally optimized to obtain the optimized point cloud of the entire survey area. The optimized point cloud of the survey area is then converted from geodetic coordinates to the standard national 2000 coordinate system and the local city coordinate system to obtain the target point cloud of the survey area. The target point cloud of the survey area is converted into point cloud data in an internationally accepted format.

[0113] Building extraction: remove noise points from the preprocessed point cloud, extract the ground point cloud, calculate the point cloud distance and normal vector, filter the low-value point cloud according to the distance and normal vector factors, delete the clustered tree point cloud around the building, and obtain the building roof point cloud data;

[0114] Segmentation and grouping of building roof laser point cloud: Randomly select initial seed points and basic points from the building roof point cloud to construct a plane by normal vector and distance determination. For the points that cannot construct a plane, use them as other points to determine the distance and normal vector between the plane and determine the inner points and outer points. After grouping the inner points, perform a cycle of inner point selection on different planes and group them in order.

[0115] Roof point cloud fitting plane equation set: Based on the outline of a single building, extract the laser point cloud of a single roof surface after segmentation and grouping. According to the grouping number of the laser point cloud of a single roof surface, select the laser point cloud of the roof surface in the same group as the plane fitting object, perform feature centering on the fitting object, establish the covariance matrix to obtain the eigenvalue and its corresponding eigenvector, and establish the plane equation set and normal vector set for the single fitting object after traversal processing;

[0116] Generate roof structure type from multiple fitting planes: Take the building outline as the range, obtain the polygonal area of ​​the building roof patch, remove the plane equations and corresponding normal vectors with an area less than the threshold, obtain a new set of plane equations and normal vectors, and determine the number of normal vectors, perform normal vector normalization, area weight analysis, and normal vector comparison analysis to generate the house structure type;

[0117] Verification of building structure type: classify the building roof structure type and record it in the corresponding building outline attribute. Use the segmented and grouped roof patch point cloud and building outline as the data basis to obtain the ridge line. Predict the structure type based on the geometric features and ridge line. Determine the predicted structure type and compare it with the structure type in the corresponding building outline attribute. Modify the inconsistent structure type.

[0118] like Figure 2 As shown, a method for automatically identifying building roof structures based on airborne point clouds, the specific steps of extracting buildings are:

[0119] Eliminate noise points: Remove extremely high and extremely low discrete points and noise points from the general format point cloud data to obtain a three-dimensional laser point cloud including topography, buildings and surrounding small pieces, trees, power lines, etc. within the target measurement area;

[0120] Extract ground point cloud: Use opening and closing operations on the three-dimensional laser points within the target measurement area to preliminarily determine the ground points, remove isolated high points and fill small depressions to obtain a relatively smooth surface model, compare the point cloud with the constructed surface model, and reclassify the ground points and non-ground points; points that are close to the surface model are determined as ground points; points that are far away are determined as non-ground points. Update the surface model based on the new classification results, and repeat the above process until the classification results are stable or a certain number of iterations are reached, and obtain a three-dimensional laser point cloud containing ground points and non-ground points;

[0121] Calculate the point cloud distance and normal vector: Calculate the distance from each point to the ground point directly below and assign it to the point cloud attribute; the calculation formula is:

[0122] Ax+By+Cz+D=0

[0123]

[0124] (A, B, C) represents the normal vector of the plane formed by nearby ground points,

[0125] D represents the distance from the origin of the coordinate axis (0,0,0) to the plane.

[0126] (x, y, z) represents the coordinates of nearby ground points in the three coordinate axis directions.

[0127] (x0, y0, z0) represents the coordinates of the point cloud in the three coordinate axis directions;

[0128] d represents the distance from the point to the ground point directly below;

[0129] Then, using the nearby ground points as reference and the neighborhood points as support, we select 1-5 times the average distance between the nearby ground points as the neighborhood radius value, calculate the normal vector from each point in the point cloud to the nearby ground points, and assign it to the point attributes; the calculation formula is:

[0130] Ax+By+Cz+D=0

[0131]

[0132] Ax+By+Cz+D=0 represents the plane equation of the ground point,

[0133] (A, B, C) represents the normal vector of the plane formed by nearby ground points,

[0134] D represents the distance from the origin of the coordinate axis (0,0,0) to the plane.

[0135] (x1, y1, z1) represents the coordinates of any point in the three coordinate axes, namely point Q.

[0136] (x, y, z) represents the coordinates of any point on the ground in the three coordinate axis directions, namely point P.

[0137] represent exist projection on

[0138] Filter low vegetation point cloud: set the distance threshold by the minimum height difference between the building roof and the nearby ground, and obtain the distance value to the nearby ground point by traversing each point in the laser point cloud. Use the distance threshold as a reference to remove the laser point cloud of low vegetation, small objects close to the ground, and wall facades to obtain the processed laser point cloud; the reference range of the distance value is between 0.5 meters and 2 meters;

[0139] Eliminate irrelevant point clouds around buildings: set the density radius to calculate the number of points around each point in the point cloud, obtain the density of the point, determine the relationship between the density of the point and the set core point threshold, and determine it as a core point when the number of points is greater than or equal to the set core point threshold. When the density of the point is less than the set core point threshold, the point is within a certain core point density radius and is determined as a boundary point to form a cluster. The cluster category is determined according to the shape and density of the trees, and the clusters belonging to the tree point cloud are eliminated to obtain the building roof point cloud data. The reference range of the density radius is between 0.5 meters and 3 meters.

[0140] like Figure 3 As shown, a method for automatically identifying a building roof structure based on an airborne point cloud, the specific steps of segmenting and grouping the building roof surface laser point cloud are:

[0141] Select the initial seed point and base point: randomly select a point from the building roof point cloud as the initial seed point, randomly select another two points from the remaining point cloud as the base points of the plane, and add them to the current plane model;

[0142] Judge the normal vector and distance and construct the plane: judge the distance and normal vector consistency between the base point and the seed point. When the distance is less than or equal to the set basic distance threshold and the normal vector consistency, construct the plane. When the distance is greater than the set basic distance threshold or the normal vector is inconsistent, discard the base point and reselect the base point and calculate the distance and normal vector consistency between the new base point and the seed point. Loop until the base point and the seed point are all constructed with the plane, and cache the points that cannot be constructed as other points; the reference range of the basic distance threshold is between 10 meters and 100 meters, and the plane equation calculation formula is:

[0143] Ax+By+Cz+D=0

[0144] D=-(Ax4+By4+Cz4)

[0145] Represents the normal vector of the plane formed by the top surface of the building roof, that is,

[0146] (x, y, z) represents the coordinates of nearby ground points in the three coordinate axis directions.

[0147] (x4, y4, z4) represents the coordinates of the seed point in the three coordinate axis directions.

[0148] (x2, y2, z2) and (x3, y3, z3) represent the coordinates of the two base points in the three coordinate axis directions, and D represents the distance from the origin of the coordinate axis (0,0,0) to the plane;

[0149] Get the inner points of the roof patch: Analyze and obtain the point cloud of the target object surface. The surface point cloud has up and down fluctuations. Set the fluctuation threshold, traverse the other points in the laser point cloud except the seed point and the basic point, obtain the distance between each other point and the constructed plane model, and determine the consistency between the normal vector of the other point and the normal vector of the constructed plane. When the distance is less than the fluctuation threshold and the normal vector of the other point is determined to be consistent with the normal vector of the constructed plane, the other point is classified as an inner point. When the distance is greater than or equal to the fluctuation threshold or the normal vector of the other point is determined to be inconsistent with the normal vector of the constructed plane, the other point is classified as an outer point. The reference range of the fluctuation threshold is between 0.10 meters and 0.30 meters.

[0150] Count the number of inliers that meet the conditions: Set a threshold for the number of inliers based on the minimum number of point clouds required for the plane model, and count the relationship between the number of inliers and the threshold. When the number of inliers is less than the threshold, all inliers are classified as non-planar discrete points. When the number of inliers is greater than or equal to the threshold, all inliers are classified as the same plane points to form the same plane point cloud. The same plane point clouds are grouped and assigned group numbers according to the order of judgment and screening.

[0151] Loop to select points in different planes: remove the remaining point cloud of the same plane points that have been assigned numbers, start from the step of selecting initial seed points and basic points to the step of counting the number of qualified internal points, select points in different planes and assign group numbers in the order of selection, loop this step until the remaining points cannot construct a plane, and end.

[0152] like Figure 4 As shown, a method for automatically identifying the roof structure of a building based on airborne point cloud, the specific steps of fitting the roof point cloud to a set of plane equations are:

[0153] Extracting laser point cloud of single roof patch: Based on the outline of a single building, extract the segmented point cloud of the building roof, and select the laser point cloud of the roof patch in the same group as the plane fitting object set according to the group number of the laser point cloud of the single roof patch. When the number in the attribute of the fitting object set is equal to 1, directly process the fitting object; when the number in the attribute of the fitting object set is greater than 1, process them as fitting objects in the order of the numbers.

[0154] Feature centering: After receiving a single fitting object, subtract the mean value M = (x5, y5, z5) of the point cloud in each dimension from the coordinate value of each dimension of the fitting object to obtain the feature-centered fitting object P. The mean value of each dimension after the transformation is zero. The calculation formula is:

[0155]

[0156] x i ,y i ,z i Represents the coordinates of the centralized point cloud;

[0157] (x5, y5, z5) represents the mean of the point cloud in this dimension, i.e., M;

[0158] p represents the feature-centered fitting object;

[0159] Establish the covariance matrix: calculate the feature-centered fitting object, calculate the covariance matrix, and obtain its eigenvalues ​​λ1, λ2, λ3 and eigenvectors V1, V2, V3. The calculation formula is:

[0160]

[0161] n represents the number of points in the fitting object,

[0162] C represents the covariance matrix of the three coordinates,

[0163] cov(x,y) represents the covariance of the x-coordinate and the y-coordinate.

[0164] cov(x,y) represents the variance of the x coordinate,

[0165] x i represents the x-coordinate of the i-th point in the fitting object,

[0166] y i Represents the y coordinate of the i-th point in the fitting object,

[0167] z i Represents the z coordinate of the i-th point in the fitting object;

[0168] Establish a set of plane equations and a set of normal vectors: select the eigenvector λ3 corresponding to the minimum eigenvalue in the covariance matrix to determine the plane equation; the eigenvector V3 corresponding to the minimum eigenvalue λ3 is the direction perpendicular to the fitting plane, V3 = (a, b, c) is the plane normal vector, the plane normal vector (a, b, c,) and the average position M of the point cloud determine the D1 value, D1 = -a*x5-b*y5-c*z5, and the final plane equation of the roof patch is: ax+by+cz+D1=0; the plane equations and normal vectors of the single fitting object after the traversal cycle are integrated to form a one-to-one corresponding set of plane equations and normal vectors.

[0169] like Figure 5 As shown, a method for automatically identifying a building roof structure based on an airborne point cloud, wherein the specific steps of generating the roof structure type by using multiple fitting planes are as follows:

[0170] Calculate the area of ​​the building roof patch: Take the building outline as the range, combine the plane equation set of the building roof patch and the corresponding numbered point cloud, obtain the coordinates of the corner points of each plane, and calculate the polygonal area of ​​the building roof patch. The formula is:

[0171]

[0172] (x i ,y i ,z i ) represents the coordinates of the i-th point on the plane;

[0173] Represents the difference product of the vectors of the three corner points of a single plane.

[0174] (n1, n2, n3) represents the coordinates of the difference product of the vectors of the three corner points of a single plane in the directions of the three coordinate axes;

[0175] n represents the number of corner points of a single plane,

[0176] Represents the plane area composed of the first point, the i-th point, and the i+1-th point of a single plane.

[0177] S represents the area of ​​a single polygon;

[0178] Construct a valid building roof normal vector set: According to the polygonal area of ​​a single plane corresponding to the normal vector set, obtain the area mean and median, select the minimum of the area mean and median as the effective area threshold, remove the plane equations and corresponding normal vectors whose areas are less than the effective area threshold, and obtain a new plane equation set and normal vector set;

[0179] Determine the number of normal vectors: Count the number of normal vectors in the new normal vector set and perform determination processing. When it is determined that the number of normal vectors is equal to 1 and the absolute value of the normal vector Z value is less than the set normal vector Z threshold, a flat roof house structure type is directly generated. When the absolute value of the normal vector Z value is greater than or equal to the set normal vector Z threshold, a single slope roof house structure type is directly generated. When it is determined that the number of normal vectors is equal to 2, an area weight analysis process is performed. When it is determined that the number of normal vectors is greater than or equal to 3, a normal vector normalization process is performed. The reference range of the normal vector Z threshold is between 0 meters and 0.1 meters.

[0180] Normalization processing of normal vectors: screen the new normal vector set in the same direction, remove the normal vectors with small polygonal areas, and simultaneously delete the single plane equation corresponding to the normal vector in the plane equation set to obtain a new normal vector set and plane equation set; when the number of new normal vectors is equal to 1, the flat roof house structure type is directly generated when the absolute value of the normal vector Z value is less than the set normal vector Z threshold, and the single slope roof house structure type is directly generated when the absolute value of the normal vector Z value is greater than or equal to the set threshold; when the number of new normal vectors is equal to 2 and the number of new normal vectors is greater than or equal to 3, process the normal vector results respectively;

[0181] Area weight analysis and processing generates house structure type: when the number of normal vectors is equal to 2, determine the relationship between the plane areas corresponding to the two normal vectors; when the difference between the plane areas corresponding to the two normal vectors is greater than or equal to the set area threshold, obtain the normal vector corresponding to the larger polygon area, and directly determine the normal vector Z value. When the absolute value of the normal vector Z value is less than the set normal vector Z threshold, directly generate the flat roof house structure type; when the absolute value of the normal vector Z value is greater than or equal to the set normal vector Z threshold, directly generate the single slope roof house structure type; when the difference between the plane areas corresponding to the two normal vectors is less than the set area threshold, calculate the difference product of the two new normal vectors. When the absolute value of the difference product Z value is less than the set normal vector Z threshold, directly generate the mansard roof house structure type; when the absolute value of the difference product Z value is greater than or equal to the set normal vector Z threshold, directly generate the single slope roof house structure type; the reference range value of the area threshold is 3km 2 -5km 2 ;

[0182] Normal vector comparison analysis generates house structure type: when the number of normalized normal vectors is greater than or equal to 3, they are sorted in descending order according to the polygon area value, and the normal vector with the largest polygon area is selected as the basic vector; when the difference between the current two areas is greater than or equal to 2 times the area threshold, the basic vector Z value is directly judged. When the absolute value of the basic vector Z is less than the set normal vector Z threshold, the flat roof house structure type is directly generated; when the absolute value of the basic vector Z value is greater than or equal to the set normal vector Z threshold, the single slope roof house structure type is directly generated; when the difference between the current two areas is less than 2 times the area threshold, when the absolute value of the basic vector Z value is greater than or equal to the set normal vector Z threshold, the difference product of other normal vectors and the basic vector is traversed, when the absolute value of the difference product Z value is less than the set normal vector Z threshold, the man-slope roof house structure type is directly generated, when the absolute value of the difference product Z value is greater than or equal to the set normal vector Z threshold, the single slope roof house structure type is directly generated; when the absolute value of the basic vector Z value is less than the set normal vector Z threshold, the flat roof house structure type is directly generated.

[0183] like Figure 6 As shown, a method for automatically identifying the roof structure of a building based on airborne point cloud, the specific steps of the house structure type verification are:

[0184] Classify the roof structure: by fitting the roof structure type of the raw building, assign 0 to the flat roof, 1 to the single-slope roof, and 2 to the gable roof, and record it in the corresponding building outline attribute;

[0185] Roof line extraction: Based on the segmented and grouped roof patch point cloud and building outline, obtain the seed point, traverse other points of the building roof patch, obtain the regional points and boundary points, and merge and fit to obtain the ridge line equation;

[0186] Judgment based on geometric features: when the ridge line equation cannot be fitted, the building structure type is a flat roof; when the ridge line equation is fitted, the point clouds on the left and right sides of the ridge line equation are judged. When there is a point on one side of the ridge line equation and the ridge line coincides with the edge of the roof, the building structure type is a single-slope roof; when there are points on both sides of the ridge line equation and the ridge line coincides with the edge of the roof, the building structure type is a human-slope roof;

[0187] Structural type modification: Compare and analyze the building structure type obtained with the corresponding building outline attribute structure type; when the two structure types are the same, directly determine that the structure type is correct; when the two structure types are different, add a building outline attribute field, fill in the structure type obtained by the verification method, and modify the structure type.

[0188] like Figure 7 As shown, a method for automatically identifying the roof structure of a building based on airborne point cloud, the specific steps of extracting the roof ridge line are:

[0189] Seed point selection: Based on the segmented roof patch point cloud and building outline, the height range of the roof point cloud [h min ,h max ], select a height close to the highest point h of the roof patch max The point is taken as the seed point, where h min Represents the lowest point of the roof patch;

[0190] Get regional points: Set the initial growth region R, including the seed point. Traverse the remaining roof patch points and calculate the distance and normal vector angle between each point and the seed point; when the distance is less than or equal to the set distance threshold dth and the normal vector angle is greater than or equal to the set normal vector angle threshold cosθth, d(Pi,Pj)≤dth and ≥cosθth, place the point in the growth region R. When the distance is greater than the set distance threshold or the normal vector angle is greater than the set normal vector angle threshold, discard the point.

[0191] Boundary point determination: When the growth area R cannot continue to increase, check the boundary point of R. Select the boundary point P b , so that P b Belongs to R, there is a point P nb does not belong to R, so that P b To P nb Distance L(P b ,P nb ) is less than or equal to d neighbor , its neighbor is a distance slightly larger than dth.

[0192] Merge fitting: Get multiple growth areas (R1, R2, R3…Rn ), calculate the adjacent growth area R i and R j The distance davg(R i ,R j ), determine davg(R i ,R j ) and the set merging threshold; when davg(R i ,R j ) is less than or equal to the merge threshold, then merge R i and R j , and fit the combined points to obtain the straight line equation of the ridge line. Its calculation formula is:

[0193]

[0194] n bi and n bj Represents R i and R j The number of boundary points,

[0195] and Represents a boundary point.

[0196] The following is a detailed description of a specific implementation case:

[0197] The collected airborne lidar data in the built-up area of ​​1,200 square kilometers were processed. A total of three flights were made in the survey area. The first flight covered 8 routes, the second flight covered 7 routes, and the third flight covered 8 routes, with a total of 23 routes. The point cloud density was greater than 16 points per square meter. Preprocess the airborne laser radar point cloud, optimize the registration of each flight strip laser point cloud, and then perform coordinate conversion and format conversion on the matched optimized point cloud; remove noise points and filter low-value point clouds from the preprocessed point cloud to obtain the building roof point cloud data; randomly select initial seed points and basic points from the building roof point cloud to determine the normal vector and distance to construct a plane, determine the inner point and outer point, group the inner points, and then perform inner point selection cycles on different planes and group them in sequence; based on the outline of a single building, extract the laser point cloud of the single roof surface after segmentation and grouping, and establish a set of plane equations and normal vectors for each building; classify the building roof structure type and record it in the corresponding building outline attribute, obtain the ridge line based on the segmented and grouped roof surface point cloud and building outline, predict the structure type based on the geometric features and the ridge line, judge the predicted structure type and compare it with the structure type in the corresponding building outline attribute, and modify the inconsistent structure type.

[0198] like Figure 1A method for automatically identifying a building roof structure based on an airborne point cloud is shown, the method comprising the following steps:

[0199] Point cloud preprocessing: Through the mutual coordination and combined operation between the laser transmitter, inertial navigation and GNSS equipment, the projection zone UTM50 laser point cloud data is obtained, and the laser point clouds of the first flight 8 routes, the second flight 7 routes, and the third flight 8 routes are matched to obtain the registration point cloud of each flight, and then the rigid registration point cloud of each flight is regionally optimized to obtain the optimized point cloud of the entire survey area of ​​1,200 square kilometers; then the optimized point cloud of the survey area is converted from geodetic coordinates to the 2000 Tianjin city coordinate system to obtain the target point cloud of the survey area; the target point cloud of the survey area is converted into the internationally accepted LAS format point cloud data;

[0200] Building extraction: remove noise points from the preprocessed point cloud, extract the ground point cloud, calculate the point cloud distance and normal vector, filter the low-value point cloud according to the distance and normal vector factors, delete the clustered tree point cloud around the building, and obtain the building roof point cloud data;

[0201] Segmentation and grouping of building roof laser point cloud: Randomly select initial seed points and basic points from the building roof point cloud to construct a plane by normal vector and distance determination. For the points that cannot construct a plane, use them as other points to determine the distance and normal vector between the plane and determine the inner points and outer points. After grouping the inner points, perform a cycle of inner point selection on different planes and group them in order.

[0202] Roof point cloud fitting plane equation set: Based on the outline of a single building, extract the laser point cloud of a single roof surface after segmentation and grouping. According to the grouping number of the laser point cloud of a single roof surface, select the laser point cloud of the roof surface in the same group as the plane fitting object, perform feature centering on the fitting object, establish the covariance matrix to obtain the eigenvalue and its corresponding eigenvector, and establish the plane equation set and normal vector set for the single fitting object after traversal processing;

[0203] Generate roof structure type from multiple fitting planes: Take the building outline as the range, obtain the polygonal area of ​​the building roof patch, remove the plane equations and corresponding normal vectors with an area less than the threshold, obtain a new set of plane equations and normal vectors, and determine the number of normal vectors, perform normal vector normalization, area weight analysis, and normal vector comparison analysis to generate the house structure type;

[0204] Verification of building structure type: classify the building roof structure type and record it in the corresponding building outline attribute. Use the segmented and grouped roof patch point cloud and building outline as the data basis to obtain the ridge line. Predict the structure type based on the geometric features and ridge line. Determine the predicted structure type and compare it with the structure type in the corresponding building outline attribute. Modify the inconsistent structure type.

[0205] like Figure 2 A method for automatically identifying building roof structures based on airborne point clouds is shown, and the specific steps of extracting the building house are:

[0206] Eliminate noise points: Remove extremely high discrete points and noise points with an absolute elevation of more than 350 meters and an absolute elevation of less than -10 meters from the general LAS format point cloud data, and obtain a three-dimensional laser point cloud including topography, buildings and surrounding small pieces, trees, power lines, etc. within the target survey area;

[0207] Extract ground point cloud: Use opening and closing operations on the three-dimensional laser points within the target measurement area to preliminarily determine the ground points, remove isolated high points and fill small depressions to obtain a relatively smooth surface model, compare the point cloud with the constructed surface model, and reclassify the ground points and non-ground points; points that are close to the surface model are determined as ground points; points that are far away are determined as non-ground points. Update the surface model based on the new classification results, and repeat the above process until the classification results are stable or a certain number of iterations are reached, and obtain a three-dimensional laser point cloud containing ground points and non-ground points;

[0208] Calculate the point cloud distance and normal vector: Calculate the distance from each point to the ground point directly below and assign it to the point cloud attribute; the calculation formula is:

[0209] Ax+By+Cz+D=0

[0210]

[0211] (A, B, C) represents the normal vector of the plane formed by nearby ground points,

[0212] D represents the distance from the origin of the coordinate axis (0,0,0) to the plane.

[0213] (x, y, z) represents the coordinates of nearby ground points in the three coordinate axis directions.

[0214] (x0, y0, z0) represents the coordinates of the point cloud in the three coordinate axis directions;

[0215] d represents the distance from the point to the ground point directly below;

[0216] Then, taking the nearby ground points as reference and the neighborhood points as support, we select twice the average distance between the nearby ground points as the neighborhood radius value, calculate the normal vector from each point in the point cloud to the nearby ground points, and assign it to the attributes of the point; the calculation formula is:

[0217] Ax+By+Cz+D=0

[0218]

[0219] Ax+By+Cz+D=0 represents the plane equation of the ground point,

[0220] (A, B, C) represents the normal vector of the plane formed by nearby ground points,

[0221] D represents the distance from the origin of the coordinate axis (0,0,0) to the plane.

[0222] (x1, y1, z1) represents the coordinates of any point in the three coordinate axes, namely point Q.

[0223] (x, y, z) represents the coordinates of any point on the ground in the three coordinate axis directions, namely point P.

[0224] represent exist projection on

[0225] Filter low vegetation point clouds: Set the distance threshold to 1.5 meters based on the minimum height difference between the building roof and the nearby ground. By traversing each point in the laser point cloud, obtain the distance value to the nearby ground point. Use the distance threshold of 1.5 meters as a reference to remove low vegetation, small pieces close to the ground, and laser point clouds of wall facades, such as Figure 8 The ground points and low values ​​around the shown building are shown by the point cloud, and the processed laser point cloud is obtained;

[0226] Eliminate irrelevant point clouds around buildings: set the density radius, calculate the number of points around each point in the point cloud, obtain the density of the point, determine the relationship between the density of the point and the set core point threshold, and determine it as a core point when the number of points is greater than or equal to the set core point threshold. When the density of the point is less than the set core point threshold, the point is within a certain core point density radius and is determined as a boundary point to form a cluster. According to the shape and density of the tree, the cluster category is determined, and the cluster belonging to the tree point cloud is eliminated to obtain the building roof point cloud data, such as Fig. 9 The point cloud of the roof patch of the building is shown, where the density radius is set to 1 meter and the core point threshold is set to 12 points per square meter;

[0227] like Figure 3A method for automatically identifying a building roof structure based on an airborne point cloud is shown, wherein the specific steps of segmenting and grouping the building roof surface laser point cloud are as follows:

[0228] Select the initial seed point and base point: randomly select a point from the building roof point cloud as the initial seed point, randomly select another two points from the remaining point cloud as the base points of the plane, and add them to the current plane model;

[0229] Judge the normal vector and distance and construct the plane: judge the distance and normal vector consistency between the base point and the seed point. When the distance is less than or equal to the set basic distance threshold and the normal vector consistency, build the plane. When the distance is greater than the set basic distance threshold or the normal vector is inconsistent, discard the base point and reselect the base point and calculate the distance and normal vector consistency between the new base point and the seed point. Loop until the base point and the seed point are all constructed with the plane, and cache the points that cannot be constructed as other points; the basic distance threshold is set to 80 meters, and the plane equation calculation formula is:

[0230] Ax+By+Cz+D=0

[0231] D=-(Ax4+By4+Cz4)

[0232] Represents the normal vector of the plane formed by the top surface of the building roof, that is,

[0233] (x, y, z) represents the coordinates of nearby ground points in the three coordinate axis directions.

[0234] (x4, y4, z4) represents the coordinates of the seed point in the three coordinate axis directions.

[0235] (x2, y2, z2) and (x3, y3, z3) represent the coordinates of the two base points in the three coordinate axis directions, and D represents the distance from the origin of the coordinate axis (0,0,0) to the plane;

[0236] Get the inner points of the roof patch: Analyze and obtain the point cloud of the target object surface. The surface point cloud has up and down fluctuations. Set the fluctuation threshold, traverse the other points in the laser point cloud except the seed point and the basic point, obtain the distance between each other point and the constructed plane model, and determine the consistency between the normal vector of the other point and the normal vector of the constructed plane. When the distance is less than the fluctuation threshold and the normal vector of the other point is determined to be consistent with the normal vector of the constructed plane, the other point is classified as an inner point. When the distance is greater than or equal to the fluctuation threshold or the normal vector of the other point is determined to be inconsistent with the normal vector of the constructed plane, the other point is classified as an outer point. The fluctuation threshold is 0.15 meters.

[0237] Count the number of qualified inliers: Based on the minimum number of point clouds required for the plane model, set the number threshold to 60, and count the relationship between the number of inliers and the number threshold. When the number of inliers is less than the number threshold, all inliers are classified as non-planar discrete points; when the number of inliers is greater than or equal to the number threshold, all inliers are classified as the same plane points to form the same plane point cloud, group the same plane point clouds, and assign group numbers according to the order of judgment and screening. The group number consists of 10 characters, 1-4 are the map number, 5-8 are the sequence number of the building outline, and 9-10 are the group sequence numbers of the laser point cloud within the building outline;

[0238] Loop to select points in different planes: remove the numbered points in the same plane from the remaining point cloud to the step of selecting the initial seed point and the basic point to the step of counting the number of qualified points to select points in different planes and assign group numbers in the order of selection. Loop this step until the remaining points cannot construct a plane, and then end. Fig.10 The laser point cloud segmentation grouping of the building roof surface is shown.

[0239] like Figure 4 A method for automatically identifying a building roof structure based on an airborne point cloud is shown, wherein the specific steps of fitting the roof point cloud to a set of plane equations are as follows:

[0240] Extracting laser point cloud of single roof patch: Based on the outline of a single building, extract the segmented point cloud of the building roof, and select the laser point cloud of the roof patch in the same group as the plane fitting object set according to the group number of the laser point cloud of the single roof patch. When the number in the attribute of the fitting object set is equal to 1, directly process the fitting object; when the number in the attribute of the fitting object set is greater than 1, process them as fitting objects in the order of the numbers.

[0241] Feature centering: After receiving a single fitting object, subtract the mean value M = (x5, y5, z5) of the point cloud in each dimension from the coordinate value of each dimension of the fitting object to obtain the feature-centered fitting object P. The mean value of each dimension after the transformation is zero. The calculation formula is:

[0242]

[0243]

[0244] x i ,y i ,z i Represents the coordinates of the centralized point cloud;

[0245] (x5, y5, z5) represents the mean of the point cloud in this dimension, i.e., M;

[0246] p represents the feature-centered fitting object;

[0247] Establish the covariance matrix: calculate the feature-centered fitting object, calculate the covariance matrix, and obtain its eigenvalues ​​λ1, λ2, λ3 and eigenvectors V1, V2, V3. The calculation formula is:

[0248]

[0249] n represents the number of points in the fitting object,

[0250] C represents the covariance matrix of the three coordinates,

[0251] cov(x,y) represents the covariance of the x-coordinate and the y-coordinate.

[0252] cov(x,y) represents the variance of the x coordinate,

[0253] x i represents the x-coordinate of the i-th point in the fitting object,

[0254] y i Represents the y coordinate of the i-th point in the fitting object,

[0255] z i Represents the z coordinate of the i-th point in the fitting object;

[0256] Establish a set of plane equations and a set of normal vectors: select the eigenvector λ3 corresponding to the minimum eigenvalue in the covariance matrix to determine the plane equation; the eigenvector V3 corresponding to the minimum eigenvalue λ3 is the direction perpendicular to the fitting plane, V3 = (a, b, c) is the plane normal vector, the plane normal vector (a, b, c,) and the average position M of the point cloud determine the D1 value, D1 = -a*x5-b*y5-c*z5, and the final plane equation of the roof patch is: ax+by+cz+D1=0; the plane equations and normal vectors of the single fitting object after the traversal cycle are integrated to form a one-to-one corresponding set of plane equations and normal vectors.

[0257] like Figure 5 A method for automatically identifying a building roof structure based on an airborne point cloud is shown, wherein the specific steps of generating the roof structure type by using multiple fitting planes are as follows:

[0258] Calculate the area of ​​the building roof patch: Take the building outline as the range, combine the plane equation set of the building roof patch and the corresponding numbered point cloud, and obtain the coordinates of the corner points of each plane, such as Fig.11 The polygonal area of ​​the building roof patch is calculated by fitting the plane and the intersection line and point shown in the figure, where the formula is:

[0259]

[0260]

[0261] (x i ,y i ,z i ) represents the coordinates of the i-th point on the plane;

[0262] Represents the difference product of the vectors of the three corner points of a single plane.

[0263] (n1, n2, n3) represents the coordinates of the difference product of the vectors of the three corner points of a single plane in the directions of the three coordinate axes;

[0264] n represents the number of corner points of a single plane,

[0265] Represents the plane area composed of the first point, the i-th point, and the i+1-th point of a single plane.

[0266] S represents the area of ​​a single polygon;

[0267] Construct a valid building roof normal vector set: According to the polygonal area of ​​a single plane corresponding to the normal vector set, obtain the area mean and median, select the minimum of the area mean and median as the effective area threshold, remove the plane equations and corresponding normal vectors whose areas are less than the effective area threshold, and obtain a new plane equation set and normal vector set;

[0268] Determine the number of normal vectors: Count the number of normal vectors in the new normal vector set and perform determination processing. When the number of normal vectors is determined to be equal to 1, the absolute value of the normal vector Z value is less than the set normal vector Z threshold, and the flat roof house structure type is directly generated. When the absolute value of the normal vector Z value is greater than or equal to the set normal vector Z threshold, the single slope roof house structure type is directly generated. When the number of normal vectors is determined to be equal to 2, area weight analysis processing is performed; when the number of normal vectors is determined to be greater than or equal to 3, normal vector normalization processing is performed; the normal vector Z threshold is 0.05;

[0269] Normalization processing of normal vectors: screen the new normal vector set in the same direction, remove the normal vectors with small polygonal areas, and simultaneously delete the single plane equation corresponding to the normal vector in the plane equation set to obtain a new normal vector set and plane equation set; when the number of new normal vectors is equal to 1, the flat roof house structure type is directly generated when the absolute value of the normal vector Z value is less than the set normal vector Z threshold, and the single slope roof house structure type is directly generated when the absolute value of the normal vector Z value is greater than or equal to the set threshold; when the number of new normal vectors is equal to 2 and the number of new normal vectors is greater than or equal to 3, process the normal vector results respectively;

[0270] Area weight analysis and processing generates house structure type: when the number of normal vectors is equal to 2, determine the relationship between the plane areas corresponding to the two normal vectors; when the difference between the plane areas corresponding to the two normal vectors is greater than or equal to the set area threshold, obtain the normal vector corresponding to the larger polygon area, and directly determine the normal vector Z value. When the absolute value of the normal vector Z value is less than the set normal vector Z threshold, directly generate the flat roof house structure type; when the absolute value of the normal vector Z value is greater than or equal to the set normal vector Z threshold, directly generate the single slope roof house structure type; when the difference between the plane areas corresponding to the two normal vectors is less than the set area threshold, calculate the difference product of the two new normal vectors. When the absolute value of the difference product Z value is less than the set normal vector Z threshold, directly generate the mansard roof house structure type; when the absolute value of the difference product Z value is greater than or equal to the set normal vector Z threshold, directly generate the single slope roof house structure type; the area threshold is 5km 2 ;

[0271] Normal vector comparison analysis generates house structure type: when the number of normalized normal vectors is greater than or equal to 3, they are sorted in descending order according to the polygon area value, and the normal vector with the largest polygon area is selected as the basic vector; when the difference between the current two areas is greater than or equal to 2 times the area threshold, the basic vector Z value is directly judged. When the absolute value of the basic vector Z is less than the set normal vector Z threshold, the flat roof house structure type is directly generated; when the absolute value of the basic vector Z value is greater than or equal to the set normal vector Z threshold, the single slope roof house structure type is directly generated; when the difference between the current two areas is less than 2 times the area threshold, when the absolute value of the basic vector Z value is greater than or equal to the set normal vector Z threshold, the difference product of other normal vectors and the basic vector is traversed, when the absolute value of the difference product Z value is less than the set normal vector Z threshold, the man-slope roof house structure type is directly generated, when the absolute value of the difference product Z value is greater than or equal to the set normal vector Z threshold, the single slope roof house structure type is directly generated; when the absolute value of the basic vector Z value is less than the set normal vector Z threshold, the flat roof house structure type is directly generated.

[0272] like Figure 6 As shown, a method for automatically identifying the roof structure of a building based on airborne point cloud, the specific steps of the house structure type verification are:

[0273] Classify the roof structure: by fitting the roof structure type of the raw building, assign 0 to the flat roof, 1 to the single-slope roof, and 2 to the gable roof, and record it in the corresponding building outline attribute;

[0274] Roof line extraction: Based on the segmented and grouped roof patch point cloud and building outline, obtain the seed point, traverse other points of the building roof patch, obtain the regional points and boundary points, and merge and fit to obtain the ridge line equation;

[0275] Judgment based on geometric features: when the ridge line equation cannot be fitted, the building structure type is a flat roof; when the ridge line equation is fitted, the point clouds on the left and right sides of the ridge line equation are judged. When there is a point on one side of the ridge line equation and the ridge line coincides with the edge of the roof, the building structure type is a single-slope roof; when there are points on both sides of the ridge line equation and the ridge line coincides with the edge of the roof, the building structure type is a human-slope roof;

[0276] Structural type modification: Compare and analyze the building structure type obtained with the corresponding building outline attribute structure type; when the two structure types are the same, directly determine that the structure type is correct; when the two structure types are different, add a building outline attribute field, fill in the structure type obtained by the verification method, and modify the structure type, such as Fig.12 Building outlines and structure types are shown.

[0277] like Figure 7 As shown, a method for automatically identifying the roof structure of a building based on airborne point cloud, the specific steps of extracting the roof ridge line are:

[0278] Seed point selection: Based on the segmented roof patch point cloud and building outline, the height range of the roof point cloud [h min ,h max ], select a height close to the highest point h of the roof patch max The point is taken as the seed point, where h min Represents the lowest point of the roof patch;

[0279] Get regional points: Set the initial growth region R, including the seed point. Traverse the remaining roof patch points and calculate the distance and normal vector angle between each point and the seed point; when the distance is less than or equal to the set distance threshold dth and the normal vector angle is greater than or equal to the set normal vector angle threshold cosθth, d(Pi,Pj)≤dth and ≥cosθth, place the point in the growth area R. When the distance is greater than the set distance threshold or the normal vector angle is greater than the set normal vector angle threshold, discard the point. The distance threshold is set to 30 meters and the normal vector angle threshold is set to 0.05.

[0280] Boundary point determination: When the growth area R cannot continue to increase, check the boundary point of R. Select the boundary point P b , so that P b Belongs to R, there is a point P nb does not belong to R, so that P b To P nb Distance L(P b ,P nb ) is less than or equal to d neighbor , its neighboris a distance slightly larger than dth, where d neighbor It is 1.2 times the set distance threshold.

[0281] Merge fitting: Get multiple growth areas (R1, R2, R3…R n ), calculate the adjacent growth area R i and R j The distance davg(R i ,R j ), determine davg(R i ,R j ) and the set merging threshold; when davg(R i ,R j ) is less than or equal to the merge threshold, then merge R i and R j , and fit the combined points to obtain the straight line equation of the ridge line. Its calculation formula is:

[0282]

[0283] n bi and n bj Represents R i and R j The number of boundary points,

[0284] and Represents a boundary point.

[0285] Glossary:

[0286] The LAS format is a binary file format designed specifically for 3D point cloud data, managed and maintained by the American Society for Photogrammetry and Remote Sensing. It is widely used in the fields of laser radar scanning and geographic information systems, and has become an industrial standard format for laser point cloud data.

[0287] GNSS, also known as the Global Navigation Satellite System, is an air-based radio navigation and positioning system that can provide users with all-weather three-dimensional coordinates, speed and time information at any location on the Earth's surface or in near-Earth space. It includes one or more satellite constellations and their augmentation systems required to support specific tasks, including China's Beidou Satellite Navigation System, the United States' Global Positioning System, Russia's GLONASS Satellite Navigation System and the European Union's Galileo Satellite Navigation System.

[0288] Projection zone UTM50: It is the 50th projection zone (Northern Hemisphere) of the Universal Transverse Mercator projection. It is a conformal transverse cylindrical projection. It divides the earth ellipsoid into multiple 6-degree wide longitude zones and performs conformal projection in each zone, thereby reducing the area deformation on the map. It is especially suitable for map making in mid-latitude areas.

[0289] The foregoing method descriptions and flow charts are provided only as illustrative examples and are not intended to require or imply that the above operations or steps of various aspects must be performed in the order given. As will be appreciated by those skilled in the art, the order of the boxes in the foregoing aspects may be performed in any order. Words such as "thereafter," "then," "next," and the like are not intended to limit the order of operations or steps; these words are only used to guide the reader through the description of the method. In addition, any singular reference to a claim element, for example, the use of the article "a," "an," or "the" is not to be construed as limiting the element to the singular.

[0290] The various illustrative logic blocks, modules, circuits, and algorithm steps described in conjunction with the aspects disclosed herein can all be implemented as electronic hardware, computer software, or a combination thereof. In order to clearly represent the interchangeability between hardware and software, the various illustrative components, blocks, modules, circuits, and steps are generally described above around their functions. Whether such functions are implemented as hardware or software depends on the specific application and the design constraints imposed on the entire method. A skilled person can implement the described functions in an alternative manner for each specific application, but such implementation decisions should not be interpreted as causing a departure from the scope of protection of the present invention.

[0291] The method in the embodiment of the present invention uses airborne radar to quickly obtain high-precision three-dimensional point clouds of urban building roof surfaces over a large area, and has the advantages of high precision and high penetration. At the same time, an improved algorithm can be used to segment building roof surfaces and fit roof plane equations. It is a high-precision and efficient data acquisition method for quickly obtaining urban three-dimensional model structure types.

[0292] The foregoing description of the disclosed aspects is provided to enable any person skilled in the art to implement or use the present invention. Various modifications to these aspects will be apparent to those skilled in the art, and the general principles defined herein may also be applied to other embodiments without departing from the spirit and scope of the present invention. Therefore, the present invention is not intended to be limited to the aspects provided herein, but rather to the widest scope consistent with the principles and novel features disclosed herein.

Claims

1. A method for automatically identifying building roof structures based on airborne point clouds, characterized in that The method comprises the following steps: laser point cloud preprocessing, building house extraction, segmentation and grouping of building roof surface laser point cloud, roof point cloud fitting plane equation set, multiple fitting planes generating roof structure type, house structure type verification, wherein: Verification of building structure type: Classify the building roof structure type and record it in the corresponding building outline attribute. Use the segmented and grouped roof patch point cloud and building outline as the data basis to obtain the ridge line. Predict the structure type based on the geometric features and ridge line. Determine the predicted structure type and compare it with the structure type in the corresponding building outline attribute. Modify the inconsistent structure type. Point cloud preprocessing: The original laser point cloud data is obtained through the mutual coordination and combined operation between the laser transmitter, inertial navigation and GNSS equipment. The laser point cloud is matched with the point cloud of a single flight and a single flight strip to obtain the registration point cloud of each flight. Then, the rigid registration point cloud of each flight is regionally optimized to obtain the optimized point cloud of the entire survey area. The optimized point cloud of the survey area is then converted from geodetic coordinates to the standard national 2000 coordinate system and the local city coordinate system to obtain the target point cloud of the survey area. The target point cloud of the survey area is converted into point cloud data in an internationally accepted format. Building extraction: remove noise points from the preprocessed point cloud, extract the ground point cloud, calculate the point cloud distance and normal vector, filter the low-value point cloud according to the distance and normal vector factors, delete the clustered tree point cloud around the building, and obtain the building roof point cloud data; Segmentation and grouping of building roof laser point cloud: Randomly select initial seed points and basic points from the building roof point cloud to construct a plane by normal vector and distance determination. For the points that cannot construct a plane, use them as other points to determine the distance and normal vector between the plane and determine the inner points and outer points. After grouping the inner points, perform a cycle of inner point selection on different planes and group them in order. Roof point cloud fitting plane equation set: Based on the outline of a single building, extract the laser point cloud of a single roof surface after segmentation and grouping. According to the grouping number of the laser point cloud of a single roof surface, select the laser point cloud of the roof surface in the same group as the plane fitting object, perform feature centering on the fitting object, establish the covariance matrix to obtain the eigenvalue and its corresponding eigenvector, and establish the plane equation set and normal vector set for the single fitting object after traversal processing; Generate roof structure type from multiple fitting planes: Take the building outline as the range, obtain the polygonal area of ​​the building roof patch, remove the plane equations and corresponding normal vectors with an area less than the threshold, obtain a new set of plane equations and normal vectors, and determine the number of normal vectors, perform normal vector normalization, area weight analysis, and normal vector comparison analysis to generate the house structure type; Verification of building structure type: classify the building roof structure type and record it in the corresponding building outline attribute. Use the segmented and grouped roof patch point cloud and building outline as the data basis to obtain the ridge line. Predict the structure type based on the geometric features and ridge line. Determine the predicted structure type and compare it with the structure type in the corresponding building outline attribute. Modify the inconsistent structure type.

2. The method for automatically identifying a building roof structure based on airborne point cloud according to claim 1, characterized in that: The specific steps of extracting the building house are as follows: Eliminate noise points: Remove extremely high and extremely low discrete points and noise points from the general format point cloud data to obtain a three-dimensional laser point cloud including topography, buildings and surrounding small pieces, trees, power lines, etc. within the target measurement area; Extract ground point cloud: Use opening and closing operations on the three-dimensional laser points within the target measurement area to preliminarily determine the ground points, remove isolated high points and fill small depressions to obtain a relatively smooth surface model, compare the point cloud with the constructed surface model, and reclassify the ground points and non-ground points; points that are close to the surface model are determined as ground points; points that are far away are determined as non-ground points. Update the surface model based on the new classification results, and repeat the above process until the classification results are stable or a certain number of iterations are reached, and obtain a three-dimensional laser point cloud containing ground points and non-ground points; Calculate the point cloud distance and normal vector: Calculate the distance from each point to the ground point directly below and assign it to the point cloud attribute; the calculation formula is: Ax+By+Cz+D=0 (A, B, C) represents the normal vector of the plane formed by nearby ground points, D represents the distance from the origin of the coordinate axis (0,0,0) to the plane. (x, y, z) represents the coordinates of nearby ground points in the three coordinate axis directions. (x0, y0, z0) represents the coordinates of the point cloud in the three coordinate axis directions; d represents the distance from the point to the ground point directly below; Then, using the nearby ground points as reference and the neighborhood points as support, we select 1-5 times the average distance between the nearby ground points as the neighborhood radius value, calculate the normal vector from each point in the point cloud to the nearby ground points, and assign it to the point attributes; the calculation formula is: Ax+By+Cz+D=0 Ax+By+Cz+D=0 represents the plane equation of the ground point, (A, B, C) represents the normal vector of the plane formed by nearby ground points, D represents the distance from the origin of the coordinate axis (0,0,0) to the plane. (x1, y1, z1) represents the coordinates of any point in the three coordinate axes, namely point Q. (x, y, z) represents the coordinates of any point on the ground in the three coordinate axis directions, namely point P. represent exist projection on Filter low vegetation point cloud: set the distance threshold by the minimum height difference between the building roof and the nearby ground, and obtain the distance value to the nearby ground point by traversing each point in the laser point cloud. Use the distance threshold as a reference to remove the laser point cloud of low vegetation, small objects close to the ground, and wall facades to obtain the processed laser point cloud; the reference range of the distance value is between 0.5 meters and 2 meters; Eliminate irrelevant point clouds around buildings: set the density radius to calculate the number of points around each point in the point cloud, obtain the density of the point, determine the relationship between the density of the point and the set core point threshold, and determine it as a core point when the number of points is greater than or equal to the set core point threshold. When the density of the point is less than the set core point threshold, the point is within a certain core point density radius and is determined as a boundary point to form a cluster. The cluster category is determined according to the shape and density of the trees, and the clusters belonging to the tree point cloud are eliminated to obtain the building roof point cloud data. The reference range of the density radius is between 0.5 meters and 3 meters.

3. The method for automatically identifying a building roof structure based on airborne point cloud according to claim 1, characterized in that: The specific steps of segmenting and grouping the building roof surface laser point cloud are: Select the initial seed point and base point: randomly select a point from the building roof point cloud as the initial seed point, randomly select another two points from the remaining point cloud as the base points of the plane, and add them to the current plane model; Judge the normal vector and distance and construct the plane: judge the distance and normal vector consistency between the base point and the seed point. When the distance is less than or equal to the set base distance threshold and the normal vector consistency, build the plane. When the distance is greater than the set base distance threshold or the normal vector is inconsistent, discard the base point and reselect the base point and calculate the distance and normal vector consistency between the new base point and the seed point. Loop until the base point and the seed point are all constructed with the plane, and cache the points that cannot construct the plane as other points. The reference range for setting the basic distance threshold is between 10 meters and 100 meters, and the calculation formula for constructing the plane equation is: Ax+By+Cz+D=0 Represents the normal vector of the plane formed by the top surface of the building roof, that is, (x, y, z) represents the coordinates of nearby ground points in the three coordinate axis directions. (x4, y4, z4) represents the coordinates of the seed point in the three coordinate axis directions. (x2, y2, z2) and (x3, y3, z3) represent the coordinates of two base points in the three coordinate axis directions. D represents the distance from the origin of the coordinate axis (0,0,0) to the plane; Get the inner points of the roof patch: Analyze and obtain the point cloud of the target object surface. The surface point cloud has up and down fluctuations. Set the fluctuation threshold, traverse the other points in the laser point cloud except the seed point and the basic point, obtain the distance between each other point and the constructed plane model, and determine the consistency between the normal vector of the other point and the normal vector of the constructed plane. When the distance is less than the fluctuation threshold and the normal vector of the other point is determined to be consistent with the normal vector of the constructed plane, the other point is classified as an inner point. When the distance is greater than or equal to the fluctuation threshold or the normal vector of the other point is determined to be inconsistent with the normal vector of the constructed plane, the other point is classified as an outer point. The reference range of the fluctuation threshold is between 0.10 meters and 0.30 meters. Count the number of inliers that meet the conditions: Set a threshold for the number of inliers based on the minimum number of point clouds required for the plane model, and count the relationship between the number of inliers and the threshold. When the number of inliers is less than the threshold, all inliers are classified as non-planar discrete points. When the number of inliers is greater than or equal to the threshold, all inliers are classified as the same plane points to form the same plane point cloud. The same plane point clouds are grouped and assigned group numbers according to the order of judgment and screening. Loop to select points in different planes: remove the remaining point cloud of the same plane points that have been assigned numbers, start from the step of selecting initial seed points and basic points to the step of counting the number of qualified internal points, select points in different planes and assign group numbers in the order of selection, loop this step until the remaining points cannot construct a plane, and end.

4. The method for automatically identifying a building roof structure based on airborne point cloud according to claim 1, characterized in that: The specific steps of fitting the roof point cloud to the plane equation set are: Extracting laser point cloud of single roof patch: Based on the outline of a single building, extract the segmented point cloud of the building roof, and select the laser point cloud of the roof patch in the same group as the plane fitting object set according to the group number of the laser point cloud of the single roof patch. When the number in the attribute of the fitting object set is equal to 1, directly process the fitting object; when the number in the attribute of the fitting object set is greater than 1, process them as fitting objects in the order of the numbers. Feature centering: After receiving a single fitting object, subtract the mean value M = (x5, y5, z5) of the point cloud in each dimension from the coordinate value of each dimension of the fitting object to obtain the feature-centered fitting object P. The mean value of each dimension after the transformation is zero. The calculation formula is: x i ,y i ,z i Represents the coordinates of the centralized point cloud; (x5, y5, z5) represents the mean of the point cloud in this dimension, i.e., M; p represents the feature-centered fitting object; Establish the covariance matrix: calculate the feature-centered fitting object, calculate the covariance matrix, and obtain its eigenvalues ​​λ1, λ2, λ3 and eigenvectors V1, V2, V3. The calculation formula is: n represents the number of points in the fitting object, C represents the covariance matrix of the three coordinates, cov(x,y) represents the covariance of the x-coordinate and the y-coordinate. cov(x,y) represents the variance of the x coordinate, x i represents the x-coordinate of the i-th point in the fitting object, y i Represents the y coordinate of the i-th point in the fitting object, z i Represents the z coordinate of the i-th point in the fitting object; Establish a set of plane equations and a set of normal vectors: select the eigenvector λ3 corresponding to the minimum eigenvalue in the covariance matrix to determine the plane equation; the eigenvector V3 corresponding to the minimum eigenvalue λ3 is the direction perpendicular to the fitting plane, V3 = (a, b, c) is the plane normal vector, the plane normal vector (a, b, c,) and the average position M of the point cloud determine the D1 value, D1 = -a*x5-b*y5-c*z5, and the final plane equation of the roof patch is: ax+by+cz+D1=0; the plane equations and normal vectors of the single fitting object after the traversal cycle are integrated to form a one-to-one corresponding set of plane equations and normal vectors.

5. The method for automatically identifying a building roof structure based on airborne point cloud according to claim 1, characterized in that: The specific steps of generating the roof structure type by using multiple fitting planes are as follows: Calculate the area of ​​the building roof patch: Take the building outline as the range, combine the plane equation set of the building roof patch and the corresponding numbered point cloud, obtain the coordinates of the corner points of each plane, and calculate the polygonal area of ​​the building roof patch. The formula is: (x i ,y i ,z i ) represents the coordinates of the i-th point on the plane; Represents the difference product of the vectors of the three corner points of a single plane. (n1, n2, n3) represents the coordinates of the difference product of the vectors of the three corner points of a single plane in the directions of the three coordinate axes; n represents the number of corner points of a single plane, Represents the plane area composed of the first point, the i-th point, and the i+1-th point of a single plane. S represents the area of ​​a single polygon; Construct a valid building roof normal vector set: According to the polygonal area of ​​a single plane corresponding to the normal vector set, obtain the area mean and median, select the minimum of the area mean and median as the effective area threshold, remove the plane equations and corresponding normal vectors whose areas are less than the effective area threshold, and obtain a new plane equation set and normal vector set; Determine the number of normal vectors: Count the number of normal vectors in the new normal vector set and perform determination processing. When it is determined that the number of normal vectors is equal to 1 and the absolute value of the normal vector Z value is less than the set normal vector Z threshold, directly generate a flat roof house structure type; when the absolute value of the normal vector Z value is greater than or equal to the set normal vector Z threshold, directly generate a single slope roof house structure type; when it is determined that the number of normal vectors is equal to 2, perform area weight analysis processing; when it is determined that the number of normal vectors is greater than or equal to 3, perform normal vector normalization processing; The reference range of the normal vector Z threshold is between 0 meters and 0.1 meters; Normal vector normalization: filter the new normal vector set in the same direction, remove the normal vectors with small polygon area, and simultaneously delete the single plane equation corresponding to the normal vector in the plane equation set to obtain a new normal vector set and plane equation set; When the number of new normal vectors is equal to 1, the flat roof house structure type is directly generated when the absolute value of the normal vector Z value is less than the set normal vector Z threshold, and the single slope roof house structure type is directly generated when the absolute value of the normal vector Z value is greater than or equal to the set threshold; when the number of new normal vectors is equal to 2 and the number of new normal vectors is greater than or equal to 3, the normal vector results are processed respectively; Area weight analysis and processing to generate house structure type: when the number of normal vectors is equal to 2, determine the relationship between the plane areas corresponding to the two normal vectors; When the difference in the plane areas corresponding to the two normal vectors is greater than or equal to the set area threshold, the normal vector corresponding to the larger polygon area is obtained, and the normal vector Z value is directly determined. When the absolute value of the normal vector Z value is less than the set normal vector Z threshold, the flat roof house structure type is directly generated. When the absolute value of the normal vector Z value is greater than or equal to the set normal vector Z threshold, the single slope roof house structure type is directly generated; when the difference in the plane areas corresponding to the two normal vectors is less than the set area threshold, the difference product of the two new normal vectors is calculated. When the absolute value of the difference product Z value is less than the set normal vector Z threshold, the human-slope roof house structure type is directly generated. When the absolute value of the difference product Z value is greater than or equal to the set normal vector Z threshold, the single slope roof house structure type is directly generated; the reference range value of the area threshold is 3km 2 -5km 2 ; Normal vector comparison analysis generates house structure type: when the number of normal vectors after normalization is greater than or equal to 3, they are sorted in descending order according to the polygon area value, and the normal vector with the largest polygon area is selected as the basic vector; when the difference between the current two areas is greater than or equal to 2 times the area threshold, the basic vector Z value is directly determined. When the absolute value of the basic vector Z value is less than the set normal vector Z threshold, the flat roof house structure type is directly generated; when the absolute value of the basic vector Z value is greater than or equal to the set normal vector Z threshold, the single slope roof house structure type is directly generated; when the difference between the current two areas is less than 2 times the area threshold, when the absolute value of the basic vector Z value is greater than or equal to the set normal vector Z threshold, the difference product of other normal vectors and the basic vector is traversed, when the absolute value of the difference product Z value is less than the set normal vector Z threshold, the human slope roof house structure type is directly generated, when the absolute value of the difference product Z value is greater than or equal to the set normal vector Z threshold, the single slope roof house structure type is directly generated; when the absolute value of the basic vector Z value is less than the set normal vector Z threshold, the flat roof house structure type is directly generated.

6. The method for automatically identifying a building roof structure based on airborne point cloud according to claim 1, characterized in that: The specific steps of the house structure type verification are: Classify the roof structure: by fitting the roof structure type of the raw building, assign 0 to the flat roof, 1 to the single-slope roof, and 2 to the gable roof, and record it in the corresponding building outline attribute; Roof line extraction: Based on the segmented and grouped roof patch point cloud and building outline, obtain the seed point, traverse other points of the building roof patch, obtain the regional points and boundary points, and merge and fit to obtain the ridge line equation; Judgment based on geometric features: when the ridge line equation cannot be fitted, the building structure type is a flat roof; when the ridge line equation is fitted, the point clouds on the left and right sides of the ridge line equation are judged. When there is a point on one side of the ridge line equation and the ridge line coincides with the edge of the roof, the building structure type is a single-slope roof; when there are points on both sides of the ridge line equation and the ridge line coincides with the edge of the roof, the building structure type is a human-slope roof; Structural type modification: Compare and analyze the building structure type and the corresponding building outline attribute structure type; when the two structure types are the same, directly determine that the structure type is correct; When two structure types are different, add a building outline attribute field, fill in the structure type obtained through the verification method, and modify the structure type.

7. The method for automatically identifying a building roof structure based on airborne point cloud according to claim 6, characterized in that: The specific steps of extracting the ridge line are: Seed point selection: Based on the segmented roof patch point cloud and building outline, the height range of the roof point cloud [h min ,h max ], select a height close to the highest point h of the roof patch max The point is taken as the seed point, where h min Represents the lowest point of the roof patch; Get regional points: Set the initial growth region R, including the seed point. Traverse the remaining roof patch points and calculate the distance and normal vector angle between each point and the seed point; when the distance is less than or equal to the set distance threshold dth and the normal vector angle is greater than or equal to the set normal vector angle threshold cosθth, d(P i ,P j )≤dth and The point is placed in the growth region R. When the distance is greater than the set distance threshold or the normal vector angle is greater than the set normal vector angle threshold, the point is discarded. Boundary point determination: When the growth area R cannot continue to increase, check the boundary point of R. Select the boundary point P b , so that P b Belongs to R, there is a point P nb does not belong to R, so that P b To P nb Distance L(P b ,P nb ) is less than or equal to d neighbor , its neighbor is a distance slightly larger than dth. Merge fitting: Get multiple growth areas (R1, R2, R3…R n ), calculate the adjacent growth area R i and R j The distance davg(R i ,R j ), determine davg(R i ,R j ) and the relationship between the set merging threshold; when davg(R i ,R j ) is less than or equal to the merge threshold, then merge R i and R j , and fit the combined points to obtain the straight line equation of the ridge line. Its calculation formula is: n bi and n bj Represents R i and R j The number of boundary points, and Represents a boundary point.

Images