A glacier elevation change monitoring method and device based on a heterogenous laser radar intersection

Abstract

This invention discloses a method and device for monitoring glacier elevation changes based on the intersection points of heterogeneous lidar systems. The method includes: acquiring altimetry data from first and second lidar systems with different orbital configurations; calculating a system deviation correction constant based on the difference in representative statistics of a stable reference surface to unify the elevation benchmarks of the two types of data; extracting elevation changes from intersection points and heterogeneous systems within a preset geographic grid using a spatial cross-pairing mechanism; obtaining prior measurement accuracy estimates for various intersection points; and constructing an inverse variance weighted model to fuse the elevation changes. This invention achieves deep integration of multi-source altimetry data through benchmark collaborative calibration, spatial cross-pairing, and a weighted fusion model, significantly improving the spatial coverage and data reliability of glacier elevation change monitoring.

Description

Technical Field

[0001] This invention belongs to the field of satellite remote sensing and geodesy technology, specifically relating to a method and device for monitoring glacier elevation changes based on the intersection of heterogeneous lidar sources. Background Technology

[0002] Glacier elevation changes are a key indicator for assessing glacier mass balance, predicting sea-level rise, and forecasting regional water resource changes. Satellite lidar technology, particularly the ICESat-2 satellite with its photon-counting laser altimetry system, can provide surface elevation data with sub-meter vertical accuracy, and has become an important tool for monitoring glacier changes. Through repeating orbit self-crosspoint analysis of ICESat-2, the time-varying sequence of glacier surface elevation can be obtained.

[0003] However, ICESat-2 employs a 91-day precisely repeating orbit design with an interval of approximately 3 kilometers between adjacent ground tracks. In complex glacier regions such as the southeastern Tibetan Plateau (SETP), characterized by steep terrain and frequent cloud and snow cover, this design results in an extremely limited number of available crossover points and systematic spatial cover gaps in glacier terminus regions with intense mass exchange. The sampling sparsity of a single data source limits the reliability of regional glacier change characterization.

[0004] To overcome the aforementioned sampling bottlenecks, existing technologies are beginning to explore the combined application of multi-source lidar data. For example, the GEDI instrument, aboard the International Space Station (ISS), uses a non-sunsynchronous orbit with a quasi-random ground trajectory, complementing the fixed orbit of ICESat-2 in space. Previous research by the inventors (Remote Sensing, 2025) based on experiments in the SETP glacier region showed that while the accuracy of GEDI data at single points is limited in steep terrain areas, its dense spatial sampling can theoretically effectively fill the coverage gaps of ICESat-2.

[0005] However, the direct fusion of multi-source and heterogeneous satellite altimetry data presents the following key technical challenges that urgently need to be addressed:

[0006] First, ICESat-2 and GEDI differ fundamentally in sensor design, measurement mechanism, and orbit characteristics. The systematic elevation deviation between the two exhibits spatiotemporal heterogeneity, and existing technologies lack a reliable unified benchmark method for complex terrains such as glaciers.

[0007] Secondly, the geometric characteristics of fixed repeating orbits (such as ICESat-2) and quasi-random orbits (such as GEDI) are completely different. There is still a lack of a unified and automated computational framework for spatial pairing, intersection extraction and quality control of heterogeneous intersection points.

[0008] Finally, the measurement accuracy of different types of intersections (self-intersections and heterogeneous intersections) varies significantly. If simple averaging or equal-weighted fusion is used, obvious local noise will be introduced. Existing technologies still lack a scientific weighting mechanism that can take into account both the coverage of multi-source data and the fusion accuracy.

[0009] In summary, existing technologies for monitoring glacier elevation changes in complex terrain areas still face the dual challenges of insufficient spatial coverage and difficulty in achieving high accuracy in multi-source data fusion. A comprehensive technical solution has yet to be developed that can systematically address heterogeneous sensor bias correction, unified intersection point extraction, and multi-precision data fusion. Summary of the Invention

[0010] To address the aforementioned technical problems, this invention provides a method and device for monitoring glacier elevation changes based on the intersection points of heterogeneous lidar systems. This invention employs a heterogeneous benchmark collaborative calibration strategy, utilizing the difference in representative statistics at stable reference surfaces to accurately obtain the system deviation correction constant, thus constructing a unified elevation benchmark path for heterogeneous lidar altimetry data. Simultaneously, it utilizes an orbital space cross-pairing mechanism to transform discrete footprint point sequences into continuous line segment sequences to capture self-intersections and heterogeneous intersection points. Furthermore, it introduces an inverse variance weighted fusion model based on the accuracy estimate of stable terrain prior measurements. By characterizing the elevation change rate of different types of intersection points with precision and dynamically allocating weights, it achieves adaptive correction of heterogeneous satellite altimetry errors and deep integration of multi-source data, thereby significantly improving the spatiotemporal continuity and data reliability of elevation change monitoring in complex glacier terrain.

[0011] A method for monitoring glacier elevation changes based on the intersection of heterogeneous lidar sources includes the following steps:

[0012] S100 acquires first and second lidar altimetry data of the target glacier region, wherein the first and second lidar altimetry data have different orbital configuration characteristics.

[0013] S200 extracts representative observation values ​​of the first lidar altimetry data and the second lidar altimetry data based on a stable reference surface, calculates and obtains the system deviation correction constant, and uses the system deviation correction constant to unify the elevation benchmark of the second lidar altimetry data.

[0014] Based on the first and second lidar altimetry data after the unified elevation benchmark, S300 performs spatial cross-pairing within a preset geographic grid, extracts the elevation estimate at the intersection point, and calculates the elevation change of various intersection points.

[0015] S400 obtains the prior measurement accuracy estimates of the various intersection points on the stable reference surface, and constructs an inverse variance weighted model based on the prior measurement accuracy estimates; the inverse variance weighted model is used to fuse the elevation changes of the various intersection points in the target glacier area to obtain the fused elevation change results of the target glacier area.

[0016] Preferably, in step S200, the specific method for calculating and obtaining the system deviation correction constant includes:

[0017] S201 extracts bare ground pixels within the target glacier area as the stable reference surface and removes footprint data with local slopes greater than or equal to a preset slope threshold.

[0018] S202 divides the stable reference surface into a regular grid of a preset scale. For each effective grid that contains both the first lidar altimeter data and the second lidar altimeter data, the representative statistical difference between the two types of observations is calculated.

[0019] S203 extracts a global robust statistic from the differences in the representative statistics of all valid grids, which is used as the system bias correction constant.

[0020] Preferably, the representative statistical difference is the difference between the median elevation of the second lidar altimetry data and the median elevation of the first lidar altimetry data within the effective grid; the global robust statistic is the global median of the representative statistical differences across all effective grids.

[0021] Preferably, step S300 specifically includes:

[0022] S301 divides the target glacier area into a geographic grid of a preset scale and determines whether the trajectory of the first lidar altimeter data and the trajectory of the second lidar altimeter data share the same geographic grid.

[0023] S302 connects the discrete footprint point sequences of each track into a continuous line segment sequence within the shared grid, and performs exhaustive spatial intersection calculation on the continuous line segment sequences from different tracks to obtain spatial intersection points.

[0024] S303 determines whether the distance along the track from the spatial intersection point to the nearest real footprint point and the local slope meet the preset threshold conditions; if they meet the conditions, the two footprint points closest to the spatial intersection point are taken along each track and distance-weighted linear interpolation is performed to obtain the elevation estimate at the intersection point, and the difference is used to obtain the elevation change.

[0025] Preferably, the various intersection points include three types: a first self-intersection point formed by the first lidar height measurement data, a second self-intersection point formed by the second lidar height measurement data, and a heterogeneous intersection point formed by the first lidar height measurement data and the second lidar height measurement data.

[0026] Preferably, in step S400, the method for constructing the inverse variance weighted model specifically includes:

[0027] S401 extracts the elevation change rate of various intersection points on the stable reference surface; calculates the normalized median absolute deviation of the elevation change rate of each type of intersection point, and uses it as an estimate of the prior measurement accuracy of that type of intersection point. ;

[0028] S402 divides the target glacier region into multiple independent analysis units, and for the various intersection points falling within the same analysis unit, a formula is used. Construct an inverse variance weighting scheme for weight allocation.

[0029] Preferably, the first lidar altimeter data is ICESat-2 satellite data with fixed repeating orbit characteristics, and the second lidar altimeter data is GEDI data with non-sun-synchronous or non-repeating orbit characteristics; the preset slope threshold is 5°, the size of the regular grid is 100 m × 100 m, and the size of the geographic grid is 0.1° × 0.1°.

[0030] The present invention also provides a glacier elevation change monitoring device, comprising:

[0031] The data acquisition module is used to acquire first lidar altimeter data and second lidar altimeter data, wherein the first lidar altimeter data and the second lidar altimeter data have different orbital configuration characteristics.

[0032] The benchmark correction module is used to extract representative observation values ​​of the first lidar altimetry data and the second lidar altimetry data based on a stable reference surface, calculate and obtain the system deviation correction constant, and use the system deviation correction constant to unify the elevation benchmark of the second lidar altimetry data.

[0033] The intersection point extraction module is used to perform spatial cross-pairing within a preset geographic grid based on the first lidar altimetry data and the second lidar altimetry data after a unified elevation benchmark, extract the elevation estimate at the intersection point, and calculate the elevation change of various intersection points.

[0034] The weighted fusion module is used to obtain the prior measurement accuracy estimates of the various intersection points on the stable reference surface, construct an inverse variance weighted model based on the prior measurement accuracy estimates, and use the inverse variance weighted model to fuse the elevation changes of the various intersection points in the target glacier area to obtain the fused elevation change results of the target glacier area.

[0035] The above technical solution has the following advantages or beneficial effects:

[0036] (1) This invention overcomes the problem of sparse sampling under complex terrain and snowy conditions by introducing second lidar data from non-sunsynchronous / quasi-random orbits and constructing a spatial cross-pairing mechanism for continuous line segment sequences. Experimental data show that, compared with a single data source, this invention increases the total number of fused cross points by about 43.5 times, expands the effective coverage area by about 6.1 times, and increases the spatial sampling density to 5.75 points / km², effectively filling the observation gaps in key areas such as glacier termini.

[0037] (2) This invention uses low-slope stable bare land as a reference surface and extracts the difference of representative statistics (median) as a system bias correction constant. This benchmark collaborative calibration strategy effectively suppresses the error amplification caused by local terrain, making the internal consistency of the annual average elevation change rate calculated from different intersection types better than 0.07 m / yr, thus ensuring the reliability of heterogeneous data fusion.

[0038] (3) In the fusion stage, this invention constructs an inverse variance weighted model based on prior measurement accuracy. Weights are dynamically assigned for different types of intersection points, effectively suppressing local noise introduced by low-precision data sources. Verified with independent gravity satellite data, the regional mass loss rate estimated by this invention is highly consistent with the benchmark data within the statistical error range (-5.80 ± 0.41 Gt / yr vs. -5.81 ± 0.42 Gt / yr), ensuring high fidelity while expanding the sample size.

[0039] (4) The “system deviation correction-unified cross pairing-inverse variance weighted fusion” technical framework constructed by this invention does not depend on specific parameters of specific sensors. For monitoring glacier elevation changes, it can be applied to a variety of spaceborne altimetry mission combinations with complementary orbital configurations, providing a standardized data processing scheme. Attached Figure Description

[0040] The accompanying drawings, which are incorporated in and form part of this specification, illustrate embodiments consistent with the invention and, together with the description, serve to explain the principles of the invention, wherein:

[0041] Figure 1 This is a flowchart illustrating the glacier elevation change monitoring method of Embodiment 1 of the present invention.

[0042] Figure 2 This is a schematic diagram of the intersection point extraction algorithm in Embodiment 1 of the present invention;

[0043] Figure 3 This is a spatial coverage pattern diagram of the three intersection types in the SETP glacier region according to Embodiment 1 of the present invention;

[0044] Figure 4 This is a distribution characteristic diagram of the glacier elevation change rate along the elevation gradient in Embodiment 1 of the present invention; Figure 4 (a)-(d) correspond to the four cases: GEDI-GEDI, IS2-IS2, GEDI-IS2, and inverse variance weighted fusion, respectively. Detailed Implementation

[0045] It should be noted that the following description is merely a preferred embodiment of the present invention and is not intended to limit the invention. Various modifications and variations can be made to the present invention by those skilled in the art. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the scope of protection of the present invention.

[0046] The overall technical approach of this invention is to effectively combine the originally sparse but high-precision ICESat-2 data with the densely sampled but relatively low-precision GEDI data with complementary orbits through a chain-like technical path of "benchmark unification - intersection point generation - accuracy-oriented adaptive fusion". By quantifying accuracy and scientifically weighting, the limitations of a single data source are overcome, and the monitoring range and accuracy are improved simultaneously.

[0047] Example 1

[0048] like Figure 1 As shown, this embodiment of the invention provides a method for monitoring glacier elevation changes based on the intersection of heterogeneous lidar sources, specifically including the following steps:

[0049] S100: Acquires altimeter data from heterogeneous lidar and performs multi-level quality control.

[0050] In this embodiment, the heterogeneous lidar altimeter data includes data from the first satellite lidar ICESat-2 ATL06 land ice product (V7 version) and data from the second satellite lidar GEDI L2A product (V2 version).

[0051] The specific preprocessing steps include:

[0052] S101 Primary Screening

[0053] Obtain ICESat-2 ATL06 land ice products and GEDI L2A products with matching time spans. Read the official quality flags of both (e.g., atl06_quality_summary=0 for ATL06, quality_flag=1 for GEDI), covering the time span from 2019 to 2025, and remove invalid observations.

[0054] S102 Cloud Cover and Outlier Filtering

[0055] An external digital elevation model (DEM) – NASADEM – is introduced as a coarse elevation reference to calculate the elevation difference between the laser footprint elevation and the reference DEM. If the absolute vertical deviation exceeds a preset gross error threshold (set at 100 m), it is considered cloud contamination or a severe multipath effect and is discarded.

[0056] S103 Fine Constraints

[0057] Threshold constraints are imposed on the elevation uncertainty and geolocation error of the first satellite lidar data, specifically: elevation uncertainty < 0.8 m, geolocation error < 4.5 m; for the second satellite lidar data, the lowest mode elevation (elev_lowestmode) is uniformly extracted as the height measurement value representing the real ground.

[0058] S200: Deviation correction for heterogeneous sensor systems based on a stable reference surface.

[0059] Due to inherent differences in track configuration, pointing accuracy, and penetration characteristics among different sensors, direct fusion would introduce significant systematic errors. This step aims to unify the elevation datum, and the specific operations are as follows:

[0060] S201 Reference Surface Extraction

[0061] The high-precision land cover product ESA WorldCover was used to extract "bare land pixels" in the study area as a stable reference surface for system bias calibration. A slope mask was then superimposed to retain only the footprint points in areas where the local slope is less than the first preset slope threshold of 5°, so as to eliminate the projection of horizontal positioning error caused by steep slope onto vertical error.

[0062] S202 Gridding Deviation Calculation

[0063] The study area was divided into a regular grid with a predetermined spatial scale of 100 m × 100 m. For each effective grid cell containing observations from both types of sensors, the difference between the medians of the two types of observations within it was calculated. ;

[0064] in, This refers to the local elevation deviation within a single grid. and These are the footprint elevation sets of GEDI and ICESat-2 within the grid, respectively.

[0065] S203 Global Calibration

[0066] The global median of all valid grid cell bias values ​​is extracted as the final system bias correction constant. Apply this correction constant to all raw second sensor data: This achieves physical unification of elevation benchmarks.

[0067] S300: Unified intersection extraction of heterogeneous unions.

[0068] Combination Figure 2 As shown, this step seamlessly extracts three types of intersection points in space. The specific operation is as follows:

[0069] S301 Dimensional Reduction Search Space

[0070] The study area was divided into a regular geographical grid of 0.1° × 0.1° (approximately 10 km × 10 km). The grid cells traversed by each satellite orbit were recorded. A potential intersection was only determined when two different orbits shared the same grid cell, thereby greatly reducing the computational overhead of matching massive amounts of heterogeneous data.

[0071] S302 Exhaustive Pairing and Intersection Calculation

[0072] Within the shared grid, the discrete footprint point sequence of each track is connected into a continuous line segment sequence. Exhaustive spatial intersection calculations are performed on all line segments within the two tracks (track A and track B) to obtain the spatial intersection points on the two-dimensional plane.

[0073] S303 Nearest Neighbor Interpolation and Quality Check

[0074] Using the extracted spatial intersection point as the center, calculate the distance along the track from it to the nearest real footprint point on both tracks. Set an effective distance threshold of 100m and a local slope threshold. A rigorous screening process is performed. Once the criteria are met, distance-weighted linear interpolation is performed on the two footprints closest to the intersection point along each track to obtain two elevation estimates at the intersection point. and The difference That is, the intersection point during the time period Elevation changes within the area.

[0075] The algorithm described above synchronously outputs three types of cross-point sets: ICESat-2 self-cross-points (IS2-IS2), GEDI self-cross-points (GEDI-GEDI), and ICESat-2 / GEDI heterogeneous cross-points (GEDI-IS2).

[0076] S400: Inverse variance weighted fusion based on robust accuracy estimation.

[0077] Because the data sources at the three types of intersections are of different quality, they cannot be directly weighted and averaged.

[0078] S401 Accuracy Prior Estimation

[0079] Again utilizing the intersections in stable bare land areas, assuming the theoretical value of bare land elevation change is 0, the normalized median absolute deviation (NMAD) of the annual change rate measured in bare land for each intersection type is calculated as a priori estimate of the measurement accuracy for that type of intersection. .

[0080] Under stable bare ground conditions, the measurement accuracy of different intersection types is dominated by sensor noise and track geometry, and their relative accuracy relationships are consistent in glacial regions. Therefore, the prior accuracy of bare ground can be directly used to guide the weight allocation of intersections in glacial regions.

[0081] In this embodiment, the measured It is 0.33 m / yr. It is 0.98 m / yr. It is 1.78m / yr.

[0082] S402 Inverse Variance Weighted

[0083] A regular geographic grid is used to spatially partition the target glacier region, with each grid serving as an independent analysis unit. For the three types of intersections falling within the same analysis unit (ICESat-2 self-intersection, GEDI self-intersection, and ICESat-2 and GEDI heterogeneous intersection), a weighted fusion is performed using the inverse square of the standard deviation (i.e., inverse variance) of the elevation change residuals of each intersection point as the weight, calculating the representative elevation change rate of the grid. Using the GEDI self-intersection point with the lowest relative accuracy as the weight reference (its weight is set to 1), the optimal fusion weight ratio for the three types of intersections is IS2-IS2 : GEDI-GEDI : IS2-GEDI ≈ 30 : 1 : 3.3. This weighting strategy effectively incorporates the advantages of GEDI's dense spatial coverage while minimizing noise pollution caused by its lower single-point accuracy.

[0084] The aforementioned weighting ratio is calculated based on prior accuracy estimates obtained from bare land areas (where the residual standard deviations of IS2-IS2, GEDI-IS2, and GEDI-GEDI are approximately σ1, σ2, and σ3, respectively, and are weighted by inverse variance, i.e., weighting ratio 1 / σ1² : 1 / σ2² : 1 / σ3²). This ratio can be adaptively adjusted according to different sensor combinations or study areas, but the weighted fusion method of this invention does not depend on this specific value.

[0085] S500: Spatiotemporal evolution analysis of glacier elevation changes.

[0086] Based on the fused intersection elevation change rate data ( Combining the glacier outline boundaries, the average mass balance of each glacier and even the entire region is calculated, and a model is generated. Figure 4 The diagram shows the distribution characteristics of the glacier elevation change rate along the altitudinal gradient, thus enabling a quantitative inversion of the glacier evolution law.

[0087] Comparative Example 1

[0088] To further verify the technical effects and substantial progress of the technical solution of the present invention, the following comparative examples are set up for comparison with Embodiment 1 of the present invention.

[0089] Test conditions and experimental setup:

[0090] Test area: Southeastern Tibetan Plateau (SETP) glacier region.

[0091] Time span: 2019 to 2025.

[0092] Comparative Example 1: The conventional single-sensor cross-point method in the existing technology is adopted, which relies only on ICESat-2 data, uses precise repeating track self-crossing extraction logic, does not introduce heterogeneous data, and does not involve multi-source weighted fusion processing.

[0093] Example 1: Using the method provided in Example 1 of this invention, ICESat-2 and GEDI heterogeneous data are combined, and heterogeneous system bias correction, unified extraction of three types of cross points, and inverse variance weighted fusion are performed in sequence.

[0094] The comparative test results and analysis are as follows:

[0095] 1. Comparative Analysis of Data Sampling Density and Spatial Coverage

[0096] In Comparative Example 1, limited by the orbital span of a single ICESat-2 sensor and the combined effects of high-altitude cloudy weather conditions, only 20,075 usable intersection points were obtained, covering only 160 large glaciers in the region, with an effective coverage area of ​​805.51 km², resulting in an overall spatial point density of approximately 0.13 points / km². Particularly in the glacier termini (ablation zone) where material dynamics change drastically, there are systematic spatial observation blind spots, making it difficult to provide spatially representative observation samples.

[0097] In Example 1, by constructing a heterogeneous joint architecture and introducing GEDI-IS2 heterogeneous intersection points, the sampling bottleneck of a single sensor was effectively overcome. Data shows that the total number of intersection points after fusion increased to 872,714 (approximately 43.5 times that of Comparative Example 1); the number of covered glaciers increased to 3,551 (an increase of approximately 22.2 times); the effective coverage area expanded to 4,895.00 km² (an increase of approximately 6.1 times); and the overall spatial point density increased to approximately 5.75 points / km². These results objectively demonstrate that the method of this invention significantly fills the observation gap at glacier terminals in complex terrain, achieving near-full coverage monitoring of regional glacier elevation changes.

[0098] 2. Cross-validation analysis of measurement accuracy and independent observation data

[0099] To assess the impact of introducing heterogeneous non-repeating orbit data on the overall measurement accuracy, the elevation change rate extracted in Example 1 was converted into regional glacier mass balance and compared with the concurrent independent observation results from the Gravity Recovery and Climate Experiment Satellite (GRACE / GRACE-FO) for verification.

[0100] Calculations showed that the glacier mass loss rate in the SETP region, obtained based on the method in Example 1, was -5.80 ± 0.41 Gt / yr; while the concurrent observation results of GRACE / GRACE-FO, serving as an independent validation benchmark, were -5.81 ± 0.42 Gt / yr. The two results are highly consistent within the statistical error range. This validation result demonstrates that, while increasing the sample size tenfold and significantly improving spatial coverage, this invention effectively controls the systematic errors and local noise that may be introduced by heterogeneous data through bare land systematic bias correction and inverse variance weighting mechanisms, ensuring high fidelity and high accuracy in regional-scale glacier mass balance assessment.

[0101] In summary, the present invention not only overcomes the technical defects of existing single altimetry technology in complex terrain areas with sparse sampling and insufficient spatial representativeness, but also achieves an order-of-magnitude increase in the number of monitoring samples and coverage without sacrificing vertical measurement accuracy.

[0102] The above description is merely a preferred embodiment of the present invention and does not limit the scope of patent protection of the present invention. Any equivalent structural or procedural transformations made based on the description and drawings of the present invention, or direct or indirect applications in other glacier evolution assessment and monitoring scenarios, are similarly included within the scope of patent protection of the present invention.

Claims

1. A method for monitoring glacier elevation changes based on the intersection of heterogeneous lidar sources, characterized in that, Includes the following steps: S100 acquires first and second lidar altimetry data of the target glacier region, wherein the first and second lidar altimetry data have different orbital configuration characteristics. S200 extracts representative observation values ​​of the first lidar altimetry data and the second lidar altimetry data based on a stable reference surface, calculates and obtains the system deviation correction constant, and uses the system deviation correction constant to unify the elevation benchmark of the second lidar altimetry data. Based on the first and second lidar altimetry data after the unified elevation benchmark, S300 performs spatial cross-pairing within a preset geographic grid, extracts the elevation estimate at the intersection point, and calculates the elevation change of various intersection points. S400 obtains the prior measurement accuracy estimates of the various intersection points on the stable reference surface, and constructs an inverse variance weighted model based on the prior measurement accuracy estimates; the inverse variance weighted model is used to fuse the elevation changes of the various intersection points in the target glacier area to obtain the fused elevation change results of the target glacier area.

2. The elevation change monitoring method according to claim 1, characterized in that, In step S200, the specific method for calculating and obtaining the system deviation correction constant includes: S201 extracts bare ground pixels within the target glacier area as the stable reference surface and removes footprint data with local slopes greater than or equal to a preset slope threshold. S202 divides the stable reference surface into a regular grid of a preset scale. For each effective grid that contains both the first lidar altimeter data and the second lidar altimeter data, the representative statistical difference between the two types of observations is calculated. S203 extracts a global robust statistic from the differences in the representative statistics of all valid grids, which is used as the system bias correction constant.

3. The elevation change monitoring method according to claim 2, characterized in that, The representative statistical difference is the difference between the median elevation of the second lidar altimetry data and the median elevation of the first lidar altimetry data within the effective grid; the global robust statistic is the global median of the representative statistical differences across all effective grids.

4. The elevation change monitoring method according to claim 1, characterized in that, Step S300 specifically includes: S301 divides the target glacier area into a geographic grid of a preset scale and determines whether the trajectory of the first lidar altimeter data and the trajectory of the second lidar altimeter data share the same geographic grid. S302 connects the discrete footprint point sequences of each track into a continuous line segment sequence within the shared grid, and performs exhaustive spatial intersection calculation on the continuous line segment sequences from different tracks to obtain spatial intersection points. S303 determines whether the distance along the track from the spatial intersection point to the nearest real footprint point and the local slope meet the preset threshold conditions; if they meet the conditions, the two footprint points closest to the spatial intersection point are taken along each track and distance-weighted linear interpolation is performed to obtain the elevation estimate at the intersection point, and the difference is used to obtain the elevation change.

5. The elevation change monitoring method according to claim 1 or 4, characterized in that, The various types of intersections include three types: a first self-intersection formed by the first lidar altimeter data, a second self-intersection formed by the second lidar altimeter data, and a heterogeneous intersection formed by the first lidar altimeter data and the second lidar altimeter data.

6. The elevation change monitoring method according to claim 1, characterized in that, In step S400, the method for constructing the inverse variance weighted model specifically includes: S401 extracts the elevation change rate of various intersection points on the stable reference surface; calculates the normalized median absolute deviation of the elevation change rate of each type of intersection point, and uses it as an estimate of the prior measurement accuracy of that type of intersection point. ; S402 divides the target glacier region into multiple independent analysis units, and for the various intersection points falling within the same analysis unit, a formula is used. Construct an inverse variance weighting scheme for weight allocation.

7. The elevation change monitoring method according to any one of claims 1 to 6, characterized in that, The first lidar altimeter data is ICESat-2 satellite data with fixed repeating orbit characteristics, and the second lidar altimeter data is GEDI data with non-sun-synchronous or non-repeating orbit characteristics; the preset slope threshold is 5°, the size of the regular grid is 100 m × 100 m, and the size of the geographic grid is 0.1° × 0.1°.

8. A device for monitoring changes in glacier elevation, characterized in that, include: The data acquisition module is used to acquire first lidar altimeter data and second lidar altimeter data, wherein the first lidar altimeter data and the second lidar altimeter data have different orbital configuration characteristics. The benchmark correction module is used to extract representative observation values ​​of the first lidar altimetry data and the second lidar altimetry data based on a stable reference surface, calculate and obtain the system deviation correction constant, and use the system deviation correction constant to unify the elevation benchmark of the second lidar altimetry data. The intersection point extraction module is used to perform spatial cross-pairing within a preset geographic grid based on the first lidar altimetry data and the second lidar altimetry data after a unified elevation benchmark, extract the elevation estimate at the intersection point, and calculate the elevation change of various intersection points. The weighted fusion module is used to obtain the prior measurement accuracy estimates of the various intersection points on the stable reference surface, construct an inverse variance weighted model based on the prior measurement accuracy estimates, and use the inverse variance weighted model to fuse the elevation changes of the various intersection points in the target glacier area to obtain the fused elevation change results of the target glacier area.

Images