Fine and reliable surface and canopy top extraction method and system applied to icesat-2 data

Abstract

The application discloses a fine and high-reliability ground surface and canopy top extraction method and system applied to ICESat-2 data, first, an adaptive threshold ellipsoid filtering method is used to roughly distinguish noise photon groups and signal photon groups, then a dense area noise removal method is used to further remove noise photon groups in the signal photon groups, then an improved local outlier factor (LOF) algorithm is used to remove discrete single noise photons adjacent to the signal photons to obtain a fine denoising result, and finally, an improved local direction center algorithm is used in combination with a uniaxial inverse distance weight statistical method to finely distinguish the ground surface and the canopy top. The application provides a new idea for extracting ground surface and canopy top height information from ICESat-2 data, and can effectively assist carbon sink calculation.

Description

TECHNICAL FIELD

[0001] The application belongs to the technical field of laser radar, and particularly relates to a fine and high-reliability ground surface and canopy top extraction method applied to ICESat-2 data. BACKGROUND

[0002] As an active remote sensing detection technology, spaceborne laser (radar) altimetry has the ability to penetrate observation operations, and is widely used in long-range ranging, environmental monitoring and three-dimensional imaging fields due to its high measurement accuracy, high spatial and temporal resolution and high vertical resolution. There are three kinds of spaceborne laser altimeters, namely discrete laser radar, full waveform laser radar and photon counting laser radar. Among them, the first application of photon counting laser radar is on the Ice, Cloud and Land Elevation Satellite-2 (ICESat-2) launched by the National Aeronautics and Space Administration (NASA) of the United States in 2018. The Advanced Topographic Laser Altimeter System (ATLAS) carried by ICESat-2 adopts a high-sensitivity single-photon detector and designs 3 groups of 6 beams, each group consisting of a strong energy beam and a weak energy beam, which can adapt to the detection of different reflectivity targets. At the same time, its laser repetition frequency reaches 10 kHz, and the footprint spacing is only 0.7 m, which has great advantages in monitoring polar ice layer changes, global vegetation conditions and clouds and aerosols. However, due to the high sensitivity of ATLAS to photons, it is easily affected by atmospheric scattering, solar radiation and instrument noise, and the obtained photon signal contains a large amount of background noise. Therefore, noise removal is needed for photon data to obtain real signal photons, and through further classification of signal photons, the ground surface and canopy top are obtained to serve scientific problems such as forest inversion, carbon calculation and height measurement. SUMMARY

[0003] In view of the photon denoising and photon classification problems of ICESat-2 forest area, a fine and high-reliability ground surface and canopy top extraction method and system are proposed, which performs multi-level denoising on photon data through adaptive threshold ellipsoid filtering, dense noise removal and improved LOF algorithm to accurately obtain signal photons, and then uses the improved local direction center algorithm to finely extract the external photons of the signal photons, and based on the adaptive size moving window combined with the single-axis inverse distance weight statistical method, the ground surface photons, the canopy and the canopy top photons are distinguished.

[0004] In order to achieve the above purpose, the technical scheme provided by the application is a fine and high-reliability ground surface and canopy top extraction method applied to ICESat-2 data, which comprises the following steps:

[0005] Step 1, using an adaptive threshold ellipsoid filtering method to filter the photon data, and roughly distinguishing noise photons and signal photons;

[0006] Step 1.1, set the threshold of the number of photons for ellipsoid filtering, calculate the radius of each axis of the ellipsoid;

[0007] Step 1.2, screen the potential filtering point set with each photon as the center, and roughly divide the photon group into noise photon group and signal photon group;

[0008] Step 2, further remove the noise photon group in the signal photon group by using the dense area denoising method;

[0009] Step 2.1, merge the signal photon groups according to the overlap rate;

[0010] Step 2.2, re-screen the unmerged signal photon groups;

[0011] Step 2.3, remove the noise photon group according to the distance between the photon groups;

[0012] Step 3, remove the discrete single noise photon close to the signal photon based on the improved LOF algorithm;

[0013] Step 4, classify the signal photons after removing the noise to obtain the ground surface and canopy top photons;

[0014] Step 4.1, distinguish the external photons and internal photons based on the improved local direction center algorithm;

[0015] Step 4.2, distinguish the ground surface photons and canopy top photons in the external photons based on the single-axis inverse distance weight statistical method.

[0016] Moreover, the fixed filtering threshold is used in step 1.1 to inversely calculate the radius of each axis of the ellipsoid to adaptively calculate the expression of the ellipsoid, so as to calculate the radius r of the along-track axis at For example, first, project the photons onto the along-track axis direction, then divide the region where the photons are located into 2n equal intervals, then sample n intervals at intervals, and calculate the average AveSP of the number of photons in the sampling region at Then, calculate the average interval distance AveDis at And the ratio of the average AveSP of the number of photons at Finally, according to the filtering photon number threshold T sph And the ellipsoid amplification coefficient R cof The radius formula r of the along-track axis is obtained at As follows:

[0017]

[0018] AveDis at =L at / (2×n) (2)

[0019]

[0020] where L at represents the length of the photon distribution along the track axis, pn t represents the number of photons in the tth segment of the sampling segment.

[0021] the radius of the track axis r ct and the radius of the elevation axis r h The calculation process is the same as the radius of the track axis r at , and the same way to obtain r ct and r h After that, the ellipsoid expression is as follows:

[0022]

[0023] where c(c at , c ct , c h ) is the coordinate of the photon at the center of the ellipsoid.

[0024] Moreover, in step 1.2, according to the principle of proximity similarity, the following rules are set: taking each photon as the center of the ellipsoid, all the photons in the ellipsoid are collected as a group, if the number of photons in the group is less than the set threshold of the filtered photon number, it is directly considered that all the photons in the group are noise photons, called noise photon group, otherwise, it is considered that all the photons are signal photons, called signal photon group, as shown in the following formula:

[0025]

[0026]

[0027] where group_ng represents the noise photon group, group_sg represents the signal photon group, sum() represents the sum of the number of photons in the set, o represents the number of candidate photons, p l_at , p l_ct and p l_h respectively represent the value of the lth candidate photon point along the track axis, the track axis and the elevation axis, r at , r ct and r h respectively represent the radius of the ellipsoid in the direction of the track axis, the track axis and the elevation axis, T sph is the set threshold of the filtered photon number.

[0028] In order to quickly and efficiently filter out most of the noise photons, the photons in the noise photon group and the signal photon group are not filtered as the center photon, the photons in the signal photon group are counted into the filtering range of the next center photon, and the photons in the noise photon group are not counted.

[0029] Further, in the step 2.1, the number of photons in the overlapping part is first obtained according to the number of overlaps of the photons in two adjacent signal photon groups, and then the ratio of the number of photons in the overlapping part to the number of photons in the signal photon group with fewer photons is taken as the overlap rate, when the overlap rate is greater than a set overlap rate threshold, the two signal photon groups are merged, the overlap rate of the merged signal photon group and its adjacent signal photon group is calculated, and whether the merging needs to be performed again is determined according to the size relationship between the overlap rate and the threshold. After the merging judgment operation is performed on all signal photon groups, a plurality of merged signal photon groups and a plurality of unmerged signal photon groups are obtained.

[0030] Further, in the step 2.2, the signal photon group with the largest number of photons after merging is selected as the baseline signal photon group, and then the signal photon group with the number of photons greater than the baseline signal photon group by a certain proportion is also regarded as the signal photon group after merging. After screening, the number of photons in the signal photon group that does not meet the condition is relatively small, which is called an outlier photon group.

[0031] Further, in the step 2.3, the noise photon group is removed by calculating the distance between the outlier photon group and the signal photon group after merging, the distance between the photon groups refers to the minimum value of the minimum distance between the photons in the two photon groups, and the calculation formula is as follows:

[0032] h(A,B)=min{min{d(a,b)}} (7)

[0033] In the formula, h(A,B) represents the distance between the photon group A and the photon group B, a and b represent the photons in the photon group A and the photon group B respectively, min represents the minimum value, d(a,b) represents the distance between a and b, and the larger h(A,B) is, the farther the two photon groups are, and the more likely the two photon groups are noise photon groups.

[0034] The threshold value T is set dh is the average value of the distance between all photons, if the distance between the outlier photon group and the signal photon group after merging is greater than the threshold value, the outlier photon group is regarded as a noise photon group.

[0035] Further, in the step 3, the LOF algorithm is improved, and an elliptical region search is used instead of a circular region search to adapt to the characteristics that the density distribution of the photons is higher in the horizontal direction, and the major axis and the minor axis of the ellipse are respectively n1r at and n2r h , wherein r at is the radius along the track axis, r h is the radius along the height axis, and n1 and n2 are constants.

[0036] Further, in the step 4.1, the local direction centrality metric DCM is used to quantitatively measure the internal points and the boundary points, and the specific calculation method is as follows:

[0037]

[0038] In the formula, α i Let represent the i-th angle formed by the i-th nearest neighbor point to the center point, π be pi, k be the number of photons within the ellipse, and m1r be the major and minor axes of the ellipse. at m2r h r at r is the radius along the track axis. h Let m1 be the radius of the elevation axis, and m2 be constants.

[0039] when The condition is met when k points are evenly distributed around the center point and all angles are equal; the minimum value of DCM is 0. The maximum value of DCM is 4(k-1)π when one angle is 2π and the other angles are 0. 2 / k 2 Therefore, the DCM is normalized to [0,1] as follows:

[0040]

[0041] In the formula, α i Let represent the i-th angle formed by the i-th nearest neighbor point to the center point, π be pi, k be the number of photons within the ellipse, and m1r be the major and minor axes of the ellipse. at m2r h r at r is the radius along the track axis. h Let m1 be the radius of the elevation axis, and m2 be constants.

[0042] After obtaining the DCM′ value of each point, they are sorted. Points with smaller DCM′ values ​​are more likely to be internal points, while points with larger DCM′ values ​​are more likely to be external points. The number of external points is adjusted multiple times, and the corresponding classification accuracy is recorded. Finally, the data with the highest classification accuracy is selected, and the surface photons and canopy top photons are identified as external photons, while the canopy photons are identified as internal photons.

[0043] Furthermore, in step 4.2, a window is constructed along the track direction with each external photon as the center, and the external photon and other signal photons in the window are classified according to the inverse distance weighted sum in the elevation direction to determine whether it is a surface photon or a canopy top photon.

[0044] Window size and classification criteria are as follows:

[0045] W at =2×w cof ×r at (10)

[0046]

[0047]

[0048] In the formula, W at represents the window width in the track direction, w cof and r at respectively represent the window size coefficient and the radius of the ellipsoid in the track direction, TOC represents that the external photon belongs to the top of the canopy, ground represents that the external photon belongs to the ground, h c represents the elevation value of the external photon to be calculated, assuming that there are qh photons in the window whose elevation is greater than h c , then the elevation of the ph photon is h ph , assuming that there are ql photons in the window whose elevation is less than or equal to h c , then the elevation of the pl photon is h pl .

[0049] The application also provides a fine and high-reliability ground and canopy top extraction system applied to ICESat-2 data, which is used to realize the fine and high-reliability ground and canopy top extraction method applied to ICESat-2 data.

[0050] Moreover, the application comprises a processor and a memory, the memory is used to store program instructions, and the processor is used to call the program instructions in the memory to execute the fine and high-reliability ground and canopy top extraction method applied to ICESat-2 data.

[0051] Compared with the prior art, the application has the following advantages:

[0052] 1) The fixed filtering threshold value is used to inversely calculate the radii of the ellipsoid in all directions to adaptively calculate the ellipsoid expression, thereby reducing the number of parameters and reducing the sensitivity of the method to parameters; 2) The photon data is subjected to multi-stage denoising through adaptive threshold ellipsoid filtering, dense noise removal and improved LOF algorithm, thereby greatly reducing the influence of noise on subsequent photon classification; 3) The improved local direction center algorithm is used to obtain accurate and dense external photons in the signal photons, and the single-axis inverse distance weighted statistical method is used to accurately classify the signal photons and obtain the ground and the canopy top. BRIEF DESCRIPTION OF DRAWINGS

[0053] Figure 1 is a flowchart of an embodiment of the application. DETAILED DESCRIPTION

[0054] The application provides a fine and high-reliability ground and canopy top extraction method and system applied to ICESat-2 data, and the technical solutions of the application are further described below with reference to the drawings and embodiments.

[0055] Example 1

[0056] like Figure 1 As shown, this invention provides a method for extracting fine and highly reliable surface and canopy top data from ICESat-2 data, comprising the following steps:

[0057] Step 1: Use the adaptive threshold ellipsoidal filtering method to filter the photon data and roughly distinguish between noise photons and signal photons.

[0058] To reduce the sensitivity of the filtering method to parameters and decrease the number of manual adjustments, while also considering noise removal in both the rail and vertical directions (noise in the vertical direction is often overlooked by researchers), an adaptive threshold ellipsoidal filtering method is proposed for coarse noise removal.

[0059] Step 1.1: Set the photon count threshold for ellipsoidal filtering and calculate the radii of each axis of the ellipsoid.

[0060] This invention adaptively calculates the ellipsoid expression by using a fixed filtering threshold to inversely calculate the radii of the ellipsoid in all directions. This is used to calculate the radius r along the orbital axis. at For example, firstly, the photons are projected onto the direction along the orbital axis, and the region where the photons are located is divided into 2n segments with a distance of 2n. Then, n segments are sampled at intervals, and the mean value of the number of photons in the sampled region, AveSP, is calculated. at Then calculate the average interval distance AveDis at and the mean number of photons AveSP at The ratio is used to obtain the unit photon distance, and finally, based on the filtered photon number threshold T... sph And the magnification factor R of the ellipsoid cof The formula for the radius r along the track axis is obtained. at as follows:

[0061]

[0062] AveDis at =L at / (2×n) (2)

[0063]

[0064] In the formula, L at pn represents the length of the photon distribution along the orbital axis. t This represents the number of photons in the t-th segment of the sampling segment.

[0065] vertical rail axis radius r ct and elevation axis radius r h Calculation process and the radius r along the track axis at The same method is used to obtain r. ct and rh After that, the ellipsoid expression is as follows:

[0066]

[0067] In the formula, c(c at ,c ct ,c h ) is the coordinate of the center photon of the ellipsoid.

[0068] Step 1.2, screening the potential filtering point set with each photon as the center, roughly dividing the photon group into a noise photon group and a signal photon group.

[0069] Since the distance between the photons is close, in order to avoid frequent calculation of similar photons and reduce the amount of subsequent calculation, according to the similarity principle, the following rules are set: taking each photon as the center of the ellipsoid, all the photons in the ellipsoid are collected as a group, if the number of photons in the group is less than the set filtering photon number threshold, it is directly considered that all the photons in the group are noise photons, called noise photon group, otherwise it is considered that all the photons are signal photons, called signal photon group, as shown in the following formula:

[0070]

[0071]

[0072] In the formula, group_ng represents the noise photon group, group_sg represents the signal photon group, sum() represents the sum of the number of photons in the set, o represents the number of candidate photon points, p l_at , p l_ct and p l_h respectively represent the value of the lth candidate photon point along the track axis, the vertical track axis and the elevation axis, r at , r ct and r h respectively represent the radius of the ellipsoid in the direction of the track axis, the vertical track axis and the elevation axis, T sph is the set filtering photon number threshold.

[0073] In order to quickly and efficiently filter out most of the noise photons, the photons in the noise photon group and the signal photon group are no longer filtered as the center photon (the center of the ellipsoid), the photons in the signal photon group are counted into the filtering range of the next center photon, and the photons in the noise photon group are not counted.

[0074] Step 2, using the dense area noise removal method to further remove the noise photon group in the signal photon group.

[0075] Due to the influence of clouds or other reasons, there may be non-uniform noise conditions, that is, the noise density in some places may be larger, close to the signal photon density, at which time the adaptive threshold ellipsoid filtering method may not be able to filter it out, and further processing is needed for this kind of situation. Therefore, this step mainly merges the signal photon groups obtained in step 1, finds the signal photon track, and then removes the noise photon groups that do not belong to the signal photon track according to the set rules.

[0076] Step 2.1, merging signal photon groups according to overlap rate.

[0077] Since step 1 divides the photons into multiple small groups, and the distribution of signal photons is generally highly concentrated, there will be some overlap between adjacent signal photon groups. In order to find the correct track of signal photons and use it as a reference for removing noise photon groups, it is necessary to merge the signal photon groups obtained in step 1. First, the number of photons in the overlapping part is obtained according to the number of overlapping photons in two adjacent ellipsoids (signal photon groups), then the ratio of the number of photons in the overlapping part to the number of photons in the ellipsoid (signal photon group) with fewer photons is taken as the overlap rate, when the overlap rate is greater than the set overlap rate threshold, the two signal photon groups are merged, the overlap rate of the merged signal photon group and its adjacent signal photon group is calculated, and according to the size relationship between the overlap rate and the threshold, it is determined whether it needs to be merged again. Since the signal photon distribution is dense, it can be considered that the photon distribution in the ellipsoid of the signal photon group is relatively uniform, and then when merging the signal photon groups, the overlap rate threshold can be represented by the volume. According to the rules set in step 1.2, the ratio of the intersection volume of two ellipsoids to the volume of the ellipsoid should be less than 0.5, and the threshold is generally set to [0.1, 0.3]. After the merging judgment operation on all signal photon groups, a number of merged signal photon groups and a number of unmerged signal photon groups are obtained.

[0078] Step 2.2, re-screening of unmerged signal photon groups.

[0079] The merged signal photon group can basically represent the track of the signal photon, but there are the following two cases that make the signal photon group that should be merged not merged: 1) the special distribution of signal photons (the distribution of photons on one side of the intersection is less) makes the overlap rate of adjacent signal photon groups less than the threshold; 2) special terrain and other ground object interference, causing the signal photon track to be interrupted and unable to be merged. Therefore, the signal photon group with the largest number of photons after merging is selected as the baseline signal photon group, and then the photon group with a number of photons greater than a certain proportion of the baseline signal photon group in the unmerged signal photon group is also regarded as a merged signal photon group. After screening, the number of photons in the signal photon group that does not meet the conditions is generally small, which is called an outlier photon group.

[0080] Step 2.3. Remove the noise photon group according to the distance between the photon groups.

[0081] Step 2.2. The outlier photon group is screened according to the number of photons in the signal photon group. At this time, the outlier photon group includes two categories: 1) the outlier photon group far away from the signal photon trajectory and consisting of dense noise; 2) the photon group consisting of signal photons or part of noise photons, which is close to or in the trajectory, but is not combined due to the distribution of photons not being within the intersection range. The noise photon group is removed by calculating the distance between the outlier photon group and the combined signal photon group. The distance between the photon groups refers to the minimum value of the minimum distance between the photons in the two photon groups. The calculation formula is as follows:

[0082] h(A, B) = min{min{d(a, b)}} (7)

[0083] In the formula, h(A, B) represents the distance between the photon group A and the photon group B, a and b represent the photons in the photon group A and the photon group B, min represents the minimum value, d(a, b) represents the distance between a and b, and h(A, B) is larger, indicating that the two photon groups are farther apart, and it is more likely to be a noise photon group.

[0084] Set the threshold value T dh as the average distance between all photons. If the distance between the outlier photon group and the combined signal photon group is greater than the threshold value, the outlier photon group is regarded as a noise photon group.

[0085] Step 3. Remove the discrete single noise photon close to the signal photon based on the improved LOF algorithm.

[0086] In order to further remove the discrete single noise photon close to the signal photon, the Local Outlier Factor (LOF) algorithm is improved. The basic idea of the LOF algorithm is to calculate the score of a point, which represents the local density between the given point and its adjacent points, and the outlier points are considered to be the points with a density level significantly lower than the threshold score. The LOF algorithm is improved in the present application, and an elliptical region search is used instead of a circular region search to adapt to the characteristics of the high density distribution of photons in the horizontal direction. The major axis and the minor axis of the ellipse are n1r at , n2r h , where r at is the radius along the track axis, r h is the radius of the elevation axis, and n1, n2 are constants. Since the photon points on the signal photon trajectory are relatively dense, and the noise points around them exhibit a certain outlier trend locally, the improved Local Outlier Factor algorithm can better remove this part of noise.

[0087] Step 4. Classify the noise-removed signal photons to obtain the ground and canopy top photons.

[0088] Step 4.1. Distinguish the external photons and internal photons based on the improved local direction centrality algorithm.

[0089] The core idea of the local direction centrality algorithm is to measure the K-Nearest Neighbor (KNN) distribution of each point to distinguish the boundary points (external points) and internal points of a group of points. The algorithm considers that the boundary points of a group of points should outline the contour of the group of points, and the internal points should be surrounded by their neighboring points in all directions, while the boundary points will only have neighboring points in certain directions. In order to quantitatively measure the internal points and boundary points, the local direction centrality metric DCM is proposed, which is calculated as follows:

[0090]

[0091] wherein αi represents the ith angle formed by the nearest i points of the center point, k is the number of nearest points of the center point, and π is the circular constant. i

[0092] For two-dimensional angles, when is established, i.e. only when k points are uniformly distributed around the center point, all angles are equal, and DCM reaches the minimum value 0; when one angle is 2π and the other angles are 0, the maximum value of DCM is 4(k-1)π 2 / k 2 ; therefore, DCM can be normalized to [0,1] as follows:

[0093]

[0094] wherein αi represents the ith angle formed by the nearest i points of the center point, k is the number of nearest points of the center point, and π is the circular constant. i

[0095] The local direction centrality algorithm uses KNN to find the surrounding points of the center point, which can better distinguish whether the center point is a boundary point or an internal point when applied to a group of clustered points. However, for signal photons, signal photons are composed of ground photons, canopy photons and canopy top photons, and when the vegetation is high, there may be a large gap between canopy photons and ground photons, and there are few signal photons in the gap area. Therefore, using KNN to find the surrounding points of the center photon will lead to misclassification of the center photon between the canopy and the ground. Therefore, the process of using KNN to find the surrounding points of the center point is modified to use an elliptical range to find the surrounding photons of the center point, i.e. k represents the number of photons within the elliptical range, and the major axis and minor axis of the ellipse are m1r at , m2r h ​​where r at is the radius along the track axis, r h is the radius along the elevation axis, and m1 and m2 are constants. This makes the center point not affected by the interval distance when searching for surrounding points, and avoids the photons between the canopy and the ground being misclassified as external photons. After obtaining the DCM' value of each point, the points are sorted according to the DCM' value. The smaller the DCM' value, the more likely the point is an internal point. The larger the DCM' value, the more likely the point is an external point. The number of external points is adjusted multiple times, and the corresponding classification accuracy is recorded. Finally, the data with the highest classification accuracy is selected. The ground photons and the canopy top photons are obtained as external photons, and the canopy photons are obtained as internal photons.

[0096] Step 4.2, distinguishing the ground photons and the canopy top photons in the external photons based on the single-axis inverse distance weighting statistical method.

[0097] Since the external photons need to be distinguished as ground photons and canopy top photons, the single-axis inverse distance weighting statistical method is used to determine whether the external photons are ground photons or canopy top photons within a certain size window. A window is constructed in the track direction with each external photon as the center. According to the inverse distance weighting summation of the external photon and other signal photons in the window in the elevation direction, the external photon is classified as a ground photon or a canopy top photon.

[0098] The window size and classification criteria are as follows:

[0099] W at = 2 x w cof x r at (10)

[0100]

[0101]

[0102] In the formula, W at represents the window width in the track direction, w cof and r at represent the window size coefficient and the radius of the ellipsoid in the track direction, respectively, TOC represents that the external photon belongs to the canopy top, ground represents that the external photon belongs to the ground, h c represents the elevation value of the external photon to be calculated, qh represents the number of photons in the window with an elevation greater than h c , h ph represents the elevation of the ph-th photon, ql represents the number of photons in the window with an elevation less than or equal to h c , and h pl represents the elevation of the pl-th photon.

[0103] Example Two

[0104] Based on the same inventive concept, the application further provides a fine and high-reliability ground and canopy top extraction system applied to ICESat-2 data, which comprises a processor and a memory, the memory is used for storing program instructions, and the processor is used for calling the program instructions in the memory to execute the fine and high-reliability ground and canopy top extraction method applied to ICESat-2 data.

[0105] In the implementation, the method provided by the technical scheme of the application can be automatically run by a computer software technology, and the system device of the method, such as a computer readable storage medium storing the corresponding computer program of the technical scheme of the application and a computer device including the running of the corresponding computer program, should also be within the protection scope of the application.

[0106] The specific embodiments described herein are merely illustrative of the spirit of the application. Those skilled in the art of the application can make various modifications or supplements to the described specific embodiments or replace them with similar ways, without departing from the spirit of the application or exceeding the scope defined by the appended claims.

Claims

1. A method for fine and reliable extraction of ground and canopy top from ICESat-2 data, characterized in that, The method comprises the following steps: Step 1, filtering photon data by using an adaptive threshold ellipsoid filtering method to roughly distinguish noise photons and signal photons; Step 1.1, setting the photon number threshold of ellipsoid filtering, and calculating the radius of each axis of the ellipsoid; Adaptive calculation of the ellipsoid expression is performed by using a fixed filtering threshold to inversely calculate the semi-axes of the ellipsoid, so as to calculate the radius along the track axis For example, first, the photons are projected along the track axis direction, and the area where the photons are located is evenly divided into 2 n segments at a distance of 1 n segment interval, and then the sampling is performed at the interval The average value of the number of photons in the sampling area is calculated The ratio of the average interval distance to the average value of the number of photons is calculated to obtain the unit photon distance, and finally, the radius formula along the track axis is obtained according to the filtering photon number threshold and the ellipsoid magnification coefficient as follows:​ （1） （2） （3） wherein denotes the length of the distribution of photons along the track axis, denotes the number of photons in the segment, t denotes the number of photons in the segment, Pitch axis radius and the elevation axis radius The calculation process is the same as along the track axis radius and the same way to obtain and After that, the ellipsoid expression is as follows: （4） In the formula, is the coordinate of the ellipsoid center photon; Step 1.2, screening a potential filtering point set with each photon as the center, and roughly dividing the photon group into a noise photon group and a signal photon group; Step 2, further removing the noise photon group in the signal photon group by using a dense area noise removal method; Step 2.1, merging the signal photon groups according to the overlap rate; Step 2.2, re-screening the unmerged signal photon groups; Step 2.3, removing the noise photon group according to the distance between the photon groups; Step 3, removing discrete single noise photons close to the signal photons based on an improved LOF algorithm; The LOF algorithm is improved by using an elliptical region search instead of a circular region search to adapt to the higher density distribution of photons in the horizontal direction. The major axis and the minor axis of the ellipse are respectively , wherein is the radius along the track axis, is the radius along the elevation axis, , is a constant; Step 4, classifying the signal photons after removing the noise to obtain ground surface and canopy top photons; Step 4.1, distinguishing external photons and internal photons based on an improved local direction centrality algorithm; The local direction centrality measure DCM is used to quantitatively measure internal points and boundary points, and the specific calculation method is as follows: （8） In the formula, the first angle formed by the nearest neighbors of the center point, i the second angle formed by the nearest neighbors of the center point, i the third angle formed by the nearest neighbors of the center point, pi, k the number of photons in the elliptical range, the major axis and the minor axis of the ellipse are , , the radius along the track axis, the radius along the elevation axis, , a constant; When holds, i.e. only when k are evenly distributed around the center point, all angles are equal, DCM takes the minimum value 0; when one of the angles is and all other angles are 0, the maximum value DCM of is obtained; thus DCM is normalized to [0, 1] as follows: （9） In the formula, the first angle formed by the nearest neighbors of the center point, i the second angle formed by the nearest neighbors of the center point, i the third angle formed by the nearest neighbors of the center point, pi, k the number of photons in the elliptical range, the major axis and the minor axis of the ellipse are , , the radius along the track axis, the radius along the elevation axis, , a constant; After the value of each point is obtained , it is sorted, , the smaller the value, the more likely the point is an interior point, , the larger the value, the more likely the point is an exterior point; the number of exterior points is adjusted multiple times, and the corresponding classification accuracy is recorded, and finally the group of data with the highest classification accuracy is selected, and finally the ground photons and the top photons of the canopy are obtained as the exterior photons, and the canopy photons are obtained as the interior photons; Step 4.2, distinguishing ground surface photons and canopy top photons in the external photons based on a single-axis inverse distance weight statistical method.

2. The method for fine and highly reliable surface and canopy top extraction applied to ICESat-2 data according to claim 1, characterized in that: In step 1.2, according to the similarity principle, the rules are as follows: taking each photon as the center of the ellipsoid, all the photons in the ellipsoid are collected as a group, if the number of photons in the group is less than the set filtering photon number threshold, it is directly considered that all the photons in the group are noise photons, which is called a noise photon group, otherwise it is considered that all the photons are signal photons, which is called a signal photon group, as shown in the following formula: （5） （6） wherein, represents the set of noise photons, represents the set of signal photons, represents summing the number of photons within the set, o represents the number of candidate photons, , and represent the value of the l th candidate photon point along the along-track, cross-track, and height axes, respectively, , and represent the radii of the ellipsoid along the along-track, cross-track, and height axes, respectively, is a set filter photon number threshold value; In order to quickly and efficiently filter out most of the noise photons, the photons in the noise photon group and the signal photon group are not used as the center photon for filtering, the photons in the signal photon group are counted into the filtering range of the next center photon, and the photons in the noise photon group are not counted.

3. The method for fine and highly reliable surface and canopy top extraction applied to ICESat-2 data of claim 1, wherein: In step 2.1, the number of photons in the overlapping part is first calculated according to the number of overlapping photons in two adjacent signal photon groups, and then the ratio of the number of photons in the overlapping part to the number of photons in the signal photon group with fewer photons is taken as the overlap rate. When the overlap rate is greater than the set overlap rate threshold, the two signal photon groups are merged, the overlap rate of the merged signal photon group and its adjacent signal photon group is calculated, and whether it needs to be merged again is determined according to the size relationship between the overlap rate and the threshold. After the merging judgment operation is performed on all signal photon groups, a plurality of merged signal photon groups and a plurality of unmerged signal photon groups are obtained.

4. The method for fine and highly reliable surface and canopy top extraction applied to ICESat-2 data of claim 1, wherein: In step 2.2, the signal photon group with the largest number of photons after merging is selected as the baseline signal photon group, and then the signal photon groups with the number of photons greater than a certain proportion of the baseline signal photon group are also regarded as the signal photon groups after merging. After screening, the number of photons in the signal photon group that does not meet the condition is relatively small, which is called an outlier photon group.

5. The method for fine and highly reliable surface and canopy top extraction applied to ICESat-2 data of claim 1, wherein: In step 2.3, the distance between the outlier photon group and the merged signal photon group is calculated to remove the noise photon group. The distance between the photon groups refers to the minimum value of the minimum distance between the photons in the two photon groups, and the calculation formula is as follows: （7） wherein represents a photon group A and a photon group B the distance between a and b represents a photon in a photon group A and a photon group B min represents a minimum value, represents a and b the distance between the greater the distance between the two photon groups, the more likely it is a noise photon group; Setting threshold T dh is the average distance between all photons, and if the distance between the outlier photon group and the combined signal photon group is greater than the threshold, the outlier photon group is considered a noise photon group.

6. The method for fine and highly reliable surface and canopy top extraction applied to ICESat-2 data of claim 1, wherein: In step 4.2, a window is constructed along the track direction, centered on each external photon, and a classification is made according to the inverse distance weighted sum of the external photon and other signal photons within the window in the elevation direction, to determine whether it is a ground photon or a canopy top photon; the window size and classification criteria are as follows: （10） （11） （12） wherein, represents the window width in the track direction, and represent the window size coefficient and the radius of the ellipsoid in the track direction, respectively, TOC represents that the external photon belongs to the top of the canopy, ground represents that the external photon belongs to the ground surface, represents the elevation value of the external photon to be calculated, assuming that there are qh photons in the window whose elevations are greater than , the elevation of the ph th photon is , assuming that there are ql photons in the window whose elevations are less than or equal to , the elevation of the pl th photon is .

7. A fine high-reliability surface and canopy top extraction system applied to ICESat-2 data, characterized by, The application comprises a processor and a memory, the memory is used for storing program instructions, and the processor is used for calling the program instructions in the memory to execute the fine and high-reliability ground and canopy top extraction method applied to ICESat-2 data according to any one of claims 1-6.

Images