A distributed glacier degree-day factor calculation method based on multi-source remote sensing data

Abstract

The application discloses a kind of distributed glacier degree-day factor calculation method based on multi-source remote sensing data, specifically including the following steps: step 1, multi-source data acquisition and pretreatment;Step 2, glacier under surface classification delimitation;Step 3, grid scale glacier ablation amount estimation;Step 4, spatial continuous positive accumulated temperature sequence construction;Step 5: distributed initial degree-day factor inversion;Step 6: spatial calibration of distributed glacier degree-day factor.The method disclosed in the application realizes the distributed calculation of glacier degree-day factor, solves the representative error caused by the traditional method "point instead of surface", and can finely depict the spatial difference of glacier internal ablation.The application fully considers the influence of moraine coverage and non-moraine coverage on ablation, improves the physical accuracy of glacier runoff simulation under complex under surface conditions, combines the high spatial coverage capability of multi-source remote sensing and the high-precision advantage of ground observation, and provides reference basis for glacier research in remote areas lacking observation data.

Description

Technical Field

[0001] This invention belongs to the field of glacier hydrology technology, and in particular relates to a distributed method for calculating the glacier day-degree factor based on multi-source remote sensing data. Background Technology

[0002] Against the backdrop of escalating global warming, glaciers, as a core component of the "Asian water tower," are rapidly melting, leading not only to sea-level rise but also profoundly impacting water resource security, ecological balance, and the frequency of natural disasters in high-altitude regions and their downstream basins. Therefore, accurate simulation and forecasting of glacier runoff changes are of significant strategic importance. Among numerous methods for simulating glacier ablation, the degree-day model has become the mainstream approach in glacier hydrology research and application due to its clear physical mechanisms, moderate data threshold, and high efficiency in large-scale distributed simulations. The degree-day factor (DDF), as the core parameter connecting positive accumulated temperature and ablation, directly determines the simulation accuracy and forecast reliability of the runoff model.

[0003] However, existing methods for determining the diurnal variation factor of glaciers have significant limitations both theoretically and practically, severely restricting the simulation accuracy of distributed models.

[0004] First, there is a severe lack of spatial representativeness. Traditional methods rely heavily on data from a very limited number of discontinuous ground observation piles. Researchers are often forced to adopt a "point-to-surface" strategy, using the DDF value of a single site to extrapolate the melting characteristics of the entire glacier, or even across regions, the melting characteristics of neighboring glaciers. This method completely ignores the complex spatial heterogeneity within glaciers—variations in altitude, aspect, slope, and local microclimate can all lead to orders-of-magnitude differences in snow and ice melt rates. This statistically significant "spatial distortion" means that models often fail to capture the true spatial pattern of glacier melt when calculating large-scale runoff.

[0005] Secondly, the coupling of underlying physical characteristics is insufficient. The glacier surface is not a single homogeneous structure, but a complex system composed of snow cover, non-moored areas (open ice areas), and moored areas. In particular, the ablation mechanisms of moored and non-moored areas are fundamentally different: moored thickness plays a drastically different regulatory role in the ablation of the underlying ice mass through thermal barrier or endothermic enhancement effects. However, most current models still tend to use uniform or simplified DDF parameters, failing to finely distinguish the thickness distribution of moored areas and their profound impact on albedo and thermal conductivity. This neglect of underlying surface heterogeneity leads to significant systematic biases in runoff simulation results in glacier regions with widespread moored development.

[0006] In summary, how to overcome the geographical limitations of traditional observations and obtain distributed time-of-day factors with high spatial resolution that can dynamically reflect the physical heterogeneity of the glacier's underlying surface (moorland and non-moorland) in the absence of high-altitude measured data has become a key problem that urgently needs to be solved in the current interdisciplinary field of glaciology and hydrology. Summary of the Invention

[0007] The purpose of this invention is to provide a distributed method for calculating the diurnal glacier factor based on multi-source remote sensing data, so as to solve the above-mentioned technical problems.

[0008] To achieve the above objectives, the present invention provides the following technical solution:

[0009] This invention discloses a distributed method for calculating the diurnal glacier factor based on multi-source remote sensing data, the method comprising the following steps:

[0010] Step 1: Multi-source data acquisition and preprocessing: Acquire multi-temporal remote sensing data of the study area, including: Landsat 8 / 9 OLI multispectral imagery, MODIS daily surface temperature products, and dual-period elevation data with time spans; acquire daily temperature data from synchronized surface meteorological stations; and preprocess the acquired data.

[0011] Step 2: Classification and demarcation of glacier underlying surface: Based on the acquired remote sensing data, the glacier area is divided into grids according to the set resolution, and the glacier grids are divided into three categories according to different underlying surfaces: snow cover area, non-moored area, and moored area, and the boundaries of different types of underlying surfaces are demarcated.

[0012] Step 3: Grid-scale Glacier Ablation Estimation: By distinguishing between snow-covered areas, non-moored areas, and moored areas, the elevation data of the corresponding grids is determined, thereby determining the elevation difference of each grid at different time periods. Finally, the remote sensing mass balance data for each grid, i.e., the glacier ablation amount, is generated. This specifically includes the following steps:

[0013] 1) Calculate elevation changes: Difference is performed between two periods of elevation data to obtain the elevation changes at the grid scale. :

[0014] (3)

[0015] In the formula, , The elevation data are for times t1 and t2, respectively. The selection of times t1 and t2 is determined according to the needs of the researchers.

[0016] 2) Dynamic thinning correction: To obtain the net ablation, the dynamic uplift or subsidence elevation caused by glacial movement must be deducted.

[0017] (4)

[0018] In the formula, This refers to changes in net elevation. This refers to the elevation change between t1 and t2; The rate of rise or fall; The time interval between t1 and t2;

[0019] 3) Estimation of ablation: Based on the glacier underlying surface classification results from step 2, different density conversion formulas are set to estimate the glacier ablation:

[0020] (5)

[0021] In the formula, This represents the amount of glacier ablation per unit grid within a specific time period. The density of ice; The density of water;

[0022] Step 4: Construction of Spatial Continuous Positive Accumulated Temperature Sequence: Matching the data with the divided glacier grid, interpolating to generate continuous meteorological driving data for each grid, and calculating the spatial continuous positive accumulated temperature sequence (PDD) based on this; the specific process includes:

[0023] First, temperature downscaling is performed: The MODIS LST data preprocessed in step 1 is downscaled to a resolution consistent with the glacier grid using geographic weighted regression or the vertical lapse rate method;

[0024] Then, positive accumulated temperature calculation is performed: the sum of the daily average temperatures above 0℃ for each grid within the ablation period (t1 to t2) is statistically analyzed to construct a spatially continuous positive accumulated temperature sequence.

[0025] (6)

[0026] In the formula, The positive accumulated temperature value for each grid where the daily average temperature is greater than 0℃; The average temperature of each grid cell on day d;

[0027] Step 5, Distributed Initial Degree-Day Factor Inversion: Based on the degree-day model principle, invert the initial degree-day factor for each grid; the inversion formula for the initial degree-day factor is:

[0028] (7)

[0029] In the formula: This is the initial degree-day factor, in mm·d⁻¹·℃⁻¹;

[0030] Step 6: Spatial calibration of distributed glacier diurnal factors: Input the daily temperature data of each glacier grid and the initial diurnal factor into the diurnal model to simulate and generate the mass balance data, i.e., ablation, for each grid. The specific formula is:

[0031] * (8)

[0032] In the formula, M represents the simulated ablation amount;

[0033] Based on this, the simulated mass balance data is compared with the remote sensing mass balance data to calibrate the initial degree-day factor and generate the final distributed degree-day factor data.

[0034] Furthermore, the time span in the dual-period elevation data with a time span described in step 1 is at least one year.

[0035] Furthermore, the preprocessing of the acquired data described in step 1 specifically includes: performing terrain and atmospheric correction on Landsat image data; performing 3D registration and deviation correction on elevation data from different sources to ensure spatial reference consistency; and performing cloud masking and missing value filling on MODIS daily surface temperature products.

[0036] Furthermore, the specific process for defining the boundaries of different types of underlying surfaces in step 2 is as follows:

[0037] First, the snow cover area is demarcated using the Normalized Diffusion Index (NDSI). The formula for calculating the NDSI is:

[0038] (1)

[0039] In the formula, and These represent the reflectance in the green light band and the shortwave infrared band, respectively.

[0040] When NDSI is greater than a certain set threshold, it indicates a snow-covered area; otherwise, it indicates a non-snow-covered area.

[0041] For the identification of till-covered areas, since the thermal inertia of till is significantly different from that of pure ice, the radiation temperature of the till surface is usually higher than that of pure ice during the ablation period; therefore, the classification discriminant is established as follows:

[0042] (2)

[0043] In the formula, For glacier grids; This is a snow-covered area; This is an area not covered by tidal flats; The area is covered by topsoil; The set NDSI threshold is 0.3 to 0.5. The temperature threshold is 0~0.5; This refers to the temperature of the surface of the tundish or ice.

[0044] Furthermore, the specific process of calibrating the initial degree-day factor in step 6 to generate the final distributed degree-day factor data is as follows: If based on the initial degree-day factor... The mass balance data (M) of each simulated grid and the remotely sensed mass balance data are... If the comparison error is within ±5%, the day factor is determined; otherwise, the day factor is further calibrated using the residuals of the simulated and remotely sensed mass balance data, and then compared with the remotely sensed mass balance data to continuously iterate and obtain a day factor that meets the error requirements. :

[0045] = (9)

[0046] Finally, distributed daily time factor data that distinguishes between snow-covered areas, non-moored areas, and moored areas are generated, thus obtaining a spatial distribution map of glacier daily time factors in the study area.

[0047] The beneficial effects of this invention are as follows: The method described in this invention enables distributed calculation of glacier diurnal factors, solving the representativeness error caused by the traditional method of "representing the whole by a point," and can accurately characterize the spatial differences in glacier ablation. This invention fully considers the influence of moraine cover and non-moraine cover on ablation, improves the physical accuracy of glacier runoff simulation under complex underlying surface conditions, and integrates the high spatial coverage capability of multi-source remote sensing with the high precision advantages of ground observation, providing a reliable means of parameter acquisition for glacier research in remote areas without monitoring stations or lacking observation data.

[0048] The present invention will now be described in further detail with reference to the accompanying drawings and specific embodiments. Attached Figure Description

[0049] Figure 1 This is a schematic diagram of the method flow described in this invention;

[0050] Figure 2 This is a schematic diagram comparing simulated mass balance data and remote sensing mass balance data in Example 1;

[0051] Figure 3 This is a spatial distribution map of the glacier diurnal factor generated in Example 1. Detailed Implementation

[0052] This invention discloses a distributed method for calculating the diurnal temperature variation factor of glaciers based on multi-source remote sensing data, such as... Figure 1As shown, the method includes the following steps:

[0053] Step 1: Acquisition and Preprocessing of Multi-Source Data: Acquire multi-temporal remote sensing data of the study area, including: Landsat 8 / 9 OLI multispectral imagery, MODIS daily surface temperature products (such as MOD11A1, with a spatial resolution of 1 km), and dual-period elevation data with a time span (at least one year) (such as SRTM and ICESat-2 altimetry data); acquire daily temperature data from synchronized ground meteorological stations.

[0054] The acquired data is preprocessed, including topographic and atmospheric correction of Landsat image data; 3D registration and deviation correction of elevation data from different sources to ensure spatial reference consistency; and cloud masking and missing value filling of MODIS daily surface temperature products.

[0055] Step 2: Glacier surface classification and demarcation: Based on the acquired remote sensing data (including Landsat image data and elevation data, etc.), the glacier area is divided into grids according to the set resolution (the resolution is generally 30m). Based on the multispectral features of the Landsat image (including vegetation, surface features, etc.), the glacier grids are divided into three categories according to different surface areas: snow cover area, non-moored area (open ice area), and moored area.

[0056] For delineating different types of underlying surfaces, the Normalized Difference Snow Index (NDSI) is first used to delineate the snow cover area. The formula for calculating the NDSI is:

[0057] (1)

[0058] In the formula, and These represent the reflectance in the green light band and the shortwave infrared band, respectively.

[0059] When NDSI is greater than a certain set threshold, it indicates a snow-covered area; otherwise, it indicates a non-snow-covered area.

[0060] For identifying till-covered areas, since the thermal inertia of till differs significantly from that of pure ice, the radiation temperature of the till surface is typically higher than that of pure ice during the ablation period. Therefore, a classification discriminant can be established as follows:

[0061] (2)

[0062] In the formula, For glacier grids; This is a snow-covered area; This is an area not covered by tidal flats; The area is covered by topsoil; The set NDSI threshold is typically 0.3 to 0.5 (which can be adjusted based on snow cover imagery and region). This is the temperature threshold, typically 0~0.5 (which can be adjusted based on snow cover images and the region). This refers to the temperature of the surface of the tundish or ice.

[0063] Step 3: Grid-scale Glacier Ablation Estimation: By distinguishing between snow-covered, non-moored, and moored glacier grids, the elevation and other data of the corresponding grids are determined. This allows for the determination of the elevation difference for each grid over different time periods, ultimately generating remote sensing mass balance data for each grid. Glacier mass balance data refers to the difference between the amount of snow and ice gained and lost by a glacier within a certain period (usually one year). The changes in elevation data of the glacier grid over two periods are considered as remote sensing mass balance data, also known as glacier ablation.

[0064] The estimation of glacier ablation at the grid scale includes the following steps:

[0065] 1) Calculate elevation change: Use the difference between two periods of elevation data to obtain the elevation change value at the grid scale (i.e., the elevation change of a certain location between two time periods). :

[0066] (3)

[0067] In the formula, , The elevation data are for times t1 and t2, respectively. The choice of time t1 and t2 can be determined according to the needs of the researchers, and the time interval should be at least one year.

[0068] 2) Dynamic thinning correction: To obtain the net ablation, the dynamic uplift or subsidence elevation caused by glacial movement must be deducted.

[0069] (4)

[0070] In the formula, This refers to changes in net elevation. This refers to the elevation change between t1 and t2; The rate of rise or fall; This is the time interval between t1 and t2.

[0071] 3) Estimation of ablation: Based on the glacier underlying surface classification results from step 2, different density conversion formulas are set to estimate the glacier ablation:

[0072] (5)

[0073] In the formula, It represents the amount of glacier ablation per unit grid within a specific time period (also known as ablation water equivalent, which is the thickness of the water layer formed per unit area after all the ice and snow have melted and been converted into water). The density of ice; This is the density of water.

[0074] Step 4: Construction of Spatial Continuous Positive Accumulated Temperature Sequence: Temperature data is the main driving meteorological element for the degree-day model. It is matched with the divided glacier grid, and interpolation is used to generate continuous meteorological driving data for each grid. The spatial continuous positive accumulated temperature sequence (PDD) calculated based on this can effectively characterize the distribution of heat available for snow and ice melt in different regions, thus more accurately simulating the spatial pattern of glacier ablation.

[0075] First, temperature downscaling is performed: The MODIS LST data preprocessed in step 1 is downscaled to a resolution consistent with the glacier grid using geographically weighted regression (GWR) or the vertical lapse rate method.

[0076] Then, positive accumulated temperature calculation is performed: the sum of the daily average temperature of each grid that is greater than 0℃ within the ablation period (i.e., from t1 to t2) is statistically analyzed to construct a spatially continuous positive accumulated temperature (PDD) sequence.

[0077] (6)

[0078] In the formula, The positive accumulated temperature value for each grid where the daily average temperature is greater than 0℃; The average temperature of each grid on day d can be obtained from MODIS LST data.

[0079] Step 5: Distributed Initial Day-Degree Factor Inversion: Based on the principle of the day-degree model, the initial day-degree factor for each grid is inverted. The day-degree model is a simplified method for estimating snow and ice melt based on temperature. Its basic principle is that snow and ice melt is linearly related to the cumulative temperature above 0°C (positive accumulated temperature, PDD). This model converts positive accumulated temperature into melt volume through the day-degree factor, thereby achieving a quantitative estimate of the glacier or snow melt process.

[0080] The inversion formula for the initial degree-day factor is:

[0081] (7)

[0082] In the formula: The initial degree-day factor is expressed in mm·d⁻¹·℃⁻¹.

[0083] Because the PDD (particulate matter density) differs between snow-covered areas, unmoored areas, and moored areas, their initial degree-day factors also differ.

[0084] Step 6: Spatial calibration of distributed glacier diurnal factors: Input the daily temperature data of each glacier grid and the initial diurnal factor into the diurnal model to simulate and generate the mass balance data, i.e., ablation, for each grid. The specific formula is:

[0085] * (8)

[0086] In the formula, M represents the ablation amount generated in the simulation.

[0087] Based on this, the simulated mass balance data is compared with the remote sensing mass balance data to calibrate the initial day-degree factor of each grid and generate the final distributed day-degree factor.

[0088] The specific calibration process is as follows: if based on the initial degree-day factor The mass balance data (M) of each simulated grid and the remotely sensed mass balance data are... If the comparison results meet the requirements (error within ±5%), then the day factor is determined; otherwise, the day factor is further calibrated using the residuals of the simulated and remote sensing mass balances, and then compared with the simulated mass balance data. This process is repeated iteratively to find a day factor that meets the error requirements. :

[0089] = (9)

[0090] Finally, distributed, high-precision diurnal factor data that distinguishes between snow-covered areas, non-moored areas (open ice areas), and moored areas are generated, thus obtaining a spatial distribution map of glacier diurnal factors in the study area.

[0091] Example 1

[0092] This embodiment takes the glacier region of the Qinghai-Tibet Plateau as the study area and uses the above method to obtain multi-year average daily temperature factor data. Figure 2 This figure shows a comparison between simulated mass balance data and remote sensing mass balance data in this embodiment. The vertical axis represents simulated mass balance data, and the horizontal axis represents remote sensing mass balance data. The figure shows that the data is concentrated, the effect is good, and it reflects the reliability of the day-degree factor calculation method. Finally, a spatial distribution map of the glacier day-degree factor in the Qinghai-Tibet Plateau glacier region is obtained, as shown below. Figure 3 As shown.

[0093] Finally, it should be noted that the above description is only used to illustrate the technical solution of the present invention and not to limit it. Although the present invention has been described in detail with reference to preferred embodiments, those skilled in the art should understand that modifications or equivalent substitutions can be made to the technical solution of the present invention without departing from the spirit and scope of the technical solution of the present invention.

Claims

1. A distributed method for calculating the diurnal temperature variation factor of glaciers based on multi-source remote sensing data, characterized in that, The method includes the following steps: Step 1: Multi-source data acquisition and preprocessing: Acquire multi-temporal remote sensing data of the study area, including: Landsat 8 / 9 OLI multispectral imagery, MODIS daily surface temperature products, and dual-period elevation data with time spans; acquire daily temperature data from synchronized surface meteorological stations; and preprocess the acquired data. Step 2: Classification and demarcation of glacier underlying surface: Based on the acquired remote sensing data, the glacier area is divided into grids according to the set resolution, and the glacier grids are divided into three categories according to different underlying surfaces: snow cover area, non-moored area, and moored area, and the boundaries of different types of underlying surfaces are demarcated. The specific process for defining the boundaries of different types of underlying surfaces is as follows: First, the snow cover area is demarcated using the Normalized Snow Cover Index (NDSI). The formula for calculating the NDSI is: （1） In the formula, and These represent the reflectance in the green light band and the shortwave infrared band, respectively. When NDSI is greater than a certain set threshold, it indicates a snow-covered area; otherwise, it indicates a non-snow-covered area. For the identification of till-covered areas, the classification discriminant is established as follows: （2） In the formula, For glacier grids; This is a snow-covered area; This is an area not covered by tidal flats; The area is covered by topsoil; The set NDSI threshold is 0.3 to 0.

5. The temperature threshold is 0~0.5; The temperature of the surface of the moraine or ice. Step 3: Grid-scale Glacier Ablation Estimation: By distinguishing between snow-covered areas, non-moored areas, and moored areas, the elevation data of the corresponding grids is determined, thereby determining the elevation difference of each grid at different time periods. Finally, the remote sensing mass balance data for each grid, i.e., the glacier ablation amount, is generated. This specifically includes the following steps: 1) Calculate elevation changes: Difference is performed between two periods of elevation data to obtain the elevation changes at the grid scale. : （3） In the formula, , The elevation data are for times t1 and t2, respectively. 2) Dynamic thinning correction: To obtain the net glacier ablation, the dynamic uplift or subsidence caused by glacier movement must be deducted. （4） In the formula, This refers to changes in net elevation. The speed at which the force surges up or sinks; The time interval between t1 and t2; 3) Glacier ablation estimation: Based on the glacier underlying surface classification results from step 2, different density conversion formulas are set to estimate glacier ablation: （5） In the formula, This represents the amount of glacier ablation per unit grid within a specific time period. The density of ice; The density of water; Step 4: Construction of Spatial Continuous Positive Accumulated Temperature Sequence: Matching the data with the divided glacier grid, interpolating to generate continuous meteorological driving data for each grid, and calculating the spatial continuous positive accumulated temperature sequence (PDD) based on this; the specific process includes: First, temperature downscaling is performed: The MODIS LST data preprocessed in step 1 is downscaled to a resolution consistent with the glacier grid using geographic weighted regression or the vertical lapse rate method. Then, positive accumulated temperature calculation is performed: the sum of the daily average temperatures above 0℃ for each grid within the ablation period (t1 to t2) is statistically analyzed to construct a spatially continuous positive accumulated temperature sequence. （6） In the formula, The positive accumulated temperature value for each grid where the daily average temperature is greater than 0℃; The average temperature of each grid cell on day d; Step 5, Distributed Initial Degree-Day Factor Inversion: Based on the degree-day model principle, invert the initial degree-day factor for each grid; the inversion formula for the initial degree-day factor is: （7） In the formula: This is the initial degree-day factor, in mm·d. -1 ·℃ -1 ; Step 6: Spatial calibration of distributed glacier diurnal factors: Input the positive accumulated temperature value of each glacier grid and the initial diurnal factor into the diurnal model to simulate and generate the mass balance data of each grid, i.e., glacier ablation. The specific formula is as follows: * （8） In the formula, M represents the simulated glacier ablation amount; Based on this, the simulated mass balance data is compared with the remote sensing mass balance data to calibrate the initial degree-day factor and generate the final distributed degree-day factor data.

2. The method for calculating the distributed glacier diurnal factor based on multi-source remote sensing data according to claim 1, characterized in that, The time span in the bi-period elevation data with a time span mentioned in step 1 is at least one year.

3. The method for calculating the distributed glacier diurnal factor based on multi-source remote sensing data according to claim 1, characterized in that, The preprocessing of the acquired data described in step 1 specifically includes: performing terrain and atmospheric correction on Landsat image data; performing 3D registration and deviation correction on elevation data from different sources to ensure spatial reference consistency; and performing cloud masking and missing value filling on MODIS daily surface temperature products.

4. The method for calculating the distributed glacier diurnal factor based on multi-source remote sensing data according to claim 1, characterized in that, The specific process of calibrating the initial degree-day factor in step 6 to generate the final distributed degree-day factor data is as follows: If based on the initial degree-day factor... The mass balance data (M) of each simulated grid and the remotely sensed mass balance data are... If the comparison error is within ±5%, the day factor is determined; otherwise, the day factor is further calibrated using the residuals of simulated and remotely sensed mass balance data. Then, simulated mass balance data is compared with remotely sensed mass balance data to iteratively determine a day factor that meets the error requirements. : = （9） Finally, distributed daily time factor data that distinguishes between snow-covered areas, non-moored areas, and moored areas are generated, thus obtaining a spatial distribution map of glacier daily time factors in the study area.

Images