Self-adaptive ICESat-2 glacier photon denoising method

Abstract

The invention discloses a self-adaptive ICESat-2 glacier photon denoising method, which comprises the following steps of: extracting a potential photon center point coordinate from ATL03 original photon data, fitting a curve by using the center point coordinate and constructing an envelope line; performing coarse denoising on the original photon data through the envelope line to obtain coarse screening photon point data; introducing an LOF strategy to construct a self-adaptive ellipse search window; calculating the slope of the curve to determine the search direction of the self-adaptive elliptical search window; a density clustering method is adopted, a grid adjustment mechanism is introduced to optimize density estimation, and photon space density is calculated; and based on the adaptive elliptical search window, the search direction and the photon space density, performing fine denoising on the coarse screening data to obtain signal photons. The method can accurately identify signal photons and reduce elevation calculation errors, so that the glacier surface height inversion precision is improved, meanwhile, the operation time can be greatly shortened, the algorithm operation efficiency is improved, and the actual requirement for large-scale data rapid processing is met.

Description

Technical Field

[0001] The present invention relates to the technical field of glacier data processing, and in particular to an adaptive ICESat-2 glacier photon denoising method. Background Art

[0002] The ICESat-2 satellite carries a multi-beam micropulse photon counting laser altimeter (ATLAS), which captures high-density photons along its track to measure Earth's surface height. Its ATL03 data product achieves an along-track resolution of 0.7 meters under clear atmospheric conditions, providing crucial data support for glacier surface morphology measurements. However, the ATLAS sensor is prone to recording a large number of noise photons, resulting in a significant amount of noise photons surrounding the captured surface signal photons in the ATL03 data product, especially on complex glacier surfaces. This makes it challenging to directly use this data to invert glacier height.

[0003] Several methods have been studied for retrieving glacier surface elevation from the ICESat-2 ATL03 data product. A representative example is the iterative fitting method (ATL06 algorithm), which is used by NASA to generate the ATL06 product for terrestrial glaciers. This algorithm employs a sliding window strategy, segmenting the along-track photon sequence into approximately 40-meter intervals. Using a local surface fitting model (typically a linear or quadratic surface), it identifies and aggregates ground signal photons, thereby estimating the surface elevation at the center of each segment. This process involves fitting the photon data for each segment to extract valid signal photons from the noisy data, ultimately determining the glacier surface elevation.

[0004] In addition, researchers have proposed a robust density estimation method for retrieving glacier height using ICESat-2 photon counting data. This method uses an elliptical neighborhood that matches the distribution of signal photons on the glacier surface to filter noise. First, the method processes the input glacier ATL03 data product using a multi-scale random sampling consistency fitting method. This method searches for and determines the optimal elliptical neighborhood (with orientation and size parameters) for filtering each photon within the ATL03 data. This step aims to delineate an appropriate region within the raw data, providing a basis for subsequent photon density estimation and signal photon extraction. Secondly, a smoothness-based weighted density and a similarity-based weighted density are calculated. The smoothness-based weighted density considers the smoothness of the spatial distribution of photons, while the similarity-based weighted density focuses on the degree of similarity between photons. These two weighted densities are combined to produce a mixed weighted photon density. This combined calculation process provides a more accurate assessment of the photon distribution density in the data. Finally, based on the mixed weighted photon density, the photons are segmented into adjacent segments, from which candidate signal photons are extracted. Then, an adaptive density threshold is determined based on the data characteristics, and further screening is performed using this threshold to ultimately extract the true signal photons. After these steps, the extracted signal photons are used to invert and output the glacier surface height.

[0005] The ATL06 product, described above, performs poorly in areas with drastic topography variations, such as steep slopes at the ice margin, fracture zones, ridges, and areas with numerous fissures. Its fixed sliding window strategy (with segmentation intervals of approximately 40 meters) and local surface fitting model (linear or quadratic) make it difficult to adapt to the complex and changing glacial terrain. This can easily lead to deviations in elevation estimation or signal omissions, resulting in unstable signal photon recognition and an inability to accurately extract the true glacier surface elevation from the signal photons.

[0006] The robust density estimation method for inverting glacier height using ICESat-2 photon counting data also has significant flaws. Although the multi-scale random sampling consensus (RANSAC) strategy can accurately determine ellipse parameters and achieve good denoising results, it is time-consuming. In large-scale glacier height inversion tasks, due to the need to process massive amounts of data, the long calculation time of the RANSAC strategy significantly prolongs the overall processing cycle, resulting in low efficiency and difficulty in meeting the needs of application scenarios such as quickly obtaining glacier height information, limiting the application and promotion of this method in large-scale practical projects. Summary of the Invention

[0007] The technical problem to be solved by the present invention is to provide an adaptive ICESat-2 glacier photon denoising method that can accurately identify signal photons and reduce elevation calculation errors, thereby improving the accuracy of glacier surface height inversion. At the same time, it can significantly shorten the calculation time, improve the algorithm operation efficiency, and meet the actual needs of rapid processing of large-scale data.

[0008] In order to solve the above technical problems, the present invention provides an adaptive ICESat-2 glacier photon denoising method, which is characterized by comprising the following steps:

[0009] The potential photon center coordinates are extracted from the ATL03 raw photon data, and the extracted center point coordinates are used to fit the curve and construct the envelope. The ATL03 raw photon data is coarsely denoised using the envelope to obtain coarsely screened photon point data.

[0010] Introducing LOF strategy to construct adaptive elliptical search window;

[0011] The slope of the solution curve is used to determine the search direction of the adaptive elliptical search window;

[0012] The density clustering method is used and a grid adjustment mechanism is introduced to optimize the density estimation and calculate the photon spatial density;

[0013] Based on the adaptive elliptical search window, search direction and photon spatial density, the coarse screening data is finely denoised to obtain signal photons.

[0014] Furthermore, the method for extracting the coordinates of the potential photon center point is as follows:

[0015] Divide all photons in the ATL03 original photon data into several grid areas along the track and elevation directions with a set step size. The grid area includes several along-track photon nodes and several height segments divided by each along-track photon node.

[0016] Count the number of photons in each height segment, and take the height segment with the largest number of photons in each photon section along the track as the potential signal photon unit;

[0017] Calculate the mean value of the photon coordinates within the potential signal photon unit and obtain the photon center coordinates.

[0018] Furthermore, the "segment_id" parameter in the ATL03 original photon data is used to block the original photons along the track by section, and then the heights of all sections are divided into blocks in turn to obtain the photon sections and height segments along the track;

[0019] The step length of block processing along the track is 20m, and the step length of block processing along the height is 30m.

[0020] Furthermore, the second highest height segment with the largest number of photons is selected from the same along-track photon section, and the second highest height segment with the largest number of photons is compared with the height segment with the largest number of photons. The height segment that meets the conditions is regarded as a potential signal photon unit. The conditions are as follows:

[0021] N max >N sec ·ρ

[0022] Among them, N max Represents the maximum number of photons in all height segments in a certain photon section along the track, N sec It represents the second largest number of photons among all height segments in a certain along-track photon section, and ρ>1 is an empirical coefficient.

[0023] Furthermore, the smoothing spline fitting method is used for curve fitting.

[0024] Furthermore, the method of constructing an adaptive elliptical search window is as follows:

[0025] First, the signal photons in the coarse-screened photon point data are divided into elevation intervals according to each along-track photon node and the number of photons is counted. The width of each elevation interval is set to 1m.

[0026] The local outlier factor algorithm is used to identify the number of outlier height segments, thereby estimating the vertical distribution range of the actual signal photons in each along-track photon node, that is, the bandwidth;

[0027] The length of the minor axis of the ellipse is calculated by calculating the mean of all bandwidths;

[0028] The axis ratio of the ellipse is set to 2 / 3, and the corresponding major axis length is derived to construct an adaptive ellipse search window;

[0029] Among them, the along-track photon nodes are obtained by dividing the photons in the coarse-screened photon point data into blocks along the track according to the "segment_id" parameter in the coarse-screened photon point data.

[0030] Furthermore, by calculating the first-order derivative of the fitting curve, the slope direction of the surface distribution corresponding to each photon is inverted, and this direction is taken as the main axis direction of the corresponding photon search ellipse, that is, the search direction of the adaptive elliptical search window is determined.

[0031] Furthermore, in the density clustering method, an adaptive MinPts calculation method based on point density estimation is adopted. First, the global average point density ρ is calculated, which is specifically defined as:

[0032]

[0033] Among them, N total is the total number of photon points in the coarse-screened photon point data, Stotal is the mesh area after coarse denoising, which is defined as follows:

[0034] S total =δx·δy

[0035] On this basis, the area of ​​the elliptical search region is set to:

[0036] S ellipse =πab

[0037] The final minimum number of photon points in the neighborhood is:

[0038] MinPts=S ellipse ·ρ

[0039] When the number of photon points N in the adaptive elliptical search window is greater than MinPts, the core point is marked as a signal photon. After traversing each signal photon, the photon points in the cluster are regarded as signal photons, and the other photon points are regarded as noise photons and removed.

[0040] Furthermore, an adaptive grid adjustment mechanism is introduced to optimize the density estimation process. When the initially calculated MinPts<1, the accuracy of local density estimation is gradually improved by recursively reducing the grid height δy until MinPts≥1 is satisfied.

[0041] Beneficial effects of the present invention:

[0042] 1. This invention enables rapid and efficient inversion of glacier surface elevation. While ensuring processing accuracy, it significantly shortens computation time and improves algorithm efficiency, meeting the practical needs of rapid large-scale data processing and providing more timely and accurate data support for glacier monitoring and other related fields.

[0043] 2. The present invention can automatically calculate the corresponding adaptive parameters based on the characteristics of the photon distribution in each ATL03 data, thereby improving adaptability to complex terrain, accurately identifying signal photons, reducing elevation calculation errors, and thus improving the accuracy of glacier surface height inversion.

[0044] 3. The present invention uses the smoothing spline fitting method to perform curve fitting on the glacier surface. Compared with traditional methods such as polynomial regression or piecewise linear fitting, smoothing spline fitting can flexibly adapt to the local changes of the glacier surface from gentle slopes to steep slopes while maintaining the smoothness of the curve.

[0045] 4. In the setting of the adaptive elliptical search window, the bandwidth adaptation strategy of the local outlier factor is adopted, which can effectively reflect the photon distribution characteristics of the ATL03 data under different conditions without relying on the empirical threshold. On the one hand, the recognition based on the LOF algorithm enhances the stability of the obtained signal photon distribution width; on the other hand, by automatically adjusting the elliptical shape according to the photon width, the algorithm has stronger adaptability under different photon distributions.

[0046] 5. In the process of adaptively determining the density threshold, since the photons in the grid may be sparse and the area S total The adaptive grid adjustment mechanism effectively improves the noise suppression capability in low signal-to-noise ratio environments and significantly enhances the stability and reliability of the algorithm. BRIEF DESCRIPTION OF THE DRAWINGS

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

[0048] Figure 2 1 is a schematic diagram comparing the denoising and height inversion results using four methods in Example 2 of the present invention;

[0049] Figure 3 1 is a schematic diagram comparing the height inversion results of the present method and the ATL06 product in Example 2 of the present invention;

[0050] Figure 4 1 is a schematic diagram comparing the denoising and height inversion results using four methods in Example 3 of the present invention;

[0051] Figure 5 1 is a schematic diagram comparing the height inversion results of the present method and the ATL06 product in Example 3 of the present invention;

[0052] Figure 6 1 is a schematic diagram comparing the denoising and height inversion results using four methods in the fourth embodiment of the present invention;

[0053] Figure 7 It is a schematic diagram comparing the height inversion results of this method and the ATL06 product in Example 4 of the present invention. DETAILED DESCRIPTION

[0054] The present invention will be further described below with reference to the accompanying drawings and specific embodiments so that those skilled in the art can better understand the present invention and implement it. However, the embodiments are not intended to limit the present invention.

[0055] Example 1:

[0056] Reference Figure 1As shown, an embodiment of the adaptive ICESat-2 glacier photon denoising method of the present invention effectively solves the problem of low accuracy caused by elevation estimation deviation and signal omission in areas with complex terrain of the ATL06 product. The present invention automatically calculates corresponding adaptive parameters based on the characteristics of the photon distribution in each ATL03 data, improves adaptability to complex terrain, accurately identifies signal photons, reduces elevation calculation errors, and thus improves the accuracy of glacier surface height inversion.

[0057] At the same time, given that the robust density estimation method in the existing technology introduces the multi-scale random sampling consensus (RANSAC) strategy, which leads to the disadvantage of a long algorithm consumption, the present invention introduces the local outlier factor detection (LOF) strategy to determine the size of the adaptive ellipse and determines the direction of the adaptive ellipse by fitting the near-ground curve of the glacier surface. While ensuring processing accuracy, the calculation time can be greatly shortened, the algorithm operation efficiency can be improved, and the actual needs of rapid processing of large-scale data can be met, providing more timely and accurate data support for related fields such as glacier monitoring.

[0058] Specifically, the algorithm in this application processes ATL03 raw photon data using a two-step process: coarse denoising and fine denoising. The coarse denoising step primarily removes long-range noise, reducing the computational effort required in subsequent steps. The fine denoising step uses an adaptive denoising algorithm to refine the noise and extract signal photons from the glacier surface.

[0059] During the coarse denoising process, signal photons in the ATL03 raw photon data, which effectively reflect the glacier surface elevation, are typically concentrated in a relatively narrow band, while noise photons can extend for hundreds or even thousands of meters, widely distributed vertically above and below the echo. Generally speaking, noise photons occupy a larger space than signal photons but have a lower density. Based on the aforementioned characteristics of signal and noise photons, the raw photons are first segmented by segment using the "segment_id" parameter in the ATL03 raw photon data (ATL03 raw photon data assigns a "segment_id" parameter to each photon, dividing photons into segments based on along-track distance, with a fixed length of 20 meters). The height of all segments is then segmented sequentially. Potential signal photon regions are then detected by counting the grid cells with the largest number of photons. The mean photon coordinates within these potential signal photon regions are then calculated to determine the coordinates of their center points. Finally, an approximate curve of the glacier surface is fitted using these extracted center points. Finally, the upper and lower envelopes of the fitting curve are constructed, and based on them, the original photon data is coarsely denoised to obtain coarsely screened photon point data.

[0060] In the process of extracting the center point of the signal photon, all photons are divided into several regions along the track and elevation directions with a certain step size. In order to accurately search for potential signal photon areas, strict constraints are imposed, and each along-track photon node is divided into several height segments with a step size of 30m. The number of photons in all grid cells is counted, and then the height segment with the largest number of points among all height segments corresponding to each along-track photon node is selected as the potential signal photon unit. The coordinate mean of all photons in this height segment is calculated as the signal photon distribution center of this along-track photon node.

[0061] Since the density of signal photons and noise photons in low signal-to-noise ratio areas is not clearly distinguishable, in order to avoid misidentification of signal intervals due to similar photon numbers, it is also necessary to calculate the height segment with the second largest number of photons and compare it with the height segment with the largest number. Finally, the height segment that meets the conditions is regarded as a potential signal photon unit. The conditions are as follows:

[0062] N max >N sec ·ρ

[0063] where N max Represents the maximum number of photons in all height segments of a certain photon node along the track, N sec It represents the number of photons in the second-highest height segment in the above-mentioned along-track photon node. ρ>1 is an empirical coefficient, and the preferred value is 1.5.

[0064] Then, it is necessary to fit the curve function of the glacier surface based on the center point of the signal photon obtained by the above method, and construct the upper and lower envelopes of the fitting curve to perform rough denoising on the original point cloud data. Although the glacier surface has significant volatility and complex terrain undulation characteristics as a whole, at a specific spatial scale, its local surface usually appears as a relatively smooth and approximately linear continuous structure. This local stability provides a theoretical basis and practical feasibility for extracting the near-ground curve of the glacier surface using function fitting methods such as smoothing splines. Therefore, this application selects the smoothing spline fitting method to perform curve fitting on the glacier surface. Compared with traditional methods such as polynomial regression or piecewise linear fitting, smoothing spline fitting can flexibly adapt to the local changes of the natural undulations of the glacier surface from gentle slopes to steep slopes while maintaining the smoothness of the curve. In addition, the spline curve is differentiable, which is convenient for calculating the slope, extracting the envelope or analyzing the trend of surface changes. It is an ideal choice for fitting the glacier surface.

[0065] During fine denoising, noise photons surrounding the signal photons are retained after coarse denoising. Fine denoising requires further identification and extraction of signal photons on the ground. Since photons are typically densely distributed along the track and relatively sparse in the vertical direction, traditional circular neighborhoods are difficult to effectively extract signal photons. Therefore, to increase identification accuracy, an elliptical search window is often used to process the original photons.

[0066] Since the distribution of photon point clouds under different terrains and times has different distribution characteristics, the main parameters such as noise threshold, major and minor axes of the elliptical filter kernel, and direction should be automatically determined according to different distribution characteristics in order to achieve the best extraction effect.

[0067] Specifically, it is necessary to determine the optimal major and minor axes of the elliptical search window:

[0068] The ellipse scale also significantly impacts recognition. If the major axis is too large, the noise photon density in areas with broken terrain may be high, making it difficult to identify. If the minor axis is too large, background noise (such as atmospheric photons) may be introduced, potentially leading to misidentification of near-ground noise photons and affecting the accuracy of subsequent elevation fitting. In areas with steep terrain, a narrow minor axis is also detrimental to capturing continuous surface signals within a range of vertical elevation differences.

[0069] To overcome the above problems, this application proposes a bandwidth adaptive strategy based on the local outlier factor. The specific method is as follows: first, the signal photons after the coarse denoising process are divided into elevation intervals according to each along-track photon node to count the number of photons. In order to perform more precise statistics, the width of each interval is set to 1m; using the characteristic of signal photons gathering in the vertical direction, the local outlier factor (LOF) algorithm is used to identify the number of outlier height segments, thereby estimating the vertical distribution range of the real signal photons in each along-track photon node, that is, the bandwidth. In order to avoid the influence of outliers, the mean of the bandwidth of all along-track photon nodes is calculated as the length of the final ellipse's minor axis, and the ellipse's axis ratio (that is, the ratio of the minor axis to the major axis) is set to 2 / 3, and the corresponding major axis length is derived, thereby constructing an adaptive elliptical search window for subsequent processing.

[0070] This strategy, without relying on empirical thresholds, can effectively reflect the photon distribution characteristics of ATL03 data under different conditions. On the one hand, the recognition based on the LOF algorithm enhances the stability of the acquired signal photon distribution width; on the other hand, by automatically adjusting the elliptical shape according to the photon width, the algorithm has greater adaptability under different photon distributions.

[0071] It is also necessary to determine the optimal search direction:

[0072] Due to the complex undulations of glacial surfaces, the spatial distribution of signal photons often deviates from the horizontal direction as the terrain changes. Traditional fixed-direction elliptical search models have limited recognition accuracy in non-horizontal areas and are unable to effectively extract signal photons from the actual surface structure.

[0073] To enhance the search model's adaptability to terrain changes, this application uses a method for adaptively determining the elliptical search direction based on the derivative of the terrain curve. By calculating the first-order derivative of the fitted surface curve (the fitting method is described in detail in the coarse denoising section), the terrain slope function is inverted. This function is used to calculate the terrain slope direction at each photon's corresponding position. This direction is the main axis direction of the elliptical search window used for the corresponding photon search, allowing the elliptical search window to dynamically align with the changing surface trend.

[0074] Finally, you need to determine the density threshold:

[0075] In density-based clustering methods, the selection of the minimum number of neighborhood points (MinPts) plays a key role in determining whether a point belongs to a high-density area. To adapt to the density distribution of different areas, an adaptive MinPts calculation method based on point density estimation is adopted. That is, a method for adaptively calculating the threshold MinPts based on the average point density in the block area. This method first calculates the global average point density ρ, which is specifically defined as:

[0076]

[0077] Among them, N total is the total number of photon points in the coarse-screened photon point data, S total is the mesh area after coarse denoising, which is defined as follows:

[0078] S total =δx·δy, where δy is the grid height and δx is the grid width.

[0079] On this basis, the area of ​​the elliptical search region is set to:

[0080] S ellipse =πab

[0081] The final minimum number of photon points in the neighborhood is:

[0082] MinPts=S ellipse ·ρ

[0083] When the number of photon points N within the search ellipse is greater than MinPts, the core point is marked as a signal photon. After traversing each photon, the photons in the cluster are considered signal photons, and the others are considered noise photons.

[0084] However, in practical applications, if the photons in a grid are sparse and the area S totalIf the value of ρ is large, the estimated average density ρ is significantly lower, which results in the expected number of points in the ellipse MinPts << (much less than) 1. In this case, even if a region contains only isolated photons, it may be mistakenly identified as a signal point, seriously affecting the accuracy of clustering.

[0085] To address these issues, this application introduces an adaptive grid adjustment mechanism to optimize the density estimation process. When the initially calculated MinPts < 1, the grid height δy is recursively reduced to gradually improve the accuracy of the local density estimate until MinPts ≥ 1 is satisfied. This strategy effectively improves noise suppression in low signal-to-noise ratio environments and significantly enhances the algorithm's stability and reliability.

[0086] Finally, by using the adaptive parameters mentioned above, the real signal photons can be accurately extracted from the ATL03 data, thereby further obtaining accurate glacier elevation information.

[0087] This invention overcomes some of the shortcomings of existing technologies, enabling rapid and efficient inversion of glacier surface elevation. While ensuring processing accuracy, it significantly reduces computation time and improves algorithm efficiency, meeting the practical needs of rapid large-scale data processing and providing more timely and accurate data support for glacier monitoring and other related fields.

[0088] To evaluate the performance of photon denoising and signal identification algorithms, this application constructed an ICESat-2 reference dataset using a manual annotation method based on visual interpretation. This dataset uses confidence labels (0–4) for photon events from the ATL03 product and ground-based photon classification results from the ATL08 product to initially screen potential signal photons. Subsequently, by combining high-resolution Google Earth imagery with DEM data, photon resolution and confirmation were performed above dense ground-based photon areas, eliminating misclassified or missed photons. Finally, the actual types of ground and non-ground signal photons were manually annotated. This reference dataset effectively reflects the actual surface conditions and has important validation significance and application value in areas lacking high-precision airborne lidar data.

[0089] This algorithm uses three common classification metrics based on reference data for evaluation: recall (R), precision (P), and the harmonic mean of the two - F-score (F). These metrics are calculated based on the correspondence between the classification results and the reference true labels and are defined as follows:

[0090]

[0091] Among them, TP (True Positive) represents the number of real signal photons correctly identified as signal photons, FP (False Positive) represents the number of noise photons mistakenly identified as signal photons, and FN (False Negative) represents the number of real signal photons mistakenly identified as noise photons. The recall rate R represents the ratio of the number of successfully detected signal photons to the total number of real signal photons, reflecting the integrity of the algorithm; the precision rate P represents the ratio of the points identified as signal photons to the actual signal photons, measuring the accuracy of the algorithm; the F value is the harmonic average of the recall rate and the precision rate, and it is an important indicator to measure the overall performance of the algorithm by comprehensively considering the balance between the two. The higher the F value, the better the algorithm extracts signal photons. Eight laser beam data were selected, and the results were obtained under the calculation of the three algorithms, as shown in the following table:

[0092]

[0093] As shown in the table above, the average classification accuracy of this algorithm is 95.99% and the average recall rate is 99.23%, which are slightly higher than those of the other two algorithms, proving that the accuracy of this algorithm is slightly improved.

[0094] In terms of running speed, this algorithm efficiently determines the parameters of the ellipse neighborhood through steps such as segmentation along the track (every 20 meters), height grouping, spline curve fitting, and outlier factor detection. Taking a data set of 350,000 photons as an example, the computational workload of this algorithm is approximately 3.25 million operations (O(3.25×10 6 ))). In contrast, the multi-scale RANSAC algorithm needs to perform multiple iterations independently for each photon (assuming 10 RANSAC iterations and 5 scales), with a total computational cost of up to 3.5 billion operations (O(3.5×10 9 ), which is approximately three orders of magnitude higher than our algorithm. This difference in complexity directly leads to a significant difference in runtime. Our algorithm runs approximately 1,000 times faster than the multi-scale RANSAC algorithm (on a single core), significantly improving processing efficiency and making it particularly suitable for real-time applications or processing large datasets.

[0095] Example 2:

[0096] Reference Figure 2 As shown in the figure, based on the dataset under weak beam conditions, it can be seen from the figure that the noise distribution is uneven. Four methods are used to process the dataset respectively, among which: Visual Results is the result of processing the dataset by the manual annotation method based on visual interpretation, Our Method is the result of processing the dataset by the method of this application, ATL03 is the result of processing the dataset by photon confidence labeling, and DBSCAN is the result of processing the dataset by the density-based clustering method;

[0097] Since the weak light beam signal is not strong, the density distribution difference between signal photons and noise photons is not obvious, making it difficult to effectively distinguish between the two. Judging from the confidence labeling results of ATL03 and the DBSCAN clustering results, there is a more obvious misjudgment of noise photons when the noise photon density difference is large. Although the algorithm proposed in this application identifies relatively sparse signal photons in this area, and there are cases where signal photons are misjudged as noise in some areas, it is still better than other methods overall and can more accurately extract the height information of the glacier surface. This advantage provides a higher accuracy guarantee for the subsequent glacier height inversion.

[0098] Specifically, from Figure 2 It can be clearly seen in areas (a1) and (a2) that when the noise photon density varies greatly, the confidence labeling results of ATL03 and the results of the DBSCAN algorithm cannot achieve satisfactory results. The main reason is that the DBSCAN algorithm relies on global fixed parameters and cannot dynamically adjust the parameter settings according to different density areas.

[0099] Again, the method of this application is compared with the results of processing the data set using the ATL06 algorithm. Since the ATL06 product does not have the along-track distance parameter, the latitude is selected as the horizontal axis to draw a comparison chart of the height inversion results, as shown in the figure. Figure 3 As shown, from Figure 3 The height inversion results for regions (b1) and (b2) show that the official ATL06 land ice height product cannot provide accurate results in some complex situations. However, this method can more accurately extract glacier surface height information, not only providing reliable glacier surface elevations but also providing valuable applications for subsequent glacier dynamic monitoring and change analysis.

[0100] Example 3:

[0101] Reference Figure 4 As shown, based on the data set under strong light beam conditions, it can be seen in the figure that the noise distribution in the data set is sparse, and four methods are used for processing respectively. The number of signal photons is much larger than the noise photons, and there is a significant difference in the spatial density between the two. Under this type of data, all three methods can achieve good denoising effects, but there are still differences in performance. It is worth noting that the adaptive threshold determination method used in this application has certain limitations in this scenario. Because this method adaptively determines the division threshold through density differences, the original photons are divided into two categories: signals and noise, and when the number of noise photons is much lower than that of signal photons, some signal photons with lower density are inevitably misidentified as noise photons, resulting in a slight loss of edge information. However, in comparison with the ATL06 land ice height product chart Figure 5It can be seen that the glacier height results obtained by the two methods are consistent, and there is no large height error.

[0102] In comparison, the DBSCAN algorithm performs best for this type of data, preserving the signal photons more completely. However, this algorithm also has significant shortcomings. Its denoising performance is highly dependent on parameter settings, requiring repeated adjustments for different data types to achieve optimal results, which reduces the algorithm's versatility and efficiency in practical applications.

[0103] Example 4:

[0104] Reference Figure 6 As shown in the figure, based on the dataset under strong beam conditions, four methods were used to process them respectively. The figure shows the denoising effect of weak beams in this dataset and a schematic diagram of height inversion. In the data type with uniform noise distribution, the weak beam band has a significantly reduced signal-to-noise ratio due to the small number of signal photons. Under this condition, although the three methods can still extract relatively complete surface information as a whole, there are obvious differences in the extraction results. Specifically, the confidence method based on the ATL03 marker performed the worst, with multiple breaks and lack of continuity in the extraction results. In contrast, the method used in this study performed similarly to the DBSCAN algorithm under weak beam conditions. Both can effectively identify the glacier surface contour and overcome the interference caused by signal sparseness to a certain extent, resulting in relatively ideal extraction results.

[0105] Reference Figure 7 As shown in the figure, compared with the ATL06 land ice height product, the glacier height results obtained by the two are consistent, and there is no large height error.

[0106] In summary, based on Examples 2 to 4, this algorithm was validated on datasets with three different noise photon distribution types. The experimental results show that this algorithm performs best in the case of uneven noise distribution, achieving significant improvements over the other two algorithms. In the cases of uniform and sparse noise distribution, the three algorithms perform similarly. Further comparisons of the inverted altitude results with NASA's ATL06 land ice height product revealed broadly consistent results, with this algorithm outperforming it in some complex scenarios.

[0107] The above embodiments are only preferred embodiments for fully illustrating the present invention, and the protection scope of the present invention is not limited thereto. Any equivalent substitution or modification made by those skilled in the art based on the present invention is within the protection scope of the present invention.

Claims

1. An adaptive ICESat-2 glacier photon denoising method, characterized by: The following steps are involved: Extract the potential photon center coordinates from the ATL03 raw photon data, fit the curve using the extracted center coordinates, and construct the envelope. The ATL03 raw photon data is coarsely denoised using the envelope curve to obtain coarsely screened photon point data; Introducing LOF strategy to construct adaptive elliptical search window; The slope of the solution curve is used to determine the search direction of the adaptive elliptical search window; The density clustering method is used and a grid adjustment mechanism is introduced to optimize the density estimation and calculate the photon spatial density; Based on the adaptive elliptical search window, search direction and photon spatial density, the coarse screening data is finely denoised to obtain signal photons.

2. The adaptive ICESat-2 glacier photon denoising method according to claim 1, wherein: The method to extract the coordinates of the potential photon center point is as follows: Divide all photons in the ATL03 original photon data into several grid areas along the track and elevation directions with a set step size. The grid area includes several along-track photon nodes and several height segments divided by each along-track photon node. Count the number of photons in each height segment, and take the height segment with the largest number of photons in each photon section along the track as the potential signal photon unit; Calculate the mean value of the photon coordinates within the potential signal photon unit and obtain the photon center coordinates.

3. The adaptive ICESat-2 glacier photon denoising method according to claim 2, wherein: The "segment_id" parameter in the ATL03 original photon data is used to block the original photons along the track by section, and then the heights of all sections are divided into blocks in turn to obtain the photon sections and height segments along the track; The step length of block processing along the track is 20m, and the step length of block processing along the height is 30m.

4. The adaptive ICESat-2 glacier photon denoising method according to claim 2, wherein: In the same along-track photon section, the second highest height segment with the largest number of photons is selected and compared with the highest height segment with the highest number of photons. The height segment that meets the conditions is regarded as a potential signal photon unit. The conditions are as follows: N max >N sec ·r Among them, N max Represents the maximum number of photons in all height segments in a certain photon section along the track, N sec It represents the second largest number of photons among all height segments in a certain along-track photon section, and ρ>1 is an empirical coefficient.

5. The adaptive ICESat-2 glacier photon denoising method according to claim 1, wherein: The smoothing spline fitting method is used for curve fitting.

6. The adaptive ICESat-2 glacier photon denoising method according to claim 2, wherein: The method to construct an adaptive elliptical search window is as follows: First, the signal photons in the coarse-screened photon point data are divided into elevation intervals according to each along-track photon node and the number of photons is counted. The width of each elevation interval is set to 1m. The local outlier factor algorithm is used to identify the number of outlier height segments, thereby estimating the vertical distribution range of the actual signal photons in each along-track photon node, that is, the bandwidth; The length of the minor axis of the ellipse is calculated by calculating the mean of all bandwidths; The axis ratio of the ellipse is set to 2 / 3, and the corresponding major axis length is derived to construct an adaptive ellipse search window; Among them, the along-track photon nodes are obtained by dividing the photons in the coarse-screened photon point data into blocks along the track according to the "segment_id" parameter in the coarse-screened photon point data.

7. The adaptive ICESat-2 glacier photon denoising method according to claim 1, wherein: By calculating the first-order derivative of the fitting curve, the slope direction of the surface distribution corresponding to each photon is inverted, and this direction is the main axis direction of the corresponding photon search ellipse, that is, the search direction of the adaptive elliptical search window is determined.

8. The adaptive ICESat-2 glacier photon denoising method according to claim 1, wherein: In the density clustering method, an adaptive MinPts calculation method based on point density estimation is adopted. First, the global average point density ρ is calculated, which is specifically defined as: Among them, N total is the total number of photon points in the coarse-screened photon point data, S total is the mesh area after coarse denoising, which is defined as follows: S total =δx·δy On this basis, the area of ​​the elliptical search region is set to: S ellipse =πab The final minimum number of photon points in the neighborhood is: MinPts=S ellipse ·r When the number of photon points N in the adaptive elliptical search window is greater than MinPts, the core point is marked as a signal photon. After traversing each signal photon, the photon points in the cluster are regarded as signal photons, and the other photon points are regarded as noise photons and removed.

9. The adaptive ICESat-2 glacier photon denoising method according to claim 8, wherein: An adaptive grid adjustment mechanism is introduced to optimize the density estimation process. When the initially calculated MinPts<1, the accuracy of local density estimation is gradually improved by recursively reducing the grid height δy until MinPts≥1 is satisfied.

Images