Glacier mass balance calculation method for generating elevation time sequence based on ASTER stereo data

Abstract

The invention discloses a glacier mass balance calculation method for generating an elevation time sequence based on ASTER stereo data. The glacier mass balance calculation method comprises the following steps of: 1, collecting a glacier ASTER L1A stereo image and preprocessing the glacier ASTER L1A stereo image; 2, the quality of the ASTER stereoscopic image pair is improved through radiation correction and strip removal processing; 3, calculating a rational polynomial coefficient RPC model of the ASTER stereo image pair, and detecting and correcting a jitter error of a satellite in a vertical orbit direction; 4, performing residual deviation correction by using the reference ASTER DEM data; 5, high-precision ASTER DEM data are obtained through abnormal value filtering; 6, generating time-continuous ASTER DEM data by adopting a Gaussian process regression method; and 7, calculating the volume change of each glacier, and polymerizing the volume change of each glacier to obtain the overall mass balance change of the regional glacier. According to the method, high-precision and time-continuous ASTER DEM data are generated through correcting satellite jitter errors by a system, and abnormal value filtering and statistical aggregation are combined, so that the mass balance change of the large-scale regional glacier is accurately calculated.

Description

Technical Field

[0001] This invention belongs to the field of remote sensing image processing technology, specifically a method for calculating glacier mass balance based on ASTER stereo data to generate elevation time series. Background Technology

[0002] Glaciers are sensitive indicators of climate change, and changes in their mass directly affect regional water supply and global sea-level rise. Glacier mass balance, as an important indicator reflecting the health of glaciers, is of great significance for understanding climate change, water resource management, and sea-level change. Therefore, accurate monitoring of glacier mass balance has important scientific value and practical application significance for assessing the impact of climate change and predicting future water resource change trends.

[0003] Traditional methods for monitoring glacier mass balance primarily rely on field measurements. These involve setting up observation facilities such as snowplows and snow pits on the glacier surface to periodically measure glacier melt and accumulation. However, due to the complex terrain and harsh environment of high-altitude glacier regions, field observations face numerous challenges: First, reaching glacier areas requires overcoming difficulties such as inconvenient transportation and high-altitude hypoxia, resulting in high safety risks and labor costs. Second, the limited number of field observation points makes it difficult to reflect the spatial variation characteristics of the entire glacier region, leading to insufficient representativeness of the observation results. Third, long-term continuous observation requires sustained human and financial investment, which is difficult to maintain in remote areas. Furthermore, for some sparsely populated glacier areas, field observation is almost impossible. These limitations make traditional methods insufficient to meet the needs of large-scale and long-term glacier monitoring.

[0004] In recent years, the rapid development of remote sensing technology has provided an effective means for large-scale glacier monitoring. Satellite remote sensing has advantages such as wide coverage, fixed revisit period, and no terrain limitations, enabling systematic monitoring of global glaciers. By generating digital elevation models (DEMs) from satellite stereo images, information on changes in glacier surface elevation can be obtained, thereby calculating changes in glacier volume and mass. It can monitor a large number of glaciers simultaneously without requiring personnel to enter dangerous areas, and can also retrospectively analyze historical image data for long-term series analysis. In particular, since its launch in 2000, the ASTER satellite has accumulated a large amount of stereo image data, providing valuable data resources for glacier change research.

[0005] However, existing DEM generation methods based on remote sensing imagery face numerous technical challenges in practical applications, limiting the accuracy and reliability of glacier mass balance calculations. Specifically: 1) Satellites inevitably experience attitude jitter during operation. This jitter consists of three components: vertical to the orbital direction, parallel to the orbital direction, and rotation around the nadir axis. In particular, jitter in the vertical orbital direction directly causes the epipolar lines of stereo image pairs to deviate from their theoretical positions, preventing the perspective rays of the two images from intersecting correctly and severely affecting the accuracy of stereo matching. Traditional DEM generation methods often neglect or simplify jitter correction. This results in significant systematic biases in the generated DEM data, especially in mountainous areas with drastic elevation changes. These errors can reach several meters or even tens of meters, severely impacting the accuracy of glacier mass balance calculations. 2) The original DEM data contains numerous outliers and noise due to factors such as cloud cover, seasonal snow cover, low image contrast, and shadow occlusion. Specifically: clouds completely block surface information, leading to matching failures; seasonal snow cover alters surface reflectivity, causing elevation misjudgments; low-contrast areas (such as glacier accumulation zones) lack sufficient texture information, resulting in poor matching reliability; and shadowed mountain areas are particularly difficult to obtain accurate texture information. Effective elevation information is needed; however, if these outliers are not effectively identified and removed, they will introduce huge errors in subsequent analysis. Although existing methods employ some simple threshold filtering techniques, they are often too coarse, either over-filtering, leading to the loss of effective data, or under-filtering, resulting in the retention of outliers, making it difficult to achieve a good balance between data integrity and quality; 3) Due to the limitation of satellite revisit cycles, it is difficult for a single sensor to obtain time-dense observation data. Taking the ASTER satellite as an example, although the ASTER satellite revisit cycle is 16 days, considering factors such as cloud cover and data quality, the actual usable image data is sparse and irregular in time. This problem of insufficient temporal resolution makes it difficult to capture the seasonal change characteristics of glaciers, to distinguish between different processes of glacier ablation and accumulation, and to accurately analyze interannual change trends. Therefore, some studies have attempted to directly use discrete DEM data to calculate the multi-year average rate of change, but this method ignores the nonlinear characteristics and seasonal fluctuations of glacier changes, and is prone to large estimation biases. It is evident that how to reconstruct a time-continuous high-order sequence from time-sparse observation data is a key issue in improving the accuracy of glacier mass balance calculation.

[0006] Furthermore, existing research often lacks a systematic processing workflow when dealing with large-scale remote sensing data. Many studies only improve one or a few technical aspects, lacking a complete technical chain from data preprocessing, error correction, outlier filtering to time series interpolation. For example, although some studies have noticed the satellite jitter problem, they only corrected some error components; some studies have adopted complex interpolation algorithms, but the quality control of the raw data is not strict enough, resulting in the problem of "garbage in, garbage out". In addition, glaciers in different regions have different topographic features and change patterns, and the parameter settings and algorithm selection of existing methods often lack specificity, versatility and adaptability. Moreover, when dealing with thousands of image data, computational efficiency is also an important issue. Although many algorithms have high accuracy, their computational complexity is too high, making it difficult to apply to large-scale monitoring at the regional scale.

[0007] In summary, current methods for calculating glacier mass balance based on remote sensing imagery have significant shortcomings in key technical aspects such as satellite jitter correction, outlier filtering, and time series interpolation. These shortcomings make it difficult to meet the needs of high-precision, large-scale, and long-term glacier monitoring. Therefore, there is an urgent need to develop a complete technical solution that can systematically correct satellite jitter errors, effectively filter outliers, generate high-order sequences with continuous time, and accurately calculate glacier mass balance, so as to provide reliable data support for glacier change research and climate change assessment. Summary of the Invention

[0008] The purpose of this invention is to provide a method for calculating the mass balance of glaciers based on ASTER three-dimensional data to generate elevation time series. By systematically correcting satellite jitter errors, high-precision and time-continuous ASTER DEM data is generated. Combined with outlier filtering and statistical aggregation, the mass balance changes of large-scale regional glaciers can be accurately calculated.

[0009] This invention is achieved through the following technical solution: A method for calculating glacier mass balance based on ASTER 3D data to generate elevation time series includes the following steps: Step 1: Collect ASTER L1A stereo images within the glacier study area. Select images with cloud cover of less than 95% and a time span covering the study period as the original ASTER stereo image pairs for generating ASTER DEM data. Step 2: Perform radiometric correction and destriping on the original ASTER stereo image pairs to obtain high-quality ASTER stereo image pairs. Step 3: Calculate the rational polynomial coefficient RPC model of the high-quality ASTER stereo image pair, and perform two-dimensional correlation analysis based on the theoretical epipolar line. Use a 7th-order polynomial to fit and correct the low-frequency error, and fit the sum of 8 sine curves to correct the high-frequency jitter error. Correct the jitter error in the vertical orbit direction of the satellite through the low-frequency error and high-frequency error to obtain the corrected ASTER stereo image pair and RPC model. Step 4: Calculate the initial ASTER DEM data using the corrected stereo image pair and RPC model. Use the reference ASTER DEM data to correct the residual bias of the initial ASTER DEM data. Then, use quadratic regression to fit the relationship between the elevation deviation and the sensor angle to correct the satellite vertical orbit system error. Finally, use sine / cosine function modeling to gradually eliminate system errors and high-frequency noise to obtain the corrected ASTER DEM data. Step 5: Use glacier contour data to generate a buffer and crop the corrected ASTER DEM data. Use reference ASTER DEM data to filter out outliers. Remove outliers through weighted least squares fitting and Gaussian process regression to obtain high-precision ASTER DEM data. Step 6: Construct a time covariance function using a linear kernel, an exponential sine square kernel, a radial basis function kernel, and a rational quadratic kernel. Then, use the Gaussian process regression method to interpolate the high-precision ASTER DEM data into an elevation time series with a one-month time step to generate time-continuous ASTER DEM data. Step 7: Calculate the volume change of each glacier based on the elevation time series, and then aggregate the volume changes of each glacier to obtain the overall mass balance change of the regional glaciers.

[0010] Furthermore, the specific process of step 1 is as follows: Step 1.1: Based on the geographical and temporal ranges, retrieve all available ASTERL1A stereo image pairs within the glacier study area; Step 1.2: Perform preliminary quality screening on the available ASTER L1A stereo image pairs: select ASTER stereo image pairs with cloud cover of less than 95% and a time span covering the study period. Step 1.3: Stitch together the ASTER stereo image pairs from every three consecutive acquisition days to obtain the original ASTER stereo image pairs used to generate ASTERDEM data.

[0011] Furthermore, the specific process of step 2 is as follows: Step 2.1: Convert the DN values ​​of the original ASTER stereo image pair into sensor radiance using radiometric calibration coefficients. , is represented as: In the formula: It is the radiance of the sensor; It is a numerical value; It is the gain factor; It is the offset; Step 2.2: Measure the sensor radiance. Converted to surface reflectance ， Represented as: In the formula: It is the surface reflectance; It is the distance between the Earth and the Sun; It is solar irradiance; It is the zenith angle of the sun; Step 2.3: Correct the striping effect using sensor response calibration data. The calibration function is a linear function composed of gain coefficient, offset coefficient and quadratic term coefficient, used to reduce the striping effect in the image and obtain high-quality ASTER stereo image pairs.

[0012] Furthermore, the specific process of step 3 is as follows: Step 3.1: Based on the high-quality ASTER stereo image pair, convert the geocentric coordinates to topographic Cartesian coordinates, and then to geodetic coordinates. Set the vertical grid density to 200 meters and define the range between -500 meters and +8850 meters. Step 3.2: Perform two-dimensional correlation analysis around the theoretical epipolar line to analyze the corresponding points of the ASTER stereo image pairs. The vibration intensity in pixels is equal to the distance between the maximum correlation point and the epipolar line, and the influence of the satellite vertical orbit jitter component is directly observed. Step 3.3: Model the parallax error caused by the jitter component in the vertical orbit direction of the satellite, and use a 7th-order polynomial to correct the low-frequency error. , is represented as: In the formula: Represents the fitting coefficients of a 7th-order polynomial. Represents the coordinates of the satellite's parallel orbit direction on the image; Step 3.4: Estimate the high-frequency jitter error by fitting the sum of eight sine curves along the axis on a 1000-pixel wide column with 90% overlap. , is represented as: In the formula: Representing the The amplitude of a sine curve represents the intensity of the jitter component at that frequency; Representing the The wavelength of a sine curve represents the spatial period of the frequency jitter component. Representing the The phase of a sine curve represents the spatial starting position of the jitter component at that frequency. The coordinates representing the direction of the satellite's parallel orbit in the image; Step 3.5: Calculate the results... and As correction value ， Applied to rear-view images, bilinear interpolation is used to resample the rear-view images to obtain corrected ASTER stereo image pairs and RPC models.

[0013] Furthermore, the specific process of step 4 is as follows: Step 4.1: Calculate the initial ASTER DEM data using the corrected ASTER stereo image pairs and the RPC model; Step 4.2: Calculate the elevation deviation between the generated initial ASTER DEM data and the reference ASTER DEM data in non-glacier areas. , is represented as: In the formula: This represents the initial ASTER DEM data; This indicates that the reference data is ASTER DEM data; Step 4.3: Use quadratic regression to fit the elevation deviation. The relationship with the sensor angle is used to correct errors in the satellite's vertical orbit system. , is represented as: In the formula: The sensor angle indicating the vertical orbital direction of the satellite; , and These are the coefficients of the quadratic regression model; Step 4.4: Fitting the low-frequency error of the satellite's parallel orbit direction , is represented as: In the formula: The coordinates of the satellite's parallel orbit direction on the image. ， 、 、 These are the coefficients of the quadratic regression model; Step 4.5: Use sine / cosine function fitting to remove high-frequency errors. , is represented as: In the formula: Indicates the first The wavelength of each frequency component; Indicates the first The amplitude coefficient of the sinusoidal term of each frequency component; It is the first The amplitude coefficient of the cosine term of each frequency component; It is an index of the frequency component; These are the coordinates of the satellite's parallel orbit direction on the image; Step 4.6, use , and The initial ASTER DEM data was progressively corrected to obtain the corrected ASTER DEM data.

[0014] Furthermore, the specific process of step 5 is as follows: Step 5.1: Use glacier contour data to generate a 10 km buffer zone around each glacier, and then crop and correct the ASTER DEM data based on the buffer zone. Step 5.2: Exclude corrected ASTER DEM data with a root mean square error greater than 20 meters from the reference ASTER DEM data on ice-free terrain, and retain corrected ASTER DEM data that meet the requirements. Step 5.3, with radius The circle represents the filtering window, excluding data where the absolute height difference between the corrected ASTER DEM data and the reference ASTER DEM data exceeds a threshold. The pixels are represented as: In the formula: The threshold representing outlier filtering; Representing the In this iteration, the corrected ASTER DEM data is in pixels. The elevation value of the location; Represents reference ASTER DEM data at the pixel level The elevation value of the location; Step 5.4: Considering terrain slope and stereo correlation quality, calculate the elevation measurement error of each pixel selected and retained in Step 5.3. , is represented as: In the formula: Indicates elevation measurement error; Indicates the slope of the terrain; Indicates the stereo correlation quality; This represents the registration error. Ideally, the registration error is calculated on pixels with a small slope and good stereo correlation quality. Specifically, the elevation measurement error is expressed as follows: In the formula: and These are empirical coefficients, which respectively reflect the degree of influence of slope and related mass on elevation error; Step 5.5: Perform two consecutive weighted least squares fitting operations to check the elevation measurement error. Perform adjustment processing, calculate robust linear elevation change rate, and remove outliers outside the 99% confidence interval of the first fit; Step 5.6: Calculate the maximum allowable linear elevation change rate for each pixel. Centered on, with radius as Within a circular neighborhood, collect the linear elevation change rate of all pixels. Calculate its 80th percentile. With 20th percentile Elevation threshold , is represented as: The elevation difference threshold represents the reasonable fluctuation range of the elevation change rate within the neighborhood. The maximum linear change rate is used to constrain the elevation difference threshold. , is represented as: In the formula: Indicates the dynamic elevation difference threshold; Indicates the basic elevation difference threshold; This represents the time difference between the observation time and the reference DEM acquisition time. Step 5.7: Set the dynamic threshold. As a filtering condition applied in step 5.3, it is used for replacement. This filters out the remaining outliers and yields high-precision ASTER DEM data.

[0015] Furthermore, in step 5.5, two consecutive weighted least squares fitting operations are performed to correct the elevation measurement error. The specific process of adjustment is as follows: Step 5.5.1, First Fit: For pixels Using elevation observations at all its time points Perform linear fitting. Represented as: In the formula: Indicates time Elevation observation values; Indicates the intercept; Indicates the linear rate of change of elevation; Represents the residual; Elevation measurement error is used in the fitting process. Calculate weights After the fitting is complete, outliers outside the 99% confidence interval of the first fitting are deleted, that is, observations that meet the following conditions are deleted: In the formula: It represents the standard error of the fit; 2.576 corresponds to the critical value at the 99% confidence level. Step 5.5.2, Second Fitting: Use the filtered data to perform weighted least squares fitting again to obtain a robust linear rate of change of elevation, and retain the value after the second fitting.

[0016] Furthermore, the specific process of step 6 is as follows: Step 6.1: Construct a time covariance function that comprehensively considers multiple characteristics of glacier elevation changes. This function is composed of the following four kernel functions: 1) The linear kernel that captures long-term linear trends is represented as: In the formula: It is the variance parameter of the linear kernel; It's a time difference; 2) The sine square kernel of the periodic index that captures seasonal variations is expressed as: In the formula: It is the variance parameter of the periodic kernel; It is a cycle; It is a length scale parameter; 3) The radial basis function kernel that reflects the smoothness of elevation over time is expressed as: In the formula: It is the variance parameter of the RBF kernel. It is a length scale parameter; 4) The product of a rational quadratic kernel and a linear kernel that captures long-term nonlinear variations is expressed as: In the formula: It is the variance parameter of a rational quadratic kernel; These are shape parameters; It is a length scale parameter; Therefore, the final time covariance function is expressed as: Step 6.2: Based on the time covariance function Using the Gaussian process regression method with a time step of one month, the high-precision ASTER DEM data was interpolated into an elevation time series. Step 6.3: Execute Step 6.2 five times consecutively. After each Gaussian process regression, remove outlier observations based on the prediction error, gradually tighten the outlier judgment criteria, and generate time-continuous ASTER DEM data.

[0017] Furthermore, the data processing procedure of the Gaussian process regression method in step 6.2 is as follows: assuming that the elevation observations at any finite number of time points follow a multivariate Gaussian distribution, expressed as: In the formula: It is a vector of observed elevation values. It is the vector of elevation values ​​to be predicted; It is a mean vector; It is determined by the kernel function The constructed covariance matrix is ​​expressed as: In the formula: It is the covariance matrix between observation time points; It is the covariance matrix between the observed time point and the predicted time point; It is the covariance matrix between the predicted time point and the observed time point; It is the covariance matrix between the predicted time points; It is the observation noise variance, i.e., the elevation measurement error calculated in step 5.4. ; It is the identity matrix; Therefore, the elevation value at the predicted time point is obtained through a conditional probability distribution, expressed as: In the formula: .

[0018] Step 6.3: Execute Step 6.2 five times consecutively. After each Gaussian process regression, eliminate abnormal observations based on the prediction error and gradually tighten the outlier judgment criteria. Specifically, perform Gaussian process regression calculations with a monthly time step, that is, make predictions on the first day of each month to generate a spatial resolution time series of elevations with a resolution of 100 meters within a 10-kilometer area around the glacier, and generate time-continuous ASTER DEM data.

[0019] Furthermore, the specific process of step 7 is as follows: Step 7.1: Based on the elevation time series, calculate the elevation change of each pixel between adjacent time steps. In time and Elevation changes between , is represented as: In the formula: Represents pixels In time The elevation values ​​are compared at 12-month intervals to calculate the annual elevation change, thereby eliminating the impact of seasonal fluctuations. Step 7.2: Combine pixel area to calculate the volume change of each glacier. Its volume change is expressed as: In the formula, Represents glaciers Volume change, Represents pixels The area; Step 7.3: Aggregate the volume changes of all glaciers to the entire region, calculate the total volume change of the region, and express it as: Step 7.4: Based on the ice and snow density conversion coefficient, convert the volume change into a mass change, expressed as: In the formula: Indicates a change in quality; This represents the ice and snow density conversion coefficient, which is used in real time. As the average density of ice and snow; Therefore, the mass balance rate is expressed as: In the formula: Indicates the mass balance rate; This indicates the total area of ​​the glacier.

[0020] The present invention has the following beneficial technical effects: 1) This invention addresses the attitude jitter problem during ASTER satellite operation by establishing a complete error detection and correction process, which includes: performing two-dimensional correlation analysis around the theoretical polar line to directly observe the influence of the satellite's vertical orbit jitter component; using a 7th-order polynomial fitting to correct low-frequency errors; and superimposing eight sine curves to correct high-frequency jitter, thus eliminating the system deviation caused by jitter in a hierarchical and multi-scale manner.

[0021] 2) This invention establishes a multi-level deviation correction process based on reference ASTER DEM data, systematically eliminating various error sources: using quadratic regression fitting to correct systematic errors related to sensor angles, eliminating the deviation trend in the vertical orbit direction of the satellite; using polynomial fitting to correct low-frequency errors in the parallel orbit direction of the satellite, solving the elevation shift problem that varies with the satellite orbit position; and using sine / cosine function modeling to remove periodic high-frequency noise, effectively suppressing the ripple error caused by satellite jitter.

[0022] 3) In response to the characteristics of sparse and irregular distribution of satellite observation data over time, this invention innovatively adopts a Gaussian process regression method with combined kernel functions: by integrating linear kernels, periodic kernels, radial basis function kernels and rational quadratic kernels, a time covariance function is constructed that can simultaneously capture long-term linear trends, seasonal fluctuations, short-term smooth changes and nonlinear evolution, thereby reflecting multiple characteristics of glacier elevation changes, improving data coverage and accuracy, and helping to improve the accuracy of glacier mass balance calculation.

[0023] 4) This invention makes full use of multi-source auxiliary data such as reference ASTER DEM and glacier catalog data to establish a data fusion and complementarity mechanism: reference ASTER DEM is used for residual deviation correction and outlier filtering, providing constraints for stable terrain; glacier catalog data is used to generate processing buffers and statistical analysis ranges; for glacier areas without data coverage, the regional average elevation change rate is used for reasonable compensation. Attached Figure Description

[0024] Figure 1 This is a flowchart of the present invention; Figure 2 This is a comparison image of the original image (a) and the image after destriping (b); Figure 3 This is a schematic diagram illustrating the effect of jitter on the direction of the perspective rays in the rearview image; Figure 4 This is a schematic diagram of the optimal height approximation; Figure 5 This is a DEM data distribution map; Figure 6 This is a schematic diagram of the changes in glacier elevation in southeastern Tibet from 2015 to 2024, obtained by this invention. Detailed Implementation

[0025] The present invention will be further described in detail below with reference to specific embodiments. These descriptions are for explanation purposes only and are not intended to limit the scope of the invention.

[0026] refer to Figure 1As shown, a method for calculating glacier mass balance based on ASTER 3D data to generate elevation time series includes the following steps: Step 1: Collect ASTER L1A stereo images within the glacier study area. Select images with cloud cover less than 95% and a time span covering the study period as the original ASTER stereo image pairs for generating ASTER DEM data. The specific process is as follows: Step 1.1: Based on the geographical and temporal ranges, retrieve all available ASTERL1A stereo image pairs within the glacier study area; Step 1.2: Perform preliminary quality screening on the available ASTER L1A stereo image pairs. Specifically, select ASTER stereo image pairs with cloud cover of less than 95% and a time span covering the study period. The image pairs should have complete data files and no missing or damaged strips. In order to improve the uniformity of the time series, select image pairs from different years and seasons to capture the complete cycle of glacier changes. Step 1.3: Stitch together the ASTER stereo image pairs from every three consecutive days of acquisition to obtain the original ASTER stereo image pairs used to generate ASTERDEM data. Stitching facilitates the subsequent joint calculation of a longer DEM strip, increases the number of stable and ice-free terrains within a single DEM strip, and thus improves the accuracy of subsequent deviation correction. Step 2: Perform radiometric correction and destriping on the original ASTER stereo image pairs: Convert the luminance values ​​(DN) to surface reflectance, and use sensor response calibration data to correct for striping effects, obtaining high-quality ASTER stereo image pairs. The specific process is as follows: Step 2.1: Convert the DN values ​​of the original ASTER stereo image pairs into sensor radiance using radiometric calibration coefficients. , is represented as: In the formula: It is the radiance of the sensor; It is a numerical value; It is the gain factor; It is the offset; Step 2.2: Measure the sensor radiance. Converted to surface reflectance ， Represented as: In the formula: It is the surface reflectance; It is the distance between the Earth and the Sun; It is solar irradiance; It is the zenith angle of the sun; Step 2.3: Correct the striping effect using sensor response calibration data. The calibration function is a linear function composed of gain coefficient, offset coefficient and quadratic term coefficient, which can effectively reduce the striping effect in the image and obtain high-quality ASTER stereo image pairs. Step 3: Calculate the rational polynomial coefficient RPC model for the high-quality ASTER stereo image pair, and perform two-dimensional correlation analysis based on the theoretical epipolar line to detect and correct the jitter error in the vertical orbit direction of the satellite, obtaining the corrected ASTER stereo image pair and RPC model. The specific process is as follows: Step 3.1: Based on the high-quality ASTER stereo image pair, convert the geocentric coordinates to topographic Cartesian coordinates, and then to geodetic coordinates. Set the vertical grid density to 200 meters and define the range between -500 meters and +8850 meters. Step 3.2: Perform two-dimensional correlation analysis around the theoretical epipolar line to analyze the corresponding points of the ASTER stereo image pairs. The vibration intensity in pixels is equal to the distance between the maximum correlation point and the epipolar line, and the influence of the satellite vertical orbit jitter component is directly observed. Step 3.3: Model the parallax error caused by the jitter component in the vertical orbit direction of the satellite, and use a 7th-order polynomial to correct the low-frequency error. , is represented as: In the formula: Represents the fitting coefficients of a 7th-order polynomial. Represents the coordinates of the satellite's parallel orbit direction on the image; Step 3.4: Estimate the high-frequency jitter error by fitting the sum of eight sine curves along the axis on a 1000-pixel wide column with 90% overlap. , is represented as: In the formula: Representing the The amplitude of a sine curve represents the intensity of the jitter component at that frequency; Representing the The wavelength of a sine curve represents the spatial period of the frequency jitter component. Representing the The phase of a sine curve represents the spatial starting position of the jitter component at that frequency. The coordinates representing the direction of the satellite's parallel orbit in the image; Step 3.5: Calculate the results... and As correction value ，Applied to rear-view images, bilinear interpolation is used to resample the rear-view images to obtain corrected ASTER stereo image pairs and RPC models; Step 4: Calculate the initial ASTER DEM data using the corrected ASTER stereo image pairs and the RPC model. Correct the residual bias of the initial ASTER DEM data using reference ASTER DEM data. Through polynomial fitting and sine / cosine function modeling, systematic errors are gradually eliminated to obtain the corrected ASTER DEM data. The specific process is as follows: Step 4.1: Calculate the initial ASTER DEM data using the corrected ASTER stereo image pairs and the RPC model; Step 4.2: Calculate the elevation deviation between the generated initial ASTER DEM data and the reference ASTER DEM data in non-glacier areas. , is represented as: In the formula: This represents the initial ASTER DEM data; This indicates that the reference data is ASTER DEM data; Step 4.3: Use quadratic regression to fit the elevation deviation. The relationship with the sensor angle is used to correct errors in the satellite's vertical orbit system. , is represented as: In the formula: The sensor angle indicating the vertical orbital direction of the satellite; , and These are the coefficients of the quadratic regression model; Step 4.4: Fitting the low-frequency error of the satellite's parallel orbit direction , is represented as: In the formula: The coordinates of the satellite's parallel orbit direction on the image. ， 、 、 These are the coefficients of the quadratic regression model; Step 4.5: Use sine / cosine function fitting to remove high-frequency errors. , is represented as: In the formula: Indicates the first The wavelength of each frequency component; Indicates the first The amplitude coefficient of the sinusoidal term of each frequency component; It is the first The amplitude coefficient of the cosine term of each frequency component; It is an index of the frequency component; These are the coordinates of the satellite's parallel orbit direction on the image; Step 4.6, use , and The initial ASTER DEM data is gradually corrected to obtain the corrected ASTER DEM data. Step 5: Generate a buffer using glacier contour data and crop and correct the ASTER DEM data. Filter outliers using reference ASTER DEM data. Remove outliers through weighted least squares fitting and Gaussian process regression to obtain high-precision ASTER DEM data. The specific process is as follows: Step 5.1: Use glacier contour data to generate a 10 km buffer zone around each glacier, and then crop and correct the ASTER DEM data based on the buffer zone. Step 5.2: Exclude corrected ASTER DEM data with a root mean square error greater than 20 meters from the reference ASTER DEM data on ice-free terrain, and retain corrected ASTER DEM data that meet the requirements. Step 5.3, with radius The circle represents the filtering window, excluding data where the absolute height difference between the corrected ASTER DEM data and the reference ASTER DEM data exceeds a threshold. The pixels are represented as: In the formula: The threshold representing outlier filtering; Representing the In this iteration, the corrected ASTER DEM data is in pixels. The elevation value of the location; Represents reference ASTER DEM data at the pixel level The elevation value of the location; Step 5.4: Considering terrain slope and stereo correlation quality, calculate the elevation measurement error of each pixel selected and retained in Step 5.3. , is represented as: In the formula: Indicates elevation measurement error; Indicates the slope of the terrain; Indicates the stereo correlation quality; This represents the registration error. Ideally, the registration error is calculated on pixels with a small slope and good stereo correlation quality. Specifically, the elevation measurement error is expressed as follows: In the formula: and These are empirical coefficients, which respectively reflect the degree of influence of slope and related mass on elevation error; Step 5.5: Perform two consecutive weighted least squares fitting operations to check the elevation measurement error. Adjustment processing is performed, robust linear elevation change rates are calculated, and outliers outside the 99% confidence interval of the first fit are removed. The specific process is as follows: Step 5.5.1, First Fit: For pixels Using elevation observations at all its time points Perform linear fitting. Represented as: In the formula: Indicates time Elevation observation values; Indicates the intercept; Indicates the linear rate of change of elevation; Represents the residual; Elevation measurement error is used in the fitting process. Calculate weights After the fitting is complete, outliers outside the 99% confidence interval of the first fitting are deleted, that is, observations that meet the following conditions are deleted: In the formula: It represents the standard error of the fit; 2.576 corresponds to the critical value at the 99% confidence level. Step 5.5.2, Second Fitting: Use the filtered data to perform weighted least squares fitting again to obtain a robust linear rate of change of elevation, and retain the value after the second fitting; Step 5.6: Calculate the maximum allowable linear elevation change rate for each pixel. Centered on, with radius as Within a circular neighborhood, collect the linear elevation change rate of all pixels. Calculate its 80th percentile. With 20th percentile Elevation threshold , is represented as: The elevation difference threshold represents the reasonable fluctuation range of the elevation change rate within the neighborhood. The maximum linear change rate is used to constrain the elevation difference threshold. , is represented as: In the formula: Indicates the dynamic elevation difference threshold; Indicates the basic elevation difference threshold; This represents the time difference between the observation time and the acquisition time of the reference ASTER DEM data; Step 5.7: Set the dynamic threshold. As a filtering condition applied in step 5.3, it is used for replacement. This filters out the remaining outliers and obtains high-precision ASTER DEM data. Step 6: Using the Gaussian process regression method, the high-precision ASTER DEM data is interpolated into an elevation time series with a one-month time step to generate time-continuous ASTER DEM data. The specific process is as follows: Step 6.1: Construct a time covariance function that comprehensively considers multiple characteristics of glacier elevation changes. This function is composed of the following four kernel functions: 1) The linear kernel that captures long-term linear trends is represented as: In the formula: It is the variance parameter of the linear kernel; It's a time difference; 2) The sine square kernel of the periodic index that captures seasonal variations is expressed as: In the formula: It is the variance parameter of the periodic kernel; It is a cycle; It is a length scale parameter; 3) The radial basis function kernel that reflects the smoothness of elevation over time is expressed as: In the formula: It is the variance parameter of the RBF kernel. It is a length scale parameter; 4) The product of a rational quadratic kernel and a linear kernel that captures long-term nonlinear variations is expressed as: In the formula: It is the variance parameter of a rational quadratic kernel; These are shape parameters; It is a length scale parameter; Therefore, the final time covariance function is expressed as: The time covariance function can simultaneously describe the linear trend, seasonal fluctuations, short-term smooth changes, and nonlinear long-term changes in glacier elevation. Step 6.2: Based on the time covariance function, use the Gaussian process regression method with a time step of one month to interpolate the high-precision ASTER DEM data; The data processing procedure of the Gaussian process regression method is as follows: It is assumed that the elevation observations at any finite number of time points follow a multivariate Gaussian distribution, expressed as: In the formula: It is a vector of observed elevation values. It is the vector of elevation values ​​to be predicted; It is a mean vector; It is determined by the kernel function The constructed covariance matrix is ​​expressed as: In the formula: It is the covariance matrix between observation time points; It is the covariance matrix between the observed time point and the predicted time point; It is the covariance matrix between the predicted time point and the observed time point; It is the covariance matrix between the predicted time points; It is the observation noise variance, i.e., the elevation measurement error calculated in step 5.4. ; It is the identity matrix; Therefore, the elevation value at the predicted time point is obtained through a conditional probability distribution, expressed as: In the formula: Step 6.3: To further filter out outliers, Step 6.2 is executed five times consecutively. After each Gaussian process regression, outlier observations are removed based on the prediction error. The outlier judgment criteria are gradually tightened so that the final retained observation data is highly consistent with the Gaussian process regression model. Specifically, Gaussian process regression calculation is performed with a monthly time step, that is, prediction is performed on the first day of each month to generate an elevation time series with a spatial resolution of 100 meters within a 10-kilometer area around the glacier, generating time-continuous ASTER DEM data. Step 7: Calculate the volume change of each glacier based on the elevation time series, and then aggregate the volume changes of each glacier to obtain the overall mass balance change of the regional glaciers. The specific process is as follows: Step 7.1: Based on the elevation time series, calculate the elevation change of each pixel between adjacent time steps. In time and Elevation changes between , is represented as: In the formula: Represents pixels In time The elevation values ​​are compared at 12-month intervals to calculate the annual elevation change, thereby eliminating the impact of seasonal fluctuations. Step 7.2: Combine pixel area to calculate the volume change of each glacier. Its volume change is expressed as: In the formula, Represents glaciers Volume change, Represents pixels The area; Step 7.3: Aggregate the volume changes of all glaciers to the entire region, calculate the total volume change of the region, and express it as: Step 7.4: Based on the ice and snow density conversion coefficient, convert the volume change into a mass change, expressed as: In the formula: Indicates a change in quality; This represents the ice and snow density conversion coefficient, which is used in real time. As the average density of ice and snow; Therefore, the mass balance rate is expressed as: In the formula: Indicates the mass balance rate; This indicates the total area of ​​the glacier.

[0027] Figure 2 (a) is the original image. Figure 2 (b) shows the image after destriping, from which it can be seen that: Figure 2 The original image in (a) shows obvious striping, which is caused by the inconsistent sensitivity of each unit in the sensor linear array; Figure 2 (b) The striping effect in the image after destripping of sensor response calibration data was significantly suppressed, and the image quality was significantly improved, providing a better data foundation for subsequent stereo matching.

[0028] Figure 3 This diagram illustrates how satellite attitude jitter affects the geometry of ASTER stereo image pairs. The figure shows the imaging geometry of the 3N band (nadir view) and the 3B band (back view), as well as the positional relationships of terrain undulations, candidate altitudes, and optimal altitudes. Blue dots represent the image point positions corresponding to candidate altitudes, and green dots represent the true matching points corresponding to optimal altitudes. As can be seen from the figure, satellite vertical orbit jitter can cause the perspective rays of the two bands to not intersect correctly. Therefore, correction is necessary to obtain accurate elevation information.

[0029] Depend on Figure 4 It can be seen that the DEM calculation using the rational polynomial coefficient model is an iterative process of determining the optimal elevation of ground points by gradually narrowing the height search range. For each point on the target grid, the algorithm calculates the normalized cross-correlation coefficient of two images at different candidate heights and selects the height with the highest correlation as the optimal elevation value of that point.

[0030] Figure 5 The figure shows the spatial distribution of 1,632 DEM bands generated in southeastern Tibet. Each black band represents the coverage area of ​​a DEM data, and the red dashed lines mark the glacier distribution range in the study area. The inset in the upper right corner shows the geographical location of the study area on the Tibetan Plateau. As can be seen from the figure, the DEM bands cover the entire glacier region of southeastern Tibet, and there is overlap between the bands. This high-density data coverage provides sufficient observation data for subsequent elevation time series interpolation. In addition, the scale and geographic coordinate range (91°E-100°E) are also marked in the figure to facilitate understanding the spatial scale of the data.

[0031] Figure 6 The time-series elevation interpolation results for the first day of each year from 2015 to 2024 are presented in grid form. Each square in the figure represents the elevation distribution of a specific glacier region at a specific time point. The color intensity indicates the elevation change, with the color scale ranging from -1.5 m / a to +1.0 m / a. This figure effectively demonstrates the temporally continuous and spatially complete elevation change data generated by the Gaussian process interpolation method, providing intuitive visualization results for glacier mass balance analysis.

Claims

1. A method for calculating glacier mass balance based on ASTER 3D data to generate elevation time series, characterized in that, Includes the following steps: Step 1: Collect ASTER L1A stereo images within the glacier study area. Select images with cloud cover of less than 95% and a time span covering the study period as the original ASTER stereo image pairs for generating ASTER DEM data. Step 2: Perform radiometric correction and destriping on the original ASTER stereo image pairs to obtain high-quality ASTER stereo image pairs. Step 3: Calculate the rational polynomial coefficient RPC model of the high-quality ASTER stereo image pair, and perform two-dimensional correlation analysis based on the theoretical epipolar line. Use a 7th-order polynomial to fit and correct the low-frequency error, and fit the sum of 8 sine curves to correct the high-frequency jitter error. Correct the jitter error in the vertical orbit direction of the satellite through the low-frequency error and high-frequency error to obtain the corrected ASTER stereo image pair and RPC model. Step 4: Calculate the initial ASTER DEM data using the corrected stereo image pair and RPC model. Use the reference ASTER DEM data to perform residual bias correction on the initial ASTER DEM data. Then, use quadratic regression to fit the relationship between elevation deviation and sensor angle to correct the satellite vertical orbit system error. Finally, use sine / cosine function modeling to gradually eliminate system errors and high-frequency noise to obtain the corrected ASTER DEM data. Step 5: Use glacier contour data to generate a buffer and crop the corrected ASTER DEM data. Use reference ASTER DEM data to filter out outliers. Remove outliers through weighted least squares fitting and Gaussian process regression to obtain high-precision ASTER DEM data. Step 6: Construct a time covariance function using a linear kernel, an exponential sine square kernel, a radial basis function kernel, and a rational quadratic kernel. Then, use the Gaussian process regression method to interpolate the high-precision ASTER DEM data into an elevation time series with a one-month time step to generate time-continuous ASTER DEM data. Step 7: Calculate the volume change of each glacier based on the elevation time series, and then aggregate the volume changes of each glacier to obtain the overall mass balance change of the regional glaciers.

2. The method for calculating glacier mass balance based on ASTER 3D data to generate elevation time series according to claim 1, characterized in that, The specific process of step 1 is as follows: Step 1.1: Based on the geographical and temporal ranges, retrieve all available ASTER L1A stereo image pairs within the glacier study area; Step 1.2: Perform preliminary quality screening on the available ASTER L1A stereo image pairs: select ASTER stereo image pairs with cloud cover of less than 95% and a time span covering the study period. Step 1.3: Stitch together the ASTER stereo image pairs from every three consecutive acquisition days to obtain the original ASTER stereo image pairs used to generate ASTER DEM data.

3. The method for calculating glacier mass balance based on ASTER 3D data to generate elevation time series according to claim 1, characterized in that, The specific process of step 2 is as follows: Step 2.1: Convert the DN values ​​of the original ASTER stereo image pair into sensor radiance using radiometric calibration coefficients. , is represented as: In the formula: It is the radiance of the sensor; It is a numerical value; It is the gain factor; It is the offset; Step 2.2: Measure the sensor radiance. Converted to surface reflectance ， Represented as: In the formula: It is the surface reflectance; It is the distance between the Earth and the Sun; It is solar irradiance; It is the zenith angle of the sun; Step 2.3: Correct the striping effect using sensor response calibration data. The calibration function is a linear function composed of gain coefficient, offset coefficient and quadratic term coefficient, used to reduce the striping effect in the image and obtain high-quality ASTER stereo image pairs.

4. The method for calculating glacier mass balance based on ASTER 3D data to generate elevation time series according to claim 1, characterized in that, The specific process of step 3 is as follows: Step 3.1: Based on the high-quality ASTER stereo image pair, convert the geocentric coordinates to topographic Cartesian coordinates, and then to geodetic coordinates. Set the vertical grid density to 200 meters and define the range between -500 meters and +8850 meters. Step 3.2: Perform two-dimensional correlation analysis around the theoretical epipolar line to analyze the corresponding points of the ASTER stereo image pairs. The vibration intensity in pixels is equal to the distance between the maximum correlation point and the epipolar line, and the influence of the satellite vertical orbit jitter component is directly observed. Step 3.3: Model the parallax error caused by the jitter component in the vertical orbit direction of the satellite, and use a 7th-order polynomial to correct the low-frequency error. , is represented as: In the formula: Represents the fitting coefficients of a 7th-order polynomial. Represents the coordinates of the satellite's parallel orbit direction on the image; Step 3.4: Estimate the high-frequency jitter error by fitting the sum of eight sine curves along the axis on a 1000-pixel wide column with 90% overlap. , is represented as: In the formula: Representing the The amplitude of a sine curve represents the intensity of the jitter component at that frequency; Representing the The wavelength of a sine curve represents the spatial period of the frequency jitter component. Representing the The phase of a sine curve represents the spatial starting position of the jitter component at that frequency. The coordinates representing the direction of the satellite's parallel orbit in the image; Step 3.5: Calculate the results... and As correction value ， Applied to rear-view images, bilinear interpolation is used to resample the rear-view images to obtain corrected ASTER stereo image pairs and RPC models.

5. The method for calculating glacier mass balance based on ASTER stereo data and time series elevation generation according to claim 1, characterized in that, The specific process of step 4 is as follows: Step 4.1: Calculate the initial ASTER DEM data using the corrected ASTER stereo image pairs and the RPC model; Step 4.2: Calculate the elevation deviation between the generated initial ASTER DEM data and the reference ASTER DEM data in non-glacier areas. , is represented as: In the formula: This represents the initial ASTER DEM data; This indicates that the reference data is ASTER DEM data; Step 4.3: Use quadratic regression to fit the elevation deviation. The relationship with the sensor angle is used to correct errors in the satellite's vertical orbit system. , is represented as: In the formula: The sensor angle indicating the vertical orbital direction of the satellite; 、 and These are the coefficients of the quadratic regression model; Step 4.4: Fitting the low-frequency error of the satellite's parallel orbit direction , is represented as: In the formula: The coordinates of the satellite's parallel orbit direction on the image. ， 、 、 These are the coefficients of the quadratic regression model; Step 4.5: Use sine / cosine function fitting to remove high-frequency errors. , is represented as: In the formula: Indicates the first The wavelength of each frequency component; surface Show the first The amplitude coefficient of the sinusoidal term of each frequency component; It is the first The amplitude coefficient of the cosine term of each frequency component; It is an index of the frequency component; These are the coordinates of the satellite's parallel orbit direction on the image; Step 4.6, use , and The initial ASTER DEM data was progressively corrected to obtain the corrected ASTER DEM data.

6. The method for calculating glacier mass balance based on ASTER stereo data and time series elevation generation according to claim 1, characterized in that, The specific process of step 5 is as follows: Step 5.1: Use glacier contour data to generate a 10 km buffer zone around each glacier, and then crop and correct the ASTER DEM data based on the buffer zone. Step 5.2: Exclude corrected ASTER DEM data with a root mean square error greater than 20 meters from the reference ASTER DEM data on ice-free terrain, and retain corrected ASTER DEM data that meet the requirements. Step 5.3, with radius The circle represents the filtering window, excluding data where the absolute height difference between the corrected ASTER DEM data and the reference ASTER DEM data exceeds a threshold. The pixels are represented as: In the formula: The threshold representing outlier filtering; Representing the In this iteration, the corrected ASTER DEM data is in pixels. The elevation value of the location; Represents reference ASTER DEM data at the pixel level The elevation value of the location; Step 5.4: Considering terrain slope and stereo correlation quality, calculate the elevation measurement error of each pixel selected and retained in Step 5.

3. , is represented as: In the formula: Indicates elevation measurement error; Indicates the slope of the terrain; Indicates the stereo correlation quality; This represents the registration error. Ideally, the registration error is calculated on pixels with a small slope and good stereo correlation quality. Specifically, the elevation measurement error is expressed as follows: In the formula: and These are empirical coefficients, which respectively reflect the degree of influence of slope and related mass on elevation error; Step 5.5: Perform two consecutive weighted least squares fitting operations to check the elevation measurement error. Perform adjustment processing, calculate robust linear elevation change rate, and remove outliers outside the 99% confidence interval of the first fit; Step 5.6: Calculate the maximum allowable linear elevation change rate for each pixel. Centered on, with radius as Within a circular neighborhood, collect the linear elevation change rate of all pixels. Calculate its 80th percentile. With 20th percentile Elevation threshold , is represented as: The elevation difference threshold represents the reasonable fluctuation range of the elevation change rate within the neighborhood. The maximum linear change rate is used to constrain the elevation difference threshold. , is represented as: In the formula: Indicates the dynamic elevation difference threshold; Indicates the basic elevation difference threshold; This represents the time difference between the observation time and the reference DEM acquisition time. Step 5.7: Set the dynamic threshold. As a filtering condition applied in step 5.3, it is used for replacement. This filters out the remaining outliers and yields high-precision ASTER DEM data.

7. The method for calculating glacier mass balance based on ASTER stereo data and time series elevation generation according to claim 6, characterized in that, In step 5.5, two consecutive weighted least squares fitting operations are performed to correct the elevation measurement error. The specific process of adjustment is as follows: Step 5.5.1, First Fit: For pixels Using elevation observations at all its time points Perform linear fitting. Represented as: In the formula: Indicates time Elevation observation values; Indicates the intercept; Indicates the linear rate of change of elevation; Represents the residual; Elevation measurement error is used in the fitting process. Calculate weights After the fitting is complete, outliers outside the 99% confidence interval of the first fitting are deleted, that is, observations that meet the following conditions are deleted: In the formula: It represents the standard error of the fit; 2.576 corresponds to the critical value at the 99% confidence level. Step 5.5.2, Second Fitting: Use the filtered data to perform weighted least squares fitting again to obtain a robust linear rate of change of elevation, and retain the value after the second fitting.

8. The method for calculating glacier mass balance based on ASTER stereo data and time series elevation generation according to claim 1, characterized in that, The specific process of step 6 is as follows: Step 6.1: Construct a time covariance function that comprehensively considers multiple characteristics of glacier elevation changes. This function is composed of the following four kernel functions: 1) The linear kernel that captures long-term linear trends is represented as: In the formula: It is the variance parameter of the linear kernel; It's a time difference; 2) The sine square kernel of the periodic index that captures seasonal variations is expressed as: In the formula: It is the variance parameter of the periodic kernel; It is a cycle; It is a length scale parameter; 3) The radial basis function kernel that reflects the smoothness of elevation over time is expressed as: In the formula: It is the variance parameter of the RBF kernel. It is a length scale parameter; 4) The product of a rational quadratic kernel and a linear kernel that captures long-term nonlinear variations is expressed as: In the formula: It is the variance parameter of a rational quadratic kernel; These are shape parameters; It is a length scale parameter; Therefore, the final time covariance function is expressed as: Step 6.2: Based on the time covariance function Using the Gaussian process regression method with a time step of one month, the high-precision ASTER DEM data was interpolated into an elevation time series. Step 6.3: Execute Step 6.2 five times consecutively. After each Gaussian process regression, remove outlier observations based on the prediction error, gradually tighten the outlier judgment criteria, and generate time-continuous ASTER DEM data.

9. The method for calculating glacier mass balance based on ASTER stereo data and time series elevation generation according to claim 8, characterized in that, The data processing procedure of the Gaussian process regression method in step 6.2 is as follows: Assuming that the elevation observations at any finite number of time points follow a multivariate Gaussian distribution, it can be expressed as: In the formula: It is a vector of observed elevation values. It is the vector of elevation values ​​to be predicted; It is a mean vector; It is determined by the kernel function The constructed covariance matrix is ​​expressed as: In the formula: It is the covariance matrix between observation time points; It is the covariance matrix between the observed time point and the predicted time point; It is the covariance matrix between the predicted time point and the observed time point; It is the covariance matrix between the predicted time points; It is the observation noise variance, i.e., the elevation measurement error calculated in step 5.

4. ; It is the identity matrix; Therefore, the elevation value at the predicted time point is obtained through a conditional probability distribution, expressed as: In the formula: Step 6.3: Execute Step 6.2 five times consecutively. After each Gaussian process regression, eliminate abnormal observations based on the prediction error and gradually tighten the outlier judgment criteria. Specifically, perform Gaussian process regression calculations with a monthly time step, that is, make predictions on the first day of each month to generate a spatial resolution time series of elevations with a resolution of 100 meters within a 10-kilometer area around the glacier, and generate time-continuous ASTER DEM data.

10. The method for calculating glacier mass balance based on ASTER stereo data and time series elevation generation according to claim 1, characterized in that, The specific process of step 7 is as follows: Step 7.1: Based on the elevation time series, calculate the elevation change of each pixel between adjacent time steps. In time and Elevation changes between , is represented as: Mode middle: Represents pixels In time The elevation values ​​are compared at 12-month intervals to calculate the annual elevation change, thereby eliminating the impact of seasonal fluctuations. Step 7.2: Combine pixel area to calculate the volume change of each glacier. ,That Volume change is expressed as: Mode middle, Represents glaciers Volume change, Represents pixels The area; Step 7.3: Aggregate the volume changes of all glaciers to the entire region, calculate the total volume change of the region, and express it as: Step 7.4: Based on the ice and snow density conversion coefficient, convert the volume change into a mass change, expressed as: In the formula ： Indicates a change in quality; This represents the ice and snow density conversion coefficient, which is used in real time. As the average density of ice and snow; Therefore, the mass balance rate is expressed as: In the formula: Indicates mass balance rate ； This indicates the total area of ​​the glacier.

Images