An Optimization Method for UAV-borne LiDAR Point Cloud Filtering

Abstract

An optimization method for filtering LiDAR point clouds carried by an unmanned aerial vehicle. The method includes establishing an index of the point cloud using a KD tree; calculating and determining the maximum search radius for identifying isolated points r ; corresponding to the neighborhood threshold k 1 and the maximum search radius r for resampling; corresponding to the neighborhood threshold k 2 and the maximum search radius r for resampling. By establishing an index of the point cloud using a KD tree, the identification and removal of isolated points are achieved based on the search for neighboring points; the local tangent plane fitting at the position of each sampling point is realized using the moving least squares method, and the resampling of the point cloud is realized based on the method of projective transformation, thereby realizing the filtering and optimization processing of the original sampling data; greatly weakening the influence of sampling errors on the point cloud data, improving the consistency between the obtained data and the collected point cloud data, enhancing the universality of the airborne LiDAR point cloud data, and broadening the application field of lightweight and economically applicable unmanned aerial vehicle borne LiDAR devices.

Description

Technical Field

[0001] The present invention relates to an optimization method for filtering LiDAR point clouds carried by unmanned aerial vehicles, and belongs to the technical field of remote sensing image processing. Background Art

[0002] With the continuous progress of sensor technology, LiDAR (Light Detection And Ranging) provides a more convenient means for people to collect ground object and geomorphic feature points. By measuring the time difference between the emitted laser and the reflected laser, the distance between the sensor and the geographical entity can be calculated. Through the effective integration with GPS (Global Positioning System) and IMU (Inertial Measurement Unit), LiDAR can directly obtain a set of three-dimensional data points that can characterize the geographical features of the earth's surface.

[0003] According to the different sensor platforms, LiDAR can be classified into airborne LiDAR (also known as ALS, Airborne Laser Scanning), vehicle-mounted LiDAR, ground LiDAR, and handheld LiDAR. Among them, airborne LiDAR has attracted attention due to its high data acquisition efficiency, wide application fields, and high operation safety. Most airborne LiDARs use manned aircraft as the sensor platform. Due to the high cost of data acquisition, the popularity of airborne LiDAR in practical applications is not high. In recent years, thanks to the progress of sensor platform technology, airborne LiDAR has been developing more and more towards the direction of lightweight and economical applicability. However, it inevitably reduces the data acquisition accuracy, limits the application fields of unmanned aerial vehicle (UAV)-borne LiDAR point clouds, and also poses higher requirements for subsequent point cloud data processing. Summary of the Invention

[0004] The object of the present invention is to provide an optimization method for filtering UAV-borne LiDAR point clouds, which can realize the filtering and optimization processing of the original sampling data, improve the consistency between the obtained data and the collected point cloud data, and broaden the application scope of lightweight and economical UAV-borne LiDAR devices.

[0005] To achieve the above object, the present invention provides an optimization method for filtering UAV-borne LiDAR point clouds, including the following steps:

[0006] S1. Establish an index of the point cloud using a KD tree;

[0007] S2. Calculate and determine the maximum search radius for identifying isolated points r ;

[0008] S3. Resample corresponding to the neighborhood threshold k 1 and the maximum search radius r for resampling;

[0009] S4. Resample corresponding to the neighborhood threshold k 2 and the maximum search radius r for resampling.

[0010] Furthermore, the specific process of S1 is as follows:

[0011] S1.1. Build a KD tree:

[0012] S1.1-1. Divide the points along the first dimension of the sampling point coordinates x direction. Select the point with the middle coordinate in the first dimension of all sampling points as the splitting point, and divide the space where the sampling points are located into two subspaces, the left and the right; for all sampling points in the left space, their first dimension coordinates are less than the median point, and for all sampling points in the right space, their coordinates are greater than or equal to the median point;

[0013] S1.1-2. For the sampling points in the left and right subspaces, select the second dimension of the sampling points y direction for further point division. For the left and right spaces, respectively select the point with the middle coordinate in the second dimension as the splitting point, and divide the left and right spaces into two subspaces again; for all sampling points in the left subspace of the left space, their second dimension coordinates are less than the median point, and for all sampling points in the right subspace, their coordinates are greater than or equal to the median point; for all sampling points in the left subspace of the right space, their second dimension coordinates are less than the median point, and for all sampling points in the right subspace, their coordinates are greater than or equal to the median point;

[0014] S1.1-3. For the four subspaces obtained by the division in S1.1-2, select the third dimension of the sampling points z direction and perform further point division according to the method of S1.1-1 or S1.1-2, and further divide the four subspaces into eight subspaces;

[0015] S1.1-4. Repeat S1.1-1 to S1.1-3 until all sampling points are added to the KD tree, and the KD tree construction process ends;

[0016] S1.2. k - Neighborhood search: Given the airborne LiDAR point cloud set , select any sampling point S in the set , and implement the search for the neighborhood points of the sampling point in two steps: Determine the search path and find the neighboring points, specifically:

[0017] S1.2-1. Determine the search path:

[0018] S1.2-1-1. Construct a stack for storing the nodes on the search path P S , and set it to an empty stack;

[0019] S1.2-1-2. Take the root node of the KD tree as the starting point, push onto the stack, and compare the coordinates of and in the first dimension. If , take the root node of the left subtree as the new starting point and search downward. Otherwise, if , search downward along the root node of the right subtree as the new starting point;

[0020] S1.2-1-3. Check if is empty. If it is, the search process ends. Otherwise, push the new root node onto the stack, and compare the coordinates of and in the second dimension. If , take the root node of the left subtree as the new starting point and search downward. Otherwise, if , search downward along the root node of the right subtree as the new starting point;

[0021] S1.2-1-4. Check if is empty. If it is, the search process ends. Otherwise, push the new root node onto the stack, and compare the coordinates of and in the third dimension. If , take the root node of the left subtree as the new starting point and search downward. Otherwise, if , search downward along the root node of the right subtree as the new starting point;

[0022] S1.2-1-5. Repeat S1.2-1-2 to S1.2-1-4 until the search path encounters a leaf node of the KD tree;

[0023] S1.2-2. Search for neighboring points: Assume the current sampling point is . After determining the search path for neighboring points, set up a neighboring point list vertList={} and a distance list distList={} to store the neighboring points and their corresponding distances respectively. Then, based on the backtracking method, determine the k neighboring points of the current sampling point. The steps are as follows:

[0024] S1.2-2-1. Check the stack Whether it is an empty stack. If so, the adjacent point search process ends; otherwise, go to S1.2-2-2;

[0025] S1.2-2-2. Pop the top element from in turn according to the principle of last in first out, calculate the distance between and . If =0, consider and to be the same point, go back to S1.2-2-1; otherwise, go to S1.2-2-3;

[0026] S1.2-2-3. If the length of vertList is less than k at this time, directly add and to the end of the corresponding list, that is, add to the list vertList, add to the list distList; if the length of vertList is equal to k , judge compared with the maximum distance stored in the list distList. If , then delete from distList and delete the corresponding adjacent point from vertList, and insert into the corresponding position of the distList list in ascending order of adjacent point distance, and insert into the corresponding position of vertList, ensuring that the subscript of in vertList is the same as the subscript in distList;

[0027] S1.2-2-4. Judge whether is a leaf node. If so, go back to S1.2-2-1; if not, go to S1.2-2-5;

[0028] S1.2-2-5. With as the center and as the radius, draw a circle, and judge whether there is a coincidence between the circle and the hyperplane where is located. If there is a coincidence, then perform adjacent point search from the other half space of the hyperplane where is located. Take as the starting point according to the path search method, determine the search path according to the method described in S1.2-1, and add all the nodes on the path to the stack and return to S1.2-2-1.

[0029] Further, the specific process of S2 is as follows:

[0030] S2.1. Calculate the maximum search radius:

[0031] Based on the constructed KD tree, select all or a part of the airborne LiDAR point cloud set as the sample for determining the maximum search radius r , denoted as ; traverse S r each sampling point in , and find the corresponding k 2 nearest neighbors in ascending order of distance, record the distance between this nearest neighbor and the current sampling point as r i . After all sampling points in the sample S r are traversed, take the average value of the maximum search radii r i of all traversed sampling points as the maximum search radius of the entire point cloud r , and the calculation formula is as follows: ;

[0032] S2.2. Isolated point identification and removal:

[0033] Based on the constructed point cloud KD tree index, traverse each sampling point in the airborne LiDAR point cloud set , set the nearest neighbor set of the current sampling point to be empty, add the sampling point to the nearest neighbor set , and perform the following processing:

[0034] S2.2-1. Sequentially find the nearest neighbors of in ascending order of distance, and handle them in the following two cases:

[0035] S2.2-1-1. If the Euclidean distance between and is less than r and the number of points in the nearest neighbor set is less than k 1, add the nearest neighbor to the nearest neighbor set , and continue to traverse the next sampling point;

[0036] S2.2-1-2. If the Euclidean distance between is greater than r , stop searching;

[0037] S2.2-2. Count the number of points in the set of neighboring points . If the number is less than 3, that is, the number of sampling points in the set of neighboring points does not meet the requirements of plane fitting, then consider as an outlier, delete it from the original point cloud, and traverse the next sampling point.

[0038] Furthermore, the process of S3 is to traverse each sampling point in , take the first neighboring points from the set of neighboring points of the current sampling point , denoted as

[0039] S3.1. Calculate the centroid of all sampling points in the set of neighboring points : ;

[0040] S3.2. Construct the covariance matrix: ;

[0041] S3.3. Calculate the eigenvalues and eigenvectors of the covariance matrix , select the eigenvector corresponding to the smallest eigenvalue as the normal of the tangent plane, denoted as ; Thus, the equation of the fitted plane is determined as: ;

[0042] where ;

[0043] Expand the above formula to get: ;

[0044] where ;

[0045] S3.4. Project the current sampling point onto the fitted tangent plane, and assign the projected point coordinates to the current point as the new coordinates of the current point, specifically as follows:

[0046] S3.4-1. Calculate the directed distance from the point d to the fitted tangent plane: ;

[0047] S3.4-2, Calculation Point The projected point coordinates on the fitted tangent plane are as follows: .

[0048] Furthermore, the process of S4 is to expand the search range of the neighborhood and traverse each sampling point in again, and accordingly implement the secondary resampling process of the point cloud data, which is specifically as follows:

[0049] S4.1. Set the neighborhood set of the current sampling point to be empty;

[0050] S4.2. Search for the neighboring points of in ascending order of distance, and handle them in the following two cases:

[0051] S4.2-1. If the Euclidean distance between and is less than r and the number of points in the neighboring point set is less than k 2, add the neighboring point to the neighboring point set , and continue to traverse the next sampling point;

[0052] S4.2-2. If the Euclidean distance between and is greater than r , stop the search;

[0053] S4.3. Based on the methods described in S3.1~S3.3, use the current point and all points in the neighboring point set to fit the tangent plane;

[0054] S4.4. Based on the method described in S3.4, project the current point onto the fitted tangent plane, and assign the coordinates of the projected point to the current point as the new coordinates of the current point;

[0055] S4.5. After the traversal is completed, output the processed point cloud data to a file.

[0056] The present invention realizes the identification and elimination of isolated points by using a KD-tree to establish an index of the point cloud and based on the search for neighboring points; uses the moving least squares method to realize the fitting of the local tangent plane at the position of each sampling point, and realizes the resampling of the point cloud based on the method of projection transformation, thereby realizing the filtering and optimization processing of the original sampling data; greatly weakens the influence of sampling errors on the point cloud data, improves the consistency between the obtained data and the collected point cloud data, and then effectively improves the quality of the DEM constructed based on the airborne LiDAR point cloud, improves the universality of the airborne LiDAR point cloud data to a certain extent, and broadens the application field of lightweight and economical unmanned airborne LiDAR equipment. Description of the Drawings

[0057] Figure 1 is the flowchart of the work of the present invention;

[0058] Figure 2 is the data of the arched roof collected in the embodiment of the present invention;

[0059] Figure 3 is the data of the flat road surface collected in the embodiment of the present invention;

[0060] Figure 4 is the data of the forest land collected in the embodiment of the present invention;

[0061] Figure 5 is the schematic diagram of the KD-tree construction in the embodiment of the present invention;

[0062] Figure 6 is the schematic diagram of the search for neighboring points in the embodiment of the present invention;

[0063] Figure 7 is the schematic diagram of the projection of a point on a plane in the embodiment of the present invention;

[0064] Figure 8 is the effect diagram of filtering and optimizing the point cloud data collected for the house in the embodiment of the present invention;

[0065] Figure 9 is the effect diagram of filtering and optimizing the point cloud data collected for the flat road surface in the embodiment of the present invention;

[0066] Figure 10 is the effect diagram of filtering and optimizing the point cloud data collected for the forest land in the embodiment of the present invention. Detailed Embodiment

[0067] The present invention will be further described below with reference to the drawings.

[0068] As Figure 1 shown, an unmanned airborne LiDAR point cloud filtering and optimizing method includes the following steps:

[0069] S1. Establish an index of the point cloud using a KD tree;

[0070] S2. Calculate and determine the maximum search radius for identifying isolated points r ;

[0071] S3. Resample corresponding to the neighborhood threshold k 1 and the maximum search radius r ;

[0072] S4. Resample corresponding to the neighborhood threshold k 2 and the maximum search radius r ;

[0073] The experimental area of this embodiment is located in Huarong District, Ezhou City. The average drop of the survey area is 30 meters, the vegetation coverage rate is about 60%, there are residential areas, straight roads, artificially planted osmanthus forests and natural forests, and the slopes and embankments in the area are relatively obvious. Based on the airborne LiDAR point cloud data of the research area collected, the algorithm proposed by the present invention is used to perform filtering and optimization processing on it. During the operation of the algorithm, two parameters need to be determined, namely k 1 and k 2, which respectively represent the number of adjacent sampling points of the current sampling point. k The value range of 1 is usually selected from 8 to 12. k The value range of 2 is usually selected as 1.5 times of k 1, that is, 12 to 18;

[0074] In order to better analyze the performance of the Haida S2 laser and the characteristics of the point cloud data collected by it, under the condition of the same time period, the same survey area, and the same carrier (DJI M350), the Haida S2 set four flight altitudes of 80m, 100m, 120m, and 140m to collect four flights of LiDAR point cloud data, and the Riegl 1850 collected 1 flight of LiDAR point cloud data at a flight altitude of 150m. Superimpose the four flights of data corresponding to the Haida S2 at flight altitudes of 80m, 100m, 120m, and 140m (represented in yellow) and the 1 flight of data of the Riegl 1850 (represented in white), and intercept the cross-section in any case. The effect is as shown in Figure 2 in (a), (b) to Figure 4 in (a), (b);

[0075] The following processing is performed on the collected LiDAR point cloud data:

[0076] (1) Select the point with the x coordinate in the middle among all sampling points as the segmentation point for dimension division. As shown in Figure 5 , select as the segmentation point for dimension division. At this time, the sampling point That is the root node of the constructed KD tree. Based on this, the space where all sampling points are located can be divided into two sub-spaces, the left and the right. For all sampling points in the left space, their x coordinates are all less than , and for all sampling points in the right space, their coordinates are all greater than .

[0077] (2) For the sampling points in the left and right sub-spaces, select the y -axis direction of the sampling point coordinates for further dimensional division. Taking the left space as an example, for the sampling point 's y coordinate is in the middle of all sampling points in the left space. Therefore, select it as the splitting point for dimensional division. At this time, is added as the left child node of the KD tree. Based on this, all sampling points in the left space are divided into two parts, the upper part contains the points and ; the lower part contains the points , and . Similarly, is added as the right child node of the KD tree, and all sampling points in the left space are also divided into two parts;

[0078] (3) Repeat steps (1) to (2) until all sampling points are added to the KD tree, and the construction process of the KD tree ends;

[0079] (4) Based on the KD tree, determine the search path of the sampling point . As shown in Figure 6 , construct a stack P S for storing the nodes on the search path and set it to an empty stack; push onto the stack and start the search with the root node as the starting point. Considering that 's x coordinate is greater than , so search downward along the right subtree of the KD tree; push onto the stack and select the root node of the right subtree as the new starting point for the search. Considering that 's y coordinate is greater than , so continue to search downward along the right subtree of the right subtree; push onto the stack and use the node as the new starting point for the search. Considering that 's x coordinate is less than , so search downward along the left subtree; push onto the stack and combine Figure 6It can be seen that at this time, the found point is itself. Therefore, the search ends. At this time, the elements contained in the stack are ;

[0080] (5) If it is necessary to search for the nearest k ( k = 3) neighboring points, first set the neighboring point list vertList = {} and the distance list distList = {}, and perform the search according to the following steps: The stack P S is not empty, directly go to the next step; according to the principle of last in first out, pop the top element P S from , and calculate the distance between the current sampling point and . At this time, the distance is 0, it is judged that the target neighboring point coincides with the current sampling point, and enter the next step; pop the top element P S of , calculate the distance between the current sampling point and . At this time, since the length of vertList is less than k , therefore, directly add and to the corresponding lists. At this time, vertList = { }, distList = { }; Considering that is a non-leaf node, therefore, it is necessary to expand the search range: with as the center and as the radius to make a circle, there is an overlap with the hyperplane where the sampling point is located. Therefore, it is necessary to search for neighboring points from the other half space of the hyperplane where is located;

[0081] Push the nodes on the other branch of the KD tree corresponding to the point onto the stack according to the path search method respectively, and push onto the stack. At this time, ; Continue to pop the top element P S of , calculate the distance between the current sampling point and d 7-9 . At this time, since the length of vertList is less than k , therefore, directly add andd 7-9 Add it to the corresponding list. At this time, ; because it is a leaf node, so there is no need to expand the search scope; continue to pop the P S top element of the stack , calculate the distance between the current sampling point and . At this time, since the length of vertList is less than , so directly add k and to the corresponding list. At this time, ; because is not a leaf node, so it is necessary to increase the search scope: with as the center and as the radius to make a circle, there is an overlap with the hyperplane where the sampling point is located. Therefore, it is necessary to find neighboring points from the other half space of the hyperplane where is located; According to the method of determining the search path, it is necessary to push the

[0082] two nodes onto the stack. At this time, ; continue to pop the top element of the stack P . At this time, the length of vertList is equal to S . As shown in k , is greater than the maximum distance Figure 6 shown, is greater than the maximum distance , so it does not meet the neighboring point requirements; continue to pop the P S top element of the stack , calculate the distance between the current sampling point and it, . Therefore, it is necessary to update vertList and distList. The updated lists are respectively: ; because is not a leaf node, so it is necessary to increase the search scope: with as the center and as the radius to make a circle, there is an overlap with the hyperplane where the sampling point is located. Therefore, it is necessary to find neighboring points from the other half space of the hyperplane where

[0083] is located; According to the method of determining the search path, it is necessary to push the two nodes onto the stack. At this time,P S The top element of the stack . Calculate the current sampling point and the distance between . At this time, the length of vertList is equal to k , as Figure 6 shown is greater than the maximum distance . Therefore, it does not meet the requirements of neighboring points; continue to pop the stack P S The top element of the stack ; Calculate the current sampling point and the distance between . Therefore, it is necessary to update vertList and distList, and the updated lists are respectively: ; At this time, the stack P S is empty, and the neighboring point search process ends. The points adjacent to k ( k = 3) from near to far are respectively ;

[0084] (6) As Figure 7 shown, assume the coordinates of the current point are . Its projected point coordinates on the plane can be calculated according to the following steps: Calculate the directed distance from the point d to the plane: ; Calculate the projected point coordinates of the point on the plane : ;

[0085] (7) Expand the search range of the neighborhood and traverse each sampling point in again, and accordingly implement the secondary resampling processing of the point cloud data: Set the neighborhood set of the current sampling point to be empty ; Search for the neighboring points of in ascending order of distance, and handle them in the following two cases:

[0086] If and the Euclidean distance between them is less than r , and the number of points in the neighboring point set is less than k 2, add the neighboring point to the neighboring point set , continue to traverse the next sampling point;

[0087] If and the Euclidean distance between them is greater than r , stop searching;

[0088] Use all the points in the current point and the neighboring point set to fit the tangent plane, specifically as follows:

[0089] Calculate the centroid of all sampling points in the neighboring point set : ;

[0090] Construct the covariance matrix: ;

[0091] Calculate the eigenvalues and eigenvectors of the covariance matrix , select the eigenvector corresponding to the smallest eigenvalue as the normal of the tangent plane, denoted as ; Thus, the equation of the fitted plane is determined as: ;

[0092] where ;

[0093] Expand the above formula to get: ;

[0094] where ;

[0095] Based on the fitted tangent plane, project the current point onto the fitted tangent plane and assign the coordinates of the projection point to the current point as the new coordinates of the current point, specifically as follows:

[0096] Calculate the directed distance from the point d to the fitted tangent plane: ;

[0097] Calculate the coordinates of the projection point of the point on the fitted tangent plane : ;

[0098] After the traversal is completed, output the point cloud data generated by resampling to a file.

[0099] The HaiDa S2 airborne LiDAR device was used to collect field data in the study area (the average elevation difference in the study area is 30 meters, the vegetation coverage rate is about 60%, there are residential areas, straight roads, artificially planted osmanthus forests and natural forests, and the slopes and embankments in the area are relatively obvious). The data collected in four flights at flight altitudes of 80m, 100m, 120m, and 140m were superimposed. Then, the data after superposition was processed using the algorithm of the present invention, and the processed data was compared with the data collected by the Riegl 1850 laser sensor at a flight altitude of 150m under the same time period, the same survey area, and the same carrier (DJI M350). The results are as follows:

[0100] As Figure 8 shown in the residential area, (a) is the position of the house section, (b) is the point cloud at the section position, and (c) is the local enlarged effect. It can be seen from the figure that the point cloud collected by the HaiDa S2 laser sensor generally maintains a high consistency with the roof surface of the house, but the distribution of the point cloud along the normal direction of the plane where the roof is located is relatively discrete, with a thickness of about 15 - 20 cm, which has an adverse effect on the post-processing and further application of the airborne LiDAR point cloud data. After filtering and optimizing the point cloud data using the algorithm proposed by the present invention, on the one hand, the data maintains the consistency with the roof surface, and more importantly, the distribution of the point cloud along the normal direction of the plane where the roof is located becomes more concentrated, and the thickness is controlled within 5 cm.

[0101] As Figure 9 shown in the plane road surface area, (a) is the position of the straight road section, (b) is the point cloud at the section position, and (c) is the local enlarged effect. It can be seen from the figure that the point cloud of the straight road surface collected by the HaiDa S2 laser sensor generally maintains a high consistency with the sampling entity, but the distribution of the point cloud along the normal direction of the plane where the road surface is located is relatively discrete, with a thickness of about 15 - 20 cm. After filtering and optimizing the point cloud data using the algorithm of the present invention, the point cloud better maintains the consistency with the road surface morphology, and the distribution of the point cloud along the normal direction of the plane where the road surface is located becomes more concentrated, and the thickness is controlled within 5 cm.

[0102] As Figure 10The forest area shown, (a) is the location of the forest land profile, (b) is the point cloud at the profile location, and (c), (d) are the local magnification effects. Compared with the roof of the house and the flat road surface, the situation of the forest land is more complex, and the distribution of the point cloud in the vertical space is affected by many factors such as the density of vegetation cover, the growth of vegetation, and the characteristics of the vegetation itself. The ground point density of the point cloud collected by Heda S2 is relatively sparse compared with the roof of the house and the flat road surface, and the thickness of this part of the point cloud in the vertical direction is also relatively small. However, by comparing the distribution of the point cloud in the vertical direction before and after filtering optimization, it can be seen that the distribution of the point cloud after filtering treatment is more concentrated in the vertical direction than before, and the influence of random error is weakened to a great extent. After processing the LiDAR point cloud from the vegetation branches and the top surface with the algorithm proposed in the present invention, the distribution of the point cloud has changed slightly and maintains a high consistency with that before processing. It is worth mentioning that compared with the original sampled point cloud, when processing the point cloud data with the algorithm in the text, some isolated points are removed, and the sampling noise of the original point cloud data is suppressed to a certain extent.

[0103] In summary, when there is a large amount of noise in the point cloud data collected by the affordable airborne LiDAR sensor, filtering and optimizing it based on the algorithm of the present invention can effectively weaken the influence of sampling error on the point cloud data while effectively retaining the surface feature points of geographical entities, thereby providing data guarantee for the reconstruction of high-fidelity real-scene three-dimensional models and the construction of high-quality DEMs. The present invention can maximize the universality of airborne LiDAR point cloud data and expand the application fields of lightweight and affordable unmanned aerial vehicle-borne LiDAR devices.

Claims

1. An optimization method for filtering LiDAR point clouds carried by an unmanned aerial vehicle, characterized in that, It includes the following steps: S1. Establish an index of the point cloud using a KD tree. The specific process is as follows: S1.

1. Construct a KD tree: S1.1-1. Along the first dimension of the sampling point coordinates x Perform point division in this direction, select the point with the middle coordinate in the first dimension of all sampling points as the splitting point, and divide the space where the sampling points are located into two sub-spaces on the left and right; for all sampling points in the left space, their first dimension coordinates are less than the median point, and for all sampling points in the right space, their coordinates are greater than or equal to the median point; S1.1-2. For the sampling points in the left and right subspaces, select the second dimension y direction to re-divide the points. For the left and right spaces, respectively select the points with the second dimension coordinate in the middle as the splitting points, and re-divide the left and right spaces into two subspaces respectively; for all sampling points in the left subspace of the left space, their second dimension coordinates are less than the median point, and for all sampling points in the right subspace, their coordinates are greater than or equal to the median point; for all sampling points in the left subspace of the right space, their second dimension coordinates are less than the median point, and for all sampling points in the right subspace, their coordinates are greater than or equal to the median point; S1.1-3. For the four subspaces obtained by dividing S1.1-2, select the third dimension of the sampling points z and perform re-division of the points according to the method of S1.1-1 or S1.1-2 in the direction, further dividing the four subspaces into eight subspaces; S1.1-4. Repeat S1.1-1 to S1.1-3 until all sampled points are added to the KD tree, and the KD tree construction process ends. S1.2, k - Neighborhood search: Given the airborne LiDAR point cloud set , select any sampling point S in the set , and implement the search for the neighborhood points of the sampling point in two steps: Determine the search path and find neighboring points, specifically: S1.2-1. Determine the search path: S1.2-1-1. Construct a stack for storing the nodes on the search path P S and set it to an empty stack; S1.2-1-2. Using the root node of the KD tree as the starting point, push onto the stack, and compare with in terms of the magnitude in the first dimension of their coordinates. If , then use the root node of the left subtree as the new starting point to search downward. Conversely, if , then use the root node of the right subtree as the new starting point to search downward; S1.2-1-3, Judgment Check if it is empty. If it is empty, the search process ends; otherwise, push the new root node onto the stack , and compare with in terms of their sizes in the second dimension. If , use the root node of the left subtree as the new starting point to search downward; conversely, if , use the root node of the right subtree as the new starting point to search downward; S1.2-1-4, Judgment Check if it is empty. If it is empty, the search process ends; otherwise, push the new root node onto the stack and compare with the magnitudes of the coordinates in the third dimension. If , use the root node of the left subtree as the new starting point to search downward; conversely, if , use the root node of the right subtree as the new starting point to search downward; S1.2-1-5. Repeat S1.2-1-2 to S1.2-1-4 until the search path encounters a leaf node of the KD tree. S1.2-2. Find adjacent points: Assume the current sampling point is . After determining the adjacent point search path, set the adjacent point list vertList = {} and the distance list distList = {} to store the adjacent points and their corresponding distances respectively. Then, based on the backtracking method, determine the k adjacent points of the current sampling point. The steps are as follows: S1.2-2-1, Judgment Stack Determine whether it is an empty stack. If so, the adjacent point search process ends; otherwise, proceed to S1.2-2-2; S1.2-2-2. Pop the top element from in sequence according to the last-in-first-out principle , calculate the distance between and . If = 0, consider and to be the same point, go back to S1.2-2-1; otherwise, enter S1.2-2-3. S1.2-2-3. If the length of vertList is less than k , directly add and to the end of the corresponding list, that is, add to the list vertList, add to the list distList; if the length of vertList is equal to k , judge and the maximum distance stored in the list distList for comparison. If , then delete from distList and delete the corresponding neighboring point from vertList, and insert into the corresponding position of the distList list in ascending order of the neighboring point distance, and insert into the corresponding position of vertList, ensuring that the subscript of in vertList is the same as the subscript in distList; S1.2-2-4, Judgment Whether it is a leaf node. If it is, go back to S1.2-2-1. If not, go to S1.2-2-5; S1.2-2-5. With as the center and as the radius, draw a circle and determine whether there is an overlap between the circle and the hyperplane where is located. If there is an overlap, then perform a neighboring point search from the other half space of the hyperplane where is located. According to the path search method, use as the starting point, determine the search path according to the method described in S1.2-1, and add all the nodes on the path to the stack and return to S1.2-2-1; S2. Calculate and determine the maximum search radius for identifying isolated points r ; The specific process is as follows: S2.

1. Calculate the maximum search radius: Based on the constructed KD tree, select all or a part of the airborne LiDAR point cloud set as a sample for determining the maximum search radius r , denoted as ; Traverse S r each sampling point in , search in ascending order of distance for the corresponding k 2 neighboring points, record this neighboring point the distance between it and the current sampling point is r i , after all sampling points in the sample S r are traversed, take the average value of the maximum search radius of all traversed sampling points r i as the maximum search radius of the entire point cloud r , the calculation formula is as follows: ; S2.

2. Identify and remove isolated points: Based on the constructed point cloud KD-tree index, traverse each sampling point in the airborne LiDAR point cloud set in the set , set the neighboring point set of the current sampling point to be empty, add the sampling point to the neighboring point set , and perform the following processing: S2.2-1. Search for the neighboring points of in ascending order of distance, and handle them in the following two cases: ​ S2.2-1-1. If and the Euclidean distance between them is less than r and the number of points in the neighboring point set is less than k 1, add the neighboring point to the neighboring point set , and continue to traverse the next sampling point; S2.2-1-2. If and the Euclidean distance between them is greater than r , stop searching; S2.2-2. Statistical adjacent point set Count the number of midpoints. If the number is less than 3, that is, the number of sampling points in the adjacent point set cannot meet the requirements of plane fitting, then it is considered to be an isolated point, delete it from the original point cloud, and traverse the next sampling point; S3. Corresponding to the neighborhood threshold k 1 and the maximum search radius r Perform resampling. The process is to traverse each sampling point in , take the first adjacent points from the set of adjacent points of the current sampling point , denoted as , and fit the corresponding tangent plane as follows: S3.

1. Calculate adjacent points The centroid of all sampling points in the set: ; S3.

2. Construct the covariance matrix: ; S3.

3. Calculate the covariance matrix of the eigenvalues and eigenvectors , select the eigenvector corresponding to the smallest eigenvalue as the normal of the tangent plane, denoted as ; thus, the equation of the fitting plane is determined as: ; Among them, ; Expanding the above formula gives: ; Among them, ; S3.

4. Project the current sampling point onto the fitted tangent plane, and assign the coordinates of the projected point to the current point as the new coordinates of the current point, as follows: S3.4-1, Calculation Point Directed distance to the fitting tangent plane d : ; S3.4-2, Calculation Point On the fitted tangent plane Projection point coordinates: ； S4. Corresponding to the neighborhood threshold k 2 and the maximum search radius r Perform resampling. The process is to expand the search range of the neighborhood and traverse again each sampling point in , and accordingly implement the secondary resampling process of the point cloud data, specifically as follows: S4.

1. Set the current sampling point whose neighborhood set is empty; S4.

2. Search in ascending order of distance for adjacent points , and handle in the following two cases: S4.2-1. If and the Euclidean distance between them is less than r , and the number of points in the set of neighboring points is less than k 2, add the neighboring point to the set of neighboring points , and continue to traverse the next sampling point; S4.2-2. If and the Euclidean distance between them is greater than r , stop searching; S4.

3. Using the method described in S3.1~S3.3, fit the tangent plane with all the points in the current point and the adjacent point set ; S4.

4. Based on the method described in S3.4, project the current point onto the fitted tangent plane and assign the coordinates of the projected point to the current point as the new coordinates of the current point. S4.

5. After the traversal is completed, output the processed point cloud data to a file.

Images