Laser elevation control point screening method considering vegetation coverage and terrain distribution

Abstract

The invention discloses a laser elevation control point screening method considering vegetation coverage and terrain distribution. The method comprises the following steps: extracting land elevation data from an ICESat 2ATL08 product, and preprocessing the land elevation data; preliminarily screening the preprocessed land elevation data according to influence factors of laser elevation quality, and extracting 100-meter unit interval data meeting quality requirements; in each 100-meter unit interval data meeting the quality requirement, extracting the midpoint position of 20-meter fitting elevation values which meet the quality requirement and are uniformly distributed under different terrains in each 100-meter unit interval by utilizing an elevation quality screening condition, and calculating the height of each 100-meter unit interval according to the 100-meter fitting elevation values and the 20-meter fitting elevation values which meet the quality requirement. Screening a final fitting elevation value of the midpoint position according to an elevation error evaluation result; all finally reserved fitting elevation values serve as elevation control points to be imported into a data set.

Description

Technical Field

[0001] This invention relates to the technical field of aerospace photogrammetry, and more particularly to a method for selecting laser elevation control points that takes into account vegetation cover and terrain distribution. Background Technology

[0002] Elevation control points, as fundamental data in surveying and mapping and geographic information, play an irreplaceable role in the production of remote sensing digital products and the construction of global geographic information resources. In recent years, global mapping missions have significantly increased the demand for high-precision, high-coverage elevation control points, especially in areas with complex terrain and difficult field control operations. This places higher demands on the accuracy, density, and production efficiency of elevation control points. Traditional elevation control point measurement work is inefficient, time-consuming, and costly, and is limited by factors such as climate and terrain, making it difficult to meet the needs of large-scale, accurate mapping. In contrast, spaceborne laser altimetry technology, with its characteristics of low atmospheric interference, strong penetration, high accuracy, and global coverage, is gradually becoming a key technology for acquiring high-precision control points in stereoscopic mapping.

[0003] Chinese invention patent CN112985358A discloses a method and system for extracting global elevation control points from ICESat 2 / ATLAS. This includes three steps: preprocessing, coarse screening, and fine screening of ICESat 2 ATL08 data products. Laser points with gentle terrain and low cloud cover are retained as ground elevation control points, and the data is finally imported into the ICESat 2 global elevation control point database. Chinese invention patent CN116467293 discloses a method and system for filtering global elevation control points based on ICESat 2 data. This includes three steps: extracting land elevation parameters from the ICESat 2 ATL08 product and performing preprocessing, initial screening, and a joint evaluation system. The evaluation results of each parameter of the obtained control points are weighted, fused, and classified. Finally, the classified elevation control point data is imported into the global elevation control point database.

[0004] The two methods described above are based on the ATL08 data product from the ICESat 2 satellite. They analyze the 100-meter fitted elevation values ​​using the data quality labels provided by ATL08 and can automatically extract high-precision global elevation control point data after hierarchical classification. However, the error factors considered in the above ATL08 data screening process are limited. They do not take into account the impact of surface type factors such as vegetation cover on the accuracy of laser elevation. Furthermore, the extracted elevation control points are only at a ground spatial resolution of 100 meters, which is insufficient to show small-scale terrain undulations. The sparse number and distribution of control points are also insufficient to show more detailed terrain features. Summary of the Invention

[0005] This invention provides a laser elevation control point selection method that takes into account vegetation cover and terrain distribution. It can automatically extract high-precision global elevation control point data with a spatial resolution of 20 meters, solving the technical problem that the spatial resolution of elevation control points extracted using ICESat 2 is low and cannot show the small-scale terrain undulation features. This provides data support for stereo mapping and product quality inspection.

[0006] This invention can be achieved through the following technical solutions:

[0007] A method for selecting laser elevation control points that considers vegetation cover and topographic distribution includes the following steps:

[0008] Step 1: Extract land elevation data from the ICESat 2ATL08 product and perform preprocessing;

[0009] Step 2: Based on the factors affecting the quality of laser elevation data, perform initial screening on the preprocessed land elevation data and extract data in 100-meter unit intervals that meet the quality requirements.

[0010] Step 3: Within each 100-meter unit interval of data that meets the quality requirements, use the elevation quality screening conditions to extract the 20-meter fitted elevation values ​​that are evenly distributed under different terrains and meet the quality requirements.

[0011] Step 4: At the midpoint of each 100-meter unit interval, combine the 100-meter fitted elevation value and the 20-meter fitted elevation value that meet the quality requirements, and select the final fitted elevation value for that midpoint position based on the elevation error assessment results.

[0012] Step 5: Import all the fitted elevation values ​​that are ultimately retained into the dataset as elevation control points.

[0013] Furthermore, in step one, through preprocessing, the retained quality labels include: 100-meter fitted elevation value longitude label, 100-meter fitted elevation value latitude label, 100-meter fitted elevation value h_te_best_fit label, 20-meter fitted elevation value longitude_20m label, 100-meter fitted elevation value latitude_20m label, 20-meter fitted elevation value h_te_best_fit_20m label, 100-meter unit interval best reference DEM elevation value dem_h label, 100-meter unit interval land cover type segment_landcover label, and 100-meter single... The tags for vegetation cover (segment_cover), terrain slope (terrain_slope), terrain inclination (h_te_skew), cloud confidence (cloud_flag_atm), total photon count (n_seg_ph), terrain photon count (n_te_photons), terrain photon distribution (subset_te_flag), median elevation (h_te_median), and elevation interpolation (h_te_interp) are specified for each unit of measurement.

[0014] Furthermore, in step three, the elevation quality screening conditions are set as gross error factor and median elevation difference;

[0015] To address gross errors, when filtering data using the 20-meter fitted elevation value h_te_best_fit_20m label against the best reference DEM elevation value dem_h for a 100-meter unit interval, 20-meter fitted elevation values ​​h_te_best_fit_20m labels with an elevation difference ΔH2 greater than the threshold 3σ are removed. The elevation difference ΔH2 is calculated using the following equation.

[0016] ΔH2=|h_te_best_fit_20m-dem_h|>3σ

[0017] In the formula, σ represents the absolute elevation accuracy threshold of the reference DEM;

[0018] For the median elevation difference, when filtering the 20-meter fitted elevation value h_te_best_fit_20m using the median elevation of the 100-meter unit interval (h_te_median), 20-meter fitted elevation values ​​h_te_best_fit_20m with an elevation difference ΔH3 greater than the threshold TMedian are removed. The elevation difference ΔH3 is calculated using the following equation.

[0019] △H3=|h_te_best_fit_20m-h_te_median|>TMedianⅠ

[0020] In the formula, the threshold TMedianⅠ represents the dynamic threshold based on the change in terrain slope within the unit interval.

[0021] Furthermore, the threshold TMedianⅠ is set to 0.5 when the terrain_slope_angle value of a 100-meter unit interval is between 0 and 2;

[0022] When the value of terrain_slope_angle in a 100-meter unit interval is between 2 and 6, the threshold TMedianⅠ is set to 0.7;

[0023] When the value of terrain_slope_angle in a 100-meter unit interval is greater than 6, the threshold TMedianⅠ is set to 1.5;

[0024] The value of terrain_slope_angle is the angle value converted from the radian value corresponding to the terrain_slope label for a 100-meter unit interval.

[0025] Furthermore, in step four, at the midpoint of each 100-meter unit interval, different filtering strategies are adopted for different situations at the midpoint location to select the final fitted elevation value corresponding to that midpoint location. The filtering strategies are as follows:

[0026] If only a 100-meter fitted elevation value h_te_best_fit that meets the quality requirements exists, and no 20-meter fitted elevation value h_te_best_fit_20m that meets the quality requirements exists, then the 100-meter fitted elevation value h_te_best_fit is used as the final elevation value for that location and is retained.

[0027] When both a 100-meter fitted elevation value h_te_best_fit and a 20-meter fitted elevation value h_te_best_fit_20m that meet the quality requirements exist, the elevation error of the two fitted elevation values ​​is evaluated using the 100-meter unit interval elevation median label h_te_median and the 100-meter unit interval elevation interpolation label h_te_interp. The fitted elevation value with the smaller error evaluation result is selected as the final elevation value for that midpoint position and retained.

[0028] Furthermore, the elevation error index E(h) used for elevation error assessment is calculated using the following formula.

[0029] E(h)=|h-h_te_median|*w1+|h-h_te_interp|*w2

[0030] In the formula, h represents two fitted elevation values ​​for elevation error assessment, and w1 and w2 are error weights set based on different terrain slopes.

[0031] Furthermore, the error weights w1 and w2 are set to 0.8 and 0.2 respectively when the terrain_slope_angle value of the 100-meter unit interval is between 0 and 2.

[0032] When the terrain_slope_angle value for a 100-meter unit interval is between 2 and 6, w1 is set to 0.6 and w2 is set to 0.4;

[0033] When the `terrain_slope_angle` value is greater than 6 for a 100-meter unit interval, w1 is set to 0.2 and w2 is set to 0.8.

[0034] The value of terrain_slope_angle is the angle value converted from the radian value corresponding to the terrain_slope label for a 100-meter unit interval.

[0035] Furthermore, after filtering at the midpoint, the 100-meter fitted elevation values ​​h_te_best_fit are further filtered using the 100-meter unit interval elevation median label h_te_median. The following 100-meter fitted elevation values ​​with elevation differences greater than the threshold TMedianⅡ are removed:

[0036] △H4=|h_te_best_fit-h_te_median|>TMedianⅡ

[0037] In the formula, the threshold TMedianⅡ represents the dynamic threshold based on the change in terrain slope within a 100-meter unit interval.

[0038] Furthermore, in step two, the influencing factors include gross error factors, surface factors, topographic factors, atmospheric factors, and photon quality factors.

[0039] Regarding gross errors, when filtering data for the optimal reference DEM elevation value dem_h label within a 100-meter unit interval, 100-meter unit intervals with an elevation difference ΔH1 greater than a threshold are excluded. The formula for calculating the elevation difference ΔH1 is as follows.

[0040] △H1=|h_te_best_fit-dem_h|>3σ

[0041] In the formula, σ represents the absolute elevation accuracy threshold of the reference DEM;

[0042] For surface factors, the segment_landcover label of the 100-meter unit interval is used to filter the vegetation coverage segment_cover label of the 100-meter unit interval, and the 100-meter unit interval with the segment_cover label above the vegetation coverage threshold TCover is removed.

[0043] For terrain factors, when using the terrain slope (terrain_slope) label to filter the terrain inclination (h_te_skew) label within a 100-meter unit interval, the terrain_slope label is first converted from radian values ​​to angle values ​​(terrain_slope_angle). Then, 100-meter unit intervals with an angle value (terrain_slope_angle) above the terrain slope threshold (TSlope) are removed, as are 100-meter unit intervals with an h_te_skew label above the inclination threshold (TSkew).

[0044] For atmospheric factors, when filtering data for the cloud_flag_atm label of cloud amount confidence level in a 100-meter unit interval, remove 100-meter unit intervals where the cloud_flag_atm label is below the cloud amount threshold Tcloud.

[0045] Regarding photon quality factors, when filtering data on the topographic photon count (n_te_photons) within a 100-meter unit interval using the total photon count (n_seg_ph) and the topographic photon distribution (subset_te_flag) within the same 100-meter unit interval, the following criteria are excluded: 100-meter unit intervals with n_seg_ph exceeding the photon count threshold (TPh); 100-meter unit intervals with the sum of data within the subset_te_flag group below the topographic photon distribution threshold (TFlag); and 100-meter unit intervals with a topographic photon proportion (ΔP) less than the threshold (TPhotos). The formula for calculating the topographic photon proportion (ΔP) is as follows:

[0046] △P=|n_te_photons / n_seg_ph| <TPhotos

[0047] In the formula, the threshold TPhotos represents the terrain photon ratio threshold.

[0048] The beneficial technical effects of this invention are as follows:

[0049] 1. First, all beam data is read from the ICESat 2ATL08 product track and preprocessed, retaining only the quality labels used for screening. This ensures both the amount of data used for screening and avoids redundancy in data quality labels, saving processing time. Second, the data is initially screened based on five factors affecting laser elevation quality, such as surface vegetation cover, extracting only 100-meter unit intervals that meet the requirements. Then, for the extracted 100-meter unit intervals, dynamic quality screening based on terrain slope is performed on the 20-meter fitted elevation values, extracting 20-meter fitted elevation values ​​that meet the elevation quality requirements, thus reducing the minimum interval of the final elevation control points from 100 meters to 20 meters. Finally, at the midpoint of each unit interval, the final elevation value at that location is selected by combining the 100-meter and 20-meter fitted elevation values ​​that meet the quality requirements under different terrain conditions and the elevation error assessment results of the two. This ensures that all fitted elevations in the data are fully considered, and the fitted elevation with the best quality is used as the final elevation control point.

[0050] 2. Compared with existing technologies, this invention has developed a more comprehensive screening strategy to address the impact of various error factors on ICESat 2 satellites, such as gross errors, terrain, surface conditions, atmosphere, and photon quality. This results in more robust extracted elevation control points under different environments. Furthermore, the use of 20-meter fitted elevation values ​​to extract elevation control points results in a higher density, meeting the requirements for elevation control point density in remote sensing data production at different resolutions and map sizes. The optimal selection of fitted elevations at overlapping locations further improves the accuracy of elevation control points, solving the problems of sparse control points and unstable elevation accuracy in previous elevation control point screening systems. This provides denser and more accurate control point data support for the production and inspection of higher-resolution remote sensing mapping products. Attached Figure Description

[0051] Figure 1 This is a schematic diagram of the overall process of the present invention. Detailed Implementation

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

[0053] like Figure 1As shown, this invention provides a method for selecting laser elevation control points that considers vegetation cover and topographic distribution. Land elevation parameters are extracted from the ICESat 2ATL08 product and preprocessed, retaining only the quality labels used in data selection to reduce data volume. Then, initial screening is performed using five types of factors influencing laser elevation quality, such as surface vegetation, to extract 100-meter unit intervals that meet quality requirements. Within each compliant unit interval, 20-meter fitted elevation values ​​with balanced distribution under different terrain conditions are extracted. At the midpoint of each unit interval, the 100-meter and 20-meter fitted elevation values ​​that meet quality requirements are combined, and the final elevation value at that location is selected based on the elevation error assessment results. All compliant fitted elevation values ​​are stored as results in the elevation control point dataset, providing a large-scale, high-precision elevation control point dataset with 20-meter intervals for the remote sensing field.

[0054] The details are as follows:

[0055] Step 1: Data Preprocessing

[0056] When traversing all beam data in the ICESat 2ATL08 product track, only the quality labels used for screening are retained to reduce data redundancy.

[0057] In this example, the data quality labels that need to be retained include: longitude label for 100-meter fitted elevation values, latitude label for 100-meter fitted elevation values, h_te_best_fit label for 100-meter fitted elevation values, longitude_20m label for 20-meter fitted elevation values, latitude_20m label for 100-meter fitted elevation values, h_te_best_fit_20m label for 20-meter fitted elevation values, dem_h label for the best reference DEM elevation value for the 100-meter unit interval, segment_landcover label for the 100-meter unit interval, and vegetation label for the 100-meter unit interval. The tags are: segment_cover (cover), terrain slope (100-meter interval), terrain skew (100-meter interval), cloud confidence (100-meter interval), total number of photons (100-meter interval), number of terrain photons (100-meter interval), terrain photon distribution (100-meter interval), median elevation (100-meter interval), and elevation interpolation (100-meter interval).

[0058] Step 2: Data extraction within a 100-meter unit interval:

[0059] All the ICESat 2ATL08 data obtained from preprocessing were used as a unit interval of 100 meters. The data were initially screened by five categories of factors affecting laser elevation quality: gross error, surface, topography, atmosphere, and photons. Only the 100-meter unit intervals that met the requirements were extracted.

[0060] The details are as follows:

[0061] 1. Gross error factor screening includes data screening of the best reference DEM elevation value dem_h label for a 100-meter unit interval, and removing unit intervals with elevation difference ΔH1=|h_te_best_fit-dem_h|>3σ, where σ represents the absolute elevation accuracy threshold of the reference DEM, which is an empirically given value such as 10.

[0062] 2. Surface factor filtering includes data filtering using the 100-meter unit interval land cover type label (segment_landcover) and the 100-meter unit interval vegetation cover label (segment_cover). Unit intervals within the land cover type set TLand (available in ATL08) are removed. Considering that the elevation control points extracted in this invention only target conventional land, the default Tland set is set to {70, 80, 200}, representing snow and ice, inland water bodies, and ocean, respectively.

[0063] Meanwhile, considering that areas with excessive vegetation cover may affect the elevation accuracy of terrain photons, the segment_cover unit intervals above the vegetation cover threshold TCover are excluded, and the default TCover is set to 65.

[0064] 3. Topographic factor filtering includes data filtering using the terrain slope (terrain_slope) and terrain inclination (h_te_skew) labels for 100-meter unit intervals. First, terrain_slope is converted from radian values ​​to angle values ​​(terrain_slope_angle). Considering that areas with significant terrain undulations can affect the elevation accuracy of control points, unit intervals with terrain_slope_angle exceeding the terrain slope threshold (TSlope), with a default setting of TSlope to 25, are removed. Similarly, unit intervals with h_te_skew exceeding the inclination threshold (TSkew), with a default setting of TSkew to 1, are also removed.

[0065] 4. Atmospheric factor screening includes using the cloud_flag_atm label of cloud amount confidence level in a 100-meter unit interval to filter data, and removing unit intervals where cloud_flag_atm is below the cloud amount threshold Tcloud. Based on previous invention experience, Tcloud is set to 2 by default.

[0066] 5. The photon mass factor involves data screening using the total number of photons n_seg_ph label in a 100-meter unit interval, the number of terrain photons n_te_photons label in a 100-meter unit interval, and the terrain photon distribution subset_te_flag label in a 100-meter unit interval. Unit intervals with n_seg_ph above the photon number threshold TPh are excluded, and TPh is default set to 500; since the quality of terrain photons has a crucial impact on the accuracy of the fitted elevation, to ensure the uniform distribution of terrain photons within a unit interval, unit intervals with the sum of the within-group data of subset_te_flag below the terrain photon distribution threshold TFlag are excluded, and TFlag is default set to 5; at the same time, to ensure that there is a sufficient proportion of terrain photons in the unit interval for elevation fitting, unit intervals with the terrain photon ratio △P = |n_te_photons / n_seg_ph| < TPhotos are excluded, where TPhotos represents the terrain photon ratio threshold, and is default set to 0.5.

[0067] Step 3: Extraction of 20-meter fitted elevation data

[0068] Within each 100-meter unit interval that meets the requirements, using the elevation quality screening conditions, the quality screening of all 20-meter fitted elevation values is continued, and 20-meter fitted elevation values with balanced distribution under different terrains are extracted and retained. Among them, the elevation quality screening conditions are set as the gross error factor and the median elevation difference.

[0069] 1. For the gross error factor

[0070] The gross error factor screening of the 20-meter fitted elevation value h_te_best_fit_20m is carried out using the best reference DEM elevation value dem_h label in a 100-meter unit interval, and 20-meter fitted elevation values with the elevation difference ΔH2 = |h_te_best_fit_20m - dem_h| > 3σ are excluded, where σ represents the absolute elevation accuracy threshold of the reference DEM, and the empirically given value is such as 10.

[0071] 2. For the median elevation difference

[0072] The median elevation difference screening of the 20-meter fitted elevation value h_te_best_fit_20m is carried out using the elevation median h_te_median label in a 100-meter unit interval, and 20-meter fitted elevation values with the elevation difference △H3 = |h_te_best_fit_20m - h_te_median| > TMedianⅠ are excluded, where the threshold TMedianⅠ represents the median elevation difference threshold.

[0073] Considering that the accuracy of elevation values ​​decreases as the terrain slope increases, in order to ensure that there are enough usable elevation control points under different terrain slopes, the threshold TMedianⅠ is set as a dynamic threshold that changes according to the terrain slope where the fitted elevation value is located. Specifically, when the terrain_slope_angle value of a 100-meter unit interval is between 0 and 2, the threshold TMedianⅠ is set to 0.5.

[0074] When the terrain_slope_angle value of a 100-meter unit interval is between 2 and 6, set the threshold TMedianⅠ to 0.7;

[0075] When the value of terrain_slope_angle in a 100-meter unit interval is greater than 6, the threshold TMedianⅠn is set to 1.5.

[0076] For the 20-meter fitted elevation value in step three that meets the elevation quality requirements but is not located at the midpoint of the unit interval, it is retained as the final elevation value for that location.

[0077] Step 4: Fitted Elevation Fusion Extraction

[0078] At the midpoint of each 100-meter interval, combining the 100-meter and 20-meter fitted elevation values ​​that meet the quality requirements, and based on the elevation error assessment results, the final fitted elevation value at that location is selected and retained. The specific process is as follows:

[0079] After obtaining the 100-meter fitted elevation value h_te_best_fit that meets the quality requirements, different filtering strategies are adopted for the following two situations regarding the elevation values ​​at the midpoint of the unit interval:

[0080] i. When there is only a 100-meter fitted elevation value h_te_best_fit that meets the quality requirements, and there is no 20-meter fitted elevation value h_te_best_fit_20m that meets the quality requirements, then the 100-meter fitted elevation value h_te_best_fit shall be used as the final elevation value for that location.

[0081] ii. When both a 100-meter fitted elevation value h_te_best_fit and a 20-meter fitted elevation value h_te_best_fit_20m that meet the quality requirements exist, the elevation error of the two fitted elevation values ​​is evaluated using the 100-meter unit interval elevation median label h_te_median and the 100-meter unit interval elevation interpolation label h_te_interp. The fitted elevation value with the smaller elevation error is selected and retained.

[0082] The elevation error index used for elevation error assessment is calculated using the following formula.

[0083] E(h)=|h-h_te_median|*w1+|h-h_te_interp|*w2

[0084] Where h represents two fitted elevation values ​​used for elevation error assessment, and w1 and w2 are error weights set based on different terrain slopes.

[0085] When the value of terrain_slope_angle in a unit interval is between 0 and 2, w1 = 0.8 and w2 = 0.2; when the value of terrain_slope_angle in a unit interval is between 2 and 6, w1 = 0.6 and w2 = 0.4; when the value of terrain_slope_angle in a unit interval is greater than 6, w1 = 0.2 and w2 = 0.8.

[0086] In addition, to ensure that the 100-meter fitted elevation value and the 20-meter fitted elevation value have the same accuracy standard, the median elevation difference of the 100-meter fitted elevation value h_te_best_fit is also used to filter the median elevation difference using the 100-meter unit interval median elevation value h_te_median. 100-meter fitted elevation values ​​with an elevation difference ΔH4 = |h_te_best_fit - h_te_median| > TMedianⅡ are removed, where TMedianⅡ represents the median elevation difference threshold.

[0087] When the value of terrain_slope_angle in a 100-meter unit interval is between 0 and 2, set the threshold TMedianⅡ to 0.5;

[0088] When the terrain_slope_angle value of a 100-meter unit interval is between 2 and 6, set the threshold TMedianⅡ to 0.7;

[0089] When the value of terrain_slope_angle in a 100-meter unit interval is greater than 6, the threshold TMedianⅡ is set to 1.5.

[0090] The fitted elevation value that remains at the midpoint of the unit interval in step four is retained as the final elevation value for that location.

[0091] Step 5: Import all the final fitted elevation values ​​into the elevation control point dataset as the result.

[0092] To address the issues of sparse control points and unstable elevation accuracy in previous methods of selecting elevation control points due to the influence of various factors such as gross errors, terrain, surface conditions, atmosphere, and photon quality on ICESat 2 satellites, this new method for selecting global elevation control points can extract higher-precision, higher-density 20-meter interval elevation control points from ICESat 2ATL08 global data products. This provides denser and more accurate control point data support for the production and testing of higher-resolution remote sensing mapping products.

[0093] While specific embodiments of the present invention have been described above, those skilled in the art should understand that these are merely illustrative examples. Various changes or modifications can be made to these embodiments without departing from the principles and essence of the present invention. Therefore, the scope of protection of the present invention is defined by the appended claims.

Claims

1. A method for selecting laser elevation control points considering vegetation cover and topographic distribution, characterized in that... Includes the following steps: Step 1: Extract land elevation data from the ICESat 2ATL08 product and perform preprocessing; Step 2: Based on the factors affecting the quality of laser elevation data, perform initial screening on the preprocessed land elevation data and extract data in 100-meter unit intervals that meet the quality requirements. Step 3: Within each 100-meter unit interval of data that meets the quality requirements, use the elevation quality screening conditions to extract the 20-meter fitted elevation values ​​that are evenly distributed under different terrains and meet the quality requirements. Step 4: At the midpoint of each 100-meter unit interval, combine the 100-meter fitted elevation value and the 20-meter fitted elevation value that meet the quality requirements, and select the final fitted elevation value for that midpoint position based on the elevation error assessment results. Step 5: Import all the fitted elevation values ​​that are ultimately retained into the dataset as elevation control points.

2. The laser elevation control point selection method considering vegetation cover and topographic distribution according to claim 1, characterized in that: In step one, through preprocessing, the retained quality labels include: 100-meter fitted elevation value (longitude), 100-meter fitted elevation value (latitude), 100-meter fitted elevation value (h_te_best_fit), 20-meter fitted elevation value (longitude_20m), 100-meter fitted elevation value (latitude_20m), 20-meter fitted elevation value (h_te_best_fit_20m), 100-meter unit interval best reference DEM elevation value (dem_h), 100-meter unit interval land cover type (segment_landcover), and 100-meter unit area... The tags for vegetation cover (segment_cover), terrain slope (terrain_slope), terrain inclination (h_te_skew), cloud confidence (cloud_flag_atm), total number of photons (n_seg_ph), number of terrain photons (n_te_photons), terrain photon distribution (subset_te_flag), median elevation (h_te_median), and elevation interpolation (h_te_interp) are defined as follows:

3. The laser elevation control point selection method considering vegetation cover and topographic distribution according to claim 2, characterized in that: In step three, the elevation quality screening conditions are set as gross error factor and median elevation difference; To address gross errors, when filtering data using the 20-meter fitted elevation value h_te_best_fit_20m label against the best reference DEM elevation value dem_h for a 100-meter unit interval, 20-meter fitted elevation values ​​h_te_best_fit_20m labels with an elevation difference ΔH2 greater than the threshold 3σ are removed. The elevation difference ΔH2 is calculated using the following equation. ΔH2=|h_te_best_fit_20m-dem_h|>3σ In the formula, σ represents the absolute elevation accuracy threshold of the reference DEM; For the median elevation difference, when filtering the 20-meter fitted elevation value h_te_best_fit_20m using the median elevation of the 100-meter unit interval (h_te_median), 20-meter fitted elevation values ​​h_te_best_fit_20m with an elevation difference ΔH3 greater than the threshold TMedian are removed. The elevation difference ΔH3 is calculated using the following equation. △H3=|h_te_best_fit_20m-h_te_median|>TMedianⅠ In the formula, the threshold TMedianⅠ represents the dynamic threshold based on the change in terrain slope within the unit interval.

4. The laser elevation control point selection method considering vegetation cover and topographic distribution according to claim 3, characterized in that: The threshold TMedianⅠ is set to 0.5 when the terrain_slope_angle value of a 100-meter unit interval is between 0 and 2; When the value of terrain_slope_angle in a 100-meter unit interval is between 2 and 6, the threshold TMedianⅠ is set to 0.7; When the value of terrain_slope_angle in a 100-meter unit interval is greater than 6, the threshold TMedianⅠ is set to 1.5; The value of terrain_slope_angle is the angle value converted from the radian value corresponding to the terrain_slope label for a 100-meter unit interval.

5. The laser elevation control point selection method considering vegetation cover and topographic distribution according to claim 1, characterized in that, In step four, at the midpoint of each 100-meter unit interval, different filtering strategies are adopted for different situations at the midpoint, and the final fitted elevation value corresponding to that midpoint is selected. The filtering strategies are as follows: If only a 100-meter fitted elevation value h_te_best_fit that meets the quality requirements exists, and no 20-meter fitted elevation value h_te_best_fit_20m that meets the quality requirements exists, then the 100-meter fitted elevation value h_te_best_fit is used as the final elevation value for that location and is retained. When both a 100-meter fitted elevation value h_te_best_fit and a 20-meter fitted elevation value h_te_best_fit_20m that meet the quality requirements exist, the elevation error of the two fitted elevation values ​​is evaluated using the 100-meter unit interval elevation median label h_te_median and the 100-meter unit interval elevation interpolation label h_te_interp. The fitted elevation value with the smaller error evaluation result is selected as the final elevation value for that midpoint position and retained.

6. The laser elevation control point selection method considering vegetation cover and topographic distribution according to claim 5, characterized in that: The elevation error index E(h) used for elevation error assessment is calculated using the following formula. E(h)=|h-h_te_median|*w1+|h-h_te_interp|*w2 In the formula, h represents two fitted elevation values ​​for elevation error assessment, and w1 and w2 are error weights set based on different terrain slopes.

7. The laser elevation control point selection method considering vegetation cover and topographic distribution according to claim 6, characterized in that: The error weights w1 and w2 are set to 0.8 and 0.2 respectively when the terrain_slope_angle value of the 100-meter unit interval is between 0 and 2. When the terrain_slope_angle value for a 100-meter unit interval is between 2 and 6, w1 is set to 0.6 and w2 is set to 0.4; When the `terrain_slope_angle` value is greater than 6 for a 100-meter unit interval, w1 is set to 0.2 and w2 is set to 0.

8. The value of terrain_slope_angle is the angle value converted from the radian value corresponding to the terrain_slope label for a 100-meter unit interval.

8. The laser elevation control point selection method considering vegetation cover and topographic distribution according to claim 5, characterized in that: After filtering at the midpoint, the 100-meter fitted elevation values ​​h_te_best_fit are further filtered using the 100-meter unit interval elevation median label h_te_median. The following 100-meter fitted elevation values ​​with elevation differences greater than the threshold TMedianⅡ are removed: △H4=|h_te_best_fit-h_te_median|>TMedianⅡ In the formula, the threshold TMedianⅡ represents the dynamic threshold based on the change in terrain slope within a 100-meter unit interval.

9. The laser elevation control point selection method considering vegetation cover and topographic distribution according to claim 2, characterized in that: In step two, the influencing factors include gross error factors, surface factors, topographic factors, atmospheric factors, and photon quality factors. Regarding gross errors, when filtering data for the optimal reference DEM elevation value dem_h label within a 100-meter unit interval, 100-meter unit intervals with an elevation difference ΔH1 greater than a threshold are excluded. The formula for calculating the elevation difference ΔH1 is as follows. △H1=|h_te_best_fit-dem_h|>3σ In the formula, σ represents the absolute elevation accuracy threshold of the reference DEM; For surface factors, the segment_landcover label of the 100-meter unit interval is used to filter the vegetation coverage segment_cover label of the 100-meter unit interval, and the 100-meter unit interval with the segment_cover label above the vegetation coverage threshold TCover is removed. For terrain factors, when using the terrain slope (terrain_slope) label to filter the terrain inclination (h_te_skew) label within a 100-meter unit interval, the terrain_slope label is first converted from radian values ​​to angle values ​​(terrain_slope_angle). Then, 100-meter unit intervals with an angle value (terrain_slope_angle) above the terrain slope threshold (TSlope) are removed, as are 100-meter unit intervals with an h_te_skew label above the inclination threshold (TSkew). For atmospheric factors, when filtering data for the cloud_flag_atm label of cloud amount confidence level in a 100-meter unit interval, remove 100-meter unit intervals where the cloud_flag_atm label is below the cloud amount threshold Tcloud. Regarding photon quality factors, when filtering data on the topographic photon count (n_te_photons) within a 100-meter unit interval using the total photon count (n_seg_ph) and the topographic photon distribution (subset_te_flag) within the same 100-meter unit interval, the following criteria are excluded: 100-meter unit intervals with n_seg_ph exceeding the photon count threshold (TPh); 100-meter unit intervals with the sum of data within the subset_te_flag group below the topographic photon distribution threshold (TFlag); and 100-meter unit intervals with a topographic photon proportion (ΔP) less than the threshold (TPhotos). The formula for calculating the topographic photon proportion (ΔP) is as follows: △P=|n_te_photons / n_seg_ph| <TPhotos In the formula, the threshold TPhotos represents the terrain photon ratio threshold.

Images