Method and system for removing airborne LiDAR noise caused by building glass

Abstract

The present invention discloses a method and system for removing airborne LiDAR noise caused by building glass. The noise removal method includes: 1. Segmenting the airborne LiDAR point cloud data to obtain a ground point set and a non-ground point set; 2. Segmenting the non-ground point set to obtain an underground noise point set and a remaining point set; 3. Segmenting the remaining point set to obtain a high point set and an unclassified point set; 4. Segmenting the high point set to obtain a building point clustering set, and the remaining points are noise points on the horizontal periphery of the building top; 5. Extracting a building roof patch point set and a non-roof patch point set from the building point clustering set; 6. Extracting high-noise points from the non-roof patch point set; 7. Extracting noise points within the building plane coverage range and with an elevation close to the ground from the unclassified point set; 8. Extracting elevation outlier points in the non-noise point set as noise points. This method can remove airborne LiDAR noise caused by building exterior wall glass.

Description

Technical Field

[0001] The present invention belongs to the technical field of airborne lidar point cloud data processing, and particularly relates to a method and system for removing airborne LiDAR noise caused by exterior wall glass of buildings in a survey area. Background Technique

[0002] In general, in the hardware design of lidar (Light Detection And Ranging, LiDAR), the transmitter and receiver are at the same station, that is, the positions of the transmitter and receiver are close in space. When the detection distance is far, it can be considered that the transmitter and receiver are at the same position. On this basis, the basic assumption of LiDAR ranging is that both the transmitted pulse and the received pulse are transmitted along the same straight line, but the transmission directions of the transmitted pulse and the received pulse are opposite, and the time difference (Time-of-flight) between the transmitted pulse and the laser pulse is used to measure the distance.

[0003] Glass is one of the commonly used materials in modern buildings. On the exterior walls of buildings, glass windows and glass curtain walls are very common, and the surface normal vectors are close to being parallel to the ground. The working wavelength of LiDAR is generally in the micrometer or sub-micrometer range, and the roughness of glass is also in this range. When a laser pulse is incident on glass, in addition to transmission, strong specular reflection or diffuse reflection will also occur. The scattered / reflected laser pulse may also encounter other scattering or reflecting surfaces, resulting in multiple scattering or reflection, thus causing the non-linear transmission of the laser pulse. This causes the laser receiver to finally receive not the echo pulse generated at the first scattering point / reflection point, and the distance calculated from the time difference between the received pulse and the transmitted pulse is greater than the distance between the laser and the first reflection point / scattering point. For example, in Figure 1 , the transmitted pulse is emitted from point P, reaches the first reflection point P wall on the exterior wall glass, and is reflected at P wall , then reaches the second reflection point / scattering point P grd on the ground, and another reflection or scattering occurs. In Figure 1 (a), a part of the pulse reflected or scattered at P grd is incident on P wall , and is reflected again and received by the laser; in Figure 1 (b), a part of the pulse reflected or scattered at P grd is directly received by the laser. Finally, the laser receives the pulse with a longer transmission path after multiple reflections or scatterings. The distance d wall from P wall to the laser is =|PP wall |, which is significantly less than the actual transmission distance d all of the received pulse (for Figure 1 (a), it is 2|PPwall |+2|P wall P grd |; For Figure 1 (b), it is |PP wall |+|P wall P grd |+|P grd P|) of half. However, when calculating the distance between the laser and the scattering surface, the available one is d all , not d wall . If the received pulse corresponds to the correct transmitted pulse, the scattering point / reflection point calculated based on the straight-line transmission assumption of the laser pulse is at a false position P0 instead of P wall . Therefore, at the outer wall of a building using glass materials, noise is prone to occur itself.

[0004] The adoption of multi-pulse technology exacerbates the above noise problem. The main goal of multi-pulse technology is to increase the pulse repetition frequency of LiDAR, so it allows multiple transmitted pulses to exist simultaneously in the air. The names of multi-pulse technologies from different manufacturers are different, including FMP (Fixed Multipulse), CMP (Continuous Multipulse) of Optech, MPiA (Multiple Pulse in Air) of Leica, MTA (Multiple-Time-Around) of Rigel, etc. For example, the Riegl VQ-1560i using MTA technology supports up to 20 laser pulses in the air simultaneously. However, multi-pulse technology also brings another problem, that is, how to determine the correspondence between the received pulse and the transmitted pulse (see Figure 2 ), which directly involves the calculation of the pulse flight time and the measurement of the distance. Figure 2 is a comparison schematic diagram of single-pulse and multi-pulse. In Figure 2 , T1 and T2 are transmitted pulses, R1 and R2 are received pulses, the horizontal axis is the time axis, recording the time of pulse transmission and reception. The solid line indicates that the transmitted pulse corresponding to the received pulse can be determined, and the dashed line indicates that the transmitted pulse corresponding to the received pulse cannot be determined. Figure 2 (a) is the case of using single-pulse technology, Figure 2 (b) is the case of using multi-pulse technology. Currently, the main solution of hardware manufacturers to this problem is to judge the transmitted pulse based on the neighborhood continuity assumption and use technologies such as variable-period measurement to amplify the distance discontinuity caused by incorrect judgment of the transmitted pulse.

[0005] Assume that due to non-straight-line transmission, the time for the received pulse to reach the receiver is ΔT′ later than the time for the pulse backscattered at P wall to reach the receiver. In the case of correct judgment of the transmitted pulse, the corresponding reflection point / scattering point of the former should be from Pwall The position of the transmitted pulse is cΔT′ / 2 distance along the transmission direction, that is, P0, where c is the speed of light. If the transmitted pulse is calculated incorrectly and the selected transmitted pulse is emitted earlier than the correct transmitted pulse, the reflection point / scattering point P - It will be lower, possibly lower than the ground elevation; if the transmit pulse is calculated incorrectly and the selected transmit pulse is later than the correct transmit pulse, then the position of the reflection point / scattering point P + It will be higher, and may be close to or even exceed the roof elevation of the upper floors.

[0006] However, building exteriors often have large areas not covered by glass. Balconies, outdoor air conditioners, and other structures are interspersed with glass. These walls and structures utilize a large number of non-glass materials, such as dry-hanging stone, paint, natural stone paint, and thermal insulation panels. Compared to glass, these materials have a higher roughness and are less susceptible to strong specular reflections. Consequently, there are many non-noise points on building exteriors. Figure 3 This diagram shows the distribution of noise and non-noise points on a building's exterior wall. Rectangular blocks represent tall buildings, solid horizontal lines represent the ground, solid black points represent reflection / scattering points on the exterior wall, and white points represent noise points. Of the four closely located points P1, P2, P3, and P4, P2 and P4 fall on the glass, generating non-linear laser pulses, while P1 and P3 experience only normal backscattering. Therefore, if the transmitted pulse corresponding to the echo pulses of P1 and P3 is correctly determined, the positions of P1 and P3 are accurate. However, for P2 and P4, even if the transmitted pulse corresponding to their echoes is correctly determined, the time difference between the transmitted and received pulses is significantly greater than that of P1 and P3. The noise points located close to the ground are the positions of P2 and P4 inferred when the transmitted pulse is correctly determined. The noise points located significantly above the ground are the positions of P2 and P4 inferred when the transmitted pulse is incorrectly selected. The noise points located significantly below the ground are the positions of P2 and P4 inferred when the transmitted pulse is incorrectly selected. Because the echoes from P2 and P4 are not isolated, and there are many similar echoes, all received by the LiDAR at similar times, the transmitted pulse judgment algorithm does not clearly violate the neighborhood continuity assumption and cannot be treated as noise removal. In this case, if the received pulse is assigned a transmitted pulse that is earlier than the correct transmitted pulse, the calculated point position may be significantly below the ground. If the received pulse is assigned a transmitted pulse that is later than the correct transmitted pulse, the calculated point position will be significantly higher than the height of the glass where the reflection occurred, approaching or even exceeding the height of the roof.

[0007] In summary, at the exterior wall glass of a building, due to multiple scattering / reflections of laser pulses resulting in non-linear transmission, and incorrect judgment of superimposed emission pulses, the airborne LiDAR point cloud noise may be distributed either below the building or even below the ground without building coverage (the elevation may be close to the ground or significantly lower than the ground), or above the building top and near the roof.

[0008] Since it is difficult to determine whether non-linear transmission has occurred for the received pulses, it is very difficult to handle this type of noise at the signal processing level. In data post-processing, since this type of noise is not consistent with the assumptions for discrete isolated noise in traditional point cloud processing, there is no ready-made noise reduction method. The reflection noise processed in Chinese Patent Application CN202110204990.0 is caused by a glass plate installed outside the LiDAR hardware, and is not the same problem as that solved in this application. Summary of the Invention

[0009] Object of the Invention: The present invention provides a method and system for removing airborne LiDAR noise caused by building glass, aiming to remove the airborne LiDAR noise caused by the exterior wall glass of a building.

[0010] Technical Solution: On the one hand, the present invention discloses a method for removing airborne LiDAR noise caused by building glass, including the following steps:

[0011] S1. Use the three-dimensional DBSCAN algorithm to segment the airborne LiDAR point cloud data S all to obtain a ground point clustering set C grd , and all the points included in C grd constitute the ground point set S grd , remove the points in S all from S grd to obtain a non-ground point set S nongrd ;

[0012] S2. Extract the points from the non-ground point set S nongrd whose elevation is lower than the ground and the elevation difference from the ground belongs to the first elevation difference range to obtain an underground noise point set S lownoise , and add S lownoise to the noise point set S noise ; S remain =S nongrd -S noise is the remaining point set; the first elevation difference range is [Δz low1 , Δz up1 , Δz low1 is the lower limit of the first elevation difference, and Δz up1 is the upper limit of the first elevation difference, both of which are negative values;

[0013] S3. Extract points from the remaining point set S remain whose elevation is higher than the ground and the elevation difference from the ground belongs to the second elevation difference range to form a high point set S high ; S unclassify = S remain - S high is an unclassified point set; the second elevation difference range is [Δz low2 , Δz up2 , Δz low2 is the lower limit of the second elevation difference, and Δz up2 is the upper limit of the second elevation difference, both being positive values;

[0014] S4. Use the two-dimensional DBSCAN algorithm to segment the high point set S high to obtain a building point clustering set C build2 , and all the points included in C build2 form a set S build2 ; Remove the points in S high from S build2 to obtain a building top horizontal peripheral noise point set S horznoise = S high - S build2 , and add the points in S horznoise to the noise point set S noise ;

[0015] S5. Use the three-dimensional DBSCAN algorithm to extract building roof patch points from the set S build2 to form a set S roof ; S offroof = S build2 - S roof is a non-roof patch point set;

[0016] S6. Extract points from the non-roof patch point set S offroof whose elevation difference from the roof is within the third elevation difference range to form a high noise point set S highnoise , and add S highnoise to the noise point set S noise ; the third elevation difference range is [Δz low3 , Δz up3 , Δz low3 is the lower limit of the third elevation difference, and Δz up3 is the upper limit of the third elevation difference, both being positive values;

[0017] S7. Extract noise points that are within the building plane coverage range and have an elevation close to the ground from the unclassified point set S unclassify and add them to the ground noise point set S grdnoise , and add S grdnoise to the noise point set S noise ;

[0018] S8. Extract elevation outliers from the non-noise point set S nonoise = S all - S noise and add them to the noise point set S noise .

[0019] Furthermore, in step S1, the unique echo points in S all are segmented using the 3D DBSCAN algorithm to extract the ground point clustering set C grd ; specifically including:

[0020] S11. The unique echo points in S all form a set S 1 / 1 . Indexes are constructed for the points in S 1 / 1 according to their three-dimensional coordinates;

[0021] S12. Traverse the points in S 1 / 1 . Let the current point be P dbscan1 . Use the index to search in S 1 / 1 for whether there are at least N dbscan1 neighborhood points within a spherical neighborhood centered at P dbscan with radius R dbscan . If so, add P dbscan1 to the first core point set S core1 ; N dbscan is the minimum neighborhood number, and N dbscan is an integer within the interval [3, 8]; R dbscan = d spacing + ε 11 , where d spacing is the average ground adjacent point spacing, and ε 11 is the preset adjacent point spacing margin, and ε 11 > 0; ρ density is the planned ground point density of the laser pulse;

[0022] S13. Segment the points in the first core point set S core1 using 3D Euclidean clustering with a search radius of R dbscan , and the resulting clusters form the set C Euclidean1 ;

[0023] S14. Traverse each cluster in C Euclidean1 . Let the current cluster be C E . The points included in C E form a set S E . Traverse the points in S E . Let the current point be P E, search in S using the index for neighborhood points within the spherical neighborhood centered at P with radius R, and add them to set S; after traversing all points in S, calculate set S = S ∪ S; 1 / 1 search for neighborhood points within the spherical neighborhood centered at P with radius R in S E and add them to set S dbscan ; after traversing all points in S neighb , calculate set S E = S cluster ∪ S E ; neighb

[0024] S15. If the number of points in S is greater than or equal to N and less than or equal to N, then the points in S form cluster C, and add C to the ground point cluster set C; N is the minimum number of points for the first cluster, N = N + ε, N is the number of points on the roof with the largest area in the survey area, N = A / ρ, A is the area of the roof with the largest area in the survey area; ε > 0 is the preset margin for the minimum number of points in the cluster; N is the maximum number of points for the first cluster, and its value is the number of points in the airborne LiDAR point cloud data S; cluster If the number of points in S is greater than or equal to N min1 and less than or equal to N max1 , then the points in S form cluster C cluster , and add C ds to the ground point cluster set C ds ; N grd is the minimum number of points for the first cluster, N min1 = N min1 + ε maxroof , N 12 is the number of points on the roof with the largest area in the survey area, N maxroof = A maxroof / ρ max ; A density is the area of the roof with the largest area in the survey area; ε max > 0 is the preset margin for the minimum number of points in the cluster; N 12 is the maximum number of points for the first cluster, and its value is the number of points in the airborne LiDAR point cloud data S max1 ; all

[0025] S16. Clear set S neighb , and jump to step S14 to process the next cluster in C until all clusters in C have been processed, and C is the finally obtained ground point cluster set. Euclidean1 Euclidean1 grd

[0026] Furthermore, the step S2 of extracting points with an elevation lower than the ground and an elevation difference from the ground belonging to the first elevation difference range from the non-ground point set S includes the steps: nongrd

[0027] S21. Establish a three-dimensional index for the points in the ground point set S; grd

[0028] S22. Traverse the points in the non-ground point set S: Let the current point be P nongrd , and use the three-dimensional index to find the point P in S that is closest to P nongrd grd nongrd ​​​​​​​​​​grd ; If P nongrd and P grd 's three - dimensional Euclidean distance D1 is greater than the lower limit of the first elevation difference Δz low1 's absolute value |Δz low1 |, then re - execute S22 to process the next point in S nongrd . Otherwise, enter S23;

[0029] S23: If the elevation coordinate value z nongrd of P nongrd and the elevation coordinate value z grd of P grd satisfy:

[0030] Δz low1 ≤z nongrd - z grd ≤Δz up1 , then add P nongrd to the set S lownoise ;

[0031] S24: Jump to step S22 to process the next point in the non - ground point set S nongrd until all points in S nongrd have been processed. The points in the set S lownoise are the extracted points.

[0032] Furthermore, in step S3, extracting points with elevations higher than the ground and elevation differences from the ground within the second elevation difference range from the remaining point set S remain specifically includes:

[0033] S31: Establish a three - dimensional index for the points in the ground point set S grd ;

[0034] S32: Traverse the points in the remaining point set S remain : Let the current point be P remain . Use the three - dimensional index to find the point P grd in S remain that is closest to P grd ; If the three - dimensional Euclidean distance D2 between P remain and P grd is greater than the upper limit of the second elevation difference Δz up2 's absolute value |Δz up2 |, then re - execute S32 to process the next point in S remain . Otherwise, enter S33;

[0035] S33: If the elevation coordinate value z remain of P remain and the elevation coordinate value z grd of Pgrd The difference satisfies:

[0036] Δz low2 ≤z remain -z grd ≤Δz up2 If so, add P remain to the set S high ;

[0037] S34. Jump to step S32 to process the next point in the remaining point set S remain until all points in S remain are processed. The points in the set S high are the extracted points.

[0038] Furthermore, in step S4, the two-dimensional DBSCAN algorithm is used to segment the high-point set S high to obtain the building point clustering set C build2 , specifically including:

[0039] S41. The points corresponding to the unique echo and the first echo in the multiple echoes in the set S high form the set S2; build a two-dimensional index for the points in S2 according to the coordinates in the XY plane;

[0040] S42. Traverse the points in S2. Let the current point be P dbscan2 . Use the index to search in S2 for whether there are at least N dbscan2 neighborhood points within a circular neighborhood with P dbscan as the center and R dbscan as the radius. If so, add P dbscan2 to the second core point set S core2 ;

[0041] S43. Segment the points in the second core point set S core2 using two-dimensional Euclidean clustering, with a search radius of R dbscan . The obtained clusters form the set C Euclidean2 ;

[0042] S44. Traverse each cluster in C Euclidean2 . Let the current cluster be C E . The points included in C E form the set S E ; traverse the points in S E . Let the current point be P E . Use the two-dimensional index to search for the neighborhood points within a circular neighborhood with P E as the center and R dbscan as the radius in S2, and add them to the set S neighb ; after traversing SE After the points within it, calculate the set S cluster = S E ∪ S neighb ;

[0043] S45. If the number of points within S cluster is greater than or equal to N min2 , and less than or equal to N max2 , then the points within S cluster form a cluster C ds . Add C ds to the building point cluster set C build2 ; N min2 is the minimum number of points for the second clustering, N min2 = N minroof - ε 21 , N minroof is the number of points on the roof with the smallest area in the survey area, N minroof = A min / ρ density , A min is the area of the roof with the smallest area in the survey area; ε 21 > 0, is the preset margin of the minimum number of points for clustering; N max2 is the maximum number of points for the second clustering, and its value is the number of points in the set S high ;

[0044] S46. Clear the set S neighb , and jump to step S44 to process the next cluster in C Euclidean2 until all clusters in C Euclidean2 have been processed.

[0045] Furthermore, in the said step S7, extract the noise points within the building plane coverage range and with elevation close to the ground from the unclassified point set S unclassify , specifically including:

[0046] S71. Perform three-dimensional Euclidean clustering on the points within S unclassify , with a search radius of R dbscan , and the obtained clusters form a set C unclassify ;

[0047] S72. Traverse the clusters within C unclassify : Let the current cluster be c unclassify ;

[0048] S73. Traverse the building clusters within the building point cluster set C build2 , let the current building cluster be c build2 , and calculate the two-dimensional vector boundary b build2 of the corresponding building based on the XY coordinates of all the constituent points of c build2, calculate c unclassify Among the constituent points of buildgrd :

[0049] Condition 1: The XY coordinates of p unclassify are within the two-dimensional vector boundary b build2 ;

[0050] Condition 2: The elevation z unclassify of p unclassify is within the interval [z grdnear -Δz, z grdnear +Δz];

[0051] p unclassify ∈c unclassify , z grdnear is the elevation of the point closest to p grd among all points in S unclassify , and Δz is the elevation difference threshold;

[0052] Calculate the proportion k buildgrd of the points that meet the above two conditions: k buildgrd = N buildgrd / N unclassify , where N unclassify is the number of points in c unclassify ;

[0053] S74. If k buildgrd ≥ k th , then judge all points in c unclassify as noise points within the building plane coverage range and with elevation close to the ground, add them to the ground noise point set S grdnoise , and jump to step S72 to process the next cluster in C unclassify until all clusters in C unclassify are processed; k th is the preset proportion threshold;

[0054] If k buildgrd <k th , then jump to step S73 to judge whether the points in the next cluster in c build2 are noise points within the building plane coverage range and with elevation close to the ground, until all clusters in c unclassify are traversed, jump to step S72 to process the next cluster in C build2 until all clusters in C unclassify are processed. unclassify

[0055] Furthermore, in step S8, for the non-noise point set S nonoiseOne of the following two methods is used to extract elevation outliers:

[0056] Method 1: Calculate S nonoise The average value Z and standard deviation σ of the elevations of the inliers z ; Traverse S nonoise For the points in , determine the points with elevations not in the interval as elevation outliers;

[0057] where B low and B high are the standard deviation multiples of the lower and upper limits of the elevation interval deviating from the average value respectively;

[0058] Method 2: Sort the points in S nonoise by elevation, and determine the first N nonoise u points with the lowest elevations as elevation outliers, and / or determine the first N nonoise v points with the highest elevations as elevation outliers;

[0059] where N nonoise is the number of points in S nonoise , and u and v are the minimum elevation proportion value and maximum elevation proportion value respectively; the value ranges of u and v are [0.001%, 0.01%].

[0060] Furthermore, the above airborne LiDAR noise removal method further includes:

[0061] S9. Recover the misclassified points in the noise point set S noise , specifically:

[0062] S91. Segment the points in the noise point set S noise using three-dimensional Euclidean clustering to obtain the clustering set C noise ;

[0063] S92. Traverse the clusters in C noise , let the current cluster be c k , calculate the minimum value z k and the maximum value z min of the elevations of the points constituting c max , calculate the elevation difference Δz k of the inliers in c k =z max -z min ;

[0064] If Δz k <Δz th , then remove all the points in c k from the set S noise ; Δz th is the preset elevation difference threshold.

[0065] On the other hand, the present invention also discloses a system for implementing the method for removing airborne LiDAR noise caused by building glass, including:

[0066] A ground point and non-ground point division module, which is used to segment the airborne LiDAR point cloud data S all by using a three-dimensional DBSCAN algorithm to obtain a ground point clustering set C grd , and all the points included in C grd constitute a ground point set S grd , and the points in S all are removed from S grd to obtain a non-ground point set S nongrd ;

[0067] A sub-ground noise point extraction module, which is used to extract points from the non-ground point set S nongrd whose elevation is lower than the ground and the elevation difference from the ground belongs to the first elevation difference range to obtain a sub-ground noise point set S lownoise , and add S lownoise to the noise point set S noise ; S remain =S nongrd -S noise is the remaining point set; the first elevation difference range is [Δz low1 , Δz up1 , Δz low1 is the lower limit of the first elevation difference, and Δz up1 is the upper limit of the first elevation difference, both of which are negative values;

[0068] A high point extraction module, which is used to extract points from the remaining point set S remain whose elevation is higher than the ground and the elevation difference from the ground belongs to the second elevation difference range to form a high point set S high ; S unclassify =S remain -S high is the unclassified point set; the second elevation difference range is [Δz low2 , Δz up2 , Δz low2 is the lower limit of the second elevation difference, and Δz up2 is the upper limit of the second elevation difference, both of which are positive values;

[0069] A building top horizontal peripheral noise point extraction module, which is used to extract building top horizontal peripheral noise points, specifically: segment the high point set S high by using a two-dimensional DBSCAN algorithm to obtain a building point clustering set C build2 , and all the points included in C build2 constitute a set S build2 ; the points in S high are removed from Sbuild2 The points in build2 are used to obtain the set S of horizontal peripheral noise points on the top of the building. horznoise Let S' = S high Let S'' = S' - S build2 Then, add the points in S'' horznoise to the set S of noise points. noise ;

[0070] Building roof patch point and non-roof patch point division module, which is used to extract the building roof patch points from the set S build2 to form the set S' roof ; Let S'' = S' offroof Let S''' = S'' - S' build2 be the set of non-roof patch points; roof

[0071] Roof noise point extraction module, which is used to extract the points with elevation difference within the third elevation difference range from the set S of non-roof patch points offroof to form the set S of high-noise points highnoise Then, add S highnoise to the set S of noise points. noise The third elevation difference range is [Δz low3 , Δz up3 , where Δz low3 is the lower limit of the third elevation difference and Δz up3 is the upper limit of the third elevation difference, both of which are positive values;

[0072] Noise point extraction module inside the building and close to the ground, which is used to extract the noise points within the building plane coverage range and with elevation close to the ground from the set S of unclassified points unclassify and add them to the set S of ground noise points grdnoise Then, add S grdnoise to the set S of noise points. noise ;

[0073] Elevation outlier extraction module, which is used to extract the elevation outliers from the set S of non-noise points nonoise Let S' = S all Let S'' = S' - S noise and add them to the set S of noise points. noise ;

[0074] Furthermore, the above system further includes a misclassified noise point recovery module, which is used to recover the misclassified points in the set S of noise points, specifically: noise ;

[0075] S91. Segment the points in the set S of noise points using three-dimensional Euclidean clustering to obtain the clustering set C noise ; noise ;

[0076] S92. Traverse C​noise For the clustering in k , set the current cluster as c k , calculate the minimum value z min and the maximum value z max of the elevations of the points that make up c k , and calculate the elevation difference Δz k of the points inside c max = z min - z

[0077] If Δz k < Δz th , then remove all the points in c k from the set S noise ; Δz th is the preset elevation difference threshold.

[0078] Advantageous effects: The method and system for removing airborne LiDAR noise caused by building glass disclosed in the present invention have the following advantages:

[0079] 1. This method and system are applicable not only to the noise caused by the glass on the exterior wall of the building, but also to the noise caused by the incorrect calculation of the correspondence between the transmitted pulse and the received pulse on the high building;

[0080] 2. This method and system do not require any prior classification of the point cloud, with reasonable assumptions and being applicable in the vast majority of cases;

[0081] 3. This method and system have a stable effect, and more than 95% of the noise points can be extracted when the parameter settings are reasonable. Description of the Drawings

[0082] Figure 1 is a schematic diagram of the measurement error caused by the glass on the exterior wall of the building;

[0083] Figure 2 is a comparison schematic diagram of using single-pulse and multi-pulse technologies;

[0084] Figure 3 is a schematic diagram of the distribution of noise points and non-noise points at the exterior wall of the building;

[0085] Figure 4 is a flowchart of the method for removing airborne LiDAR noise caused by building glass disclosed in the present invention;

[0086] Figure 5 is a schematic diagram of the composition of the system for removing airborne LiDAR noise caused by building glass disclosed in the present invention. Detailed Embodiments

[0087] The present invention will be further clarified below in conjunction with the drawings and specific embodiments.

[0088] The present invention discloses a method for removing airborne LiDAR noise caused by building glass, as Figure 4 shown, which includes the following steps:

[0089] S1. Use the three-dimensional DBSCAN algorithm to segment the airborne LiDAR point cloud data S all to obtain a ground point clustering set C grd , and all the points included in C grd constitute the ground point set S grd . Remove the points in S all from S grd to obtain a non-ground point set S nongrd ;

[0090] Since the echo type of ground points is mainly the unique echo, that is, there is only one echo. To reduce the search space and improve the execution efficiency, only segment the set S all composed of points with the unique echo type in S 1 / 1 to extract the ground point clustering set C grd ; specifically including:

[0091] S11. The unique echo points in S all constitute the set S 1 / 1 . Build an index for the points in S 1 / 1 according to the three-dimensional coordinates; in this embodiment, build a KD tree or octree index based on the X, Y, and Z coordinates of the points in S 1 / 1 ;

[0092] S12. Traverse the points in S 1 / 1 . Let the current point be P dbscan1 . Use the index to search in S 1 / 1 to check whether there are at least N dbscan1 neighboring points within a spherical neighborhood with P dbscan as the center and R dbscan as the radius. If so, add P dbscan1 to the first core point set S core1 ; N dbscan is the minimum neighborhood number, which can take an integer within the interval [3, 8]. In this embodiment, the value is 4; R dbscan = d spacing + ε 11 , where d spacing is the average ground neighboring point spacing, and ε 11 is the preset neighboring point spacing margin, ε 11 > 0, that is, the value of the search radius R dbscan is a value slightly larger than the average ground neighboring point spacing, and ε 11 can be valued according to experience. ρdensity is the planned ground point density of the laser pulse;

[0093] S13. Segment the points in the first core point set S core1 using three-dimensional Euclidean clustering with a search radius of R dbscan to obtain clusters that form the set C Euclidean1 ;

[0094] S14. Traverse each cluster in C Euclidean1 Let the current cluster be C E , and the points included in C E form the set S E ; Traverse the points in S E Let the current point be P E , and use the index to search for the neighborhood points within the spherical neighborhood with P 1 / 1 as the center and R E as the radius in S dbscan , and add them to the set S neighb ; After traversing the points in S E , calculate the set S cluster = S E ∪ S neighb ;

[0095] S15. If the number of points in S cluster is greater than or equal to N min1 , and less than or equal to N max1 , then the points in S cluster form the cluster C ds , and add C ds to the ground point cluster set C grd ; N min1 is the minimum number of points in the first cluster, N min1 = N maxroof + ε 12 , N maxroof is the number of points on the roof with the largest area in the survey area, N maxroof = A max / ρ density , A max is the area of the roof with the largest area in the survey area, N min1 The purpose of such a setting is to avoid identifying the roof of a large building as a ground point; ε 12 > 0 is the preset margin of the minimum number of points in the cluster; N max1 is the maximum number of points in the first cluster, and can be simply taken as the number of points in the airborne LiDAR point cloud data S all ;

[0096] S16. Clear the set S neighb , and jump to step S14 to process C Euclidean1the next cluster in until C Euclidean1 until all clusters in C are processed grd which is the finally obtained ground point cluster set.

[0097] S2. Extract points from the non-ground point set S nongrd whose elevation is lower than the ground and the elevation difference from the ground belongs to the first elevation difference range, to obtain the underground noise point set S lownoise , and add S lownoise to the noise point set S noise ; S remain = S nongrd - S noise is the remaining point set; the first elevation difference range is [Δz low1 , Δz up1 , where Δz low1 is the lower limit of the first elevation difference and Δz up1 is the upper limit of the first elevation difference, both of which are negative values;

[0098] Extracting points from the non-ground point set S nongrd whose elevation is lower than the ground and the elevation difference from the ground belongs to the first elevation difference range includes the steps of:

[0099] S21. Establish a three-dimensional index for the points in the ground point set S grd ;

[0100] Similarly, in this embodiment, a KD tree or octree index is established according to the X, Y, and Z coordinates of the points in S grd ;

[0101] S22. Traverse the points in the non-ground point set S nongrd : Let the current point be P nongrd , and use the three-dimensional index to find the point P grd in S nongrd that is closest to P grd ; if the three-dimensional Euclidean distance D1 between P nongrd and P grd is greater than the absolute value Δz low1 of the lower limit Δz low1 of the first elevation difference, then re-execute S22 to process the next point in S nongrd , otherwise enter S23; the limitation of the three-dimensional Euclidean distance between P nongrd and P grd here is to limit the search range and avoid using ground points that are too far away to judge the elevation difference;

[0102] S23. If the elevation coordinate value z nongrd of P nongrd and the elevation coordinate value z grd of Pgrd The difference satisfies:

[0103] Δz low1 ≤z nongrd -z grd ≤Δz up1 If so, add P nongrd to the set S lownoise ;

[0104] S24. Jump to step S22 to process the next point in the non-ground point set S nongrd until all points in S nongrd are processed. The points in the set S lownoise are the extracted points.

[0105] In step S2, noise points that are lower than the ground and whose elevation difference from the ground is within [Δz low1 , Δz up1 are extracted, that is, Figure 1 the P_points in. The first elevation difference lower limit Δz low1 can take a relatively large negative value, and the specific value is set according to requirements. For example, if the lowest noise point is 95 m lower than the ground, then Δz low1 can take 100 m. Δz up1 is negative; to avoid classifying river embankments, stairs leading to the underground, etc. as noise, Δz up1 can be within the range of [-1.0, -5.0].

[0106] S3. Extract points from the remaining point set S remain whose elevation is higher than the ground and whose elevation difference from the ground belongs to the second elevation difference range to form a high point set S high ; S unclassify = S remain - S high is the unclassified point set; the second elevation difference range is [Δz low2 , Δz up2 , Δz low2 is the second elevation difference lower limit, and Δz up2 is the second elevation difference upper limit, both of which are positive values;

[0107] Extract points from the remaining point set S remain whose elevation is higher than the ground and whose elevation difference from the ground belongs to the second elevation difference range, specifically including:

[0108] S31. Establish a three-dimensional index for the points in the ground point set S grd ;

[0109] S32. Traverse the points in the remaining point set S remain : Let the current point be P remain, find S using three-dimensional indexing grd The mid-distance P remain The nearest point P grd ; If P remain and P grd 's three-dimensional Euclidean distance D2 is greater than the absolute value |Δz up2 | of the upper limit of the second elevation difference, then re-execute S32 to process the next point in S up2 |, otherwise enter S33; remain

[0110] S33, If the difference between the elevation coordinate value z remain of P remain and the elevation coordinate value z grd of P grd satisfies:

[0111] Δz low2 ≤z remain -z grd ≤Δz up2 , then add P remain to the set S high ;

[0112] S34, Jump to step S32 to process the next point in the remaining point set S remain until all points in S remain are processed. The points in the set S high are the extracted points.

[0113] Since high-noise points are generally at a high altitude, beyond the elevation range of common trees, so extract points that are high enough above the ground (including points on buildings and high-noise points above the ground) to find candidate noise points from them. To ensure that higher building points and the noise on them are extracted, Δz low2 can empirically take a positive value, such as 50m, and Δz up2 can take a larger positive value, higher than the building height, and in this embodiment, the value is 200m.

[0114] S4, Use the two-dimensional DBSCAN algorithm to segment the high-point set S high to obtain the building point clustering set C build2 , and all the points included in C build2 constitute the set S build2 ; Remove the points in S high from S build2 to obtain the building top horizontal perimeter noise point set S horznoise =S high -S build2 , and add the points in S horznoise to the noise point set S noise ;​

[0115] The main purpose of step S4 is to extract non-noise points located on the building and noise points on the horizontal perimeter of the building in the XY plane. Since the building roof is generally impenetrable and its echo types are mainly single echo and the first echo of multiple echoes, in order to reduce the search space, step S4 only selects points with single echo or the first echo of multiple echoes from S high for segmentation, so as to obtain the building point clustering set C build2 , specifically including:

[0116] S41. The points corresponding to the single echo and the first echo of multiple echoes in set S high constitute set S2; a two-dimensional index is constructed for the points in S2 according to the XY coordinates, and a KD tree or a quadtree can be used.

[0117] S42. Traverse the points in S2, and let the current point be P dbscan2 . Use the two-dimensional index to search in S2 whether there are at least N dbscan2 neighboring points within the circular neighborhood with P dbscan as the center and R dbscan as the radius. If so, add P dbscan2 to the second core point set S core2 ;

[0118] S43. Segment the points in the second core point set S core2 using two-dimensional Euclidean clustering, with a search radius of R dbscan , and the obtained clusters constitute set C Euclidean2 ;

[0119] S44. Traverse each cluster in C Euclidean2 , and let the current cluster be C E . The points included in C E constitute set S E ; Traverse the points in S E , and let the current point be P E . Use the two-dimensional index to search for neighboring points within the circular neighborhood with P E as the center and R dbscan as the radius in S2, and add them to set S neighb ; After traversing the points in S E , calculate set S cluster = S E ∪S neighb ;

[0120] S45. If the number of points in S cluster is greater than or equal to N min2 , and less than or equal to N max2 , then Scluster The points within form cluster C ds , add C ds to the building point cluster set C build2 ; N min2 is the minimum number of points for the second cluster, N min2 = N minroof - ε 21 , N minroof is the number of points on the roof with the smallest area in the survey area, N minroof = A min / ρ density , A min is the area of the roof with the smallest area in the survey area; ε 21 > 0, is the preset margin of the minimum number of points for clustering; N max2 is the maximum number of points for the second cluster, and the value is the number of points in set S high ;

[0121] S46. Clear set S neighb , and jump to step S44 to process the next cluster in C Euclidean2 until all clusters in C Euclidean2 have been processed.

[0122] The set S of horizontal peripheral noise points on the top of the building obtained in step S4 horznoise , that is Figure 1 part of the points of P + in.

[0123] S5. Extract the building roof patch points from set S build2 to form set S roof ; S offroof = S build2 - S roof is the non-roof patch point set;

[0124] Step S5 uses the 3D DBSCAN algorithm to extract the building roof patch points from set S build2 to form set S roof , and the specific steps are similar to those of extracting the ground point set S 1 / 1 from set S grd in S1, including:

[0125] S51. Build an index for the points in S build2 according to the XYZ coordinates;

[0126] S52. Traverse the points in S build2 , let the current point be P dbscan3 , and use the index to search in S build2 for the points with P dbscan3 as the center and R dbscanWithin a spherical neighborhood with a radius, does there exist at least N dbscan neighborhood points? If so, add P dbscan3 to the third core point set S core3 ;

[0127] S53. Segment the points in the third core point set S core3 using three-dimensional Euclidean clustering with a search radius of R dbscan , and the resulting clusters form the set C Euclidean3 ;

[0128] S54. Traverse each cluster in C Euclidean3 . Let the current cluster be C E , and the points included in C E form the set S E ; Traverse the points in S E . Let the current point be P E , and use the index to search for neighborhood points within a spherical neighborhood centered at P build2 with a radius of R E in S dbscan , and add them to the set S neighb ; After traversing the points in S E , calculate the set S cluster = S E ∪ S neighb ;

[0129] S55. If the number of points in S cluster is greater than or equal to N min3 , and less than or equal to N max3 , then the points in S cluster form the cluster C ds , and add C ds to the roof patch point cluster set C roof ; N min3 is the minimum number of points in the third cluster, N min3 = N minroofseg - ε 31 , N minroofseg is the number of points on the roof with the smallest area in the survey area, N minroofseg = A minseg / ρ density , A minseg is the area of the roof with the smallest area in the survey area; ε 31 > 0 is the preset margin of the minimum number of points in the third cluster; N max3 is the maximum number of points in the third cluster, and its value is the number of points in S build2 ;

[0130] S56. Clear the set S neighb , and jump to step S54 to process C Euclidean3the next cluster in until C Euclidean3 until all clusters in C are processed roof is the finally obtained set of roof patch point clusters, C roof the points in form the set S roof .

[0131] S6. Extract from the non-roof patch point set S offroof the points with elevation difference from the roof within the third elevation difference range to form the high-noise point set S highnoise , and add S highnoise to the noise point set S noise ; the third elevation difference range is [Δz low3 , Δz up3 , where Δz low3 is the lower limit of the third elevation difference, and Δz up3 is the upper limit of the third elevation difference, both being positive values; the specific steps for extracting S highnoise include:

[0132] S61. Establish a three-dimensional index for the points in the building roof patch point set S roof .

[0133] S62. Traverse the points in the non-roof patch point set S offroof : Let the current point be P offroof , and use the three-dimensional index to find the point P roof in S offroof that is closest to P roof ; if the three-dimensional Euclidean distance D3 between P offroof and P roof is greater than the absolute value |Δz up3 | of the upper limit Δz up3 of the third elevation difference, then re-execute S62 to process the next point in S offroof , otherwise enter S63; the value of the upper limit Δz up3 of the third elevation difference should be less than or equal to the minimum distance between adjacent high-rise buildings, so as to only use the building roofs near directly below the noise points.

[0134] S63. If the difference between the elevation coordinate value z offroof of P offroof and the elevation coordinate value z roof of P roof satisfies: Δz low3 ≤ z offroof - z roof ≤ Δz up3 , then add P offroof to the high-noise point set S highnoise ;

[0135] The lower limit Δz of the third elevation differencelow3 Set as the elevation difference between the rooftop accessory and the roof. The rooftop accessory is an accessory structure such as a parapet wall, a drying rack, a light well, etc. set on the roof. In this embodiment, Δz low3 is set to 2m. Δz up3 Can take a relatively large positive value to include the noise point with the largest elevation difference from the roof.

[0136] S64. Jump to step S62 to process the next point in the non-roof patch point set S offroof until all points in S offroof are processed. The points in the set S highnoise are the extracted points. After that, add S highnoise to the noise point set S noise .

[0137] The noise points in the set S highnoise extracted in step S6 are Figure 1 part of the points in P + .

[0138] S7. Extract the noise points within the building plane coverage range and close to the ground elevation from the unclassified point set S unclassify , and add them to the ground noise point set S grdnoise . Add S grdnoise to the noise point set S noise ; specifically including:

[0139] S71. Perform three-dimensional Euclidean clustering on the points in S unclassify with a search radius of R dbscan to obtain the clusters that form the set C unclassify ; The minimum number of points for three-dimensional Euclidean clustering can be set to 1, and the maximum number of points can empirically take a value greater than or equal to 10 and less than or equal to 100;

[0140] S72. Traverse the clusters in C unclassify : Let the current cluster be c unclassify ;

[0141] S73. Traverse the building clusters in the building point cluster set C build2 . Let the current building cluster be c build2 . Based on the XY coordinates of all the constituent points of c build2 , calculate the two-dimensional vector boundary b build2 of its corresponding building. Calculate the number N unclassify of the points that simultaneously satisfy the following two conditions among the constituent points of c buildgrd :

[0142] Condition 1: p unclassifyThe XY coordinates of are located within the two-dimensional vector boundary b build2 ;

[0143] Condition 2: p unclassify has an elevation z unclassify within the interval [z grdnear -Δz, z grdnear +Δz];

[0144] p unclassify ∈c unclassify , z grdnear is the elevation of the point closest to p grd among all points in S unclassify . Δz is the elevation difference threshold, and a value within [2, 5] can be empirically taken; in specific implementation, methods such as the raster-vector conversion method, the Alpha-shape method, and the boundary tracing algorithm based on Delaunay triangulation can be used to calculate the two-dimensional vector boundary.

[0145] Calculate the proportion k buildgrd of the points that meet the above two conditions: k buildgrd = N buildgrd / N unclassify , where N unclassify is the number of points in c unclassify ;

[0146] S74. If k buildgrd ≥ k th , then judge all points in c unclassify as noise points within the building plane coverage range and with an elevation close to the ground, and jump to step S72 to process the next cluster in C unclassify until all clusters in C unclassify are processed; k th is a preset proportion threshold, with a value greater than 0.5 and less than or equal to 1;

[0147] If k buildgrd < k th , then jump to step S73 to judge whether the points in c build2 are noise points within the building plane coverage range and with an elevation close to the ground for the next cluster in c unclassify until all clusters in c build2 are traversed, then jump to step S72 to process the next cluster in C unclassify until all clusters in C unclassify are processed, and add the ground noise point set S grdnoise to S noise .

[0148] The noise points extracted in step S7 that are within the building plane coverage range and have an elevation close to the ground areFigure 1 Point P0

[0149] S8. Extract elevation outliers from the non-noise point set S nonoise = S all - S noise and add them to the noise point set S noise .

[0150] Elevation outliers refer to points with extremely high or low elevations, which may be caused by the glass of building facades, or by flying birds or instrument noise, etc. For the non-noise point set S nonoise extracting elevation outliers can be done using a method based on the normal distribution or a fixed ratio, specifically one of the following two methods:

[0151] Method 1: Calculate the mean Z and standard deviation σ of the elevations of the inlier points in S nonoise ; Traverse the points in S z ; Determine the points with elevations not in the interval nonoise as elevation outliers;

[0152] where B low and B high are the multiples of the standard deviation by which the lower and upper limits of the elevation interval deviate from the mean, respectively, and can take values greater than or equal to 3.0 according to the elevation distribution of the point cloud.

[0153] Method 2: Sort the points in S nonoise by elevation, and determine the first N nonoise u points with the lowest elevations as elevation outliers, and / or determine the first N nonoise v points with the highest elevations as elevation outliers;

[0154] where N nonoise is the number of points in S nonoise , and u and v are the minimum elevation ratio value and the maximum elevation ratio value, respectively; the value ranges of u and v are [0.001%, 0.01%].

[0155] Through the above steps S1 - S8, the noise points in the airborne LiDAR point cloud data S all are extracted to obtain the noise point set S noise , but there may be misclassified points in it, that is, points that are not noise are classified as noise points, and the misclassified points can be restored through step S9.

[0156] S9. Restore the misclassified points in the noise point set S noise , specifically as follows:

[0157] S91. Restore the misclassified points in the noise point set S noise ​Points within are segmented using three-dimensional Euclidean clustering to obtain a clustering set C noise ; The maximum number of points for three-dimensional Euclidean clustering can be set to S noise The number of points, the minimum number of points can empirically be set to an integer greater than or equal to 10, and the search radius can be set to R dbscan ;

[0158] S92. Traverse the clusters in C noise Let the current cluster be c k , and calculate the minimum value z k of the elevations of the points that make up c min and the maximum value z max , and calculate the elevation difference Δz k of the points within c k = z max - z min ;

[0159] If Δz k < Δz th , then remove all the points in c k from the set S noise ; Δz th is a preset elevation difference threshold. In this embodiment, the value of Δz th is the absolute value of the first elevation difference upper limit Δz up1 in step S23, that is, Δz th = |Δz up1 |.

[0160] The present invention also discloses a system for implementing the above method for removing airborne LiDAR noise caused by building glass, as Figure 5 shown, including:

[0161] A ground point and non-ground point division module 1, configured to segment the airborne LiDAR point cloud data S all using a three-dimensional DBSCAN algorithm to obtain a ground point clustering set C grd , and all the points included in C grd constitute a ground point set S grd , and the points in S all are removed from S grd to obtain a non-ground point set S nongrd ;

[0162] A below-ground noise point extraction module 2, configured to extract points from the non-ground point set S nongrd whose elevations are lower than the ground and whose elevation differences from the ground belong to a first elevation difference range to obtain a below-ground noise point set S lownoise , and add S lownoise to the noise point set S noise ; S remain = Snongrd -S noise is the remaining point set; the first elevation difference range is [Δz low1 , Δz up1 , where Δz low1 is the lower limit of the first elevation difference, and Δz up1 is the upper limit of the first elevation difference, both of which are negative values;

[0163] The high point extraction module 3 is used to extract points from the remaining point set S remain whose elevation is higher than the ground and the elevation difference from the ground belongs to the second elevation difference range, and form a high point set S high ; S unclassify = S remain - S high is the unclassified point set; the second elevation difference range is [Δz low2 , Δz up2 , where Δz low2 is the lower limit of the second elevation difference, and Δz up2 is the upper limit of the second elevation difference, both of which are positive values;

[0164] The building top horizontal peripheral noise point extraction module 4 is used to extract the building top horizontal peripheral noise points. Specifically: use the two-dimensional DBSCAN algorithm to segment the high point set S high to obtain the building point clustering set C build2 , and all the points included in C build2 form the set S build2 ; Remove the points in S high from S build2 to obtain the building top horizontal peripheral noise point set S horznoise = S high - S build2 , and add the points in S horznoise to the noise point set S noise ;

[0165] The building roof patch point and non-roof patch point division module 5 is used to extract the building roof patch points from the set S build2 to form the set S roof ; S offroof = S build2 - S roof is the non-roof patch point set;

[0166] The roof noise point extraction module 6 is used to extract points with an elevation difference from the roof within the third elevation difference range from the non-roof patch point set S offroof to form the high noise point set S highnoise , and add S highnoise to the noise point set S noise ; The third elevation difference range is [Δzlow3 , Δz up3 , Δz low3 is the lower limit of the third elevation difference, Δz up3 is the upper limit of the third elevation difference, both being positive values;

[0167] The noise point extraction module 7 inside the building and close to the ground is used to extract the noise points within the building plane coverage range and with an elevation close to the ground from the unclassified point set S unclassify and add them to the ground noise point set S grdnoise , and add S grdnoise to the noise point set S noise ;

[0168] The elevation outlier extraction module 8 is used to extract the elevation outliers from the non-noise point set S nonoise = S all - S noise and add them to the noise point set S noise .

[0169] The misclassified noise point recovery module 9 is used to recover the misclassified points in the noise point set S noise according to steps S91 - S92.

Claims

1. A method for removing airborne LiDAR noise caused by building glass, characterized in that: including the following steps: S1. Use the 3D DBSCAN algorithm to segment the airborne LiDAR point cloud data S all to obtain the ground point clustering set C grd , where all the points included in C grd constitute the ground point set S grd . Remove the points in S all from S grd to obtain the non-ground point set S nongrd ; S2. Extract points from the non-ground point set S nongrd with elevations lower than the ground and elevation differences from the ground within the first elevation difference range to obtain the underground noise point set S lownoise , and add S lownoise to the noise point set S noise ; S remain = S nongrd - S noise is the remaining point set; the first elevation difference range is [Δz low1 , Δz up1 , where Δz low1 is the lower limit of the first elevation difference and Δz up1 is the upper limit of the first elevation difference, both being negative values; S3. Extract points from the remaining point set S remain whose elevation is higher than the ground and the elevation difference from the ground elevation belongs to the second elevation difference range to form a high point set S high ; S unclassify = S remain - S high is the unclassified point set; the second elevation difference range is [Δz low2 , Δz up2 , Δz low2 is the lower limit of the second elevation difference, and Δz up2 is the upper limit of the second elevation difference, both of which are positive values; S4. Use the two-dimensional DBSCAN algorithm to segment the high-point set S high to obtain the building point clustering set C build2 , and all the points included in C build2 constitute the set S build2 ; Remove the points in S high from S build2 to obtain the set S horznoise of horizontal peripheral noise points at the top of the building = S high - S build2 , and add the points in S horznoise to the noise point set S noise ; S5. Use the 3D DBSCAN algorithm to extract the building roof patch points from the set S build2 to form the set S roof ; S offroof = S build2 - S roof is the non-roof patch point set; S6. Extract points from the non-roof patch point set S offroof whose elevation difference from the roof is within the third elevation difference range to form a high-noise point set S highnoise , and add S highnoise to the noise point set S noise ; the third elevation difference range is [Δz low3 , Δz up3 , where Δz low3 is the lower limit of the third elevation difference and Δz up3 is the upper limit of the third elevation difference, both being positive values; S7. Extract noise points from the unclassified point set S unclassify that are within the building's planar coverage and have an elevation close to the ground, and add them to the ground noise point set S grdnoise . Add S grdnoise to the noise point set S noise ; S8. Extract elevation outliers from the non-noise point set S nonoise = S all - S noise and add them to the noise point set S noise .

2. The method for removing airborne LiDAR noise caused by building glass according to claim 1, characterized in that: In the step S1, the only echo point in S all is segmented by using the three-dimensional DBSCAN algorithm to extract the ground point clustering set C grd ; specifically including: S11, S all The only echo points in it form a set S 1 / 1 , for S 1 / 1 Build an index for the points in it according to the three-dimensional coordinates; S12. Traverse the points in S 1 / 1 Let the current point be P dbscan1 , and use the index to search in S 1 / 1 for whether there are at least N dbscan1 neighborhood points within the spherical neighborhood centered at P dbscan with radius R dbscan . If so, add P dbscan1 to the first core point set S core1 ; N dbscan is the minimum number of neighborhoods, and N dbscan is an integer within the interval [3, 8]; R dbscan = d spacing + ε 11 , where d spacing is the average ground adjacent point spacing, and ε 11 is the preset adjacent point spacing margin, and ε 11 > 0; ρ density is the planned ground point density of the laser pulse; S13. Segment the points in the first core point set S core1 using three-dimensional Euclidean clustering with a search radius of R dbscan , and the resulting clusters form the set C Euclidean1 ; S14. Traverse C Euclidean1 For each cluster within C, let the current cluster be C E , where C E contains points that form a set S E ; Traverse the points within S E , let the current point be P E , and use the index to search for neighborhood points within the spherical neighborhood centered at P 1 / 1 with radius R E in S dbscan , and add them to the set S neighb ; After traversing the points within S E , calculate the set S cluster = S E ∪ S neighb ; S15. If the number of points in S cluster is greater than or equal to N min1 and less than or equal to N max1 , then the points in S cluster form a cluster C ds . Add C ds to the ground point cluster set C grd . N min1 is the minimum number of points for the first cluster, N min1 = N maxroof + ε 12 . N maxroof is the number of points of the roof with the largest area in the survey area. N maxroof = A max / ρ density . A max is the area of the roof with the largest area in the survey area. ε 12 > 0 is the preset margin of the minimum number of points for clustering. N max1 is the maximum number of points for the first cluster, and its value is the number of points in the airborne LiDAR point cloud data S all . S16. Clear the set S neighb , jump to step S14 to process the next cluster in C Euclidean1 until all clusters in C Euclidean1 have been processed. C grd is the finally obtained ground point cluster set.

3. The method for removing airborne LiDAR noise caused by building glass according to claim 1, characterized in that The step S2 extracts points from the non-ground point set S nongrd with elevations lower than the ground and elevation differences from the ground within a first elevation difference range, including the steps of: S21. Build a three-dimensional index for the points in the ground point set S grd ; S22, for non-ground point set S nongrd Traverse the points inside: Set the current point to be P nongrd , use the three-dimensional index to find S grd Middle distance P nongrd The nearest point P grd ; If P nongrd With P grd The three-dimensional Euclidean distance D1 is greater than the first elevation difference lower limit Δz low1 The absolute value of |Δz low1 |, then re-execute S22 and process S nongrd The next point in, otherwise go to S23; S23. If P nongrd 's elevation coordinate value z nongrd and P grd 's elevation coordinate value z grd 's difference satisfies: Δz low1 ≤z nongrd -z grd ≤Δz up1 If so, then add P nongrd to set S lownoise ; S24. Jump to step S22 to process the next point in the non-ground point set S nongrd until all points in S nongrd have been processed. The points in the set S lownoise are the extracted points.

4. The method for removing airborne LiDAR noise caused by building glass according to claim 1, wherein In step S3, points with an elevation higher than the ground and a difference in elevation from the ground within the second elevation difference range are extracted from the remaining point set S remain Specifically, it includes: S31, ground point set S grd Establish a three-dimensional index for the points in; S32, for the remaining point set S remain Traverse the points inside: Set the current point to be P remain , use the three-dimensional index to find S grd Middle distance P remain The nearest point P grd ; If P remain With P grd The three-dimensional Euclidean distance D2 is greater than the second elevation difference upper limit Δz up2 The absolute value of |Δz up2 |, then re-execute S32 and process S remain The next point in, otherwise go to S33; S33. If P remain 's elevation coordinate value z remain and P grd 's elevation coordinate value z grd 's difference satisfies: Δz low2 ≤z remain -z grd ≤Δz up2 , then P remain Add to collection S high ; S34. Jump to step S32 to process the next point in the remaining point set S remain until all points in S remain have been processed. The points in the set S high are the extracted points.

5. The method for removing airborne LiDAR noise caused by building glass according to claim 1, wherein, In the step S4, the two-dimensional DBSCAN algorithm is used to segment the high point set S high to obtain the building point clustering set C build2 , which specifically includes: S41. Set S high The points corresponding to the only echo and the first echo among the multiple echoes in the set form set S2; two-dimensional indexes are constructed for the points in S2 according to the coordinates in the XY plane; S42. Traverse the points in S2, and let the current point be P dbscan2 , and use the index to search in S2 for a circular neighborhood with P dbscan2 as the center and R dbscan as the radius to check if there are at least N dbscan neighborhood points. If so, add P dbscan2 to the second core point set S core2 ; S43. Segment the points in the second core point set S core2 using two-dimensional Euclidean clustering with a search radius of R dbscan , and the resulting clusters form the set C Euclidean2 ; S44. Traverse C Euclidean2 For each cluster within C, let the current cluster be C E , C E The points it contains form the set S E ; Traverse the points within S E , let the current point be P E , use the two-dimensional index to search in S2 for the neighborhood points within the circular neighborhood centered at P E with radius R dbscan , and add them to the set S neighb ; After traversing the points within S E , calculate the set S cluster = S E ∪ S neighb ; S45. If the number of points in S cluster is greater than or equal to N min2 and less than or equal to N max2 , then the points in S cluster form a cluster C ds . Add C ds to the building point cluster set C build2 ; N min2 is the minimum number of points for the second cluster, N min2 = N minroof - ε 21 , N minroof is the number of points on the roof with the smallest area in the survey area, N minroof = A min / ρ density , A min is the area of the roof with the smallest area in the survey area; ε 21 > 0 is the preset margin of the minimum number of points for clustering; N max2 is the maximum number of points for the second cluster, and its value is the number of points in the set S high ; S46. Clear the set S neighb , and jump to step S44 to process the next cluster in C Euclidean2 until all clusters in C Euclidean2 have been processed.

6. The method for removing airborne LiDAR noise caused by building glass according to claim 1, characterized in that, In the step S7, noise points that are within the building plane coverage range and have an elevation close to the ground are extracted from the unclassified point set S unclassify , specifically including: S71. Perform three-dimensional Euclidean clustering on the points in S unclassify with a search radius of R dbscan to obtain clusters that form a set C unclassify ; S72. Traverse the clusters within C unclassify Let the current cluster be c unclassify ; S73. Cluster the building points in set C build2 Traverse the building clusters in C, and let the current building cluster be c build2 , and based on c build2 Calculate the two-dimensional vector boundary b of the corresponding building based on the XY coordinates of all the constituent points build2 , calculate the number N of points that satisfy the following two conditions simultaneously among the constituent points of c unclassify : buildgrd ​ Condition 1: p unclassify The XY coordinates of build2 are within the two-dimensional vector boundary b Condition 2: p unclassify with elevation z unclassify in the interval [z grdnear -Δz, z grdnear +Δz]; p unclassify ∈c unclassify ,z grdnear is S grd Among all the points in, the elevation of the point closest to p unclassify is the closest, and Δz is the elevation difference threshold; Calculate the proportion k of the points that meet the above two conditions buildgrd : k buildgrd = N buildgrd / N unclassify , where N unclassify is the number of midpoints of c unclassify ; S74. If k buildgrd ≥ k th , then all the points within c unclassify are determined to be noise points within the building plane coverage range and with elevations close to the ground, and are added to the ground noise point set S grdnoise , and jump to step S72 to process the next cluster in C unclassify until all the clusters in C unclassify have been processed; k th is a preset ratio threshold; If k buildgrd <k th , then jump to step S73, for c build2 The next cluster judgment c in unclassify Is the point in the noise point within the building plane coverage and close to the ground elevation? build2 All clusters in are traversed, jump to step S72 to process C unclassify The next cluster in C unclassify All clusters in are processed.

7. The method for removing airborne LiDAR noise caused by building glass according to claim 1, wherein, The step S8 is for the non-noise point set S nonoise One of the following two methods is adopted to extract elevation outliers: Method 1: Calculate S nonoise The average elevation of the interior points and the standard deviation σ z ; Traverse the points in S nonoise and determine the points with elevations not within the interval as elevation outliers; where B low and B high are the standard deviation multiples of the lower and upper limits of the elevation interval deviating from the average value, respectively; Method 2: Sort the points in S nonoise by elevation, and determine the first N nonoise u points with the lowest elevation as elevation outliers, and / or determine the first N nonoise v points with the highest elevation as elevation outliers; Where N nonoise For S nonoise The number of points in , u and v are the minimum elevation ratio value and the maximum elevation ratio value respectively; the value range of u, v is [0.001%, 0.01%].

8. The method for removing airborne LiDAR noise caused by building glass according to claim 1, wherein further including: S9. Recover the misclassified points in the noise point set S noise Specifically, S91, for the noise point set S noise The points in the cluster are segmented using three-dimensional Euclidean clustering, and the resulting cluster set C noise ; S92. Traverse C noise For the clusters in k , set the current cluster as c k , calculate the minimum value z min and the maximum value z max of the elevations of the points that make up c k , and calculate the elevation difference Δz k of the points inside c max = z min - z If Δz k <Δz th , then all points in c k are removed from the set S noise ; Δz th is a preset elevation difference threshold.

9. An airborne LiDAR noise removal system caused by building glass, characterized in that, including: The ground point and non-ground point segmentation module is used to use the 3D DBSCAN algorithm to analyze the airborne LiDAR point cloud data S all Perform segmentation to obtain the ground point cluster set C grd , C grd All points included constitute the ground point set S grd , S all Remove the S grd Points in the, get the non-ground point set S nongrd ; Below-ground noise point extraction module, used to extract points from the non-ground point set S nongrd whose elevation is lower than the ground and the elevation difference from the ground belongs to the first elevation difference range, to obtain the underground noise point set S lownoise , and add S lownoise to the noise point set S noise ; S remain = S nongrd - S noise is the remaining point set; the first elevation difference range is [Δz low1 , Δz up1 , Δz low1 is the lower limit of the first elevation difference, and Δz up1 is the upper limit of the first elevation difference, both of which are negative values; A high point extraction module for extracting points from the remaining point set S remain whose elevation is higher than the ground and the elevation difference from the ground elevation belongs to the second elevation difference range, to form a high point set S high ; S unclassify = S remain - S high is the unclassified point set; the second elevation difference range is [Δz low2 , Δz up2 , Δz low2 is the lower limit of the second elevation difference, and Δz up2 is the upper limit of the second elevation difference, both of which are positive values; Building top horizontal perimeter noise point extraction module, which is used to extract the building top horizontal perimeter noise points. Specifically: Use the two-dimensional DBSCAN algorithm to segment the high point set S high to obtain the building point clustering set C build2 , and all the points included in C build2 constitute the set S build2 ; Remove the points in S high from S build2 to obtain the building top horizontal perimeter noise point set S horznoise = S high - S build2 , and add the points in S horznoise to the noise point set S noise ; A building roof patch point and non-roof patch point division module, which is used to extract building roof patch points from the set S build2 to form the set S roof ; S offroof = S build2 - S roof is the non-roof patch point set; Roof noise point extraction module is used to extract the noise from the non-roof patch point set S offroof Extract the points whose elevation difference with the roof is within the third elevation difference range from the original point to form a high noise point set S highnoise , S highnoise Add noise point set S noise The third elevation difference range is [Δz low3 ,Δz up3 ], Δz low3 is the lower limit of the third elevation difference, Δz up3 is the upper limit of the third elevation difference, all of which are positive; A noise point extraction module inside the building and close to the ground, which is used to extract noise points within the building plane coverage range and with an elevation close to the ground from the unclassified point set S unclassify and add them to the ground noise point set S grdnoise , add S grdnoise to the noise point set S noise ; The elevation outlier extraction module is used to extract the outlier points from the non-noise point set S nonoise =S all -S noise Extract the elevation outliers and add the noise point set S noise middle.

10. The airborne LiDAR noise removal system caused by building glass according to claim 9, characterized in that, It further includes a misclassified noise point recovery module for recovering misclassified points in the noise point set S noise specifically as follows: S91. Segment the points in the noise point set S noise using three-dimensional Euclidean clustering to obtain the clustering set C noise ; S92. Traverse C noise For the clusters in k , assume the current cluster is c k . Calculate the minimum value z min and the maximum value z max of the elevations of the points that make up c k . Calculate the elevation difference Δz k of the points inside c max = z min - z If Δz k <Δz th , then all points in c k are removed from the set S noise ; Δz th is a preset elevation difference threshold.

Images