A hemispherical-scale snow depth estimation method and system based on spaceborne lidar

Abstract

This invention relates to the field of lidar snow depth estimation technology, specifically providing a hemispherical-scale snow depth estimation method and system based on spaceborne lidar. The method includes: downloading and preprocessing snow cover-related data collected by the spaceborne lidar within a preset hemispherical scale; obtaining snow cover period and snowless period light segments within the preset hemispherical scale based on the preprocessed snow cover-related data; obtaining the intersection points of the snow cover period and snowless period light segments using a dynamic discrimination rule based on a slope-limited intersection point distance threshold; obtaining snow depth values ​​based on the intersection points and a snow depth estimation model; and correcting the snow depth values ​​based on a snow depth reference error model to obtain corrected snow depth values, thus completing the hemispherical-scale snow depth estimation based on spaceborne lidar. This invention is not limited by factors such as altitude or cost when generating snow depth data, and can compensate for sparse station snow depth data. It plays an important role in verifying the authenticity of snow depth products, constructing snow cover models, and measuring snow depth in high-altitude mountainous areas.

Description

Technical Field

[0001] This invention belongs to the field of lidar snow depth estimation technology, specifically relating to a hemispherical scale snow depth estimation method and system based on spaceborne lidar. Background Technology

[0002] Snow cover is the most widely distributed element in the cryosphere. With high albedo and low thermal conductivity, it plays a vital role in surface energy flow and the regulation of regional and global ecosystems. Simultaneously, snowmelt also plays a crucial role in the global water cycle and water balance. Snow depth is one of the most important and fundamental parameters in snow cover research, and it is also a key factor in forecasting, monitoring, and warning of snowmelt-related flooding disasters.

[0003] Currently, long-term snow depth data at global or hemispherical scales mainly come from passive microwave inversion or atmospheric reanalysis data. However, the spatial resolution of these data is low; existing mainstream passive microwave inversion algorithms tend to saturate in brightness-temperature difference when the snow depth reaches about 50 cm; in high-altitude mountainous areas, snow depth inversion involves significant uncertainty, and suitable observational data is lacking for verification. Direct manual measurement or observation using meteorological stations is the most accurate and intuitive way to obtain snow depth data, and this data is widely used for the verification of snow depth products. However, due to limitations imposed by topography, altitude, and natural conditions, observation efficiency is low, and the data is discontinuous and uneven. Especially in remote high-altitude mountainous areas, direct snow depth measurements are very rare, and the spatial heterogeneity of snow cover is particularly significant due to factors such as topography and climate. Therefore, it is crucial to develop a dataset with wide coverage, high spatial resolution, an inversion process unaffected by snow depth thresholds, and the ability to verify snow depth in high-altitude mountainous areas.

[0004] By comparing the elevation differences at the ICESat orbital intersections before and after snow cover, the depth of new snow on the ice sheet can be estimated, and the results have been proven reliable. While ICESat and high-precision DEMs can be used to determine snow depth on small-scale mountain glaciers, ICESat has only one laser beam, resulting in a sparse footprint, making it difficult to apply to snow depth measurements in non-polar regions. ICESat-2 has six laser beams, a spot footprint size of 17m, and an interval of 0.7m between adjacent spots, increasing the density of cross-track sampling and achieving a surface elevation accuracy of 0.1m, enabling large-scale, wide-area dense observation of snow cover. The standard ICESat-2 products ATL06 and ATL08 have the potential to estimate snow depth in mid-to-high latitude regions; most snow depth measurements are performed using ICESat-2 elevations during the snow cover period and accurate digital terrain models (DTMs) during the snowless period. Related studies have shown that the accuracy of ICESat-2 snow depth estimation depends on terrain slope and forest canopy density. However, these methods rely on high-precision ground reference elevations for registration in the validation areas, making large-scale snow depth estimation impossible. Using only the elevation difference between trajectory intersections in the ICESat-2 ATL08 product can accurately estimate snow depth in flat areas, but it introduces larger errors in mountainous regions. Furthermore, the validation used only a small number of intersections, leading to significant randomness. Using the SlideRule Earth project, the ATL06 product was improved, and snow depth was validated across different vegetation covers, elevations, and slopes. The results indicate that ICESat-2 SlideRule may provide a new snow depth dataset spanning the western United States and even the entire globe.

[0005] ICESat-2 has good snow depth estimation capabilities, but the selected areas are limited, the snow cover type is singular, and the method has poor universality, making it only applicable to snow depth estimation under specific background conditions. Summary of the Invention

[0006] To address the problems existing in the prior art, this invention proposes a hemispherical-scale snow depth estimation method and system based on spaceborne lidar. The snow depth value is obtained by subtracting the intersection points of the lidar segments during snowy and snowless periods. This method takes into account the relationship between terrain and the geometric structure of the lidar segments, dynamically determines the distance between intersection points using terrain slope, and corrects for snow depth estimation errors caused by terrain undulations. Verification of the results obtained by this invention shows that this method can provide a new technical approach for acquiring snow depth data.

[0007] A hemispherical-scale snow depth estimation method based on spaceborne lidar, the method comprising:

[0008] Download snow-related data within a preset hemispherical scale collected by the spaceborne lidar, and preprocess the snow-related data;

[0009] Based on the preprocessed snow cover-related data, the light segments during the snow cover period and the light segments during the snowless period are obtained within a preset hemispherical scale;

[0010] A dynamic discrimination rule based on a slope-limited intersection distance threshold is adopted to obtain the intersection points of the light segment during the snowy period and the light segment during the snowless period; wherein, the dynamic discrimination rule based on a slope-limited intersection distance threshold includes: within a preset hemispherical scale, setting a threshold for judging the center distance between intersection points based on a preset terrain slope threshold;

[0011] Based on the intersection point and the snow depth estimation model, the snow depth value within a preset hemispherical scale is obtained; wherein, the snow depth estimation model is the difference between the elevation values ​​of the light segment during the snow-covered period and the light segment during the snowless period at the intersection point;

[0012] Based on the terrain slope of the light segment during the snowless period and the horizontal distance between the light segment during the snow-covered period and the light segment during the snowless period when the intersection is formed, a snow depth reference error model is constructed.

[0013] The snow depth value is corrected based on the snow depth reference error model to obtain the corrected snow depth value, thus completing the hemispherical scale snow depth estimation based on spaceborne lidar.

[0014] Preferably, the snow-related data includes location data, elevation data, slope data, land use type data, and snow cover data within a preset hemispherical scale;

[0015] The method for preprocessing the snow cover-related data includes:

[0016] Remove light segments with a slope greater than 5°, as well as light segments of water bodies, buildings, and ice surfaces from the snow-related data.

[0017] Preferably, the snow cover period includes: for the period between 30°N and 55°N, the months with snow cover are defined as December, January, February, and March, totaling 4 months; and for the period between 55°N and 70°N, the months with snow cover are defined as December, January, February, March, and April, totaling 5 months.

[0018] Preferably, the method for obtaining the snow-free period of the light segment includes:

[0019] Import the MOD10A1 snow cover product into the GEE software and set the NDSI of the MOD10A1 snow cover product to 0.2. When the NDSI of the MOD10A1 snow cover product is less than 0.2, it is determined that there is no snow cover; otherwise, it is determined that there is snow cover.

[0020] Import the snow cover map into ArcGIS software and divide it into a 5°×5° grid to obtain the snow cover situation of each 5° cell within the preset date. Change the date to obtain the first and last days of snow cover for each cell.

[0021] Based on the first day of snow accumulation and the last day of snow accumulation, the snow-free period of the light segment is obtained, i.e., the snow-free period.

[0022] The snowless period is compared with the corresponding time period's segment_landcover field in the preprocessed snow cover data. If the segment_landcover field in the preprocessed snow cover data is 1, the snowless period time period is correctly identified, and the snowless period light segment is obtained. If the segment_landcover field is 2, the snowless period time period is incorrectly identified, and the segment_landcover field comparison is repeated on a different snowless date until the segment_landcover field is 1, thus obtaining the snowless period light segment. Here, segment_landcover is a snow cover identifier; a segment_landcover field of 1 represents snowless land, and a segment_landcover field of 2 represents land with snow.

[0023] Preferably, the dynamic discrimination rule for the intersection distance threshold based on slope specifically includes:

[0024] The preset terrain slope threshold is 2°. When the terrain slope is greater than or equal to 2°, the threshold for determining the center distance between intersections is set to 2m; when the terrain slope is less than 2°, the threshold for determining the center distance between intersections is set to 17m.

[0025] Preferably, when calculating the snow depth value, the elevation difference between the light segments during the snow cover period and the snowless period at the intersection is screened. When the elevation difference is negative or greater than the maximum snow depth in the current area, the corresponding light segment is removed. When the segment_landcover field of the light segment during the snow cover period and the snowless period are both 0, 1 and 3, the elevation difference is not calculated and the snow depth value is directly set to 0. Among them, a segment_landcover field of 0 represents water body and a segment_landcover field of 3 represents ice surface.

[0026] Preferably, the method for obtaining the horizontal distance between the light segment during the snow-covered period and the light segment during the snowless period when the intersection point is formed includes:

[0027] Convert the latitude and longitude corresponding to the intersection of the light segment during the snow-covered period and the light segment during the snowless period into radians, and calculate the difference between the radians of longitude and latitude respectively.

[0028] Based on the aforementioned radian difference, the initial distance between the light segment during the snow cover period and the light segment during the snowless period is obtained using the Haversine formula;

[0029] Calculate the great circle angle based on the initial distance between the light segment during the snow-covered period and the light segment during the snowless period;

[0030] Based on the great circle angle and the Earth's radius, the horizontal distance between the light segment during the snowy season and the light segment during the snowless season is obtained.

[0031] This invention also provides a hemispherical-scale snow depth estimation system based on spaceborne lidar for implementing the method, comprising:

[0032] The data acquisition module is used to download snow-related data within a preset hemispherical scale collected by the spaceborne lidar, and to preprocess the snow-related data.

[0033] The light segment segmentation module is used to obtain the light segments during the snow-covered period and the light segments during the snowless period within a preset hemispherical scale based on the preprocessed snow-related data.

[0034] The intersection point acquisition module is used to obtain the intersection points of the light segment during the snowy season and the light segment during the snowless season by adopting a dynamic discrimination rule based on the intersection point distance threshold limited by the slope; wherein, the dynamic discrimination rule based on the intersection point distance threshold limited by the slope includes: setting a threshold for judging the center distance between intersection points based on a preset terrain slope threshold within a preset hemispherical scale;

[0035] The snow depth calculation module is used to obtain the snow depth value within a preset hemispherical scale based on the intersection point and the snow depth estimation model; wherein, the snow depth estimation model is the difference between the elevation values ​​of the light segment during the snow-covered period and the light segment during the snowless period at the intersection point;

[0036] The error model construction module is used to construct a snow depth reference error model based on the terrain slope of the light segment during the snowless period and the horizontal distance between the light segment during the snow-covered period and the light segment during the snowless period when the intersection is formed.

[0037] The error correction module is used to correct the snow depth value based on the snow depth reference error model, obtain the corrected snow depth value, and complete the hemispherical scale snow depth estimation based on the spaceborne lidar.

[0038] Compared with the prior art, the beneficial effects of the present invention are as follows: The present invention utilizes ICESat-2 spaceborne lidar data, and through the relationship between terrain and optical segment geometry, combined with dynamic discrimination of intersection distance thresholds limited by slope, calculates the difference between the optical segment elevations during snowy and snowless periods to obtain snow depth values. Furthermore, it uses methods such as elevation difference filtering and terrain geometry correction to reduce snow depth estimation errors, and finally obtains a seasonal snow depth dataset for the Northern Hemisphere from 2018 to 2020.

[0039] This invention provides a new approach and method for snow depth data acquisition. This method is not limited by factors such as altitude and cost when generating snow depth data. It can make up for sparse snow depth data from stations and plays an important role in verifying the authenticity of snow depth products, constructing snow cover models, and measuring snow depth in high-altitude mountainous areas. Attached Figure Description

[0040] To more clearly illustrate the technical solution of the present invention, the drawings used in the embodiments are briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0041] Figure 1 This is a schematic diagram of the ICESat-2 satellite altimetry process according to an embodiment of the present invention;

[0042] Figure 2 This is a schematic diagram of the processing flow of the ICESat-2 snow depth estimation method according to an embodiment of the present invention;

[0043] Figure 3 This is a schematic diagram of the ICESat-2 trajectory intersection formation method according to an embodiment of the present invention; (a) is the intersection formed by repeating tracks of different periods; (b) is the intersection formed by different tracks;

[0044] Figure 4 This is a schematic diagram of snow depth reference error according to an embodiment of the present invention;

[0045] Figure 5 This is a statistical chart showing the classification of snow depth estimation deviations in an embodiment of the present invention.

[0046] Figure 6 This is a statistical chart showing the classification of snow depth estimation deviations for different snow-covered months according to an embodiment of the present invention.

[0047] Figure 7 This is a scatter plot showing the deviation between the snow depth at the site and the ICESat-2 snow depth estimation in an embodiment of the present invention.

[0048] Figure 8 This is a statistical chart showing the error levels of snow depth estimation in different ranges according to an embodiment of the present invention.

[0049] Figure 9 This is a graph showing the relationship between terrain slope and snow depth estimation errors in an embodiment of the present invention (data set with slope ≤ 5°).

[0050] Figure 10 This is a diagram showing the relationship between terrain slope and snow depth estimation errors in an embodiment of the present invention.

[0051] Figure 11 This is a graph showing the relationship between altitude and snow depth estimation errors in an embodiment of the present invention.

[0052] Figure 12 This is a graph showing the relationship between the distance between intersections and the snow depth estimation deviation in an embodiment of the present invention;

[0053] Figure 13This is a graph showing the relationship between the distance between stations and intersections and the snow depth estimation deviation in an embodiment of the present invention;

[0054] Figure 14 This is a graph showing the relationship between snowless months and snow depth estimation deviation in an embodiment of the present invention;

[0055] Figure 15 This is a heatmap showing the factors affecting the relative error of ICESat-2 snow depth estimation in an embodiment of the present invention. Detailed Implementation

[0056] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0057] To make the above-mentioned objects, features and advantages of the present invention more apparent and understandable, the present invention will be further described in detail below with reference to the accompanying drawings and specific embodiments.

[0058] The software platform and products used in the embodiments of this invention are described as follows:

[0059] Google Earth Engine (GEE) is a cloud-based geospatial analytics platform designed to provide global-scale geospatial data storage, processing, and analysis services. GEE boasts a petabyte-scale catalog of publicly available and free geospatial datasets, allowing users instant access to years of historical data and the latest real-time data. Leveraging cloud infrastructure, GEE automatically subdivides and allocates computational tasks, enabling massively parallel processing and supporting high-throughput analysis. The platform integrates powerful visualization tools, supporting real-time viewing of analysis results for intuitive user understanding. GEE provides JavaScript and Python APIs, allowing users to quickly implement complex geospatial analyses using familiar programming languages. Application scenarios include environmental monitoring, climate change research, and geospatial applications.

[0060] MOD10A1 is one of the MODIS (Moderate Resolution Imaging Spectroradiometer) satellite remote sensing data products, belonging to the snow cover dataset, providing daily global snow cover information at a 500-meter resolution from the Terra satellite. The data includes: Snow cover, expressed as a percentage of the snow cover in each pixel; Snow albedo, providing albedo information for the snow surface; Quality Assessment (QA) data, containing quality assessment information for the snow cover data to help users determine its reliability; and the Normalized Difference Snow Index (NDSI), used to identify snow-covered areas. MOD10A1 products are widely used in meteorology, climate research, water resource management, and environmental monitoring.

[0061] ArcGIS is a powerful Geographic Information System (GIS) software widely used in various fields such as map creation, spatial analysis, data management, and geographic information visualization.

[0062] Example 1

[0063] like Figure 1 , Figure 2 As shown, a hemispherical-scale snow depth estimation method based on spaceborne lidar is presented, which includes:

[0064] S1: Batch download snow-related data within a preset hemispherical scale collected by the spaceborne lidar, and preprocess the snow-related data; batch download data, and read lat, lon, slope, h_te_best_fit(H fit The information includes landcover, snowcover, date, etc., where lat represents latitude, lon represents longitude, slope represents slope, date represents date, and the index H... fit This indicates the best-fit elevation of the light segment in 100-meter segments, landcover indicates the land use type, snowcover indicates whether there is snow, and date indicates the date of snow cover.

[0065] When downloading and reading snow-related data in batches, the 30°N-70°N area is divided into blocks of 5°×5° and processed one by one. In a further implementation, the snow-related data includes location data, elevation data, slope data, land use type data, and snow cover data within a preset hemispherical scale. In this embodiment, the snow depth proxy index is also constructed using the snow-related data. The snow depth proxy index includes snow depth data and snow depth reference error data.

[0066] Methods for preprocessing snow cover-related data include:

[0067] Remove light segments with a slope greater than 5° from snow-related data, as well as light segments from water bodies, buildings, and ice surfaces.

[0068] S2: Based on the preprocessed snow cover data, obtain the light segments during the snow cover period and the light segments during the snowless period within the preset hemispherical scale.

[0069] A further implementation method is that the snow cover period includes: the snow cover months in the 30°N-55°N range are determined to be 4 months, namely December, January, February, and March; and the snow cover months in the 55°N-70°N range are determined to be 5 months, namely December, January, February, March, and April.

[0070] A further implementation method for obtaining snowless light segments includes:

[0071] To make the snowless period closer to the first or last day of snow cover, this invention uses the MOD10A1 product on the GEE platform to determine the snowless period. First, binarization is performed using the NDSI band with a threshold of 0.2 to obtain data on whether a region has snow or no snow on a specific date. Specifically, the MOD10A1 snow cover product is imported into the GEE software, and the NDSI band threshold of the MOD10A1 snow cover product is set to 0.2. When the NDSI band threshold of the MOD10A1 snow cover product is less than 0.2, it is determined to be snowless; otherwise, it is determined to have snow cover, thus obtaining a binarized snow cover map.

[0072] The binarized snow cover map is imported into ArcGIS software and divided into 5°×5° grids to obtain the snow cover situation of each 5° cell within a preset date. The date is changed to obtain the first and last days of snow cover for each cell. In this embodiment, the first or last day of snow cover is pushed forward or backward by 3-5 days to obtain the snowless period, i.e., the snowless period, and the total number of snowless days is guaranteed to be 40 days.

[0073] Based on the first day of snow cover and the last day of snow cover, the light segment during the snowless period is obtained.

[0074] The snowless period is compared with the corresponding time period's segment_landcover field in the preprocessed snow cover data. If the segment_landcover field in the preprocessed snow cover data is 1, the snowless period time period is correctly identified, and the snowless period light segment is obtained. If the segment_landcover field is 2, the snowless period time period is incorrectly identified. The segment_landcover field comparison is repeated on a different snowless date until the segment_landcover field is 1, thus obtaining the snowless period light segment. Here, segment_landcover is the snow cover identifier; a segment_landcover field of 1 represents snowless land, and a segment_landcover field of 2 represents land with snow.

[0075] Repeat the above cross-comparison process until the segment_landcover field of MOD10A1 image and ICESat-2 both show no snow on the same day and within the same processing unit. Then, the day can be determined as the earliest snow-free date to appear in the frame. This method can make the time interval between the snow-free period and the snow-covered period shorter.

[0076] S3: Employing a dynamic discrimination rule based on a slope-defined threshold for intersection point distance, the intersection points of the light segments during the snowy season and the light segments during the snowless season are obtained. The dynamic discrimination rule based on a slope-defined threshold for intersection point distance includes: setting a threshold for judging the center distance between intersection points within a preset hemispherical scale, based on a preset terrain slope threshold. For example... Figure 3 As shown, Figure 3 (a) Represents the intersection of repeating orbits with different periods: the satellite returns to this point after a complete revisit cycle, and the two trajectories appear to be parallel but slightly offset; Figure 3 (b) Represents intersections formed by different orbits: satellites may not have completed a full revisit cycle, and the intersections may have been formed at any time.

[0077] A further implementation method is that the dynamic discrimination rule for the intersection distance threshold based on slope limitation specifically includes:

[0078] In this embodiment, there are two ways in which ICESat-2 footprint intersections are formed. The first is through intersections formed by repeating orbits with different periods, where the two orbits appear to be parallel. This method results in a large number of dense and continuous intersections. The second is through intersections formed by different orbits, resulting in fewer and sparser intersections. In low- and mid-latitude regions, satellite orbits are sparse, and revisits will cause orbital shifts, resulting in mostly discrete intersections formed by different orbits.

[0079] The footprint intersections are not perfectly coincident. Since the diameter of the ICESat-2 light spot footprint is 17m, a distance between the center points of the footprints within the range of 0-17m can be considered an intersection. When the distance is 0, the two light segments are considered to be completely coincident. Although a smaller distance indicates that the two light segments represent the same geographical location, it results in too few intersections meeting the criteria. Setting the distance threshold to 17m causes too much error in areas with significant terrain undulations. Figure 4 As shown.

[0080] This invention dynamically selects a distance threshold based on the relationship between snow depth error and slope. When the slope is greater than 5 degrees, the snow depth error exceeds 30cm, therefore this invention only considers areas with a slope of less than 5 degrees. The preset terrain slope threshold is 2°. When the slope is greater than or equal to 2°, it is considered a gentle slope (2-5°), and the terrain will have significant undulations. To reduce the error, the threshold for determining the center distance between intersections is set to 2m. When the terrain slope is less than 2°, the terrain is relatively flat, and the threshold for determining the center distance between intersections is set to 17m.

[0081] The dynamic discrimination rule for the intersection distance threshold limited by slope can generate more intersections under the condition that the snow depth estimation error is small.

[0082] S4: Based on the intersection point and snow depth estimation model, obtain the snow depth value within the preset hemispherical scale; where the snow depth estimation model is the difference between the intersection point elevation values ​​of the light segments during the snow cover period and the snowless period; in the ICESat-2ATL08 data, the "dem_h" field represents the elevation value of a specific location extracted from the external digital elevation model (DEM). Since this value is constant over time, it cannot be used to calculate the daily snow depth by subtracting the intersection points; "h_te_best_fit"(H fit The H field represents the best-fit elevation for the midpoint of each 100-meter segment. This segment elevation is obtained by interpolation of the midpoint of the 100-meter segment using a selected best-fit method (including linear, cubic, and quartic polynomials). fit The difference between the values ​​can represent the snow depth; in this embodiment, the snow depth is calculated as follows:

[0083] SD=H fit2 -H fit1 ,

[0084] Where SD represents the snow depth value, and H fit2 H represents the snow surface elevation when there is snow cover. fit1 It indicates the elevation of the ground surface when there is no snow cover.

[0085] A further implementation method includes filtering the elevation difference between the light segments during the snow-covered and snowless periods at the intersection when calculating the snow depth value. If the elevation difference is negative or greater than the maximum snow depth in the current area, the corresponding light segment is removed. Specifically, if SD < 0, it indicates that the calculated snow depth is negative, which may be due to errors caused by a long time interval between the snowless period and the first or last day of snow cover, or by the growth of vegetation or the height of other objects on the ground exceeding the snow depth. In this case, the light segment is deleted.

[0086] If SD > SDmax, where SDmax is the maximum snow depth in the defined area, the light segment should also be deleted.

[0087] When the segment_landcover field of the light segment during the snow cover period and the light segment during the snowless period are both 0, 1 and 3, the elevation difference is not calculated and the snow depth value is directly set to 0; where segment_landcover field of 0 represents water body and segment_landcover field of 3 represents ice surface.

[0088] In low-latitude regions, there is often no snow during the snow cover period, meaning the segment_snowcover field is always 1. In this case, no elevation difference calculation is performed, and SD is directly set to 0.

[0089] S5: Snow depth estimation is achieved through the difference between crosspoints. Typically, when subtracting between two light segments, the crosspoints are not perfectly aligned. When forming snow depth crosspoints, there is a 0-17m gap between the light segments during the snow-covered and snow-free periods. Topographical undulations along this distance affect the snow depth estimation accuracy. The ICESat-2 calibration specifies a vertical elevation error of 10cm, but it also has a horizontal positioning error of approximately 3.5±2.1m. In areas with significant topographical undulations, the horizontal error will propagate vertically; therefore, the elevation error caused by topographical undulations must be taken into account.

[0090] A snow depth reference error model is constructed based on the terrain slope of the light segment during the snowless period and the horizontal distance between the light segment during the snow-covered period and the light segment during the snowless period when the intersection point is formed.

[0091] A further implementation method involves obtaining the horizontal distance between the light segment during the snow-covered period and the light segment during the snowless period when an intersection point is formed, including:

[0092] Convert the latitude and longitude of the intersection point into radians, and calculate the radian difference between the longitude and latitude of the light segment during the snowy period and the light segment during the snowless period. Specifically, the specific latitude and longitude values ​​between the light segments during the snowy and snowless periods are known. Let the latitude and longitude of the point with snow cover be (lat1, lon1), and the latitude and longitude of the point without snow cover be (lat2, lon2). The Earth's radius R is usually taken as 6,371,000 meters. The conversion formula is:

[0093] Radians=degrees×(π / 180),

[0094] Radians represent radians, and degrees represent degrees.

[0095] lat1 and lat2 need to be converted to radians:

[0096] lat'1=lat1×(π / 180),

[0097] lat'2=lat2×(π / 180),

[0098] Similarly, lon1 and lon2 need to be converted to radians:

[0099] lon'1 = lon1 × (π / 180),

[0100] lon'2 = lon2 × (π / 180).

[0101] Calculate the difference in radians between the longitude and latitude of two points:

[0102] Calculate the difference between the longitude and latitude of two points:

[0103] Δlat=lat'2-lat'1,

[0104] Δlon=lon'2-lon'1.

[0105] Based on the radian difference, the initial distances between the light segments during the snow-covered period and the light segments during the snowless period are obtained using the Haversine formula; the specific formula is as follows:

[0106]

[0107] The great circle angle c is calculated based on the initial distance between the light segment during the snow-covered period and the light segment during the snowless period; the specific calculation formula is as follows:

[0108]

[0109] Here, atan2 is a function that returns radians and is used to calculate the great circle angle between two points.

[0110] Based on the great circle angle and the Earth's radius, the horizontal distance between the light segment during the snowy season and the light segment during the snowless season is obtained. Furthermore, the actual distance between the two points is obtained by multiplying the great circle angle c by the Earth's radius (R):

[0111] d = R·c.

[0112] d is the horizontal distance between the light segments during the snowy and snowless periods, and its unit is the same as R, which is meters.

[0113] Based on the horizontal distance, the terrain slope of the snowless, clear section is known and is represented by the slope field. The vertical change over this distance can be calculated using the tan function, i.e., the snow depth reference error is calculated as follows:

[0114] Reference error=distance×tan(slope),

[0115] Among them, Reference error is the snow depth reference error; distance is the horizontal distance between the two light segments when the intersection point is formed; slope is the terrain slope of the light segment during the snowless period. A slope of "+" indicates an uphill slope along the track direction, and a slope of "-" indicates a downhill slope along the track direction.

[0116] S6: Correct the snow depth value based on the snow depth reference error model to obtain the corrected snow depth value, and complete the hemispherical scale snow depth estimation based on the spaceborne lidar.

[0117] In this embodiment, principles and indicators for quality inspection of the ICESat-2ATL08 snow depth dataset are also provided.

[0118] During the data quality inspection process, it must be ensured that the snow depth at the site and the ICESat-2 snow depth are on the same day.

[0119] The elevation difference between the station and ICESat-2 should not exceed 50m, and they should be kept on the same horizontal plane as much as possible.

[0120] The horizontal distance between the station and ICESat-2 is within 7km;

[0121] Specifically using Correlation, R 2 Evaluation metrics such as RMSE, MAE, and MedAE are defined as follows:

[0122]

[0123] In the formula, r xy It is the Pearson correlation coefficient between variables x and y, ranging from -1 to 1; x i and y i These are the estimated snow depth values ​​for the i-th snow depth in the sample; and These are the sample means of x and y, respectively, which are the average snow depths at each station; is the estimated value, and n is the sample size.

[0124] In the snow depth distribution density map of the dataset, the ICESat-2 ATL08 snow depth dataset contains a total of 1,133,709 data points, including 239,058 in January, 187,325 in February, 126,866 in March, 279,316 in April, and 301,144 in December. Snow depth points are densely distributed between 55°N and 70°N, while the snow depth is mostly zero between 30°N and 40°N.

[0125] Figure 5 A statistical chart classifies snow depth estimation deviations into five levels: >20cm, 5-20cm, -5-5cm, -20-(-5)cm, and <(-20)cm. Figure 5 It can be seen that the percentage of errors greater than 20cm is 19.3%, 5-20cm is 11.4%, -5-5cm is 60.5%, -20-(-5)cm is 6.1%, and <(-20)cm is 2.7%. Overall, the percentage of snow depth estimation errors >5cm is 30.7%, and the percentage of errors less than -5cm is 8.8%, indicating that the number of overestimations is much greater than the number of underestimations.

[0126] Figure 6 A statistical chart classifies the snow depth estimation deviation for different snow cover months. The validation data is divided into three parts with a snow depth estimation deviation of 5 cm as the threshold. Figure 6 The percentage of each part is indicated by numbers. Throughout the snow cover period, the portion with the largest deviation in snow depth estimation being between -5 and 5 cm accounted for the largest proportion, with the highest percentage in December at 70.2%, while the percentage was lower in March and April, reaching 50%. From Figure 6 It can be observed that in December and January, the proportion of snow depth estimation deviations greater than 5cm was 28.6% and 27.4%, respectively, while the proportion of deviations less than -5cm was almost zero. However, as the snow cover period progressed, the proportion of snow depth estimation deviations less than -5cm gradually increased, reaching 25.5% in April, exceeding the proportion of deviations greater than 5cm.

[0127] Figures 7-8 To illustrate the relationship between snow depth and estimation error across different ranges, the accuracy of ICESat-2 snow depth estimation is related to the actual snow thickness; different ranges of snow depth have different estimation errors. The relationship between the two is as follows: Figure 7The gray bars represent the buffer zone along the y=0 line, where the snow depth estimation error is within the range of -5 to 5 cm. When the snow depth is less than 5 cm, most snow depths are overestimated; when the snow depth is between 5 and 40 cm, both overestimation and underestimation occur; and when the snow depth is greater than 40 cm, most snow depths are underestimated. To quantitatively describe the relationship between the two, based on the approximate snow depth range in the ICESat-2 ATL08 dataset, 0-5 cm is defined as thin snow, 5-20 cm as moderate snow, 20-40 cm as deep snow, and greater than 40 cm as heavy snow. The number of snow depth estimation errors in different ranges was statistically analyzed, such as... Figure 8 It can be observed that the proportions of snow depth estimation errors greater than 20cm, 5-20cm, -5-5cm, -20-(-5)cm, and <(-20cm) are as follows: for thin snow, 0.2%, 14.0%, 85.8%, 0, and 0, respectively; for moderate snow, 2.5%, 20.3%, 66.9%, 10.3%, and 0, respectively; for deep snow, 1.2%, 9.6%, 51.4%, 36.9%, and 0.9%, respectively; and for heavy snow, 1.3%, 2.7%, 50.4%, 41.4%, and 4.2%, respectively.

[0128] Figures 9-10 To illustrate the relationship between terrain slope and snow depth estimation errors, the dataset created in this invention centralizes the relationship between slope and snow depth estimation errors as follows: Figure 9 It can be observed that as the slope increases from 0° to 5°, the snow depth estimation error does not change significantly, and the difference in R between the two remains constant. 2 The value is 0. This indicates that the dynamic discrimination rule for the intersection distance threshold limited by the slope in the above embodiment has almost no impact on the snow depth estimation when the central slope of the dataset is less than 5°. To further intuitively verify the relationship between slope and snow depth estimation error, the experiment retained the light segment with a terrain slope greater than 5°, and selected the light segment with consistent underlying surface, snow cover date, snowless date, intersection distance, distance between verification point and station, altitude, and other conditions. Only the influence of slope on snow depth measurement error was considered. Figure 10 The study found a positive correlation between the absolute errors in slope and snow depth estimation, with a correlation coefficient of 0.85 and an R-value of [missing value]. 2 The error is 0.72, MAE is 8.39cm, and RMSE is 12.58cm. For every 10-degree increase in slope, the error increases by approximately 50cm; when the slope is less than 5°, the error is less than 25cm, and when the slope is greater than 20°, the error exceeds 1m.

[0129] Figure 11The graph shows the relationship between altitude and snow depth estimation error. The experiment analyzed the ICESat-2 snow depth estimation results at different altitudes. A total of 10268 ICESat-2 snow depth points were used for validation. Figure 11 It can be observed that at the same altitude, there will be different snow depth estimation errors, and the error does not change significantly with increasing altitude. The correlation between altitude and snow depth estimation error is -0.01, R0. 2 The correlation between the two values ​​is 0, indicating a very poor correlation. Altitude does not affect the estimation of snow depth.

[0130] Figure 12 The graph shows the relationship between the distance between intersection points and the relative error in snow depth estimation. To achieve a balance between snow depth estimation error and the number of intersection points, the experiment adopted a dynamic discrimination rule for the intersection point distance threshold limited by slope: when the slope is ≥2°, the distance is set to 2m; when the slope is <2°, the distance is set to 17m. The diameter of the ICESat-2 spot foot point is 17m; therefore, the center distance between intersection points in the dataset should be within the range of 0-17m. The relationship between the distance between snow depth intersection points and the snow depth estimation error is shown below. Figure 12 It can be observed that the error changes irregularly as the distance increases, indicating that the dynamic discrimination rule of the intersection distance threshold limited by the slope has a good effect.

[0131] Figure 13 The relationship between the distance between the station and the intersection point and the relative error of snow depth estimation was discussed. Different colors represent the difference between the station's elevation and the elevation of the snow depth intersection point, with each level divided into 10m increments, for a total of 5 levels. The smaller the elevation difference, the more likely the snow depth intersection point and the station are at the same level, and the more representative the station is of the snow depth value at the surrounding intersection points. Therefore, the smaller the elevation difference, the more reliable the result at that point. Figure 13 It can be seen that as the distance increases, the relative error of the snow depth deviation on the fitted line increases very little, indicating that in most cases, the snow depth at a station can represent the snow depth value of ICESat-2 ATL08 within that range. In the dataset created in this invention, we have removed non-flat areas with a slope greater than 5°. In flat areas, when the difference between the station elevation and the light segment elevation is small (less than 50m), the snow depth at the station within a 7km range can verify the snow depth at the ICESat-2 snow depth intersection point.

[0132] Figure 14 The relationship between snowless months and snow depth estimation error was discussed. The ICESat-2 snow depth value is the difference between the light range during the snow-covered period and the snowless period. The choice of snowless dates directly affects the snow depth estimation. Figure 14It can be observed that snow depth is generally overestimated. The closer the snowless period is to the snow cover period, the fewer outliers appear in the dataset, such as in November 2019 and April 2020. Specifically, the medians for October 2019, November 2019, April 2020, and May 2020 are all greater than 0, with a mean error greater than 0, and the error distribution is right-skewed, indicating that snow depth is overestimated. Conversely, the medians for June 2020 and August 2020 are less than 0, with a mean error less than 0, and the error distribution is left-skewed, indicating that snow depth is underestimated. In August, the Northern Hemisphere is in summer, with vigorous vegetation growth, leading to an overestimation of surface elevation and indirectly causing the most severe underestimation of snow depth. Between 50°N and 70°N, the snow cover months selected in the dataset are from December 2019 to April 2020. November 2019 and May 2020 are closest to the snow cover period, and the boxes for these two months are the flattest, with the most concentrated snow depth errors. Therefore, the smaller the time difference between the non-snow cover period and the snow cover period, the smaller the change in surface elevation, and the more accurate the snow depth estimation.

[0133] Figure 15 To comprehensively analyze the factors that cause errors in snow depth estimation, the degree of influence on snow depth estimation, from greatest to least, is: slope, underlying surface type, different months with snow cover, snow depth range, snowless months, distance between the measured station and the snow depth intersection point, distance between snow depth intersection points, and altitude.

[0134] Seasonal snowmelt plays a vital role in the global water cycle and water balance. ICESat-2 can be used to generate large-scale, high-precision snow depth datasets, which are not limited by factors such as altitude, climate, or cost. This data can compensate for sparse snow depth data from individual stations and provide a reference for snow cover model construction and snow depth measurement in high-altitude mountainous areas. This embodiment utilizes the ICESat-2 ATL08 product, and by analyzing the relationship between topography and light segment geometry, the snow depth value is obtained by subtracting the light segment elevation between snowy and snowless periods, thus acquiring a seasonal snow depth dataset for the Northern Hemisphere from 2018 to 2020.

[0135] The ICESat-2 ATL08 snow depth dataset was validated using snow depth data from over 43,000 stations in the Northern Hemisphere. The results showed that slope was the most significant factor causing snow depth estimation errors. Estimation results were relatively reliable when the slope was less than 5°; however, the greater the slope, the higher the error, with an increase of approximately 50 cm for every 10° increase in slope. Significant differences in snow depth estimation performance were observed across different underlying surfaces, including bare land, herbaceous vegetation, herbaceous wetlands, cultivated land, shrubs, open forests, and closed forests. 2The values ​​were 0.88, 0.79, 0.70, 0.68, 0.51, 0.21, and 0.18, respectively. Surfaces with sparse vegetation and minimal elevation changes during freeze-thaw cycles are more suitable for ICESat-2 snow depth estimation. Different snow cover months affect snow depth estimation. During the accumulation period, the snow is shallower, and slight elevation changes caused by topography and surface differences often exceed the snow depth, resulting in overestimation. During the ablation period, snow particles are larger and have wider gaps, allowing laser signals to penetrate part of the snow layer, leading to underestimation of snow depth. The ICESat-2 snow depth dataset created in this invention exhibits strong spatial heterogeneity; overestimation or underestimation of snow depth may occur simultaneously at different points within a small region.

[0136] Example 2

[0137] This invention also provides a hemispherical-scale snow depth estimation system based on spaceborne lidar, and a method for implementing this system, including:

[0138] The data acquisition module is used to download snow-related data within a preset hemispherical scale collected by the spaceborne lidar and to preprocess the snow-related data.

[0139] The light segment segmentation module is used to obtain the light segments during the snow-covered period and the light segments during the snowless period within a preset hemispherical scale based on the preprocessed snow-related data.

[0140] The intersection point acquisition module is used to obtain the intersection points of the light segments during the snowy season and the light segments during the snowless season by adopting a dynamic discrimination rule based on the intersection point distance threshold that is limited by the slope. The dynamic discrimination rule based on the intersection point distance threshold that is limited by the slope includes: setting a threshold for judging the center distance between intersection points based on a preset terrain slope threshold within a preset hemispherical scale.

[0141] The snow depth calculation module is used to obtain the snow depth value within a preset hemispherical scale based on the intersection point and the snow depth estimation model; the snow depth estimation model is to subtract the elevation values ​​of the light segment during the snow cover period and the snowless period at the intersection point.

[0142] The error model construction module is used to construct a snow depth reference error model based on the terrain slope of the light segment during the snowless period and the horizontal distance between the light segment during the snow-covered period and the light segment during the snowless period when they form an intersection.

[0143] The error correction module is used to correct the snow depth value based on the snow depth reference error model, obtain the corrected snow depth value, and complete the hemispherical scale snow depth estimation based on the spaceborne lidar.

[0144] The embodiments described above are merely preferred embodiments of the present invention and are not intended to limit the scope of the present invention. Various modifications and improvements made to the technical solutions of the present invention by those skilled in the art without departing from the spirit of the present invention should fall within the protection scope defined by the claims of the present invention.

Claims

1. A method for estimating snow depth at a hemispherical scale based on spaceborne lidar, characterized in that, The method includes: Download snow-related data within a preset hemispherical scale collected by the spaceborne lidar, and preprocess the snow-related data; Based on the preprocessed snow cover-related data, the light segments during the snow cover period and the light segments during the snowless period are obtained within a preset hemispherical scale; A dynamic discrimination rule based on a slope-limited intersection distance threshold is adopted to obtain the intersection points of the light segment during the snowy period and the light segment during the snowless period; wherein, the dynamic discrimination rule based on a slope-limited intersection distance threshold includes: within a preset hemispherical scale, setting a threshold for judging the center distance between intersection points based on a preset terrain slope threshold; Based on the intersection point and the snow depth estimation model, the snow depth value within a preset hemispherical scale is obtained; wherein, the snow depth estimation model is the difference between the elevation values ​​of the light segment during the snow-covered period and the light segment during the snowless period at the intersection point; Based on the terrain slope of the light segment during the snowless period and the horizontal distance between the light segment during the snow-covered period and the light segment during the snowless period when the intersection is formed, a snow depth reference error model is constructed. The snow depth value is corrected based on the snow depth reference error model to obtain the corrected snow depth value, thus completing the hemispherical scale snow depth estimation based on the spaceborne lidar. The method for obtaining the horizontal distance between the light segment during the snow-covered period and the light segment during the snowless period when the intersection point is formed includes: Convert the latitude and longitude corresponding to the intersection point into radians, and calculate the radian difference between the longitude and latitude of the light segment during the snowy period and the light segment during the snowless period at the intersection point; Based on the aforementioned radian difference, the initial distance between the light segment during the snow cover period and the light segment during the snowless period is obtained using the Haversine formula; Calculate the great circle angle based on the initial distance between the light segment during the snow-covered period and the light segment during the snowless period; Based on the great circle angle and the Earth's radius, the horizontal distance between the light segment during the snow-covered period and the light segment during the snowless period is obtained; Based on the horizontal distance, the terrain slope of the bare section during the snowless period is known and is represented by the slope field. The vertical change over this distance is calculated using the tan function, i.e., the snow depth reference error is calculated as follows: Reference error=distance×tan(slope), Among them, Reference error is the snow depth reference error; distance is the horizontal distance between the two light segments when the intersection point is formed; slope is the terrain slope of the light segment during the snowless period. A slope of "+" indicates an uphill slope along the track direction, and a slope of "-" indicates a downhill slope along the track direction.

2. The method according to claim 1, characterized in that, The snow-related data includes location data, elevation data, slope data, land use type data, and snow cover data within a preset hemispherical scale. The method for preprocessing the snow cover-related data includes: Remove light segments with a slope greater than 5°, as well as light segments related to water bodies, buildings, and ice surfaces from the snow accumulation-related data.

3. The method according to claim 1, characterized in that, The snow cover period includes: for the period between 30°N and 55°N, the snow cover months are defined as December, January, February, and March, totaling 4 months; and for the period between 55°N and 70°N, the snow cover months are defined as December, January, February, March, and April, totaling 5 months.

4. The method according to claim 1, characterized in that, The method for obtaining the snowless period light segment includes: Import the MOD10A1 snow cover product into the GEE software, set the band threshold of the Normalized Difference Snow Index (NDSI) of the MOD10A1 snow cover product to 0.

2. When the NDSI of the MOD10A1 snow cover product is less than 0.2, it is determined that there is no snow cover; otherwise, it is determined that there is snow cover. Import the snow cover map into ArcGIS software and divide it into a 5°×5° grid to obtain the snow cover situation of each 5° cell within the preset date. Change the date to obtain the first and last days of snow cover for each cell. Based on the first day of snow accumulation and the last day of snow accumulation, the snowless period is obtained; The snowless period is compared with the corresponding time period's segment_landcover field in the preprocessed snow cover data. If the segment_landcover field in the preprocessed snow cover data is 1, the snowless period is correctly identified, and the snowless period light segment is obtained. If the segment_landcover field is 2, the snowless period is incorrectly identified, and the segment_landcover field comparison is repeated on a different snowless date until the segment_landcover field is 1, thus obtaining the snowless period light segment. Here, segment_landcover is a snow cover identifier; a segment_landcover field of 1 represents snowless land, and a segment_landcover field of 2 represents land with snow.

5. The method according to claim 1, characterized in that, The dynamic discrimination rule based on the intersection distance threshold limited by slope specifically includes: The preset terrain slope threshold is 2°. When the terrain slope is greater than or equal to 2°, the threshold for determining the center distance between intersections is set to 2m; when the terrain slope is less than 2°, the threshold for determining the center distance between intersections is set to 17m.

6. The method according to claim 1, characterized in that, This also includes filtering the elevation difference between the light segments during the snow cover period and the snowless period at the intersection when calculating the snow depth value. When the elevation difference is negative or greater than the maximum snow depth in the current area, the corresponding light segment is removed. When the segment_landcover field of the light segment during the snow cover period and the light segment during the snowless period are both 0, 1 and 3, the elevation difference is not calculated and the snow depth value is directly set to 0. Among them, a segment_landcover field of 0 represents water body and a segment_landcover field of 3 represents ice surface.

7. A hemispherical-scale snow depth estimation system based on spaceborne lidar, used to implement the method described in any one of claims 1-6, characterized in that, include: The data acquisition module is used to download snow-related data within a preset hemispherical scale collected by the spaceborne lidar and to preprocess the snow-related data. The light segment segmentation module is used to obtain the light segments during the snow cover period and the light segments during the snowless period within a preset hemispherical scale based on the preprocessed snow cover-related data. The intersection point acquisition module is used to obtain the intersection points of the light segment during the snowy season and the light segment during the snowless season by adopting a dynamic discrimination rule based on the intersection point distance threshold limited by the slope; wherein, the dynamic discrimination rule based on the intersection point distance threshold limited by the slope includes: setting a threshold for judging the center distance between intersection points based on a preset terrain slope threshold within a preset hemispherical scale; The snow depth calculation module is used to obtain the snow depth value within a preset hemispherical scale based on the intersection point and the snow depth estimation model; wherein, the snow depth estimation model is the difference between the elevation values ​​of the light segment during the snow-covered period and the light segment during the snowless period at the intersection point; The error model construction module is used to construct a snow depth reference error model based on the terrain slope of the light segment during the snowless period and the horizontal distance between the light segment during the snow-covered period and the light segment during the snowless period when the intersection is formed. The error correction module is used to correct the snow depth value based on the snow depth reference error model, obtain the corrected snow depth value, and complete the hemispherical scale snow depth estimation based on the spaceborne lidar.

Images