Hemisphere scale snow depth estimation method and system based on satellite-borne laser radar

Abstract

The invention relates to the technical field of laser radar snow depth estimation, and particularly provides a hemisphere scale snow depth estimation method and system based on a satellite-borne laser radar, and the method comprises the steps: downloading and preprocessing accumulated snow related data collected by the satellite-borne laser radar in a preset hemisphere scale; based on the pre-processed accumulated snow related data, a snow accumulation period light section and a snow-free period light section in a preset hemispherical scale are obtained; adopting a gradient limitation-based intersection distance threshold dynamic discrimination rule to obtain an intersection of the snow accumulation period light section and the snow-free period light section; obtaining a snow depth value based on the intersection point and a snow depth estimation model; and correcting the snow depth value based on the snow depth reference error model to obtain the corrected snow depth value, and completing the hemispherical scale snow depth estimation based on the satellite-borne laser radar. The method is not limited by factors such as altitude and cost when making the snow depth data, can compensate for sparse site snow depth data, and plays an important role in authenticity inspection of snow depth products, construction of a snow accumulation model and measurement of the snow depth of a high-altitude mountainous area.

Description

Technical Field

[0001] The present invention belongs to the technical field of snow depth estimation using laser radar, and in particular relates to a hemispherical scale snow depth estimation method and system based on satellite-borne laser radar. Background Art

[0002] Snow is the most widely distributed element in the cryosphere. It has high albedo and low thermal conductivity, and is of great significance to the flow of surface energy and the regulation of regional and global ecosystems. At the same time, snowmelt also plays an important role in the global water cycle and water balance. Snow depth is one of the most important and basic parameters in snow research, and is also an important factor in predicting, monitoring, and warning of snowmelt flood disasters.

[0003] At present, global or hemispheric scale, long-term series of snow depth data mainly come from passive microwave inversion or atmospheric reanalysis data. However, the spatial resolution of the data is low; the existing mainstream passive microwave inversion algorithm tends to saturate when the snow depth reaches about 50cm; in high-altitude mountainous areas, there is great uncertainty in the inversion of snow depth, and there is a lack of suitable observational data to verify it. Direct manual measurement or observation using meteorological stations is the most accurate and intuitive way to obtain snow depth data, and this part of the data is also widely used in the verification of snow depth products. However, due to the limitations of terrain, altitude and natural conditions, the observation efficiency is low, and the data is discontinuous and uneven. Especially in high-altitude remote mountainous areas, direct snow depth measurements are very rare, and the spatial heterogeneity of snow accumulation is particularly significant due to factors such as terrain and climate. Therefore, it is crucial to develop a dataset with wide coverage, high spatial resolution, inversion process not limited by snow depth threshold, and the ability to verify snow depth in high-altitude mountainous areas.

[0004] By comparing the elevation difference at the intersection of ICESat tracks before and after snow accumulation, the depth of new snow on the ice sheet can be estimated, and the results have been proven to be reliable. ICESat and high-precision DEM can be used to obtain the snow depth on small-scale mountain glaciers, but ICESat has only one laser beam with sparse footprints, which makes it difficult to apply to snow depth measurement in non-polar regions. ICESat-2 has a total of 6 laser beams, with a spot footprint size of 17m and a spacing of 0.7m between adjacent spots, which increases the density of cross-track sampling. The surface elevation accuracy can reach 0.1m, which can achieve large-scale and large-scale intensive observation of snow accumulation. ICESat-2's standard products ATL06 and ATL08 have the potential to estimate snow depth in mid- and high-latitude regions. Most of them measure snow depth through the elevation of ICESat-2 during the snow accumulation period and the accurate digital terrain model (DTM) during the snow-free period. Related studies have shown that the accuracy of ICESat-2 snow depth estimation depends on terrain slope and forest canopy density, but these methods have high-precision ground reference elevations for alignment in the verified areas, and cannot achieve snow depth estimation over a large range. Using only the elevation difference between the intersection points of the ICESat-2 ATL08 product track can accurately estimate the snow depth in flat areas, but there will be greater errors in mountainous areas. However, only a small number of intersection points were used in the verification, and the randomness is too high. Using the SlideRule Earth project, the ATL06 product was improved and the snow depth at different vegetation coverage, altitudes, and slopes was verified. The results show that ICESat-2 SlideRule may provide a new snow depth dataset across the western United States and even the global land.

[0005] ICESat-2 has good snow depth estimation capabilities, but the selected area is limited, the snow cover type is single, and the method has poor universality, and is only suitable for snow depth estimation under specific background conditions. Summary of the invention

[0006] In order to solve the problems existing in the prior art, the present invention proposes a hemispherical-scale snow depth estimation method and system based on satellite-borne laser radar, which obtains the snow depth value by subtracting the intersection points of the light segments during the snowy and snowless periods; the method takes into account the relationship between the terrain and the geometric structure of the light segments, uses the terrain slope to dynamically determine the distance between the intersection points, and corrects the snow depth estimation error caused by the terrain undulation. By testing the results obtained by the present invention, it can be found that the method can provide a new technical approach for obtaining snow depth data.

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

[0008] Downloading snow accumulation related data within a preset hemispherical scale collected by a satellite-borne laser radar, and preprocessing the snow accumulation related data;

[0009] Based on the pre-processed snow-related data, the light segments of the snow-covered period and the snow-free period within the preset hemispherical scale are obtained;

[0010] The intersection point of the light segment in the snow accumulation period and the light segment in the snow-free period is obtained by adopting a dynamic discrimination rule of the intersection point distance threshold defined by the slope; wherein the dynamic discrimination rule of the intersection point distance threshold defined by the slope comprises: setting a threshold for discriminating the center distance between intersection points based on a preset terrain slope threshold within a preset hemispherical scale;

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

[0012] A snow depth reference error model is constructed based on the terrain slope of the snow-free period light segment and the horizontal distance between the snow-covered period light segment and the snow-free period light segment when the intersection is formed;

[0013] The snow depth value is corrected based on the snow depth reference error model to obtain a corrected snow depth value, thereby completing hemispherical-scale snow depth estimation based on satellite-borne laser radar.

[0014] Preferably, the snow accumulation 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 accumulation related data includes:

[0016] The light segments with a slope greater than 5° and the light segments of water bodies, buildings and ice surfaces in the snow accumulation related data are eliminated.

[0017] Preferably, the snow accumulation period light segment includes: the snow accumulation months in the light segment determined at 30°N-55°N are December, January, February, and March, a total of 4 months; the snow accumulation months in the light segment determined at 55°N-70°N are December, January, February, March, and April, a total of 5 months.

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

[0019] Import the MOD10A1 snow product into the GEE software, set the NDSI of the MOD10A1 snow product to 0.2, and when the NDSI of the MOD10A1 snow product is less than 0.2, it is judged as no snow cover, otherwise it is judged as snow cover;

[0020] The snow cover map was imported into ArcGIS software and divided into 5°×5° grids to obtain the snow cover situation on each 5° pixel within the preset date. The date was changed to obtain the first and last days of snow cover for each pixel.

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

[0022] The snow-free period is compared with the segment_landcover field of the corresponding time period in the pre-processed snow accumulation related data. If the segment_landcover field of the corresponding time period in the pre-processed snow accumulation related data is 1, the snow-free period time period is correctly judged, and the snow-free period light segment is obtained; if the segment_landcover field is 2, the snow-free period time period is incorrectly judged, then the snow-free date is changed and the segment_landcover field comparison is repeated until the segment_landcover field is 1, and the snow-free period light segment is obtained; wherein, segment_landcove is a snow cover identifier, a segment_landcover field of 1 represents snow-free land, and a segment_landcover field of 2 represents snowy land.

[0023] Preferably, the dynamic discrimination rule of intersection distance threshold based on slope limitation 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 in the snow accumulation period and the snow-free period at the intersection is screened, and when the elevation difference is negative or greater than the maximum snow depth in the current area, the corresponding light segment is eliminated; when the segment_landcover fields of the light segments in the snow accumulation period and the snow-free period are both 0, 1, and 3, the elevation difference calculation is not performed, and the snow depth value is directly set to 0; wherein, a segment_landcover field of 0 represents a water body, and a segment_landcover field of 3 represents an ice surface.

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

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

[0028] Based on the arc difference, the initial distance between the light segment during the snow accumulation period and the light segment during the snow-free period is obtained by using the Haversine formula;

[0029] The great circle angle is calculated based on the initial distance between the light segment during the snowy period and the light segment during the snowless period;

[0030] Based on the great circle angle and the radius of the earth, the horizontal distance between the light segment during the snow accumulation period and the light segment during the snow-free period is obtained.

[0031] The present invention also provides a hemispherical scale snow depth estimation system based on a satellite-borne laser radar, which is used to implement the method, comprising:

[0032] A data acquisition module, used to download snow accumulation related data within a preset hemispherical scale collected by a satellite-borne laser radar, and pre-process the snow accumulation related data;

[0033] A light segment division module is used to obtain the light segments of the snow accumulation period and the light segments of the snow-free period within a preset hemispherical scale based on the preprocessed snow accumulation related data;

[0034] The intersection acquisition module is used to obtain the intersection of the light segment during the snow accumulation period and the light segment during the snow-free period by adopting a dynamic discrimination rule of the intersection distance threshold defined by the slope; wherein the dynamic discrimination rule of the intersection distance threshold defined by the slope includes: setting a threshold for discriminating the center distance between intersections based on a preset terrain slope threshold within a preset hemispherical scale;

[0035] A snow depth value 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 to make a difference between the elevation values ​​of the light segment during the snow accumulation period and the light segment during the snow-free period at the intersection point;

[0036] An error model building module is used to build a snow depth reference error model based on the terrain slope of the snow-free period light segment and the horizontal distance between the snow-covered period light segment and the snow-free period light segment 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 to obtain the corrected snow depth value and complete the hemispherical scale snow depth estimation based on the satellite-borne laser radar.

[0038] Compared with the prior art, the beneficial effects of the present invention are as follows: the present invention utilizes the ICESat-2 data of the satellite-borne lidar, through the relationship between the terrain and the geometric structure of the light segment, combined with the dynamic discrimination of the intersection distance threshold limited by the slope, and the light segment elevations during the snowy period and the snow-free period are subtracted to obtain the snow depth value, and uses means such as elevation difference screening and terrain geometric structure correction to reduce the snow depth estimation error, and finally obtains the seasonal snow depth data set in the Northern Hemisphere from 2018 to 2020.

[0039] The present invention provides a new idea and method for obtaining snow depth data. This method is not restricted by factors such as altitude and cost when producing snow depth data, and can make up for sparse snow depth data at stations. It plays an important role in the authenticity verification of snow depth products, the construction of snow accumulation models, and the measurement of snow depth in high-altitude mountainous areas. BRIEF DESCRIPTION OF THE DRAWINGS

[0040] In order to more clearly illustrate the technical solution of the present invention, the following briefly introduces the drawings required for use in the embodiments. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying creative labor.

[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 It is a schematic diagram of a processing flow of an ICESat-2 snow depth estimation method according to an embodiment of the present invention;

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

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

[0045] Figure 5 A classified statistical diagram of the number of snow depth estimation deviations according to an embodiment of the present invention;

[0046] Figure 6 A classified statistical diagram of snow depth estimation deviations in different snow accumulation months according to an embodiment of the present invention;

[0047] Figure 7 A scatter plot of the snow depth of the site and the estimated deviation of the snow depth of ICESat-2 according to an embodiment of the present invention;

[0048] Figure 8 A statistical diagram of snow depth estimation error levels in different ranges according to an embodiment of the present invention;

[0049] Fig. 9 This is a relationship diagram between terrain slope and snow depth estimation error in an embodiment of the present invention (slope≤5° in the data set);

[0050] Fig.10 A relationship diagram between terrain slope and snow depth estimation error according to an embodiment of the present invention;

[0051] Fig.11 A relationship diagram of altitude and snow depth estimation error according to an embodiment of the present invention;

[0052] Fig.12 A relationship diagram of the distance between intersection points and the snow depth estimation deviation according to an embodiment of the present invention;

[0053] Fig.13A graph showing the relationship between the distance between the station and the intersection and the snow depth estimation deviation according to an embodiment of the present invention;

[0054] Fig.14 A relationship diagram of snow-free months and snow depth estimation deviation according to an embodiment of the present invention;

[0055] Fig.15 This is a heat map of factors that affect the relative error of ICESat-2 snow depth estimation according to an embodiment of the present invention. DETAILED DESCRIPTION

[0056] The following will be combined with the drawings in the embodiments of the present invention to clearly and completely describe the technical solutions in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0057] In order to make the above-mentioned objects, features and advantages of the present invention more obvious and easy to understand, the present invention is 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 the present invention are described as follows:

[0059] Google Earth Engine (GEE) is a cloud computing-based geospatial analysis platform designed to provide geospatial data storage, processing, and analysis services on a global scale. GEE has a PB-level catalog of open and free geospatial datasets, and users can instantly access years of historical data and the latest real-time data. Based on cloud infrastructure, GEE can automatically subdivide and distribute computing tasks, achieve large-scale parallel processing, and support high-throughput analysis. The platform integrates powerful visualization tools and supports real-time viewing of analysis results, which is convenient for users to understand intuitively. GEE provides JavaScript and Python APIs, and users can quickly implement complex geospatial analysis through 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. It belongs to the snow cover dataset and provides daily global 500-meter resolution snow cover information from the Terra satellite. The data content includes: Snow cover, which expresses the proportion of snow in each pixel as a percentage. Snow albedo: provides albedo information on the snow surface. Quality assessment (QA) data: contains quality assessment information of snow cover data to help users judge the reliability of the data. NDSI index: Normalized difference snow index, used to identify snow-covered areas. MOD10A1 products are widely used in meteorology, climate research, water resources management, and environmental monitoring.

[0061] ArcGIS is a powerful geographic information system (GIS) software that is widely used in many fields such as map production, spatial analysis, data management, and geographic information visualization.

[0062] Embodiment 1

[0063] like Figure 1 , Figure 2 As shown in FIG. 1 , a hemispherical-scale snow depth estimation method based on space-borne lidar includes:

[0064] S1: Batch download snow data collected by satellite-borne lidar within a preset hemispherical scale, and preprocess the snow data; batch download data, read lat, lon, slope, h_te_best_fit (H fit ), landcover, snowcover, date and other information, where lat represents latitude, lon represents longitude, slope represents slope, date represents date, and index H fit It represents the best-fit elevation of the light segment on the 100-meter segment, landcover represents the land use type, snowcover represents whether there is snow, and date represents the snow accumulation date.

[0065] When downloading and reading snow-related data in batches, 30°N-70°N is divided into 5°×5° blocks one by one; a further implementation method is that the snow-related data includes position data, elevation data, slope data, land use type data and snow cover data within a preset hemispherical scale; in this embodiment, it also includes using snow-related data to construct a snow depth proxy indicator, and the snow depth proxy indicator includes snow depth data and snow depth reference error data.

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

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

[0068] S2: Based on the preprocessed snow-related data, the light segments of the snow-covered period and the snow-free period are obtained within a preset hemispherical scale.

[0069] A further implementation method is that the snow accumulation period light segment includes: the snow accumulation months in the light segment determined at 30°N-55°N are December, January, February, and March, a total of 4 months; the snow accumulation months in the light segment determined at 55°N-70°N are December, January, February, March, and April, a total of 5 months.

[0070] A further embodiment is that the method for obtaining the snow-free period includes:

[0071] In order to make the snow-free period closer to the first day of snow accumulation or the last day of snow accumulation, the present invention uses the MOD10A1 product on the GEE platform to determine the snow-free period; first, the NDSI band is used for binarization processing, and the threshold is set to 0.2, so as to obtain data on whether there is snow or no snow in the determined area on the determined date. Specifically, the MOD10A1 snow product is imported into the GEE software, and the band threshold of the MOD10A1 snow product NDSI is set to 0.2. When the band threshold of the MOD10A1 snow product NDSI is less than 0.2, it is determined to be without snow cover, otherwise it is determined to be covered with snow, and a binary snow cover map can be obtained;

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

[0073] Based on the first day of snow accumulation and the last day of snow accumulation, the snow-free period is obtained.

[0074] The snow-free period is compared with the segment_landcover field of the corresponding time period in the preprocessed snow-related data. If the segment_landcover field of the corresponding time period in the preprocessed snow-related data is 1, the snow-free period is correctly judged and the snow-free period light segment is obtained; if the segment_landcover field is 2, the snow-free period is incorrectly judged, then the snow-free date is changed and the segment_landcover field comparison is repeated until the segment_landcover field is 1 to obtain the snow-free period light segment; wherein, segment_landcove is the snow cover identifier, the segment_landcover field is 1 for snow-free land, and the segment_landcover field is 2 for snow-covered land.

[0075] Repeat the above cross-comparison process until the segment_landcover field of the MOD10A1 image and the ICESat-2 shows no snow on the same day and in the same processing unit. This day can be determined as the earliest snow-free date that appears in the box. This method can shorten the time interval between the snow-free period and the snow accumulation period.

[0076] S3: Adopting the dynamic discrimination rule of intersection distance threshold based on slope limitation to obtain the intersection point of the light segment during the snow accumulation period and the light segment during the snow-free period; wherein, the dynamic discrimination rule of intersection distance threshold based on slope limitation includes: setting the threshold for discriminating the center distance between intersections based on the preset terrain slope threshold within the preset hemispherical scale. Figure 3 As shown, Figure 3 (a) The intersection point formed by repeated orbits with different periods: the satellite returns to this point again after a complete revisit cycle. The two tracks appear to be parallel but slightly offset. Figure 3 (b) Represents the intersection points formed by different orbits: the satellite does not necessarily go through a complete revisit cycle, and the intersection points may be formed at any time.

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

[0078] In this embodiment, there are two ways to form intersections of ICESat-2 footprints. The first is the intersection formed by repeated orbits of different periods, that is, the two orbits appear to be in a parallel state. This method will form more dense and continuous intersections. The second is the intersection formed by different orbits, which form fewer and sparse intersections. In low and medium latitudes, satellite tracks are sparse, and track deviations will occur during revisits, which are basically discrete intersections formed by different tracks.

[0079] The intersection points of footprints are not completely overlapped. Since the diameter of the ICESat-2 spot foot is 17m, the foot center distance within the range of 0-17m can be considered to be an intersection. When the distance is equal to 0, the two light segments are considered to be completely overlapped. Although the smaller the distance, the more it can be said that the two light segments represent the same geographical location, it will result in too few intersection points that meet the conditions; when the distance threshold is set to 17m, the error in areas with large terrain undulations is too large. Figure 4 shown.

[0080] The present invention dynamically selects the distance threshold based on the relationship between the snow depth error and the slope. When the slope is greater than 5 degrees, the snow depth error exceeds 30cm, so the present invention only considers areas with a slope less than 5 degrees; the preset terrain slope threshold is 2°. When the slope is greater than or equal to 2°, it is a gentle slope (2-5°), and the terrain will have large ups and downs. In order 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 of intersection distance threshold limited by slope can generate more intersections under the condition of smaller snow depth estimation error.

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

[0083] SD=H fit2 -H fit1 ,

[0084] Among them, SD represents the snow depth, H fit2 Indicates the snow surface elevation when covered with snow, H fit1 Represents the elevation of the ground surface when there is no snow cover.

[0085] A further implementation method is to screen the elevation difference between the light segments in the snow accumulation period and the snow-free period at the intersection when calculating the snow depth value, and remove the corresponding light segment when the elevation difference is negative or greater than the maximum snow depth in the current area. Specifically, if SD<0, it means that the calculated snow depth is negative, which may be due to the error caused by the long time interval between the snow-free period and the first day of snow accumulation or the last day of snow accumulation, and the growth of vegetation or the height of other objects on the surface has exceeded the snow depth, so the light segment is deleted;

[0086] If SD>SDmax, SDmax is the maximum snow depth value in the determined area, and the light segment is also deleted.

[0087] When the segment_landcover fields of the snow-covered period and the snow-free period are both 0, 1, and 3, no elevation difference calculation is performed and the snow depth is directly set to 0; among them, the segment_landcover field is 0 for water bodies, and the segment_landcover field is 3 for ice surfaces.

[0088] In low-latitude areas, there is no snow during the snow accumulation period most of the time, that is, when the segment_snowcover field is 1, no height difference calculation is performed in this part, and SD is directly set to 0.

[0089] S5: Snow depth estimation is achieved by making a difference between intersection points. Usually, when making a difference between two light segments, the intersection points of footprints are not completely overlapped. When forming a snow depth intersection point, there will be a 0-17m interval between the light segments during the snow accumulation period and the snow-free period. The terrain undulation over this distance will affect the snow depth estimation effect. The vertical altimetry error calibrated by ICESat-2 is 10cm, but it also has a horizontal positioning error of about 3.5±2.1m. In areas with large terrain undulations, the horizontal error will be transmitted in the vertical direction. Therefore, the altimetry error caused by terrain undulations needs to be taken into account.

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

[0091] A further implementation method is that the method for obtaining the horizontal distance between the light segment during the snow accumulation period and the light segment during the snow-free period when the intersection is formed includes:

[0092] Convert the longitude and latitude corresponding to the intersection point into radians, and calculate the radian difference between the longitude and latitude of the snowy period light segment and the snow-free period light segment at the intersection point. Specifically, the specific longitude and latitude values ​​between the snowy period light segment and the snow-free period light segment are known, and the longitude and latitude of the snowy period point are recorded as: (lat 1 ,lon 1 ), the latitude and longitude of the snow-free period point are: (lat 2 ,lon2 ), the radius of the earth R is usually taken as 6371000 meters; the conversion formula is:

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

[0094] Radians represent radians and degrees represent degrees.

[0095] lat 1 and lat 2 Need to convert to radians:

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

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

[0098] Similarly, lon 1 and lon 2 Need to convert to radians:

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

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

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

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

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

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

[0105] Based on the arc difference, the initial distance between the light segment during the snow-covered period and the light segment during the snow-free period is obtained by 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 accumulation period and the light segment during the snow-free period. The specific calculation formula is as follows:

[0108]

[0109] 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 radius of the earth, the horizontal distance between the light segment during the snowy period and the light segment during the snowless period is obtained. Furthermore, the actual distance between the two points is obtained by multiplying the great circle angle c by the radius of the earth (R):

[0111] d = R·c.

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

[0113] Based on the horizontal distance, the terrain slope of the light section in the snow-free period is known, which is the slope field. The vertical change over this distance can be calculated using the tan function, that is, the calculation method of the snow depth reference error is as follows:

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

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

[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 satellite-borne 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 check, it is necessary to ensure that the snow depth of the site and ICESat-2 belong to the same day;

[0119] The altitude difference between the station and ICESat-2 should not exceed 50m, and they should be kept at the same level as much as possible;

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

[0121] Specific use of Correlation, R 2 ,RMSE,MAE,MedAE and other evaluation indicators. The definitions of each indicator are as follows:

[0122]

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

[0124] In the snow depth point distribution density map of the dataset, the ICESat-2 ATL08 snow depth dataset obtained 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. The snow depth points are densely distributed between 55-70°N, and the snow depth is mostly zero between 30-40°N.

[0125] Figure 5 This is a statistical chart of the number of snow depth estimation deviations. The snow depth estimation deviations are divided into the following five levels: >20cm, 5-20cm, -5-5cm, -20-(-5)cm, <(-20)cm. Figure 5 It can be seen that the number of deviations greater than 20cm accounts for 19.3%, 5-20cm accounts for 11.4%, -5-5cm accounts for 60.5%, -20-(-5)cm accounts for 6.1%, and <(-20)cm accounts for 2.7%. In general, the proportion of snow depth estimation deviations > 5cm accounts for 30.7%, and less than -5cm accounts for 8.8%. The number of overestimates is much greater than the number of underestimates.

[0126] Figure 6 This is a statistical diagram of the snow depth estimation deviation classification in different snow months. The verification data is divided into three parts with the snow depth estimation deviation equal to 5cm as the threshold. Figure 6 The percentage of each part is marked with numbers. During the entire snow accumulation period, the snow depth estimation deviation of -5-5cm accounts for the largest proportion, with the highest in December, reaching 70.2%, and the proportion in March and April is relatively low, reaching 50%. Figure 6 It can be found that in December and January, the proportion of snow depth estimation deviations greater than 5 cm accounted for 28.6% and 27.4% respectively, while the proportion of snow depth estimation deviations less than -5 cm was almost 0. However, as the snow accumulation period progressed, the proportion of snow depth estimation deviations less than -5 cm gradually increased, and in April, the proportion of deviations less than -5 cm exceeded the proportion of deviations greater than 5 cm, reaching 25.5%.

[0127] Figure 7-Figure 8The relationship between snow depth and estimation error in different ranges. The accuracy of ICESat-2 snow depth estimation is related to the thickness of the snow itself. Snow depth in different ranges has different estimation errors. The relationship between the two is as follows: Figure 7 , the gray strip represents the buffer zone along the y=0 line where the snow depth estimation error is within the range of -5-5cm. When the snow depth range is less than 5cm, most of them are overestimated; when the snow depth is 5-40cm, there are both overestimations and underestimates; when the snow depth is greater than 40cm, most of them are underestimated. In order to quantitatively describe the relationship between the two, according to the approximate range of snow depth in the ICESat-2 ATL08 dataset, 0-5cm is defined as thin snow, 5-20cm as moderate snow, 20-40cm as deep snow, and greater than 40cm as heavy snow, and the number of snow depth estimation errors in different ranges is statistically analyzed by level, such as Figure 8 It can be found that the proportion of snow depth estimation errors greater than 20cm, 5-20cm, -5-5cm, -20-(-5)cm and <(-20cm) in the corresponding categories is 0.2%, 14.0%, 85.8%, 0, 0 respectively for thin snow; 2.5%, 20.3%, 66.9%, 10.3%, 0 respectively for moderate snow; 1.2%, 9.6%, 51.4%, 36.9%, 0.9% respectively for deep snow; and 1.3%, 2.7%, 50.4%, 41.4%, 4.2% respectively for heavy snow.

[0128] Figure 9-10 is the relationship between terrain slope and snow depth estimation error. The relationship between slope and snow depth estimation error in the dataset prepared by the present invention is as follows: Fig. 9 It can be found that as the slope increases from 0 to 5°, the snow depth estimation error does not change significantly, and the R 2 0. This shows that the dynamic discrimination rule of the intersection distance threshold defined by the slope in the embodiment has almost no effect on the snow depth estimation when the slope in the data set is less than 5°. In order to further intuitively verify the relationship between the slope and the snow depth estimation error, the experiment retained the light segments with a terrain slope greater than 5°, and selected the light segments with the same other conditions such as the underlying surface of the ICESat-2 light segment, snow accumulation date, snow-free date, intersection distance, distance between the verification point and the station, and altitude. Only the influence of the slope on the snow depth measurement error was considered, as shown in Figure 2. Fig.10 It was found that the absolute error of slope and snow depth estimation was positively correlated, with a Correlation of 0.85 and R 2The error is 0.72, MAE is 8.39cm, and RMSE is 12.58cm. The error increases by about 50cm for every 10-degree increase in slope; when the slope is less than 5°, the error is less than 25cm; when the slope is greater than 20°, the error exceeds 1m.

[0129] Fig.11 The figure is the relationship between altitude and snow depth estimation error. The experiment analyzes the ICESat-2 snow depth estimation at different altitudes. There are 10268 ICESat-2 snow depth points involved in the verification, such as Fig.11 It can be found that at the same altitude, there are different snow depth estimation errors, and the errors do not change significantly with increasing altitude. The correlation between altitude and snow depth estimation error is -0.01, R 2 The correlation between the two is very poor, and the altitude will not affect the estimation of snow depth.

[0130] Fig.12 The figure below is the relationship between the distance between intersections and the relative error of snow depth estimation. In order to strike a balance between the snow depth estimation error and the number of intersections, the experiment adopts a slope-limited intersection distance threshold dynamic discrimination rule, that is, when the slope is ≥2°, the distance is set to 2m; when the slope is <2°, the distance is set to 17m. The foot point diameter of ICESat-2 light spot is 17m, so the center distance between intersections in the dataset should be in the range of 0-17m. The relationship between the distance between snow depth intersections and the snow depth estimation deviation is shown in Figure 2. Fig.12 ,It can be found that as the distance increases, the error changes without regularity, which shows that ,the dynamic discrimination rule of the intersection distance threshold limited by the slope has a good ,effect.

[0131] Fig.13 The relationship between the distance between the station and the intersection and the relative error of the snow depth estimation deviation is discussed. Different colors represent the difference between the station altitude and the snow depth intersection altitude. There are 5 levels with every 10m. The smaller the altitude difference, the snow depth intersection and the station are at the same level, and the station can better represent the snow depth value at the surrounding intersections. Therefore, the smaller the altitude difference, the more reliable the result of the point. Fig.13 It can be seen that as the distance increases, the relative error of the snow depth deviation on the fitting line increases very little, indicating that in most cases, the snow depth of the station can represent the snow depth value of ICESat-2 ATL08 within this range. In the data set produced by the present invention, we have eliminated non-flat areas with a slope greater than 5°. In flat areas, when the difference between the station altitude and the light band altitude is small (less than 50m), the snow depth of the station within 7km can verify the snow depth of the ICESat-2 snow depth intersection point.

[0132] Fig.14The relationship between snow-free months and snow depth estimation errors is discussed. The ICESat-2 snow depth value is the difference between the light periods during the snow accumulation period and the snow-free period. The choice of snow-free dates will directly affect the snow depth estimation, such as Fig.14 . It can be found that the snow depth is generally overestimated. When the snow-free period is closer to the snow accumulation period, the fewer outliers in the data set, such as 2019-11 and 2020-04. Among them, the medians of 2019-10, 2019-11, 2020-04, and 2020-05 are all greater than 0, the mean error is greater than 0, and the error distribution is biased to the right, that is, the snow depth is overestimated; while the medians of 2020-06 and 2020-08 are less than 0, the mean error is less than 0, and the error distribution is biased to the left, that is, the snow depth is underestimated; in August, the northern hemisphere is in summer, vegetation grows vigorously, and the surface elevation is overestimated, which indirectly leads to the most serious underestimation of snow depth. Between 50°N and 70°N, the snow-covered months selected in the dataset are from 2019-12 to 2020-04. 2019-11 and 2020-05 are closest to the snow-covered period, and the boxes of these two months are the flattest, with the most concentrated snow depth errors. Therefore, the smaller the time difference between the non-snow-covered period and the snow-covered period, the smaller the change in surface elevation, and the more accurate the snow depth estimation.

[0133] Fig.15 In order to conduct a comprehensive analysis of the factors that cause snow depth estimation errors, the factors that affect snow depth estimation are, from large to small, slope, underlying surface type, different months of snow accumulation, snow depth range, snow-free months, the distance between the measuring station and the snow depth intersection point, the distance between the snow depth intersection points, and altitude.

[0134] Seasonal snowmelt plays an important role in the global water cycle and water balance. ICESat-2 can be used to prepare a large-scale, high-precision snow depth dataset, which is not limited by factors such as altitude, climate, and cost. It can make up for the sparse snow depth data of stations and provide a reference for the construction of snow models and the measurement of snow depth in high-altitude mountainous areas. This embodiment uses the ICESat-2ATL08 product, and through the relationship between the terrain and the geometric structure of the light segment, the light segment elevation in the snowy period and the snowless period is subtracted to obtain the snow depth value, and the seasonal snow depth dataset of the Northern Hemisphere from 2018 to 2020 is obtained.

[0135] The ICESat-2 ATL08 snow depth dataset was verified using snow depths at more than 43,000 stations in the Northern Hemisphere. It was found that slope is the main factor causing snow depth estimation errors. When the slope is less than 5°, the estimation result is relatively reliable. The greater the slope, the higher the error. For every 10° increase in slope, the error increases by about 50 cm. The snow depth estimation results on different underlying surfaces are significantly different. On bare land, herbaceous vegetation, herbaceous wetlands, cultivated land, shrubs, open forests, and closed forests, R 2They are 0.88, 0.79, 0.70, 0.68, 0.51, 0.21, and 0.18 respectively; the underlying surface with less vegetation and smaller elevation changes during the freeze-thaw period is more suitable for ICESat-2 to estimate snow depth. Different snow accumulation months will have an impact on snow depth estimation. During the accumulation period, the snow is shallower, and the slight elevation difference caused by terrain and underlying surface differences is often greater than the snow depth, which is mainly manifested as overestimation; during the ablation period, the snow particle size is larger and the gap is wider. The laser signal may penetrate part of the snow layer, causing the snow depth to be underestimated. The ICESat-2 snow depth dataset produced by the present invention has strong spatial heterogeneity, and the snow depth may be overestimated or underestimated at different points in a small area.

[0136] Embodiment 2

[0137] The present invention also provides a hemispherical scale snow depth estimation system based on a satellite-borne laser radar, which is used to implement a method, including:

[0138] A data acquisition module is used to download snow accumulation related data within a preset hemispherical scale collected by the satellite-borne laser radar and pre-process the snow accumulation related data;

[0139] A light segment division module is used to obtain the light segments of the snow accumulation period and the light segments of the snow-free period within a preset hemispherical scale based on the preprocessed snow accumulation related data;

[0140] The intersection acquisition module is used to obtain the intersection of the light segment during the snow accumulation period and the light segment during the snow-free period by adopting a dynamic discrimination rule of the intersection distance threshold defined by the slope; wherein the dynamic discrimination rule of the intersection distance threshold defined by the slope includes: setting a threshold for discriminating the center distance between intersections based on a preset terrain slope threshold within a preset hemispherical scale;

[0141] A snow depth value 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 to make a difference between the elevation values ​​of the light segments during the snow accumulation period and the snow-free period at the intersection point;

[0142] An error model building module is used to build a snow depth reference error model based on the terrain slope of the snow-free period light segment and the horizontal distance between the snow-covered period light segment and the snow-free period light segment when forming 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 satellite-borne lidar.

[0144] The embodiments described above are only descriptions of the preferred embodiments of the present invention and are not intended to limit the scope of the present invention. Without departing from the design spirit of the present invention, various modifications and improvements made to the technical solutions of the present invention by ordinary technicians in this field should all fall within the protection scope determined by the claims of the present invention.

Claims

1. A hemispherical-scale snow depth estimation method based on satellite-borne lidar is characterized by: The method comprises: Downloading snow accumulation related data within a preset hemispherical scale collected by a satellite-borne laser radar, and preprocessing the snow accumulation related data; Based on the pre-processed snow-related data, the light segments of the snow-covered period and the snow-free period within the preset hemispherical scale are obtained; The intersection point of the light segment in the snow accumulation period and the light segment in the snow-free period is obtained by adopting a dynamic discrimination rule of the intersection point distance threshold defined by the slope; wherein the dynamic discrimination rule of the intersection point distance threshold defined by the slope comprises: setting a threshold for discriminating the center distance between intersection points based on a preset terrain slope threshold within a preset hemispherical scale; Based on the intersection and the snow depth estimation model, a snow depth value within a preset hemispherical scale is obtained; wherein the snow depth estimation model is to make a difference between the elevation values ​​of the light segment during the snow accumulation period and the light segment during the snow-free period at the intersection; A snow depth reference error model is constructed based on the terrain slope of the snow-free period light segment and the horizontal distance between the snow-covered period light segment and the snow-free period light segment when the intersection is formed; The snow depth value is corrected based on the snow depth reference error model to obtain a corrected snow depth value, thereby completing hemispherical-scale snow depth estimation based on satellite-borne laser radar.

2. The method according to claim 1, characterized in that The snow accumulation 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 accumulation related data includes: The light segments with a slope greater than 5° and the light segments of water bodies, buildings and ice surfaces in the snow accumulation related data are eliminated.

3. The method according to claim 1, characterized in that The snow accumulation period light period includes: the snow accumulation months in the light period determined at 30°N-55°N are December, January, February, and March, a total of 4 months; the snow accumulation months in the light period determined at 55°N-70°N are December, January, February, March, and April, a total of 5 months.

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

2. When the NDSI of the MOD10A1 snow product is less than 0.2, it is judged as no snow cover, otherwise it is judged as snow cover; The snow cover map was imported into ArcGIS software and divided into 5°×5° grids to obtain the snow cover situation on each 5° pixel within the preset date. The date was changed to obtain the first and last days of snow cover for each pixel. Based on the first day of snow accumulation and the last day of snow accumulation, a snow-free period is obtained; The snow-free period time period is compared with the segment_landcover field of the corresponding time period in the pre-processed snow accumulation related data. If the segment_landcover field of the corresponding time period in the pre-processed snow accumulation related data is 1, the snow-free period time period is correctly judged, and the snow-free period light segment is obtained; if the segment_landcover field is 2, the snow-free period time period is incorrectly judged, then the snow-free date is changed and the segment_landcover field comparison is repeated until the segment_landcover field is 1, and the snow-free period light segment is obtained; wherein, segment_landcove is a snow cover identifier, a segment_landcover field of 1 represents snow-free land, and a segment_landcover field of 2 represents snowy land.

5. The method according to claim 1, characterized in that The dynamic discrimination rules of intersection distance threshold based on slope limitation specifically include: 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 It also includes screening the elevation difference between the light segments in the snow accumulation period and the snow-free 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 eliminated. When the segment_landcover fields of the light segments in the snow accumulation period and the snow-free period are both 0, 1, and 3, the elevation difference calculation is not performed, and the snow depth value is directly set to 0. Among them, the segment_landcover field of 0 represents water body, and the segment_landcover field of 3 represents ice surface.

7. The method according to claim 1, characterized in that The method for obtaining the horizontal distance between the light segment during the snow accumulation period and the light segment during the snow-free period when the intersection is formed includes: Convert the longitude and latitude corresponding to the intersection point into radians, and calculate the radian differences between the longitude and latitude of the light segment during the snow-covered period and the light segment during the snow-free period at the intersection point; Based on the arc difference, the initial distance between the light segment in the snowy period and the light segment in the snowless period is obtained by the Haversine formula; The great circle angle is calculated based on the initial distance between the light segment during the snow-covered period and the light segment during the snow-free period; Based on the great circle angle and the radius of the earth, the horizontal distance between the light segment during the snow accumulation period and the light segment during the snow-free period is obtained.

8. A hemispherical-scale snow depth estimation system based on satellite-borne laser radar, used to implement the method according to any one of claims 1 to 7, characterized in that: include: A data acquisition module, used to download snow accumulation related data within a preset hemispherical scale collected by a satellite-borne laser radar, and pre-process the snow accumulation related data; A light segmentation module is used to obtain the light segments of the snow accumulation period and the light segments of the snow-free period within a preset hemispherical scale based on the preprocessed snow accumulation related data; The intersection acquisition module is used to obtain the intersection of the light segment in the snow accumulation period and the light segment in the snow-free period by adopting the dynamic discrimination rule of the intersection distance threshold defined by the slope; wherein the dynamic discrimination rule of the intersection distance threshold defined by the slope includes: setting the threshold for discriminating the center distance between intersections based on the preset terrain slope threshold within the preset hemispherical scale; A snow depth value 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 to make a difference between the elevation values ​​of the light segment during the snow accumulation period and the light segment during the snow-free period at the intersection point; An error model building module is used to build a snow depth reference error model based on the terrain slope of the snow-free period light segment and the horizontal distance between the snow-covered period light segment and the snow-free period light segment 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 to obtain the corrected snow depth value and complete the hemispherical scale snow depth estimation based on the satellite-borne laser radar.

Images