A method for filtering airborne point cloud based on fracture line constraint
By employing an airborne point cloud filtering method based on fracture line constraints, combined with fracture point extraction based on curvature and elevation differences, coarse and fine filtering are performed to solve the problem of accurate extraction of ground points under complex terrain, thus forming a high-precision digital elevation model.
Patent Information
- Authority / Receiving Office
- CN · China
- Patent Type
- Patents(China)
- Current Assignee / Owner
- NANJING UNIV OF INFORMATION SCI & TECH
- Filing Date
- 2022-08-31
- Publication Date
- 2026-04-28
AI Technical Summary
Existing airborne laser point cloud filtering methods struggle to effectively distinguish between ground and non-ground points in complex terrains and scenarios, often leading to the erroneous filtering of ground points. This is especially true under conditions of abrupt elevation changes and noise, where the filtering effect is poor.
A fast airborne point cloud filtering method based on fracture line constraints is adopted. By extracting fracture points and combining curvature values and elevation differences, coarse and fine filtering are performed. Slope and flatness are calculated using equilateral triangles, and non-ground points are filtered out by combining quadratic surface fitting method.
It improves the completeness and accuracy of ground point extraction, is applicable to a variety of complex scenarios, overcomes the effects of elevation changes and noise, and forms a high-precision digital elevation model.
Smart Images
Figure CN115294302B_ABST
Abstract
Description
Technical Field
[0001] This invention relates to the technical field of intelligent processing of laser point cloud measurement data and application research of geographic information data, specifically to a fast airborne point cloud filtering method based on fracture line constraints. Background Technology
[0002] Airborne laser point cloud filtering, as a means of accurately separating and identifying ground points, is a crucial step in building high-precision digital terrain models. Accurate and high-precision point cloud filtering plays a vital role in applications such as ground information extraction, road planning, and improving target recognition accuracy. Laser scanning systems can obtain azimuth information through scanning, providing a fast and accurate solution for large-area mapping. In the point cloud data processing process, ground information is extremely important; however, the sheer number of ground points and their integrity are easily affected by noise, occlusion, and other factors, leading to slow data processing speeds and difficulties in target classification and identification.
[0003] Currently, there are three main types of airborne laser point cloud filtering methods: The first type is based on slope variation, calculating the slope between each point and its neighbors, and filtering out features by comparing the actual slope between two points with a slope threshold. This type of method tends to incorrectly filter out a large number of ground points when dealing with steep slopes or areas with abrupt elevation changes, and the steeper the slope, the lower the filtering accuracy. The second type is mathematical morphology methods, which perform erosion and dilation on the terrain point cloud to construct an opening operation. Utilizing the large elevation changes of non-ground points after the opening operation, points with elevation differences greater than a threshold are classified as non-ground points and filtered out. The filtering results obtained by this method depend on the choice of the optimal window size; a smaller window size will preserve larger buildings, while a larger window will smooth out terrain details. The third type is surface-based methods, which iteratively select ground measurements from the original dataset to gradually approximate the ground surface, creating an approximate bare Earth surface. This type of method works well in areas with gentle terrain changes, but it often misclassifies the top surfaces of large buildings as terrain surfaces.
[0004] In recent years, the theory of LiDAR scan lines has been widely applied to airborne LiDAR point cloud filtering. Scan line-based filtering models can consider the terrain features of the entire scanned area, enabling fast and efficient processing of large areas and large data volumes, providing a new approach for filtering ground features from airborne LiDAR point clouds. Yin et al. used spherical coordinates instead of Cartesian coordinates. Their algorithm does not assume roads are at the lowest point, using azimuth and radial distance spherical coordinates as criteria. Ground points and non-ground points are segmented by the breakpoints and inflection points of the radial distance curve in the spherical coordinates. This algorithm avoids over-segmentation and under-segmentation. Zheng et al. processed point clouds using variable radius annulus and B-spline fitting to obtain ground point clouds. For each scan line, a variable radius annulus was used to roll over the points on the scan line, retaining those points. B-spline fitting was then used to fit the surface, traversing each point in the region. The difference between the fitted elevation and the actual elevation was used to retain ground points and filter out ground features. JM Sánchez et al. proposed a ground filtering algorithm based on scan line processing. The core of this algorithm is to use a one-dimensional iterative spline interpolation algorithm to improve ground point extraction.
[0005] However, the filtering models based on one-dimensional terrain features mentioned above are mostly suitable for extracting ground points when the terrain is flat. They are difficult to apply to areas with rugged terrain and complex scenes. They cannot avoid situations where abrupt changes in elevation caused by fault lines lead to the incorrect filtering of a large number of ground points. Summary of the Invention
[0006] Purpose of the invention: To address the above shortcomings, this invention provides a fast airborne point cloud filtering method based on fracture line constraints. It can build a scan line filtering model that considers fracture line constraints based on the advantages of scan line filtering. It is applicable to a variety of complex scenarios and takes into account the influence of elevation changes and noise points, thereby improving the completeness and accuracy of ground point extraction.
[0007] Technical Solution: To solve the above problems, this invention discloses an airborne point cloud fast filtering method based on fracture line constraints, specifically including the following steps:
[0008] (1) Collect airborne laser point cloud data within the measurement area, and extract the fracture points on the fracture line based on the average curvature and elevation values of each point cloud and its corresponding neighboring point clouds.
[0009] (2) Determine the abscissa of each scan line according to the range of the measurement area and the preset scan line spacing; obtain the abscissa and ordinate of each point cloud, calculate the distance of each point cloud from each scan line in the horizontal axis direction, assign the point cloud located within the set distance to the corresponding scan line, and then arrange the point clouds on the scan line in ascending order according to the ordinate value.
[0010] (3) Perform coarse filtering on each scan line to filter out point clouds that are not ground points and ground features on the scan line and extract the scan center point and reference point; the scan center point is the point cloud at the leftmost end of the scan line after coarse filtering, and the reference point is the point cloud closest to the scan center point;
[0011] (4) Perform fine filtering for each scan line, wherein the fine filtering specifically includes:
[0012] (4.1) Obtain the point cloud adjacent to the reference point B according to the point cloud order on the scan line as the test point, and determine whether the test point is a break point. If it is a break point, retain the point cloud; if it is not a break point, construct a directed vector with the test point and calculate the angle between the vectors. ,like If the value is not greater than the threshold, the test point is retained; otherwise, the test point is filtered out. The specific formula is:
[0013]
[0014] In the formula, B is the reference point. The point to be measured; The center point of the scan;
[0015] (4.2) If the above test point is retained, the retained test point is replaced with the reference point, and step (4.1) is repeated until all point clouds on the scan line are traversed; if the above test point is filtered out, the reference point B is not replaced, and step (4.1) is repeated until all point clouds on the scan line are traversed.
[0016] (5) Use the quadratic surface fitting method to filter out point clouds that are not ground points in the measurement area and obtain the final filtered point cloud.
[0017] Furthermore, the extraction of the break point in step (1) specifically includes the following steps:
[0018] (1.1) Calculate the average curvature value H of each point cloud in the measurement area and its neighboring point clouds; the specific formula is:
[0019]
[0020] In the formula, , , , , , ; ; , , , , All The partial derivatives; The expression for fitting a surface to any point cloud and its neighborhood points using the least squares method and then parameterizing it is given by a, b, and c, where a, b, and c are the fitting parameters.
[0021] (1.2) For each point cloud in the measurement area, calculate the ratio G of the number of neighboring points with elevations lower than that point to the total number of points in the entire area;
[0022] (1.3) Calculate the eigenvalue D of each point cloud based on H and G, and retain the point clouds with eigenvalue D greater than the preset value as breakpoints; the specific formula is as follows:
[0023]
[0024] In the formula, Weights for Type I breakpoints. The weights are for Type II breakpoints.
[0025] Furthermore, the coarse filtration described in step (3) specifically includes the following steps:
[0026] (3.1) Define the leftmost point cloud on each scan line as the initial point. Select the point cloud adjacent to the initial point according to the point cloud order on the scan line, calculate the height difference between the point cloud and the initial point, and if the height difference is less than the set height difference threshold, then the point cloud is a ground point and is denoted as [missing information]. Otherwise, filter out the point and continue selecting subsequent new points according to the point cloud sequence on the scan line, repeating the above operation until the ground point is determined. ;
[0027] (3.2) Selection and Adjacent point clouds are denoted as ,by , Construct equilateral triangles of different lengths. equilateral triangle ;like If the angle is less than the included angle threshold, then it is determined that... The point is a ground point; if If the angle is greater than the included angle threshold, then it is determined that... The point is a non-ground point, so filter it out. And continue selecting subsequent new points as new points according to the point cloud sequence on the scan line. Reconstruct the equilateral triangle and obtain the included angle. Then Compare with the included angle threshold until a ground point is obtained. point; Take 1...n;
[0028] (3.3) Let Increment by 1, repeat step (3.2) until the next new point is the last point cloud on the scan line, filter out non-ground points on the scan line, and determine all ground points on the scan line;
[0029] (3.4) Calculate the flatness of each point cloud on the scan line and retain the point clouds whose flatness is greater than the flatness threshold.
[0030] Furthermore, step (5) specifically includes the following steps:
[0031] (5.1) Divide the original point cloud data into grids based on the range of horizontal and vertical coordinates of each point cloud, forming overlapping and non-overlapping regions between different grids; obtain m ground points on the filtered scan line falling within the cell, and use the least squares principle to calculate the sum of squared errors. Surface parameters of the cell to reach the minimum , , , , , The specific formula is:
[0032]
[0033] In the formula, The coordinates of the point cloud;
[0034] (5.2) After determining the parameters of each grid surface, determine the grid position of each point cloud in the original point cloud: if the original point cloud falls into the non-overlapping area of the grid, then calculate the fitted elevation of the point cloud in a single grid. If the original point cloud falls within the grid overlap region, then calculate the fitted elevation of the point cloud on different grid surfaces. The average value; The formula obtained is:
[0035]
[0036] (5.3) Compare the difference between the fitted elevation value and the actual elevation of each point cloud, filter out the point clouds with a difference greater than the preset threshold, and obtain the final filtered point cloud.
[0037] Furthermore, the coarse filtration described in step (3) also includes:
[0038] (3.5) For any point cloud in the point cloud, obtain the ratio of the number of neighboring points with an elevation greater than that point to the total number of neighboring points, and retain point clouds with a ratio less than a set threshold.
[0039] Furthermore, in step (3.2), the geometric relations of an equilateral triangle are calculated to obtain... The specific formula for defining coordinates is as follows:
[0040]
[0041]
[0042] In the formula, Represented as ( , , ), Represented as ( , , ), Represented as ( , , ), Represented as ( , , ).
[0043] Furthermore, the formula for calculating the flatness of the point cloud in step (3.4) is as follows:
[0044]
[0045]
[0046] In the formula, Represented as the neighborhood of the point cloud. ; Represented as a covariance matrix; The center point of the spatial neighborhood; Point cloud flatness; , , These are the eigenvalues of the covariance matrix D in three-dimensional space.
[0047] Furthermore, the present invention also provides a computer device, including a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor, when executing the computer program, implements the steps of any of the above-described methods. A computer-readable storage medium storing a computer program thereon, wherein the computer program, when executed by a processor, implements the steps of any of the above-described methods.
[0048] Beneficial Effects: Compared with existing technologies, the airborne point cloud fast filtering method based on fracture line constraints described in this invention has the following significant advantages: The method combines curvature values and elevation differences to obtain fracture points, making it suitable for filtering different types of survey areas; for airborne laser point cloud data in various scenarios, it considers constructing equilateral triangles to calculate slope values and point cloud flatness for coarse filtering; on the other hand, it incorporates fracture point constraints to overcome the influence of elevation abrupt changes and noise points, achieving fine filtering; finally, it uses a quadratic surface fitting method to filter out non-ground points within the measurement area, improving the completeness and accuracy of ground point extraction and forming a digital elevation model. Attached Figure Description
[0049] Figure 1 The diagram shown is a flowchart of the method described in this invention;
[0050] Figure 2 The diagram shown illustrates different fracture lines in the method described in this invention. Figure 2 (a) is a schematic diagram of a Type I fracture line. Figure 2 (b) is a schematic diagram of a Type II fracture line;
[0051] Figure 3 The image shown is a result of extracting the fracture line points using the method described in this invention.
[0052] Figure 4 The diagram shown is a schematic diagram of the scan line coarse filtering method of the present invention;
[0053] Figure 5 The diagram shown is a schematic of the vectors constructed in the method of the present invention and the angle between the vectors;
[0054] Figure 6 The diagram shown is a comparison of the results before and after fine filtering in the method described in this invention. Figure 6 (a) is a schematic diagram before fine filtering. Figure 6 (b) is a schematic diagram after fine filtering;
[0055] Figure 7 The diagram shown is a schematic representation of an example of grid division of the measurement area in the method described in this invention;
[0056] Figure 8 The diagram shows the ground point cloud filtering results in different types of complex scenarios in the method described in this invention. Figure 8 (a) is a schematic diagram after terrain 1 filtering. Figure 8 (b) is a schematic diagram after terrain 2 filtering. Figure 8 (c) is a schematic diagram of the terrain after 3-stage filtering. Figure 8 (d) is a schematic diagram of the terrain after 4-filtering. Detailed Implementation
[0057] The technical solution of the present invention will be further described below with reference to the accompanying drawings.
[0058] like Figure 1 As shown, the present invention provides a fast airborne point cloud filtering method based on fracture line constraints, which specifically includes the following steps:
[0059] Step 1: Collect airborne laser point cloud data within the measurement area, and extract the fracture points on the fracture line by combining elevation differences and curvature values.
[0060] like Figure 2 As shown in (a), fault lines include Class I fault lines with continuous curved surfaces at ridgelines, valley lines, etc.; such as Figure 2 As shown in (b), the fault lines also include Type II fault lines where the curved surface is discontinuous at artificial structures, steep slopes, etc. Since the fault points on Type I fault lines are continuous, Type I fault points can be directly extracted based on the absolute value of the average curvature of the point cloud. The absolute value of the average curvature reflects the degree of surface undulation. When this value is large, it indicates that the local surface undulation is large, and vice versa, it indicates that the ground is relatively flat. Points with larger absolute values of average curvature can be extracted as potential fault points. Since the ground points at Type II fault lines are discontinuous and there are significant elevation differences, directly extracting this type of fault line based on curvature is not very effective. Type II fault points can be extracted by statistically analyzing the ratio of the number of neighboring points with elevations lower than the point to the total number of points in the entire region. If the ratio is greater than a given threshold, the point is judged as a potential fault point.
[0061] Because different methods are used to extract fracture points on different types of fracture lines, this method combines elevation difference and curvature value to extract fracture points in order to be applicable to the filtering of measurement areas containing different types of fracture lines. Specifically, it includes the following steps:
[0062] (1) For any point cloud and its neighborhood point clouds, fit the surface using the least squares method and parameterize it to obtain: a, b, and c are all fitting parameters; The partial derivatives are respectively , , , , ; Calculate the average curvature value H of each point cloud in the measurement area and its neighboring point clouds; the specific formula is:
[0063]
[0064] In the formula, , , , , , L, M, and N are the second fundamental invariants of the surface; E, F, and G are the second fundamental invariants of the surface. The direction of the principal normal to the surface;
[0065] (2) For each point cloud in the measurement area, calculate the ratio G of the number of neighboring points with an elevation lower than that point to the total number of points in the entire area;
[0066] (3) Calculate the eigenvalues D of each point cloud based on H and G, using the following formula:
[0067]
[0068] In the formula, Weights for Type I breakpoints. The weights are for Type II breakpoints.
[0069] The feature value D of each point cloud is compared with a preset value. If the feature value D is greater than the preset value (set to 0.6 in this embodiment), the point cloud is considered a break point; otherwise, it is considered a non-break point. Based on the above steps, the break lines formed by each break point are finally obtained as follows: Figure 3 As shown in the figure, the results show that the fracture line basically conforms to the actual situation and can meet the filtering constraint requirements.
[0070] Step 2: Extract one-dimensional terrain feature points from the point cloud of the entire area to obtain the x-coordinate and y-coordinate of each point cloud. Determine the x-coordinate of each scan line based on the measurement area and the preset scan line spacing; calculate the distance of each point cloud from each scan line in the horizontal direction, assign the point clouds within the set distance to the corresponding scan line, and arrange the point clouds on the scan line in ascending order of their y-coordinate values.
[0071] Step 3: Perform coarse filtering on each scan line to remove point clouds that are not ground points or ground features, and select the scan center. This coarse filtering includes various techniques, such as determining whether adjacent points are non-ground points by judging whether their elevation difference and slope are within a preset threshold range. In this method, the coarse filtering first uses the points on the scan line to construct an equilateral triangle to calculate the slope and filter most ground features, then removes potential building points based on ground flatness, and finally filters non-ground points based on elevation differences. Specific steps include:
[0072] (1) Use the points on the scan line to construct an equilateral triangle to calculate the slope and filter most ground features.
[0073] (1.1) Define the leftmost point cloud on each scan line as the initial point. Select the point cloud adjacent to the initial point according to the point cloud order on the scan line, calculate the height difference between the point cloud and the initial point, and if the height difference is less than the set height difference threshold h, then the point cloud is a ground point and is denoted as h. Otherwise, filter out the point and continue selecting subsequent new points according to the point cloud sequence on the scan line, repeating the above operation until the ground point is determined. ;
[0074] (1.2) Selection and Adjacent point clouds are denoted as ,by , Construct equilateral triangles of different lengths. equilateral triangle ;like If the angle is less than the included angle threshold, then it is determined that... The point is a ground point; if If the angle is greater than the included angle threshold, then it is determined that... The point is a non-ground point, so filter it out. And continue selecting subsequent new points as new points according to the point cloud sequence on the scan line. Reconstruct the equilateral triangle and obtain the included angle. Then Compare with the included angle threshold until a ground point is obtained. point; Take 1...n;
[0075] Specifically, such as Figure 4 As shown, select the ground point Adjacent point clouds are denoted as ,by , Construct an equilateral triangle using length With equilateral triangle ;like If the angle is less than the included angle threshold, then it is determined that... The point is a ground point; if If the angle is greater than the included angle threshold, then it is determined that... The point is a non-ground point, so filter it out. And continue selecting subsequent new points as new points according to the point cloud sequence on the scan line. Reconstruct an equilateral triangle to obtain the included angle. And compare it with the included angle threshold until a ground point is obtained. point;
[0076] Among them, it is obtained by calculation from the geometric relations of equilateral triangles. The coordinates can be calculated similarly. Coordinates. The specific formula is:
[0077]
[0078]
[0079] In the formula, the current point ( , , ), ( , , ), ( , , ), ( , , ).
[0080] (1.3) Let Increment by 1, repeat step (3.2) until the next new point is the last point cloud on the scan line. Traverse all points on the scan line and filter out non-ground points on the scan line using this method to determine all ground points on the scan line.
[0081] (2) Eliminate potential building points based on ground flatness. Calculate the point cloud flatness for each point cloud along the scan line, and retain point clouds with a flatness greater than the flatness threshold. The specific formula for calculating point cloud flatness is as follows:
[0082]
[0083]
[0084] In the formula, Represented as the neighborhood of the point cloud. ; Represented as a covariance matrix; The center point of the spatial neighborhood; Point cloud flatness; , , ( ≤ ≤ ) are the eigenvalues of the covariance matrix D in three-dimensional space.
[0085] (3) After filtering out vegetation, vehicles, people and other ground features on the scan line, the influence of non-ground points with flat local neighborhoods, such as buildings, on the selection of the scan center needs to be considered. In order to remove buildings on the scan line, for any point in the point cloud on the scan line, the ratio of the number of neighboring points with an elevation greater than that point to the total number of neighboring points is obtained. The smaller the ratio, the greater the possibility that the point is a ground point. Therefore, point clouds with a ratio less than the set threshold are retained.
[0086] After coarse filtering for each scan line, the leftmost endpoint of the scan line is taken as the scan center O, and the point closest to the scan center is taken as the reference point B.
[0087] Step 4: Consider the influence of the break point on the fracture line and perform fine filtering on the scan line. Specifically, this includes:
[0088] (1) Obtain the point cloud adjacent to the reference point B according to the point cloud order on the scan line as the test point, and determine whether the test point is a break point. If it is a break point, retain the point cloud; if it is not a break point, construct a directed vector with the test point and calculate the angle between the vectors. ,like If the value is not greater than the threshold, the test point is retained; otherwise, the test point is filtered out. The specific formula is:
[0089]
[0090] In the formula, B is the reference point. The point to be measured; The center point of the scan;
[0091] (2) If the above test point is retained, the retained test point is replaced with the reference point, and step (1) is repeated until all point clouds on the scan line are traversed; if the above test point is filtered out, the reference point B is not replaced, and step (1) is repeated until all point clouds on the scan line are traversed.
[0092] like Figure 5 As shown, the current reference point is Points, construct directed vectors , Obtain the angle between vectors If the angle between the vectors is less than the threshold, it is considered a non-ground point and is filtered out. If a point is identified as a non-ground point, the reference point is kept unchanged, and new points are selected to reconstruct the directed vector. The determination of whether a new point to be measured is a ground point is then repeated until the determination of all point clouds along the scan line is completed. Figure 6 As shown, Figure 6 (a) is the original point cloud data image. Figure 6 (b) is the point cloud image data after fine filtering, which can distinguish between ground points and non-ground points in the image.
[0093] Step 5: Use the quadratic surface fitting method to filter out point clouds that are not ground points within the measurement area. That is, filter the point clouds between the scan lines to obtain the final filtered point cloud and form a digital elevation model.
[0094] (1) Divide the original point cloud data into grids based on the preset horizontal and vertical coordinate ranges of each point cloud. The grid cell size should, in principle, be larger than the side length of the largest feature. Overlapping and non-overlapping areas are formed between different grids; for example... Figure 7 As shown, after all grids are divided within the measurement area, there are non-overlapping areas, two-grid overlapping areas, and four-grid overlapping areas.
[0095] Obtain the m ground points on the filtered scan line that fall within the cell, and calculate the sum of squared errors using the least squares principle. Approximate surface parameters of the cell to reach the minimum , , , , , The specific formula is:
[0096]
[0097] In the formula, The coordinates of the point cloud;
[0098] (2) After determining the parameters of each grid surface, determine the grid position of each point cloud in the original point cloud: if the original point cloud falls into the non-overlapping area of the grid, calculate the fitted elevation of the point cloud in a single grid. If the original point cloud falls within the grid overlap region, then calculate the fitted elevation of the point cloud on different grid surfaces. The average value; The formula obtained is:
[0099]
[0100] (3) Compare the difference between the fitted elevation value and the actual elevation of each point cloud, filter out the point clouds with a difference greater than the preset threshold, and obtain the final filtered point cloud belonging to the ground points.
[0101] like Figure 8 As shown, the method described in this invention is used to extract ground points for different types of measurement areas; for example... Figure 8 (a) 8 (b) Figure 8 (c) and (d) are the filtered results for four different terrains. According to the figures, point cloud data of all ground points in the observation area can be extracted.
Claims
1. A fast airborne point cloud filtering method based on fracture line constraints, characterized in that, Includes the following steps: (1) Collect airborne laser point cloud data within the measurement area, and extract the fracture points on the fracture line based on the average curvature and elevation values of each point cloud and its corresponding neighboring point clouds. (2) Determine the abscissa of each scan line according to the range of the measurement area and the preset scan line spacing; obtain the abscissa and ordinate of each point cloud, calculate the distance of each point cloud from each scan line in the horizontal axis direction, assign the point cloud located within the set distance to the corresponding scan line, and then arrange the point clouds on the scan line in ascending order according to the ordinate value. (3) Perform coarse filtering on each scan line to filter out point clouds that are not ground points and ground features on the scan line and determine the scan center point and reference point; the scan center point is the point cloud at the leftmost end of the scan line after coarse filtering, and the reference point is the point cloud closest to the scan center point; (4) Perform fine filtering for each scan line, wherein the fine filtering specifically includes: (4.1) Obtain the point cloud adjacent to the reference point according to the point cloud order on the scan line as the test point, and determine whether the test point is a break point. If it is a break point, retain the point cloud; if it is not a break point, construct a directed vector with the test point and calculate the angle between the vectors. ,like If the value is not greater than the threshold, the test point is retained; otherwise, the test point is filtered out. The specific formula is: ； In the formula, B is the reference point. The point to be measured; The center point of the scan; (4.2) If the above test point is retained, the retained test point is replaced with the reference point, and step (4.1) is repeated until all point clouds on the scan line are traversed; if the above test point is filtered out, the reference point B is not replaced, and step (4.1) is repeated until all point clouds on the scan line are traversed. (5) Use the quadratic surface fitting method to filter out point clouds that are not ground points in the measurement area and obtain the final filtered point cloud.
2. The airborne point cloud fast filtering method based on fracture line constraints according to claim 1, characterized in that, Extracting the break point in step (1) specifically includes the following steps: (1.1) Calculate the average curvature value H of each point cloud in the measurement area and its neighboring point clouds; the specific formula is: ； In the formula, , , , , , ; ; , , , , All The partial derivatives; The expression for fitting a surface to any point cloud and its neighborhood points using the least squares method and then parameterizing it is given by a, b, and c, where a, b, and c are the fitting parameters. (1.2) For each point cloud in the measurement area, calculate the ratio of the number of neighboring points with elevations lower than that point to the total number of neighboring points. ; (1.3) Based on the H and Calculate the feature value D for each point cloud, and retain the point clouds with feature value D greater than a preset value as breakpoints; the specific formula is as follows: ； In the formula, Weights for Type I breakpoints. The weights are for Type II breakpoints.
3. The airborne point cloud fast filtering method based on fracture line constraints according to claim 1, characterized in that, The coarse filtration described in step (3) specifically includes the following steps: (3.1) Define the leftmost point cloud on each scan line as the initial point. Select the point cloud adjacent to the initial point according to the point cloud order on the scan line, calculate the height difference between the point cloud and the initial point, and if the height difference is less than the set height difference threshold, then the point cloud is a ground point and is denoted as [missing information]. Otherwise, filter out the point and continue selecting subsequent new points according to the point cloud sequence on the scan line, repeating the above operation until the ground point is determined. ; (3.2) Selection and Adjacent point clouds are denoted as ,by , Construct equilateral triangles of different lengths. equilateral triangle ;like If the angle is less than the included angle threshold, then it is determined that... The point is a ground point; if If the angle is greater than the included angle threshold, then it is determined that... The point is a non-ground point, so filter it out. And continue selecting subsequent new points as new points according to the point cloud sequence on the scan line. Reconstruct the equilateral triangle and obtain the included angle. Then Compare with the included angle threshold until a ground point is obtained. point; Take values from 1 to n, where n is a positive integer; (3.3) Let Increment by 1, repeat step (3.2) until the next new point is the last point cloud on the scan line, filter out non-ground points on the scan line, and determine all ground points on the scan line; (3.4) Calculate the flatness of each point cloud on the scan line and retain the point clouds whose flatness is greater than the flatness threshold.
4. The airborne point cloud fast filtering method based on fracture line constraints according to claim 1, characterized in that, Step (5) specifically includes the following steps: (5.1) Divide the original point cloud data into grids based on the range of horizontal and vertical coordinates of each point cloud, forming overlapping and non-overlapping regions between different grids; obtain m ground points on the filtered scan line falling within the cell, and use the least squares principle to calculate the sum of squared errors. Surface parameters of the cell to reach the minimum , , , , , The specific formula is as follows: ； In the formula, The coordinates of the point cloud; (5.2) After determining the parameters of each grid surface, determine the grid position of each point cloud in the original point cloud: if the original point cloud falls into the non-overlapping area of the grid, then calculate the fitted elevation of the point cloud in a single grid. If the original point cloud falls within the grid overlap region, then calculate the fitted elevation of the point cloud on different grid surfaces. The average value; The formula obtained is: ； (5.3) Compare the difference between the fitted elevation value and the actual elevation of each point cloud, filter out the point clouds with a difference greater than the preset threshold, and obtain the final filtered point cloud.
5. The airborne point cloud fast filtering method based on fracture line constraints according to claim 3, characterized in that, The coarse filtration described in step (3) also includes: (3.5) For any point cloud, calculate and count the ratio of the number of neighboring points with an elevation greater than that point to the total number of neighboring points, and retain point clouds with a ratio less than a set threshold.
6. The airborne point cloud fast filtering method based on fracture line constraints according to claim 3, characterized in that, In step (3.2), the geometric relations of an equilateral triangle are calculated to obtain... The specific formula for defining coordinates is as follows: ； ； In the formula, Represented as ( , , ), Represented as ( , , ), Represented as ( , , ), Represented as ( , , ).
7. The airborne point cloud fast filtering method based on fracture line constraints according to claim 3, characterized in that, The formula for calculating the flatness of the point cloud in step (3.4) is: ； ； In the formula, Represented as the neighborhood of the point cloud. ; Represented as a covariance matrix; The center point of the spatial neighborhood; Point cloud flatness; , , These are the eigenvalues of the covariance matrix D in three-dimensional space.
8. A computer device, comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, characterized in that, When the processor executes the computer program, it implements the steps of the method according to any one of claims 1 to 7.
9. A computer-readable storage medium having a computer program stored thereon, characterized in that, When the computer program is executed by a processor, it implements the steps of the method according to any one of claims 1 to 7.