Method and device for screening elevation control points based on ICESat-2 data

Abstract

The invention provides a method and device for screening elevation control points based on ICESat-2 data, and relates to the technical field of data processing.The method comprises the steps that ATL03 and ATL08 of ICESat-2 are combined, all original ground photon points are extracted from the ATL03, all the original ground photon points in every 100-meter segments in the ATL03 are converted to a horizontal datum plane with every 100-meter segments in the ATL08 in the rail direction as the unit, and all the original ground photon points in every 100-meter segments in the ATL03 are extracted; the method comprises the following steps: obtaining ground photon points after terrain removal, constructing an elevation frequency histogram and calculating skewness of the elevation frequency histogram and ground photon vertical aggregation degree of each 100-meter segment, screening potential elevation control points in the ATL08 based on attribute parameters of each potential elevation control point in the ATL08 and the corresponding ground photon vertical aggregation degree and skewness, and obtaining the potential elevation control points in the ATL08. Therefore, the extraction precision of the elevation control points is effectively improved.

Description

Technical Field

[0001] The present invention relates to the technical field of data processing, and in particular to a method and device for screening elevation control points based on ICESat-2 data. Background Art

[0002] Elevation control points are mainly used for fine terrain mapping, and their extraction accuracy directly affects the accuracy of the final terrain mapping results. If the accuracy of the elevation control points is not high, the terrain fitting error will increase, thereby reducing the reliability of the mapping results.

[0003] The existing ICESat-2 elevation control point screening method mainly relies on the ATL08 product, and the elevation control points are screened by setting empirical thresholds. However, since only ATL08 parameters are used, it is difficult to effectively eliminate the anomalies of the vertical discrete distribution of ground points (such as multiple scattering, complex terrain, and confusion between low vegetation and ground points), resulting in limited accuracy of elevation control points. Summary of the invention

[0004] The object of the present invention is to provide a method and device for screening elevation control points based on ICESat-2 data, so as to alleviate the technical problem of poor elevation control point extraction accuracy existing in the existing ICESat-2 elevation control point screening method.

[0005] In a first aspect, the present invention provides a method for screening elevation control points based on ICESat-2 data, comprising: combining the ATL03 data product and the ATL08 data product of ICESat-2 to extract all original ground photon points from the ATL03 data product; taking the ATL08 data product in the along-track direction of each 100-meter segment as a unit, converting all original ground photon points in each 100-meter segment of the ATL03 data product to a horizontal reference plane to obtain all ground photon points after removing the terrain in each 100-meter segment; constructing an elevation frequency histogram for each 100-meter segment based on all ground photon points after removing the terrain in each 100-meter segment, and calculating the skewness of the elevation frequency histogram and the elevation frequency of each 100-meter segment. The vertical concentration of ground photons; wherein, the horizontal axis of the elevation frequency histogram represents the elevation interval, and the vertical axis represents the proportion of the number of ground photon points falling into the elevation interval to the total number of ground photon points in the 100m segment to which it belongs; the vertical concentration of ground photons represents the ratio of the peak of the elevation frequency histogram to the height standard deviation of the ground photon points in the 100-meter segment; based on the attribute parameters of each potential elevation control point in the ATL08 data product, the vertical concentration of ground photons in the 100-meter segment to which it belongs and the skewness of the corresponding elevation frequency histogram, the potential elevation control points in the ATL08 data product are screened to obtain the elevation control point screening results; wherein, the potential elevation control point represents the ATL08 ground photon point used for elevation control screening.

[0006] In an alternative embodiment, the ATL03 data product and the ATL08 data product of ICESat-2 are combined to extract all the original ground photon points from the ATL03 data product, including: extracting the attribute parameters of each photon point in the ATL03 data product and the ATL08 data product; spatially matching the ATL03 data product and the ATL08 data product based on the attribute parameters; using the classification information of the target photon points in the ATL08 data product as the classification label of the photon points in the ATL03 data product that are spatially matched with them; where the target photon points represent any one of all the photon points in the ATL08 data product; determining all the original ground photon points based on the classification labels of all the photon points in the ATL03 data product.

[0007] In an alternative embodiment, taking every 100-meter segment along the track direction of the ATL08 data product as a unit, all the original ground photon points within each 100-meter segment in the ATL03 data product are converted to the horizontal datum plane, including: obtaining the along-track distance and elevation value of each original ground photon point within the target 100-meter segment in the ATL03 data product; where the target 100-meter segment represents any one of all the 100-meter segments in the ATL03 data product; fitting the along-track distance and elevation value of the original ground photon points within the target 100-meter segment to obtain the terrain slope of the target 100-meter segment; constructing a two-dimensional rotation matrix based on the terrain slope; and converting each original ground photon point within the target 100-meter segment to the horizontal datum plane based on the two-dimensional rotation matrix to obtain the corresponding ground photon points after removing the terrain.

[0008] In an alternative embodiment, fitting the along-track distance and elevation value of the original ground photon points within the target 100-meter segment to obtain the terrain slope of the target 100-meter segment includes: constructing a unary linear regression equation for the elevation and along-track distance of the original ground photon points; where the elevation of the original ground photon points is the dependent variable of the unary linear regression equation, the along-track distance of the original ground photon points is the independent variable of the unary linear regression equation, and the terrain slope of the segment to which the original ground photon points belong is the slope of the unary linear regression equation; substituting the along-track distance and elevation value of the original ground photon points within the target 100-meter segment into the unary linear regression equation to obtain a target system of equations; and solving the target system of equations using the least squares method to obtain the terrain slope of the target 100-meter segment.

[0009] In an optional embodiment, based on all ground photon points after removing the terrain in each 100-meter segment, an elevation frequency histogram of each 100-meter segment is constructed, and the skewness of the elevation frequency histogram and the vertical concentration of ground photons in each 100-meter segment are calculated, including: grouping all ground photon points after removing the terrain in the target 100-meter segment at preset elevation intervals, and counting the number of ground photon points falling into each elevation interval; wherein the target 100-meter segment represents any segment of all 100-meter segments in the ATL03 data product; taking the elevation interval as the horizontal axis, the ground photon points falling into the elevation interval are counted. The ratio of the number of photon points to the total number of ground photon points in the 100m segment to which it belongs is taken as the vertical axis, and the elevation frequency histogram of the target 100-meter segment is constructed; the peak of the elevation frequency histogram of the target 100-meter segment is determined; the height standard deviation of all ground photon points after removing the terrain in the target 100-meter segment is calculated; based on the height standard deviation and the elevation of all ground photon points after removing the terrain in the target 100-meter segment, the skewness of the elevation frequency histogram of the target 100-meter segment is calculated; the ratio of the peak and the height standard deviation is calculated, and the ratio result is used as the vertical concentration of ground photons in the target 100-meter segment.

[0010] In an optional implementation, based on the attribute parameters of each potential elevation control point in the ATL08 data product, the ground photon vertical concentration of the 100-meter segment to which it belongs, and the skewness of the corresponding elevation frequency histogram, the potential elevation control points in the ATL08 data product are screened, including: calculating the difference between the elevation of the target potential elevation control point and its reference elevation; wherein the target potential elevation control point represents any potential elevation control point in the ATL08 data product; and eliminating the target potential elevation control point when it is determined that the target potential elevation control point satisfies any of the following preset conditions; wherein the preset condition Including: the absolute value of the difference is greater than the first threshold, the cloud cover identifier is greater than the second threshold, the surface type identifier is equal to the third threshold, the number of ground photon points in the 100-meter segment is less than or equal to the fourth threshold, the total number of photon points in the 100-meter segment is greater than or equal to the fifth threshold, the proportion of ground photon points in the 100-meter segment is less than or equal to the sixth threshold, the point discrete distribution identifier is equal to 1, the sum of the ground photon distribution identifiers of the 100-meter segment is not equal to 5, the vertical concentration of ground photons in the 100-meter segment is less than the seventh threshold, and the absolute value of the skewness of the elevation frequency histogram of the 100-meter segment is greater than the eighth threshold.

[0011] In a second aspect, the present invention provides a device for screening elevation control points based on ICESat-2 data, comprising: an extraction module for combining the ATL03 data product and the ATL08 data product of ICESat-2 to extract all original ground photon points from the ATL03 data product; a conversion module for converting all original ground photon points in each 100-meter segment of the ATL08 data product in the along-track direction to a horizontal reference plane, so as to obtain all ground photon points after removing the terrain in each 100-meter segment; a calculation module for constructing an elevation frequency histogram of each 100-meter segment based on all ground photon points after removing the terrain in each 100-meter segment, and calculating the skewness and The vertical concentration of ground photons in each 100-meter segment; wherein, the horizontal axis of the elevation frequency histogram represents the elevation interval, and the vertical axis represents the proportion of the number of ground photon points falling into the elevation interval to the total number of ground photon points in the 100m segment to which it belongs; the vertical concentration of ground photons represents the ratio of the peak of the elevation frequency histogram to the height standard deviation of the ground photon points in the 100-meter segment; the screening module is used to screen the potential elevation control points in the ATL08 data product based on the attribute parameters of each potential elevation control point in the ATL08 data product, the vertical concentration of ground photons in the 100-meter segment to which it belongs, and the skewness of the corresponding elevation frequency histogram, to obtain the elevation control point screening results; wherein, the potential elevation control point represents the ATL08 ground photon point used for elevation control screening.

[0012] In an optional embodiment, the extraction module is specifically used to: extract attribute parameters of each photon point in the ATL03 data product and the ATL08 data product; spatially match the ATL03 data product and the ATL08 data product based on the attribute parameters; use the classification information of the target photon point in the ATL08 data product as the classification label of the photon point in the ATL03 data product that matches it spatially; wherein the target photon point represents any photon point among all photon points in the ATL08 data product; and determine all original ground photon points based on the classification labels of all photon points in the ATL03 data product.

[0013] In a third aspect, the present invention provides an electronic device, including a memory and a processor, wherein the memory stores a computer program that can be run on the processor, and when the processor executes the computer program, the method for screening elevation control points based on ICESat-2 data described in any one of the aforementioned embodiments is implemented.

[0014] In a fourth aspect, the present invention provides a computer-readable storage medium storing computer instructions, which, when executed by a processor, implement the method for screening elevation control points based on ICESat-2 data as described in any one of the aforementioned embodiments.

[0015] The present invention provides a method for screening elevation control points based on ICESat-2 data. First, the ATL03 data product and the ATL08 data product of ICESat-2 are combined to extract all original ground photon points from the ATL03 data product. Then, taking every 100-meter segment along the track direction of the ATL08 data product as a unit, all the original ground photon points within each 100-meter segment in the ATL03 data product are converted to the horizontal reference plane, obtaining all the ground photon points after terrain removal within each 100-meter segment. Furthermore, an elevation frequency histogram that can accurately reflect the photon aggregation degree is constructed for each 100-meter segment, and the skewness of the elevation frequency histogram and the vertical aggregation degree of ground photons in each 100-meter segment are calculated. Finally, based on the attribute parameters of each potential elevation control point in the ATL08 data product, the vertical aggregation degree of ground photons in its affiliated 100-meter segment, and the skewness of the corresponding elevation frequency histogram, the potential elevation control points in the ATL08 data product are screened to obtain the screening result of elevation control points. By combining the two data products of ATL03 and ATL08, the present invention proposes two new parameters for screening elevation control points: the vertical aggregation degree of ground photons and skewness. Compared with screening only using the ATL08 parameters, it can effectively improve the extraction accuracy of elevation control points and provide technical support for global elevation control point extraction and fine terrain mapping. BRIEF DESCRIPTION OF THE DRAWINGS

[0016] In order to more clearly illustrate the specific embodiments of the present invention or the technical solutions in the prior art, the following will briefly introduce the drawings required for use in the description of the specific embodiments or the prior art. Obviously, the drawings in the following description are some embodiments of the present invention. For those of ordinary skill in the art, without creative efforts, other drawings can also be obtained based on these drawings.

[0017] Figure 1 It is a flowchart of a method for screening elevation control points based on ICESat-2 data provided by an embodiment of the present invention; Figure 2 It is a frequency distribution histogram of the ground elevation error between the original data and the elevation control points in the flat area; Figure 3 It is a frequency distribution histogram of the ground elevation error between the original data and the elevation control points in the hilly area; Figure 4 It is a frequency distribution histogram of the ground elevation error between the original data and the elevation control points in the mountainous area; Figure 5 It is a comparison result diagram of the elevation control point accuracy verification for the ATL08 original data, without adding new parameters, and after adding new parameters; Figure 6A result comparison diagram between the solution provided by the embodiment of the present invention and the elevation control point screening solution proposed by Li et al., (2021); Figure 7 A functional module diagram of a device for screening elevation control points based on ICESat-2 data provided by an embodiment of the present invention; Figure 8 A schematic diagram of an electronic device provided by an embodiment of the present invention. Specific embodiments

[0018] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings in the embodiments of the present invention. Obviously, the described embodiments are some, but not all, of the embodiments of the present invention. Usually, the components of the embodiments of the present invention described and illustrated herein can be arranged and designed in various different configurations.

[0019] Therefore, the following detailed description of the embodiments of the present invention provided in the drawings is not intended to limit the scope of the claimed present invention, but merely represents selected embodiments of the present invention. All other embodiments obtained by those of ordinary skill in the art based on the embodiments of the present invention without creative efforts fall within the scope of protection of the present invention.

[0020] The following will describe in detail some embodiments of the present invention with reference to the accompanying drawings. Without conflict, the following embodiments and the features in the embodiments can be combined with each other.

[0021] Embodiment 1 Figure 1 A flowchart of a method for screening elevation control points based on ICESat-2 data provided by an embodiment of the present invention, as Figure 1 shown, the method specifically includes the following steps: Step S102, combining the ATL03 data product and the ATL08 data product of ICESat-2 to extract all original ground photon points from the ATL03 data product.

[0022] ICESat-2 / ATLAS provides 23 standard data products (ATL00~ATL23, except ATL05). The embodiments of the present invention mainly conduct research on ATL03 and ATL08. Among them, ATL03 calculates the precise geolocation information (i.e., latitude, longitude, and elevation) of each received photon through the photon round-trip time, laser position, and attitude angle. Based on the preprocessing of ATL03, ATL08 removes a large number of noise photons and classifies the remaining photons, including noise points, ground points, canopy points, and top canopy points.

[0023] The ATL08 data product also includes ground photon points (only one point per 100m segment) for elevation control screening, which are referred to as potential elevation control points in the embodiments of the present invention. The method provided by the embodiments of the present invention is to eliminate abnormal elevation control points from these potential elevation control points to obtain high-precision elevation control points, so as to provide technical support for subsequent fine topographic mapping.

[0024] At the beginning of the execution of the method of the embodiments of the present invention, first, the ATL03 data product and the ATL08 data product of ICESat-2 are combined, and the photon points in the above two data products are spatially matched and classified and labeled, so as to assign the classification information of ATL08 to the ATL03 photon points, thereby constructing photon-level classified point cloud data, and then all the original ground photon points in the ATL03 data product are determined according to the classification information.

[0025] Step S104: Taking each 100-meter segment along the track direction of the ATL08 data product as a unit, convert all the original ground photon points within each 100-meter segment of the ATL03 data product to the horizontal reference plane to obtain all the ground photon points after removing the terrain within each 100-meter segment.

[0026] It is known that the degree of aggregation of ground photon points in the vertical direction is closely related to their accuracy. The more dispersed the elevation distribution of ground photon points, the higher the probability that the ground photon points deviate from the true ground. In order to quantify the vertical aggregation degree of the original ground photon points, the embodiments of the present invention combine the ATL08 data product and the ATL03 data product, and take each 100-meter segment along the track direction of the ATL08 data product as a unit, and respectively calculate the elevation frequency histogram of the original ground photon points within each 100-meter segment of the ATL03 data product, extract its peak and elevation standard deviation, and use the peak / elevation standard deviation as the aggregation degree index.

[0027] However, due to the influence of terrain undulation on the distribution of the original ground photon points, if the elevation standard deviation is directly calculated based on the original ground photon points or the elevation frequency histogram is constructed, the photon aggregation degree cannot be accurately reflected. To eliminate the influence of terrain, the embodiments of the present invention convert all the original ground photon points within each 100-meter segment of the ATL03 data product to the horizontal reference plane, that is, remove the influence of terrain undulation, to obtain the ground photon points after removing the terrain. Optionally, a linear regression is used to fit the terrain slope, and a rotation matrix is constructed to transform the ground photon points, so as to obtain the photon distribution after removing the influence of terrain undulation, so as to improve the accuracy of photon aggregation degree evaluation.

[0028] Step S106: Based on all the ground photon points after removing the terrain within each 100-meter segment, construct the elevation frequency histogram of each 100-meter segment, and calculate the skewness of the elevation frequency histogram and the vertical aggregation degree of the ground photons in each 100-meter segment.

[0029] Among them, the horizontal axis of the elevation frequency histogram represents the elevation interval, and the vertical axis represents the proportion of the number of ground photon points falling within the elevation interval to the total number of ground photon points within its affiliated 100m segment; the vertical aggregation degree of ground photons represents the ratio of the peak of the elevation frequency histogram to the height standard deviation of the ground photon points within the 100m segment.

[0030] According to the definitions of the horizontal and vertical axes of the above-mentioned elevation frequency histogram, to construct the elevation frequency histogram for each 100m segment, it is necessary to divide the elevations in this segment at fixed intervals, count the number of ground photon points falling within each elevation interval, and then calculate the proportion to the total number of ground photon points within its affiliated 100m segment.

[0031] In the embodiment of the present invention, the vertical dispersion degree of the ground photon points within the 100m segment is characterized by calculating the height standard deviation of the ground photon points within the 100m segment. The peak of the elevation frequency histogram represents the highest frequency value in the elevation frequency histogram (that is, the maximum value among the proportions of the number of ground photon points falling within the elevation interval to the total number of ground photon points), which is also the main peak of the elevation frequency distribution, corresponding to the elevation value where the ground photons are most concentrated. The skewness of the elevation frequency histogram is used to measure the symmetry of the elevation distribution of the ground photon points, reflecting the skewness of the distribution of the ground photon points. A positive skewness value indicates a right-skewed distribution, and a negative skewness value indicates a left-skewed distribution.

[0032] Step S108, based on the attribute parameters of each potential elevation control point in the ATL08 data product, the vertical aggregation degree of ground photons in its affiliated 100m segment, and the skewness of the corresponding elevation frequency histogram, screen the potential elevation control points in the ATL08 data product to obtain the screening result of the elevation control points.

[0033] Among them, the potential elevation control point refers to the ATL08 ground photon point used for elevation control screening.

[0034] In the prior art, only based on the attribute parameters of the potential elevation control points in the ATL08 data product, the elevation control points are screened by setting empirical thresholds. However, this screening method cannot effectively eliminate abnormal situations where the vertical distribution of ground points is discrete (such as multiple scattering, complex terrain, confusion between low vegetation and ground points). Therefore, the embodiment of the present invention proposes to calculate two parameters, namely, the vertical aggregation degree of ground photons in the 100m segment to which the potential elevation control point belongs and the skewness of the corresponding elevation frequency histogram, and further screen the potential elevation control points in the ATL08 data product based on the above two parameters, so as to eliminate potential elevation control points affected by abnormal elevations, clouds, water bodies, ground object reflections, solar radiation, atmospheric interference, and the vertical dispersion degree of ground photons, achieving the effect of improving the extraction accuracy of elevation control points.

[0035] The embodiment of the present invention provides a method for screening elevation control points based on ICESat-2 data. First, the ATL03 data product and the ATL08 data product of ICESat-2 are combined to extract all original ground photon points from the ATL03 data product. Then, with the ATL08 data product in the along-track direction as a unit, all original ground photon points in each 100-meter segment of the ATL03 data product are converted to a horizontal reference plane to obtain all ground photon points after removing the terrain in each 100-meter segment. Then, an elevation frequency histogram that can accurately reflect the photon concentration is constructed for each 100-meter segment, and the skewness of the elevation frequency histogram and the vertical concentration of ground photons in each 100-meter segment are calculated. Finally, based on the attribute parameters of each potential elevation control point in the ATL08 data product, the vertical concentration of ground photons in the 100-meter segment to which it belongs, and the skewness of the corresponding elevation frequency histogram, the potential elevation control points in the ATL08 data product are screened to obtain the elevation control point screening result. The embodiment of the present invention proposes two new parameters for screening elevation control points by combining the two data products ATL03 and ATL08: ground photon vertical concentration and skewness. Compared with screening using only ATL08 parameters, the extraction accuracy of elevation control points can be effectively improved, providing technical support for global elevation control point extraction and fine terrain mapping.

[0036] In an optional implementation, the above step S102 combines the ATL03 data product and the ATL08 data product of ICESat-2 to extract all original ground photon points from the ATL03 data product, specifically including the following steps: Step S1021, extracting the attribute parameters of each photon point in the ATL03 data product and the ATL08 data product.

[0037] Step S1022: spatially match the ATL03 data product and the ATL08 data product based on the attribute parameters.

[0038] Specifically, firstly, the attribute parameters of each photon point are extracted from the ATL03 data product and the ATL08 data product. The following Table 1 shows the attribute parameters of the photon point in the ATL03 data product, and Table 2 shows the attribute parameters of the photon point in the ATL08 data product.

[0039] Taking the method of spatially matching any classified photon point in the ATL08 data product with the photon point in the ATL03 data product as an example, the above step S1022 will be described. Specifically, to spatially match the photon points in the ATL03 data product and the ATL08 data product, first connect each segmented photon using the ATL03 attribute parameters in Table 1 to obtain the global index of all photons in ATL03. Then, obtain the classed_pc_indx and ph_segment_id corresponding to any classified photon point (hereinafter referred to as the classified photon point) in the ATL08 data product. Finally, calculate the global index position of the classified photon point in the original photon sequence of ATL03 based on its classed_pc_indx and ph_segment_id parameters and the ph_index_beg and segment_id parameters provided by ATL03, and determine the original photon point in the ATL03 data product that matches the classified photon point. Thus, the spatial matching of a pair of photon points is completed.

[0040] Table 1 Attribute parameters of photon points in the ATL03 data product

[0041] Table 2 Attribute parameters of photon points in the ATL08 data product

[0042] Step S1023, use the classification information of the target photon point in the ATL08 data product as the classification label of the photon point in the ATL03 data product that is spatially matched with it; where the target photon point represents any photon point among all photon points in the ATL08 data product.

[0043] In the embodiment of the present invention, the classification information classed_pc_flag of each photon point in the ATL08 data product is respectively assigned to the ATL03 photon that is spatially matched with it, that is, the transmission of the classification information from the ATL08 data product to the ATL03 data product is completed. Furthermore, most photon points in the ATL03 data product have classification labels (that is, classification information), and the classification labels include: noise points, ground points, canopy points, and top canopy points.

[0044] Step S1024, determine all original ground photon points based on the classification labels of all photon points in the ATL03 data product.

[0045] That is to say, extract the photon points in the ATL03 data product with the classification label of ground points, that is, obtain all original ground photon points, as the basic data for subsequent data processing steps.

[0046] In an alternative embodiment, in step S104 above, taking every 100-meter segment along the track direction of the ATL08 data product as a unit, all the original ground photon points within each 100-meter segment in the ATL03 data product are converted to the horizontal datum plane, which specifically includes the following steps: Step S1041: Obtain the along-track distance and elevation value of each original ground photon point within the target 100-meter segment in the ATL03 data product; where the target 100-meter segment represents any one of all the 100-meter segments in the ATL03 data product.

[0047] Step S1042: Fit the along-track distance and elevation value of the original ground photon points within the target 100-meter segment to obtain the terrain slope of the target 100-meter segment.

[0048] As can be seen from the above description, in the embodiment of the present invention, taking every 100-meter segment along the track direction of the ATL08 data product as a unit, the original ground photon points within each 100-meter segment in the ATL03 data product are processed. Specifically, first, obtain the along-track distance and elevation value of each original ground photon point within the target 100-meter segment in the ATL03 data product, and then use linear regression to fit the along-track distance and elevation value of the original ground photon points to estimate the terrain slope along the track direction of the target 100-meter segment.

[0049] Step S1043: Construct a two-dimensional rotation matrix based on the terrain slope.

[0050] In the embodiment of the present invention, the terrain slope is the tangent value of the slope angle , that is, . Therefore, after obtaining the terrain slope, rotate around the X-axis in the slope direction to construct the two-dimensional rotation matrix: .

[0051] Step S1044: Based on the two-dimensional rotation matrix, convert each original ground photon point within the target 100-meter segment to the horizontal datum plane to obtain the corresponding ground photon points after removing the terrain.

[0052] For each original ground photon point within the target 100-meter segment, applying the two-dimensional rotation matrix transformation can obtain the ground photon points after removing the slope influence (that is, removing the terrain). The above transformation process is expressed as: ; where represents the coordinate of the original ground photon point, represents the along-track distance of the original ground photon point, represents the elevation of the original ground photon point, represents the coordinate of the ground photon point after removing the terrain, represents the along-track distance of the ground photon point after removing the terrain, Represents the elevation of the ground photon points after removing the terrain.

[0053] After the above steps, all the original ground photon points within the target 100-meter segment will be transformed to a horizontal datum plane, eliminating the slope influence, such that the elevation frequency histogram and the photon distribution characteristics can more accurately reflect the vertical aggregation degree of the ground photon points.

[0054] In an alternative embodiment, in step S1042 above, fitting the along-track distance and elevation values of the original ground photon points within the target 100-meter segment to obtain the terrain slope of the target 100-meter segment specifically includes the following steps: Step S10421, constructing a univariate linear regression equation for the elevation and along-track distance of the original ground photon points; wherein, the elevation of the original ground photon points is the dependent variable of the univariate linear regression equation, the along-track distance of the original ground photon points is the independent variable of the univariate linear regression equation, and the terrain slope of the segment to which the original ground photon points belong is the slope of the univariate linear regression equation.

[0055] That is, the univariate linear regression equation is expressed as: , where represents the coordinates of the original ground photon points, represents the terrain slope of the segment to which the original ground photon points belong, represents the intercept.

[0056] Step S10422, substituting the along-track distance and elevation values of the original ground photon points within the target 100-meter segment into the univariate linear regression equation to obtain a target system of equations.

[0057] Step S10423, solving the target system of equations using the least squares method to obtain the terrain slope of the target 100-meter segment.

[0058] Specifically, the least squares method (Least Squares Method, LSM) is used to calculate the regression parameters and of the univariate linear regression equation to minimize the sum of squared errors, that is: , where represents the total number of the original ground photon points within the target 100-meter segment.

[0059] Taking the derivatives of a and b and setting the derivatives to zero, the optimal solutions can be obtained: .

[0060] In an optional implementation, the above step S106, based on all ground photon points after removing the terrain in each 100-meter segment, constructs an elevation frequency histogram for each 100-meter segment, and calculates the skewness of the elevation frequency histogram and the vertical concentration of ground photons for each 100-meter segment, specifically includes the following steps: Step S1061, group all ground photon points after removing the terrain in the target 100-meter segment according to the preset elevation interval, and count the number of ground photon points falling into each elevation interval; wherein the target 100-meter segment represents any segment of all 100-meter segments in the ATL03 data product.

[0061] Step S1062, with the elevation interval as the horizontal axis and the ratio of the number of ground photon points falling into the elevation interval to the total number of ground photon points in the 100m segment to which it belongs as the vertical axis, construct an elevation frequency histogram of the target 100-meter segment.

[0062] Specifically, all ground photon points after removing the terrain in the target 100-meter segment are grouped at fixed elevation intervals (exemplarily 0.1m), and the number of photons in each elevation interval is counted. The proportion of the number of ground photon points in each elevation interval to the total number of ground photon points in the 100m segment to which it belongs is calculated, and then an elevation frequency histogram is constructed. The height of the bars in the elevation frequency histogram corresponds to the proportion of the number of ground photon points in each elevation interval to the total number of ground photon points, reflecting the frequency of occurrence of the elevation.

[0063] Step S1063, determine the peak of the elevation frequency histogram of the target 100-meter segment. That is, determine the frequency value where the ground photon points are most concentrated in the segment. .

[0064] Step S1064, calculate the standard deviation of the heights of all ground photon points after removing the terrain within the target 100-meter segment.

[0065] Specifically, the height standard deviation is calculated using the following formula: Calculation: ,in, It represents the average elevation of all ground photon points after removing the terrain within the target 100-meter segment, that is, .

[0066] Step S1065, based on the height standard deviation and the elevations of all ground photon points after removing the terrain within the target 100-meter segment, the skewness of the elevation frequency histogram of the target 100-meter segment is calculated.

[0067] The calculation formula of skewness is: .

[0068] Step S1066, calculate the ratio of the peak value to the height standard deviation, and use the ratio result as the vertical concentration of ground photons in the target 100-meter segment.

[0069] In order to quantify the vertical distribution concentration of ground photon points, the embodiment of the present invention proposes to use the ratio of the wave crest to the elevation standard deviation as the ground photon vertical concentration, which is expressed as: Obviously, the vertical concentration of ground photons can characterize the degree of concentration of ground photon points. The larger the value of the vertical concentration of ground photons, the more concentrated the ground photon points are at a certain elevation and the flatter the terrain surface. The smaller the value of the vertical concentration of ground photons, the more dispersed the distribution of ground photon points is, which may be greatly affected by multiple scattering, low vegetation or local undulations.

[0070] In an optional implementation, the above step S108, based on the attribute parameters of each potential elevation control point in the ATL08 data product, the ground photon vertical concentration of the 100-meter segment to which it belongs, and the skewness of the corresponding elevation frequency histogram, screens the potential elevation control points in the ATL08 data product, specifically comprising the following steps: Step S1081, calculating the difference between the elevation of the target potential elevation control point and its reference elevation; wherein the target potential elevation control point represents any potential elevation control point in the ATL08 data product.

[0071] Step S1082: When it is determined that the target potential elevation control point satisfies any of the following preset conditions, the target potential elevation control point is eliminated.

[0072] Among them, the preset conditions include: the absolute value of the difference is greater than the first threshold, the cloud cover identifier is greater than the second threshold, the surface type identifier is equal to the third threshold, the number of ground photon points in the 100-meter segment is less than or equal to the fourth threshold, the total number of photon points in the 100-meter segment is greater than or equal to the fifth threshold, the proportion of ground photon points in the 100-meter segment is less than or equal to the sixth threshold, the point discrete distribution identifier is equal to 1, the sum of the ground photon distribution identifiers of the 100-meter segment is not equal to 5, the vertical concentration of ground photons in the 100-meter segment is less than the seventh threshold, and the absolute value of the skewness of the elevation frequency histogram of the 100-meter segment is greater than the eighth threshold.

[0073] Based on the above steps, it can be seen that in order to improve the accuracy of the elevation control point screening results, the embodiment of the present invention adopts a multiple screening strategy to eliminate potential elevation control points affected by abnormal elevation, cloud layer, water body, ground object reflection, solar radiation, atmospheric interference and vertical dispersion of ground photons. The embodiment of the present invention does not specifically limit the values ​​of the first to eighth thresholds, and users can configure them according to actual conditions.

[0074] Exemplarily, the screening scheme of elevation control points is as follows: Step 1. Remove outlier points: Calculate the difference between the elevation of the target potential elevation control point (h_te_best_fit) and its reference elevation (dem_h), and remove outliers with an absolute value greater than 50m.

[0075] Step 2. Eliminate potential elevation control points affected by clouds: Based on the cloud coverage reflected by the cloud flag cloud_flag_atm, eliminate ground points with cloud_flag_atm greater than 2 to avoid high-density clouds or misclassified points interfering with the accuracy of terrain extraction.

[0076] Step 3. Remove water interference: In the water area, since the laser echo is mainly specular reflection, the pulse expansion is minimal and the energy loss is small. The instrument will detect multiple surface echoes, so that the received photon signals are concentrated at lower elevations, which in turn affects the accuracy of terrain extraction. Therefore, the embodiment of the present invention removes the potential elevation control points with the surface type identifier segment_landcover=80 (water body) to reduce the impact of the water surface.

[0077] Step 4. Eliminate ground reflection, solar radiation and atmospheric influence: Use the number of ground photon points (n_te_photons), the total number of photon points (n_seg_ph) and the ratio of ground photon points to signal photon number (te_radio) for screening, and retain potential elevation control points that meet the following conditions: ①n_te_photons>40. Based on this condition, in order to ensure sufficient ground photon signals, the target potential elevation control points whose number of ground photon points in the 100-meter segment is less than or equal to the fourth threshold are eliminated as described above.

[0078] ②n_seg_ph<500, based on this condition, the abnormally high-density photon areas are removed, corresponding to the above removal of the target potential elevation control points whose total photon points in the 100-meter segment are greater than or equal to the fifth threshold.

[0079] ③te_radio>0.5, based on this condition to ensure the effective surface photon ratio, where, , Indicates the number of ground photon points contained in the 100m segment. Indicates the number of canopy photon points contained in the 100m segment, Indicates the number of photons at the top of the canopy contained in the 100m segment. This inequality corresponds to the above elimination of target potential elevation control points whose ground photon point ratio in the 100m segment is less than or equal to the sixth threshold.

[0080] Step 5. Eliminate potential elevation control points with large vertical geolocation errors: The point discrete distribution identifier psf_flag indicates the total vertical geolocation error caused by ranging errors and local terrain slopes. psf_flag has two values: 0 and 1, indicating that the total vertical geolocation error of the potential elevation control point is less than or equal to 1 and greater than 1, respectively. In the embodiments of the present invention, potential elevation control points with psf_flag = 1 are eliminated to improve the accuracy of elevation control points.

[0081] Step 6. Identify and eliminate potential elevation control points with uneven ground photon distribution or large surface fitting errors: The ground photon distribution identifier subset_te_flag reflects the source situation of the original ground photons used to calculate the 100-meter segmented statistical data. The values of subset_te_flag, -1, 0, and 1, indicate no photon points, no ground photon points, and ground photon points in each 20-meter ATL03 segment, respectively. If the sum of subset_te_flag is 5, it means that all 5 20-meter ATL03 segments contain ground photon points, ensuring the spatial uniformity of the data, thereby reducing the influence caused by uneven photon distribution or surface fitting errors. In the embodiments of the present invention, only potential elevation control points with the sum of subset_te_flag equal to 5 are retained.

[0082] Step 7. Eliminate potential elevation control points affected by multiple scattering, complex terrain, and low vegetation confusion: Use and to further screen potential elevation control points, eliminate abnormal points with large vertical discreteness, and retain potential elevation control points that meet the following conditions: ① , corresponding to the target potential elevation control point in the above text that eliminates the vertical aggregation degree of ground photons in the 100-meter segment to which it belongs is less than the seventh threshold.

[0083] ② , corresponding to the target potential elevation control point in the above text that eliminates the absolute value of the skewness of the elevation frequency histogram in the 100-meter segment to which it belongs is greater than the eighth threshold.

[0084] In summary, the present invention combines two data products, ATL03 and ATL08, comprehensively considers the influence of factors such as laser characteristics, clouds, water bodies, ground object reflections, solar radiation, atmosphere, and the vertical aggregation degree of ground photons on the elevation accuracy of potential elevation control points, and proposes the above-mentioned multiple screening strategies for high-precision elevation control points, thereby effectively improving the extraction accuracy of elevation control points, enhancing the reliability of the data, and providing technical support for global elevation control point extraction and fine terrain mapping.

[0085] To verify the functionality of the method provided in the embodiments of the present invention, the following experiments were carried out. Given that the accuracy requirements for elevation control points are different in flat, hilly and mountainous regions (0.5m, 0.7m and 1.5m respectively), before the experiment, potential elevation control points were divided into flat (slope < 2°), hilly (2° ≤ slope < 6°) and mountainous (6° ≤ slope < 25°) regions according to the terrain slope, and their accuracies were verified respectively.

[0086] The accuracy of the elevation control points extracted by the present invention was verified using the DEM products extracted from airborne lidar data. Figure 2 is the frequency distribution histogram of the ground elevation error between the original data and the elevation control points in the flat region, Figure 3 is the frequency distribution histogram of the ground elevation error between the original data and the elevation control points in the hilly region, Figure 4 is the frequency distribution histogram of the ground elevation error between the original data and the elevation control points in the mountainous region. Figures 2 to 4 The results show that through the screening of elevation control points, the proportion of potential elevation control points in the low error range has increased significantly, while the proportion of potential elevation control points with high errors has decreased significantly, especially in the hilly and mountainous regions.

[0087] Figure 5 shows the comparison results of the accuracy verification of the elevation control points for the ATL08 original data, without adding new parameters (vertical aggregation degree and skewness of ground photons), and after adding new parameters. According to Figure 5 It can be seen that in flat, hilly and mountainous regions, the proportions of ground points where the original data meet the accuracy requirements are only 83.92%, 68.54% and 48.5% respectively; the corresponding elevation errors RMSE are 1.69m, 3.35m and 6m respectively, all exceeding the accuracy requirements of elevation control points. After the screening in steps 1-6 of the above screening scheme, the proportions of potential elevation control points meeting the accuracy requirements are increased to 93.5%, 92.88% and 94.23% respectively; the elevation errors RMSE are reduced to 0.3m, 0.43m and 0.81m respectively. After adding new parameters (step 7), the proportions of potential elevation control points meeting the accuracy requirements are further increased to 94.47%, 96.44% and 97.01% respectively; the elevation errors RMSE are reduced to 0.27m, 0.33m and 0.56m respectively. The results show that the elevation control point screening method provided in the embodiments of the present invention can significantly improve the extraction accuracy of elevation control points, effectively eliminate potential elevation control points that do not meet the accuracy conditions, and the addition of new parameters further improves the accuracy, especially in mountainous regions.

[0088] Furthermore, compared with the elevation control point screening scheme proposed by Li et al. (2021) (Li, B; Xie H; Tong X; Tang H; Wang X. A method of extracting high-accuracy elevation control points from ICESat-2 altimetry data. Photogrammetric Engineering & Remote Sensing. 2021, 87(11): 821-830), the results are as Figure 6 shown. Compared with the method of Li et al. (2021), in the embodiments of the present invention, by adding new screening parameters (vertical aggregation degree and skewness of ground photons), potential elevation control points with large discrete distributions of elevation are effectively removed. In flat, hilly, and mountainous areas, the proportions of elevation control points meeting the accuracy requirements are increased by 1.05%, 2.95%, and 2.67% respectively, and the RMSE errors are reduced by 0.04 m, 0.09 m, and 0.21 m respectively.

[0089] Embodiment 2 The embodiments of the present invention also provide a device for screening elevation control points based on ICESat-2 data. This device is mainly used to execute the method for screening elevation control points based on ICESat-2 data provided in the above Embodiment 1. The following is a specific introduction to the device for screening elevation control points based on ICESat-2 data provided in the embodiments of the present invention.

[0090] Figure 7 is a functional module diagram of a device for screening elevation control points based on ICESat-2 data provided in the embodiments of the present invention. As Figure 7 shown, this device mainly includes: an extraction module 10, a conversion module 20, a calculation module 30, and a screening module 40, where: The extraction module 10 is used to combine the ATL03 data product and the ATL08 data product of ICESat-2 to extract all original ground photon points from the ATL03 data product.

[0091] The conversion module 20 is used to take every 100-meter segment along the track direction of the ATL08 data product as a unit, and convert all the original ground photon points within each 100-meter segment in the ATL03 data product to the horizontal datum plane to obtain all the ground photon points after removing the terrain within each 100-meter segment.

[0092] A calculation module 30 is configured to construct an elevation frequency histogram for each 100-meter segment based on all ground photon points after terrain removal within each 100-meter segment, and calculate the skewness of the elevation frequency histogram and the vertical aggregation degree of ground photons for each 100-meter segment. Wherein, the horizontal axis of the elevation frequency histogram represents the elevation interval, and the vertical axis represents the proportion of the number of ground photon points falling within the elevation interval to the total number of ground photon points within its affiliated 100-meter segment. The vertical aggregation degree of ground photons represents the ratio of the peak of the elevation frequency histogram to the height standard deviation of the ground photon points within the 100-meter segment.

[0093] A screening module 40 is configured to screen potential elevation control points in the ATL08 data product based on the attribute parameters of each potential elevation control point in the ATL08 data product, the vertical aggregation degree of ground photons of its affiliated 100-meter segment, and the skewness of the corresponding elevation frequency histogram, to obtain an elevation control point screening result. Wherein, the potential elevation control point refers to the ATL08 ground photon points used for elevation control screening.

[0094] An embodiment of the present invention provides a device for screening elevation control points based on ICESat-2 data. First, the ATL03 data product and the ATL08 data product of ICESat-2 are combined to extract all original ground photon points from the ATL03 data product. Then, taking each 100-meter segment along the track direction of the ATL08 data product as a unit, all original ground photon points within each 100-meter segment in the ATL03 data product are converted to the horizontal reference plane to obtain all ground photon points after terrain removal within each 100-meter segment. Furthermore, an elevation frequency histogram that can accurately reflect the photon aggregation degree is constructed for each 100-meter segment, and the skewness of the elevation frequency histogram and the vertical aggregation degree of ground photons for each 100-meter segment are calculated. Finally, based on the attribute parameters of each potential elevation control point in the ATL08 data product, the vertical aggregation degree of ground photons of its affiliated 100-meter segment, and the skewness of the corresponding elevation frequency histogram, the potential elevation control points in the ATL08 data product are screened to obtain an elevation control point screening result. The embodiment of the present invention proposes two new parameters for screening elevation control points by combining the ATL03 and ATL08 data products: the vertical aggregation degree of ground photons and the skewness. Compared with screening only using ATL08 parameters, it can effectively improve the extraction accuracy of elevation control points and provide technical support for global elevation control point extraction and fine terrain mapping.

[0095] Optionally, the extraction module 10 is specifically configured to: Extract the attribute parameters of each photon point in the ATL03 data product and the ATL08 data product.

[0096] Perform spatial matching on the ATL03 data product and the ATL08 data product based on the attribute parameters.

[0097] Use the classification information of the target photon points in the ATL08 data product as the classification label of the photon points in the ATL03 data product that are spatially matched to them; wherein, the target photon points represent any one of all the photon points in the ATL08 data product.

[0098] Determine all the original ground photon points based on the classification labels of all the photon points in the ATL03 data product.

[0099] Optionally, the conversion module 20 includes: An acquisition unit for acquiring the along-track distance and elevation value of each original ground photon point within a target 100-meter segment in the ATL03 data product; wherein, the target 100-meter segment represents any one of all the 100-meter segments in the ATL03 data product.

[0100] A fitting unit for fitting the along-track distance and elevation value of the original ground photon points within the target 100-meter segment to obtain the terrain slope of the target 100-meter segment.

[0101] A construction unit for constructing a two-dimensional rotation matrix based on the terrain slope.

[0102] A conversion unit for converting each original ground photon point within the target 100-meter segment to a horizontal reference plane based on the two-dimensional rotation matrix to obtain the corresponding ground photon points after terrain removal.

[0103] Optionally, the fitting unit is specifically configured to: Construct a unary linear regression equation for the elevation and along-track distance of the original ground photon points; wherein, the elevation of the original ground photon points is the dependent variable of the unary linear regression equation, the along-track distance of the original ground photon points is the independent variable of the unary linear regression equation, and the terrain slope of the segment to which the original ground photon points belong is the slope of the unary linear regression equation.

[0104] Substitute the along-track distance and elevation value of the original ground photon points within the target 100-meter segment into the unary linear regression equation to obtain a target system of equations.

[0105] Solve the target system of equations using the least squares method to obtain the terrain slope of the target 100-meter segment.

[0106] Optionally, the calculation module 30 is specifically configured to: Group all the ground photon points after terrain removal within the target 100-meter segment at a preset elevation interval and count the number of ground photon points falling into each elevation interval; wherein, the target 100-meter segment represents any one of all the 100-meter segments in the ATL03 data product.

[0107] Taking the elevation interval as the horizontal axis and the ratio of the number of ground photon points falling within the elevation interval to the total number of ground photon points within its corresponding 100m segment as the vertical axis, construct the elevation frequency histogram for the target 100m segment.

[0108] Determine the peak of the elevation frequency histogram for the target 100m segment.

[0109] Calculate the height standard deviation of all ground photon points after terrain removal within the target 100m segment.

[0110] Based on the height standard deviation and the elevations of all ground photon points after terrain removal within the target 100m segment, calculate the skewness of the elevation frequency histogram for the target 100m segment.

[0111] Calculate the ratio of the peak to the height standard deviation, and use the ratio result as the vertical aggregation degree of ground photons for the target 100m segment.

[0112] Optionally, the screening module 40 is specifically configured to: Calculate the difference between the elevation of the target potential elevation control point and its reference elevation; wherein, the target potential elevation control point represents any potential elevation control point in the ATL08 data product.

[0113] When it is determined that the target potential elevation control point meets any of the following preset conditions, eliminate the target potential elevation control point.

[0114] Among them, the preset conditions include: the absolute value of the difference is greater than the first threshold, the cloud cover identifier is greater than the second threshold, the surface type identifier is equal to the third threshold, the number of ground photon points within its corresponding 100m segment is less than or equal to the fourth threshold, the total number of photon points within its corresponding 100m segment is greater than or equal to the fifth threshold, the ratio of ground photon points within its corresponding 100m segment is less than or equal to the sixth threshold, the point discrete distribution identifier is equal to 1, the sum of the ground photon distribution identifiers for its corresponding 100m segment is not equal to 5, the vertical aggregation degree of ground photons for its corresponding 100m segment is less than the seventh threshold, and the absolute value of the skewness of the elevation frequency histogram for its corresponding 100m segment is greater than the eighth threshold.

[0115] Embodiment III See Figure 8 , this embodiment of the present invention provides an electronic device, which includes: a processor 60, a memory 61, a bus 62, and a communication interface 63, and the processor 60, the communication interface 63, and the memory 61 are connected through the bus 62; the processor 60 is used to execute an executable module stored in the memory 61, such as a computer program.

[0116] Among them, the memory 61 may include high-speed random access memory (RAM), or may also include non-volatile memory, such as at least one disk memory. The communication connection between this system network element and at least one other network element is realized through at least one communication interface 63 (which can be wired or wireless), and the Internet, wide area network, local area network, metropolitan area network, etc. can be used.

[0117] The bus 62 can be an ISA bus, a PCI bus, an EISA bus, etc. The bus can be divided into an address bus, a data bus, a control bus, etc. For the sake of easy representation, Figure 8 only a bidirectional arrow is used in the figure, but it does not mean that there is only one bus or one type of bus.

[0118] Among them, the memory 61 is used to store a program. After receiving an execution instruction, the processor 60 executes the program. The method executed by the device defined by any embodiment of the foregoing embodiments of the present invention can be applied to or implemented by the processor 60.

[0119] The processor 60 may be an integrated circuit chip with signal processing capabilities. In the implementation process, the steps of the above method can be completed by the integrated logic circuit in the hardware of the processor 60 or the instructions in the form of software. The above-mentioned processor 60 can be a general-purpose processor, including a central processing unit (CPU for short), a network processor (NP for short), etc.; it can also be a digital signal processor (DSP for short), an application specific integrated circuit (ASIC for short), a field-programmable gate array (FPGA for short), or other programmable logic devices, discrete gate or transistor logic devices, discrete hardware components. It can implement or execute the various methods, steps and logic block diagrams disclosed in the embodiments of the present invention. The general-purpose processor can be a microprocessor or the processor can also be any conventional processor, etc. The steps of the method disclosed in combination with the embodiments of the present invention can be directly embodied as being executed and completed by a hardware decoding processor, or executed and completed by a combination of hardware and software modules in the decoding processor. The software module can be located in a mature storage medium in the art such as random access memory, flash memory, read-only memory, programmable read-only memory, or electrically erasable programmable memory, register, etc. This storage medium is located in the memory 61, and the processor 60 reads the information in the memory 61 and combines its hardware to complete the steps of the above method.

[0120] A computer program product of a method and device for screening elevation control points based on ICESat-2 data provided by an embodiment of the present invention includes a computer-readable storage medium storing non-volatile program code executable by a processor. The instructions included in the program code can be used to execute the method described in the foregoing method embodiment. For specific implementation, reference can be made to the method embodiment and will not be elaborated here.

[0121] In addition, in each embodiment of the present invention, each functional unit can be integrated in a processing unit, or each unit can exist physically alone, or two or more units can be integrated in one unit.

[0122] If the above-mentioned function is implemented in the form of a software functional unit and sold or used as an independent product, it can be stored in a non-volatile computer-readable storage medium executable by a processor. Based on such understanding, the technical solution of the present invention, in essence, or the part that contributes to the prior art, or a part of this technical solution, can be embodied in the form of a software product. This computer software product is stored in a storage medium and includes several instructions to enable a computer device (which can be a personal computer, a server, or a network device, etc.) to execute all or part of the steps of the method described in each embodiment of the present invention. The foregoing storage medium includes: various media such as USB flash drives, mobile hard disks, read-only memories (ROM, Read-Only Memory), random access memories (RAM, Random Access Memory), magnetic disks, or optical discs that can store program code.

[0123] It should be noted that similar reference numerals and letters indicate similar items in the following drawings. Therefore, once an item is defined in one drawing, it does not need to be further defined and explained in subsequent drawings.

[0124] In the description of the present invention, it should be noted that the orientation or positional relationship indicated by the terms "center", "upper", "lower", "left", "right", "vertical", "horizontal", "inner", "outer", etc. is based on the orientation or positional relationship shown in the drawings, or the orientation or positional relationship in which the product of the present invention is usually placed during use. It is only for the convenience of describing the present invention and simplifying the description, rather than indicating or implying that the device or element referred to must have a specific orientation, be constructed and operated in a specific orientation, and therefore cannot be understood as a limitation to the present invention. In addition, the terms "first", "second", "third", etc. are only used for distinguishing descriptions and cannot be understood as indicating or implying relative importance.

[0125] In addition, terms such as "horizontal", "vertical", "hanging", etc. do not mean that the components are required to be absolutely horizontal or hanging, but can be slightly inclined. For example, "horizontal" only means that its direction is more horizontal relative to "vertical", and does not mean that the structure must be completely horizontal, but can be slightly inclined.

[0126] In the description of the present invention, it should also be noted that, unless otherwise clearly specified and defined, the terms "arranged", "installed", "connected", and "coupled" should be understood in a broad sense. For example, it can be a fixed connection, a detachable connection, or an integral connection; it can be a mechanical connection or an electrical connection; it can be directly connected or indirectly connected through an intermediate medium, and it can be the communication inside two components. For those of ordinary skill in the art, the specific meanings of the above terms in the present invention can be understood according to specific situations.

[0127] Finally, it should be noted that: the above embodiments are only used to illustrate the technical solutions of the present invention, rather than to limit it; although the present invention has been described in detail with reference to the foregoing embodiments, those of ordinary skill in the art should understand that: they can still modify the technical solutions recorded in the foregoing embodiments, or perform equivalent replacements on some or all of the technical features; and these modifications or replacements do not make the essence of the corresponding technical solutions deviate from the scope of the technical solutions of the embodiments of the present invention.

Claims

1. A method for screening elevation control points based on ICESat-2 data, characterized in that include: Combine the ATL03 data product and the ATL08 data product of ICESat-2 to extract all raw ground photon points from the ATL03 data product; Taking each 100-meter segment of the ATL08 data product in the along-track direction as a unit, all original ground photon points in each 100-meter segment of the ATL03 data product are converted to a horizontal reference plane to obtain all ground photon points after removing the terrain in each 100-meter segment; Based on all the ground photon points after removing the terrain in each 100-meter segment, an elevation frequency histogram of each 100-meter segment is constructed, and the skewness of the elevation frequency histogram and the vertical concentration of ground photons in each 100-meter segment are calculated; wherein the horizontal axis of the elevation frequency histogram represents the elevation interval, and the vertical axis represents the proportion of the number of ground photon points falling into the elevation interval to the total number of ground photon points in the 100-meter segment to which it belongs; the vertical concentration of ground photons represents the ratio of the peak of the elevation frequency histogram to the height standard deviation of the ground photon points in the 100-meter segment; Based on the attribute parameters of each potential elevation control point in the ATL08 data product, the vertical concentration of ground photons in the 100-meter segment to which it belongs, and the skewness of the corresponding elevation frequency histogram, the potential elevation control points in the ATL08 data product are screened to obtain the elevation control point screening result; wherein the potential elevation control point represents the ATL08 ground photon point used for elevation control screening.

2. The method for screening elevation control points based on ICESat-2 data according to claim 1, wherein Combine ICESat-2's ATL03 data product with the ATL08 data product to extract all raw ground photon points from the ATL03 data product, including: Extracting attribute parameters of each photon point in the ATL03 data product and the ATL08 data product; spatially matching the ATL03 data product and the ATL08 data product based on the attribute parameters; The classification information of the target photon point in the ATL08 data product is used as the classification label of the photon point in the ATL03 data product that matches the space thereof; wherein the target photon point represents any photon point among all the photon points in the ATL08 data product; Based on the classification labels of all photon points in the ATL03 data product, all original ground photon points are determined.

3. The method for screening elevation control points based on ICESat-2 data according to claim 1, wherein Taking the ATL08 data product as a unit of each 100-meter segment along the track, all original ground photon points in each 100-meter segment of the ATL03 data product are converted to a horizontal reference plane, including: Obtaining the along-track distance and elevation value of each original ground photon point in a target 100-meter segment in the ATL03 data product; wherein the target 100-meter segment represents any segment of all 100-meter segments in the ATL03 data product; Fitting the along-track distance and elevation value of the original ground photon point within the target 100-meter segment to obtain the terrain slope of the target 100-meter segment; constructing a two-dimensional rotation matrix based on the terrain slope; Based on the two-dimensional rotation matrix, each original ground photon point within the 100-meter segment of the target is converted to a horizontal reference plane to obtain a corresponding ground photon point after removing the terrain.

4. The method for screening elevation control points based on ICESat-2 data according to claim 3, wherein The along-track distance and elevation value of the original ground photon point in the target 100-meter segment are fitted to obtain the terrain slope of the target 100-meter segment, including: Constructing a univariate linear regression equation of the elevation and along-track distance of the original ground photon point; wherein the elevation of the original ground photon point is the dependent variable of the univariate linear regression equation, the along-track distance of the original ground photon point is the independent variable of the univariate linear regression equation, and the terrain slope of the segment to which the original ground photon point belongs is the slope of the univariate linear regression equation; Substituting the along-track distance and elevation value of the original ground photon point within the target 100-meter segment into the univariate linear regression equation to obtain a target equation group; The target equation group is solved by the least square method to obtain the terrain slope of the target 100-meter segment.

5. The method for screening elevation control points based on ICESat-2 data according to claim 1, wherein Based on all ground photon points after removing the terrain in each 100-meter segment, an elevation frequency histogram of each 100-meter segment is constructed, and the skewness of the elevation frequency histogram and the vertical concentration of ground photons in each 100-meter segment are calculated, including: Group all ground photon points after removing terrain in the target 100-meter segment at preset elevation intervals, and count the number of ground photon points falling into each elevation interval; wherein the target 100-meter segment represents any segment of all 100-meter segments in the ATL03 data product; With the elevation interval as the horizontal axis and the ratio of the number of ground photon points falling into the elevation interval to the total number of ground photon points in the 100m segment to which it belongs as the vertical axis, an elevation frequency histogram of the target 100m segment is constructed; Determine the peak of the elevation frequency histogram of the target 100-meter segment; Calculate the height standard deviation of all ground photon points after removing the terrain within the 100-meter segment of the target; Calculating the skewness of the elevation frequency histogram of the target 100-meter segment based on the height standard deviation and the elevations of all ground photon points after removing the terrain within the target 100-meter segment; The ratio of the wave peak to the height standard deviation is calculated to use the ratio result as the ground photon vertical concentration of the target 100-meter segment.

6. The method for screening elevation control points based on ICESat-2 data according to claim 1, wherein Based on the attribute parameters of each potential elevation control point in the ATL08 data product, the vertical concentration of ground photons in the 100-meter segment to which it belongs, and the skewness of the corresponding elevation frequency histogram, the potential elevation control points in the ATL08 data product are screened, including: Calculate the difference between the elevation of a target potential elevation control point and its reference elevation; wherein the target potential elevation control point represents any potential elevation control point in the ATL08 data product; When it is determined that the target potential elevation control point meets any of the following preset conditions, the target potential elevation control point is eliminated; Among them, the preset conditions include: the absolute value of the difference is greater than the first threshold, the cloud cover identifier is greater than the second threshold, the surface type identifier is equal to the third threshold, the number of ground photon points in the 100-meter segment is less than or equal to the fourth threshold, the total number of photon points in the 100-meter segment is greater than or equal to the fifth threshold, the proportion of ground photon points in the 100-meter segment is less than or equal to the sixth threshold, the point discrete distribution identifier is equal to 1, the sum of the ground photon distribution identifiers of the 100-meter segment is not equal to 5, the vertical concentration of ground photons in the 100-meter segment is less than the seventh threshold, and the absolute value of the skewness of the elevation frequency histogram of the 100-meter segment is greater than the eighth threshold.

7. An apparatus for screening elevation control points based on ICESat-2 data, characterized in that, include: An extraction module, for combining the ATL03 data product and the ATL08 data product of ICESat-2 to extract all original ground photon points from the ATL03 data product; A conversion module is used to convert all original ground photon points in each 100-meter segment of the ATL03 data product to a horizontal reference plane, taking each 100-meter segment of the ATL08 data product in the along-track direction as a unit, to obtain all ground photon points after removing the terrain in each 100-meter segment; A calculation module is used to construct an elevation frequency histogram for each 100-meter segment based on all ground photon points after removing the terrain in each 100-meter segment, and calculate the skewness of the elevation frequency histogram and the vertical concentration of ground photons in each 100-meter segment; wherein the horizontal axis of the elevation frequency histogram represents the elevation interval, and the vertical axis represents the proportion of the number of ground photon points falling into the elevation interval to the total number of ground photon points in the 100-meter segment to which it belongs; the vertical concentration of ground photons represents the ratio of the peak of the elevation frequency histogram to the height standard deviation of the ground photon points in the 100-meter segment; A screening module is used to screen the potential elevation control points in the ATL08 data product based on the attribute parameters of each potential elevation control point in the ATL08 data product, the vertical concentration of ground photons in the 100-meter segment to which it belongs, and the skewness of the corresponding elevation frequency histogram, so as to obtain the elevation control point screening result; wherein the potential elevation control point represents the ATL08 ground photon point used for elevation control screening.

8. The device for screening elevation control points based on ICESat-2 data according to claim 7, wherein The extraction module is specifically used for: Extracting attribute parameters of each photon point in the ATL03 data product and the ATL08 data product; spatially matching the ATL03 data product and the ATL08 data product based on the attribute parameters; The classification information of the target photon point in the ATL08 data product is used as the classification label of the photon point in the ATL03 data product that matches the space thereof; wherein the target photon point represents any photon point among all the photon points in the ATL08 data product; Based on the classification labels of all photon points in the ATL03 data product, all original ground photon points are determined.

9. An electronic device, comprising a memory and a processor, wherein a computer program capable of running on the processor is stored on the memory, and is characterized in that When the processor executes the computer program, the method for screening elevation control points based on ICESat-2 data according to any one of claims 1 to 6 is implemented.

10. A computer-readable storage medium, characterized in that, The computer-readable storage medium stores computer instructions, and when the computer instructions are executed by a processor, the method for screening elevation control points based on ICESat-2 data according to any one of claims 1 to 6 is implemented.

Images