A method for detecting the storage volume of mountain glaciers

Abstract

This invention provides a method for detecting the reserves of mountain glaciers, relating to the field of glacier reserve detection technology. The method includes: acquiring multi-source remote sensing data of a target glacier region and performing collaborative preprocessing on the multi-source remote sensing data to obtain glacier boundary information and ice surface velocity field; based on a glacier dynamics physical model, fusing the glacier boundary information and ice surface velocity field, and calculating the glacier thickness distribution using an improved thickness inversion model; calculating the total glacier volume through spatial integration based on the glacier thickness distribution and glacier boundary information, and converting this to the total glacier reserve by combining the average density of the glacier ice; verifying the inverted glacier thickness distribution using measured thickness data independent of the modeling data, and optimizing the model parameters based on the verification results. This invention, through the deep fusion and verification optimization closed loop of multi-source remote sensing data and an improved glacier dynamics physical model, enables efficient and accurate detection of mountain glacier reserves.

Description

Technical Field

[0001] This invention relates to the field of glacier reserve detection technology, and specifically to a method for detecting the reserve of mountain glaciers. Background Technology

[0002] Currently, the detection of mountain glacier reserves mainly relies on three types of technical means: first, estimating glacier area by manually interpreting optical remote sensing images; second, using field drilling or ground-penetrating radar to measure glacier thickness at single points or in small areas, and then estimating reserves through spatial extrapolation; and third, directly applying empirical area-volume formulas established based on observation data of specific regions for large-scale estimation.

[0003] However, the aforementioned existing technologies have obvious drawbacks: First, field measurement methods are labor-intensive and dangerous, making it difficult to extend to vast and complex high-altitude areas; second, remote sensing interpretation and empirical formulas are limited by the accuracy of data sources, the subjectivity of manual interpretation, and the inherent regional limitations of empirical models, generally resulting in low automation of glacier boundary extraction, insufficient accuracy of reserve estimation, and poor universality, failing to meet the needs of large-scale, high-precision, and automated dynamic monitoring of glaciers. Summary of the Invention

[0004] In order to solve the technical problems in related technologies, the present invention provides a method for detecting the storage capacity of mountain glaciers.

[0005] To achieve the above objectives, the technical solution adopted by the present invention is as follows:

[0006] A method for detecting the storage capacity of mountain glaciers includes:

[0007] Step S1: Acquire multi-source remote sensing data of the target glacier area and perform collaborative preprocessing on the multi-source remote sensing data to obtain glacier boundary information and ice surface velocity field;

[0008] Step S2: Based on the glacier dynamics physical model, the glacier boundary information and the ice surface velocity field are fused together, and the glacier thickness distribution is calculated using an improved thickness inversion model;

[0009] Step S3: Based on the glacier thickness distribution and the glacier boundary information, calculate the total glacier volume through spatial integration, and combine it with the average density of glacier ice to obtain the total glacier reserves;

[0010] Step S4: Validate the glacier thickness distribution retrieved in Step S2 using measured thickness data independent of the modeling data, and optimize the model parameters based on the validation results.

[0011] Optionally, in step S1, acquiring multi-source remote sensing data and performing collaborative preprocessing specifically includes:

[0012] Step S1-1: Acquire multispectral optical remote sensing images of the target glacier area. Based on the glacier's reflectance spectral characteristics and normalized snow cover index, automatically extract the glacier boundary using a machine learning algorithm to obtain the glacier boundary information, wherein the glacier boundary information includes glacier area information.

[0013] Step S1-2: Acquire synthetic aperture radar image pairs of the same target glacier area with short time baselines, invert the two-dimensional motion field of the glacier surface through pixel offset tracking or interferometry, extract the ice surface velocity value along the main glacier flow line, and obtain the ice surface velocity field.

[0014] Optionally, the improved thickness inversion model, by integrating the ice surface velocity value, glacier surface slope and glacier cross-sectional characteristics, and introducing topographic adjustment factors and regional slip coefficients, obtains the glacier thickness distribution based on glacier dynamics principles.

[0015] Optionally, the topographic adjustment factor is a coefficient determined based on the surface curvature and used to correct the influence of topography on ice flow.

[0016] Optionally, the calibration method for the regional slip coefficient C is as follows:

[0017] Obtain at least one set of independent measured glacier thickness data in the target glacier area or in neighboring glacier areas with similar climate and topography.

[0018] Substitute the ice surface velocity, surface slope angle, and stream function corresponding to the measured data into the calculation formula of the improved thickness inversion model, and use the measured thickness as the target value to back-calculate the corresponding C value.

[0019] Statistical analysis is performed on the multiple C values ​​obtained by back-calculation, and the average or median value is taken as the final value of the slip coefficient C of the region.

[0020] Optionally, step S4 specifically includes:

[0021] Step S4-1: Obtain glacier thickness measurement data obtained by ground penetrating radar or airborne thickness measurement methods, which are not involved in the calibration of the slip coefficient C of the region;

[0022] Step S4-2: Compare the model inversion thickness at the measured point coordinates with the measured thickness, and calculate the root mean square error and coefficient of determination to evaluate the accuracy;

[0023] Step S4-3: If the accuracy does not reach the preset threshold, adjust the empirical constant β in the terrain adjustment factor and / or recalibrate the regional slip coefficient C, and repeat steps S2 to S4-2 until the model accuracy meets the requirements.

[0024] Optionally, in step S3, the glacier thickness distribution is spatially integrated within the area determined by the glacier boundary information to obtain the total volume of the glacier.

[0025] Optionally, the total glacier reserves are obtained by multiplying the total glacier volume by the average density of the glacier ice.

[0026] Optionally, the machine learning algorithm is a random forest classifier or a U-Net convolutional neural network.

[0027] Beneficial effects:

[0028] 1. Through the above technical solution, firstly, the method of the present invention adopts a multi-source remote sensing data collaborative preprocessing approach, which significantly reduces the reliance on a single data source and high-intensity manual intervention, thereby enabling the method of the present invention to achieve automated and efficient large-scale glacier reserve surveys. Specifically, the data source of the method of the present invention is multi-source remote sensing data, thus overcoming the inherent defects of single optical images in boundary identification; at the same time, the structured glacier boundary information and ice surface velocity field output after collaborative preprocessing can transform the original, heterogeneous multi-source remote sensing data into standardized inputs that can be directly used for physical model calculations, thus making it suitable for automated batch processing, thereby enabling the method of the present invention to achieve automated and efficient large-scale glacier reserve surveys.

[0029] Secondly, the method of this invention employs a glacier dynamic physical model to fuse information and calculate thickness. This fundamentally transcends reliance on regional empirical relationships, enabling the method to be universally applicable to different regions and types of glaciers, and yielding more reliable thickness inversion results. Specifically, by introducing a glacier dynamic physical model as the core of the calculation, the thickness inversion is based on universal physical laws rather than specific statistical experience, thus enhancing the theoretical foundation and applicability of the method. Simultaneously, fusing glacier boundary information with the ice surface velocity field allows the thickness inversion to be subject to both geometric and kinematic constraints, resulting in thickness distribution results that are more reliable and physically meaningful than those obtained through extrapolation of single information or empirical formulas.

[0030] Third, the method of this invention employs a clear calculation chain of spatial integration and density conversion, making the reserve estimation process deterministic and traceable. This avoids the uncertainty of results caused by fuzzy coefficients and opaque processes in traditional empirical formula methods. Specifically, the method of this invention calculates volume through spatial integration based on thickness distribution and boundary information, and then converts reserves by average density, forming a mathematical and physical process that is entirely based on prior outputs and has clearly defined steps. In this way, the entire reserve result generation chain is clear and verifiable, thereby effectively improving the reliability and credibility of the final reserve data.

[0031] Fourth, the method of this invention establishes a closed-loop step of verification and parameter optimization using independent measured data. This enables the method to self-calibrate based on actual observations, thereby ensuring the robustness of the technical solution in practical applications and the controllability of the final output accuracy. Specifically, the method of this invention verifies the model inversion results by introducing measured thickness independent of the modeling data, which objectively assesses the accuracy of the core steps. Then, it optimizes the model parameters based on the verification results, essentially using the "true value" information of the target region to perform localized calibration of the method. This makes the method an adaptive system that can continuously approximate the real situation, thus ensuring its reliability and accuracy in different application scenarios from a mechanism perspective.

[0032] 2. Other beneficial effects or advantages of the present invention will be described in detail in the specific embodiments. Attached Figure Description

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

[0034] in:

[0035] Figure 1 This is a flowchart illustrating the steps of a method for detecting the storage capacity of mountain glaciers provided in an exemplary embodiment of the present invention;

[0036] Figure 2 This is a comparison table of estimation results between existing technologies and our method. Detailed Implementation

[0037] To facilitate a clearer and more accurate understanding of the technical solutions of this invention by those skilled in the art, the existing related technologies and their technical problems will be described in more detail below.

[0038] Accurately quantifying the ice storage capacity of mountain glaciers is crucial for assessing water resource security in Gansu Province and the Hexi Corridor, and for studying the impacts of regional climate warming. Currently, this field relies on three main technical approaches, the specific applications of which and their inherent limitations in glacier research in Gansu Province are as follows:

[0039] Example 1: Glacier boundary interpretation and area estimation based on optical remote sensing images.

[0040] This method relies on satellite imagery from sources such as Landsat. For example, when conducting long-term monitoring of glaciers in regions like the Qilian and Altun Mountains in Gansu Province, researchers typically use band ratio thresholding methods (such as calculating the NDSI index using shortwave infrared and green light bands) to initially distinguish snow and ice from the images. However, the moraine (rock debris) covering the glacier terminus has a similar spectrum to the surrounding mountains, making it difficult for automated algorithms to distinguish. This necessitates visual comparison and manual correction by personnel on high-resolution imagery such as Google Earth. The entire process is highly subjective and inefficient. For instance, analyzing glacier changes across multiple image periods for a single region can take weeks of manual boundary delineation and inspection, and the results from different interpreters may vary significantly, failing to meet the demands for rapid, automated, and operational monitoring.

[0041] Example 2: Direct detection and spatial interpolation based on ground / airborne geophysical measurements.

[0042] To obtain the core parameter of glacier thickness, scientific research in Gansu Province has widely adopted ground-penetrating radar and emerging airborne remote sensing technologies, including ground-penetrating radar measurements and airborne ice radar measurements.

[0043] For ground-penetrating radar (GPR) measurements, taking the Qilian Mountains' July 1st Glacier as an example, researchers organized a field expedition in the summer of 2015, using GPR to measure along a pre-defined profile and obtain the glacier's thickness distribution. Based on this, they created an ice bed topographic map and calculated the glacier's annual volume to be approximately 0.1129 billion cubic meters. However, this method has significant limitations: field operations are restricted by extreme weather and complex terrain, posing high risks and costs; for glaciers covering several square kilometers, the limited measurement profiles (usually only a few main lines) cannot fully reflect their complex three-dimensional thickness distribution, and subsequent reliance on mathematical methods such as Kriging interpolation to generalize sparse point data to the entire glacier area introduces considerable uncertainty.

[0044] In 2024, my country conducted its first airborne ice radar penetration survey of the Laohugou No. 12 Glacier, the Qiyi Glacier, and the Ningchanhe No. 3 Glacier in Gansu Province, based on an "airborne remote sensing system." This technology, using radar mounted on an aircraft to penetrate the ice, efficiently acquired large-scale, high-precision subglacial topographic data, and was praised as "the first realization of ice thickness measurement in complex valley glaciers under complex terrain conditions." However, its application is still limited by extremely high economic costs and complex airspace coordination, and currently can only be used for a very small number of "typical glaciers," unable to be applied as a universal method for operational monitoring of all glaciers in Gansu Province.

[0045] Example 3: Indirect estimation based on physical models or empirical formulas (to avoid the difficulties of direct measurement, researchers often use models for indirect estimation). This mainly includes physical model methods and empirical formula methods.

[0046] For physical modeling methods, for example, glacier thickness can be inverted using parameters such as glacier surface velocity and slope based on glacier dynamics theories such as the Shallow Ice Approximation (SIA). However, studies have shown that ice storage results calculated based on the SIA model are generally smaller than those obtained by ground-penetrating radar (GPR), and their accuracy is greatly affected by the quality of digital elevation model (DEM) data and the division of calculation sections. Many physical parameters in the model (such as rheological parameters and slip coefficients) can introduce significant biases when local measured data calibration is lacking.

[0047] For empirical formula methods, in regional-scale assessments, "area-volume" statistical empirical formulas based on global or specific regional data are often used directly for rapid estimation. However, glacier morphology varies greatly depending on local climate and topography. Directly applying formulas fitted to glaciers in other regions (such as the Alps) to glaciers in the Qilian Mountains or Altun Mountains, which have vastly different topographical and climatic conditions, will lead to systematic errors and insufficient reliability.

[0048] Overall, the existing technological system faces a fundamental contradiction in addressing the monitoring needs of glaciers in Gansu Province. Field / aerial measurement methods that can provide direct and accurate thickness data (such as Example 2) are severely limited in their spatial and temporal scalability due to cost, security, and accessibility constraints. Meanwhile, remote sensing interpretation and model estimation methods with large-scale observation capabilities (such as Examples 1 and 3) suffer from low automation, reliance on human experience or simplified physical mechanisms, and poor regional applicability of parameters, making it difficult to guarantee the accuracy and reliability of their results.

[0049] In view of this, the present invention provides a novel solution, namely, a method for detecting the reserves of mountain glaciers. The technical concept of the present invention lies in addressing the fundamental contradiction in existing glacier reserve detection methods—the difficulty in simultaneously achieving scalability, automation, and regional universality—by proposing a novel technical approach: "adaptive inversion of a physical model driven by multi-source remote sensing data."

[0050] Specifically, this invention no longer relies on a single, expensive, or highly manual technology. Instead, it combines the advantages of large-scale, low-cost data acquisition with the physical mechanisms of glacier dynamics by synergistically utilizing readily available optical imagery (for automated glacier boundary extraction) and synthetic aperture radar imagery (for ice surface movement inversion). This invention introduces a topographic adjustment factor to quantify the obstructive effect of complex terrain on ice flow and establishes a regional slip coefficient to integrate and calibrate localized physical processes such as base slip. This allows the model to adapt to the complex and varied terrain and climate conditions in areas such as the Qilian Mountains in Gansu Province. Finally, by constructing a complete computational process from automated data input to physical model inversion and parameter iterative optimization, this invention aims to achieve efficient, high-precision, and regionally adaptable operational estimation of large-scale mountain glacier reserves in a standardized and repeatable manner, thereby providing more reliable technical support for water resource management.

[0051] The technical solution of the present invention will be described in detail below with reference to the accompanying drawings.

[0052] like Figure 1 As shown, the present invention provides a method for detecting the storage capacity of mountain glaciers, comprising the following steps:

[0053] Step S1: Acquire multi-source remote sensing data of the target glacier area and perform collaborative preprocessing on the multi-source remote sensing data to obtain glacier boundary information and ice surface velocity field;

[0054] Step S2: Based on the glacier dynamics physical model, the glacier boundary information and the ice surface velocity field are fused together, and the glacier thickness distribution is calculated using an improved thickness inversion model;

[0055] Step S3: Based on the glacier thickness distribution and the glacier boundary information, calculate the total glacier volume through spatial integration, and combine it with the average density of glacier ice to obtain the total glacier reserves;

[0056] Step S4: Validate the glacier thickness distribution retrieved in Step S2 using measured thickness data independent of the modeling data, and optimize the model parameters based on the validation results.

[0057] Through the above technical solution, firstly, the method of the present invention employs a multi-source remote sensing data collaborative preprocessing approach, which significantly reduces the reliance on a single data source and high-intensity manual intervention. This enables the method to achieve automated and highly efficient large-scale glacier reserve surveys. Specifically, the data source of the method is multi-source remote sensing data, thus overcoming the inherent limitations of single optical images in boundary identification. Simultaneously, the structured glacier boundary information and ice surface velocity field output after collaborative preprocessing can transform the raw, heterogeneous multi-source remote sensing data into standardized inputs that can be directly used for physical model calculations. This makes it suitable for automated batch processing, thereby enabling the method of the present invention to achieve automated and highly efficient large-scale glacier reserve surveys.

[0058] Secondly, the method of this invention employs a glacier dynamic physical model to fuse information and calculate thickness. This fundamentally transcends reliance on regional empirical relationships, enabling the method to be universally applicable to different regions and types of glaciers, and yielding more reliable thickness inversion results. Specifically, by introducing a glacier dynamic physical model as the core of the calculation, the thickness inversion is based on universal physical laws rather than specific statistical experience, thus enhancing the theoretical foundation and applicability of the method. Simultaneously, fusing glacier boundary information with the ice surface velocity field allows the thickness inversion to be subject to both geometric and kinematic constraints, resulting in thickness distribution results that are more reliable and physically meaningful than those obtained through extrapolation of single information or empirical formulas.

[0059] Third, the method of this invention employs a clear calculation chain of spatial integration and density conversion, making the reserve estimation process deterministic and traceable. This avoids the uncertainty of results caused by fuzzy coefficients and opaque processes in traditional empirical formula methods. Specifically, the method of this invention calculates volume through spatial integration based on thickness distribution and boundary information, and then converts reserves by average density, forming a mathematical and physical process that is entirely based on prior outputs and has clearly defined steps. In this way, the entire reserve result generation chain is clear and verifiable, thereby effectively improving the reliability and credibility of the final reserve data.

[0060] Fourth, the method of this invention establishes a closed-loop step of verification and parameter optimization using independent measured data. This enables the method to self-calibrate based on actual observations, thereby ensuring the robustness of the technical solution in practical applications and the controllability of the final output accuracy. Specifically, the method of this invention verifies the model inversion results by introducing measured thickness independent of the modeling data, which objectively assesses the accuracy of the core steps. Then, it optimizes the model parameters based on the verification results, essentially using the "true value" information of the target region to perform localized calibration of the method. This makes the method an adaptive system that can continuously approximate the real situation, thus ensuring its reliability and accuracy in different application scenarios from a mechanism perspective.

[0061] The method of the present invention will be further described below with reference to an exemplary embodiment.

[0062] It should be noted that this exemplary embodiment uses the Qiyi Glacier in the Qilian Mountains of Gansu Province as an example for calculation. At the same time, the goal of this exemplary embodiment is to automatically and accurately estimate the ice reserves of the Qiyi Glacier in a specific year (e.g., 2023) without relying on large-scale field exploration.

[0063] Step 1: Collaborative acquisition and preprocessing of multi-source remote sensing data.

[0064] 1. Acquisition of optical remote sensing images and automated extraction of glacier boundaries:

[0065] Data Acquisition: Acquire Sentinel-2 multispectral images covering the July 1 Glacier at the end of the ablation period in 2023 (September).

[0066] Automated processing: The U-Net convolutional neural network model (or random forest classifier) ​​is used for glacier boundary extraction. This model has been trained using glacier images from multiple locations around the world (including the Qilian Mountains) and manually labeled ground truth values, and can effectively distinguish between clean ice, surface moraine-covered ice, and surrounding bare rock.

[0067] Output: After processing, the digital boundary of the July 1st Glacier is obtained, and its area is calculated. This boundary vector diagram contains information about the glacier boundary.

[0068] For example, when acquiring multispectral optical remote sensing images of a target glacier region, the glacier boundary can be automatically extracted using machine learning algorithms based on the glacier's reflectance spectral characteristics and normalized snow cover index, thus obtaining glacier boundary information (including glacier area information). Furthermore, by acquiring synthetic aperture radar image pairs of the same target glacier region with short time baselines, the two-dimensional motion field of the glacier surface can be inverted using pixel offset tracking or interferometry techniques, and the ice surface velocity value along the main glacier flow direction can be extracted to obtain the ice surface velocity field.

[0069] 2. Radar remote sensing data acquisition and ice surface velocity field inversion:

[0070] Data acquisition: Acquire two Sentinel-1 SAR image pairs for the same region during the winter of 2023 (December-January) with a time span of 12 days.

[0071] Collaborative preprocessing: Pixel offset tracking technology was used to process the image pair, inverting the two-dimensional displacement field of the glacier surface. Geometric correction was then performed using a digital elevation model (DEM), ultimately generating ice surface velocity field data along the main glacier flow line. Its units have been standardized to the International System of Units (m / s).

[0072] In this step, it should be noted that, firstly, since the multi-source remote sensing data in the method of this invention includes optical remote sensing images (high resolution) and radar remote sensing data (synthetic aperture radar image pairs with short time baselines), this allows the method to simultaneously acquire two key and complementary types of physical information about the target glacier: its morphology and movement. This overcomes the information blind spots and identification defects inherent in relying on a single type of data. Optical remote sensing images (e.g., Sentinel-2, Landsat, etc.) are sensitive to surface materials and are the most efficient and mature data source for identifying and extracting the spatial extent (i.e., morphology) of glaciers. Simultaneously, synthetic aperture radar image pairs, through interferometry or migration tracking techniques, are extremely sensitive to minute surface deformations or movements, providing a large-scale, non-contact remote sensing method for acquiring the velocity (i.e., movement) of glacier surface movement.

[0073] Secondly, the processing of SAR images in the above embodiments allows the acquisition of ice surface flow velocity to be based on mature and accurate radar remote sensing physical observation principles, thereby providing direct and quantitative dynamic constraints for glacier thickness inversion. (Pixel offset tracking technology measures velocity by tracking the positional movement of pixels in images at different times using cross-correlation algorithms, while interferometry technology detects subwavelength-level distance changes by analyzing radar wave phase differences. Both technologies are based on the physical principles of radar wave propagation and scattering, enabling precise measurement of two-dimensional surface motion.)

[0074] Step 2: Glacier thickness inversion based on improved physical model.

[0075] 1. Input data preparation:

[0076] Extracting glacier surface slope angles from publicly available DEM data (such as ALOS DSM). .

[0077] Based on the glacier boundary and DEM obtained in the first step, the cross-sectional flow function related to the glacier width is calculated. .

[0078] 2. Apply the core thickness inversion formula:

[0079] An improved physical model was used to calculate the ice thickness. For each pixel point on the glacier, the ice thickness H was calculated using the following formula:

[0080] The formula for calculating ice thickness H is:

[0081]

[0082] In the formula, For the inverted glacier thickness, The extracted ice surface flow velocity value, It is the exponent of the law of flow of glacial ice (a dimensionless constant). The density of glacial ice, It is the acceleration due to gravity. This refers to the local slope angle of the glacier surface calculated based on a digital elevation model. The stream function, which is related to the glacier's cross-sectional area, is jointly determined by the glacier boundary information and the digital elevation model. is the cross-sectional shape factor (dimensionless constant). The topographic adjustment factor is the surface curvature. A function is used to correct the effect of terrain on ice flow. The regional slip coefficient is a parameter calibrated from measured data to account for the glacier base slip effect. (This parameter is a comprehensive physical parameter characterizing the contribution of local physical processes such as glacier base slip to the overall flow velocity; its dimensions are...) or (T represents time, M represents mass).

[0083] For example, the value can be: the exponent of the law of glacial ice flow. Density of glacial ice gravitational acceleration Cross-sectional shape factor (Setting an approximate parabolic cross-section for the July 1st Glacier).

[0084] 3. Determine the terrain adjustment factor:

[0085] According to the formula In the formula, These are empirical constants determined based on the topographical complexity of glacier regions. The surface curvature is calculated from the digital elevation model (DEM).

[0086] For the glacial topography of the Qilian Mountains valleys, empirical constants are used. This makes it possible in flat ice areas At steep icefalls or where the terrain changes drastically at the tip of an ice tongue, To correct the obstructive effect of topography on ice flow.

[0087] 4. Calibrate the slip coefficient of the region.

[0088] Regional slip coefficient The calibration method can be as follows: Obtain at least one set of independent measured glacier thickness data in the target glacier area or in a neighboring glacier area with similar climate and topographic conditions; substitute the ice surface velocity, surface slope angle, and stream function corresponding to the measured data into the calculation formula of the improved thickness inversion model, and use the measured thickness as the target value to back-calculate the corresponding thickness. Value; multiple values ​​obtained by reverse calculation Statistical analysis is performed on the values, and the average or median is taken as the regional slip coefficient. The final value.

[0089] In other words, the regional slip coefficient needs to be calibrated using measured data. .

[0090] For example, the measured thickness data of the July 1st Glacier from the 2024 "Gansu Province Typical Glacier Ground Penetrating Radar Measurement" project are used as the calibration basis.

[0091] Parameter back-calculation and determination: Five measured thickness points of GPR were selected from the upper, middle and lower reaches of the glacier, and their coordinates were determined. , , and Substituting into the above formula, we can find the 5 results. The initial values ​​were then used. The geometric mean of these values ​​was taken to ultimately determine the regional slip coefficient applicable to the glaciers on the northern slope of the Qilian Mountains. (A C value of this magnitude indicates that the contribution of slip at the base of the July 1st Glacier is relatively small, with internal deformation being the primary factor.) This allows for the transformation of limited and expensive field measurement data into a generalizable regional physical parameter.

[0092] 5. Calculate the thickness distribution:

[0093] Input all the above parameters and the prepared spatial variable field into the model, calculate pixel by pixel, and finally obtain the ice thickness distribution map of the entire basin of the July 1st Glacier. .

[0094] In this step, it should be noted that...

[0095] First, the formula for calculating ice thickness H ( In this regard, it is based on the physical laws of glacier flow and the principle of mass conservation, by integrating the velocity profile on the glacier cross section and introducing a topographic adjustment factor. The result is obtained after considering the region slip coefficient C, where, The term represents the shear stress that drives the ice to slide down. This refers to the surface flow velocity exhibited under this stress. (Number of 3 is commonly used) is the rheological index describing the non-Newtonian fluid properties of ice. The (n+2) factor is derived from the integral of the velocity profile over the cross-section of the glacier. The term is used to characterize the influence of glacier cross-sectional shape on flow rate. The basic physical relationship of the formula is: under the same stress, the greater the thickness (H), the greater the flow velocity (u); conversely, observing a specific flow velocity allows us to deduce the corresponding thickness. This elevates thickness estimation from purely statistical area experience to a physical inversion based on mechanical principles.

[0096] Meanwhile, topographic adjustment factor The introduction of this method enables the method of this invention to quantify and correct the local effects of complex terrain (such as icefalls, steep slopes, and cirques) on the standard ice flow equation, thereby effectively improving the accuracy of thickness inversion for real mountain glaciers, especially in areas with dramatic topographic relief. (This is because the standard shallow ice approximation formula assumes a gentle ice bed slope and smooth terrain, while the introduced method...) To directly address this limitation, when the terrain changes drastically (large |S| value), This leads to an increase in the denominator of the formula, resulting in a corresponding decrease in the calculated thickness H. Thus, it can be physically correlated with abrupt changes in terrain, such as icefalls, where the additional obstruction effect on ice flow means that the actual ice thickness at these locations is usually less than the estimate under the assumption of smooth terrain.

[0097] Furthermore, the introduction of the regional slip coefficient C allows this method to comprehensively absorb and characterize the influence of key physical processes not explicitly expressed in the formula, such as base slip conditions, ice temperature, and ice crystal structure, through a calibrable parameter. This enables the same physical formula framework to be flexibly adapted to glacier regions under different climatic and geological conditions through parameter adjustments. (Glacier movement is contributed by both internal deformation and base slip; the standard formula primarily describes internal deformation. C, as a comprehensive regional calibration parameter, directly and proportionally affects the final inverted thickness value when its value increases or decreases. For temperate glaciers with significant base slip, such as maritime glaciers, a smaller C value can be obtained through calibration, allowing the model to calculate a larger thickness H at the same observed flow velocity u, because some movement is contributed by slip. Conversely, for polar glaciers with frozen bases, a larger C value can be calibrated. The introduction of C is the core mechanism for achieving regional universality in this formula; it allows the model to approximate the real physical conditions of different regions through parameter localization without changing its physical core.)

[0098] Second, regarding terrain adjustment factors The calculation formula ( In this regard, it establishes a standardized method that maps continuously changing terrain curvature to adjustment coefficients, making the quantification process of terrain impact deterministic and repeatable, thereby eliminating the arbitrariness and ambiguity in handling terrain effects in the model. Specifically, the formula defines terrain impact with a linear relationship. Here, "1" is the baseline term, representing that when the terrain is flat or has a uniform slope (S≈0), the terrain does not produce additional corrections to the basic physical formula. This is a correction term, where β is a positive empirical constant and |S| is the absolute value of the surface curvature. As the curvature of the terrain increases (|S| increases), regardless of whether the terrain is convex or concave (as guaranteed by the absolute value), the correction term increases linearly, leading to… Within the linear framework, the amount of correction to the model for each unit increase in terrain curvature is constant and predictable (determined by β).

[0099] Third, regarding the regional slip coefficient Regarding the calibration method, firstly, it requires "obtaining at least one set of independent measured glacier thickness data in the target glacier area or neighboring glacier areas with similar climate and topographic conditions." This ensures that the parameter calibration process is based on real, regionally representative observations, thereby guaranteeing that the obtained regional slip coefficient C accurately reflects the physical nature of the glacier in that region, rather than theoretical assumptions. Secondly, when all other variables in the formula for calculating ice thickness H ( , , When data (such as the measured thickness, etc.) can be obtained through remote sensing or DEM, the only unknown, C, can be solved by using the measured thickness as the output target H in the formula. This maximizes the value of each expensive field measurement, making it not just represent information from a single point, but a "benchmark point" for calibrating the entire regional model. Finally, individual measured points may be affected by local anomalies (such as ice fissures or bedrock uplifts). The method of this invention obtains "at least one set" of data and back-calculates multiple C values, then determines the final value by taking the average or median. This utilizes the robustness of statistical data to filter out local noise and extract regional commonalities, ensuring that the final C value is not a fragile value dependent on a specific measurement point, but a statistically significant representative parameter that robustly reflects the overall slip characteristics of the region. This effectively improves the stability and reliability of the model in the region.

[0100] Third, in complex mountain glacial environments, the standard shallow ice approximation formula will produce systematic biases due to abrupt changes in terrain and significant differences in bottom slip conditions. This invention addresses this by introducing a terrain adjustment factor. To quantify the terrain barrier effect, and to absorb localized bottom physical processes through the regional slip coefficient C, a dual correction mechanism is achieved, making the model well applicable to complex mountain glacier environments.

[0101] Step 3: Calculation of glacier volume and reserves.

[0102] 1. Calculate the total volume using spatial integration:

[0103] The thickness distribution obtained in the second step At the glacier boundary determined in the first step Perform a double integral within the inner quadrant to calculate the total volume V of the glacier. The formula is: In the formula, This refers to the glacier boundary, that is, the glacier region determined by glacier boundary information.

[0104] In practice, the thickness values ​​of all pixels in the raster-format thickness map can be multiplied by their pixel areas (e.g., 30m × 30m) and then summed to calculate the total volume of the July 1st Glacier. (Based on this total volume, further calculations show that the average thickness of the July 1st Glacier during this period was approximately 45 meters, which is consistent with the understanding of similar glaciers in the region.)

[0105] 2. Ice storage conversion:

[0106] The formula for calculating the total reserves of glaciers is: ,in, The average density of glacier ice can be, for example, taken as 0.85 g / cm³.

[0107] Total reserves (Approximately 117 million tons).

[0108] In this step, it should be noted that the formula for calculating the total volume V of the glacier ( ,or In this regard, it defines and implements the volume summation of the thickness field in continuous space, so that the total volume calculation is completely determined by the known thickness distribution (H) and boundary range (Ω), thereby completely avoiding the additional errors and subjectivity caused by extrapolation based on a small number of point thicknesses or the use of statistical coefficients in traditional methods.

[0109] Step 4: Accuracy verification and iterative optimization.

[0110] 1. Independent verification:

[0111] Use coefficients not involved in step two The data was validated using another set of independent GPR measured thickness points (e.g., data obtained from a tributary on the other side of the glacier).

[0112] Compare the thickness values ​​retrieved by the model at these points with the measured values.

[0113] The calculated root mean square error is Meter, coefficient of determination The value reached 0.89. This indicates that the model-inverted thickness is highly correlated with the measured thickness, and the absolute error is within an acceptable range.

[0114] 2. Parameter optimization:

[0115] The verification results showed that the simulated thickness of the accumulation zone in the upper part of the glacier was systematically too thin.

[0116] Analysis: This may be due to the different physical properties of the snow on the upper part compared to ice.

[0117] Optimization methods can be:

[0118] Step S4-1: Obtain the slip coefficient of the area not involved in the above-mentioned area by ground penetrating radar or airborne thickness measurement. Calibrated glacier thickness measurement data;

[0119] Step S4-2: Compare the model inversion thickness at the measured point coordinates with the measured thickness, and calculate the root mean square error and coefficient of determination to evaluate the accuracy;

[0120] Step S4-3: If the accuracy does not reach the preset threshold, adjust the empirical constant in the terrain adjustment factor. and / or recalibrate the slip coefficient of the region Repeat steps S2 to S4-2 until the model accuracy meets the requirements.

[0121] Based on the above optimization method, the empirical constants in the terrain adjustment factor are... Make fine adjustments, using a smaller [size] in the accumulation area. The value was adjusted (from 0.08 to 0.05) and recalibrated. value.

[0122] Iteration: Repeat step two and this verification step until the error at all verification points meets the preset threshold (e.g., RMSE < 10 meters). After one round of optimization, the final estimated volume is corrected to... .

[0123] In this step, it's important to note that the parameter optimization method firstly uses "measured glacier thickness data obtained through ground-penetrating radar or airborne thickness measurements, which were not used in the calibration of the regional slip coefficient C" as the verification benchmark. This ensures that the model's performance evaluation is based on objective and independent ground truth values, thereby guaranteeing the authority and impartiality of the verification results and avoiding potential optimistic biases from self-assessment. Secondly, "calculating the root mean square error and coefficient of determination to evaluate accuracy" is used as statistical indicators. These two indicators comprehensively and quantitatively diagnose model performance from two dimensions: the magnitude of the deviation and the degree of trend agreement. This not only provides a binary judgment of "whether it meets the standard" but also accurately identifies whether the model is systematically overestimating / underestimating (RMSE reflection) or failing to capture spatial variation patterns. This reflects the fact that the system can provide a clear direction for subsequent optimization to a certain extent. Finally, "if the accuracy does not reach the preset threshold, adjust the empirical constant β in the terrain adjustment factor and / or recalibrate the regional slip coefficient C, and repeat steps S2 to S42 until the model accuracy meets the requirements" allows the system to actively approach the target accuracy through parameter adjustment, thereby ensuring that the final output results of the method of the present invention can reach a controllable and reliable quality standard under different regions and data conditions.

[0124] To visually demonstrate the effectiveness of this invention, the estimation results of this method are compared with those of traditional methods. For a detailed comparison, please refer to the appendix. Figure 2 .

[0125] The above are merely specific embodiments of the present invention, but the scope of protection of the present invention is not limited thereto. Any variations or substitutions within the technical scope disclosed in the present invention should be included within the scope of protection of the present invention. Therefore, the scope of protection of the present invention should be determined by the scope of the claims.

Claims

1. A method for detecting the storage capacity of mountain glaciers, characterized in that, include: Step S1: Acquire multi-source remote sensing data of the target glacier area and perform collaborative preprocessing on the multi-source remote sensing data to obtain glacier boundary information and ice surface velocity field; Step S2: Based on the glacier dynamics physical model, the glacier boundary information and the ice surface velocity field are fused together, and the glacier thickness distribution is calculated using an improved thickness inversion model; Step S3: Based on the glacier thickness distribution and the glacier boundary information, calculate the total glacier volume through spatial integration, and combine it with the average density of glacier ice to obtain the total glacier reserves; Step S4: Validate the glacier thickness distribution retrieved in Step S2 using measured thickness data independent of the modeling data, and optimize the model parameters based on the validation results.

2. The method for detecting the storage capacity of mountain glaciers according to claim 1, characterized in that, In step S1, acquiring multi-source remote sensing data and performing collaborative preprocessing specifically includes: Step S1-1: Acquire multispectral optical remote sensing images of the target glacier area. Based on the glacier's reflectance spectral characteristics and normalized snow cover index, automatically extract the glacier boundary using a machine learning algorithm to obtain the glacier boundary information, wherein the glacier boundary information includes glacier area information. Step S1-2: Acquire synthetic aperture radar image pairs of the same target glacier area with short time baselines, invert the two-dimensional motion field of the glacier surface through pixel offset tracking or interferometry, extract the ice surface velocity value along the main glacier flow line, and obtain the ice surface velocity field.

3. The method for detecting the storage capacity of mountain glaciers according to claim 2, characterized in that, The improved thickness inversion model, by integrating the ice surface velocity value, glacier surface slope and glacier cross-sectional characteristics, and introducing topographic adjustment factors and regional slip coefficients, obtains the glacier thickness distribution based on glacier dynamics principles.

4. The method for detecting the storage capacity of mountain glaciers according to claim 3, characterized in that, The topographic adjustment factor is a coefficient determined based on the surface curvature and used to correct the influence of topography on ice flow.

5. The method for detecting the storage capacity of mountain glaciers according to claim 3, characterized in that, The region slip coefficient The calibration method is as follows: Obtain at least one set of independent measured glacier thickness data in the target glacier area or in neighboring glacier areas with similar climate and topography. Substituting the measured ice surface velocity, surface slope angle, and stream function into the calculation formula of the improved thickness inversion model, and using the measured thickness as the target value, the corresponding thickness is calculated. value; The multiple results obtained from the reverse calculation Statistical analysis is performed on the values, and the average or median value is taken as the slip coefficient of the region. The final value.

6. The method for detecting the storage capacity of mountain glaciers according to claim 4 or 5, characterized in that, Step S4 specifically includes: Step S4-1: Obtain the slip coefficient of the area not involved in the above-mentioned area by ground penetrating radar or airborne thickness measurement. Calibrated glacier thickness measurement data; Step S4-2: Compare the model inversion thickness at the measured point coordinates with the measured thickness, and calculate the root mean square error and coefficient of determination to evaluate the accuracy; Step S4-3: If the accuracy does not reach the preset threshold, adjust the empirical constant in the terrain adjustment factor. and / or recalibrate the slip coefficient of the region Repeat steps S2 to S4-2 until the model accuracy meets the requirements.

7. The method for detecting the storage capacity of mountain glaciers according to claim 3, characterized in that, In step S3, the glacier thickness distribution is spatially integrated within the area defined by the glacier boundary information to obtain the total volume of the glacier.

8. The method for detecting the storage capacity of mountain glaciers according to claim 7, characterized in that, The total glacier reserves are calculated by multiplying the total glacier volume by the average density of the glacier ice.

9. The method for detecting the storage capacity of mountain glaciers according to claim 1, characterized in that, The machine learning algorithm is either a random forest classifier or a U-Net convolutional neural network.

Images