Glacier surface moraine identification method based on multi-source remote sensing parameters

Abstract

The invention relates to the technical field of glacier information processing, in particular to a glacier surface moraine identification method based on multi-source remote sensing parameters, which comprises the following steps: performing radiation-atmosphere joint correction and geometric-reprojection on optical, thermal infrared, synthetic aperture radar and digital elevation data to generate physical consistency data; deducing the surface temperature under the driving of short-wave albedo, downward radiation and turbulent flux by using a cross-modal visual conversion network and a symmetric neural operator, and obtaining a thickness-existence degree posterior through a normalized flow; an uncertainty field is constructed according to information entropy and variance, a representative pixel unmanned aerial vehicle radar sampling and value measurement backflow fine tuning model is selected by means of quantum annealing, and then a spatial continuous thickness field is output through anisotropic diagram diffusion. According to the method, the future thickness is integrated according to mass conservation in combination with day-by-day weather forecast, and the future thickness is coupled with the existence degree index to form a burst risk grid, so that high-precision and iterable moraine dam early warning is realized.

Description

Technical Field

[0001] The present invention relates to the technical field of glacier information processing, and in particular to a glacier surface moraine identification method based on multi-source remote sensing parameters. Background Art

[0002] The front of alpine glaciers is often covered by a debris-ice mixture (surface moraine). Debris thickness controls ice surface energy absorption and local ablation rates, and is also a trigger for moraine dam erosion and breaches. Accurately understanding its spatial distribution and multi-year evolution is crucial for disaster warning, regional water resource assessment, and Quaternary environmental reconstruction. Due to the rugged glacial topography and severe cold climate, drilling or trenching methods are costly and have limited coverage, making remote sensing technology the only viable way to conduct macroscopic monitoring of surface moraines. However, a single sensor can often only capture one aspect of the debris layer: the visible-near infrared band reflects color differences but can easily confuse coarse-grained debris with thin snowpack; thermal infrared records temperature but is affected by diurnal variations and shadows; and synthetic aperture radar is sensitive to roughness and humidity but lacks spectral discrimination. Therefore, multi-source remote sensing fusion is considered a way to break through existing bottlenecks. Summary of the Invention

[0003] To address the numerous issues with the aforementioned existing technologies, the present invention provides a glacier moraine identification method based on multi-source remote sensing parameters. This method physically unifies multi-domain remote sensing imagery, extracts cross-modal features using a two-branch visual transformation network, and infers instantaneous temperature using a symmetric neural operator coupled with shortwave albedo and turbulent flux. A posteriori for thickness and presence are then output using affine coupled flows. The uncertainty field is sampled using quantum annealing to select representative pixels. The measured thickness is then used to fine-tune the model, and a continuous thickness field is generated using anisotropic diffusion. Finally, the future thickness is integrated under daily meteorological constraints to construct a risk index.

[0004] A method for identifying glacier moraines based on multi-source remote sensing parameters includes the following steps:

[0005] Synchronously acquire optical and thermal infrared remote sensing images, synthetic aperture radar remote sensing images and digital elevation data, perform radiometric-atmospheric joint correction and geometric-reprojection processing, and generate physically consistent data;

[0006] The physical consistency data is input into a cross-modal encoder, and combined with the physical driving field composed of shortwave albedo, downwelling radiation and turbulent flux, the surface temperature is inferred using a symmetric neural operator, and the temperature and joint characteristics are used as conditions to input into a differentiable probability model to obtain the posterior data of surface moraine thickness and surface moraine presence;

[0007] An uncertainty field is constructed based on the posterior data, high entropy and spatially correlated pixels are selected through quantum optimization to perform field sampling, the measured thickness is used to fine-tune the differentiable probability model and the symmetric neural operator, and a spatially continuous thickness field is formed through anisotropic graph diffusion regularization;

[0008] The symmetric neural operator with the spatially continuous thickness field and frozen weights is used to fuse the daily weather forecast to calculate the net radiation of the pixel, and the future thickness is obtained by integrating it to a set number of years according to the law of conservation of mass, and the burst risk index raster data is generated based on the presence of surface moraines and the future thickness.

[0009] Preferably, the cross-modal encoder adopts a visual conversion network, the optical and thermal infrared remote sensing images are divided into non-overlapping patches of fixed length at the input end and linearly mapped to form a first feature sequence, the synthetic aperture radar remote sensing images and digital elevation data are divided into the same length to form a second feature sequence, and the first feature sequence and the second feature sequence are respectively processed by multi-layer self-attention modules and spliced into a joint feature vector on the hidden dimension.

[0010] Preferably, the shortwave albedo is calculated based on the normalized difference between the reflectivity of the green light band and the reflectivity of the shortwave infrared band, the downlink radiation is interpolated through the meteorological forecast field to a spatial resolution consistent with the physical consistency data, and the turbulent flux is calculated based on the near-surface meteorological elements through the friction velocity method.

[0011] Preferably, the symmetric neural operator is composed of multiple layers of Fourier integral operators connected in series, and the convolution kernels of each layer are set to a symmetrical structure according to the real part and the imaginary part in the frequency domain.

[0012] Preferably, the differentiable probability model is composed of multiple levels of affine coupled transformation units connected in series, each level of affine coupled transformation unit uses an alternating mask method to separate input channels, and the scale function and the displacement function are implemented through a multi-layer perceptron.

[0013] Preferably, the uncertainty field is obtained by linearly combining the information entropy of the pixel surface presence and the pixel thickness variance according to preset weights.

[0014] Preferably, quantum optimization uses quantum annealing to solve the quadratic disordered Boolean optimization model, sets the pixel uncertainty as the linear term coefficient, and sets the symmetric neural operator convolution kernel spectral similarity as the quadratic term coupling coefficient to obtain a set of pixels with high entropy and strong spatial correlation.

[0015] Preferably, on-site active sampling collects the pixel set in sequence through a traveling salesman path constructed by a drone. The drone is equipped with a surface wave radar to obtain a radar profile and invert the thickness according to predetermined electromagnetic velocity parameters. At the same time, it is equipped with a dual-band thermal infrared sensor to obtain brightness temperature information. The radar inversion thickness is used to fine-tune the differentiable probability model and the symmetric neural operator.

[0016] Preferably, the anisotropic graph diffusion regularization model establishes an undirected graph structure on a uniform resolution grid, the main diagonal elements of the diffusion tensor are set according to the weighted values of the pixel slope and the pixel uncertainty, and the explicit difference format is iterated until the average relative difference between two adjacent thickness fields is less than a preset threshold, and the spatially continuous thickness field is output.

[0017] Preferably, the collapse risk index is calculated based on the exponential relationship between the presence of moraines on the pixel surface and the thickness at a future time. The thickness at a future time is obtained by integrating the mass conservation equation daily with the preset average density of the debris layer and the latent heat constant as parameters. The risk index is output in a raster form for risk grading.

[0018] Compared with the prior art, the advantages and beneficial effects of the present invention are:

[0019] This invention uses a cross-modal visual conversion network + symmetric neural operators to achieve deep fusion of optical, thermal infrared, radar and terrain information, solving the problem that single-band methods cannot finely identify thin debris.

[0020] The present invention achieves a posteriori inference of thickness under energy budget constraints through a closed loop of physical driving field and differentiable probability model, overcoming the defect that the empirical coefficient fails with the scene.

[0021] The present invention uses quantum annealing active sampling and anisotropic graph diffusion regularization to output a spatially continuous and noise-suppressed thickness field, avoiding banding artifacts and significantly reducing field workload. BRIEF DESCRIPTION OF THE DRAWINGS

[0022] Figure 1 is a flow chart of the present invention;

[0023] Figure 2 This is a data processing flow chart of the present invention;

[0024] Figure 3 This is a structural diagram of the symmetric neural operator and physical drive fusion in the present invention. DETAILED DESCRIPTION

[0025] Hereinafter, embodiments of the present disclosure will be described with reference to the accompanying drawings. However, it should be understood that these descriptions are merely exemplary and are not intended to limit the scope of the present disclosure. In the following detailed description, for ease of explanation, many specific details are set forth to provide a comprehensive understanding of the embodiments of the present disclosure.

[0026] like Figure 1-Figure 2 As shown in FIG, a glacier moraine identification method based on multi-source remote sensing parameters includes the following steps:

[0027] Synchronously acquire optical and thermal infrared remote sensing images, synthetic aperture radar remote sensing images and digital elevation data, perform radiometric-atmospheric joint correction and geometric-reprojection processing, and generate physically consistent data;

[0028] Unified preprocessing of multi-source remote sensing parameters is the starting point of the entire glacier moraine identification process. Its goal is to map raw data with different observation phases, imaging mechanisms, spatial resolutions, and coordinate bases into the same physical semantic space, thereby laying a comparable numerical foundation for subsequent cross-modal feature encoding and coupled inference. This paper introduces a two-stage strategy of "radiance-atmosphere joint correction" and "geometry-reprojection unified processing" in this step. Compared with the traditional pipelined, item-by-item correction method, this can eliminate the radiation amplitude drift and spatial distortion superposition errors between multi-source data at one time, significantly reduce the cross-registration errors between optical, thermal infrared, and synthetic aperture radar, and mitigate the irreversible loss of image sharpness caused by repeated resampling.

[0029] In the radiation-atmosphere joint correction stage, the top radiance normalization of the optical and thermal infrared remote sensing images is first performed based on the orbital epoch and the solar altitude angle, and the digital quantization value of the radiometer is converted into TOA (top atmosphere) radiance; then the radiation transfer model is used to remove the atmospheric scattering and absorption effects from the TOA radiance, and directly push it back to the surface reflectivity layer. Since the atmospheric pressure in high-altitude glacier areas is relatively low and the water vapor content is thin, if the empirical formula of low-altitude areas is still used, systematic over-compensation errors will occur. The present invention introduces an atmospheric pressure and water vapor profile parameter extrapolation module in the solution of the atmospheric transmission coefficient, and calibrates the real-time water vapor column content by integrating microwave radiometer data, so that the correction coefficient is more suitable for the high-altitude and thin atmosphere. For synthetic aperture radar images, the present invention adopts a radiation-geometry integrated calibration process to integrate thermal noise removal, Beta0 calibration and slant range-ground range correction into one operation to avoid the edge stripe effect caused by multi-level resampling.

[0030] The unified geometry-reprojection processing stage uses the digital elevation model generated by laser altimetry and stereo optical images as the terrain benchmark to perform terrain correction on optical and radar images respectively. Optical images use polynomial image equations to solve the heading error and perform adjustment, while radar images use the time-space baseline-elevation fitting method to compensate for terrain distortion. All data are finally reprojected to a unified UTM coordinate system and rasterized to a resolution of 20 meters. In order to prevent frequency domain aliasing introduced by repeated interpolation, optical images use a cubic convolution kernel and radar images use a sinc-Lanczos finite kernel to ensure that the high-frequency texture in the rectangular coordinate grid still has sufficient energy. After this stage, all observations are organized into a sequence of row pixels and stored in X phys Matrix, the column-wise features of the matrix are arranged strictly in the order of band number, radiation amount, backscatter coefficient and terrain index, which is convenient for direct slicing and reading in the cross-modal encoder.

[0031] In terms of effectiveness, the combined radiation-atmosphere correction keeps the difference in visible and near-infrared reflectance observed in multiple temporal phases on the same day within 5%, and the thermal infrared brightness temperature drift within one Kelvin. After unified geometry-reprojection processing, the average deviation between pixel misalignment caused by mountain radar image distortion and optical pixel boundaries in SAR images is reduced to less than one pixel. Using the glaciers on the eastern slope of Mount Gongga as an experimental area, the preprocessing method employed achieved a mean square error (MSE) between the simultaneously acquired Sentinel-2MSI green light reflectance and the ground-based reflectance measured by the airborne spectrometer, less than half that of the traditional six-stage correction process. The SAR backscatter coefficient is highly consistent with the measured incident angle scattering response curve based on ground angular reflectors, providing a stable baseline for subsequent frequency-domain estimation of debris roughness.

[0032] In the embodiment, the research team selected a typical glacier front area on the Qinghai-Tibet Plateau. They first deployed corner reflectors and water surface diffuse reflectors in the field to evaluate the effects of geometric and radiometric dual corrections. They then acquired Sentinel-2MSI, LandsatOLI, Sentinel-1IWSLC, and ICESat-2ATL06 data with a time interval of no more than six hours. After applying the joint correction process of the present invention, the optical and thermal infrared bands were uniformly converted to surface reflectivity, and the SAR images were slant range-to-ground range corrected and calibrated to σ0. Finally, X phys The matrix is input into the cross-modal encoder, which shows significant diagonal mutual activation of spectral-radar information in the first layer self-attention weight heat map, proving that high consistency preprocessing effectively improves the coupling degree of different modalities in the feature space.

[0033] like Figure 3 As shown, the physical consistency data is input into a cross-modal encoder, combined with the physical driving field composed of shortwave albedo, downwelling radiation and turbulent flux, and the surface temperature is inferred using a symmetric neural operator. The temperature and joint features are used as conditions to input into a differentiable probability model to obtain the posterior data of surface moraine thickness and surface moraine presence.

[0034] After completing the radiation-atmosphere joint correction and geometry-reprojection, the physical consistency data is organized into the matrix X phys . The order of arrangement of the matrix column features is: green light reflectance, near-infrared reflectance, short-wave infrared reflectance, thermal infrared brightness temperature, dual-polarization backscatter coefficient, digital elevation, slope and curvature. This order ensures that subsequent slicing operations can extract spectral subsets or microwave-topographic subsets with constant time complexity. Each row corresponds to a single pixel, and the row index strictly corresponds one-to-one to the 20-meter resolution grid. Then it enters the core computing link of cross-modal encoding-physical coupling-probabilistic inference.

[0035] The cross-modal encoder utilizes a two-branch visual transformation network. The first branch receives a first feature sequence formed by optical and thermal infrared remote sensing imagery; the second branch receives a second feature sequence formed by synthetic aperture radar backscatter coefficients and digital elevation metrics. Both branches employ a fixed-length patch splitting strategy. The patch vectors are linearly expanded to a uniform dimension before being fed into a multi-layer self-attention module. The key-query-value operations within the self-attention module are all performed on the same scale space, enabling a global interaction between spectral information and backscatter texture. The latent vectors output by the two branches are concatenated in the hidden dimension to form a joint feature vector z. The cross-modal encoder is trained using a contrastive loss. Positive samples consist of bimodal features from the same pixel, while negative samples are randomly paired within the batch. The contrastive loss ensures that the cross-modal feature embedding space remains aligned.

[0036] The joint eigenvector needs to be concatenated with the pixel-level physical driving field before entering the physical inference phase. The physical driving field includes shortwave albedo α, downwelling shortwave radiation S↓, downwelling longwave radiation L↓, sensible heat flux H, and latent heat flux LE. α is calculated using the following formula:

[0037]

[0038] R G Indicates the surface reflectance in the green light band, R SWIR Represents the surface reflectivity in the shortwave infrared band. S↓ and L↓ are obtained by interpolating the short-beam top radiation and long-wave downlink radiation of the meteorological model. H and LE are calculated from ground meteorological elements according to the friction velocity method. All components of the physical driving field are spatially related to X phys Alignment, using the nearest moment in time.

[0039] The symmetric neural operator is responsible for converting (z,D phys ) is mapped to the instantaneous surface temperature field. The operator network consists of four layers of Fourier integral operators connected in series. Each layer first performs a fast Fourier transform to map the input to the frequency domain; then the convolution kernel κ(f) is applied in the frequency domain. The convolution kernel satisfies the symmetry constraint of κ(f) = κ(-f), ensuring that the time domain operator matrix obtained after the inverse transform is symmetric and positive definite. This design conforms to the physical properties of the heat diffusion equation. Network training simultaneously minimizes information alignment loss and physical residual loss. The physical residual loss is based on the energy balance relationship:

[0040] Q net =(1-α)S↓+L↓-εσT 4 -H-LE

[0041] Where T is the predicted surface temperature, ε is the surface emissivity towards the band, and ε is the Stefan constant. The second-norm residual of is used to make the network output consistent with the heat diffusion-radiation balance coupling equation.

[0042] The differentiable probability model adopts the form of multi-level affine coupled normalized flow. The input vector is the concatenated [z,T]. The normalized flow starts from the prior distribution (lognormal-Bernoulli mixture) and is mapped to the posterior distribution through a reversible affine transformation. The affine coupling block uses an alternating mask to divide the channels into two groups. The first group remains unchanged, and the second group is affine transformed according to the scale function s(·) and displacement function t(·) mapped by the perceptron of the first group, thereby ensuring that the Jacobian determinant is easy to calculate. The model maximizes the lower bound of evidence and jointly optimizes the contrast loss, physical residual loss and variational loss. The final output is the posterior mean μ d ,variance and the existence probability p c .

[0043] The motivation for introducing a posteriori uncertainty in this link is to provide a quantitative basis for active sampling selection. Information entropy p c log p c and variance The pixel uncertainty field is comprehensively constructed. Relying on the spatial correlation matrix generated by the convolution kernel spectrum of a symmetric neural operator, a quadratic unordered Boolean optimization model can be solved on a quantum annealing machine to select a subset of pixels that is both informative and spatially representative. This active sampling strategy significantly reduces fieldwork workload and significantly reduces the noise in the posterior thickness field.

[0044] Example 1 was carried out in the extremely cold region of the eastern Pamir Plateau. The study area used physical consistency data constructed by Sentinel-2, MODIS thermal infrared, Sentinel-1, ALOS-2 and ICESat-2 as input. Under the cross-validation framework, after visualizing the self-attention weights, it was found that compared with the ordinary parallel convolutional network without physical field constraints, the weight matrix output by the symmetric neural operator showed smoother radial attenuation in the low-frequency and high-frequency regions, indicating that the model has intrinsically learned the stationary characteristics of the debris layer-temperature field. After testing at 500 GPS-GPR common points, the root mean square error of the thickness was reduced from 3 meters of the traditional empirical formula regression to 1.6 meters.

[0045] Example 2 targets the southern slope glacier front with significant seasonal freeze-thaw. Remote sensing data from April and August were input into the link of the present invention respectively, and the thickness difference was inferred and compared with the results of the rod ablation monitoring. The results showed that the correlation coefficient between the thickness change trend output by the model and the measured ablation amount reached 0.89, verifying that the symmetric neural operator can still respond accurately under conditions of drastic changes in seasonal radiation flux. The same experiment was compared without adding a physical driving field, and the correlation coefficient dropped to 0.65, confirming the role of the physical driving field in improving the generalization ability of the model.

[0046] Preferably, the cross-modal encoder adopts a visual conversion network, the optical and thermal infrared remote sensing images are divided into non-overlapping patches of fixed length at the input end and linearly mapped to form a first feature sequence, the synthetic aperture radar remote sensing images and digital elevation data are divided into the same length to form a second feature sequence, and the first feature sequence and the second feature sequence are respectively processed by multi-layer self-attention modules and spliced into a joint feature vector on the hidden dimension.

[0047] The cross-modal encoder is designed to fuse the spectral and thermal radiation information carried by optical and thermal infrared remote sensing imagery with the structural, roughness, and topographic information carried by synthetic aperture radar remote sensing imagery and digital elevation data within a unified feature space. Remote sensing imagery differs in its native resolution and imaging mechanism. Direct splicing can easily lead to gradient updates dominated by one modality, resulting in information masking. This paper employs a visual transformer network as a cross-modal encoder, achieving scale normalization and semantic alignment of the two modalities through unified patch segmentation and self-attention mechanisms.

[0048] First, a fixed-length two-dimensional non-overlapping patch is created in the optical and thermal infrared channels. The patch side length is an integer multiple of the 20-meter grid unit to ensure that the patch coverage is consistent with the subsequent self-attention window. The original pixels are expanded in row-major order and filled into the matrix, and the absolute position code is inserted in the patch dimension. Each patch is mapped to the dimension d by a linear transformation. model The vectors of are used to form the first characteristic sequence. The synthetic aperture radar backscatter coefficients and digital elevation data undergo the same partitioning and mapping process to form the second characteristic sequence. The patch sizes of the two sequences at the input are consistent with the linear mapping matrix, fundamentally eliminating the scale mismatch caused by resolution differences.

[0049] The sequence enters a multi-layer self-attention module. Each layer of the self-attention module consists of multiple query-key-value mappings. The query vector With key vector After doing the dot product operation, the scaling factor And get the weight matrix W through Softmax, the weight matrix and the value vector Multiplication yields the context vector C. After concatenating the multi-head outputs in the channel dimension, they are linearly mapped back to the original dimension and the residual is added. Compared to the convolution kernel architecture, which has a fixed receptive field, the self-attention module automatically allocates long-range and short-range dependencies based on the scene during training, allowing spectral-thermal and microwave-topographic information to interact at different scales.

[0050] The present invention introduces contrast loss in the cross-modal alignment stage. The vector output by the optical-thermal infrared branch is denoted as z opt , the vector output by the radar-terrain branch is recorded as z sar When the batch size is B, define:

[0051]

[0052] sim(·,·) is the cosine similarity function, and τ is the temperature hyperparameter. By introducing this loss, bimodal features of the same pixel are brought closer together in the embedding space, while bimodal features of different pixels are pushed further apart, achieving consistent cross-modal embedding. Through backpropagation, the discrimination boundaries between optical-thermal and radar-terrain features gradually converge to a common subspace, enabling subsequent symmetric neural operators to absorb information from all modalities within a single vector space.

[0053] The joint eigenvector z is directly concatenated in the hidden dimension by the outputs at the ends of the two branches. After concatenation, it is compressed to the same length as the physical driving field dimension through a fully connected layer. The concatenation operation retains the important characteristic dimensions of each modality, while the compression operation ensures the consistency of the input structure of the subsequent symmetric neural operators. To avoid excessive differences in the gradient update speeds of the optical-thermal branch and the radar-terrain branch, the present invention applies a large forward scaling factor to the radar-terrain branch at the beginning of training to keep its gradient norm on the same order of magnitude as the optical-thermal branch. After 5,000 optimization steps, the scaling factor is removed when the loss converges stably to ensure the long-term stability of the model.

[0054] An additional advantage of the visual transformation network is the explicit interpretability provided by the self-attention weight matrix. In experiments at the glacier front, the first-layer attention head generated long-range coupling between the backscatter coefficient of the slope region and the green reflectance of the shadowed region above, reflecting the intrinsic relationship between debris accumulation and slope aspect and illumination conditions. This interpretability is used in the active sampling phase to construct the space-frequency correlation matrix, allowing the quantum annealing process to incorporate physical prior information rather than relying solely on statistically independent entropy screening.

[0055] Example 1 was conducted on the main ice ridge of the West Kunlun Mountains. The research team uniformly extracted one thousand pixels from the true value map as the evaluation set, and used the traditional convolutional neural network as a control. After using the cross-modal encoder of the present invention, the root mean square error of the thickness on the evaluation set dropped from 3.2 meters to 1.9 meters; the active sampling points selected after the information entropy-variance uncertainty mapping reduced the field workload by 60% compared with random sampling, and at the same time reduced the noise patch area ratio of the final spatial continuous thickness field from 12% to 4%. Example 2 was conducted on the northern slope of the Altun Mountains, where cloud shadows and deep snow cover existed in the experimental area. The present invention performed global self-attention alignment on the spectral features of the cloud shadows and the roughness features of the radar spots, keeping the thickness error within two meters on the test set, while the error of the control method increased to more than four meters.

[0056] The cross-modal encoder also captures directional roughness information through multi-head attention. In a field-measured reflectivity-roughness test site, the optical-thermal branch focused on patches corresponding to high-value areas of longwave radiation in the debris layer, while the radar-topography branch focused on areas with maximum directional backscattering. Fusion of these two branches in a high-dimensional space significantly improved the ability to discriminate between coarse-grained and fine-grained mixed scenes, providing sufficient conditions for the subsequent joint thickness-presence distribution of the probabilistic model. The presence of long diagonal activation bands in the self-attention heatmap also indicates that the model tends to treat pixels in the same slope aspect or elevation band as global context, consistent with the physical law that glacial moraines are distributed in zonal patterns along topography.

[0057] Preferably, the shortwave albedo is calculated based on the normalized difference between the reflectivity of the green light band and the reflectivity of the shortwave infrared band, the downlink radiation is interpolated through the meteorological forecast field to a spatial resolution consistent with the physical consistency data, and the turbulent flux is calculated based on the near-surface meteorological elements through the friction velocity method.

[0058] Shortwave albedo, downwelling radiation, and turbulent fluxes jointly determine the instantaneous energy budget of the debris layer-ice surface system. This invention provides a unified and interconnected acquisition path for these three physical quantities within a multi-source remote sensing framework, avoiding the inefficiency of point-by-point field measurements while ensuring a precise correspondence with the high-spatial-resolution remote sensing grid. The following describes this approach from three perspectives: principle, implementation details, and results.

[0059] Shortwave albedo is the ratio of the surface's instantaneous scattering of solar shortwave radiation. The surface of glacial debris has a strong absorption of visible light, but shows a higher scattering of shortwave infrared. This invention uses the green light band reflectivity R G and the shortwave infrared band reflectivity R SWIR The normalized difference of establishes an approximate parameter-free albedo estimation formula:

[0060]

[0061] Where α is the shortwave albedo, which has no unit; R G With R SWIR are the green and shortwave infrared reflectances of the surface layer after combined radiation-atmosphere correction. The formula is derived based on the empirical relationship that the difference in scattering-absorption between the two bands is inversely proportional to the integrated scattering ratio under the principle of energy conservation. Compared to polynomial models that require multi-band regression, this formula relies solely on two high-signal-to-noise ratio bands, maintaining stable output even in the presence of cloud cover or partial sensor failure. In the examples, the α values calculated from the green and shortwave infrared reflectances of the Sentinel-2 multispectral instrument using this formula were compared with the integrating sphere observations from a ground-based spectrometer, yielding a mean absolute error of less than 0.03.

[0062] Downward radiation consists of a shortwave component S↓ and a longwave component L↓, which are derived from solar radiation and atmospheric thermal radiation. The present invention uses numerical weather forecast products with a resolution better than three thousand meters, and projects the forecast field to a twenty-meter grid through bilinear interpolation. In order to avoid overestimation of shortwave radiation caused by mountain shadows, the algorithm calculates the solar incidence angle correction coefficient for each pixel according to the slope-solar azimuth after interpolation; if the corrected incidence angle is less than the horizon, the shortwave radiation is automatically set to zero. In the scene on the south slope of Mount Gongga in Example 1 of the specification, this operation reduces the estimation error of S↓ in the shadow area from an uncorrected 400 watts per square meter to less than 100 watts per square meter. Longwave radiation is not sensitive to terrain shadows and can be kept consistent with the interpolation value. After interpolation, the radiation field completely corresponds to the row and column coordinates of the physical consistency data, ensuring that the symmetric neural operator can read all physical quantities at once.

[0063] Turbulent flux refers to the sensible heat flux H and latent heat flux LE, which characterize the sensible heat and water vapor exchange between air and debris-ice surface. The present invention uses the friction velocity method to calculate. First, the near-surface temperature T is extracted based on the weather forecast product. a Specific humidity q a The horizontal wind speed U is combined with the dual-polarization synthetic aperture radar backscatter coefficient to estimate the surface roughness length z0: the backscattering in the debris roughness area is stronger, corresponding to a larger z0; the backscattering in the bare ice area is weaker, corresponding to a smaller z0. The friction velocity u is then solved using the classic Prandtl-von-Kármán relationship. * The expressions for sensible heat flux and latent heat flux are:

[0064] H=ρc p u * (T s -T a )

[0065] LE=ρL v u * (q s -q a )

[0066] Where ρ is the air density, c p is the specific heat capacity at constant pressure, L v is the latent heat of evaporation of water vapor, T s is the surface temperature predicted by the symmetric neural operator, q s is the saturated specific humidity. The flux and temperature field are iteratively coupled: the first iteration uses remote sensing brightness temperature to estimate T s , after getting the flux, input the neural operator to update T s , and then fed back to the flux formula until convergence. The example shows that after two iterations, T s The convergence error is less than 0.2°C, which meets the accuracy requirements of energy balance calculation.

[0067] By incorporating shortwave albedo, radiation field, and turbulent flux into the joint features, the model explicitly perceives the energy state of the debris-ice surface. The symmetric neural operator consists of multiple layers of Fourier integral operators, and the frequency-domain convolution kernels are constrained in a mirror-symmetric manner, structurally ensuring that the heat diffusion matrix is positive definite. The self-attention-symmetric neural operator combination captures multimodal spatial interactions while enforcing local energy conservation. The training loss function consists of a contrastive loss, an energy balance residual, and a variational lower bound; the weights are determined through cross-validation, ensuring that the model maintains physical consistency while also maximizing generalization.

[0068] The root mean square error of the thickness inference of the experimental glacier group on the Pamir Plateau of the present invention is reduced by more than 30% compared with the traditional empirical formula. In areas with significant debris coverage, the change in debris thickness is sensitive to the change in albedo. After the short-wave albedo is added, the model can automatically amplify this sensitivity, and the thickness error decreases relatively more. In addition, the coupled turbulent flux model of the present invention introduces radar estimation of roughness length to achieve dynamic regulation of heat exchange at the ice-debris interface; in the case of a sudden increase in wind speed, the flux peak output by the model and the ground eddy-related observation error are controlled within 15%, while the model error without using radar roughness information is close to 40%. The results of Example 2 show that the energy budget closure of a tributary glacier with a three-hour resolution is increased from 0.79 to 0.93, verifying the effectiveness of the physical driving field parameterization.

[0069] Preferably, the symmetric neural operator is composed of multiple layers of Fourier integral operators connected in series, and the convolution kernels of each layer are set to a symmetrical structure according to the real part and the imaginary part in the frequency domain.

[0070] The symmetric neural operator assumes the core responsibility of the mapping of "physical driving field → surface temperature" in this invention. Its design principle is to let the network structure itself reflect the symmetric positive definite characteristics of the heat diffusion-radiation coupling equation, so as to automatically satisfy the energy conservation constraint during the end-to-end training process without relying on external hard projection or post-processing. The operator uses multi-layer Fourier integral operators as the basic unit. Each layer first completes the convolution operation in the frequency domain, and then returns to the space-time domain through a fast inverse transformation; the convolution kernel is set according to the mirror image of the real and imaginary parts to ensure that the time domain kernel obtained after the inverse transformation has an even function property, thereby keeping the overall convolution matrix symmetric in a discrete sense. The following content is developed in sequence around the principle, implementation and effect.

[0071] The debris-ice-air system can be approximately described as a three-dimensional unsteady thermal diffusion-radiation coupled process at a fine scale. If horizontal flow is neglected and the vertical temperature distribution is simplified to a linear one, the one-dimensional diffusion-source term equation can be obtained:

[0072]

[0073] The symbol T represents the instantaneous surface temperature (unit K), and k represents the effective thermal diffusivity (unit m 2 s-1 ), ρ represents the average density of the debris layer (unit: kgm -3 ), c represents the specific heat capacity (unit Jkg -1 K -1 ), Q net represents the net radiation-turbulent flux source term (unit: Wm -2 ). If the equation is Fourier transformed in the spatial direction, we can get:

[0074]

[0075] It shows that in the frequency domain the diffusion operator presents -κ|f| 2 The multiplier form is proportional to the square of the frequency and is an even function. Based on this, the present invention introduces the Fourier integral operator and defines the convolution kernel directly in the frequency domain, so that the amplitude and sign of the convolution kernel for f and -f satisfy:

[0076]

[0077] Where κ(f) is the frequency domain convolution kernel, denotes the imaginary part. This mirror constraint ensures that the time-domain kernel k(x) obtained by the inverse transform satisfies both k(x) = k(-x) (an even function) and ∫k(x)dx = 0. This ensures that the convolution matrix is symmetric in discrete form and its row sum is zero, consistent with the mathematical characteristics of the heat diffusion matrix. This symmetry ensures that the neural operator naturally adheres to energy conservation when predicting temperature gradients, avoiding the occurrence of physically meaningless negative diffusion.

[0078] In practical applications, the symmetric neural operator is composed of four layers of Fourier integral operators connected in series. Let the input be the vector The operation process of the lth layer (l=1-4) is as follows: 1. Fast Fourier transform: Frequency domain convolution: Where ⊙ is the Hadamard multiplication 3. Inverse transform: 4. Nonlinear mapping: u (l) =σ(v (l) +b (l) ), σ is a linear rectification function. The convolution kernel κ of each layer (l) After being randomly sampled at the beginning of training, the symmetry is enforced according to the above mirror rule; during updating, only the amplitude is allowed to be adjusted synchronously on both sides and the phase is adjusted in opposite directions on both sides to ensure that the symmetry is not destroyed during training.

[0079] In order to balance deep learnability and numerical stability, the present invention adds batch normalization to the output of each layer to make the activation value have a mean of zero and a variance of one, and (l)An upper bound on the amplitude is introduced to prevent high-frequency channels from being amplified infinitely. Unlike conventional convolutional networks, this structure performs convolution in the frequency domain first, allowing it to capture long-range dependencies with O(n log n) complexity. Symmetric constraints further reduce the number of trainable parameters and the risk of overfitting.

[0080] The network input consists of two parts: one is the joint feature output by the cross-modal encoder, and the other is the physical driving vector of the shortwave albedo, downwelling radiation and turbulent flux stack. The two are mapped to the same length as the input grid through the fully connected layer and then added element by element to form u (0) The output is the temperature field T pred and the energy balance source term Q net Compute the physical residual:

[0081]

[0082] Neural operators with global loss:

[0083]

[0084] Joint optimization, and The two objectives are derived from the cross-modal contrast objective and the variational objective of the differentiable probabilistic model. The bidirectional gradient flow requires the symmetric neural operator to satisfy energy conservation while providing sufficient temperature information for the thickness-existence posterior solution.

[0085] The symmetric convolution kernel of the symmetric neural operator ensures that the predicted temperature field is spatially smooth and free of unphysical peaks. Experiments show that compared with ordinary Fourier layers without symmetry restrictions, the number of pixels with negative diffusion (gradient reversal) in the temperature gradient field is reduced by 90%, and the training convergence steps are shortened by about one-third. In scenarios with highly variable thickness of the debris layer in mountainous areas, the model can adaptively enhance the high-frequency kernel amplitude to reduce the impact of thin high-albedo debris layers on excessive temperature rise; in bare ice areas in river valleys, the low-frequency kernel amplitude is enhanced to achieve overall smoothness of the temperature field.

[0086] Example: A complete link of the proposed method was run on a 20-meter grid at the leading edge of the Karakoram Main Ridge glacier, comparing it to a traditional convolution-long short-term memory coupled model. Using a validation set of 360 independent temperature, thickness, and radiation flux samples, the symmetric neural operator scheme achieved a root mean square error (RMS) of 0.94 Kelvin for temperature, compared to over 1 Kelvin for the traditional scheme. The coefficient of determination of the thickness posterior for the measured thickness was improved to 0.87. Due to the symmetry of the convolution kernel, the power spectral density (PSD) revealed a symmetrical and stable heat diffusion pattern, without the spectral leakage often associated with high-order convolution kernels in traditional models.

[0087] Preferably, the differentiable probability model is composed of multiple levels of affine coupled transformation units connected in series, each level of affine coupled transformation unit uses an alternating mask method to separate input channels, and the scale function and the displacement function are implemented through a multi-layer perceptron.

[0088] In the context of glacial moraine identification, thickness and presence are continuous-discrete mixed random variables, requiring both pixel-level point estimation and a reliable uncertainty measure. This paper employs a differentiable probabilistic model to achieve this goal. The core concept is to use a reversible transformation to map a simple prior distribution, which is easy to sample and calculate density, into a complex posterior distribution. The transformation chain is continuous and differentiable, facilitating backpropagation with the preceding symmetric neural operator.

[0089] The model adopts a normalized flow architecture composed of a series of multi-level affine coupled transformation units. Let y be the target random vector, element 1 corresponds to the pixel thickness, and element 2 corresponds to the continuous representation of the moraine existence after logarithmic probability transformation. Define the prior z (0) Satisfies the lognormal-Bernoulli mixed distribution of independent dimensions. The k-th level affine coupling transformation is denoted as:

[0090]

[0091] m is a fixed binary mask, ⊙ represents element-wise multiplication, s k (·) and t k (·) are scale function and displacement function respectively. They are implemented by multi-layer perceptron. The perceptron input is the value of the masked reserved channel and outputs a vector with the same length as the transformed channel. The activation function uses a continuous differentiable rectified linear unit to ensure that the overall transformation Jacobian is only equal to s k Related without t k , which makes it easy to analyze and calculate. The alternating mask method exchanges m and 1-m in adjacent units, so that all channels are used as conditional branches at least once in the entire flow. After connecting L levels in series, a reversible mapping is obtained. The posterior log density is calculated as follows:

[0092]

[0093] in is the transformation channel index, p0 is the prior density, s k,i is the scaled output of the i-th channel at level k. Since the determinant is the product of the exponentials of the diagonal elements, the complexity is linearly proportional to the number of channels and is independent of the spatial dimension, making it suitable for batch inference of large-scale remote sensing rasters.

[0094] To close the loop with the surface temperature and joint features, the differentiable probabilistic model adopts a variational inference framework. The parameters of the mapping f are jointly optimized with the previous encoding-physical network to minimize the evidence lower bound:

[0095]

[0096] Here, x represents the cross-modal observation, θ represents the generative network parameters, and φ represents the flow model parameters. The first term uses the conditional likelihood of the symmetric neural operator output, and the second term is the Kullback-Leibler divergence between the posterior and the prior. Because f is reversible, the KL term can be quickly evaluated in closed-form analysis, ensuring efficient training.

[0097] Model output thickness posterior mean μ h and variance And the existence probability p c Through closed-form gradients, the measured thickness samples obtained by active field sampling can directly fine-tune φ, completing the a priori contraction. Due to the reversible nature of the flow, sampling only requires reverse reasoning of the prior to generate thickness-presence pairs, achieving uncertainty propagation.

[0098] The embodiment uses the western Karakoram glaciers as the test field. A 20-meter resolution grid with a total of 300,000 pixels is used as the training set, and 10% of the pixels are randomly selected as the test set. The true thickness value comes from the ground trench-radar joint profile. The experimental group applies a differentiable probability model; the control group uses a multi-layer perceptron of the same size to directly regress the thickness. The results show that the experimental group's prediction mean error is reduced by 20%, and the linear correlation coefficient between the prediction variance and the absolute error is 0.8, indicating that the variance has a good characterization of the error; the control group's variance-error correlation coefficient is less than 0.3. Further, in the active sampling stage, 1,000 pixels are selected for field measurement and sent back for fine-tuning. The root mean square error of the thickness of the experimental group is reduced by another 10%, and the control group is reduced by less than 5 percentage points, verifying the advantages of the flow model in the incremental learning scenario.

[0099] Another practical benefit of the differentiable probabilistic model is that it provides spatially correlated uncertainty for subsequent risk mapping. The uncertainty field, formed by the weighted combination of thickness variance and presence probability, plays a key role in constructing the objective function for quantum annealing optimization. This allows sampling points to be concentrated in critical areas with high thickness variance and presence probabilities close to 0.5, significantly improving field efficiency.

[0100] An uncertainty field is constructed based on the posterior data, high entropy and spatially correlated pixels are selected through quantum optimization to perform field sampling, the measured thickness is used to fine-tune the differentiable probability model and the symmetric neural operator, and a spatially continuous thickness field is formed through anisotropic graph diffusion regularization;

[0101] The thickness expectation, variance, and presence probability derived from multi-source remote sensing inversion are pixel-level random quantities, whose reliability is affected by training sample imbalance, radiation noise, and model extrapolation error. To minimize this uncertainty while limiting fieldwork costs, this paper proposes a three-step closed loop: "posterior drive-quantum optimization-anisotropic diffusion." First, the posterior statistics are converted into a measurable uncertainty field. A quantum annealing optimization model is then used to select the most informative and spatially representative sampling pixels across the entire domain. Finally, the sampling results are fed back into a differentiable probabilistic model and symmetric neural operators, and the thickness field is smoothed using anisotropic graph diffusion, suppressing local noise while preserving detail. This three-step closed loop ensures that the sampling-update-regularization process is completed in one go, eliminating multiple rounds of fieldwork.

[0102] The uncertainty field is constructed, and the pixel thickness posterior distribution is output by the differentiable probability model, whose parameter is the thickness mean μ h ,variance and existence probability p c The uncertainty index U needs to reflect the size of the variance and also needs to consider the characteristics of the probability of existence being far away from zero or the maximum amount of information at a certain moment. The present invention sets:

[0103]

[0104] The first term is information entropy, and the second term is variance. The two dimensions are directly added after normalization. c It represents the probability that the pixel is a surface moraine, dimensionless; represents the thickness variance in square meters. The gridded uncertainty field {U i}, i is the pixel index. Experiments have shown that using only variance ignores the region at the classification boundary where the model lacks confidence when the probability approaches 0.5, while using only information entropy cannot distinguish differences in thickness fluctuation amplitudes. The combination of the two can account for both continuous and discrete uncertainty.

[0105] Quantum annealing sampling point screening: Traditional greedy or heuristic methods often only consider the uncertainty amplitude and ignore spatial redundancy. This invention incorporates spatial correlation into the Boolean quadratic disordered optimization model and solves it in one go on the quantum annealing machine. The formula is:

[0106]

[0107] Among them, s i ∈{0,1} is a binary variable indicating whether the pixel is selected; a i =U i is the linear coefficient, and pixels with high uncertainty tend to be selected; b ij =-λC ij ,λ is a positive weight, C ijThe larger the amplitude of the convolution kernel spectrum of the symmetric neural operator at the pixel pair (i, j), the stronger the temperature coupling. It is hoped that the two will not be selected at the same time to avoid information redundancy.

[0108] Quantum annealing obtains s by finding the global minimum solution on the energy function H i Configuration. Compared with classical hill climbing or simulated annealing, quantum tunneling mechanism can skip local minima with shorter annealing time in exponential combination space. Get the target set S = {i|s i =1}, the shortest flight path of the UAV is generated and the flight is carried out according to the pixel center.

[0109] On-site sampling and model fine-tuning, UAV equipped with surface wave radar and dual-band thermal infrared sensor. The radar profile is inverted by time-depth transformation to measure the thickness h * , thermal infrared brightness temperature is used to verify the radiation correction. *} as a new supervision pair, fine-tuning the differentiable probability model and the symmetric neural operator with a small learning rate: the symmetric neural operator receives h * Update the top temperature boundary condition and use h for normalized flow * Recalculate the weights of each scale-shift perceptron. Thanks to the end-to-end differentiability of the model, fine-tuning does not destroy the symmetric structure of the existing convolution kernel and converges with only a limited number of iterations.

[0110] The anisotropic graph diffusion regularization, after fine-tuning, still contains random noise in the new thickness field, requiring spatial regularization. The thickness of the glacier moraine is smooth in the direction of slope and drainage channel, and has a large gradient in the direction perpendicular to the ice flow line. The present invention constructs an undirected graph G = (V, E), where the node V corresponds to the pixel and the edge weight is:

[0111]

[0112] Among them H i is the node altitude, r ij is the pixel spacing. Diffusion tensor:

[0113]

[0114] The symbol η is the uncertainty weight, which increases the diffusion coefficient in the direction of high uncertainty. Using explicit difference:

[0115]

[0116] Iterate until the average absolute difference between the two adjacent steps is less than the threshold. Symbol represents the thickness at iteration t, and λ is the total variation weight, which is used to suppress stair-step artifacts. This diffusion-total variation coupling preserves cross-section abrupt changes in high-slope regions and smoothes noise in low-slope regions, ultimately outputting a spatially continuous thickness field.

[0117] In this example, the Karola Glacier was used as a demonstration area, with a total of 320,000 grid cells selected. The mean of the posterior uncertainty field was 1, and the variance was 1.49. The quantum annealing model assumed 2,000 variables and 1,500 constraints, returning the Boolean vector with the minimum energy in a single annealing operation. The number of sampling points accounted for 1% of the total pixels; after the thickness was measured in the field, the model was fine-tuned three times, and the root mean square error of the thickness was reduced from 2.1 meters to 1.7 meters. After diffusion regularization, the area proportion of small-scale noise collapse patches was compressed from 8% to 3%. Compared with the entropy threshold + random sampling scheme, the quantum optimization method of the present invention reduced the thickness error by an additional 12% with the same field workload.

[0118] At the source of the Golmud River in Qinghai Province, the experimental team compared the smoothing with a classic isotropic Gaussian filter. The results showed that the Gaussian filter resulted in excessive blurring of the ice tongue boundary, with boundary errors exceeding five pixels. The anisotropic diffusion method of the present invention prevented the cross-slope diffusion of thickness values in the slope direction, keeping the boundary error within one pixel, validating the effectiveness of the slope-uncertainty coupled tensor design.

[0119] Through the collaboration of uncertainty measurement, quantum annealing optimization and anisotropic diffusion, the present invention achieves posterior shrinkage and thickness field smoothing under the premise of limited field sampling costs, which not only fully utilizes high-uncertainty information but also suppresses spatial redundancy; and through end-to-end fine-tuning, the temperature-thickness-energy balance link is ensured to be closed, which significantly improves the credibility and spatial continuity of glacier surface moraine identification results.

[0120] Preferably, the uncertainty field is obtained by linearly combining the information entropy of the pixel surface presence and the pixel thickness variance according to preset weights.

[0121] The pixel-level glacier moraine identification output is a continuous-discrete mixed posterior: the thickness h follows the conditional probability density p(h|x), the presence c follows the Bernoulli distribution Ber(p c To prioritize locations with the least information while limiting field costs, this binary posteriori needs to be mapped into a single, spatially measurable metric field. The uncertainty field proposed in this paper is precisely such a metric. It unifies the uncertainty of discrete classification and the uncertainty of continuous regression into the same dimension and can indicate the direction of the steepest "missing information" in a gradient sense.

[0122] First, calculate the existence information entropy. For the Bernoulli variable c, the entropy is defined as:

[0123] H c =-p c ln p c -(1-p c )ln(1-p c )

[0124] where p cIt represents the probability that the pixel belongs to the surface moraine category, dimensionless. c = 0.5, the maximum value is obtained, reflecting the position where the model lacks the most classification confidence; when p c When it approaches 0 or 1, the entropy decays to zero, and the corresponding classification is stable.

[0125] Second, the thickness posterior variance is directly read from the differentiable probability model It measures the dispersion of thickness estimates under given observations. It usually increases, indicating that the model is still uncertain in the regression sense.

[0126] The entropy term and the variance term have different dimensions. Directly adding them together will cause one term to dominate. In order to align the scales, the present invention first performs minimum-maximum normalization on the variance:

[0127]

[0128] in and is the extreme value of the variance of the entire image. After normalization It is consistent with the numerical range of the entropy term.

[0129] The uncertainty field U is given by:

[0130]

[0131] ω H With ω σ is the preset weight, satisfying ω H +ω σ =1. The present invention takes equal weights by default. When the test set shows that the classification error is significantly greater than the thickness error, ω can be appropriately increased through cross-validation. H , on the contrary, increase ω σ .

[0132] This linear superposition has two important properties. First, In both the entropy-dominated and variance-dominated regions, it points to the pixels with the fastest information gain, which is conducive to continuous optimization. Secondly, the choice of ω will not change the sparsity pattern of rank(U), and the sampling order is stable.

[0133] At the implementation level, the computational process and posterior inference share GPU tensors: for each batch of pixels, the normalized flow is decoded to obtain μ h , p c , and then calculate U by tensor addition. The whole process only contains logarithmic, exponential and normalization operations, and remains differentiable. When interacting with the quantum annealing platform, it is only necessary to project U into the first-order coefficient {a i}, and use the convolution spectrum to generate the secondary coupling coefficient {bij}, a quadratic unordered Boolean optimization model can be constructed.

[0134] Example 1: After generating a U field in a clastic-rich area in the western Kunlun Mountains, three sampling schemes were compared with equivalent field time budgets: random sampling, variance-only sampling, and entropy-variance fusion. The random scheme achieved an average thickness root mean square error of 2.3 m; the variance-only scheme reduced this to 1.9 m; and the entropy-variance fusion scheme further reduced this to 1.6 m, while the classification F1 score improved by 0.07. This demonstrates the benefits of multi-source uncertainty fusion for both continuous and discrete targets.

[0135] Example 2: Fifty drilling points were sampled for both entropy-dominated and variance-dominated pixels in the interlaced shadow and bare ice region of Mount Gongga. After updating the entropy-dominated samples, the endpoints of the existence probability converged, and the classification miss rate decreased by 40%. After updating the variance-dominated samples, the median thickness variance decreased by half. The uncertainty field exhibited a spatially mixed "banded-patchy" pattern, consistent with the aspect-trough network, demonstrating that it inherits energy-topography control while accurately exposing model blind spots.

[0136] The uncertainty field of this invention also contributes to the construction of an anisotropic diffusion tensor. By setting the vertical diffusion coefficient to 1 + ηU and the lateral coefficient to 1, diffusion is automatically enhanced in locations with steep slopes and high U, eliminating noise peaks. In low-uncertainty regions, the uncertainty field remains unchanged to preserve details. Spectral analysis shows that the introduction of adaptive uncertainty diffusion reduces the system's largest eigenvalue and reduces the number of iterative convergence steps by approximately one-third.

[0137] In summary, the entropy-variance linear fusion provides an uncertainty measurement that emphasizes both physics and statistics: the entropy term captures the fuzziness of classification boundaries, and the variance term characterizes the divergence of thickness regression; normalization and weight adjustment ensure a balance between the two; and quantum annealing and diffusion regularization achieve a closed loop of sampling-updating-smoothing on this basis, greatly improving the credibility, spatial continuity and field efficiency of glacier moraine identification results.

[0138] Preferably, quantum optimization uses quantum annealing to solve the quadratic disordered Boolean optimization model, sets the pixel uncertainty as the linear term coefficient, and sets the symmetric neural operator convolution kernel spectral similarity as the quadratic term coupling coefficient to obtain a set of pixels with high entropy and strong spatial correlation.

[0139] Quantum annealing has recently gained widespread recognition as a hardware-accelerated approach for solving large-scale binary combinatorial optimization problems. This paper transforms the multi-source remote sensing-driven glacier moraine sampling task into a Quadratic Unconstrained Binary Optimization (QUBO) model. Leveraging quantum annealing hardware, it outputs an optimal set of pixels in one go, achieving the dual objectives of "maximizing uncertainty" and "minimizing spatial redundancy." The following describes this approach from four perspectives: modeling principles, coefficient construction, quantum mapping, interpretation, and field results.

[0140] Modeling principle: Assume that the grid to be evaluated has N pixels. Define binary variables:

[0141] s i ∈{0,1},i=1,…,N

[0142] where s i =1 means the i-th pixel is selected into the field sampling set, s i =0, it is not selected. The target Hamiltonian is written as:

[0143]

[0144] a i is the linear term coefficient, corresponding to the pixel uncertainty; b ij is the quadratic coupling coefficient, corresponding to the spatial correlation of pixel pairs; H is the energy function to be minimized.

[0145] The quantum annealing hardware searches for the ground state on the energy surface through quantum tunneling, and the returned {s i} makes H obtain the minimum value. If the coefficient of the first term is positive, the energy minimization tends to choose a i If the quadratic coupling coefficient is negative, it tends to avoid selecting highly correlated pixels at the same time, thereby reducing spatial redundancy.

[0146] Coefficient construction, linear coefficient a i Directly take the uncertainty field:

[0147] a i =U i

[0148] in

[0149]

[0150] H c,i is the information entropy of the existence of moraine in the pixel table, is the normalized thickness variance. Weight ω H With ω σ The quadratic coupling coefficient b is determined in advance in the validation set and kept fixed.ij The convolution kernel spectrum similarity of the symmetric neural operator is obtained:

[0151] b ij =-λC ij

[0152] Among them, λ>0 is the spatial penalty factor; C ij is the spectrum similarity, which is calculated as follows:

[0153]

[0154] F SNO Represents the frequency-domain amplitude matrix of the convolution kernel of the symmetric neural operator, normalized to the range [0, 1]. Large spectral amplitudes indicate strong temperature coupling; setting a negative sign causes quantum annealing to disperse the sampling of pixels with strong coupling to avoid repeated acquisition of similar information.

[0155] Quantum mapping and hardware solution, the QUBO coefficient matrix Q is through the relationship:

[0156] Q ii =a i ,Q ij =b ij (i <j)

[0157] Fill in the weight format accepted by the quantum annealer. The present invention adopts the following strategy when embedding:

[0158] Weight compression: a i with b ij The maximum absolute value item is uniformly linearly scaled to the dynamic range allowed by the hardware.

[0159] Chain strength adaptation: For long-distance couplings where the physical adjacency graph is insufficient, virtual chains are used to connect logical variables, and the chain strength is set to twice the maximum value of the secondary coupling to prevent chain breakage.

[0160] Annealing schedule: select single annealing time t ann The quantum temperature window with the minimum injection noise is empirically shown to be ann The base state hit rate is highest when the hardware limit is reached.

[0161] After the hardware completes one annealing, it outputs multiple candidate solutions and selects the solution with the lowest energy as the final sampling solution. If the number of sampling points is constrained by the budget, a penalty term can be added to H or the returned solutions can be sorted by scores, retaining the top K highest s i The corresponding pixels can be used.

[0162] Post-processing and closed-loop correction, after field work is completed, high-quality thickness measurement Written into the training set:

[0163] Differentiable Probabilistic Model Fine-tuning— Append to the stream model input and continue backpropagation for 10 rounds with a small learning rate;

[0164] Symmetric neural operator correction - Update to the boundary condition, recalculate the convolution kernel spectrum, and update b accordingly ij At this point, the uncertainty field is reassessed and the next round of sampling begins; however, in actual projects, the error has significantly converged after one sampling, and the loop can be terminated when the threshold is met.

[0165] Example, experimental area: 45 square kilometers of glaciers on the Pamir Plateau.

[0166] Comparison schemes included random sampling, greedy maximum entropy sampling, and spatial k-means-based partitioned sampling. Metrics included thickness root mean square error (RMSE), classification F1 score, and field track length. The experimental results are shown in Table 1.

[0167] Table 1

[0168] plan Number of sampling points Track length km Thickness RMSEm F1 value random 400 65 2.4 0.78 Maximum Entropy 400 72 2.1 0.81 Partitioned k-means 400 60 2.0 0.80 The present invention 400 58 1.7 0.85

[0169] The results show that, under the premise of the shortest track length, the scheme of the present invention reduces the thickness RMSE by 0.7 meters compared with random sampling and by 0.4 meters compared with maximum entropy sampling; the classification F1 value is also significantly improved, indicating that the screening strategy of entropy and spectrum coupling can better take into account continuous and discrete uncertainties than simple information entropy or spatial partitioning.

[0170] Example A: Debris accumulation area on a high slope. The uncertainty peak is concentrated in the transition zone from the foot of the slope to the top of the slope. After quantum annealing, radar measurements confirmed that about 70% of the pixels are mixed layers of debris and bare ice. After the model was updated, the local H c a 60% decrease, A 45% drop.

[0171] Example B: Low-slope bare ice area. Spectral coupling matrix C ij The error is higher in the low-frequency channel, and the annealing results tend to be sparse and scattered. After the update, the thickness error is low at multiple ice tongue peaks. The control group without the treatment of the present invention still shows stripe noise at the same number of points.

[0172] Preferably, on-site active sampling collects the pixel set in sequence through a traveling salesman path constructed by a drone. The drone is equipped with a surface wave radar to obtain a radar profile and invert the thickness according to predetermined electromagnetic velocity parameters. At the same time, it is equipped with a dual-band thermal infrared sensor to obtain brightness temperature information. The radar inversion thickness is used to fine-tune the differentiable probability model and the symmetric neural operator.

[0173] The fundamental task of the active sampling process is to introduce a small number of high-precision in situ thickness observations, without sacrificing the macroscopic coverage advantage of remote sensing, to shrink the posterior distribution and calibrate the temperature-energy inference link. This invention integrates the field work process into the entire differentiable pipeline: first, quantum annealing is used to generate a representative set of pixels, then a traveling salesman path is used to generate the shortest path, and finally, the measured thickness is written back into the differentiable probabilistic model and symmetric neural operator as a pseudo-label that can be propagated through the gradient. This "task-driven-path optimization-differentiable reflow" design significantly reduces the quadratic error caused by the traditional separation of field and office work.

[0174] The UAV platform uses a multi-rotor model with a sufficient thrust-to-load ratio and is equipped with two types of sensors: a ground-penetrating radar (GPR) and a dual-band thermal infrared sensor (TIR). The GPR is responsible for obtaining the round-trip time delay between debris, ice surface, and ice base; the TIR is responsible for synchronously recording the surface brightness temperature, which serves as the measured boundary condition for the symmetric neural operator. The trajectory planning adopts the traveling salesman problem (TSP) model in the graph theory sense, where the nodes are the pixel center coordinates output by quantum annealing, and the edge weights are geodetic distances. After obtaining the Hamiltonian shortest circuit through integer linear programming, the UAV flight instruction set is generated. The instructions include longitude and latitude, a preset altitude, and hovering time to ensure that the sensor coverage window matches the radar transducer pulse repetition rate.

[0175] Radar thickness inversion is based on the propagation constant of electromagnetic waves in the debris-ice layer. The two-way delay Δt is automatically picked up by the radar controller as the time difference between the first and third interface echoes, and is output in real time by the onboard computing unit after image stabilization filtering. The thickness d is calculated according to the formula:

[0176]

[0177] where v is the velocity of electromagnetic waves in the measured medium; the effective dielectric constant ε of the conventional debris-ice mixed layer r The speed estimation formula is obtained through the pre-calibration experiment. (c is the speed of light in a vacuum) is hard-coded into the firmware. This approach avoids multiple manual adjustments in the field and improves the consistency of thickness calculations. To reduce the diffuse reflection mismatch at the ice-debris interface, the radar antenna attitude controller maintains a small inclination angle relative to the surface normal during flight, correcting for delay based on the real-time attitude angle.

[0178] The TIR sensor simultaneously collects brightness temperatures in two atmospheric window bands, calibrates them radiometrically, and writes them to the temperature stack for immediate feedback into the convolutional layer of the symmetric neural operator. If the TIR brightness temperature deviates from the remote sensing temperature estimate for the same pixel by more than a threshold, the system automatically marks that pixel for entry into the next round of quantum annealing, achieving closed-loop bootstrapping.

[0179] After the sampled data is returned to the internal server, two types of fine-tuning are triggered. The first type: Differentiable probability model with For new supervised pairs, the early affine coupling blocks are frozen, and only the weights of the last two levels of scaling and displacement functions are adjusted. The learning rate is set to 1% of the initial value to avoid disrupting the overall distribution. The second type: a symmetric neural operator uses the measured thickness as the depth boundary, recalculates the surface temperature-energy residual, and performs a small gradient descent on the convolution kernel amplitude to make the local diffusion-radiation coupling more accurate to the site.

[0180] The thickness grid obtained after fine-tuning often shows sudden drops or rises near the sampling points. To suppress interpolation ringing, the present invention introduces anisotropic graph diffusion regularization. The main diagonal elements of the diffusion tensor are selected as the weighted sum of the slope and uncertainty, and the lateral component remains at the base value. In the explicit difference format, the iteration step size is adaptively determined by the spectral radius, and the iteration termination condition is that the mean absolute difference of the pixels is less than a given threshold. Because the diffusion tensor has direction dependence, the diffusion in the high slope area is suppressed, and the ice tongue and sidewall details are preserved; the diffusion in the high uncertainty area is increased, and the noise peak tends to be degenerate.

[0181] Example: A 28-kilometer drone flight path was deployed in the Nyainqentanglha Mountains, along a high-risk section of ice-locked lakes. Quantum annealing outputted 960 sample pixels, and after TSP optimization, the track length was reduced to 23 kilometers, an 18% reduction compared to manual zoning. This was accomplished with just two battery changes for a single drone. Radar-derived thickness was verified through trenching, achieving a mean absolute error of 0.7 meters. After fine-tuning, the root mean square error (RMS) of thickness was reduced from a baseline of 2 meters to 1.4 meters. After 30 iterations of diffusion regularization, the noise patch ratio was reduced from 10% to 3%, and the ice edge curve and the digital elevation model of the high-resolution stereo image pair had an error of less than one pixel.

[0182] Compared to the classic "fixed grid drilling-kriging interpolation" process, this method increases the number of verification points by an order of magnitude and reduces errors by nearly half, while maintaining the same field hours. Compared to airborne optical-radar remote sensing inversion alone, it avoids systematic underestimation in areas with a strong debris-ice mixture. This demonstrates the efficiency and reliability of the quantum annealing-TSP-anisotropic diffusion closed loop: the annealing model ensures optimal sampling points in both information and space dimensions; path planning conserves flight paths; radar and thermal infrared sensors complement each other in measurement; and micro-adjustment and diffusion smoothing allow field information to be accurately integrated into the global posterior in the form of gradients, achieving iterative convergence for remote sensing-field integration.

[0183] Preferably, the anisotropic graph diffusion regularization model establishes an undirected graph structure on a uniform resolution grid, the main diagonal elements of the diffusion tensor are set according to the weighted values of the pixel slope and the pixel uncertainty, and the explicit difference format is iterated until the average relative difference between two adjacent thickness fields is less than a preset threshold, and the spatially continuous thickness field is output.

[0184] After active sampling is completed and the measured thickness is fed back to the model, the pixel-level thickness field still exhibits three types of local defects: the first type is random noise peaks, which are common at locations where the debris roughness changes sharply; the second type is strip-like oscillations, which are caused by systematic errors in the directionality of the radar inversion path; the third type is step-like fractures, which often occur in slope break zones where the ice-debris thickness changes suddenly and the model is not well fitted. If isotropic smoothing (such as Gaussian filtering or conventional total variation) is directly used, the real ice edge and ice surface microtopography will be overly blurred, affecting the subsequent melting evolution and risk classification. The present invention proposes an anisotropic graph diffusion regularization model to perform post-processing on the thickness field, so that the noise is weakened while retaining the spatial details controlled by the terrain.

[0185] The principle comes from the anisotropic diffusion equation:

[0186]

[0187] τ represents the pseudo time step; d represents the pixel thickness, in meters; A represents the second-order diffusion tensor; λ represents the total variational weight, dimensionless. This achieves direction-dependent smoothing, and the second term suppresses high-frequency oscillations when the gradient modulus is less than the noise threshold. The implementation steps are as follows:

[0188] 1. Graph structure construction: for each pixel v on a uniform 20-meter resolution grid i Create a graph node with eight-neighborhood connections. Edge weight:

[0189]

[0190] H i With H j are the altitudes of the two nodes respectively; r ij is the Euclidean distance between pixels. In this way, diffusion is suppressed between pixels with large elevation differences or long distances.

[0191] 2. Definition of diffusion tensor, for node v i Construct a diagonal tensor:

[0192]

[0193] U i is the uncertainty field value, unitless; η is the uncertainty magnification factor. The longitudinal component is introduced into U iAllowing high uncertainty pixels to release noise faster in the longitudinal direction and keeping the lateral component at the base value of one can avoid excessive lateral blurring. i If it is higher than the set threshold, then multiply the scaling factor k=S i / S max Suppress upslope diffusion on steep slopes.

[0194] 3. Explicit difference discretization, the five-point difference format is written as:

[0195]

[0196] w is the row width. The time step satisfies the stability condition:

[0197]

[0198] An empirical safety factor of 0.45 is usually taken.

[0199] Fourth, the iteration termination criterion is defined as the average relative difference between two adjacent steps:

[0200]

[0201] δ is a very small positive number to prevent the denominator from being zero. If ε (t) <ε stop That stops the iteration and outputs d # =d (t+1) Threshold ε stop Take 10 -3 It removes noise while maintaining slope break details.

[0202] In this example, in the glacier tongue region of the Qilian Mountains, the unnormalized thickness field exhibited striped noise along the radar track. After running 30 steps of anisotropic diffusion, the peak noise amplitude was reduced by 90%, and the glacier tongue boundary error was kept within one pixel. In contrast, using an isotropic Gaussian filter with a three-pixel FWHM processing also reduced noise, but the boundary blur width was expanded to five pixels. This method achieves a compromise between boundary clarity and noise suppression.

[0203] Another experiment was conducted on a glacier front in the eastern Pamir Plateau, a region with intense debris coverage. Applying slope suppression preserved the thickness gradient of the high-slope stratified wall. However, setting η to zero resulted in pseudo-smoothing along the high slopes, increasing the thickness shaving error by 30%. This demonstrates the effectiveness of the slope-uncertainty coupled diffusion tensor for complex terrain.

[0204] The symmetric neural operator with the spatially continuous thickness field and frozen weights is used to fuse the daily weather forecast to calculate the net radiation of the pixel, and the future thickness is obtained by integrating it to a set number of years according to the law of conservation of mass, and the burst risk index raster data is generated based on the presence of surface moraines and the future thickness.

[0205] After diffusion of the anisotropic map, the spatially continuous thickness field d0(i) eliminates random noise and retains the slope details at the grid scale, which can be used as the initial condition for thermal-mass evolution. The subsequent steps use the symmetric neural operator with frozen weights. The daily meteorological driving force is mapped to the pixel net radiation, and the future thickness is then integrated according to the conservation of mass. A burst risk index grid is constructed based on the presence of surface moraines. This process not only complies with the ice-debris energy budget equation, but also avoids the regional extrapolation error of the traditional empirical daily factor. The net radiation is calculated for each pixel i and daily time step t:

[0206]

[0207] Among them, α i is the shortwave albedo; and are the downlink shortwave and longwave radiation of the meteorological model, respectively; ε is the longwave emissivity of the debris-ice mixing surface; σ is the Stefan constant; is the instantaneous surface temperature output by the symmetric neural operator;

[0208] The friction velocity method is used to solve the near-surface meteorological factors. The convolution kernel spectrum of the symmetric neural operator is frozen in the previous step, and the parameters are no longer updated to ensure that the physical consistency constraints are not violated by subsequent fine-tuning.

[0209] when Ablation occurs when the latent heat consumed by melting and the thickness change are consistent with:

[0210]

[0211] ρ represents the equivalent density of the debris layer, L f represents the latent heat of melting of ice, and Δt is the length of one day. = net freezing, and the thickness remains unchanged, because the debris-ice mixed layer dissipates heat mainly through radiation at low temperatures and has no significant source of thickening. Let the daily cycle index t = 1,…,T (T corresponds to 1 year, 5 years or 10 years), then the future thickness is:

[0212]

[0213] The algorithm outputs a snapshot TIFF every thirty time steps so that the calculation can be resumed when unexpected long-period integration occurs.

[0214] The probability of existence of surface moraine p c,i The risk index is derived from the final posterior of the differentiable probability model. The risk index is constructed based on the debris-ice dam failure mechanism: the thinner the debris layer and the higher the presence, the greater the possibility of a thermal channel forming and initiating a failure. The index form:

[0215]

[0216] satisfy When it approaches 1, p c,i →0 tends to 0. The continuous index is better than the hard threshold classification. The color band can be directly rendered in GIS and superimposed with the administrative boundary to achieve hierarchical inspection.

[0217] Preferably, the collapse risk index is calculated based on the exponential relationship between the presence of moraines on the pixel surface and the thickness at a future time. The thickness at a future time is obtained by integrating the mass conservation equation daily with the preset average density of the debris layer and the latent heat constant as parameters. The risk index is output in a raster form for risk grading.

[0218] Moraine dam instability often begins with a process of "penetration-thinning" of the debris layer: as the surface moraine layer thins annually, internal thermal conductivity increases, pore ice is exposed, and hydrothermal pathways penetrate to the backwater surface of the dam, triggering erosion and rapid collapse. Therefore, risk assessment must consider both the presence of debris and its sufficient thinness. After completing multi-year thermal-mass integration, this method maps these two types of information into a single-valued raster—the burst risk index R—which directly serves as a GIS warning classification tool.

[0219] Thickness d at the future moment T Using the mass conservation equation:

[0220]

[0221] Among them, d t Indicates the current thickness of the pixel, in m; Represents the net radiation-turbulent flux, unit Wm -2 ρ represents the debris-ice mixture equivalent density, unit: kgm -3 ;L f Indicates the latent heat of melting of ice, unit Jkg -1 For a full year T1 step, five years T5 steps or ten years T 10 Step by step, integrate daily to get the target thickness d T .

[0222] Surface moraine existence degree p c , comes from the posterior of the differentiable probability model, and its value range is [0,1]. c →1 means the moraine is definitely present, p c →0 means it definitely does not exist.

[0223] Exponential mapping, based on fragility theory, considers thin-layer debris with high presence as a dangerous area and uses an exponential decay function:

[0224]

[0225] where d minis the minimum thickness constant to prevent division by zero. Function properties: when d T Decrease or p c As p increases, R approaches 1 rapidly; when p c →0 or d T When is large, R approaches 0.

[0226] In practical applications, data alignment is achieved: the future thickness grid and presence grid will have the same 20-meter resolution, UTM coordinate system, and consistent pixel indexing. Parallel computing is achieved by using the Tensor framework for one-time vectorized operations, completing millions of pixels in milliseconds on a GPU.

[0227] Thresholds are graded and divided into three default levels: green (R<0.3): routine inspections; orange (0.3≤R<0.6): seasonal inspections; and red (R≥0.6): on-site inspections by dedicated teams. The thresholds can be fine-tuned based on historical incident reviews.

[0228] Raster output: GeoTIFF 32-bit floating point with color scale table, which can be directly loaded into QGIS and ArcGIS.

[0229] The practical application effects of the present invention are shown in Table 2:

[0230] Table 2

[0231] index Experience-Day Model The present invention's integral-index link Thickness RMSE (m) 2.1 1.4 Historical breaking point hit rate 53％ 81％ Calculation duration (10 years) 2h 12min

[0232] The results show that: the thermal-mass integral combined with the physical driving field significantly reduces the thickness error; the exponential mapping amplifies the weight of thin-layer and high-presence pixels, greatly improving the hit rate; GPU batch processing makes the ten-year rolling forecast take only more than ten minutes, which can meet the annual update needs.

[0233] The above are merely embodiments of the present application and are not intended to limit the present application. For those skilled in the art, the present application may have various modifications and variations. Any modifications, equivalent replacements, improvements, etc. made within the spirit and principles of the present application should be included within the scope of the claims of the present application.

Claims

1. A method for identifying glacial moraines based on multi-source remote sensing parameters, characterized in that: The following steps are involved: Synchronously acquire optical and thermal infrared remote sensing images, synthetic aperture radar remote sensing images and digital elevation data, perform radiometric-atmospheric joint correction and geometric-reprojection processing, and generate physically consistent data; The physical consistency data is input into a cross-modal encoder, and combined with the physical driving field composed of shortwave albedo, downwelling radiation and turbulent flux, the surface temperature is inferred using a symmetric neural operator, and the temperature and joint characteristics are used as conditions to input into a differentiable probability model to obtain the posterior data of surface moraine thickness and surface moraine presence; An uncertainty field is constructed based on the posterior data, high entropy and spatially correlated pixels are selected through quantum optimization to perform field sampling, the measured thickness is used to fine-tune the differentiable probability model and the symmetric neural operator, and a spatially continuous thickness field is formed through anisotropic graph diffusion regularization; The symmetric neural operator with the spatially continuous thickness field and frozen weights is used to fuse the daily weather forecast to calculate the net radiation of the pixel, and the future thickness is obtained by integrating it to a set number of years according to the law of conservation of mass, and the burst risk index raster data is generated based on the presence of surface moraines and the future thickness.

2. The method according to claim 1, characterized in that The cross-modal encoder adopts a visual transformation network. The optical and thermal infrared remote sensing images are divided into non-overlapping patches of fixed length at the input and linearly mapped to form a first feature sequence. The synthetic aperture radar remote sensing images and digital elevation data are divided into the same length to form a second feature sequence. The first and second feature sequences are processed by multi-layer self-attention modules respectively and then spliced into a joint feature vector in the hidden dimension.

3. The method according to claim 1, characterized in that The shortwave albedo is calculated based on the normalized difference between the reflectivity of the green light band and the reflectivity of the shortwave infrared band. The downlink radiation is interpolated from the meteorological forecast field to a spatial resolution consistent with the physically consistent data. The turbulent flux is calculated based on the friction velocity method based on the near-surface meteorological elements.

4. The method according to claim 1, wherein The symmetric neural operator is composed of multiple layers of Fourier integral operators connected in series, and the convolution kernels of each layer are set to a symmetrical structure according to the real and imaginary parts in the frequency domain.

5. The method according to claim 1, wherein The differentiable probabilistic model is composed of multiple levels of affine coupled transformation units connected in series. Each level of affine coupled transformation unit uses an alternating mask method to separate the input channels, and the scaling function and displacement function are implemented through a multi-layer perceptron.

6. The method according to claim 1, characterized in that The uncertainty field is obtained by linearly combining the information entropy of the pixel surface moraine existence and the variance of the pixel thickness according to the preset weights.

7. The method according to claim 1, characterized in that Quantum optimization uses quantum annealing to solve the quadratic disordered Boolean optimization model, sets the pixel uncertainty as the linear coefficient, and sets the spectral similarity of the symmetric neural operator convolution kernel as the quadratic coupling coefficient to obtain a set of pixels with high entropy and strong spatial correlation.

8. The method according to claim 1, characterized in that Active on-site sampling collects the pixel set in sequence through a traveling salesman path constructed by a drone. The drone is equipped with a surface wave radar to obtain the radar profile and invert the thickness according to predetermined electromagnetic velocity parameters. At the same time, it is equipped with a dual-band thermal infrared sensor to obtain brightness temperature information. The radar inverted thickness is used to fine-tune the differentiable probability model and the symmetric neural operator.

9. The method according to claim 1, characterized in that The anisotropic graph diffusion regularization model establishes an undirected graph structure on a uniform resolution grid. The main diagonal elements of the diffusion tensor are set according to the weighted values of the pixel slope and the pixel uncertainty. The explicit difference format is used to iterate until the average relative difference between two adjacent thickness fields is less than a preset threshold, and the spatially continuous thickness field is output.

10. The method according to claim 1, characterized in that The burst risk index is calculated based on the exponential relationship between the presence of moraines at the pixel surface and their thickness at future times. The thickness at future times is obtained by daily integration of the mass conservation equation with the preset average density of the debris layer and the latent heat constant as parameters. The risk index is output in raster form for risk grading.

Images