A method for enhancing the integrity of LiDAR point cloud buildings based on an intensity voxel model

Abstract

This invention provides a method for enhancing the integrity of buildings in LiDAR point clouds based on an intensity voxel model, belonging to the field of remote sensing data processing technology. This method can fill small, localized voids on building surfaces caused by factors such as occlusion and uneven point cloud density. The algorithm first regularizes the LiDAR point cloud data into an intensity voxel model; then, it successively detects the building roof and facade using 3D connected component construction theory and intensity analysis theory under roof spatial constraints; finally, for the defects of small, localized voids caused by target occlusion and uneven point cloud density, it fills the voids using the three-dimensional morphological void-filling theory of rotatable planar structural elements. This method provides an effective scheme for filling small, localized voids on the surface of LiDAR point clouds, contributing to the development of LiDAR point cloud data processing and applications based on intensity voxel model theory.

Description

Technical Field

[0001] This invention belongs to the field of remote sensing data processing technology, specifically relating to a method for enhancing the integrity of LiDAR point cloud buildings based on an intensity voxel model. Background Technology

[0002] Building 3D reconstruction technology has wide applications in digital cities, real estate, urban planning and construction, and is one of the current research hotspots. Airborne and vehicle-mounted LiDAR (Light Detection and Ranging) are two mainstream methods for acquiring the roofs and facades of buildings in large-scale urban scenes. By comprehensively utilizing airborne and vehicle-mounted LiDAR point clouds, relatively complete basic data of buildings can be provided for 3D reconstruction. However, due to factors such as occlusion, complete transmission or absorption of ground objects causing no reflection echo, there are often unavoidable data deficiencies such as local data holes on the surface of building LiDAR point clouds, which are difficult to avoid even with multiple scans. Existing LiDAR point cloud enhancement methods are mainly divided into: (1) Enhancement based on the self-information of LiDAR point clouds. Starting from the point cloud itself, the deficiencies in the data and the relationship of the point cloud around the deficiencies are often used to fill in the data holes. That is, the neighboring point information, topological structure information, etc. are mainly used to fill in the data holes through data synthesis algorithms such as interpolation and surface fitting. However, this type of method is difficult to guarantee the objectivity of filling in some complex areas without repeating texture features. (2) Enhancement is achieved based on the complementary information of LiDAR point cloud and image. By using photogrammetry to collect entity image data around the hole location in the point cloud, the image data is solved into original image point cloud data using the Structure From Motion (SFM) algorithm and Patch-based Multi-view Stereo (PMVS) based on the patch model. After the registration of the two types of point cloud data, the deletion of the original image point cloud data in the overlapping area and the selection of the optimal density of the image point cloud data, the fusion between the two types of data is completed, thereby completing the repair. However, this type of method requires a lot of auxiliary work and has low work efficiency. (3) Enhancement based on complementary information from airborne and vehicle-mounted LiDAR point clouds: Since neither airborne nor vehicle-mounted LiDAR systems can acquire complete building point clouds independently, and the data collected by both systems are point clouds, airborne and vehicle-mounted LiDAR data have better compatibility compared to the fusion of aerial imagery and spectral information data sources. Based on this, point clouds of corresponding parts of similar buildings can be used to fill local data gaps. However, the existence of similar features is a prerequisite for this filling process; if these features are not present, filling cannot be achieved. Therefore, this invention proposes a method for enhancing the integrity of LiDAR point clouds based on an intensity voxel model. Summary of the Invention

[0003] The technical problem to be solved by the present invention is to provide a method for enhancing the integrity of LiDAR point cloud buildings based on the strength voxel model, which addresses the shortcomings of the prior art.

[0004] To solve the above-mentioned technical problems, the technical solution adopted by the present invention is as follows:

[0005] A method for enhancing the integrity of LiDAR point cloud buildings based on an intensity voxel model includes the following steps:

[0006] Step 1: Read the raw LiDAR point cloud data to form the raw LiDAR point cloud dataset;

[0007] Step 2: Regularize the original LiDAR point cloud dataset into an intensity voxel model;

[0008] Step 2.1: Remove elevation anomalies and intensity anomalies from the original LiDAR point cloud data to obtain the anomaly-removed dataset;

[0009] Step 2.1.1: Count the frequency of the elevation values ​​of each laser point in the original LiDAR point cloud data, and visualize the statistical results in the form of a histogram;

[0010] Step 2.1.2: Determine the highest elevation threshold T corresponding to the actual terrain. hi and minimum elevation threshold T li ;

[0011] Step 2.1.3: For each laser point in the original LiDAR point cloud data, if its elevation value is higher than the highest elevation threshold T... hi Or below the minimum elevation threshold T li If the laser point is positive, it is considered abnormal data and is removed; otherwise, the laser point is retained, and the final dataset with removed elevation anomalies is obtained.

[0012] Step 2.1.4: Count the frequency of intensity values ​​of each laser point in the dataset that removes elevation anomalies, and visualize the statistical results in the form of a histogram;

[0013] Step 2.1.5: Based on the frequency histogram results of the intensity values, visually determine the highest intensity threshold T corresponding to the actual terrain and features. he and minimum intensity threshold T le ;

[0014] Step 2.1.6: For each laser point in the dataset that has been removed due to elevation anomalies, if its intensity value is higher than the highest intensity threshold T... he or below the minimum intensity threshold T leIf the laser point is positive, it is considered an intensity anomaly and is removed; otherwise, the laser point is retained. This process yields a dataset of removed elevation and intensity anomalies.

[0015] Step 2.2: Regularize the outlier dataset into an intensity voxel model;

[0016] Step 2.2.1: Represent the 3D spatial extent using the axial bounding boxes (AABB type bounding boxes are selected for the axial parallel bounding boxes, where AABB refers to the smallest hexahedron with all sides parallel to the coordinate axes and containing all point clouds) of the removed outlier dataset;

[0017] Step 2.2.2: Determine the resolution (Δx, Δy, Δz) of the voxel in the x, y, and z directions, i.e., the voxel size, based on the average point spacing of the laser points in the dataset after removing outliers;

[0018] Step 2.2.3: Divide the axial bounding box according to the voxel resolution (△x, △y, △z) to obtain a 3D voxel grid, which is represented by a 3D voxel array (r rows, c columns and l layers). Each 3D voxel grid unit is a voxel.

[0019] Step 2.2.4: Map each laser point in the removed outlier dataset to a 3D voxel grid, and then assign values ​​to each voxel based on the intensity attributes of the laser points contained in the 3D voxel grid to obtain the intensity voxel model.

[0020] The specific process of assigning values ​​to each voxel based on the intensity properties of the laser points contained in the 3D voxel grid is as follows:

[0021] Voxels containing laser points are assigned the average intensity of the laser points, and voxels without laser points are assigned the value 0. The voxel values ​​are further discretized to {0, ..., 255} to obtain the voxel values.

[0022] Step 3: Based on the characteristics of the building's roof and facade, detect the 3D connected regions formed by the building's roof and facade in sequence to complete the LiDAR point cloud building detection and individualization based on the intensity voxel model;

[0023] Step 3.1 Detect building roofs based on the theory of 3D connected components. The specific steps are as follows:

[0024] Step 3.1.1: Based on the elevation jump characteristics of building edge points, search for building edge voxels from the strength voxel model as seed voxels G. s , where s = 1, 2, ...;

[0025] Step 3.1.2: Label the 3D connected regions of the seed voxels. For any unlabeled seed voxel G... s A depth-first search strategy is used to search for the seed voxel G in the intensity voxel model.s 3D connectivity and intensity difference less than intensity difference threshold T g All unlabeled voxels, and labeled as L k Where k is the index of the label, k = 1, 2, ... until all unlabeled seed voxels' 3D connected regions are labeled, resulting in several 3D connected regions. The specific steps are as follows:

[0026] 3.1.2.1: Initialize an empty stack, and put G... s Store it in the stack;

[0027] 3.1.2.2: Pop the top element of the stack and obtain elements that are 3D connected to the top element and whose intensity difference is less than the intensity difference threshold T. g All unlabeled voxels are labeled as building roof voxels and stored in the stack;

[0028] 3.1.2.3: Determine if the stack is empty. If it is, all building roof voxels in the intensity voxel model are marked. Otherwise, return to step 3.1.2.2.

[0029] Step 3.1.3: Based on area and strength characteristics, optimize the voxel set of the building roof to complete the building roof inspection. The specific steps are as follows:

[0030] 3.1.3.1: Eliminate 3D connected regions that are not building roofs based on area characteristics;

[0031] 3.1.3.2: Eliminate 3D connected regions of non-building roofs based on strength characteristics.

[0032] Step 3.2: Complete the individualization of building rooftops by clustering based on topological analysis of the connected regions of the rooftop space. The specific steps are as follows:

[0033] Step 3.2.1: Describe the topological relationships of the connected regions of each rooftop based on the topology matrix;

[0034] The topology matrix is ​​a matrix that represents the inclusion, adjacency or other topological relationships between rooftop connected regions. It can be represented by a two-dimensional matrix, where the rows and columns represent the extracted rooftop connected regions.

[0035] Step 3.2.2: Cluster the connected regions of each rooftop with topological adjacency and containment relationships into individual buildings. The specific process is as follows:

[0036] Let S and O correspond to two connected roof regions, with O having a larger horizontal projected area. Let the voxel sets within S and O be denoted as s and o, respectively. If the two are determined to be in an inclusive relationship and belong to the same building, then the roof of the building is represented by O. Furthermore, if the distance d(S, O) between S and O is less than a certain threshold (the threshold can be set according to the minimum distance between buildings specified in the standard), then the two are determined to be adjacent and belong to the same building, and the roof of the building is represented by S and O.

[0037] Step 3.3: Based on buffer analysis, detect the 3D connected regions of the building facade to complete the building facade detection;

[0038] Step 3.3.1: Detect the outline of the roof of each building;

[0039] Step 3.3.2: On the horizontal plane, establish a buffer zone with the roof outline of any building as the center and the width of one voxel on both the inner and outer sides;

[0040] Step 3.3.3: For non-zero voxels located within the buffer zone with an elevation direction density (i.e., the number of voxels with the same planar coordinates) greater than a certain density threshold, if their reflection intensity value is within ±2 standard deviations of the mean reflection intensity value of voxels within the buffer zone, then the voxel is classified as a facade voxel; otherwise, it is classified as a non-facade voxel. The density threshold is 0.22.

[0041] Step 4: Morphological local fine void filling based on rotatable planar structural elements;

[0042] Step 4.1: Centering each building vole on the c-axis and using the c-axis as the normal vector v i Construct a target template F with a range of 9*9;

[0043] Step 4.2: Rotate the target template F in the vertical and horizontal directions based on Rodrigues' rotation formula.

[0044] v rot =cosθv i +(1-cosθ)(v i ·k)+sinθ×v i

[0045] Where k is the rotation axis, θ is the rotation angle of the target template, and v rot The normal vector of the new target template is obtained by rotating the target template;

[0046] Step 4.3: Based on neighborhood search, find the set of neighborhood voxels I within a certain spatial neighborhood centered on the current voxel, and determine the positional relationship between all neighborhood voxels and the target template F in each direction. If F = 0, then determine that the voxel is located on the target template.

[0047] Step 4.4: Apply a 3*3 dimensional structural element (SE = [1,1,1; 1,1,1; 1,1,1]) to each target template to perform a shape closure operation, thereby filling small local voids in the building model.

[0048] The beneficial effects of adopting the above technical solution are as follows: The method proposed in this invention is based on voxel structure design, which makes good use of the implicit neighborhood relationship between voxels in 3D voxel data; the designed rotatable structural element can adapt to the detection and filling of voids on the surface of buildings in all directions, and the final result helps the development of LiDAR point cloud data processing and application based on voxel modeling theory. Attached Figure Description

[0049] Figure 1 This is a flowchart of a LiDAR point cloud building integrity enhancement method based on an intensity voxel model, as described in a specific embodiment of the present invention.

[0050] Figure 2 This refers to the original airborne and vehicle-mounted LiDAR point cloud data in a specific embodiment of the present invention;

[0051] Among them, (a) is airborne LiDAR point cloud data, and (b) is vehicle-mounted LiDAR point cloud data;

[0052] Figure 3 This is a flowchart illustrating the process of regularizing raw airborne and vehicle-mounted LiDAR point cloud data into an intensity voxel model in a specific embodiment of the present invention.

[0053] Figure 4 This is a side view of the intensity voxel model obtained by regularizing LiDAR point cloud data in a specific embodiment of the present invention.

[0054] Among them, (a) is the intensity voxel model corresponding to airborne LiDAR, and (b) is the intensity voxel model corresponding to vehicle-mounted LiDAR;

[0055] Figure 5 This is a flowchart illustrating the building inspection process in a specific embodiment of the present invention;

[0056] Figure 6 A top view of the buffer zone setup in a specific embodiment of the present invention;

[0057] Among them, dark gray voxels represent the extracted building roof outline, and light gray voxels represent the inner and outer buffer zones;

[0058] Figure 7 This is a specific embodiment of the present invention to complete the LiDAR point cloud building detection results based on the intensity voxel model;

[0059] Among them, (a) is the building detection result of airborne LiDAR, and (b) is the building detection result of vehicle-mounted LiDAR;

[0060] Figure 8 This is a schematic diagram of the rotating planar structural elements designed in a specific embodiment of the present invention;

[0061] Figure 9 This is the result of filling small morphological cavities in a specific embodiment of the present invention;

[0062] Among them, (a) is the filling result of airborne LiDAR, and (b) is the filling result of vehicle-mounted LiDAR. Detailed Implementation

[0063] The specific embodiments of the present invention will be described in further detail below with reference to the accompanying drawings and examples. The following examples are for illustrative purposes only and are not intended to limit the scope of the invention.

[0064] In this embodiment, the method was implemented using MATLAB 9.0.0 platform on an Intel(R) Core(TM) i5-6400 CPU@2.70GHz 2.71GHz, 4GB memory, and Windows 10 Ultimate 64-bit system. The effectiveness of the method was further verified by evaluating its accuracy.

[0065] A method for enhancing the integrity of LiDAR point cloud buildings based on an intensity voxel model, such as... Figure 1 As shown, it includes the following steps:

[0066] Step 1: Read the raw LiDAR point cloud data to form the raw LiDAR point cloud dataset.

[0067] In this embodiment, airborne and vehicle-mounted LiDAR point clouds from different scenarios are used to verify the effectiveness and feasibility of the proposed algorithm. Building roof data is acquired using an airborne laser scanning system, with an average point cloud density of 8.7 points / m². 2 There are a total of 87,542 data points. This scene includes isolated low-rise residential areas, some of which are surrounded by vegetation, such as... Figure 2 As shown in (a), the building facade data was acquired using a vehicle-mounted laser scanning system, with an average point cloud density of 1651 points / m². 2 There are a total of 1,066,651 data points. This scene includes building facades, street trees, and vehicles, among other features. Figure 2 As shown in (b).

[0068] In this embodiment, the original LiDAR point cloud dataset P = {p i (x i ,y i ,zi ), i = 1, ..., n}, where i is the index of the original LiDAR point cloud data, n is the number of original LiDAR point cloud data, p i This is the i-th original LiDAR point cloud data, with coordinates (x, y). i ,y i ,z i ).

[0069] Step 2: Regularize the raw LiDAR point cloud data into an intensity voxel model. The specific process is as follows: Figure 3 As shown;

[0070] Step 2.1: Remove elevation and intensity anomalies from the original LiDAR point cloud data to obtain the anomaly-removed dataset;

[0071] Step 2.1.1: Count the frequency of the elevation values ​​of each laser point in the original LiDAR point cloud data, and visualize the statistical results in the form of a histogram;

[0072] Step 2.1.2: Determine the highest elevation threshold T corresponding to the actual terrain. hi and minimum elevation threshold T li ;

[0073] Step 2.1.3: For each laser point in the original LiDAR point cloud data, if its elevation value is higher than the highest elevation threshold T... hi Or below the minimum elevation threshold T li If the laser point is positive, it is considered abnormal data and is removed; otherwise, the laser point is retained, and the final dataset with removed elevation anomalies is obtained.

[0074] Step 2.1.4: Count the frequency of intensity values ​​of each laser point in the dataset that removes elevation anomalies, and visualize the statistical results in the form of a histogram;

[0075] Step 2.1.5: Determine the highest intensity threshold T corresponding to the actual terrain and features. he and minimum intensity threshold T le ;

[0076] Step 2.1.6: For each laser point in the dataset that has been removed due to elevation anomalies, if its intensity value is higher than the highest intensity threshold T... he or below the minimum intensity threshold T le If the laser point is positive, it is considered an intensity anomaly and is removed; otherwise, the laser point is retained. This process yields a dataset of removed elevation and intensity anomalies.

[0077] In this embodiment, the removal of abnormal datasets is denoted as Q = {q} i’ (x i’ ,y i’ ,zi’ ), i'=1,…,t}, where i' is the index of the data in the dataset to be removed from the abnormal dataset, t is the number of data in the dataset to be removed from the abnormal dataset, and q i’ It removes the i'th data point from the abnormal dataset, whose coordinates are (x, y, y). i’ ,y i’ ,z i’ ).

[0078] In this embodiment, the highest elevation threshold T hi and minimum elevation threshold T li This is a constant, and its value needs to be determined based on the spatial distribution of the original LiDAR point cloud data. The highest elevation threshold T for airborne LiDAR point cloud data... hi Take 285m and the lowest elevation threshold T li Take 260m; the highest elevation threshold T for vehicle-mounted LiDAR point cloud data. hi Take 270m and the lowest elevation threshold T li Take 260m.

[0079] Step 2.2: Regularize the outlier dataset into an intensity voxel model;

[0080] Step 2.2.1: Represent the 3D spatial extent using the axial bounding boxes of the excluded outlier dataset;

[0081] In this embodiment, the spatial extent of the outlier dataset Q can be determined by its axially aligned bounding box (AABB). AABB = {x, y, z | x min ≤x≤x max ,y min ≤y≤y max ,z min ≤z≤z max}, where (x min ,y min ,z min ) and (x max ,y max ,z max () represents the minimum and maximum values ​​of the x, y, and z coordinates in the dataset after removing outliers, respectively. 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},z max =max{zi’ ,i'=1,…,t},z min =min{z i’ ,i'=1,…,t}.

[0082] Step 2.2.2: Determine the resolution (Δx, Δy, Δz) of the voxel in the x, y, and z directions, i.e., the voxel size, based on the average point spacing of the laser points in the outlier dataset Q;

[0083] In this embodiment, the formulas for calculating the resolutions Δx, Δy, and Δz of the voxels in the x, y, and z directions are shown in equation (1):

[0084]

[0085] Among them, S xy ={(x i’ ,y i’ Let {i'=1,…,t} be the two-dimensional point set obtained by projecting the outlier dataset Q onto the XOY plane, and let C(S) be the point set. xy ) is the point set S xy A(C(S) convex shell xy )) is a convex shell C(S) xy The area of ​​).

[0086] Step 2.2.3: Divide the axial bounding box according to the voxel resolution (△x, △y, △z) to obtain a 3D voxel grid, which is represented by a 3D voxel array (r rows, c columns and l layers). Each 3D voxel grid unit is a voxel.

[0087] In this embodiment, the axial bounding box can be divided into a 3D voxel grid based on the voxel resolution (△x, △y, △z), which is represented by a 3D voxel array. Let G be the set of voxels in the 3D voxel array, as shown in equation (2):

[0088] G={g j (r j ,c j ,l j ),j=1,…,m} (2)

[0089] Where j is the voxel index, m is the number of voxels, g is the voxel value, and (r,c,l) are the coordinates (row, column, and layer number) of the voxel in the 3D 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. R, C, and L are given by formula (3):

[0090]

[0091] in, The rounding up 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},z max =max{z i’ ,i'=1,…,t},z min =min{z i’ ,i'=1,…,t}.

[0092] Therefore, the number of voxels m can be derived as shown in equation (4):

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

[0094] Step 2.2.4: Map each laser point in the outlier dataset Q to 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 the intensity voxel model.

[0095] In this embodiment, each laser point in the abnormal dataset Q is mapped to a 3D voxel grid, and then each voxel is assigned a value based on the intensity attribute of the laser points contained in the 3D voxel grid. Specifically, voxels containing laser points are assigned the average intensity of the laser points, and voxels without laser points are assigned 0, as shown in equation (5).

[0096]

[0097] in, This is the floor function operator. The voxels are further discretized to {0,…,255} to obtain individual voxel values. This yields the intensity voxel model, completing the regularization of the outlier removal dataset.

[0098] In this embodiment, the airborne and vehicle-mounted LiDAR point clouds are divided into 3D arrays of 0.375m×0.375m×0.25m and 0.125m×0.125m×0.125m, respectively. The two 3D arrays contain a total of 64,783 and 85,489 target voxels, respectively. Figure 4 As shown.

[0099] Step 3: Based on the characteristics of the building's roof and facade, sequentially detect the 3D connected regions formed by the building's roof and facade to complete the LiDAR point cloud building detection based on the intensity voxel model. The specific process is as follows: Figure 5 As shown;

[0100] Step 3.1 Detect building roofs based on the theory of constructing 3D connected regions. The specific steps are as follows:

[0101] Step 3.1.1: Based on the elevation jump characteristics of building edge points, search for building edge voxels from the strength voxel model as seed voxels G. s , where s = 1, 2, ...;

[0102] In this embodiment, for any non-valued voxel in the intensity voxel model, if the elevation difference between the valued voxel and the valued voxels in its 8-level neighborhood is greater than the elevation difference threshold T, e (Constant, 2m), then the voxel is determined to be the building edge voxel, and it is used as the seed voxel G. s .

[0103] Step 3.1.2: Label the 3D connected regions of the seed voxels. For any unlabeled seed voxel G... s A depth-first search strategy is used to search for the seed voxel G in the intensity voxel model. s 3D connectivity and intensity difference less than intensity difference threshold T g All unlabeled voxels, and labeled as L k Where k is the index of the label, k = 1, 2, ... until all unlabeled seed voxels' 3D connected regions are labeled, resulting in several 3D connected regions. The specific steps are as follows:

[0104] 3.1.2.1: Initialize an empty stack, and put G... s Store it in the stack;

[0105] 3.1.2.2: Pop the top element of the stack and obtain elements that are 3D connected to the top element and whose intensity difference is less than the intensity difference threshold T. g All unlabeled voxels are labeled as building roof voxels and stored in the stack;

[0106] 3.1.2.3: Determine if the stack is empty. If it is, all building roof voxels in the intensity voxel model are marked. Otherwise, return to step 3.1.2.2.

[0107] Step 3.1.3: Based on area and strength characteristics, optimize the voxel set of the building roof to complete the building roof inspection. The specific steps are as follows:

[0108] In this embodiment, the characteristics of the building roof are: it has a certain area and its spatial distribution differs in intensity from that of other targets (such as vegetation). Based on the aforementioned area and intensity characteristics of the building roof, the 3D connected region of the building roof obtained above is optimized, that is, individual 3D connected regions that are not part of the building roof are removed.

[0109] 3.1.3.1: Eliminate 3D connected regions that are not building roofs based on area characteristics;

[0110] In this embodiment, let the minimum building area in the original LiDAR point cloud dataset P be A. min Let the area of ​​the largest building in the original LiDAR point cloud dataset P be A. max .

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

[0112] 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 If the condition is met, the 3D connected region is determined to be the roof of a building and is retained; otherwise, the 3D connected region is discarded.

[0113] In this embodiment, A min and A max This is a constant, defined by the user based on the given raw LiDAR point cloud data.

[0114] 3.1.3.2: Eliminate 3D connected regions of non-building roofs based on strength characteristics.

[0115] In this embodiment, for any 3D connected region, if its intensity is less than a given intensity threshold T... s If the 3D connected region is identified as the roof of a building, it will be retained; otherwise, the 3D connected region will be removed.

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

[0117] Step 3.2: Complete the individualization of building rooftops by clustering based on topological analysis of the connected regions of the rooftop space. The specific steps are as follows:

[0118] Step 3.2.1: Describe the topological relationships of the connected regions of each rooftop based on the topology matrix;

[0119] The topology matrix is ​​a matrix that represents the inclusion, adjacency or other topological relationships between rooftop connected regions. It can be represented by a two-dimensional matrix, where the rows and columns represent the extracted rooftop connected regions.

[0120] Step 3.2.2: Cluster the connected regions of each rooftop with topological adjacency and containment relationships into individual buildings. The specific process is as follows:

[0121] Let S and O correspond to two connected roof regions, with O having a larger horizontal projected area. Let the voxel sets within S and O be denoted as s and o, respectively. If the two are determined to be in an inclusive relationship and belong to the same building, then the roof of the building is represented by O. Furthermore, if the distance d(S, O) between S and O is less than a certain threshold (the threshold can be set according to the minimum distance between buildings specified in the standard), then the two are determined to be adjacent and belong to the same building, and the roof of the building is represented by S and O.

[0122] Step 3.3: Based on buffer analysis, detect the 3D connected regions of the building facade to complete the building facade detection;

[0123] Step 3.3.1: Detect the outline of the roof of each building;

[0124] Step 3.3.2: On the horizontal plane, establish a buffer zone with the roof outline of any building as the center and a width of one voxel on both the inner and outer sides, as shown below. Figure 6 As shown;

[0125] Step 3.3.3: The consistency criterion for the building's reflectance intensity is determined by the range of ±2 standard deviations of the mean reflectance intensity of voxels within the buffer zone whose elevation density (i.e., the number of voxels with the same planar coordinates) is greater than a certain density threshold. If the voxel value of a voxel within the buffer zone falls within this range, the voxel is classified as a facade voxel; otherwise, it is classified as a non-facade voxel. The density threshold is 0.22.

[0126] In this embodiment, after building detection processing based on the intensity voxel model, 9199 and 46451 building voxels were detected from the airborne and vehicle-mounted LiDAR point clouds, respectively. The building detection results show that all buildings in the scene data were successfully detected. The results are as follows... Figure 7 As shown, where, Figure 7 (a) shows the results of airborne LiDAR building detection. Figure 7 (b) shows the results of building detection using a vehicle-mounted LiDAR system.

[0127] Step 4: Morphological local fine void filling based on rotatable planar structural elements;

[0128] Step 4.1: Centering each building vole on the c-axis and using the c-axis as the normal vector v i Construct a target template F with a range of 9*9;

[0129] Step 4.2: Based on Rodriguez's rotation formula, rotate the target template F in the vertical and horizontal directions, v rot =cosθv i +(1-cosθ)(v i ·k)+sinθ×v i

[0130] Where k is the rotation axis, θ is the rotation angle of the target template, and v rot This is the normal vector of the new target template obtained after rotating the original target template;

[0131] In this embodiment, to achieve rotation of the target template along the vertical and horizontal directions, the rotation axis k is taken as the unit vector of the vertical center axis and the unit vector of the horizontal center axis of the target template, respectively, where the rotation angle θ is 30°, thus achieving the construction of the target template in 11 directions, as follows. Figure 8 As shown;

[0132] Step 4.3: Based on neighborhood search, find the set of neighborhood voxels I within a certain spatial neighborhood centered on the current voxel, and determine the positional relationship between all neighborhood voxels and the target template F in each direction. If F = 0, then determine that the voxel is located on the target template.

[0133] Step 4.4: Apply a 3*3 dimensional structural element (SE = [1,1,1; 1,1,1; 1,1,1]) to each target template to perform a shape closure operation, thereby filling small local voids in the building model.

[0134] In this embodiment, after morphological local small void filling processing based on rotatable planar structural elements, 4838 and 765 voxels were filled in the airborne and vehicle-mounted LiDAR voxel models, respectively. Figure 9 As shown, where, Figure 9 (a) shows the filling results from the airborne LiDAR. Figure 9 (b) shows the filling results from the vehicle-mounted LiDAR.

[0135] In this embodiment, to quantitatively evaluate the point cloud filling results, corresponding experiments were designed to verify the effectiveness of the algorithm. For filling morphologically based local small holes, the filling rate E was calculated. p (%) (Total number of voxels filled C) T Z, representing the total number of voids in the original voxel model T The accuracy of the algorithm is evaluated by the proportion of the algorithm's fill-in points.

[0136] Table 1 shows the accuracy index of the cavity filling results in this embodiment, which is based on voxel data from airborne and vehicle-mounted LiDAR experiments and performs morphological-based filling of local small cavities.

[0137] Table 1. Accuracy evaluation of morphological-based filling of local fine cavities.

[0138]

[0139] As shown in Table 1, the filling rates of airborne and vehicle-mounted LiDAR point clouds using the morphology-based local small hole filling method are 95.49% and 94.21%, respectively. This verifies the effectiveness of the method proposed in this invention.

[0140] This invention provides a method for enhancing the integrity of LiDAR point cloud structures based on a strength voxel model. Based on 3D connected component construction theory, and utilizing a search and labeling method based on voxel spatial neighborhood relationships and geometric features, it effectively leverages the implicit neighborhood relationships between voxels in 3D voxel data. The designed rotatable planar structural element can adaptively detect and fill voids in all directions. This method contributes to the development of LiDAR point cloud data processing and applications based on voxel theory. The filling rate of this method for small, localized voids can reach over 90%.

[0141] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, and not to limit them; although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some or all of the technical features therein; and these modifications or substitutions 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. A method for enhancing the integrity of LiDAR point cloud buildings based on an intensity voxel model, characterized in that, Includes the following steps: Step 1: Read the raw LiDAR point cloud data to form the raw LiDAR point cloud dataset; Step 2: Regularize the original LiDAR point cloud dataset into an intensity voxel model; Step 3: Based on the characteristics of the building's roof and facade, detect the 3D connected regions formed by the building's roof and facade in sequence to complete the LiDAR point cloud building detection and individualization based on the intensity voxel model; Step 3 specifically includes the following steps: Step 3.1: Detect building roofs based on the theory of constructing three-dimensional connected regions; Step 3.1 specifically includes the following steps: Step 3.1.1: Based on the elevation jump characteristics of building edge points, search for building edge voxels from the strength voxel model as seed voxels G. s , where s = 1, 2, ...; Step 3.1.2: Label the 3D connected regions of the seed voxels. For any unlabeled seed voxel G... s A depth-first search strategy is used to search for the seed voxel G in the intensity voxel model. s 3D connectivity and intensity difference less than intensity difference threshold T g All unlabeled voxels, and labeled as L k , where k is the index of the label, k = 1, 2, ... until all unlabeled seed voxels' 3D connected regions are labeled, resulting in several 3D connected regions; Step 3.1.3: Based on area and strength characteristics, optimize the voxel set of the building roof to complete the building roof inspection; Step 3.2: Complete the individualization of building rooftops based on topological clustering of rooftop spatial connectivity regions; Step 3.2 specifically includes the following steps: Step 3.2.1: Describe the topological relationships of the connected regions of each rooftop based on the topology matrix; The topology matrix is ​​a matrix representing the inclusion, adjacency, or other topological relationships between rooftop connected regions. It is shown as a two-dimensional matrix, where the rows and columns represent the extracted rooftop connected regions. Step 3.2.2: Cluster the connected roof regions with topological adjacency and containment relationships into individual buildings; Step 3.3: Based on buffer analysis, detect the 3D connected regions of the building facade to complete the building facade detection; Step 3.3 specifically includes the following steps: Step 3.3.1: Detect the outline of the roof of each building; Step 3.3.2: On the horizontal plane, establish a buffer zone with the roof outline of any building as the center and the width of one voxel on both the inner and outer sides; Step 3.3.3: For non-zero voxels with an elevation direction density greater than a certain density threshold located within the buffer zone, if their reflection intensity value is within ±2 standard deviations of the mean reflection intensity value of voxels within the buffer zone, then the voxel is judged as a facade voxel; otherwise, it is a non-facade voxel. The density threshold is 0.22, and the elevation direction density within the buffer zone is the number of voxels with the same planar coordinates. Step 4: Morphological local fine void filling based on rotatable planar structural elements; Step 4 specifically includes the following steps: Step 4.1: Centering each building vole on the c-axis and using the c-axis as the normal vector v i Construct a target template F with a range of 9*9; Step 4.2: Rotate the target template F in the vertical and horizontal directions based on Rodrigues' rotation formula. ； Where k is the axis of rotation. v is the rotation angle of the target template. rot The normal vector of the new target template is obtained by rotating the target template; Step 4.3: Based on neighborhood search, find the set of neighborhood voxels I within a certain spatial neighborhood centered on the current voxel, and determine the positional relationship between all neighborhood voxels and the target template F in each direction. If F=0, then determine that the voxel is located on the target template. Step 4.4: Apply a 3*3 shaped structural element SE=[1,1,1;1,1,1;1,1,1] to each target template to perform a shape closure operation, thereby filling small local voids in the building model.

2. The method for enhancing the integrity of LiDAR point cloud buildings based on an intensity voxel model according to claim 1, characterized in that, Step 2 specifically includes the following steps: Step 2.1: Remove elevation anomalies and intensity anomalies from the original LiDAR point cloud data to obtain the anomaly-removed dataset; Step 2.2: Regularize the outlier dataset into an intensity voxel model.

3. The method for enhancing the integrity of LiDAR point cloud buildings based on an intensity voxel model according to claim 2, characterized in that, Step 2.1 specifically includes the following steps: Step 2.1.1: Count the frequency of the elevation values ​​of each laser point in the original LiDAR point cloud data, and visualize the statistical results in the form of a histogram; Step 2.1.2: Determine the highest elevation threshold T corresponding to the actual terrain. hi and minimum elevation threshold T li ; Step 2.1.3: For each laser point in the original LiDAR point cloud data, if its elevation value is higher than the highest elevation threshold T... hi Or below the minimum elevation threshold T li If the laser point is positive, it is considered abnormal data and is removed; otherwise, the laser point is retained, and the final dataset with removed elevation anomalies is obtained. Step 2.1.4: Count the frequency of intensity values ​​of each laser point in the dataset that removes elevation anomalies, and visualize the statistical results in the form of a histogram; Step 2.1.5: Based on the frequency histogram results of the intensity values, visually determine the highest intensity threshold T corresponding to the actual terrain and features. he and minimum intensity threshold T le ; Step 2.1.6: For each laser point in the dataset that has been removed due to elevation anomalies, if its intensity value is higher than the highest intensity threshold T... he or below the minimum intensity threshold T le If the laser point is positive, it is considered an intensity anomaly and should be removed; otherwise, the laser point should be retained. This process yields a dataset of removed elevation and intensity anomalies.

4. The method for enhancing the integrity of LiDAR point cloud buildings based on an intensity voxel model according to claim 2, characterized in that, Step 2.2 specifically includes the following steps: Step 2.2.1: Represent the 3D spatial extent using the axial bounding boxes of the excluded outlier dataset; Step 2.2.2: Determine the resolution Δx, Δy, and Δz of the voxel in the x, y, and z directions, i.e., the voxel size, based on the average point spacing of the laser points in the dataset after removing outliers; Step 2.2.3: Divide the axial bounding box according to the voxel resolution Δx, Δy, Δz to obtain a 3D voxel grid. Each 3D voxel grid unit in r rows, c columns and l layers is represented by a 3D voxel array; Step 2.2.4: Map each laser point in the removed outlier dataset to a 3D voxel grid, and then assign values ​​to each voxel based on the intensity attributes of the laser points contained in the 3D voxel grid to obtain the intensity voxel model.

5. The method for enhancing the integrity of LiDAR point cloud buildings based on an intensity voxel model according to claim 1, characterized in that, Step 3.1.2 specifically includes the following steps: 3.1.2.1: Initialize an empty stack, and put G... s Store it in the stack; 3.1.2.2: Pop the top element of the stack and obtain elements that are 3D connected to the top element and whose intensity difference is less than the intensity difference threshold T. g All unlabeled voxels are labeled as building roof voxels and stored in the stack; 3.1.2.3: Determine if the stack is empty. If it is, all building roof voxels in the intensity voxel model are marked. Otherwise, return to step 3.1.2.2.

Images