Airborne LiDAR Data Building Extraction Method Based on Adaptive Local Spatial-Spectral Consistency

Abstract

The present invention proposes a method for extracting buildings from airborne LiDAR data under adaptive local spatial-spectral consistency, which relates to the technical field of remote sensing data processing. The method first eliminates abnormal data from the original airborne LiDAR point cloud data and regularizes it into an intensity voxel dataset; secondly, seed voxels are selected based on the elevation jump and edge straight-line characteristics of buildings, and clustering based on spatial connectivity is performed on them to establish a spectral consistency measure; then, building roof extraction and optimization are carried out under the constraint of adaptive local spatial-spectral consistency; finally, the building facade is extracted by combining the spatial constraint of the roof edge and local intensity consistency. The method proposes a construction scheme for the intensity voxel model of airborne LiDAR data and a building extraction scheme for complex urban scenarios under adaptive local spatial-spectral consistency on this basis, meeting the requirements of the diversity of building shapes, materials, etc. in complex urban scenarios, and contributing to the development of airborne LiDAR point cloud data processing and applications based on the intensity voxel model theory.

Description

Technical Field

[0001] The present invention relates to the technical field of remote sensing data processing, and particularly to a method for extracting buildings from airborne LiDAR data under adaptive local spatial-spectral consistency. Background Art

[0002] With the gradual definition and deepening of the concept of real-scene three-dimensional city, this has become a hot research topic in the field of surveying and mapping remote sensing. Among them, urban construction, as the core content of real-scene three-dimensional city construction, naturally cannot be ignored. Therefore, the automatic, high-precision and rapid extraction of building targets has naturally become a research hotspot. Airborne Light Detection And Ranging (LiDAR) technology can provide dense, accurate, and georeferenced true three-dimensional (3D) point cloud data, and contains the intensity information of echo signals. Therefore, airborne LiDAR data is particularly suitable for 3D target extraction. However, due to the limitations of complex building shapes and irregular point cloud distributions in urban areas, there are still technical difficulties in the research of extracting urban buildings using airborne LiDAR point cloud data. Currently, there are the following deficiencies: First, in terms of data structure, airborne LiDAR point cloud data has the characteristics of being discrete and disordered, and it is difficult to express the topological relationship between points. When 3D point cloud data is mapped into a raster image, data information loss will occur. Second, in existing building extraction research, only specific types of relatively regular buildings can be extracted, and the robustness for extracting polymorphic buildings in urban areas is insufficient. Third, there are various building materials in urban areas. When using intensity information to assist in building extraction, it is often based on global statistical characteristics, which is not suitable for building extraction in complex urban scenarios. Summary of the Invention

[0003] The technical problem to be solved by the present invention is to provide a method for extracting buildings from airborne LiDAR data under adaptive local spatial-spectral consistency in view of the above-mentioned deficiencies of the prior art. The present invention combines the voxel data structure with building target extraction. The voxel data structure can clearly express the spatial structure of buildings and is a true 3D data structure. Using it to express airborne LiDAR point cloud data will not cause information loss, and at the same time, it avoids the problem of misclassification of buildings caused by factors such as same-spectrum different objects and same-object different spectra, and can meet the requirements of the diversity of building shapes, materials, etc. in complex urban scenarios.

[0004] In order to achieve the above object of the present invention, the technical solution adopted by the method for extracting buildings from airborne LiDAR data under adaptive local spatial-spectral consistency of the present invention is as follows:

[0005] Step 1: Read the original airborne LiDAR point cloud data to form an original airborne LiDAR point cloud data set;

[0006] Step 2 regularizes the original airborne LiDAR point cloud dataset into an intensity voxel dataset;

[0007] The steps specifically include:

[0008] Step 2.1 removes abnormal data from the original airborne LiDAR point cloud dataset to obtain a dataset with abnormal data removed;

[0009] The steps specifically include:

[0010] Step 2.1.1 renders and displays the original airborne LiDAR point cloud dataset by elevation, visually determines the true elevation range of the ground objects from the side view, and sets the highest elevation threshold T h , the lowest elevation threshold T l ;

[0011] Step 2.1.2 for each laser point in the original airborne LiDAR point cloud dataset, if its elevation value is higher than the highest elevation threshold T h or lower than the lowest elevation threshold T l , then this laser point is abnormal data and is removed, otherwise this laser point is retained, and finally a dataset with elevation anomalies removed is obtained;

[0012] Step 2.1.3 for the dataset with elevation anomalies removed obtained in the above steps, counts the frequency of the intensity values of each laser point therein, and visually determines the intensity threshold corresponding to the true ground object and records it as I d ;

[0013] Step 2.1.4 determines the laser points with intensity values higher than I d in the dataset with elevation anomalies removed as intensity abnormal data and removes them, and finally obtains a dataset with abnormal data removed.

[0014] Step 2.2 regularizes the dataset with abnormal data removed into an intensity voxel dataset;

[0015] The steps specifically include:

[0016] Step 2.2.1 represents the three-dimensional space range with an axis-aligned bounding box of the dataset with abnormal data removed;

[0017] Step 2.2.2 determines the resolution (Δx, Δy, Δz) of the voxels in the x, y, and z directions according to the average point spacing of the laser points in the dataset with abnormal data removed;

[0018] The calculation formulas for the resolution Δx, Δy, and Δz of the voxels in the x, y, and z directions are as follows:

[0019]

[0020] where, S xyIt is a two-dimensional point set obtained by projecting the abnormal data set onto the XOY plane. C(S xy ) is the convex hull of the point set S xy , and A(C(S xy )) is the area of the convex hull C(S xy ).

[0021] Step 2.2.3 divides the axis-aligned bounding box according to the voxel resolution (Δx, Δy, Δz) to obtain a 3D voxel grid, and each 3D voxel grid unit is a voxel;

[0022] The division of the axis-aligned bounding box into a 3D voxel grid based on the voxel resolution (Δx, Δy, Δz) is represented by a 3D voxel array. As shown below:

[0023] V = {v j (r j , c j , l j ), j = 1, …, m}, (2)

[0024] where V is the set of voxels in the 3D voxel array, j is the voxel index; m is the number of voxels; v j is the voxel value of the j-th voxel; (r j , c j , l j ) is the coordinate of the j-th voxel in the voxel array, r j is the row number, c j is the column number, and l j is the layer number.

[0025] The number of voxels in the X direction is R, the number of voxels in the Y direction is C, and the number of voxels in the Z direction is L.

[0026] The calculation formulas for R, C, and L are as follows:

[0027]

[0028] where is the ceiling operator, xmax = max{xi', i' = 1, ..., t}, x min = min{x i' , i' = 1, ..., t}, y max = max{y i' , i' = 1, ..., t}, y min = min{y i' , i' = 1, ..., t}, i' is the index of the data in the abnormal data set removed, t is the number of data in the abnormal data set removed, and the coordinate of the i'-th data removed from the abnormal data set is (x i' , yi' , z i' ).

[0029] The calculation formula for the number of voxels m is as follows:

[0030] m = R * C * L (4)

[0031] Step 2.2.4 Map each laser point in the abnormal data removed dataset to the 3D voxel grid, and then assign values to each voxel according to the intensity attributes of the laser points contained in the 3D voxel grid to obtain the intensity voxel dataset.

[0032] The specific process of assigning values to each voxel according to the intensity attributes of the laser points contained in the 3D voxel grid is as follows:

[0033] Assign the voxel containing the laser point as the average value of the laser point intensity, assign the voxel without the laser point as -1, and further discretize the non - negative voxel values to {0, …, 255} to obtain each voxel value.

[0034] Step 3 Search for the seed voxels of the building monomers and establish their spectral consistency measures;

[0035] The steps specifically include:

[0036] Step 3.1 Search for the building seed voxels in the non - negative intensity voxel data based on the elevation jump and edge straight - line characteristics of the building;

[0037] Step 3.2 Cluster the building seed voxels under spatial connectivity constraints, and use the clustering result as the seed voxels of each building monomer;

[0038] The steps specifically include:

[0039] Step 3.2.1 Traverse all the seed voxels that are spatially connected to the k - th seed voxel using the depth - first strategy according to the spatial connectivity characteristics and mark them as L l , where l is the index of the marking label, l = 1, 2,...;

[0040] Step 3.2.2 Continue to scan the unmarked seed voxels in the intensity voxel dataset until all the seed voxels are marked, obtaining several 3D connected regions composed of seed voxels, and use the seed voxels marked with the same label as the seed voxels of the same building monomer.

[0041] Step 3.3 Statistically analyze the spectral characteristics of the clustered seed voxels and use them as the spectral consistency measures adapted to each building monomer;

[0042] Step 4 Extract and optimize the building roofs;

[0043] The steps specifically include:

[0044] Step 4.1 Extraction of 3D connected regions of building roofs under adaptive local spatial-spectral consistency constraints;

[0045] The steps specifically include:

[0046] Step 4.1.1 Initialize an empty stack, and put the set of seed voxels labeled with label L that have not been put into the stack for traversal into the stack; l and the seed voxels into the stack;

[0047] Step 4.1.2 Pop the top element from the stack, and search for voxels that are spatially connected and locally intensity-consistent with it in the non-negative voxel data of the intensity voxel dataset, label them with the same label L as this top element of the stack l and put them into the stack;

[0048] The voxels with locally consistent intensity are two voxels whose voxel value difference is less than the intensity difference threshold.

[0049] Step 4.1.3 If the stack is empty, terminate the program; otherwise, go to Step 4.1.2 until all seed voxels have been put into the stack for traversal, and thus obtain several 3D connected regions of building roofs.

[0050] Step 4.2 Optimization of the extraction results of 3D connected regions of building roofs based on building area, local normal vector consistency, and density characteristics;

[0051] The steps specifically include:

[0052] Step 4.2.1 Remove 3D connected regions that are not building roofs based on area characteristics;

[0053] The removal condition is: the minimum building area in the original airborne LiDAR point cloud dataset is A min , and the maximum building area in the original airborne LiDAR point cloud dataset is A max , for any 3D connected region, if its horizontal projected area is greater than or equal to A min and less than or equal to A max , then this 3D connected region is determined to be a building roof and is retained, otherwise this 3D connected region is removed.

[0054] The building area in the original airborne LiDAR point cloud dataset refers to the projected area of the building on the horizontal plane.

[0055] Step 4.2.2 Remove 3D connected regions that are not building roofs based on local normal vector consistency;

[0056] For any 3D connected region, all voxels have normal vector values calculated according to the principal component analysis algorithm (PCA). Based on the principle of consistent local normal vectors, the normal vector angle threshold is T. θ If the angle between the normal vector of a single voxel and its neighboring voxels within this region is within the threshold, then this voxel is retained; otherwise, it is deleted. Finally, if the number of voxels with consistent normal vector angles remaining in this 3D connected region exceeds a set proportion of the total number of voxels in this 3D connected region, then this region is determined to be a building roof and is retained; otherwise, this 3D connected region is excluded.

[0057] Step 4.2.3 Exclude 3D connected regions that are not building roofs based on density characteristics.

[0058] Set the density threshold to T. d For any 3D connected region, if its density is greater than the given density threshold T, d then this 3D connected region is determined to be a building roof and is retained; otherwise, this 3D connected region is excluded.

[0059] Step 5 Extract building facades by combining the edge space constraints of building roofs and local intensity consistency.

[0060] The specific steps include:

[0061] Step 5.1 Project the optimized 3D connected regions of each building roof onto the XY plane, and then detect the edge contour voxels of each building roof.

[0062] Step 5.2 On the horizontal plane, with the edge contour voxels of a single building roof as the center, set the edge space constraint conditions of the building roof with a width of one voxel on both the inner and outer sides.

[0063] For any non - negative voxel located under the edge space constraints of the building roof, if the difference between its reflected intensity value and the mean intensity of the building edge voxels is within plus or minus two standard deviations, then this voxel is determined to be a building facade voxel.

[0064] The beneficial effects produced by adopting the above - mentioned technical solutions are as follows:

[0065] The present invention proposes a method for extracting buildings from airborne LiDAR data under adaptive local spatial - spectral consistency. This method proposes a construction scheme for the intensity voxel model of airborne LiDAR data and a building extraction scheme for complex urban scenes under adaptive local spatial - spectral consistency, avoiding the problem of misclassification of buildings caused by factors such as same - spectrum different objects and same - object different spectra, meeting the requirements of the diversity of building shapes, materials, etc. in complex urban scenes, and contributing to the development of airborne LiDAR point cloud data processing and applications based on the intensity voxel model theory. Brief Description of the Drawings

[0066] Figure 1 It is a flowchart of the method for extracting buildings from airborne LiDAR data under adaptive local spatial-spectral consistency in the specific embodiment of the present invention;

[0067] Figure 2 It is the original airborne LiDAR point cloud data in the specific embodiment of the present invention;

[0068] Among them, (a) is the point cloud data of Area2; (b) is the point cloud data of Area3; (c) is the image corresponding to the point cloud data of Area2; (d) is the image corresponding to the point cloud data of Area3;

[0069] Figure 3 It is a flowchart of regularizing the original airborne LiDAR point cloud data into an intensity voxel dataset in the specific embodiment of the present invention;

[0070] Figure 4 It is a top view of the intensity voxel data obtained by regularizing the LiDAR point cloud data in the specific embodiment of the present invention;

[0071] Among them, (a) is the voxel data corresponding to Area2; (b) is the voxel data corresponding to Area3;

[0072] Figure 5 It is a flowchart of searching for seed voxels of building monomers and establishing their spectral consistency measures in the specific embodiment of the present invention;

[0073] Figure 6 It is the LABEL(k,u,L l ) program process in the specific embodiment of the present invention;

[0074] Figure 7 It is a top view of the seed clustering result in the specific embodiment of the present invention;

[0075] Among them, (a) is the seed clustering result of Area2; (b) is the seed clustering result of Area3;

[0076] Figure 8 It is a flowchart of extracting and optimizing the building roof under the constraint of adaptive local spatial-spectral consistency in the specific embodiment of the present invention;

[0077] Figure 9 It is the LABEL2(u,L l ) program process in the specific embodiment of the present invention;

[0078] Figure 10 It is a top view of the 3D connected region construction result of the intensity voxel data in the specific embodiment of the present invention;

[0079] Among them, (a) is the clustering result of Area2; (b) is the clustering result of Area3;

[0080] Figure 11 This is the flowchart for extracting the building facade in the specific embodiment of the present invention;

[0081] Figure 12 This is the schematic diagram for establishing the spatial constraint of the building roof edge in the specific embodiment of the present invention;

[0082] Figure 13 This is the building extraction result of the airborne LiDAR data under the adaptive local spatial-spectral consistency in the specific embodiment of the present invention;

[0083] Among them, (a) is the extraction result of Area2; (b) is the extraction result of Area3. Specific embodiment

[0084] Next, in combination with the drawings and embodiments, the specific embodiments of the present invention will be further described in detail. The following embodiments are used to illustrate the present invention, but are not used to limit the scope of the present invention.

[0085] In this embodiment, this method is programmed and implemented on the Intel(R) Core(TM) i7-7700 CPU@3.60GHz, with 32GB of memory and the Windows 10 Enterprise Edition system, using the MATLAB(R2016b) 9.1.0 platform, and the effectiveness of the method is further verified through the accuracy evaluation of the method.

[0086] A method for extracting buildings from airborne LiDAR data under adaptive local spatial-spectral consistency, as Figure 1 shown, includes the following steps:

[0087] Step 1: Read the original airborne LiDAR point cloud data to form the original airborne LiDAR point cloud data set;

[0088] In this embodiment, two sets of urban sample data (Area2 and Area3) provided by the International Society for Photogrammetry and Remote Sensing (ISPRS) Working Group III / 4, which are specifically used for testing target classification algorithms, as Figure 2 shown, are used as experimental data to test the effectiveness and feasibility of the method. The experimental data is acquired by a Leica ALS50 airborne laser scanner (flight altitude 500 meters, field of view angle 45°). These two sets of data contain high-rise urban residential buildings surrounded by trees and residential areas with small attached buildings. The density of the point cloud data is 4 points / m2 。

[0089] In this embodiment, the original airborne LiDAR point cloud dataset P = {p i (x i , y i , z i ), i = 1,..., n} is defined, where i is the index of the original airborne LiDAR point cloud data, n is the number of the original airborne LiDAR point cloud data, and p i is the i-th original airborne LiDAR point cloud data, and its coordinates are (x i , y i , z i ).

[0090] Step 2 Regularize the original airborne LiDAR point cloud dataset into an intensity voxel dataset, and the specific process is as Figure 3 shown.

[0091] Step 2.1 Remove the abnormal data from the original airborne LiDAR point cloud dataset to obtain a dataset with abnormal data removed;

[0092] The specific steps include:

[0093] Step 2.1.1 Render and display the original airborne LiDAR point cloud dataset by elevation, visually judge the true elevation range of the ground objects from the side view, and set the highest elevation threshold T h and the lowest elevation threshold T l ;

[0094] Step 2.1.2 For each laser point in the original airborne LiDAR point cloud data, if its elevation value is higher than the highest elevation threshold T h or lower than the lowest elevation threshold T l , then this laser point is abnormal data and is removed, otherwise this laser point is retained, and finally a dataset with elevation abnormal data removed is obtained;

[0095] Step 2.1.3 For the dataset with elevation abnormal data removed obtained in the above step, count the frequency of the intensity values of each laser point therein, and visually display the statistical results in the form of a histogram, and visually determine the intensity threshold corresponding to the true ground object and denote it as I d ;

[0096] Step 2.1.4 Determine the laser points with intensity values higher than I d in the dataset with elevation abnormal data removed as intensity abnormal data and remove them, and finally obtain a dataset with abnormal data removed.

[0097] In this embodiment, the dataset with abnormal data removed is denoted as Q = {qi ' (x i' , yi ', z i' ), i' = 1, ..., t}, where i' is the index of the data removed from the abnormal data set, t is the number of data removed from the abnormal data set, q i' is the i'-th data in the abnormal data set removed, and its coordinates are (x i' , yi ' , z i' ).

[0098] In this embodiment, the highest elevation threshold T h and the lowest elevation threshold T l are constants, and their values need to be determined according to the spatial distribution of the original airborne LiDAR point cloud data. The intensity threshold I d is a constant, and its value is determined according to the intensity frequency distribution of the original airborne LiDAR point cloud data.

[0099] Step 2.2 Regularize the abnormal data set removed into an intensity voxel data set.

[0100] Step 2.2.1 Represent the three-dimensional space range with an axis-aligned bounding box of the abnormal data set removed.

[0101] In this embodiment, the axis-aligned bounding box of the abnormal data set removed Q is a cuboid, and its bottom size can be determined by finding the minimum circumscribed rectangle of the projection of the abnormal data set removed Q on the XY plane, and its height can be determined by (z max - z min ).

[0102] Among them, z max is the maximum value of the z coordinates of the laser points in the abnormal data set removed Q, z min is the minimum value of the z coordinates of the laser points in the abnormal data set removed Q, z max = max{z i' , i' = 1, ..., t}, z min = min{z i' , i' = 1, ..., t}.

[0103] Step 2.2.2 Determine the resolution (Δx, Δy, Δz) of the voxels in the x, y, and z directions according to the average point spacing of the laser points in the abnormal data set removed Q.

[0104] In this embodiment, the calculation formulas for the resolution Δx, Δy, and Δz of the voxels in the x, y, and z directions are shown in Equation (5):

[0105]

[0106] Among them, S xy = {(x i ', y i'), i' = 1, ..., t} is the two-dimensional point set obtained by projecting the abnormal data set Q onto the XOY plane, and C(S xy ) is the point set S xy 's convex hull, and A(C(S xy )) is the area of the convex hull C(S xy ).

[0107] Step 2.2.3 divides the axis-aligned bounding box according to the voxel resolution (Δx, Δy, Δz) to obtain a 3D voxel grid, and each 3D voxel grid unit is a voxel.

[0108] In this embodiment, based on the voxel resolution (Δx, Δy, Δz), the axis-aligned bounding box can be divided into a 3D voxel grid, which is represented by a 3D voxel array. Let V be the set of voxels in the 3D voxel array, as shown in Equation (6):

[0109] V = {v j (r j , c j , l j ), j = 1, …, m}, (6)

[0110] where j is the voxel index; m is the number of voxels; v j is the voxel value of the j-th voxel; (r j , c j , l j ) is the coordinate (row, column, and layer number) of the j-th voxel in the voxel array. The number of voxels in the X direction is R, the number of voxels in the Y direction is C, and the number of voxels in the Z direction is L. Among them, R, C, and L are as shown in Equation (7):

[0111]

[0112] where is the ceiling operator, x max = max{x i' , i' = 1, ..., t}, x min = min{x i' , i' = 1, ..., t}, y max = max{y i' , i' = 1, ..., t}, y min = min{y i' , i' = 1, ..., t}.

[0113] From this, it can be obtained that the number of voxels m is as shown in Equation (8):

[0114] m = R * C * L (8)

[0115] Step 2.2.4 Map each laser point in the abnormal data set removed to the 3D voxel grid, and then assign values to each voxel according to the intensity attribute of the laser points contained in the 3D voxel grid to obtain an intensity voxel data set.

[0116] The specific process of assigning values to each voxel according to the intensity attribute of the laser points contained in the 3D voxel grid is as follows:

[0117] Assign the voxel containing the laser point to the average value of the laser point intensity, assign the voxel without the laser point to -1, and further discretize the non-negative voxel values to {0, …, 255} to obtain each voxel value.

[0118] In this embodiment, map each laser point in the abnormal data set Q removed to the 3D voxel grid, and then assign values to each voxel according to the average intensity of the laser points contained in the 3D voxel grid, and assign the voxel without the laser point to -1, as shown in Equation (9):

[0119]

[0120] Among them, is the floor function operator. Further discretize the non-negative voxel values to {0, …, 255} to obtain each voxel value. Thus, an intensity voxel data set is obtained, and the regularization operation of the abnormal data set removed is completed.

[0121] In this embodiment, the top view of the voxel data obtained by point cloud regularization is as Figure 4 shown. The intensity value represents different voxel intensity values. It can be seen that the intensity values of the building roof are relatively consistent and different from those of other targets.

[0122] Step 3 Search for the seed voxels of building monomers and establish their spectral consistency measures. The specific process is as Figure 5 shown;

[0123] Step 3.1 Search for the building seed voxels in the non-negative intensity voxel data according to the elevation jump and edge straight line characteristics of the building, and scan the seed voxels in turn until the kth unmarked seed voxel, where k = 1, 2, ….

[0124] Step 3.2 Cluster the building seed voxels under spatial connectivity constraints, and use the clustering result as the seed voxels of each building monomer;

[0125] Step 3.2.1 Traverse all the seed voxels spatially connected to the kth seed voxel according to the spatial connectivity characteristics by using the depth-first strategy and mark them as L l , where l is the index of the mark label, l = 1, 2, ….

[0126] In this embodiment, for the voxel value vk = the k-th voxel of u. Traverse all the seed voxels that are 3D-connected to the k-th seed voxel according to the spatial connectivity characteristics by using the depth-first strategy. As Figure 6 shown, call the program LABEL(k, u, L l ) to label the seed voxels that are 3D-connected to the k-th seed voxel with the label L l .

[0127] Step 3.2.2 Continue to scan the unlabeled seed voxels in the intensity voxel dataset until all the seed voxels are labeled, obtaining several 3D-connected regions composed of seed voxels. Take the seed voxels labeled with the same label as the seed voxels of the same building monomer.

[0128] In this embodiment, the goal of the algorithm is to merge the spatially connected seed voxels into a 3D-connected region. Assume that the 3D seed voxel array V contains a total of l 3D-connected regions. The task of the algorithm is to assign l labels to each seed voxel in V so that the seed voxels belonging to the same 3D-connected region have the same label, while the seed voxels belonging to different 3D-connected regions have different labels. In the specific implementation process, based on the idea that voxels with spatial connectivity belong to the same category, use the 3D connected component labeling algorithm to divide and label V into l 3D-connected regions. The clustering result is as Figure 7 shown. It can be seen that most of the edge seed voxels belonging to the same building are clustered into the same 3D-connected region.

[0129] Step 3.3 Statistically analyze the spectral characteristics of the clustered seed voxels and use them as the spectral consistency measure adapted to each building monomer;

[0130] In this embodiment, statistically analyze the reflection intensity information of all the seed voxels in the set of clustered seed voxels. Based on the intensity values of each voxel set, form the spectral consistency measure of each building monomer, providing an intensity constraint condition for the subsequent search for building roof voxels relying on the seed voxels.

[0131] Step 4 Building roof extraction and optimization, the specific process is as Figure 8 shown;

[0132] Step 4.1 Extract the 3D-connected region of the building roof under the adaptive local spatial-spectral consistency constraint;

[0133] Step 4.1.1 Initialize an empty stack and put the set of seed voxels labeled with the label L l that have not been put into the stack for traversal into the stack;

[0134] Step 4.1.2 Pop the top element from the stack, and search for voxels in the non - negative voxels of the intensity voxel dataset that are spatially connected and have the same local intensity as it using the depth - first traversal strategy, and mark them with the same label L as this top - element of the stack l and push them onto the stack;

[0135] Step 4.1.3 If the stack is empty, terminate the program; otherwise, go to Step 4.1.2 until all seed voxels are put into the stack for traversal.

[0136] The above - mentioned voxels with the same local intensity are two voxels whose voxel value difference is less than the intensity difference threshold.

[0137] In this embodiment, the specific algorithm process is as follows: Assume that the previous clustering result of seed voxels is n seed - voxel clusters marked as L l (where l is the index of the marked label, l = 1,..., n, and n represents the number of clustering results), and use V l to represent the set of seed voxels after clustering. Assume that the intensity mean of the seed voxels in each cluster is μ and the standard deviation is σ. First, take all the seed voxels V l1 in the first label L1, and scan these seed voxels in turn. Among all the non - negative voxels of the intensity voxel model, starting from the voxels in the corresponding neighborhood (such as 26 - neighborhood or other neighborhood scales) of the seed voxels in the first label L1, as Figure 9 shown, call the program LABEL2(V l , μ, σ, L l ) to search for non - negative voxels that are spatially connected and have the same local intensity as them according to the depth - first traversal strategy and mark them as L1. After completion, continue to scan the seed voxels in other labels until all the seed voxels and all the non - negative voxels that are spatially connected and meet the local intensity consistency are marked. Finally, n 3D connected regions can be obtained again.

[0138] In this embodiment, the goal of this algorithm is to adaptively cluster the voxels belonging to the same building into the same connected region in a complex urban scene with the characteristics of "same object, different spectra; same spectra, different objects". Among them, "adaptive" means that based on the intensity values of the seed voxels in each seed - voxel set that has completed clustering statistics, establish a local spectral consistency measure for different building monomers, and mark the non - negative voxels that meet the spatial connectivity and local spectral consistency as the same connected region. Thus, the result that the voxels belonging to the same monomer building are clustered into the same 3D connected region can be obtained, and the clustering result is as Figure 10 shown. The above - mentioned local spectral consistency measure is the intensity difference threshold T g .

[0139] In this embodiment, different clustering results will be obtained and the subsequent building extraction results will be affected by applying different spatial neighborhood scales and different intensity difference thresholds during the above marking process. The optimal spatial neighborhood scale will be determined in the experiment. The optimal intensity difference threshold is determined by the following scheme: Based on the intensity values of the seed voxel blocks after all clustering is completed, the mean μ and standard deviation σ of the seed intensities of each block are obtained. Through experiments, 3σ is used as the optimal intensity difference threshold T g . The reason for choosing the multiplier 3 is that there are also slight deviations in intensity within the same building area, and using 3σ can extract buildings more accurately.

[0140] Step 4.2 Optimization of the building roof extraction result based on building area, local normal vector consistency, and density characteristics;

[0141] In this embodiment, the characteristics of the building roof are as follows: It has a certain area; the roof is generally flat, and the local normal vectors tend to be consistent; there is a density difference in the spatial distribution from other objects (such as vegetation, etc.). According to the above area, local normal vector consistency, and density characteristics, a scheme for searching for the 3D connected region corresponding to the building roof from n 3D connected regions is established:

[0142] Step 4.2.1 Eliminate the 3D connected regions that are not building roofs based on the area characteristics.

[0143] In this embodiment, let the minimum building area in the original airborne LiDAR point cloud dataset P be A min , and let the maximum building area in the original airborne LiDAR point cloud dataset P be A max .

[0144] The building area in the original airborne LiDAR point cloud dataset P refers to the projected area of the building on the horizontal plane.

[0145] For any 3D connected region, if its horizontal projected area is greater than or equal to A min and less than or equal to A max , then this 3D connected region is determined to be a building roof and is retained; otherwise, this 3D connected region is eliminated.

[0146] In this embodiment, A min and A max are constants, which are defined by the user according to the given situation of the original airborne LiDAR point cloud data.

[0147] Step 4.2.2 Eliminate the 3D connected regions that are not building roofs based on local normal vector consistency.

[0148] For any 3D connected region, all the voxels in it have normal vector values calculated according to the principal component analysis algorithm (PCA). Since the building roof is more planar compared to non-building objects such as trees, the angles between its normal vectors are smaller and tend to be consistent. Based on the principle of local normal vector consistency, let the threshold of the normal vector angle be T θ , if the angle between the normal vector of a single voxel in this region and the normal vectors of its neighboring voxels is within the threshold, then this voxel is retained; if not, this voxel is deleted. Finally, if the number of voxels with consistent normal vector angles remaining in this 3D connected region accounts for more than a certain proportion of the total number of voxels in this 3D connected region, then this region is determined to be a building roof and is retained; otherwise, this 3D connected region is excluded.

[0149] In this embodiment, the threshold T of the normal vector angle θ and the ratio are constants, and their values can be determined according to the degree of local normal vector consistency of each 3D connected region.

[0150] Step 4.2.3 Exclude 3D connected regions of non-building roofs based on density characteristics.

[0151] For any 3D connected region, if its density is greater than the given density threshold T d , then this 3D connected region is determined to be a building roof and is retained; otherwise, this 3D connected region is excluded;

[0152] In this embodiment, the density threshold T d is a constant, and its value can be determined according to the density distribution of each 3D connected region.

[0153] Step 5 Extract the building facades by combining the edge space constraints of the building roofs and the local intensity consistency. The specific process is as Figure 11 shown.

[0154] In this embodiment, the characteristics of the building facades are: perpendicular to the building roof contour, located within a certain range around the building roof contour, and their intensities conform to the local intensity consistency. According to the above characteristics, determine the scheme for searching for the 3D connected regions corresponding to the building facades:

[0155] Step 5.1 Project each detected 3D connected region of the building roof onto the XY plane, and then detect the edge contour voxels of each building roof;

[0156] Step 5.2 On the horizontal plane, with the edge contour voxels of the single building roof as the center, set the edge space constraint conditions of the building roof with a width of one voxel on both the inner and outer sides;

[0157] Step 5.3 For any non - negative voxel under the spatial constraint at the edge of the building roof, if the difference between its reflection intensity value and the mean intensity of the building edge voxels is within plus or minus two standard deviations, then this voxel is determined to be an exterior wall voxel of this building.

[0158] In this embodiment, the reference data uses the building standard data provided by the ISPRS Working Group III / 4 (experimental data accurately classified into building point sets and non - building point sets) to quantitatively evaluate the calculation accuracy of the method of the present invention.

[0159] The result of the building extraction method proposed by the present invention is represented in the form of voxels, while the buildings in the reference data are expressed in discrete airborne LiDAR point cloud data. To compare with the reference data to evaluate the accuracy of the method proposed by the present invention, first, count the number of original airborne LiDAR point cloud data contained in the building voxels extracted by this method, and then compare with the reference data. Furthermore, use type - I error (the proportion of building laser points misclassified as non - building laser points), type - II error (the proportion of non - building laser points misclassified as building laser points), total error (the proportion of misclassified building laser points), correct rate (the proportion of correctly extracted building laser points in the total number of building laser points in the extraction result), completeness (the proportion of correctly extracted building laser points in the total number of building laser points in the standard data), quality, and Kappa coefficient to quantitatively evaluate the effectiveness of the building extraction method proposed by the present invention.

[0160] Table 1 shows the accuracy indicators of the building extraction results corresponding to the building extraction operations under adaptive local spatial - spectral consistency based on the intensity voxel data of two experimental data when the spatial neighborhood scales are 6, 18, 26, 56, 80, and 124 in this embodiment. The data in this table aims to examine the influence of different neighborhood scales on the building extraction results and thereby determine the optimal neighborhood scale.

[0161] Table 1 Accuracy indicators of building extraction results with different neighborhood scales

[0162]

[0163] As can be seen from Table 1, the average Kappa coefficients obtained at neighborhood sizes of 6, 18, 26, 56, 80, and 124 are 38.48%, 76.99%, 82.87%, 94.50%, 90.38%, and 84.75% respectively. This shows that: (1) The 56-neighborhood corresponds to the largest Kappa coefficient. Therefore, from the perspective of the Kappa coefficient index, the 56-neighborhood is the optimal neighborhood scale. (2) The increase in neighborhood scale does not necessarily mean an improvement in extraction accuracy. The idea of this algorithm is that building information can be propagated based on the spatial connectivity and local intensity consistency defined in the voxel array. Taking the 6-neighborhood as an example, building information can only be propagated in 6 directions: up, down, left, right, front, and back. As a result, only the voxels on flat-roof buildings can be merged into a 3D connected region and correctly extracted, while the voxels on spire-roof buildings may be divided into multiple 3D connected regions and thus misjudged by subsequent characteristics such as area and density. This can explain why the highest type I error is generated using the 6-neighborhood (see the type I error corresponding to different neighborhood scales in the fourth column of Table 1). As the neighborhood scale increases, the number of propagation directions increases, and more voxels are classified as buildings. Using neighborhoods of 18, 26, 56, 80, and 124 produces better results than the 6-neighborhood. However, if the neighborhood scale is too large, some non-building voxels may be misjudged as buildings, resulting in an increase in type II error (see the type II error using the 80-neighborhood), which can explain why the accuracy of the 80-neighborhood is lower than that of the 56-neighborhood. (3) The average total errors of 6, 18, 26, 56, 80, and 124 are 24.43%, 10.08%, 7.44%, 2.49%, 4.47%, and 7.39% respectively. This shows that the 56-neighborhood corresponds to the smallest total error. Therefore, from the perspective of the total error index, the 56-neighborhood is also the optimal neighborhood scale.

[0164] The building extraction results obtained by applying the method of the present invention are as Figure 13 shown. Among them, the building extraction results are in the form of building voxels, and the size of a single voxel is 0.4m × 0.4m × 0.4m, as Figure 12 shown by the black cubes in

[0165] Table 2 is a quantitative evaluation of the building extraction accuracy at the 56-neighborhood scale of 2 experimental data with reference data as the standard in this embodiment.

[0166] Table 2 Accuracy of Building Extraction Results

[0167]

[0168] As can be seen from Table 2, the average completeness rate, average accuracy rate, and average quality score for building extraction are 95.89%, 97.09%, and 93.22%, respectively. This verifies the effectiveness of the method proposed in the present invention.

[0169] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, rather than to limit them; although the present invention has been described in detail with reference to the foregoing embodiments, those of ordinary skill in the art should understand that they can still modify the technical solutions described in the foregoing embodiments, or perform equivalent replacements for some or all of the technical features; and these modifications or replacements do not cause the essence of the corresponding technical solutions to deviate from the scope defined by the claims of the present invention.

Claims

1. An airborne LiDAR data building extraction method based on adaptive local spatial-spectral consistency, characterized in that It includes the following steps: Step 1: Read the original airborne LiDAR point cloud data to form an original airborne LiDAR point cloud data set; Step 2: Regularize the original airborne LiDAR point cloud data set into an intensity voxel data set; Step 3: Search for seed voxels of building monomers and establish their spectral consistency measures; Step 3.1: Search for building seed voxels in the non-negative intensity voxel data based on the elevation jump and edge straight-line characteristics of buildings; Step 3.2: Cluster the building seed voxels under spatial connectivity constraints, and use the clustering results as the seed voxels of each building monomer; Step 3.3: Statistically analyze the spectral characteristics of the clustered seed voxels and use them as the spectral consistency measures adapted to each building monomer; Step 4: Extract and optimize the building roofs; Step 4.1: Extract the 3D connected regions of building roofs under adaptive local spatial-spectral consistency constraints; Step 4.2: Optimize the extraction results of the 3D connected regions of building roofs based on building area, local normal vector consistency, and density characteristics; Step 5: Extract the building facades by combining the spatial constraints of the building roof edges and local intensity consistency; Step 5.1 Project the optimized 3D connected regions of each building roof onto a plane, and then detect the edge contour voxels of each building roof; Step 5.2: On the horizontal plane, set the spatial constraint conditions for the building roof edges with the edge contour voxels of the single-building roof as the center and one voxel width on both the inner and outer sides; Step 5.3: For any non-negative voxel under the spatial constraints of the building roof edges, if the difference between its reflection intensity value and the mean intensity of the building edge voxels is within the range of plus or minus two standard deviations, then determine this voxel as a building facade voxel.

2. An airborne LiDAR data building extraction method based on adaptive local spatial-spectral consistency according to claim 1, characterized in that The specific steps of Step 2 include the following steps: Step 2.1: Remove abnormal data from the original airborne LiDAR point cloud data set to obtain a data set with abnormal data removed; Step 2.2: Regularize the data set with abnormal data removed into an intensity voxel data set.

3. An airborne LiDAR data building extraction method based on adaptive local spatial-spectral consistency according to claim 2, characterized in that The specific steps of Step 2.1 include the following steps: Step 2.1.1 Render and display the original airborne LiDAR point cloud dataset according to elevation, visually judge the true elevation range of the ground objects from the side view, and set the highest elevation threshold , the lowest elevation threshold ; Step 2.1.2 For each laser point in the original airborne LiDAR point cloud dataset, if its elevation value is higher than the maximum elevation threshold or lower than the minimum elevation threshold , then this laser point is abnormal data and is excluded. Otherwise, this laser point is retained, and finally, a dataset with elevation anomalies excluded is obtained; Step 2.1.3 For the dataset with elevation anomalies removed obtained in the above steps, count the frequencies of the intensity values of each laser point therein, and visually display the statistical results in the form of a histogram, and visually determine the intensity threshold corresponding to the true ground object and denote it as ; Step 2.1.4 Determine the laser points with intensity values higher than in the dataset after removing the elevation anomaly as intensity anomaly data and remove them, finally obtaining the dataset after removing anomalies.

4. An airborne LiDAR data building extraction method based on adaptive local spatial-spectral consistency according to claim 2, characterized in that The specific steps of Step 2.2 include the following steps: Step 2.2.1: Represent the three-dimensional space range with the axis-parallel bounding box of the data set with abnormal data removed; Step 2.2.2 Determine the resolution of the voxel in the x, y, and z directions according to the average point spacing of the laser points in the abnormal data set removed ; The resolution of the voxel in the x, y, and z directions The calculation formula is as follows: ； Among them, is the two-dimensional point set obtained by projecting the abnormal data set onto a plane, is the point set 's convex hull, is the convex hull 's area; Step 2.2.3 According to the voxel resolution Divide the axially parallel bounding box to obtain a 3D voxel grid, and each 3D voxel grid cell is a voxel; The voxel resolution-based The axially parallel bounding box is divided into a 3D voxel grid and represented by a 3D voxel array as follows: ； Where V is the set of voxels in the 3D voxel array, is the voxel index; is the number of voxels; It is voxel value of a voxel; It is The coordinates of the voxel in the voxel array, For the row number, For column number, is the layer number; The number of voxels in the direction is The number of voxels in the direction is The number of voxels in the direction; , , The calculation formula is as follows: ； Among them, is the ceiling operator, , , , , is the index for removing data from the abnormal dataset, is the number of data removed from the abnormal dataset. When removing the th data from the abnormal dataset, its coordinates are ; Voxel number The calculation formula is as follows: Step 2.2.4: Map each laser point in the data set with abnormal data removed into a 3D voxel grid, and then assign values to each voxel according to the intensity attributes of the laser points contained in the 3D voxel grid to obtain an intensity voxel data set; The specific process of assigning values to each voxel according to the intensity attributes of the laser points contained in the 3D voxel grid is: Assign the voxel containing the laser point with the average value of the laser point intensity, assign the voxel without the laser point with -1, and further discretize the non-negative voxel values to , to obtain the voxel values of each.

5. An airborne LiDAR data building extraction method under adaptive local spatial-spectral consistency according to claim 1, characterized in that The specific steps of Step 3.2 include the following steps: Step 3.2.1 Traverse all the seed voxels that are spatially connected to the th seed voxel according to the spatial connectivity characteristic by using the depth-first strategy, and mark them as , where is the index of the marking label, ; Step 3.2.2: Continue to scan the unmarked seed voxels in the intensity voxel data set until all seed voxels are marked, obtaining several 3D connected regions composed of seed voxels, and using the seed voxels marked with the same label as the seed voxels of the same building monomer.

6. An airborne LiDAR data building extraction method under adaptive local spatial-spectral consistency according to claim 1, characterized in that The specific steps of Step 4.1 include the following steps: Step 4.1.1 Initialize an empty stack and put the set of seed voxels labeled with the label that have not been put into the stack for traversal into the stack; Step 4.1.2 Pop the top element from the stack, search for voxels that are spatially connected and have the same local intensity in the non-negative voxel data of the intensity voxel dataset using the depth-first traversal strategy, and label them with the same label as this top element of the stack and push them into the stack; The voxels with local intensity consistency are two voxels with a difference in voxel values less than the intensity difference threshold; Step 4.1.3 If the stack is empty, terminate the program; otherwise, go to Step 4.1.2 until all seed voxels are put into the stack for traversal, thereby obtaining several 3D connected regions of building roofs.

7. An airborne LiDAR data building extraction method under adaptive local spatial-spectral consistency according to claim 1, characterized in that The specific steps of Step 4.2 are as follows: Step 4.2.1 Eliminate 3D connected regions of non-building roofs based on area characteristics; The rejection conditions are as follows: the minimum building area in the original airborne LiDAR point cloud dataset is , and the maximum building area in the original airborne LiDAR point cloud dataset is . For any 3D connected region, if its horizontal projected area is greater than or equal to and less than or equal to , then this 3D connected region is determined to be a building roof and is retained; otherwise, this 3D connected region is rejected. The building area in the original airborne LiDAR point cloud dataset refers to the projected area of the building on the horizontal plane; Step 4.2.2 Eliminate 3D connected regions of non-building roofs based on local normal vector consistency; For any 3D connected region where all voxels have normal vector values calculated according to the principal component analysis algorithm PCA, based on the principle of consistent local normal vectors, the normal vector angle threshold is , if the normal vector angle between a single voxel and its neighboring voxels within this region is within the threshold, then this voxel is retained; if it is not within the threshold, then this voxel is deleted. Finally, if the number of voxels with consistent normal vector angles remaining in this 3D connected region exceeds a set proportion of the total number of voxels in this 3D connected region, then this region is determined to be a building roof and is retained; otherwise, this 3D connected region is excluded; Step 4.2.3 Eliminate 3D connected regions of non-building roofs based on density characteristics; Set the density threshold to , for any 3D connected region, if its density is greater than the given density threshold , then the 3D connected region is determined to be a building roof and retained, otherwise the 3D connected region is removed.

Images