Unmanned aerial vehicle-mounted LiDAR point cloud filtering optimization method

Abstract

The invention discloses an unmanned aerial vehicle-mounted LiDAR point cloud filtering optimization method. The method comprises the following steps: establishing an index of a point cloud by using a KD tree; calculating and determining a maximum search radius r for identifying isolated points; performing resampling corresponding to the neighborhood threshold k1 and the maximum search radius r; and resampling corresponding to the neighborhood threshold k2 and the maximum search radius r. An index of a point cloud is established by using a KD tree, and identification and elimination of isolated points are realized based on searching of adjacent points; local tangent plane fitting at each sampling point position is realized by using a moving least square method, and resampling of point cloud is realized based on a projection transformation method, so that filtering optimization processing of original sampling data is realized; the method greatly weakens the influence of the sampling error on the point cloud data, improves the consistency between the obtained data and the collected point cloud data, improves the universality of the airborne LiDAR point cloud data, and widens the application field of lightweight, economical and applicable unmanned aerial vehicle-mounted LiDAR equipment.

Description

Technical Field

[0001] The present invention relates to a method for optimizing the filtering of LiDAR point clouds carried by an unmanned aerial vehicle, belonging 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 feature points of ground objects and landforms. 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, etc. Among them, airborne LiDAR has received 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 increasingly developing 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 a method for optimizing the filtering of 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 a method for optimizing the filtering of UAV-borne LiDAR point clouds, including the following steps: S1. Establish an index of the point cloud using a KD tree; S2. Calculate and determine the maximum search radius for identifying isolated points r ; S3. Corresponding to the neighborhood threshold k 1 and the maximum search radius rPerform resampling; S4. Corresponding to the neighborhood threshold k 2 and the maximum search radius r Perform resampling.

[0006] Furthermore, the specific process of S1 is as follows: S1.1. Construct a KD tree: 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 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 sub-spaces, select the second dimension of the sampling points y direction for re-partitioning the points. For the left and right spaces, respectively select the point with the middle coordinate in the second dimension as the splitting point, and re-partition the left and right spaces into two sub-spaces respectively; for all sampling points in the left sub-space of the left space, their second dimension coordinates are less than the median point, and for all sampling points in the right sub-space, their coordinates are greater than or equal to the median point; for all sampling points in the left sub-space of the right space, their second dimension coordinates are less than the median point, and for all sampling points in the right sub-space, their coordinates are greater than or equal to the median point; S1.1-3. For the four sub-spaces obtained by the partitioning in S1.1-2, select the third dimension of the sampling points z direction and re-partition the points according to the method of S1.1-1 or S1.1-2, and further divide the four sub-spaces into eight sub-spaces; 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; 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: S1.2-1. Determine the search path: S1.2-1-1. Construct a stack P for storing the nodes on the search path and set it to an empty stack; S1.2-1-2. Take the root node of the KD tree as the starting point, push into the stack, and compare with The size of the coordinates of in the first dimension. If , then take the root node of the left subtree as the new starting point and search downward. Conversely, if S1.2-1-3. Judge whether it is empty. If it is empty, the search process ends. Otherwise, push the new root node onto the stack, and compare with in the second dimension of the coordinates. If , then take the root node of the left subtree as the new starting point and search downward. Conversely, if , then take the root node of the right subtree as the new starting point and search downward; S1.2-1-4. Judge whether it is empty. If it is empty, the search process ends. Otherwise, push the new root node onto the stack, and compare with in the third dimension of the coordinates. If , then take the root node of the left subtree as the new starting point and search downward. Conversely, if , then take the root node of the right subtree as the new starting point and search downward; S1.2-1-5. Repeat S1.2-1-2 to S1.2-1-4 until the search path encounters the leaf node of the KD tree; S1.2-2. Search for the nearest neighbor: Assume the current sampling point is . After determining the search path for the nearest neighbor, set up the nearest neighbor list vertList={} and the distance list distList={} to store the nearest neighbors and their corresponding distances respectively. Then, based on the backtracking method, determine the k nearest neighbors of the current sampling point. The steps are as follows: S1.2-2-1. Judge whether the stack is an empty stack. If it is, the nearest neighbor search process ends. Otherwise, go to S1.2-2-2; S1.2-2-2. According to the principle of last in first out, pop the top element from in turn, calculate the distance between and . If =0, consider and to be the same point, and go back to S1.2-2-1. Otherwise, go to 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, is added to the vertList, is added to the distList; if the length of vertList is equal to k , judge and the maximum distance stored in the distList for comparison. If , then delete from the distList and delete 's corresponding neighboring point from the vertList. Insert into the corresponding position of the distList in ascending order of the neighboring point distance, and insert into the corresponding position of the vertList to ensure that 's subscript in the vertList is the same as its subscript in the distList; S1.2-2-4. Judge whether is a leaf node. If it is, go back to S1.2-2-1. If not, enter S1.2-2-5; S1.2-2-5. Take as the center and as the radius to make a circle, and judge whether there is an overlap between the circle and the hyperplane where is located. If there is an overlap, then perform neighboring point search in 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 go back to S1.2-2-1.

[0007] Further, the specific process of S2 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 the sample for determining the maximum search radius r , denoted as ; traverse each sampling point S r in , and search for the corresponding to the k2 adjacent points, record this adjacent point The distance between and the current sampling point is r i , waiting for the sample S r After all the sampling points in are traversed, take the maximum search radius of all the traversed sampling points r i The average value of is used as the maximum search radius of the entire point cloud r , and the calculation formula is as follows: ; S2.2. Isolated point identification and removal: Based on the constructed point cloud KD tree index, traverse the airborne LiDAR point cloud set Each sampling point in , set the adjacent point set of the current sampling point to be empty, and add the sampling point to the adjacent point set and do the following processing: S2.2-1. Search for the adjacent points of in ascending order of distance, and deal with them in the following two cases: The adjacent points of , and deal with them in the following two cases: S2.2-1-1. If the Euclidean distance between and is less than r and the number of points in the adjacent point set is less than k 1, add the adjacent point to the adjacent point set , and continue to traverse the next sampling point; S2.2-1-2. If the Euclidean distance between and is greater than r , stop searching; S2.2-2. Count the number of points in the adjacent point set . 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 that is an isolated point, delete it from the original point cloud, and traverse the next sampling point.

[0008] Furthermore, the process of S3 is to traverse each sampling point in , take the first adjacent points from the adjacent point set of the current sampling point , denoted as , and fit the corresponding tangent plane, specifically 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 eigenvalues of the covariance matrix and the 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: ; where, ; Expand the above formula to get: ; where, ; S3.4. Project the current sampling 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: S3.4-1. Calculate the directed distance from the point d to the fitted tangent plane: ; S3.4-2. Calculate the coordinates of the projection point of the point on the fitted tangent plane : .

[0009] 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, specifically as follows: S4.1. Set the neighborhood set of the current sampling point to be empty; S4.2. Search for the adjacent points of in ascending order of distance, and handle them in the following two cases: S4.2-1. If the Euclidean distance between and is less than r and the number of points in the adjacent point set is less than k 2, add the adjacent point to the adjacent point set , and continue to traverse the next sampling point; S4.2-2. If and The Euclidean distance between r is greater than S4.3. Based on the methods described in S3.1~S3.3, use all the points in the current point and the set of neighboring points to fit the tangent plane; 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.

[0010] Through the use of a KD tree to establish an index of the point cloud, the present invention realizes the identification and elimination of isolated points based on neighboring point search; uses the moving least squares method to realize the local tangent plane fitting 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, to a certain extent improves the universality of the airborne LiDAR point cloud data, and broadens the application field of lightweight and economical unmanned aerial vehicle borne LiDAR equipment. BRIEF DESCRIPTION OF THE DRAWINGS

[0011] Figure 1 is the flowchart of the work of the present invention; Figure 2 is the data of the arched roof collected in the embodiment of the present invention; Figure 3 is the data of the flat road surface collected in the embodiment of the present invention; Figure 4 is the data of the forest land collected in the embodiment of the present invention; Figure 5 is the schematic diagram of the KD tree construction in the embodiment of the present invention; Figure 6 is the schematic diagram of neighboring point search in the embodiment of the present invention; Figure 7 is the schematic diagram of the projection of a point on a plane in the embodiment of the present invention; 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; 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; 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 DESCRIPTION OF THE INVENTION

[0012] The present invention will be further described below in conjunction with the accompanying drawings.

[0013] As Figure 1 shown, a method for optimizing UAV-borne LiDAR point cloud filtering includes the following steps: S1. Establish an index of the point cloud using a KD tree; S2. Calculate and determine the maximum search radius for identifying isolated points r ; S3. Resample corresponding to the neighborhood threshold k 1 and the maximum search radius r ; S4. Resample corresponding to the neighborhood threshold k 2 and the maximum search radius r ;

[0014] The experimental area of this embodiment is located in Huarong District, Ezhou City. The average drop in the surveyed 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 collected airborne LiDAR point cloud data of the research area, 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 k 1 is usually selected from 8 to 12, k The value range of 2 is usually selected as 1.5 times of Figure 2 1, that is, 12 to 18; Figure 4 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 surveyed 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. The four flights of data corresponding to the Haida S2 at 80m, 100m, 120m, and 140m flight altitudes (represented in yellow) and the 1 flight of data of the Riegl 1850 (represented in white) are superimposed, and the cross-section in any case is intercepted, and the effect is as shown in (a), (b) in (1) Select the point with the x coordinate in the middle among all sampling points as the segmentation point for dimension division. As Figure 5 shown, select As the splitting point for dimensional division, at this time, the sampling point 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 less than , and for all sampling points in the right space, their coordinates are greater than ; (2) For the sampling points in the left and right sub-spaces, select the y axis direction of the sampling point coordinates for re-division of the dimension. Taking the left space as an example, the coordinate of the sampling point y 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 points and ; the lower part contains 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; (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; (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 , 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 , 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 , 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 , combined with Figure 6 it can be seen that at this time, the point found is itself. Therefore, the search ends. At this time, the elements contained in the stack are ; (5) If it is necessary to search for the 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 the , and calculate the distance between the current sampling point and . At this time, the distance is 0, judge that the target neighboring point coincides with the current sampling point, and enter the next step; Pop the top element P S of the , 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 increase the search range: With as the center and as the radius to make a circle, there is a coincidence between the circle and 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; 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 the , 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 and d 7-9 to the corresponding lists. At this time,​​ ; Since is a leaf node, there is no need to expand the search range; 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 , directly add k and to the corresponding list. At this time, ; Because is not a leaf node, it is necessary to increase the search range: with as the center and as the radius to draw 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; According to the method of determining the search path, it is necessary to push the two nodes onto the stack. At this time, and ; continue to pop the top element of the stack P S . At this time, the length of vertList is equal to , as shown in k , Figure 6 is greater than the maximum distance , so it does not meet the requirements of neighboring points; continue to pop the top element of the stack P S and calculate the distance between the current sampling point and it. , so it is necessary to update vertList and distList. The updated lists are respectively: ; Since is not a leaf node, it is necessary to increase the search range: with as the center and as the radius to draw 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; According to the method of determining the search path, it is necessary to push the two nodes onto the stack. At this time, ; continue to pop the P S 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 , so it does not meet the requirements of neighboring points; continue to pop the stack P S the top element of the stack ; calculate the distance between the current sampling point and the distance between , so, it is necessary to update vertList and distList, and the updated lists are respectively: ; at this time the stack P S is empty, the neighboring point search process ends, and the points neighboring are k ( k = 3) from near to far are respectively ; (6) As Figure 7 shown, assume the coordinates of the current point are , and its projected point coordinates on the plane can be calculated as follows: calculate the directed distance from the point d to the plane: ; calculate the projected point coordinates of the point on the plane : ; (7) 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: 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: If and the Euclidean distance between 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; If and the Euclidean distance between is greater than r , stop searching; Using the current point and the set of neighboring points to fit the tangent plane with all the points therein, specifically as follows: Calculate the centroid of all the sampling points in the set of neighboring points : ; Construct the covariance matrix: ; Calculate the eigenvalues and eigenvectors of the covariance matrix , and 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: ; wherein, ; Expand the above formula to obtain: ; wherein, ; Based on the fitted tangent plane, 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, specifically as follows: Calculate the directed distance from the point d to the fitted tangent plane: ; Calculate the coordinates of the projected point of the point on the fitted tangent plane : ; After the traversal is completed, output the resampled point cloud data to a file.

[0015] The research area (the average drop of the research 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.) was used for field data collection with the Hydra S2 airborne LiDAR device, and the data of four flight runs collected at flight altitudes of 80m, 100m, 120m, and 140m were superimposed. Then, the data after superimposition 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: As Figure 8The residential area shown, (a) is the section position of the house, (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 Heda S2 laser sensor generally maintains a high consistency with the top 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 top surface of the house, 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.

[0016] As Figure 9 The plane road surface area shown, (a) is the section position of the straight road surface, (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 Heda 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 along the normal direction of the plane where the road surface is located becomes more concentrated, and the thickness is controlled within 5 cm.

[0017] As Figure 10 The forest area shown, (a) is the section position of the forest land, (b) is the point cloud at the section position, and (c), (d) are the local enlarged effects. Compared with the top surface of the house and the straight 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 density of ground points of the point cloud collected by Heda S2 is relatively sparse compared with the top surface of the house and the plane 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 and optimization, it can be seen that the distribution of the filtered point cloud in the vertical direction is more concentrated than before, and the influence of random errors is weakened to a great extent. For the LiDAR point cloud from the vegetation branches and the top surface, after processing using the algorithm proposed by 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 using 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.

[0018] In summary, when there is significant noise in the point cloud data collected by an affordable airborne LiDAR sensor, filtering and optimizing it based on the algorithm of the present invention can effectively weaken the influence of sampling errors 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 3D 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 (UAV)-borne LiDAR devices.

Claims

1. A method for optimizing the filtering of UAV-borne LiDAR point cloud, characterized in that: The steps include: S1. Use KD tree to build the point cloud index; S2. Calculate and determine the maximum search radius for identifying isolated points r ; S3, corresponding to the neighborhood threshold k 1 and the maximum search radius r Perform resampling; S4, corresponding to the neighborhood threshold k 2 and the maximum search radius r Resample.

2. The UAV-mounted LiDAR point cloud filtering optimization method according to claim 1, characterized in that: The specific process of S1 is as follows: S1.

1. Construct KD tree: S1.1-1, along the first dimension of the sampling point coordinates x The points are divided in the direction, and the point with the first dimension coordinate of all sampling points in the middle is selected as the segmentation point, and the space where the sampling points are located is divided into two subspaces on the left and right. The first dimension coordinates of all sampling points in the left space are smaller than the median point, and the coordinates of all sampling points in the right space 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 of the sampling points y The points are divided again in the direction. For the left and right spaces, the points with the second dimension coordinates in the middle are selected as the split points, and the left and right spaces are divided into two subspaces respectively; the second dimension coordinates of all sampling points in the left subspace in the left space are less than the median point, and the coordinates of all sampling points in the right subspace are greater than or equal to the median point; the second dimension coordinates of all sampling points in the left subspace in the right space are less than the median point, and the coordinates of all sampling points in the right subspace are greater than or equal to the median point; S1.1-3, for the four subspaces obtained by S1.1-2, select the third dimension of the sampling point z Direction: The points are divided again according to the method of S1.1-1 or S1.1-2, and the four subspaces are further divided into eight subspaces; 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 is completed; S1.2, k -Neighborhood search: Given a collection of airborne LiDAR point clouds , select the collection S Any sampling point in , the sampling point is realized in two steps Searching for neighboring points: Determine the search path and search for neighboring points, specifically: S1.2-1. Determine the search path: S1.2-1-1. Construct a stack for storing nodes on the search path P S , and set it to an empty stack; S1.2-1-2, the root node of the KD tree As a starting point, Push into stack and compare and The size of the coordinates in the first dimension is , then the root node of the left subtree is used as the new starting point to search downwards. Conversely, if , then search downward along the root node of the right subtree as the new starting point; S1.2-1-3. Judgment Is it empty? If it is empty, the search process ends. Otherwise, a new root node is pushed into the stack. , and compare and The size of the coordinates in the second dimension is , then the root node of the left subtree is used as the new starting point to search downwards. Conversely, if , then search downward along the root node of the right subtree as the new starting point; S1.2-1-4. Judgment Is it empty? If it is empty, the search process ends. Otherwise, a new root node is pushed into the stack. , and compare and The size of the coordinates of in the third dimension is , then the root node of the left subtree is used as the new starting point to search downwards. Conversely, if , then search downward along the root node of the right subtree as the new starting point; 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 the adjacent point: Assume that the current sampling point is After determining the neighboring point search path, set the neighboring point list vertList={} and the distance list distList={} to store the neighboring points and their corresponding distances respectively, and then determine the current sampling point based on the backtracking method. k The steps are as follows: S1.2-2-1, judgment stack Is the stack empty? If so, the neighboring point search process ends. Otherwise, go to S1.2-2-2. S1.2-2-2, follow the principle of last in, first out Pop the top elements of the stack one by one ,calculate and The distance between ,if =0, think and If it is the same point, go back to S1.2-2-1, otherwise, go to S1.2-2-3; S1.2-2-3, if the length of vertList is less than k , directly and Add 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 The maximum distance stored in the list distList For comparison, if , then delete it from distList And remove from vertList The corresponding neighboring points are arranged in order from small to large according to the distance between the neighboring points. Insert into the corresponding position of the distList list and Insert into the corresponding position of vertList, make sure Subscripts in vertList Same subscript as in distList; S1.2-2-4, Judgment Is it a leaf node? If yes, go back to S1.2-2-1; if not, go to S1.2-2-5; S1.2-2-5, is the center of the circle, Draw a circle with radius Is there any overlap between the hyperplanes? If so, The other half of the hyperplane space is used to search for neighboring points, and the path search method is used to find the neighboring points. As a starting point, determine the search path according to the method described in S1.2-1 and add all nodes on the path to the stack and return to S1.2-2-1.

3. The UAV-borne LiDAR point cloud filtering optimization method according to claim 1, characterized in that: The specific process of S2 is: S2.

1. Calculate the maximum search radius: Based on the constructed KD tree, select the airborne LiDAR point cloud collection All or part of it is used to determine the maximum search radius r The sample is denoted as ; Traversal S r Each sampling point in , search in order of distance from small to large The corresponding k 2 adjacent points, record this adjacent point The distance from the current sampling point is r i , waiting for sample S r After all sampling points in are traversed, the maximum search radius of all traversed sampling points is taken r i The average value is taken as the maximum search radius of the entire point cloud r , the calculation formula is as follows: ; S2.2, Identification and elimination of isolated points: Based on the constructed point cloud KD tree index, traverse the airborne LiDAR point cloud collection Each sampling point in , set the current sampling point The set of neighboring points If it is empty, the sampling point Add to neighboring point set And do the following processing: S2.2-1. Search in order of distance from smallest to largest Neighboring points , divided into the following two situations: S2.2-1-1, if and The Euclidean distance between r , and the set of neighboring points The number of midpoints is less than k 1, the adjacent point Add to neighboring point set , continue to traverse the next sampling point; S2.2-1-2, if and The Euclidean distance between r , stop searching; S2.2-2. Counting neighboring point sets The number of midpoints. If the number is less than 3, that is, the number of sampling points in the neighboring point set cannot meet the requirements of plane fitting, then it is considered If it is an isolated point, delete it from the original point cloud and traverse the next sampling point.

4. The UAV-borne LiDAR point cloud filtering optimization method according to claim 1, characterized in that: The S3 process is to traverse Each sampling point in , from the current sampling point Take the previous point from the neighboring point set neighboring points, denoted as , fitting the corresponding tangent plane, as follows: S3.

1. Calculate neighboring points The centroid of all sample points in the set: ; S3.

2. Construct the covariance matrix: ; S3.

3. Calculate the covariance matrix The eigenvalue of With the eigenvector , select the eigenvector corresponding to the minimum eigenvalue as the normal of the tangent plane, denoted as ; Thus, the equation of the fitting plane is determined as: ; in, ; Expand the above formula, and we get: ; in, ; S3.4, the current sampling point Project to the fitted tangent plane and convert the projection point coordinates Assign values ​​to the current point as the new coordinates of the current point, as follows: S3.4-1. Calculation points Signed distance to the fitted tangent plane d : ; S3.4-2. Calculation points In the fitting tangent plane The coordinates of the projected point on : .

5. The UAV-borne LiDAR point cloud filtering optimization method according to claim 1, characterized in that: The process of S4 is to expand the search range of the neighborhood and traverse again Each sampling point in , and based on this, the secondary resampling of point cloud data is realized, as follows: S4.

1. Set the current sampling point Neighborhood set of is empty; S4.

2. Search in order of distance from smallest to largest Neighboring points , divided into the following two situations: S4.2-1, if and The Euclidean distance between r , and the set of neighboring points The number of midpoints is less than k 2, the adjacent point Add to neighboring point set , continue to traverse the next sampling point; S4.2-2, if and The Euclidean distance between r , then stop searching; S4.3, based on the method described in S3.1~S3.3, using the current point and the set of neighboring points Fitting of tangent planes to all points in 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, the processed point cloud data is output to a file.

Images