Rough mountain area satellite-borne photon counting laser radar point cloud data denoising method

Abstract

The invention discloses a rugged mountain area satellite-borne photon counting laser radar point cloud data denoising method. Firstly, coarse denoising is carried out, point cloud data is divided into a plurality of fragments, dense center points of each fragment are found, and an elevation confidence interval is set according to the dense center points to remove noise. Then fine denoising is carried out, point cloud fragments are divided again with dense center points as the reference, the slope angle of each fragment is calculated, and proper fragment combination is carried out; the Euclidean distance between two point clouds is decomposed into distances (DPSL and DVSL) parallel and perpendicular to the gradient direction, local distance statistics is performed by using the DVSL to distinguish signal photons and noise photons, a denoising threshold is searched based on a DVSL frequency histogram to remove noise, and elevation difference denoising is performed by using a ground point cloud fitting ground curve to further remove noise. The method is suitable for rugged bare land and vegetation areas, has a good denoising effect, provides a basis for point cloud data classification, earth surface extraction, canopy height estimation and the like, and has a wide application prospect.

Description

Technical Field

[0001] The present invention relates to the field of satellite-borne laser radar technology, and in particular to a method for denoising point cloud data of satellite-borne photon counting laser radar in a rugged mountainous area. Background Art

[0002] In recent years, with the rapid development of remote sensing technology, space-borne photon counting laser radar has become an important tool for earth observation. It emits tens of thousands of laser pulses per second and measures the time it takes to reflect back from the earth's surface, thereby accurately calculating the change in surface height. It is an active remote sensing technology that can directly, quickly and accurately obtain three-dimensional spatial information of the detection target. It has been widely used in topographic mapping, forestry surveys, power inspections, and the establishment of three-dimensional urban models. Space-borne photon counting laser radar has unique advantages in obtaining three-dimensional surface information in large areas due to its high orbit and large observation range.

[0003] Compared with traditional lidar, the laser pulse signals emitted and detected by photon counting lidar are weaker, so it is difficult to effectively distinguish between the pulses reflected from the surface of the target object, the light scattered by the atmosphere, the solar radiation, and the noise of the instrument itself. The system generates point clouds for all recorded photon events, resulting in a large number of noisy point clouds in the data. Therefore, the key to data processing and quantitative application is to effectively remove the noise point cloud and extract the real signal point cloud.

[0004] The existing technology for denoising current spaceborne photon counting lidar point cloud data still has the following shortcomings: (1) The noise density in spaceborne photon counting lidar point cloud data is not uniform and varies with time and space. The existing algorithms are prone to misjudgment in high noise density areas, mistaking noise for signal, or cannot completely remove noise in low noise density areas; (2) When dealing with different terrain boundaries, such as the junction of land and water, and the junction of vegetation and bare land, the existing algorithms may be affected by boundary effects. For example, when the local window crosses different terrain types, the statistical characteristics of the signal and noise may change significantly, making it difficult for the algorithm to correctly distinguish between the signal and noise. Summary of the invention

[0005] To solve the above problems, the present invention proposes a method for denoising point cloud data of spaceborne photon counting lidar in rugged mountainous areas;

[0006] To implement the above technical solution, the specific steps are as follows:

[0007] A method for denoising point cloud data of a spaceborne photon counting lidar in a rugged mountainous area comprises the following steps:

[0008] S1, obtain point cloud data of rugged mountainous areas through open source data;

[0009] The open source data is: select the ICESat-2 / ATLAS Level 2 data ATL03 from the official website of the National Snow and Ice Data Center of the United States; the required data are the elevation data (h_ph) of the point cloud and the along-track distance (along_track) of the point cloud, and establish a coordinate system with the along-track distance as the horizontal coordinate and the elevation data as the vertical coordinate, corresponding to the point cloud coordinate information (x,h);

[0010] S2. Perform the first denoising based on the dense center points of the point cloud data. The specific steps are as follows:

[0011] S2.1, segment the point cloud data at a fixed distance along the track and record the coordinates of the dense center point in each segment into data set A;

[0012] For large-scale planting areas in plains, the maximum value can be set to 200m due to its gentle terrain and small fluctuations; for forests in southwestern mountainous areas, the minimum value can be set to 50-100m due to its rugged terrain and large changes; because the mountain vegetation coverage area is relatively rugged, the fixed along-track distance selected by the present invention is set to 100m;

[0013] Furthermore, the fixed along-track distance is used as the horizontal coordinate, and the point cloud data is divided along the vertical direction of the horizontal coordinate. The area between two adjacent along-track distances is a segment, which is divided into M segments;

[0014] Dense center points are calculated as follows:

[0015]

[0016] In the formula, x i 、h i and x j 、h j They are point clouds p i and point cloud p j The horizontal and vertical coordinates of dist(p i ,p j ) is the point cloud p i With point cloud j The Euclidean distance between them; R is the radius of the neighborhood, num(p i ) is the point cloud p i The number of neighboring point clouds in the neighborhood, num max is num(p i ), i.e., the dense center point in the point cloud segment, and N is the number of point clouds in the point cloud data;

[0017] S2.2, using the elevation data of the dense center point in each segment as a reference to set the elevation confidence interval, retaining the point cloud markers within the elevation confidence interval, and removing the point cloud markers outside the elevation confidence interval;

[0018] The confidence interval for elevation is [h_ph-200,h_ph+200];

[0019] S3. Perform a second denoising on the point cloud data obtained after the first denoising to obtain the final denoising result. The steps are as follows:

[0020] S3.1, re-segment the point cloud according to the along-track distance of the dense center points in data set A;

[0021] Taking the distance along the track of each dense center point as the horizontal coordinate, the area between two consecutive dense center points in the vertical direction of the horizontal coordinate constitutes a segment, and the point cloud data is further divided into M segments;

[0022] S3.2, calculating the slope angle according to the horizontal and vertical coordinates of the two dense center points, and the obtained slope angle is the slope angle of this point cloud segment; calculating the slope angles of all segments in turn and merging the segments according to the slope angle values;

[0023] The slope angle calculation formula is as follows:

[0024]

[0025] In the formula, x q 、h q and x q+1 、h q+1 They are the dense center points p q and the dense center point p q+1 The horizontal and vertical coordinates of , M is the number of point cloud fragments;

[0026] If the slope angles of two consecutive point cloud segments are the same and the value difference is less than 3°, they will be merged into one point cloud segment;

[0027] S3.3. Calculate the slope angle θ after merging the segments w ;

[0028] The slope angle θ of the merged point cloud segment w is the mean slope angle of the two segments before merging;

[0029] S3.4, setting the number K of neighboring point clouds of all point clouds in the merged segment;

[0030] The K value is 20-30 for plain areas; the K value is 40-60 for mountainous forest areas, and the value of the present invention is 50;

[0031] S3.5, decomposing the Euclidean distance between adjacent point clouds into the distance parallel to the terrain slope and the distance perpendicular to the terrain slope, and then performing local distance statistics;

[0032] The expression of distance decomposition is as follows:

[0033]

[0034] i∈(1,N)

[0035] j∈(1,N)

[0036] i≠j

[0037] In the formula, x i 、h i and x j 、h j are the point clouds p in the fragment respectively i and point cloud p j The horizontal and vertical coordinates, dist(p i ,p j ) is the point cloud p i With point cloud j The Euclidean distance between i ,p j ) and DVSL(p i ,p j ) are dist(p i ,p j ) is the distance decomposed in the direction parallel to the terrain slope and in the direction perpendicular to the terrain slope, θ w is the slope angle of the corresponding point cloud segment, and N is the number of point clouds in the point cloud data;

[0038] Local distance statistics sum DVSL The calculation formula is:

[0039]

[0040] In the formula, DVSL(p i ,p j ) is the point cloud p i and point cloud p i The Euclidean distance between them is decomposed in the direction perpendicular to the terrain slope;

[0041] S3.6, establishing a length frequency histogram according to the distances of all point clouds in the segment in the direction perpendicular to the slope angle, forming a Gaussian fitting curve by fitting a Gaussian function, and performing a second denoising;

[0042] The horizontal coordinate corresponding to the peak position of the ground peak Gaussian fitting curve is set as the ground point cloud screening threshold;

[0043] The horizontal coordinate corresponding to the intersection of the noise peak Gaussian fitting curve and the nearest signal peak fitting curve is set as the denoising threshold;

[0044] The sum of each point cloudDVSL Compare with the denoising threshold. If it is greater than the threshold, it is noise and will be removed. If it is less than the threshold, it is a signal and will be retained.

[0045] The signal point cloud obtained by the denoising threshold is compared with the fitted ground line in elevation. If it is lower than the fitted ground line by more than 3m, it is considered to be a misjudged noise point cloud and is removed. The remaining point clouds are regarded as signals and retained.

[0046] The fitting ground line selects the ground point cloud according to the ground point cloud screening threshold, repeats the process of finding dense center points, finds the dense center points in the ground point cloud, and performs cubic spline interpolation with the dense center points in the ground point cloud as data;

[0047] Beneficial effects of the present invention

[0048] Compared with the traditional local distance statistics denoising algorithm, the present invention uses two-stage denoising. The process of first performing coarse denoising and then fine denoising reduces the amount of computation for subsequent fine denoising and improves the efficiency of the algorithm. At the same time, the point cloud is segmented based on dense center points and the segments are merged according to the terrain slope, which reduces the misjudgment problem caused by unstable noise density and improves the adaptability to different application scenarios. The final point cloud segmentation has a certain degree of adaptability to the slope, density and surface type, which reduces the influence of boundary effects on the algorithm results. Finally, due to the use of DVSL for local distance statistics, the influence of long-distance point clouds distributed along the slope direction on the algorithm accuracy is reduced, the denoising accuracy is improved and the false deletion rate is reduced. BRIEF DESCRIPTION OF THE DRAWINGS

[0049] Figure 1 is a flow chart of the present invention;

[0050] Figure 2 It is a schematic diagram of the overall method of the present invention;

[0051] Figure 3 The original ATL03 data point cloud of the present invention;

[0052] Figure 4 It is a rough denoising flow chart of the dense center point elevation search method of the present invention;

[0053] Figure 5 It is a schematic diagram of the division of the fixed along-track distance of the present invention;

[0054] Figure 6 A rough denoising result diagram of the dense center point elevation search method of the present invention;

[0055] Figure 7 It is a flow chart of the local distance statistics denoising algorithm based on terrain slope direction decomposition of the present invention;

[0056] Figure 8 A point cloud segmentation graph of the present invention that is segmented according to the distance along the track of the dense center points;

[0057] Fig. 9 The point cloud segmentation map after merging according to the slope angle of the present invention;

[0058] Fig.10 It is the Gaussian function fitting result diagram of the present invention;

[0059] Fig.11 This is a denoising result diagram of the elevation difference of the fitted ground line of the present invention;

[0060] Fig.12 This is the final result diagram of the present invention. DETAILED DESCRIPTION

[0061] The specific implementation method of the present invention is fully and clearly described below in conjunction with the accompanying drawings of the present invention. The examples described below are only some examples rather than all examples.

[0062] like Figure 1 , Figure 2 As shown, a method for denoising point cloud data of spaceborne photon counting lidar in a rugged mountainous area comprises the following steps:

[0063] S1, obtain point cloud data of rugged mountainous areas through open source data;

[0064] like Figure 3 As shown, in this embodiment, the ICESat-2 / ATLAS Level 2 data ATL03 of Shangri-La City is selected through the official website of the National Snow and Ice Data Center of the United States; the required data are the elevation data (h_ph) of the point cloud and the along-track distance (along_track) of the point cloud, and a coordinate system is established with the along-track distance as the horizontal coordinate and the elevation data as the vertical coordinate, corresponding to the point cloud coordinate information (x, h);

[0065] S2, such as Figure 4 As shown, the first denoising is performed based on the dense center points of the point cloud data; the specific steps are as follows:

[0066] S2.1, segment the point cloud data at a fixed distance along the track and record the coordinates of the dense center point in each segment into data set A;

[0067] For large-scale planting areas in plains, the maximum value of 200m can be set due to the gentle terrain and small fluctuations; for forests in southwestern mountainous areas, the minimum value of 50-100m can be set due to the rugged terrain and large changes; because the mountain vegetation coverage area is relatively rugged, the fixed along-track distance selected in this embodiment is set to 100m;

[0068] Specifically, the fixed along-track distance is used as the horizontal coordinate, and the point cloud data is divided along the vertical direction of the horizontal coordinate. The area between two adjacent along-track distances is a segment. The point cloud data is divided into M segments. In this embodiment, M is 27. Figure 5 As shown;

[0069] Furthermore, the dense center point of each fragment is calculated as follows:

[0070]

[0071] i∈(1,N)

[0072] j∈(1,N)

[0073] i≠j

[0074] In the formula, x i 、h i and x j 、h j They are point clouds p i and point cloud p j The horizontal and vertical coordinates of dist(p i ,p j ) is the point cloud p i With point cloud j The Euclidean distance between them; R is the radius of the neighborhood range, which is set to 10 in this embodiment; num(p i ) is the point cloud p i The number of neighboring point clouds in the neighborhood, num max is num(p i ), i.e., the dense center point in the point cloud segment, N is the number of point clouds in the point cloud data, and in this embodiment, N is 23353;

[0075] S2.2, using the elevation data of the dense center point in each segment as a reference to set the elevation confidence interval, retaining the point cloud markers within the elevation confidence interval, and removing the point cloud markers outside the elevation confidence interval;

[0076] In this embodiment, considering the steep slopes and possible trees in the rugged mountains, the confidence interval is set to [h_ph-200, h_ph+200]; the rough denoising result is as follows Figure 6 As shown;

[0077] S3, such as Figure 7 As shown, according to the point cloud data obtained after the first denoising, the second denoising is performed to obtain the final denoising result;

[0078] S3.1, re-segment the point cloud according to the along-track distance of the dense center points in data set A;

[0079] like Figure 8 As shown, the distance along the track of each dense center point is taken as the horizontal coordinate, and the area between two consecutive dense center points in the vertical direction of the horizontal coordinate constitutes a segment, and the point cloud data is further divided into M segments. In this embodiment, M is 27;

[0080] S3.2, calculating the slope angle according to the horizontal and vertical coordinates of the two dense center points, and the obtained slope angle is the slope angle of this point cloud segment; calculating the slope angles of all segments in turn and merging the segments according to the slope angle values;

[0081] The slope angle calculation formula is as follows:

[0082]

[0083] In the formula, x q 、h q and x q+1 、h q+1 They are the dense center points p q and the dense center point p q+1 , M is the number of point cloud fragments. In this embodiment, M is 27;

[0084] If the slope angle is too large, the segments that should not be merged will be merged together, reducing the accuracy of subsequent processing; if the slope angle is too small, too many segments will be divided, reducing the efficiency of the algorithm;

[0085] This embodiment uses empirical values ​​to determine the slope angles of two consecutive point cloud segments. If the slope angles of two consecutive point cloud segments are the same and the difference between the slope angles is less than 3°, the two segments are merged into one point cloud segment. The merged point cloud segment is as follows: Fig. 9 As shown;

[0086] S3.3. Calculate the slope angle θ after merging the segments w ;

[0087] Specifically, the slope angle θ of the merged point cloud fragment is w is the mean slope angle of the two segments before merging;

[0088] S3.4, setting the number of neighboring point clouds of all point clouds in the merged segment;

[0089] Specifically, the trial and error method is used to determine the number K of neighboring point clouds involved in the calculation of local distance statistics. K is essentially a smoothing factor, which is used to smooth the spatial heterogeneity within the same category of point clouds (non-noise points);

[0090] If the K value is too small, the smoothing effect is not obvious; if the K value is too large, the spatial distribution difference between the noisy point cloud and the non-noise point cloud is smoothed, making point cloud denoising difficult;

[0091] The K value is 20-30 for plain areas; the K value is 40-60 for mountainous forest areas, and the value in this embodiment is 50;

[0092] S3.5, decomposing the Euclidean distance between adjacent point clouds into: the distance parallel to the terrain slope and the distance perpendicular to the terrain slope;

[0093] Specifically, for each point cloud in the segment, the KNN algorithm is used to find the K nearest point clouds, the Euclidean distance between each point cloud and its nearest neighboring point cloud is calculated, and the Euclidean distance between the two point clouds is decomposed into: the distance parallel to the slope angle (DPSL) and the distance perpendicular to the slope angle (DVSL);

[0094] Use DVSL to replace the Euclidean distance between neighboring point clouds to perform local distance statistics;

[0095] The expression of distance decomposition is as follows:

[0096]

[0097] i∈(1,N)

[0098] j∈(1,N)

[0099] i≠j

[0100] In the formula, x i 、h i and x j 、h j are the point clouds p in the fragment respectively i and point cloud p j The coordinates of dist(p i ,p j ) is the point cloud p i With point cloud j The Euclidean distance between i ,p j ) and DVSL(p i ,p j ) are dist(p i ,p j ) is the distance decomposed in the direction parallel to the terrain slope and in the direction perpendicular to the terrain slope, θ w is the slope angle of the corresponding point cloud segment, N is the number of point clouds in the point cloud data, and in this embodiment, N is 23353;

[0101] Local distance statistics sum DVSL The calculation formula is:

[0102]

[0103] In the formula, DVSL(p i ,p j ) is the point cloud p i and point cloud p i The Euclidean distance between them is decomposed in the direction perpendicular to the terrain slope;

[0104] S3.6, establishing a length frequency histogram according to the distances of all point clouds in the segment in the direction perpendicular to the slope angle, forming a Gaussian fitting curve by fitting a Gaussian function, and performing a second denoising;

[0105] Specifically, for the obtained length frequency histogram, a Gaussian fitting curve is formed by fitting a Gaussian function, such as Fig.10 As shown;

[0106] The number of Gaussian fitting curves varies according to the different types of regional surface. If the surface is bare land, there are two Gaussian fitting curves, namely the ground peak Gaussian fitting curve and the noise peak Gaussian fitting curve; if the surface is forest, there are three Gaussian fitting curves, namely the ground peak Gaussian fitting curve, the canopy peak Gaussian fitting curve and the noise peak Gaussian fitting curve. The ground peak Gaussian fitting curve and the canopy peak Gaussian fitting curve are both signal peak Gaussian fitting curves.

[0107] The horizontal coordinate corresponding to the peak position of the ground peak Gaussian fitting curve is the ground point cloud screening threshold. In this embodiment, the ground point cloud screening threshold is 37.7784;

[0108] The horizontal coordinate corresponding to the intersection of the noise peak Gaussian fitting curve and the canopy peak fitting curve is the denoising threshold. In this embodiment, the denoising threshold is 139.57;

[0109] The sum of each point cloud DVSL Compare with the denoising threshold. If it is greater than the threshold, it is noise and will be removed. If it is less than the threshold, it is a signal and will be retained.

[0110] The ground point cloud is screened out according to the ground point cloud screening threshold, and the process of finding dense center points is repeated to find the dense center points in the ground point cloud. The dense center points in the ground point cloud are used as data for cubic spline interpolation to obtain the fitting ground line. The results are as follows: Fig.11 As shown;

[0111] The signal point cloud obtained by the denoising threshold is compared with the fitting ground line in elevation. If it is lower than the fitting ground line by more than 3m, it is considered to be a misjudged noise photon and is removed. The remaining photons are regarded as signals and retained. Fig.12 This is the final denoising result.

Claims

1. A method for denoising point cloud data of spaceborne photon counting lidar in rugged mountainous areas, characterized in that: The following steps are involved: S1, obtain point cloud data of rugged mountainous areas through open source data; S2, first denoising based on the dense center points of the point cloud data; S3, performing a second denoising on the point cloud data obtained after the first denoising to obtain the final denoising result; The second denoising is performed by establishing a length frequency histogram based on the distance in the direction perpendicular to the slope angle and forming a Gaussian fitting curve by fitting a Gaussian function.

2. The method for denoising point cloud data of spaceborne photon counting lidar in rugged mountainous areas according to claim 1 is characterized in that: The steps of performing the first denoising based on the dense center point of the point cloud data are as follows: S2.1, segment the point cloud data at a fixed distance along the track and record the coordinates of the dense center point in each segment into data set A; Dense center points are calculated as follows: In the formula, x i 、h i and x j 、h j They are point clouds p i and point cloud p j The horizontal and vertical coordinates of dist(p i ,p j ) is the point cloud p i With point cloud p j The Euclidean distance between them; R is the radius of the neighborhood, num(p i ) is the point cloud p i The number of neighboring point clouds in the neighborhood, num max is num(p i ), i.e., the dense center point in the point cloud segment, and N is the number of point clouds in the point cloud data; S2.

2. Set the elevation confidence interval based on the elevation data of the dense center point in each segment, retain the point cloud markers within the elevation confidence interval, and remove the point cloud markers outside the elevation confidence interval.

3. The method for denoising point cloud data of spaceborne photon counting lidar in rugged mountainous areas according to claim 1 is characterized in that: The specific steps of the second denoising are as follows: S3.1, re-segment the point cloud according to the along-track distance of the dense center points in data set A; S3.2, calculating the slope angle according to the horizontal and vertical coordinates of the two dense center points, and the obtained slope angle is the slope angle of this point cloud segment; calculating the slope angles of all segments in turn and merging the segments according to the slope angle values; The slope angle calculation formula is as follows: In the formula, x q 、h q and x q+1 、h q+1 They are the dense center points p q and the dense center point p q+1 The horizontal and vertical coordinates of , M is the number of point cloud fragments; If the slope angles of two consecutive point cloud segments are the same and the value difference is less than 3°, they will be merged into one point cloud segment; S3.

3. Calculate the slope angle θ after merging the segments w ; The slope angle θ of the merged point cloud segment w is the mean slope angle of the two segments before merging; S3.4, setting the number K of neighboring point clouds of all point clouds in the merged segment; S3.5, decomposing the Euclidean distance between adjacent point clouds into the distance parallel to the terrain slope and the distance perpendicular to the terrain slope, and then performing local distance statistics; The expression of distance decomposition is as follows: In the formula, x i 、h i and x j 、h j They are point clouds p i and point cloud p j The horizontal and vertical coordinates, dist(p i ,p j ) is the point cloud p i With point cloud p j The Euclidean distance between i ,p j ) and DVSL(p i ,p j ) are dist(p i ,p j ) is the distance decomposed in the direction parallel to the terrain slope and in the direction perpendicular to the terrain slope, θ w is the slope angle of the corresponding point cloud segment, and N is the number of point clouds in the point cloud data; The calculation formula for local distance statistics is: In the formula, DVSL(p i ,p j ) is the point cloud p i and point cloud p i The Euclidean distance between them is decomposed in the direction perpendicular to the terrain slope; S3.6, establishing a length frequency histogram according to the distances of all point clouds in the segment in the direction perpendicular to the slope angle, forming a Gaussian fitting curve by fitting a Gaussian function, and performing a second denoising; The Gaussian fitting curves include: a ground peak Gaussian fitting curve, a noise peak Gaussian fitting curve, and a canopy peak Gaussian fitting curve; The horizontal coordinate corresponding to the peak position of the ground peak Gaussian fitting curve is the ground point cloud screening threshold. The local distance statistics of each point cloud are compared with the denoising threshold. If it is greater than the threshold, it is noise and is removed. If it is less than the threshold, it is a signal and is retained. The signal point cloud obtained by the denoising threshold is compared with the fitting ground line in elevation. If it is lower than the fitting ground line by more than 3m, it is considered to be a misjudged noise photon and is removed. The remaining photons are regarded as signals and retained. The fitting ground line is obtained by filtering out the ground point cloud according to the ground point cloud screening threshold, repeating the process of finding dense center points, finding the dense center points in the ground point cloud, and performing cubic spline interpolation with the dense center points in the ground point cloud as data.

Images