A method for identifying wind erosion pits based on point cloud data

Abstract

The application discloses a kind of wind erosion pit identification methods based on point cloud data, obtains the point cloud data of the place to the laser radar sensor acquisition, extracts ground category point cloud;The ground category point cloud is converted into raster data, and the missing part of raster data is completed, to generate digital elevation model;Fill in the recessed area in the digital elevation model;Difference is obtained by the digital elevation model before and after filling recessed area, the difference part is identified wind erosion pit.The application can improve the wind erosion pit identification effect under grassland environment.

Description

Technical Field

[0001] This invention belongs to the fields of remote sensing and geological geography, and specifically relates to a method for identifying wind erosion pits based on point cloud data. Background Technology

[0002] Wind erosion pits are saucer-shaped, bowl-shaped, or trough-shaped depressions formed by the erosion of underlying sandy deposits by wind after surface vegetation has been damaged. The formation of wind erosion pits contributes to a series of problems, including grassland desertification. Preventing further deterioration of wind erosion pits is one of the key tasks in grassland desertification control, and how to achieve automatic identification and mapping of wind erosion pits is the primary issue in grassland desertification control.

[0003] Currently, research on wind erosion pits by scholars both domestically and internationally focuses on their dynamic processes, morphological evolution, surface airflow, formation mechanisms, grain size characteristics, morphological changes, and impacts. However, there are relatively few remote sensing image recognition methods for wind erosion pits in grasslands. Due to the unique environment of grasslands, and the lack of three-dimensional information in remote sensing images, which are also significantly affected by weather conditions, the accuracy of identifying wind erosion pits using traditional remote sensing images is not high. Extensive supplementary calculations are often required in post-processing, and even then, it is difficult to completely and effectively label wind erosion pits. Summary of the Invention

[0004] In order to overcome the shortcomings of the prior art, the purpose of this invention is to provide a wind erosion pit identification method based on point cloud data, so as to improve the wind erosion pit identification effect in grassland environment.

[0005] To achieve the above objectives, the technical solution adopted by the present invention is as follows:

[0006] A method for identifying wind erosion pits based on point cloud data includes the following steps:

[0007] Step 1: Obtain point cloud data acquired by the lidar sensor and extract point clouds of ground categories;

[0008] Step 2: Convert the point cloud of the ground category into raster data, and fill in the missing parts of the raster data to generate a digital elevation model;

[0009] Step 3: Fill in the concave areas in the digital elevation model;

[0010] Step 4: By subtracting the digital elevation models before and after filling the depression area, the difference is obtained, which is the identified wind erosion pit.

[0011] In one embodiment, step 1 involves using a clustering algorithm to classify the acquired point cloud in order to extract point clouds of the ground category.

[0012] In one embodiment, a clustering algorithm is used to cluster the acquired point cloud data according to its color attribute RGB, dividing it into three categories. The R and B categories of the clustering results are modified to the ground category in the LAS data, and the G category of the clustering results is modified to the vegetation category in the LAS data. The point clouds of the ground category are saved, and the other point clouds are deleted.

[0013] In one embodiment, converting the point cloud of the ground category into raster data involves using code to read the LAS file, mapping the point cloud into a grid, and describing the grid using a 2D array. If multiple points are mapped to the same grid, the average value of the multiple points is taken as the value of that grid.

[0014] In one embodiment, step 2 uses the neighborhood averaging method to complete the missing parts of the raster data. That is, the average value of the eight adjacent non-empty raster cells around the missing raster is taken as the value of the missing raster, and the missing cells are searched from the top left to the bottom right of the entire raster network to complete all missing raster cells.

[0015] In one embodiment, step 3, the filling method is as follows:

[0016] Step 31: Convert the DEM into a 2D array elevation, where each row of the elevation array is a 1D array, and search each row of the elevation array.

[0017] Step 32: Search point by point from left to right to determine the current point X. k Elevation value e(X) k Is it less than or equal to the next point X? k+1 Elevation value e(X) k+1 If so, increment the value of k by 1 and continue the judgment, repeating this process until e(X) is found. k )>e(X k+1 ), and place point X at this time k Let it be labeled, and the corresponding elevation value is e(label). The value of k is incremented by 1, where k is the index of the element in the 1D array.

[0018] Step 33, determine e(X) k If the condition is true, increment the value of k by 1 and continue the check until e(X) < e(label). k If )≥e(label), then the point label is compared with the current point X. k The elevation values ​​between them are all modified to e(label), that is, to fill in the DEM, and the value of k is increased by 1;

[0019] Step 34: Return to step 32 until the last point of the current line in the elevation is reached;

[0020] Step 35: Search point by point from right to left to determine the current point X. k Elevation value e(X) k Is it less than or equal to the next point X? k-1 Elevation value e(X) k-1 If so, decrement the value of k by 1 and continue the judgment, repeating this process until e(X) is found. k )>e(X k-1 ), and place point X at this time k Let it be labeled 2, and the corresponding elevation value is e(label2), and the value of k is reduced by 1;

[0021] Step 36, determine e(X) k If the condition ) < e(label2) is true, then the value of k is decremented by 1, and the process continues until e(X) < e(label2). k If )≥e(label2), then the point label2 and the current point X are compared. k The elevation values ​​between them are all modified to e(label2), which means filling in the DEM, and the value of k is reduced by 1;

[0022] Step 37: Return to step 35 until the first point of the current line in the elevation is reached;

[0023] Step 38: Name the resulting 2D array of the filled DEM as elevation_filled.

[0024] In one embodiment, step 4 involves subtracting each element in the array elevation_filled from each element in the array elevation to obtain another 2D array named elevation_difference, which represents the filling done on the DEM. The non-zero areas in elevation_difference are the locations of the wind erosion pits.

[0025] The present invention also provides a wind erosion pit identification system based on point cloud data, including a lidar sensor, a memory, a processor, and a display; the lidar sensor is mounted on a navigation device and is used to acquire point cloud data; the memory is used to store a computer program, and the processor is used to call and run the computer program stored in the memory to execute the method described in the present invention; the display shows the execution result of the processor.

[0026] The present invention also provides a computer-readable storage medium for storing a computer program that causes a computer to perform the methods described in the present invention.

[0027] The present invention also provides a computer program product containing instructions that, when executed on a computer, cause the computer to perform the method described in the present invention.

[0028] Compared with the prior art, the beneficial effects of the present invention are:

[0029] Compared to traditional remote sensing data, such as satellite imagery and aerial photographs, the point cloud data used in this invention provides detailed three-dimensional information, making object identification and classification, such as buildings, trees, roads, and water bodies, easier. It also allows for consistent terrain data across different weather conditions and seasons. This invention's wind erosion pit identification method based on point cloud data combines initial automatic classification using algorithms with manual review for refinement. This approach requires minimal workload and yields good classification results. Filling DEM depressions involves two rounds of searching, from left to right and from right to left, simplifying the process and avoiding the problem of not being able to fill depressions at the DEM edges. It accurately identifies DEM depressions, thus effectively marking wind erosion pits. Attached Figure Description

[0030] Figure 1 This is a flowchart of the identification method of the present invention.

[0031] Figure 2 It is the process of completing the missing data in the DEM.

[0032] Figure 3 It is the process of filling and subtracting specific latitude depressions in a DEM.

[0033] Figure 4 The results are obtained by identifying wind erosion pits using the method proposed in this invention. Detailed Implementation

[0034] The embodiments of the present invention will now be described in detail with reference to the accompanying drawings and examples.

[0035] As mentioned earlier, the grassland surface is relatively flat and uniform in color, and the climate is also changeable, making it difficult to identify wind erosion pits using traditional remote sensing images.

[0036] With the increasing prevalence of lidar equipment and data, the application range of lidar has expanded, providing convenience for this invention. Digital terrain elevation information can be obtained from point cloud data, and wind erosion pits can be identified by analyzing terrain changes and terrain attributes.

[0037] like Figure 1 As shown, the present invention provides a method for identifying wind erosion pits based on point cloud data, which mainly includes the following steps:

[0038] Step 1: Perform point cloud classification.

[0039] The present invention first requires acquiring point cloud data obtained by a lidar sensor and extracting point clouds of ground categories from it.

[0040] In one embodiment of the present invention, a clustering algorithm is used to classify the acquired point cloud, and the classification results are checked and processed. Finally, only the point cloud of the ground category is saved, and other point clouds are deleted, thereby extracting the point cloud of the ground category. Specifically, the clustering algorithm is used to cluster the acquired point cloud data according to its color attribute RGB, dividing it into three categories. The R and B categories of the clustering results are modified to the "ground" category in the LAS data, and the G category of the clustering results is modified to the "vegetation" category in the LAS data. The automatic point cloud classification results can be manually checked and improved. Finally, only the point cloud of the ground category is saved, and other point clouds are deleted.

[0041] Step 2: Convert the point cloud of the ground category into raster data, and fill in the missing parts of the raster data to generate a Digital Elevation Model (DEM).

[0042] In one embodiment of the present invention, the conversion of the point cloud of the ground category into raster data is achieved by using code to read the LAS file, mapping the point cloud into a grid, and describing the grid using a 2D array. If multiple points are mapped to the same grid, the average value of the multiple points is taken as the value of the grid.

[0043] Subsequently, the missing parts of the raster data are completed using the neighborhood averaging method. This involves taking the average value of the eight adjacent non-empty raster cells surrounding the missing raster as the value of that missing raster, and searching for missing cells from the top left to the bottom right of the entire raster network to complete all missing raster cells. Existing technologies primarily use software to directly generate DEMs, and the interpolation methods used to complete raster data are mainly inverse distance weighting, kriging, and natural neighbor interpolation, which differ from this invention. When the terrain or topography is complex, the distance differences between points become larger, potentially leading to significant errors in the interpolation results. To address the problems associated with inverse distance weighting, this invention employs the neighborhood averaging method, a simplified version of the inverse distance weighting method that only considers points at a distance of 1, reducing the impact of terrain undulations.

[0044] In this step, point cloud data is converted to raster. When multiple point clouds exist within a raster, the average height of these point clouds is taken. Because non-ground points are deleted, some point clouds are missing. For raster cells without point clouds, the average height of the eight surrounding non-empty raster cells is taken as the value of that raster cell to fill in the missing DEM. For example... Figure 2 As shown.

[0045] Step 3: Fill in the concave areas in the digital elevation model.

[0046] For example, the filling method used in this invention is as follows:

[0047] Step 31: Convert the digital elevation model into a 2D array `elevation`, where each row of `elevation` is a 1D array. Search each row of `elevation`; this row corresponds to the elevation data at a specific latitude of the actual study area. Plot the elevation data at this latitude as a curve, such as... Figure 3 As shown in (a).

[0048] Step 32: Search point by point from left to right to determine the current point X. k Elevation value e(X) k Is it less than or equal to the next point X? k+1 Elevation value e(X) k+1 If so, it means we are currently in an uphill phase. Increment the value of k by 1, which means moving the current point one point to the right. Continue this process until e(X) is reached. k )>e(X k+1 ),like Figure 3 Point A, as shown in (a), is now entering the downhill phase. Point X at this point... k Let k be denoted as label, and the corresponding elevation value is e(label). This point is a maximum point, and the value of k is incremented by 1. Here, k is the index of the element in the one-dimensional array.

[0049] Step 33, set the current elevation value e(X) k Compare with the elevation maximum point e(label): that is, determine e(X) k If k < e(label), then increment the value of k by 1. Figure 3 The region between points A and B shown in (a) is further evaluated, and this loop is repeated until e(X) is reached. k If ) ≥ e(label), it means that a point higher than the point label has appeared, such as... Figure 3 Point B is shown in (a). At this point, the area between points A and B is the location of the wind erosion pit. Compare point label (i.e., point A) with the current point X. k The elevation values ​​between points (i.e., point B) are all modified to the elevation value e(label) of point A, i.e., filling the DEM, as shown by the dashed line AB, and the value of k is increased by 1.

[0050] Step 34: Return to step 32 and determine the current point X. k Elevation value e(X) k Is it less than or equal to the next point X? k+1 Elevation value e(X) k+1 (, until the last point of the current line in the elevation is reached, such as...) Figure 3 Point F is shown in (a).

[0051] Step 35: The situation between points C and D is similar to that between points A and B, but at point E, the right side is about to reach the edge of the DEM, and there are no points higher than point E. The above method cannot fill this point, but as shown in the diagram, it should be filled. Therefore, after searching point by point from left to right, we then search point by point from right to left to determine the current point X. k Elevation value e(X) k Is it less than or equal to the next point X? k-1 Elevation value e(X) k-1 If so, decrement the value of k by 1 and continue the judgment, repeating this process until e(X) is found. k )>e(X k-1 ), and place point X at this time k Let this be labeled 2, with the corresponding elevation value being e(label2), and the value of k minus 1; in this embodiment, since... Figure 3 Point F, as shown in (a), is initially in a downhill position, satisfying e(X). k )>e(X k-1 Let point F be labeled as label2.

[0052] Step 36, during the subsequent leftward movement, determine e(X) k If the condition ) < e(label2) is true, then the value of k is decremented by 1, and the process continues until e(X) < e(label2). k )≥e(label2), that is, until Figure 3 Point G is shown in (a). At this point, point label2 (i.e., point F) is compared with the current point X. k The elevation values ​​between points (i.e., G) are all modified to e(label2) of point F, i.e., filling the DEM, as shown by the dashed line FG, and the value of k is reduced by 1.

[0053] Step 37: Return to step 35 until the first point of the current line in the elevation is reached;

[0054] Step 38: Perform the above processing on all rows in elevation to obtain a 2D array of the corresponding filled DEM, named elevation_filled.

[0055] Step 4: DEM differential analysis to identify wind erosion pits.

[0056] In this step, the difference between the digital elevation models before and after filling the depression area is obtained, and this difference is the identified wind erosion pit. Specifically, in step 4, the difference between each element in the array `elevation_filled` and each element in the array `elevation` is obtained to obtain another 2D array named `elevation_difference`, which represents the filling performed on the DEM, corresponding to... Figure 3 The difference results for mid-(a) latitude are as follows Figure 3 As shown in (b), the regions where elevation_difference is non-zero are the locations of wind erosion pits, and the results can be visualized at the end.

[0057] To further verify the method of the present invention, the present invention also provides the following experiment, the steps of which are as follows:

[0058] Step 1, perform point cloud classification:

[0059] The Hunshandake Grassland in Xilingol League, Inner Mongolia, was selected as the study area, with coordinates 42.7631°N115.9153°E. The obtained point cloud data was processed for point cloud classification.

[0060] The point cloud data was clustered using a Python language according to its color attribute (RGB), resulting in three categories. Vegetation points were initially separated from other points by clustering them into the green category (G). MicroStation v8 software was used to modify the classification numbers, and point clouds with incorrect classifications were manually checked and improved. Finally, other point clouds were deleted, and only the ground category point clouds were saved.

[0061] Step 2: Convert the ground point cloud into raster data, fill in the missing parts, and generate a DEM.

[0062] Use the LASpy library to read the LAS file from step 1, and map the points in the LAS to a grid. The grid is selected as a 1m*1m square. When there are two or more points in the same grid, the average height of all points is taken as the elevation value of the grid. For example, if there are three points in a grid with heights of 1305.62, 1305.66, and 1305.64, the average height of 1305.64 is taken as the elevation value of the grid. In step 2, non-ground points were deleted, resulting in missing point clouds. For grid cells without point clouds, the average value of the non-empty grid cells in the surrounding 8 grid cells is taken as the value of this grid cell to complete the DEM. For example, if a grid cell is empty, and the data of its surrounding 8 grid cells are 1305.72, 1305.69, 1305.70, empty, empty, empty, empty, 1305.69, then the average value of these 4 grid cells (1305.72, 1305.69, 1305.70, 1305.69) is taken as its elevation value (1305.70). The final generated DEM size is 515*1142 pixels.

[0063] Step 3: Fill in the recessed areas in the DEM.

[0064] The DEM from step 2 is read using Python code, and the result is stored in a 2D array `elevation`, corresponding to the dimensions of the DEM. This array has 1142 rows and 515 columns, with empty edge data. Its simplified format is as follows:

[0065]

[0066] 1309.03076 is located in row 53 and column 449.

[0067] Change all empty values ​​in the elevation field to the background value -32767, as follows:

[0068]

[0069] Then, two rounds of searching were performed on each row, for a total of 1142 times.

[0070] First, the system searches point by point from left to right, checking if the elevation value of the current point is less than or equal to the elevation value of the next point. If so, it shifts one pixel to the right and continues the check, repeating this process until the elevation value of the current point is greater than that of the next point. This point is then marked as the label. The system continues shifting to the right, constantly comparing the elevation value of the current point with that of the label point. If the current point's elevation value is greater than the label point's, the elevation values ​​between the label point and the current point are changed to the label point's elevation value, effectively filling the DEM. The label marker is then cleared, and the system continues shifting to the right, checking and filling the DEM, until the current row in the elevation is searched. After completing the left-to-right point-by-point search, a right-to-left point-by-point search is performed, following the same principle but in reverse direction.

[0071] For example, in the elevation data, line 309 [X0,X1,…,X…] 448 ,X 449 ,…,X 482 ,…,X 514 ], represented as:

[0072] [-32767,-32767,…,1310.65332,1310.75366,1310.79297,1310.79053,1310.72205,1310.74756,1310.79004,1310.79626,1310.80945,1310.84595,1310.79089,1310.71838,1310.69678,1310.65295,1310.58215,1310.52832,1310.50134,131 0.46191,1310.43713,1310.38171,1310.31494,1310.16675,1309.99866,1309.82202,1309.63379,1309.45081,1309.67700,1310.04321,1310.30176,1310.48682,1310.61670,1310.69617,1310.72449,1310.70044,1310.62524,…,-32767], where X 448 =1310.65332, X 449 =1310.75366, and so on, and X 448 Points on the left: X0, X1, ..., X 447 The value of X shows an increasing trend. 482 Point X on the right 483 ,X 484 ,…,X 514The value of is decreasing.

[0073] The search proceeds from left to right, with k gradually increasing from zero. From X0 to X... 449 That is, when k < 450, e(X) is always satisfied. k )≤e(X k+1 No filling is performed when k = 450 (X). 450 =1310.79297, X 451 =1310.79053), so e(X) k )>e(X k+1 ), and place point X at this time 450 Let this be denoted as label, and its corresponding elevation value is e(label) = 1310.79297. Then, compare the elevation values ​​of each point with the elevation value of label. For any point where 450 < k < 455, e(X) will always be present. k ) < e(label), while when k = 455 (X 455 =1310.79626), satisfying e(X) k If )≥e(label), then the point label(X) will be... 450 ) and X 455 The elevation values ​​between them are all modified to e(label), that is, X 451 ,X 452 ,X 453 and X 454 The values ​​were all changed to 1310.79297, and the data in row 309 became:

[0074] [X0,X1,…,X 448 ,X 449 ,…,X 482 ,…,X 514 ]=[-32767,-32767,…,1310.65332,1310.75366,1310.79297, 1310.79297,1310.79297,1310.79297,1310.79297,1310.79626,1310.80945,1310.84595,1310.79089,1310.71838,1310.69678,1310.65295,1310.58215,1310.52832,1310.50134,1310.46191,1310.43713,1310.38171,1310.31494,1310.166 75,1309.99866,1309.82202,1309.63379,1309.45081,1309.67700,1310.04321,1310.30176,1310.48682,1310.61670,1310.69617,1310.72449,1310.70044,1310.62524,…,-32767], that is, filling the DEM. Then continue searching to the right, when k=456 (X 456 =1310.80945, X 457 =1310.84595) satisfies e(X) k )≤e(X k+1 No filling is performed. When k = 457 (X) 457 =1310.84595, X 458 =1310.79089) has e(X k )>e(X k+1 ), and place point X at this time 457 Let this be the new label point, but when k > 457, e(X) always holds. k If ) < e(label), the filling condition is not met, and the search continues from left to right until line 309 is completed.

[0075] Based on the results of the left-to-right search and filling process described above, a second search is performed from right to left, with k gradually decreasing from 514. From X... 514 To X 481 That is, when k > 480, e(X) is always satisfied. k )≤e(X k-1 No filling is performed when k = 480 (X). 480 =1310.72449, X 479 =1310.69617), we have e(X) k )>e(X k-1 ), and place point X at this time 480 Let this be labeled 'label2', with an elevation value e(label2) = 1310.72449. Then, compare the elevation values ​​of each point with the elevation value of 'label2'. For any point where 458 < k < 480, e(X) will always be present. k ) < e(label2), while when k = 458 (X458 =1310.79089), satisfying e(X) k If )≥e(label2), then the point label2(X) will be set to 0. 480 ) and X 458 The elevation values ​​between them are all modified to e(label2), that is, X 479 ,X 478 ,…,X 459 The values ​​were all changed to 1310.72449, and the data in row 309 became:

[0076] [X0,X1,…,X 448 ,X 449 ,…,X 482 ,…,X 514 ]=[-32767,-32767,…,1310.65332,1310.75366,1310.79297,1310.79297,1310.79297,1310.79297,1310.79297,1310.79626,1310.80945,1310.84595,1310.79089, 1310.72449,1310.72449,1310.72449,1310.72449, 1310.72449,1310.72449,1310.72449,1310.72449,1310.72449,1310.72449,1310.72449, 1310.72449,1310.72449,1310.72449,1310.72449,1310.72449,1310.72449,1310.72449, 1310.72449,1310.72449,1310.72449 [,1310.72449,1310.70044,1310.62524,…,-32767], that is, filling in the DEM. Then continue searching to the left, when k=457 (X 457 =1310.84595, X 456 =1310.80945), we have e(X) k )>e(X k-1 ), and place point X at this time 457 Let this be the new label2 point. However, as we know from the previous left-to-right search process, when k < 457, e(X) always holds. k If ) < e(label2), the filling condition is not met, and the search continues from right to left until the 309th row is completed.

[0077] After all 1142 rows in the elevation array have been searched, a 2D array corresponding to the filled DEM will be obtained, named elevation_filled. In the example above, the modified row 309 corresponds to row 309 in elevation_filled.

[0078] Step 4: DEM differential analysis to identify wind erosion pits.

[0079] The subtraction between each element in the `elevation_filled` array and each element in the `elevation` array yields another 2D array named `elevation_difference`, representing the filling applied to the DEM. The non-zero regions in `elevation_difference` indicate the locations of wind erosion pits. Converting `elevation_difference` into an image output yields the following result: Figure 4 As shown in the figure, the black areas represent the identified wind erosion pits, corresponding to areas where the elevation_difference is not zero. These wind erosion pits have different shapes, sizes, and depths, and their shapes include circular, elliptical, and irregular shapes.

[0080] Taking line 309 as an example, the data in line 309 of elevation_filled is as follows:

[0081] [X0,X1,…,X 448 ,X 449 ,…,X 482 ,…,X 514 ] = [-32767, -32767, ..., 1310.65332, 1310.75366, 1310.79297, 1310.79297, 1310.79297, 1310.79297, 1310.79297, 1310.79626, 1310.80945, 1310.84595, 1310.79089, 1310.72449, 1310.72449, 1310.72449, 1310.72449, 1310.72449, 1310.72449, 1310.72449, 1310.72449 [,1310.72449,1310.72449,1310.72449,1310.72449,1310.72449,1310.72449,1310.72449,1310.72449,1310.72449,1310.72449,1310.72449,1310.72449,1310.72449,1310.72449,1310.70044,1310.62524,…,-32767], and the data in line 309 of elevation:

[0082] [X0,X1,…,X 448 ,X 449 ,…,X 482 ,…,X 514]=[-32767,-32767,…,1310.65332,1310.75366,1310.79297,1310.79053,1310.72205,1310.74756,1310.79004,1310.79626,1310.80945,1310.84595,1310.79089,1310.71838,1310.69678,1310.65295,1310.58215,1310.52832,1310.50134,131 Subtracting the values ​​from [0.46191,1310.43713,1310.38171,1310.31494,1310.16675,1309.99866,1309.82202,1309.63379,1309.45081,1309.67700,1310.04321,1310.30176,1310.48682,1310.61670,1310.69617,1310.72449,1310.70044,1310.62524,…,-32767], we get:

[0083] [X0,X1,…,X 448 ,X 449 ,…,X 482 ,…,X 514 ] = [0,0,…,0,0,0, 0.00244,0.07092, 0.04541,0.00293 ,0,0,0,0, 0.00610,0.02771,0.07153,0.14233,0.19617,0.22314, 0.26257,0.28735,0.34277,0.40955,0.55774,0.72583,0.90247,1.09070,1.27368, 1.04749,0.68127,0.42273,0.23767,0.10779,0.02832 The non-zero regions in this row [,0,0,0,…,0] represent the locations of wind erosion pits. All 1142 rows are calculated using this method, ultimately yielding `elevation_difference`.

[0084] To enable the application of the method of the present invention, the present invention also provides a wind erosion pit identification system based on point cloud data, including a lidar sensor, a memory, a processor, and a display; the lidar sensor is mounted on a navigation device for acquiring point cloud data, for example, the navigation device can be an unmanned or manned aircraft; the memory is used to store computer programs and can also be expanded to store the acquired point cloud data; the processor is used to call and run the computer program stored in the memory to execute the method described in the present invention; the display shows the execution result of the processor. For example, the display and the processor may not be in the same space.

[0085] The above process embodiments and accompanying drawings are descriptions of the functions of the present invention and are not intended to limit the present invention. The variable names, symbols, and software used in the present invention are not considered as limiting conditions of the present invention. Any equivalent transformations made based on the present invention are included within the scope of patent protection of the present invention.

Claims

1. A method for identifying wind erosion pits based on point cloud data, comprising the following steps: Step 1: Obtain point cloud data acquired by the lidar sensor and extract point clouds of ground categories; Step 2: Convert the point cloud of the ground category into raster data, and fill in the missing parts of the raster data to generate a digital elevation model; Step 3: Fill in the concave areas in the digital elevation model; Step 4: By subtracting the digital elevation models before and after filling the depression area, the difference is obtained, and the difference is the identified wind erosion pit. The characteristic feature is that, in step 3, the filling method is as follows: Step 31: Convert the digital elevation model into a 2D array elevation, where each row of the elevation array is a 1D array, and search each row of the elevation array. Step 32: Search point by point from left to right and determine the current point. elevation value Is it less than or equal to the next point? elevation value If so, then Increment the value by 1, continue the judgment, and repeat until... , at this point Let this be labeled, and the corresponding elevation value is... , The value is incremented by 1, where The index of an element in a one-dimensional array; Step 33, Determine Is it true? If it is true, then... The value is incremented by 1, the judgment is continued, and the loop repeats until... At this point, the point label is compared with the current point. The elevation values ​​between them were all modified to That is, filling in the digital elevation model. Increment the value by 1; Step 34: Return to step 32 until the last point of the current line in the elevation is reached; Step 35: Search point by point from right to left and determine the current point. elevation value Is it less than or equal to the next point? elevation value If so, then Decrement the value by 1, continue the judgment, and repeat until... , at this point Let it be labeled as label2, and the corresponding elevation value is , Decrease the value by 1; Step 36, Determine Is it true? If it is true, then... The value is decremented by 1, and the judgment continues, repeating this cycle until... At this point, the point label2 is compared with the current point. The elevation values ​​between them were all modified to That is, filling in the digital elevation model. Decrease the value by 1; Step 37: Return to step 35 until the first point of the current line in the elevation is reached; Step 38: Name the resulting 2D array of the corresponding filled digital elevation model as elevation_filled.

2. The wind erosion pit identification method based on point cloud data according to claim 1, characterized in that, In step 1, a clustering algorithm is used to classify the acquired point cloud in order to extract the point cloud of the ground category.

3. The wind erosion pit identification method based on point cloud data according to claim 2, characterized in that, The acquired point cloud data was clustered according to its color attribute RGB using a clustering algorithm, resulting in three categories. The R and B categories of the clustering results were modified to match the ground category in the LAS data, and the G category of the clustering results was modified to match the vegetation category in the LAS data. The point cloud of the ground category was saved, and the other point clouds were deleted.

4. The wind erosion pit identification method based on point cloud data according to claim 1, characterized in that, The process of converting the point cloud of the ground category into raster data involves using code to read the LAS file, mapping the point cloud into a grid, and describing the grid using a 2D array. If multiple points are mapped to the same grid, the average value of the multiple points is taken as the value of that grid.

5. The wind erosion pit identification method based on point cloud data according to claim 1 or 4, characterized in that, In step 2, the missing parts of the raster data are filled in using the neighborhood averaging method. That is, the average value of the eight non-empty rasters surrounding the missing raster is taken as the value of the missing raster. The missing cells are searched from the top left to the bottom right of the entire raster network to fill in all the missing rasters.

6. The wind erosion pit identification method based on point cloud data according to claim 1, characterized in that, In step 4, the difference between each element in the array elevation_filled and each element in the array elevation is calculated to obtain another 2D array named elevation_difference, which represents the filling of the digital elevation model. The non-zero areas in elevation_difference are the locations of the wind erosion pits.

7. A wind erosion pit identification system based on point cloud data, characterized in that, The device includes a lidar sensor, a memory, a processor, and a display; the lidar sensor is mounted on a navigation device for acquiring location cloud data; the memory stores a computer program, and the processor calls and runs the computer program stored in the memory to perform the method according to any one of claims 1-6. The display shows the execution results of the processor.

8. A computer-readable storage medium, characterized in that, Used to store a computer program that causes a computer to perform the method as described in any one of claims 1-6.

9. A computer program product containing instructions, characterized in that, When the instructions are executed on a computer, the computer performs the method according to any one of claims 1-6.

Images