Multi-feature graph cut building extraction method for large-scale airborne LiDAR data

Abstract

The invention relates to a multi-feature graph cut building extraction method for large-scale airborne LiDAR data, and belongs to the technical field of laser radar point cloud. The method comprises the following specific steps: 1) preprocessing airborne LiDAR point cloud data, and dividing the point cloud data into ground points and non-ground points; 2) calculating point features and surface features of the point cloud data, and normalizing the features; 3) constructing a graph for the point cloud data by using the thought of a graph cut algorithm, representing the graph by using an energy function, and finding the minimum cut of the graph by using a maximum flow / minimum cut algorithm to obtain an initial building extraction result; and 4) performing post-processing on an initial building extraction result to realize stable extraction of the building. The method can effectively solve the problems that when large-scale airborne LiDAR point cloud data is oriented, high-precision building extraction is difficult to achieve through a single feature, and many features are complex in calculation and high in calculation cost.

Description

Technical Field

[0001] The present invention belongs to the technical field of laser radar point cloud, and in particular relates to a multi-feature graph cutting building extraction method for large-scale airborne LiDAR data. Background Art

[0002] Buildings are one of the important features in cities, and the extraction of buildings is a key step in many applications, such as urban planning, disaster assessment, cadastral management, smart cities, and many other aspects. At present, airborne LiDAR technology has the advantages of high efficiency, high measurement accuracy, and little interference from the external environment, and has become one of the conventional means of urban data acquisition.

[0003] Building extraction algorithms are mainly divided into deep learning methods and traditional methods. The deep learning method uses a neural network model to learn features from LiDAR data through a large amount of training, and uses labeled data for verification to obtain a building extraction network. At present, this type of method has achieved good results on small data sets, but this type of method is very dependent on a large amount of data to train the model, and this process is very time-consuming. At the same time, due to the complexity of the urban environment, model training requires a large number of labeled samples. Therefore, this type of method is difficult to apply on a large range of data sets.

[0004] The traditional method is to extract buildings by mining the features of point cloud data itself and fusing the features. This method does not require data annotation and works well on large-scale data sets. However, this method still has two difficulties and challenges: (1) When facing large-scale data, many features are complex to calculate and require high computational costs. (2) Using only a single feature often does not achieve the desired effect. How to fuse multiple effective features to achieve higher-precision pre-earthquake building extraction.

[0005] In view of the above problems, the present invention discloses a multi-feature graph cutting building extraction method for large-scale airborne LiDAR data. Summary of the invention

[0006] The technical problem to be solved by the present invention is that it is difficult to extract buildings with high precision using a single feature when facing large-scale airborne LiDAR point cloud data, and the calculation of many features is complex and costly.

[0007] The advantages of the present invention are that the data items in the energy function are reconstructed by fusing point features and surface features, the smoothing items in the energy function are constructed by using height values ​​while considering the contextual information of the point neighborhood, and the maximum flow / minimum cut algorithm is used for robust extraction of buildings. In order to achieve the above technical objectives, the present invention adopts the following technical solutions.

[0008] The technical solution of the present invention is: a multi-feature graph cutting building extraction method for large-scale airborne LiDAR data, and the specific steps are as follows: Step 1: Preprocessing Firstly, statistical filtering is used to remove noise points in the airborne LiDAR point cloud, then the point cloud is resampled using raster resampling, and finally the point cloud data is divided into ground points and non-ground points using a ground filtering algorithm.

[0009] :Feature calculation and normalization First, the point features of the point cloud data are calculated, an initial tensor value is assigned to each point, and then the influence of the neighboring points within a certain radius of each point on the point is summarized to obtain the tensor value of the point after summary, and the summarized tensor value is decomposed to obtain the eigenvalue and eigenvector, and the point feature of the point is calculated using the eigenvalue. Then the surface features of the point cloud data are calculated. The present invention uses the homogeneity features of the gray-level co-occurrence matrix as the surface features of each point. Finally, considering that the value ranges of point features and surface features are significantly different, they are normalized using a logical function.

[0010] : Building extraction Based on the idea of ​​graph cut algorithm, a graph is constructed for point cloud data and represented by energy function. The data items of energy function are reconstructed by point features and surface features, and the height value of each point is used to construct the smoothing term in the energy function. The maximum flow / minimum cut algorithm is used to find the minimum cut of the graph to obtain the initial building extraction result.

[0011] : Post-processing Connectivity analysis was performed on the initial building extraction results to remove some noise data caused by vegetation, and morphological filters were used to fill holes and optimize the building extraction results.

[0012] The beneficial effects of the present invention are as follows: the data items of the energy function are reconstructed using point features and surface features, and then the height value of each point is used to construct the smoothing term in the energy function. The maximum flow / minimum cut algorithm can be used to obtain building data at one time, and only a small number of parameters need to be set to realize the extraction of buildings from large-scale airborne LiDAR point cloud data. BRIEF DESCRIPTION OF THE DRAWINGS

[0013] Figure 1 is a flow chart of the present invention; Figure 2 It is a data display diagram of the present invention; wherein, Figure 2 (a) is the location of the city of Wechingen, Figure 2 (b) is an image of the city of Weichingen; Figure 2 (c) Airborne LiDAR data of the city of Weichingen; Figure 3 The following are LiDAR point cloud data displays for three typical regions of the present invention; Figure 3 (a) is the elevation coloring display of area 1. Figure 3 (b) is the elevation coloring display of area 2. Figure 3 (c) Elevation coloring display of area 3; Figure 4 These are three typical regional result diagrams after preprocessing of the present invention; wherein, Figure 4 (a) Coloring the non-ground point elevations in region 1. Figure 4 (b) Coloring the non-ground point elevations in region 2, Figure 4 (c) Coloring the elevation of non-ground points in region 3; Figure 5 The three typical regional point feature coloring diagrams of the present invention are shown in FIG. 1 , wherein: Figure 5 (a) PIV coloring for region 1, Figure 5 (b) PIV coloring for region 2, Figure 5 (c) PIV coloring for region 3; Figure 6 The three typical regional surface feature coloring diagrams of the present invention are shown in FIG. 1 , wherein: Figure 6 (a) Coloring the GLCM homogeneity features for region 1, Figure 6 (b) Coloring the GLCM homogeneity features for region 2, Figure 6 (c) Coloring the GLCM homogeneity features of region 3; Figure 7 These are the three typical regional building extraction result maps of the present invention; wherein, Figure 7 (a) is the building extraction result of area 1. Figure 7 (b) is the building extraction result of area 2. Figure 7 (c) is the building extraction result of area 3; Figure 8 The point feature colored map of the entire city of Weixingen of the present invention; Fig. 9 This is a building extraction result map of the entire city of Weixingen according to the present invention. DETAILED DESCRIPTION

[0014] Combination Figure 1 The specific steps of the multi-feature graph cutting building extraction method for large-scale airborne LiDAR data of the present invention are described. Figure 2 The data used in the present invention are described below. The present invention uses the commonly used Vaihingen urban LiDAR dataset provided by the International Society for Photogrammetry and Remote Sensing for experiments. The dataset contains image data and airborne LiDAR data. Figure 2(b) Description: The experimental data is located in the main urban area, which contains more than 1,800 buildings of various types and sizes. The vegetation coverage in different building areas is inconsistent. In addition, there are special features such as bridges. It is a typical urban area, which can greatly verify the robustness of the algorithm. Figure 2 (c) The area of ​​the entire LiDAR point cloud data is about 2.0 km 2 The point cloud data contains 5 flight strips, with an average flight strip overlap of 30% and an average point density of 4-7 points / m 2 The data has an accuracy of about 0.1 m in both horizontal and vertical directions. Figure 3 Note: This example first conducts building extraction experiments on the three most commonly used typical areas, and finally extracts buildings from the entire area.

[0015] : Preprocessing Combination Figure 4 Description: Preprocess the airborne LiDAR point cloud data to divide the point cloud data into ground points and non-ground points. Specifically, first use statistical filtering to remove noise points in the airborne LiDAR point cloud, then use raster resampling to resample the point cloud, and finally use the ground filtering algorithm to classify the point cloud data into ground points. Specifically, mark the ground point for each point category and record the point cloud data as P ={ x , y , z , l},in x , y , z , are the spatial coordinate values ​​of the point cloud, l Indicates the category of the point. If the category of the point is a ground point, then l The value of is 2, otherwise l The value of is 0.

[0016] :Feature calculation and normalization For the preprocessed point cloud data, calculate its point features and surface features, and normalize the features; the specific steps are as follows: (1) Calculation point features Assume that the total number of point cloud data is n , assign a second-order tensor value T to each point in the point cloud data, denoted as P i ={ x , y , z , T}( i =1,2,…, n )Tensor T It is defined by the following formula: (1) Among them is Eigenvalues, Considering that the initial features of each point are not obvious, the eigenvalues ​​of the initial tensor value of each point are set to 1 The eigenvectors are , and record the initial tensor value as T ball , called the spherical tensor.

[0017] Take any point in the point cloud p i ( i =1,2,…, n ) as an example, the radius of the point is calculated as r In the vicinity m Neighborhood points q j ( j =1,2,…, m ) on its tensor value, and obtain the tensor value of the point , which is specifically expressed by the following formula: (2) in, Represents the calculation of neighborhood points p j Point p i The tensor value after the influence of Respectively indicate points p i and neighboring points q j The initial tensor value of T ball , is the attenuation function, which is defined as follows: (3) in d is the Euclidean distance between two points, is a scaling factor that controls the attenuation strength, usually 1.2 times the radius. w is the weight factor, expressed as , r is the radius of the search.

[0018] For tensor values Decompose to get eigenvalues ​​and eigenvectors, and use eigenvalues ​​to calculate the point p i The point eigenvalue of is denoted as PIV and is expressed by the following formula: (4) According to the above steps, by traversing each point in the point cloud P, the point feature value of each point can be calculated, recorded as P i ={ x , y , z , PIV , l} ( i =1,2,…, n ).

[0019] Combination Figure 5 Note: After calculating the point features, each point is displayed by color mapping.

[0020] (2) Calculating surface features The number of pixels in the image is determined based on the range of the x and y values ​​of the point cloud data and the set step size. The elevation value of the point within each pixel is mapped to the grayscale value of the image, where the grayscale value of the pixel at the ground point is assigned to 0, and the point cloud data is projected into a two-dimensional height image.

[0021] Then, the gray level co-occurrence matrix of the two-dimensional height image at four directions of 0°, 45°, 90°, and 135° is calculated respectively, and the homogeneity feature is calculated on this basis, which is recorded as Hom , expressed using the following formula: (5) in m and n is the gray value of the height image, G is the maximum gray value, is the joint conditional probability density, d is the step value, is the direction of generating the gray-level co-occurrence matrix.

[0022] The average values ​​of the four directions are taken as the homogeneity feature value of the two-dimensional height image. Since the grid sampling is performed in advance in Step 1, there is only one point in each pixel. The two-dimensional height image after calculating the homogeneity feature is inversely transformed, and each point in the point cloud can obtain its homogeneity feature value, that is, its surface feature, recorded as P i ={ x , y , z , PIV , Hom , l} ( i =1,2,…, n ).

[0023] Combination Figure 6 Note: After calculating the surface features, each point is displayed by color mapping.

[0024] (3) Feature normalization Considering that the value ranges of point features and surface features are significantly different, the features are normalized using a logical function, which is expressed by the following formula: in x 0 is the feature threshold, k Controls the steepness of the logistic function curve.

[0025] The normalized point features and surface features are recorded as f PIV and f Hom , point cloud data is represented as P={ x , y , z , f PIV , f Hom , l}.

[0026] : Building extraction The idea of ​​graph cut algorithm is used to construct a graph for point cloud data. Specifically, the graph has background points and foreground points. These two nodes are connected to each point in the point cloud data. In addition, each point is also connected to its adjacent points. The graph is represented by energy function, which is expressed as , defined using the following formula: (6) Among them, the first is a data item, the second item is the smoothing term. is the weight of the smoothing term.

[0027] Data Item Indicate point p i Assigning Tags l The penalty is constructed using point features and surface features, which represents the connection between each point of the point cloud data in the figure and the background point and foreground point, and is defined using the following formula: (7) in w 1 and w 2 Represents the weights of point features and surface features, satisfying w 1 + w 2 =1.

[0028] Smooth term Indicate pointp i Its neighboring points q j The penalty for discontinuity between is constructed using the height value of each point, representing the connection between adjacent points in the graph, and is defined using the following formula: (8) in Indicates the height of the point, It affects the steepness of the curve, and its value is related to the height of the building in the LiDAR data.

[0029] The maximum flow / minimum cut algorithm is used to segment the graph to achieve building extraction; specifically, the maximum flow / minimum cut algorithm is used to find the minimum cut of the graph, and a building label is assigned to each point in the point cloud data, thereby obtaining the initial building extraction result.

[0030] Post-processing Connectivity analysis is performed on the initial building extraction results to remove some noise data caused by vegetation, and morphological filters are used to fill holes and optimize the building extraction results. Figure 7 Illustrate the building extraction results for three typical areas.

[0031] This embodiment reconstructs the data items in the energy function by fusing point features and surface features, and uses the height value to construct the smoothing item in the energy function while considering the context information of the point neighborhood. The maximum flow / minimum cut algorithm is used to segment the graph at one time to achieve robust extraction of buildings. Figure 8 This shows the point feature calculation results for the entire area. Fig. 9 The experimental results show that the proposed method can robustly extract buildings from large-scale airborne LiDAR point cloud data.

[0032] The above embodiments are only used to illustrate the present invention, but not to limit the present invention. A person skilled in the art may make various changes and modifications without departing from the spirit and scope of the present invention. Therefore, all equivalent technical solutions also belong to the scope of the present invention, and the patent protection scope of the present invention shall be defined by the claims.

[0033] The technical contents not described in detail in the present invention are all known technologies.

Claims

1. A multi-feature graph cutting building extraction method for large-scale airborne LiDAR data, characterized by: Step 1: Preprocess the airborne LiDAR point cloud data and divide the point cloud data into ground points and non-ground points; Step 2: Calculate the point features and surface features of the point cloud data and normalize the features; Step 3: Use the idea of ​​graph cut algorithm to construct a graph for point cloud data and represent it using energy function. Use maximum flow / minimum cut algorithm to find the minimum cut of the graph and obtain the initial building extraction result. Step 4: Post-process the initial building extraction results to achieve robust extraction of buildings.

2. The multi-feature graph cut building extraction method for large-scale airborne LiDAR data according to claim 1, characterized in that: In the above Step 3, the idea of ​​graph cut algorithm is used to construct a graph for point cloud data. The graph has background points and foreground points. These two nodes are connected to each point in the point cloud data. In addition, each point is also connected to its adjacent points. The graph is represented by energy function, which is recorded as , defined using the following formula: The first item is the data item , indicating a point p i Assigning Tags l The penalty is constructed using point features and surface features, which represents the connection between each point of the point cloud data in the figure and the background point and foreground point, and is defined using the following formula: in w 1 and w 2 Represents the weights of point features and surface features, satisfying w 1 + w 2 =1; In the second item is the weight of the smoothing term, is a smoothing term, indicating that the point p i Its neighboring points q j The penalty for discontinuity between is constructed using the height value of each point, representing the connection between adjacent points in the graph, and is defined using the following formula: in Indicates the height of the point, It affects the steepness of the curve, and its value is related to the height of the building in the LiDAR data.

Images